[cyln]

Hello,

I am running a transient simulation using pimpleFoam. I have been trying to optmize the convergence rate for a while, however, my residuals are not converging. The photos of residuals from two different runs are attached:
1. 1 time step, 600 pimple iterations
2. 5 time steps, 200 pimple iterations

I would be grateful if you could suggest me anything. My system files are as follows:

My controlDict file:

application     pimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         0.00015;

deltaT          0.00015;

writeControl    timeStep;

writeInterval   1;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable yes;

My fvSchemes:

ddtSchemes
{
    default	    CrankNicolson 0.9;
}

gradSchemes
{
    default         Gauss linear;
     grad(p)         cellLimited Gauss linear 1;
     grad(U)         cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;

	default          		none;
	div(phi,U)       		 Gauss LUST grad(U);
	div(phi,k)       		 Gauss linear ;
	div(phi,omega)   		 Gauss linear ;
	div(phi,R)	 		Gauss linear ;
       div((nuEff*dev2(T(grad(U))))) Gauss linear;   //version 4.1
}

laplacianSchemes
{
    //default         Gauss linear corrected 0.33;
    default         Gauss linear limited corrected 0.33;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited corrected 0.33;
}

fluxRequired
{
    default         no;
    p;
}

wallDist
{
    method meshWave;
}

My fvSolution file:

solvers
{
    p
    {

        solver           GAMG;
        tolerance        0;
        relTol           1e-2;
        smoother         GaussSeidel;
        nPreSweeps       0;
        nPostSweeps      3;
        cacheAgglomeration on;
        agglomerator     faceAreaPair;
        nCellsInCoarsestLevel 50000;
        mergeLevels      1;

    }


    pFinal
    {
        solver           GAMG;
        tolerance        1e-4;
        relTol           0;
        smoother         GaussSeidel;
        nPreSweeps       0;
        nPostSweeps      3;
        cacheAgglomeration on;
        agglomerator     faceAreaPair;
        nCellsInCoarsestLevel 50000;
        mergeLevels      1;
    }

    "(U|k|omega|nut)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance        0;
        relTol           1e-3; //
        nSweeps          4;

    }

    "(U|k|omega|nut)Final"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance        1e-3; //
        relTol           0;
        nSweeps          4;

    }

}



PIMPLE
{
    nNonOrthogonalCorrectors 0;
    nCorrectors		     2;
    nOuterCorrectors        600;
    pRefCell                 0;
    pRefValue                0;
	
	 residualControl
    {
        "(p|U)"
        {
            tolerance 1e-6;
            relTol 0;
        }
    }
}

relaxationFactors
{
    fields
    {
        p               0.3;//0.3
	pFinal          1.0;
    }
    equations
    {
        U               0.7;
        k               0.7;
        omega           0.7;
	UFinal		1.0;
	omegaFinal	1.0;
	kFinal		1.0;
    }
}

[matejfor]

Hi.

a) your endTime = timeStep or I need glasses.

b) what is your mesh quality? You are trying to be as 2nd order as possible in discretisation, but still you are limiting the diffusion flux, correction is not enough?

c) you have two PISO loops and 600 outer loops? and you have a single time step. and you have problems to converge within this time step.... and you are underrelaxing ... your case is not that transient as I expected.

Can you explain why you have such a setup? What is the application and what are you trying to achieve?

[cyln]

Hi matejfor, thanks for your reply.

a. The simulation runs itself when I am away. So I check the convergence at one time step only because the residuals do not fall below that in following time steps. Once I get the first time step converged, I will fix the endTime.

b. My mesh quality is below. It is not because I want to limit the diffusion flux. What would you recommend?

Mesh stats
    points:           2341634
    faces:            20711284
    internal faces:   20259420
    cells:            9802885
    faces per cell:   4.17945
    boundary patches: 10
    point zones:      0
    face zones:       2
    cell zones:       2

Overall number of cells of each type:
    hexahedra:     0
    prisms:        1759164
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    8043721
    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                  
    INLET1              3389     2815     ok (non-closed singly connected)  
    INLET2              7406     4610     ok (non-closed singly connected)  
    CYL                 25066    12744    ok (non-closed singly connected)  
    OUTLET              16585    8496     ok (non-closed singly connected)  
    SURFACE2              30361    16233    ok (non-closed singly connected)  
    SURFACE1              29804    15940    ok (non-closed singly connected)  
    DISCHARGE           12985    8327     ok (non-closed singly connected)  
    PLUG                6812     3534     ok (non-closed singly connected)  
    SIDE11              159185   83626    ok (non-closed singly connected)  
    SIDE22              160271   84167    ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-0.2085 0 -0.749957) (2.0434 1.0606 0.749957)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (4.31244e-16 4.57242e-14 -2.94413e-14) OK.
    Max cell openness = 4.94922e-15 OK.
    Max aspect ratio = 104.056 OK.
    Minimum face area = 2.14489e-09. Maximum face area = 0.00100591.  Face area magnitudes OK.
    Min volume = 1.39841e-13. Max volume = 1.0076e-05.  Total volume = 1.98826.  Cell volumes OK.
    Mesh non-orthogonality Max: 85.2692 average: 17.155
   *Number of severely non-orthogonal (> 70 degrees) faces: 15785.
    Non-orthogonality check OK.
  <<Writing 15785 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
    Max skewness = 3.76284 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

c. At the beginning, I was running without PISO but then wanted to see if it would get better with PISO. As for under-relaxation, I used them to achieve convergence within this time step (like a steady solution) but it did not work. I removed the under-relaxations later. No, it is a transient solution


It is very weird that especially my momentum equations are not converging. :/

[matejfor]

Hi ,

it§s not weird your case is not converging. Your mesh is not good. Your max. non-orthogonality is way to high. It should be below 70 ideally. 85 is too high.
Your aspect ratio should be within double digit number ideally.

What to do:
(1) go to upwind and generally first order differencing should help.
But the results will be colourful pictures only.
(2) get back to meshing and try a better mesh. if you are running latest OpenFOAM version 1706, you may run chekcMesh -writeAllFields to get scalar fields depicting the mesh quality parameters in paraview.

On other topics:
if the case is transient, your time step should be a fraction of your overall time. if you persist on CrankNicolson 0.9 you need to keep your maxCo below 1. But if you are happy with blending factor 0.5 you may grow higher.

The good start would be: nCorrectors 2; nOuterCorrectors 2; and see whether your pressure equation is converging nicely. If not, you may need to test increasing one or the other parameters. Your relaxation factors should be set to one.

Now the solvers. For pressure and for other variables as well, your:
tolerance 0;
relTol 1e-2;

should have a different values.
tolerance is an absolute tolerance which is the algebraic solver trying to reach. 0 is very bad number. 1e-8 is more likely.
relTol should be 0.1 or 0.01 OK. For velocity 0.1 is very reasonable value.
for pFinal, set relTol to 0.

What is actually the application?
And finally your discretisation is set to quite tight 2nd order setup. If you have a problem with convergence, it is a good idea to go towards the 1st order. In velocity you may try linearUpwind. For k and omega you really cannot use linear scheme. try upwind or linearUpwind.


Now I will fabulate from the setup you have shown, that you are either solving really specific problem, or you are not very solid in CFD with openfoam.

I would suggest looking at the collection of the tutorials here: https://wiki.openfoam.com/Tutorials
for a start.

hope this helps

Matej

[cyln]

Hi Matej,

Thanks for your recommendations. I am regenerating the mesh. I will retry with a better mesh quality, along with your other suggestions. I was kind of avoiding a new mesh and trying to get around the convergence problem with other changes because I map some experimental data at the boundaries whenever I use a new mesh and it takes time to map data in openfoam. I am using ICEM and then converting the mesh to openfoam.

I was setting tolerance to zero because I wanted OpenFOAM to use relTol before the pressure equation is solved. I defined the residual control additionally as 1e-6. Anyway, I set tolerance to 1e-8;

The discretization is purposely set to second-order due to accuracy reasons. The application is shear layer mixing and the flow is highly three dimensional. Certain vortical structures need to be captured accurately.

I will write back when I see my residuals after the changes

[matejfor]

Hi cyln,

second order is OK, but your original settings was a complete overkill especially with the mesh quality you have.

ICEM hexa? or ICEM tet?

OpenFOAM tends to work better with hexahedral based mesh. (if you have tets, you may try to dualise them into the fully polyhedral mesh)

Mat

[cyln]

Hi Matej,


My mesh includes tetrahedral cells and prism layers (to capture the BL). Unfortunately, my domain is not really suitable for an hexahedral mesh.

Let me see what I can do.

[cyln]

Hi Matej,

I improved my mesh quality with a max non-orhoganility of 74. I can see the improvement from the reduced number of iterations required to solve momentum equations.

However, there is another problem I have had since the beginning. When I run my simulation, I was getting " Floating point exception(core dumped) " error. I checked my boundary conditions but I could not find anything. The error is as follows:

PIMPLE: max iterations = 100
    field "(p|U)"	: relTol 0, tolerance 1e-06

Reading field p

AMI: Creating addressing and weights between 215211 source faces and 215315 target faces
AMI: Patch source sum(weights) min/max/average = 0, 1, 0.999643
AMI: Patch target sum(weights) min/max/average = 0, 1, 0.999657
Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type RAS
Selecting RAS turbulence model kOmegaSST
Selecting patchDistMethod meshWave
bounding k, min: 0 max: 2.4895 average: 0.9636
bounding omega, min: 0 max: 3.06e+07 average: 89.61
kOmegaSSTCoeffs
{
    alphaK1         0.85;
    alphaK2         1;
    alphaOmega1     0.5;
    alphaOmega2     0.856;
    gamma1          0.555556;
    gamma2          0.44;
    beta1           0.075;
    beta2           0.0828;
    betaStar        0.09;
    a1              0.31;
    b1              1;
    c1              10;
    F3              false;
}

No MRF models present

No finite volume options present


Starting time loop

Courant Number mean: 2.63784e-06 max: 2.31795
Time = 0.00015

PIMPLE: iteration 1
--> FOAM Warning : 
    From function void Foam::lduMatrix::operator-=(const Foam::lduMatrix&)
    in file matrices/lduMatrix/lduMatrix/lduMatrixOperations.C at line 286
    Unknown matrix type combination
    this : diagonal:0 symmetric:0 asymmetric:1
    A    : diagonal:0 symmetric:0 asymmetric:0
   Normalisation factor = 0.876324
DILUPBiCG:  Iteration 0 residual = 1
DILUPBiCG:  Iteration 1 residual = 0.562062
DILUPBiCG:  Iteration 2 residual = 0.160759
DILUPBiCG:  Iteration 3 residual = 0.0585505
DILUPBiCG:  Iteration 4 residual = 0.0224445
DILUPBiCG:  Iteration 5 residual = 0.00881576
DILUPBiCG:  Solving for Ux, Initial residual = 1, Final residual = 0.00881576, No Iterations 5
   Normalisation factor = 1.70916e-07
DILUPBiCG:  Iteration 0 residual = 1
DILUPBiCG:  Iteration 1 residual = 0.266825
DILUPBiCG:  Iteration 2 residual = 0.108854
DILUPBiCG:  Iteration 3 residual = 0.0295189
DILUPBiCG:  Iteration 4 residual = 0.00901428
DILUPBiCG:  Solving for Uy, Initial residual = 1, Final residual = 0.00901428, No Iterations 4
   Normalisation factor = 1.88245e-07
DILUPBiCG:  Iteration 0 residual = 1
DILUPBiCG:  Iteration 1 residual = 0.381921
DILUPBiCG:  Iteration 2 residual = 0.115358
DILUPBiCG:  Iteration 3 residual = 0.0400545
DILUPBiCG:  Iteration 4 residual = 0.0117323
DILUPBiCG:  Iteration 5 residual = 0.00325518
DILUPBiCG:  Solving for Uz, Initial residual = 1, Final residual = 0.00325518, No Iterations 5
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::divide(Foam::Field<double>&, double const&, Foam::UList<double> const&) at ??:?
#4  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::dimensioned<double> const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) in "/opt/openfoam4/platforms/linux64GccDPInt32Opt/bin/pimpleFoam"
#5  ? in "/opt/openfoam4/platforms/linux64GccDPInt32Opt/bin/pimpleFoam"
#6  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#7  ? in "/opt/openfoam4/platforms/linux64GccDPInt32Opt/bin/pimpleFoam"

[1]+  Floating point exception(core dumped) pimpleFoam > log &

Then I started to use the following code:

unset FOAM_SIGFPE

This let the simulation continue. With the new mesh and setup, After Pimple Iteration 3, my residuals start to diverge. Before that, everything goes well. I think that my setup and mesh were so bad before, that this divergence was delayed.

I am running a turbulent case. I tried the laminar case as well to check if it was due to my turbulence BCs. I got the same error in both cases.

I am quite sure this floating point error is not related to my mesh since I have got the same error with 3 different meshes.


Have you ever experienced this error?

[sheaker]

Hello. Maybe try:
nNonOrthogonalCorrectors 5;
in Your fvSolution, pimple. It looks like this options fits perfectly for Your mesh.
Have a nice day.
Sheaker

[matejfor]

Excellent progress. You are having FPE - floating point exception which happens 99% cases in situation when you are dividing by zero. I do not believe increasing number of correctors will help here, what my guess would be is your initial conditions are not good.

Also it si recommended, if your mesh is more tetrahedral then hexahedral to use for gradients not Gauss integration but leastSquares Look here for comparison https://www.openfoam.com/documentati...t-example.html.


But your main problem I believe are the initial conditions. Do not start from zero, start from potentialFoam solution to see any difference. You can also use FO to dump excessive velocity in the domain https://www.openfoam.com/releases/op...ityConstraints

[cyln]

Well, I am having a division by zero error (FPE) for potentialFoam as well. Since we are solving the laplace equation, this must be something to do with my velocity field or ICs.

As far as I know, initializing flow field with potentialFoam is quite a new feature in OF. I do not know how long it has been though. Since I did it for the first time, I show you how I did it below. I added the following to fvSolution file.

Phi
{
solver  GAMG;
smoother DIC;

tolerance  1e-6;
relTol     0.01;
}

potentialFlow
{
  nNonOrthogonalCorrectors 10;
}

Phi
{
$p;
}

My velocity field (0/U) is defined as follows:
Note that there is mapped data at INLET1 and I shortened here.

internalField   uniform (1 1 1);

boundaryField
{
    INLET1
    {
       type    fixedValue;
       value  nonuniform List<vector>
        2673
        (
    ( 17.4495   0.0000   0.0000 )
               ...
               ...
    (  4.7022   0.0000   0.0000 )       
        )
    ;
    }

    INLET2
    {
         type            zeroGradient;
    }

    SIDE11
    {
        type              cyclicAMI;
    }

    SIDE22
    {
        type              cyclicAMI;
    }

    CYL
    {
        type            slip;
    }

    OUTLET
    {
        type            zeroGradient;
    }

    SURFACE1
    {
       type            noSlip;
    }

    SURFACE2
    {
        type            noSlip;
    }

    DISCHARGE
    {
        type            noSlip;
    }

    PLUG
    {
        type            noSlip;
    }

}

My pressure field (0/p) is as follows:

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

internalField   uniform 1.0;

boundaryField
{
    INLET1
    {
        type            zeroGradient;
    }

    INLET2
    {
        type            totalPressure;
        p0              uniform 100000;
        value           uniform 100000;
        gamma            1.4;
    }

    SIDE11
    {
        type            cyclicAMI;
    }

    SIDE22
    {
        type            cyclicAMI;
    }

    CYL
    {
        type            zeroGradient;


    }

    OUTLET
    {
        type            fixedValue;
        value           uniform 1.0;
    }

    SURFACE1
    {
       type            zeroGradient;

    }

    SURFACE2
    {
       type            zeroGradient;

    }

    DISCHARGE
    {
       type            zeroGradient;

    }

    PLUG
    {
        type            zeroGradient;
    
    }

}

What am I missing here?

[cyln]

I changed my cyclicAMI boundary condition to symmetry just to check it was the problem. It did not give me any FPE error, however, momentum equations still diverge after the Pimple Iteration 3 as it did when i did "unset FOAM_SIGFPE" with cyclicAMI.

I dont know what this means. :/

[matejfor]

Well the first thing I've noticed is the initial velocity field.
It is set to (1 1 1) but the mapped inlet seems to set velocities only in X direction.
There might be a reason for that I do not see, but it caught my eye.

Inlet 2 with totalPressure and zeroGradient is potentially dangerous as depending on your pressure you might have got a micro-pumping (pure numerical one) across the inlet. So I would suggest to use pressureInletOutletVelocity which is basically a zeroGradient, but allows only normal direction inlet to make solution more stable.

mat

[cyln]

Hi Matej,

You were right about the zeroGradient pressure BC. Honestly, i set it as pressureInletOutletVelocity at the beginning and then I changed it to zeroGradient since I was trying to find the error.

Now I switched back to pressureInletOutletVelocity BC and I got rid of the divergence after pimple iteration 3. Now i reduced the number of problems to 1 only which is FPE. SIDE11 and SIDE22 boundaries are set to symmetry BC at this situation and that is how I got rid of the FPE error. These boundaries must be set to cyclicAMI boundary condition as before. However, when I do that, I am receiving FPE error. There are are people who had the same problem with cyclicAMI BC. See the links:

Changing to cyclic gives Floating point Exeptions

simpleFoam: floating point exception

Also, yea, internal field was set to (0 0 0), then set to (1 1 1). At the ccurrent situation, I set it back to (0 0 0) and my residuals seem to converge well with symmetric BCs on the sides.



Do you or someone else have experience with cyclicAMI boundary condition? To give more detail, I set these BCs as follows:

- Since I am generating the mesh using ICEM, I was setting SIDE11 and SIDE22 to periodic BC (called cyclic in openfoam). When I convert my mesh to foam, the periodic BC was giving error. So some other people set the boundary as symmetry on ICEM, converted the mesh to foam and then converted the BC from symmetry to cyclicAMI. That is how I did it.

This is boundary file after converting from symmetry to cyclicAMI

10
(
    INLET1
    {
        type            patch;
        nFaces          2673;
        startFace       20264662;
    }
    INLET2
    {
        type            patch;
        nFaces          13206;
        startFace       20267335;
    }
    CYL
    {
        type            wall;
        inGroups        1(wall);
        nFaces          65099;
        startFace       20280541;
    }
    OUTLET
    {
        type            patch;
        nFaces          47929;
        startFace       20345640;
    }
    SURFACE2
    {
        type            wall;
        inGroups        1(wall);
        nFaces          12944;
        startFace       20393569;
    }
    SURFACE1
    {
        type            wall;
        inGroups        1(wall);
        nFaces          13221;
        startFace       20406513;
    }
    DISCHARGE
    {
        type            wall;
        inGroups        1(wall);
        nFaces          436;
        startFace       20419734;
    }
    PLUG
    {
        type            wall;
        inGroups        1(wall);
        nFaces          4376;
        startFace       20420170;
    }
    SIDE11
    {
        type            cyclicAMI;
        inGroups        1(cyclicAMI);
        nFaces          215211;
        startFace       20424546;
        matchTolerance  0.001;
        transform       rotational;
        neighbourPatch  SIDE22;
        rotationAxis    (1 0 0);
        rotationCentre  (0 0 0);
    }
    SIDE22
    {
        type            cyclicAMI;
        inGroups        1(cyclicAMI);
        nFaces          215315;
        startFace       20639757;
        matchTolerance  0.001;
        transform       rotational;
        neighbourPatch  SIDE11;
        rotationAxis    (1 0 0);
        rotationCentre  (0 0 0);
    }
)

[cyln]

I kept it running for a while (with symmetry BCs on the sides). I see that momentum equations converge really fast, however, the pressure residuals are getting crazy. They keep building up through pimple iterations. The most interesting thing is that after the pressure residuals are built up, I still get an FPE (with symmetric BCs set).

I increased nNonOrthogonalCorrectors, tried different numbers. The general trend in the residuals is attached.

Crazy questions in my mind:
1- Is cyclicAMI causing the FPE error?
2- Why am I getting an FPE error when symmetry BCs are set on SIDE11 and SIDE22? This may suggest cyclicAMI is not the reason for FPE error.

[cyln]

Hello,

the control parameters I changed one by one are as follows:

1. Mesh: After optimizing non-orthoganality, I tried different meshes. My simulation still started diverging, leading to FPE later.

2. Discretization schemes are as set to second-order accuracy. Any change in the discretization do not stop the divergence, the simulation started diverging at pimple iteration 3 again.

3. My boundary, U and p files are as shown above. They are straigtforward. I dont see any issue with them.

4. Switching from cyclicAMI to symmetry still does not help. As far as I can understand, the divergence issue is not related to cyclicAMI.

5. I changed the number of non-orthogonal correctors and ncorrectors one by one. It did not help

6. I used potentialFoam to initialize my domain. It did not help



Whatever I change, my simulation starts diverging at pimple iteration 3, leading to a FPE error. This must be suggesting something. Does anyone see the problem here?

[cyln]

Hi Matej,

I found the route of the problem. There were two things:

1. Relaxation factors: Pimple algorithm is a combination of piso and simple. Therefore, correct implementation of pimple requires the usage of under-relaxation factors. Otherwise, the simulation diverges quickly as I experienced.

2. Courant number. Even though it is said to be able to work with high Co numbers, it does not tolerate very high Co numbers. Its value still should be kept below certain value.

I figured these out when I set my SIDE11 and SIDE22 boundaries to symmetry and the simulation is laminar. My residuals are attached.

As a next move, i converted symmetry BCs to cyclicAMI and I had a FPE error during the solution of pressure equation. I hope to figure it out