[knuckles]

I'm running LES of a turbulent jet using ~1.6 million cells and a Courant Number of 0.3. My solver is customized, but all of the non-standard features are turned off and it should be behaving identically to reactingFoam with all reactions turned off. At some point, an instability arises: at first, I get a few warnings that T is outside of the allowable range (300-5000 K) - say 5100 K - but then these warnings get more frequent and more problematic, until there are thousands at each time step and the associated temperatures range from -1e16 to +1e16.

My best guess for what's going on is that I've failed to appropriately model the pressure-momentum coupling, and that I need to add more pimple loops - is this the correct fix? I have seen an article [https://openfoamwiki.net/index.php/O...hm_in_OpenFOAM] which suggests, among other things:

• Set nOuterCorrectors up to a very high number (~ 50 to 200) and control your pimple loop with the residual control
• Set nCorrectors to a small number (1 to 3; you have to play around)

...however, I've never seen this first suggestion implemented in a tutorial other than the one described at the link above.

My question: What settings should I be changing to avoid this instability?

For reference:

The output of checkMesh:

Mesh stats
    points:           1617892
    faces:            4824873
    internal faces:   4797927
    cells:            1603800
    faces per cell:   6
    boundary patches: 5
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     1603800
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology                  
    JET                 405      424      ok (non-closed singly connected)  
    PILOT               864      900      ok (non-closed singly connected)  
    COFLOW              1404     1440     ok (non-closed singly connected)  
    OUTERRADIUS         21600    21636    ok (non-closed singly connected)  
    OUTLET              2673     2692     ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-0.07172601826 -0.07172601826 0) (0.07172601826 0.07172601826 0.36)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (2.5360889850782e-16 2.502408869494276e-16 3.295433281148142e-16) OK.
    Max cell openness = 4.165398281849697e-16 OK.
    Max aspect ratio = 22.5683629471026 OK.
    Minimum face area = 5.055728864375342e-08. Maximum face area = 3.708057695313814e-05.  Face area magnitudes OK.
    Min volume = 1.872805105800212e-11. Max volume = 4.104974063197774e-08.  Total volume = 0.003402726340535457.  Cell volumes OK.
    Mesh non-orthogonality Max: 22.23980401213051 average: 2.729797707905458
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.5165905340833714 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

My system/fvSchemes

ddtSchemes
{
    default         backward;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    div(phi,rho)    Gauss limitedLinear phi 1;
    default         Gauss limitedLinear 1;
    div(phi,U)      Gauss limitedLinearV 1;
    div(phi,Yi_h)   Gauss multivariateSelection 
    {
        CO2               limitedLinear01 1; 
        H2O               limitedLinear01 1;
        NO                limitedLinear01 1;
        N2                limitedLinear01 1;
        h                 limitedLinear   1;
    };
    div((muEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p;
}

My system/fvsolution:

"(p|rho)Final" 
    { 
        $p; 
        relTol          0; 
    } 
 
    "(U|h|Z|varZ|Yi|k|epsilon)" 
    { 
        solver          smoothSolver; 
        smoother        symGaussSeidel; 
        tolerance       1e-06; 
        relTol          0.01; 
    } 
 
    "(U|h|Z|varZ|Yi|k|epsilon)Final" 
    { 
        $U; 
        relTol          0; 
    } 
} 
 
PIMPLE 
{ 
    momentumPredictor yes; 
    nOuterCorrectors  2; 
    nCorrectors     3; 
    nNonOrthogonalCorrectors 2; 
    turbOnFinalIterOnly false; 
} 
 
relaxationFactors 
{ 
    fields 
    { 
    } 
    equations 
    { 
        ".*"    1; 
    } 
}

[knuckles]

I believe I found the problem. I had a turbulentInlet BC and the magnitude of the fluctuations was far too large - presumably this caused this instability. For reference, a fluctuationScale of 0.1 or larger is very big.