[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. The solver works fine with second-order spatial schemes and first-order (Euler) time advance, but when I try to switch to second-order (backward) time advance, the pressure field develops non-physical "beads" of high and low pressure. If I leave the simulation for long enough, many other non-physical things happen to the pressure field, and the simulation often diverges or stalls (delta t goes to zero).

The evolution of the p field is illustrated here: http://imgur.com/a/fcF2I [note: the colour scale varies between images].

The test case: a turbulent jet with laminar coflow

The boundary conditions:

• Lower plane: turbulentInlet for U, zeroGradient for p
• Outer radius: inlet which may also behave as an outlet if flow conditions require (inletOutlet for U, outletInlet for p; expect flow predominantly in to domain)
• Upper plane: outlet which may also behave as an inlet if flow conditions require (inletOutlet for U, outletInlet for p; expect flow predominantly out of domain)

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.

system/fvSchemes:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.1.0                                 |
|   \\  /    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
{
    div(phi,rho)    Gauss limitedLinear phi 1;
    default         Gauss limitedLinear 1;
    div(phi,U)      Gauss filteredLinear2V 0.2 0;
    div(phi,Yi_h)   Gauss multivariateSelection 
    {
        CO2               limitedLinear01 1; 
        H2O               limitedLinear01 1;
        NO                limitedLinear01 1;
        N2                limitedLinear01 1;
        h                 filteredLinear2 0.2 0;
    };
    div((muEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p;
}


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

system/fvSolution:

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

solvers
{
    "(p|rho)"
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-08;
        relTol       0;
        // relTol          0.01;
    }

    "(p|rho)Final"
    {
        $p;
        relTol          0;
    }

    "(U|h|Z|varZ|Yi|k|epsilon)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-06;
        relTol          0;
        // relTol          0.01;
    }

    "(U|h|Z|varZ|Yi|k|epsilon)Final"
    {
        $U;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor yes; 
    nOuterCorrectors  50;
    nCorrectors     3;
    nNonOrthogonalCorrectors 1;
    turbOnFinalIterOnly false;

    residualControl
    {
        U
        {
            tolerance 1e-4;
            relTol       0;
        }
        p
        {
            tolerance 1e-4;
            relTol       0;
        }
    }
}

relaxationFactors
{
    fields
    {
        p    0.3;
        pFinal 1;
    }
    equations
    {
        "(U|k|epsilon)"                  0.7; 
        "(U|k|epsilon)Final"               1; 
    }
}


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

As suggested at https://openfoamwiki.net/index.php/O...hm_in_OpenFOAM, I have used residualControl to determine when the outer PIMPLE loops are finished;with these settings, it seems to take 5-6 iterations.

My question: Why is the p field becoming non-physical, and what can I do to fix it while maintaining second-order accuracy? I have tried tightening the solver tolerances and adding adding more inner and non-orthogonal correctors, but they do not seem to help. I was thinking that I should add more PIMPLE iterations, but 5-6 already seems like a lot to me. Should I be reducing the Courant number? I would be happy to do this for a transition period, but if I have to keep Co very low for the remainder of my simulation then I think that this could become prohibitively expensive.

[Artur]

Hi,

Have you tried the Crank-Nicolson time scheme? It also gives you an option to blend in a bit of first-order time which could help with your issue. Also, have you tried running the solver with a single outer loop, i.e. in PISO mode, or is that not an option for your case? I found that the latter sometimes helps with limiting problems with convergence and reduces run times, although may require you to reduce the time step a bit to keep the residuals under control. Furthermore, I take it you are using adjustable time step? If so, maybe try switching it off and fixing delta t to something which you already know will keep the Co<0.5 or whatever criterion you are following. I've found the adjustable time step with backward scheme to yield very similar looking issues with one of the old OF versions when using a multi-phase solver.

All the best,

A

[knuckles]

Hi Artur,

Thanks for the reply. These are all new ideas to me, so I will try them out!

Why would you expect that switching from PIMPLE to PISO would be more stable? Is the idea to switch from (PIMPLE with bigger time step) to (PISO with smaller time step), with the expectation that the two changes will kind of cancel out, or maybe even lead to improvements in both simulation time and stability?

[Artur]

Not sure exactly what your solver is doing, but for some of my cases with cavitation I've found that big time steps with pimple may cause trouble in other transport equations, which then pollutes your solution and eventually causes divergence. You said you switched all the extra pieces off but just maybe there is still something there which is causing issues, hence my idea. Could be entirely wrong though

[knuckles]

Hi Artur,

I tried some of your suggestions. Here are my results so far:

• Keeping PIMPLE with delta t = 2e-6 (which should give a Co around 0.55) gives a somewhat similar pressure field. PIMPLE typically loops 6 times to respect the error tolerances I've set
• Switching PISO with a delta t of 5e-7 (which should give a Co around 0.14) causes instabilities large enough to fully halt the solver with a nan
• Switching to PISO with a delta t of 2.5e-7 (which should give a Co around 0.07) causes the emergence of similar pressure waves, but with a much larger amplitude

My thinking is that 6 PIMPLE loops should have a similar computational complexity to 6 PISO time steps, and therefore the first case and the last (with 1/8th the time step) are approximately comparable in terms of cost. Given that the first result seems less incorrect than the last, I think my most promising route forward is to try PIMPLE again with a smaller, fixed delta t. Do you have any other suggestions?

[Artur]

I agree with you in terms of rough cost equivalence of the two algorithms. It's probably a good idea to try PIMPLE with reduced time step just to get a full picture of how the problem behaves in different configurations. If it still persists there are always the standard questions of BCs and whether the mesh is sufficiently dense. Perhaps it'd also be worth to run the case with a built-in pimple/piso solver using the same grid and setup just to verify there isn't anything in the custom one you're using that's causing the issue?

[knuckles]

Hi Artur,

My newest results:

• PIMPLE with delta t = 1e-6 (which should give a max Co around 0.27, and actually gives a max Co around 0.19) has the same p problem
• Switching to reactingFoam [the most similar base solver] with PIMPLE and delta t = 1e-6 has the same p problem

The discrepancy between the max Co ("should give" vs "actually gives") is because my "should give" calculations are based on a back-of-the-envelope analysis and OpenFoam is using the actual instantaneous data.


I can keep reducing the time step, but I feel like Co = 0.19 is all ready pretty low, no?

[Artur]

Yes it is very low from a flow point of view, although from reaction point of view I have no idea to be honest. I suppose this brings the checklist to the BCs and mesh... After that I'm really out of ideas since your schemes and solvers look OK to me (I normally wouldn't use smoothSolver myself but that should give these strange results I reckon). Did you base your BCs on an existing case or a tutorial?

[knuckles]

My BCs are based on the "pilotedDiffusionFlame" tutorial case from the "flameletfoam" solver: https://openfoamwiki.net/index.php/E...n/flameletFoam. The geometry in the flameletfoam test case is quite similar to mine, but they run RANS (and thus use a 2D domain with wedge boundaries) whereas I run LES (and thus use a 3D domain with those boundaries removed). Their tutorial files are on github at https://github.com/flameletFoam/flam...ionFlame/ras/0. For reference, my jet Reynolds number (rho * v_bulk * d_jet / mu) is 22,400, which I believe is much larger than the value for the flameletfoam test case.

I made the mesh myself. I don't consider myself a meshing expert, but I assume that it isn't horribly broken given the fact that the checkMesh utility has no complaints about it. I can post the blockMeshDict if that might help.

[Artur]

The mesh is probably fine in terms of quality, I'm no expert in non-wall-bounded flows so not sure what kind of resolution you need to resolve turbulent stresses to a level sufficient for the sgs model, so it's probably worth checking that with the literature unless you're confident in the values you've chosen.
On thing which caught my attention is the proximity of the outlet and the fact that you (the tutorial) use inletOutlet on U and fixedValue on p while from your images one can tell there's quite a lot of turbulence reaching the outlet. I realise this is quite expensive but maybe moving the outlet further and switching the U BC to pressureInletOutletVelocity would help? I found it to work very well with fV on p. Optionally, I think you might also look into using a convective outlet BC on U which afaik is supposed to be design for the exact kind of situation like what you're having now (turbulence reaching edge of the domain), there's some code and more discussion here:
http://www.cfd-online.com/Forums/ope...ady-flows.html
I'm sorry but if none of these help then I'm officially out of ideas

A

[Alish1984]

Originally Posted by knuckles

My BCs are based on the "pilotedDiffusionFlame" tutorial case from the "flameletfoam" solver: https://openfoamwiki.net/index.php/E...n/flameletFoam. The geometry in the flameletfoam test case is quite similar to mine, but they run RANS (and thus use a 2D domain with wedge boundaries) whereas I run LES (and thus use a 3D domain with those boundaries removed). Their tutorial files are on github at https://github.com/flameletFoam/flam...ionFlame/ras/0. For reference, my jet Reynolds number (rho * v_bulk * d_jet / mu) is 22,400, which I believe is much larger than the value for the flameletfoam test case.

I made the mesh myself. I don't consider myself a meshing expert, but I assume that it isn't horribly broken given the fact that the checkMesh utility has no complaints about it. I can post the blockMeshDict if that might help.

Hi,

What is ur diameter, velocity and level of turbulence intensity of the jet?

[Hen Cruise]

Hi knuckles,
did you solve the non-physical p fields problem? i recently meet a similar pressure issue when i use the 2nd time scheme in reactingfoam to calculate jet flame, or can you please give me some advice?