[BendikS]

Hi,
I'm running a case in OpenFOAM (v1606+) where I have trouble with solution convergence and bounding of k and omega.

The case is a oscillating cover (looks like a 2D cut of a bowler hat) in infinite fluid. There is no inflow or outflow, just still water and the cover oscillating in the y-direction. The simulation is obviously transient.

For the mesh motion I'm using dynamicMotionSolverFvMesh with the displacementLaplacian solver and InverseDistance diffusivity model.

For turbulence I'm using the k-omega SST model, without wall functions.

For boundary conditions I have zeroGradient for U and P on all the domain sides (except for front and back which are empty). On the cover I have zeroGradient for P and movingWallVelocity with uniform value (0,0,0) for U.

k and omega also have zeroGradient on all sides, but with fixedValue uniform 1e-11; on the cover. nut have the same BC on the cover, while the other fields are set to calculated.

I have been trying to have second order schemes, but with central schemes I get bounding problems with k and omega and the solution crashes.

I just got the solution to run by switching to the following scheme:

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         cellMDLimited Gauss linear 0.5;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss inearUpwind cellMDLimited Gauss linear 0.5;

    div(phi,k)      Gauss linearUpwind cellMDLimited Gauss linear 0.5;
    div(phi,omega)  Gauss linearUpwind cellMDLimited Gauss linear 0.5;

    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default        Gauss linear limited corrected 0.5;

    laplacian(diffusivity,cellMotionU) Gauss linear limited corrected 0.5;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited corrected 0.5;
}
wallDist
{
    method meshWave;
}
fluxRequired
{
    default         no;
    pcorr           ;
    p               ;
}

My solution dict:

solvers
{
    pcorr
    {
        solver              GAMG;
        tolerance           0.02;
        relTol              0;
        smoother            GaussSeidel;
        nPreSweeps          0;
        nPostSweeps         2;
        cacheAgglomeration     no;
        nCellsInCoarsestLevel     10;
        agglomerator        faceAreaPair;
        mergeLevels         1;
    }

    p
    {
        $pcorr;
        tolerance           1e-7;
        relTol              0.01;
    }

    pFinal
    {
        $p;
        tolerance           1e-07;
        relTol              0;
    }

    U
    {
        solver              PBiCG;
    preconditioner        DILU;
        tolerance         1e-08;
        relTol              0;//0.1;
    minIter            3;
    }

    UFinal
    {
        $U;
        //tolerance           1e-05;
        //relTol              0;
    }

    "cellDisplacement.*"
    {
        solver              GAMG;
        tolerance           1e-8;
        relTol              0;
        smoother            GaussSeidel;
        cacheAgglomeration     true;
        nCellsInCoarsestLevel     10;
        agglomerator        faceAreaPair;
        mergeLevels         1;
    }

    "(k|omega)"
    {
    solver            PBiCG;
    preconditioner        DILU;
    tolerance         1e-08;
    reltol             0;
    minIter            3;
    }
    "(k|omega)Final"
    {
    solver            PBiCG;
    preconditioner        DILU;
    tolerance         1e-08;
    reltol             0;
    minIter            3;
    }
}

PIMPLE
{
    correctPhi              yes;
    nOuterCorrectors        2; 
    nCorrectors             1; 
    nNonOrthogonalCorrectors     1; 
    pRefCell             1001; 
    pRefValue             0; 
}

relaxationFactors
{
    fields
    {
    }
    equations
    {
        "U.*"           1;
       "k.*"           1;
    “omega.*”         1;
    }
}

But I still have issues with bounding values of k and omega. I have dt = 1e-4 which results in max Courant numbers of around 0.4. The mesh is generated with SHM and the checkMesh results are listed below:

Time = 0

Mesh stats
    points:           64712
    internal points:  0
    faces:            127108
    internal faces:   63212
    cells:            31584
    faces per cell:   6.02584
    boundary patches: 6
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     30474
    prisms:        202
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     908
    Breakdown of polyhedra by number of faces:
        faces   number of cells
            7   798
            8   110

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
                   inlet       40       82  ok (non-closed singly connected)
                  outlet       40       82  ok (non-closed singly connected)
            topAndBottom       80      164  ok (non-closed singly connected)
                   front    31584    32356  ok (non-closed singly connected)
                    back    31584    32356  ok (non-closed singly connected)
                   cover      568     1136  ok (non-closed singly connected)

Checking geometry...
    Overall domain bounding box (-0.625 -0.625 0.00625) (0.625 0.625 0.02125)
    Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
    Mesh has 2 solution (non-empty) directions (1 1 0)
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (-7.57749e-19 9.55767e-20 7.24336e-15) OK.
    Max cell openness = 4.09121e-16 OK.
    Max aspect ratio = 24.2247 OK.
    Minimum face area = 4.86616e-08. Maximum face area = 0.000985343.  Face area magnitudes OK.
    Min volume = 7.29923e-10. Max volume = 1.47801e-05.  Total volume = 0.0234301.  Cell volumes OK.
    Mesh non-orthogonality Max: 35.522 average: 3.79143
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 1.90833 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

Any suggestions?

Thanks!

[matejfor]

Hi,

not still sure about the BCs although I read the description twice, but anyways...

a] I am not very sure about the simulation with k-omega at the walls without wall functions employed. I am afraid you simply need to have a wall function involved by the definition of the model. Please refer to the list of functions here: http://openfoam.com/documentation/cp...Functions.html

b] Second order - the linearUpwind is actually second order scheme, would you like to use more precise schemes? Your mesh is actually of a great quality so I would not worry much about the limitation of non-orthogonality. Not needed for your case.
I wonder you'd be able to go to more precise schemes on turbulence in OpenFOAM.
Actually what would happen if your omega would be upwinded only? Would the bounding problem remain?

c] Instead of PBiCG you may try PBiCGStab - it is a new stabilised version

d] general though - you are going to have a problem with k and omega if your U gradients are unphysical. Therefore your velocity field should be well converged for each time step. Is it the case?

e] Maybe others will spot different problem.

[BendikS]

Thanks for the answer Matej.

Sorry for being a bit unclear on the BCs, here are my BCs for k, nut, omega, p and U.

edit: note that the left and right sides of the domain is called inlet and outlet even though there's now flow through them.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 1e-11;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    topAndBottom
    {
        type            zeroGradient;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    cover
    {
        type            fixedValue;
        value           uniform 1e-11;
    }
}

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 1e-11;

boundaryField
{
    inlet
    {
        type            calculated;
        value           uniform 1e-11;
    }
    outlet
    {
        type            calculated;
        value           uniform 1e-11;
    }
    topAndBottom
    {
        type            calculated;
        value           uniform 1e-11;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    cover
    {
        type            fixedValue;
        value           uniform 1e-11;
    }
}

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      omega;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 -1 0 0 0 0];

internalField   uniform 1;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    topAndBottom
    {
        type            zeroGradient;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    cover
    {
        type            fixedValue;
        value           uniform 1e-11;
    }
}

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    topAndBottom
    {
        type            zeroGradient;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    cover
    {
        type            zeroGradient;
    }
}

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "0";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    topAndBottom
    {
        type            zeroGradient;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    cover
    {
        type            movingWallVelocity;
        value           uniform (0 0 0);
    }
}

a) Im not using any wall functions as described here http://www.wolfdynamics.com/images/O...ng/module8.pdf on slide 34. I just switched from using wall functions, described in the slides above slide 34, but I had the same problem there aswell.

b) I do know that linearUpwind is second order, so Im happy to settle with it. Ideally I would prefer central schemes, but I have not managed to get it to run, so the next step was to switch to linearUpwind and check if that helped. It did, as the simulations now are running without crashing, but I still have the boundedness problems as mentioned. I have not tried to use first order upwind schemes for omega, but if I nothing else helps I guess that's the only solution. Why just for omega btw? and not for both omega and k?

c) I'll try PBiCGStab instead, thanks.

d) I have no problems with convergence in U. The initial residual is reasonably low, and it drops down to ~10^-15 in a couple of iterations.

[matejfor]

Hi,
one more thing, I do not know where did you got these lines:
div(phi,U) Gauss inearUpwind cellMDLimited Gauss linear 0.5;
div(phi,k) Gauss linearUpwind cellMDLimited Gauss linear 0.5;
div(phi,omega) Gauss linearUpwind cellMDLimited Gauss linear 0.5;

but it does not make much sense any more.

The current style is:
div(phi,U) Gauss linearUpwind grad(U);
div(phi,k) Gauss linearUpwind grad(k);
div(phi,omega) Gauss uwind;

If you wish to use purely linear scheme, as it is the holy grail of second order, you are always going to have problems with boundedness, as the scheme is bounded only in 1D situation. So if you manage to keep your Co<0.1 approx and Pe<=2 you should be OK. But it is hardly ever achievable. Even for LES we use LUST, which blends linear with linearUpwind.

Second thing: if you insist on accuracy of your spacial discretisation you should be demanding in the same way on your time discretisation. Using Euler (first order - very dispersive on your oscilations), you should use Crank-Nicolson scheme.

[BendikS]

Hi,

I believed those two statements were the same, but in the latest simulations I have been using the "good old" version.

I switched to Upwind for omega now, and the bounding issue for omega stopped, but k is still having some "oscillatory" bounding issues, i.e. its fine for maybe ten iterations, but then k is negative for one iterations or two, before it's behaving well again.

Regarding central schemes, I guess it will be to hard to get stable solutions with it.

yup, switching to crank nicolson aswell, just wanted to figure out the spatial schemes first.

[BendikS]

Ok so I switched to
ddt: CrankNicolson 0.9;
grad: default cellMDLimited Gauss Upwind 0.5;
div(U): linearUpwind grad(U);
div(k): upwind;
div(omega): upwind;
and gauss linear corrected for laplacian schemes.

The simulations runs, and I'm not having the bounding problems anymore, but suddenly around the t = 8s, the pressure residual increases and behaves a bit strangely considering that the cover is just continuing to oscillate. There is nothing "new" happening in the simulations physically.

I attached the residual plot.

Any thoughts?

Thanks

[minh khang]

I got the same issue.

[Artur.Ant]

From my experience:

sometimes chosing the second order scheme for pressure can give some issued with k an omega too (even if upwind is used for turbulence).
I tried a 3D case compressible, supersonic with shockWaves within the domain:
using MUSCL for div(phid, p) in some cases, not always, I got some troubles with bounding k and omega. In some other cases also p was bounded in small regions of the domain.

changing the scheme for div(phid, p) to LUST everything stabilised very quickly and work smoothly.

I also got the impression that the choice of the numerical scheme is depend on the type of processor. Somebody can confirm it?