[C-L]

Hi all,

I have recently been using LES to simulate a turbulent flow at Re=180 on a hex mesh. The results have compared really well to DNS results for rms and mean profiles but I'm having issues when extending the problem to polyhedral meshes. The figures below show the profile (left is on the polymesh and right is the hexmesh).

The results look reasonably similar in the ZY plane but there are large spikes in velocity in the XZ plane on the forced cyclic boundary. These spikes don't exist anywhere else so I think it is due to the fvOptions function that I am using. The mesh was generated using starccm+ and converted to OF using ccm26ToFoam and createPatch to match the two cyclic boundaries in Z and X.

Has anyone else experienced something similar/ has any suggestions as to why this is occurring?

My case file is set up almost exactly as the channel395 tutorial but I have changed the viscosity to 1/180 and the Ubar to 15.55556 to make Re_b = 2800, Re_tau = 180(ish) and the pressure gradient roughly equal to 1. I am also using the Lagrangian Dynamic LES model.

Thanks in advance,
Charlie


FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
momentumSource
{
type meanVelocityForce;
active yes;

meanVelocityForceCoeffs
{
selectionMode all;

fieldNames (U);
Ubar (0 0 15.5555556);
}
}

[Santiago]

Several could be the reasons why the 《red spots》 in the flow normal planes. The first is whether the two faces are topologically identical. Other if the volumes that share these faces are hexas or otherwise; if the face-normals of the internal faces are very slanted wrt the mean flow direction, there's a tendency to increase mass flux through the boundary face. Finally try to use the cubeVolRoot filter as the second filter, you need to be sure that the grid is not producing excessive oscillations over the bc's. What are you using as discretization for the convective term?

Sent from my GT-I8190L using CFD Online Forum mobile app

[C-L]

Thanks for the informative reply. I will provide more details when I am in the office tomorrow but I can answer some of those questions now:

-the mesh faces are topologically identical - the mesh nodes/faces match on both boundaries to a higher tolerance. (Default is 10e-3 I think?).

- the volumes are likely to be mixed hex and poly but 95% of the mesh is poly. All the cells have a reasonable quality/skewness.

- the face angles on the internal faces of the boundary is not something that I have looked at but are almost certainly not normal to the flow direction. I think this would be tough to specify in the meshing process? I will certainly investigate this tomorrow.

- I will detail the fvSchemes tomorrow but I'm fairly sure the convective term is Gaussian linear. I'm not sure where the second filter is specified?

Thanks for the reply

Sent from my SM-G920F using CFD Online Forum mobile app

[Santiago]

The filter that you specify in LESProperties is, in fact, the second filter if dynamic procedures are used. Have in mind that FOAM is a library thought to be used for MILES, thus the primary filter is the grid itself

Sent from my GT-I8190L using CFD Online Forum mobile app

[C-L]

Originally Posted by Santiago

Several could be the reasons why the 《red spots》 in the flow normal planes. The first is whether the two faces are topologically identical. Other if the volumes that share these faces are hexas or otherwise; if the face-normals of the internal faces are very slanted wrt the mean flow direction, there's a tendency to increase mass flux through the boundary face. Finally try to use the cubeVolRoot filter as the second filter, you need to be sure that the grid is not producing excessive oscillations over the bc's. What are you using as discretization for the convective term?

Sent from my GT-I8190L using CFD Online Forum mobile app

I have just had a look at my case files and the matchTolerance for the cyclic boundary conditions is 0.0001. I am using Gauss linear for the advective terms in the U equation but Gauss limitedLinear 1 for the transport equations in the LES model (flm and fmm for dynamic Lagrangian). I'm still not sure where to find the cubeVolRoot filter: Is this specified as the delta function in the LES model? If so I am currently using vanDriest damping.

[C-L]

Hi, I just got that message after I posted that reply. Should I swap vanDriest to cubeRootVol? Thanks for the help

[Santiago]

You need to have TDV/NVD algorithms for the convective term if you want to avoid oscillations due to the grid. Yes, delta defines the filter and 'delta' for the convolution

Sent from my GT-I8190L using CFD Online Forum mobile app

[C-L]

Ah okay, so should I specify Gauss filteredLinear as the div(phi,U) term?

I have now checked the mesh face orthogonality and it appears that most of the large spikes do occur at on cells that have faces not normal to the mean flow direction (see the attached photo). If this is the main reason why I am getting issues on the PBC I'm not sure I can fix it using meshing software! Eventually I will be simulating more complex geometries so polyhedral cells would be most useful for this. Is there a numerical trick I can use to correct these jumps? Perhaps if I specified a constant pressure gradient instead of the iterative one?

Also, would you be able to provide a reference for the issues with non-orthogonality in PBC so I can read about it further? Or is it more of a speculative issue?

Thanks in advance,
Charlie

[Santiago]

Ummm... On top of my head I have no specific reference in which they could mention this issue, in general good unstructured meshing is more of an "art". Anyway, I could recommend you to read "Finite-Volume CFD procedure and Adaptive Error Control Strategy for grids of arbitrary topology" of Muzaferija & Gosman as a starting point, then you could delve a bit on the references therein mentioned. As for how to "smooth-out" those, you could either do a filtering of the resulting field (low-pass filter on post-processing) or when running the simulations using dissipative (upwind) schemes for the convective terms; in the end it all depend on the degree of accuracy you wish to get, both options will "regularize" the solution in a way that you'll end up seeing the solution of a slightly different NS equation.

[C-L]

Thanks for the reference, I will be sure to go through it. I think adding dissipation into the numerics would reduce the accurately too much for this work so I think I will have to investigate the meshing further.

Would another solution be defining a constant body force to the Ueq in order to fix the mean pressure gradient? This would mean there would be no iterative forcing on the boundary so would this remove the boundary errors?

[Santiago]

Mmmm, I honestly don't think that's the problem. The calculation of the driving force given a bulk flux/velocity returns a scalar, not a field. If you refer to what adjustPhi does, well, mass conservation on the boundaries must be satisfied at all times at each iteration after finishing the pressure correction step of PISO; but this function is only called when pressure on a given boundary are defined as zeroGradient but in your case those BC's should be 'cyclic'.

[C-L]

I have decided to run the problem with a constant forcing term just for a sanity check but I expect the continuity problem will still persist for the reasons you mentioned. I think it is time to investigate snappyHexMesh as an alternative for complex geometries!

Thanks for the advice.

[Santiago]

Originally Posted by C-L

I have decided to run the problem with a constant forcing term just for a sanity check but I expect the continuity problem will still persist for the reasons you mentioned. I think it is time to investigate snappyHexMesh as an alternative for complex geometries!

Thanks for the advice.

I would instead try to generate prisms on the boundaries (extrude the boundary faces) and then let the internal domain to be whatever you had at the beginning. Using snappy will make no difference, stick to whatever meshing software you know best.

[C-L]

As you suggested I have extruded the faces on both boundaries in order to align them normal to the flow direction. However, after running the code again I am still getting spikes on the boundary (see the image below).

The only issue that I can see from my case files are that I have not used the 'preservePatches' function but I have not used it for any of the other cases which have all worked well. Is there anything obviously wrong with the files below?

I have also attached a few timesteps from the pimpleFoam log - the only issue I can see is that the bounding flm and fmm values for the LES model are very large, corresponding to the spikes on the boundary.

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     pimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         15000;

deltaT          0.00001;

maxCo		1;

adjustTimeStep	yes;

writeControl    timeStep;

writeInterval   500;

purgeWrite      0;

writeFormat     ascii;

writePrecision  8;

writeCompression on;

timeFormat      general;

timePrecision   8;

runTimeModifiable true;


// ************************************************************************* //

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         backward;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss filteredLinear;
    div(phi,flm)    Gauss filteredLinear;
    div(phi,fmm)    Gauss filteredLinear;
    div(phi,nuTilda) Gauss filteredLinear;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method meshWave;
}


// ************************************************************************* //

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          GAMG;
        tolerance       0;
        relTol          0.05;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    pFinal
    {
        $p;
        smoother        DICGaussSeidel;
        tolerance       1e-08;
	relTol		0;
    }

    "(U|k|nuTilda|flm|fmm)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-08;
    }

    "(U|k|nuTilda|flm|fmm)Final"
    {
        $U;
    }
}

PIMPLE
{
    nOuterCorrectors 1;
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefPoint	    (0 0 0);
    pRefValue       0;
}


// ************************************************************************* //

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    note        "mesh decomposition control dictionary";
    object      decomposeParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

numberOfSubdomains  64;

//- Keep owner and neighbour on same processor for faces in zones:
// preserveFaceZones (heater solid1 solid3);

//- Keep owner and neighbour on same processor for faces in patches:
//  (makes sense only for cyclic patches)
//preservePatches (cyclic_half0 cyclic_half1);

//- Keep all of faceSet on a single processor. This puts all cells
//  connected with a point, edge or face on the same processor.
//  (just having face connected cells might not guarantee a balanced
//  decomposition)
// The processor can be -1 (the decompositionMethod chooses the processor
// for a good load balance) or explicitly provided (upsets balance).
//singleProcessorFaceSets ((f0 -1));


//- Use the volScalarField named here as a weight for each cell in the
//  decomposition.  For example, use a particle population field to decompose
//  for a balanced number of particles in a lagrangian simulation.
// weightField dsmcRhoNMean;

// method          scotch;
// method          hierarchical;
method          simple;
// method          metis;
// method          manual;
// method          multiLevel;
// method          structured;  // does 2D decomposition of structured mesh

multiLevelCoeffs
{
    // Decomposition methods to apply in turn. This is like hierarchical but
    // fully general - every method can be used at every level.

    level0
    {
        numberOfSubdomains  64;
        //method simple;
        //simpleCoeffs
        //{
        //    n           (2 1 1);
        //    delta       0.001;
        //}
        method scotch;
    }
    level1
    {
        numberOfSubdomains  4;
        method scotch;
    }
}

// Desired output

simpleCoeffs
{
    n           (4 4 4);
    delta       0.001;
}

hierarchicalCoeffs
{
    n           (1 2 1);
    delta       0.001;
    order       xyz;
}

metisCoeffs
{
 /*
    processorWeights
    (
        1
        1
        1
        1
    );
  */
}

scotchCoeffs
{
    //processorWeights
    //(
    //    1
    //    1
    //    1
    //    1
    //);
    //writeGraph  true;
    //strategy "b";
}

manualCoeffs
{
    dataFile    "decompositionData";
}

structuredCoeffs
{
    // Patches to do 2D decomposition on. Structured mesh only; cells have
    // to be in 'columns' on top of patches.
    patches     (bottomPatch);
}

//// Is the case distributed? Note: command-line argument -roots takes
//// precedence
//distributed     yes;
//// Per slave (so nProcs-1 entries) the directory above the case.
//roots
//(
//    "/tmp"
//    "/tmp"
//);

// ************************************************************************* //

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

simulationType LES;

LES
{
    LESModel        dynamicLagrangian;

    turbulence      on;

    printCoeffs     on;

    delta           vanDriest;
    
    dynamicLagrangianCoeffs
	{
	filter       simple;
	ce         1.048;
	theta     1.5;
	}


    cubeRootVolCoeffs
    {
        deltaCoeff      1;
    }

    PrandtlCoeffs
    {
        delta           cubeRootVol;
        cubeRootVolCoeffs
        {
            deltaCoeff      1;
        }

        smoothCoeffs
        {
            delta           cubeRootVol;
            cubeRootVolCoeffs
            {
                deltaCoeff      1;
            }

            maxDeltaRatio   1.1;
        }

        Cdelta          0.158;
    }

    vanDriestCoeffs
    {
        delta           cubeRootVol;
        cubeRootVolCoeffs
        {
            deltaCoeff      1;
        }

        smoothCoeffs
        {
            delta           cubeRootVol;
            cubeRootVolCoeffs
            {
                deltaCoeff      1;
            }

            maxDeltaRatio   1.1;
        }

        Aplus           26;
        Cdelta          0.158;
    }

    smoothCoeffs
    {
        delta           cubeRootVol;
        cubeRootVolCoeffs
        {
            deltaCoeff      1;
        }

        maxDeltaRatio   1.1;
    }
}


// ************************************************************************* //

Courant Number mean: 0.045302183 max: 1.0263954
deltaT = 8.5202102e-05
Time = 4.069965095738113

PIMPLE: iteration 1
smoothSolver:  Solving for Ux, Initial residual = 0.0032301636, Final residual = 7.6381712e-09, No Iterations 6
smoothSolver:  Solving for Uy, Initial residual = 0.0034618555, Final residual = 9.157908e-09, No Iterations 10
smoothSolver:  Solving for Uz, Initial residual = 0.00032816168, Final residual = 7.7291454e-09, No Iterations 4
Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92385471
GAMG:  Solving for p, Initial residual = 0.040015976, Final residual = 0.00082851291, No Iterations 2
time step continuity errors : sum local = 4.160954e-09, global = -6.2520222e-10, cumulative = -7.3492066e-05
Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92407577
GAMG:  Solving for p, Initial residual = 0.0031098961, Final residual = 7.7192676e-09, No Iterations 23
time step continuity errors : sum local = 6.2523534e-10, global = -6.2520222e-10, cumulative = -7.3492691e-05
Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92426602
smoothSolver:  Solving for flm, Initial residual = 0.029661143, Final residual = 5.2469417e-09, No Iterations 5
bounding flm, min: -667100.81 max: 2002938.7 average: 9.1001611
smoothSolver:  Solving for fmm, Initial residual = 0.013253498, Final residual = 3.8549385e-09, No Iterations 5
bounding fmm, min: -34953685 max: 2.3359435e+09 average: 13780.572
ExecutionTime = 172032.13 s  ClockTime = 172379 s

Courant Number mean: 0.044137162 max: 1.0137825
deltaT = 8.4043768e-05
Time = 4.070049139505807



PIMPLE: iteration 1
smoothSolver:  Solving for Ux, Initial residual = 0.0030180317, Final residual = 7.48811e-09, No Iterations 3
smoothSolver:  Solving for Uy, Initial residual = 0.0032379585, Final residual = 9.068158e-09, No Iterations 3
smoothSolver:  Solving for Uz, Initial residual = 0.00030659573, Final residual = 5.8584048e-09, No Iterations 2
Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.9253623
GAMG:  Solving for p, Initial residual = 0.036520884, Final residual = 0.00062184272, No Iterations 2
time step continuity errors : sum local = 2.953574e-09, global = -5.9219333e-10, cumulative = -7.3506524e-05
Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92528735
GAMG:  Solving for p, Initial residual = 0.0026411813, Final residual = 8.5933446e-09, No Iterations 25
time step continuity errors : sum local = 5.9222588e-10, global = -5.9219333e-10, cumulative = -7.3507117e-05
Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92539805
smoothSolver:  Solving for flm, Initial residual = 0.041895575, Final residual = 9.7797156e-10, No Iterations 6
bounding flm, min: -848496.88 max: 904626.11 average: 7.4115014
smoothSolver:  Solving for fmm, Initial residual = 0.012833257, Final residual = 2.5986852e-09, No Iterations 5
bounding fmm, min: -29135461 max: 1.2818239e+09 average: 12415.778
ExecutionTime = 172069.39 s  ClockTime = 172416 s

Courant Number mean: 0.041233053 max: 0.98170613
deltaT = 8.1079342e-05
Time = 4.071023715504507

PIMPLE: iteration 1
smoothSolver:  Solving for Ux, Initial residual = 0.0030767242, Final residual = 9.0609907e-09, No Iterations 3
smoothSolver:  Solving for Uy, Initial residual = 0.0032985252, Final residual = 5.9680932e-10, No Iterations 4
smoothSolver:  Solving for Uz, Initial residual = 0.00031225705, Final residual = 6.2052363e-09, No Iterations 2
Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92538027
GAMG:  Solving for p, Initial residual = 0.040505436, Final residual = 0.00061931897, No Iterations 2
time step continuity errors : sum local = 3.0254994e-09, global = -6.0425345e-10, cumulative = -7.3507721e-05
Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92538532
GAMG:  Solving for p, Initial residual = 0.0029228641, Final residual = 9.70328e-09, No Iterations 23
time step continuity errors : sum local = 6.0429122e-10, global = -6.0425345e-10, cumulative = -7.3508325e-05
Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92541771
smoothSolver:  Solving for flm, Initial residual = 0.041116896, Final residual = 1.3082178e-09, No Iterations 6
bounding flm, min: -981020.44 max: 809933.1 average: 7.3560911
smoothSolver:  Solving for fmm, Initial residual = 0.01373234, Final residual = 3.4180578e-09, No Iterations 5
bounding fmm, min: -52698474 max: 1.3968272e+09 average: 12280.821
ExecutionTime = 172072.28 s  ClockTime = 172419 s

Courant Number mean: 0.042001431 max: 0.97433007
deltaT = 8.3215477e-05
Time = 4.071106930981983