[robyTKD]

Hi Foamers,

I am working on a pimpleFoam simulation with MRF option activated. I have problem in making the simulation stable to perform 6 seconds. As you can see below, deltaT suddenly decreases and after a lot of iterations the simulation stops with an error. Moreover I encounter a bounding omega warning (with high number of iteration to solve omega equation 1001) before the decrease of deltaT. Is it related to the increasing number of iteration to solve p?

Courant Number mean: 0.219033 max: 10.0022
deltaT = 0.000256884
Time = 1.69945

DILUPBiCG:  Solving for Ux, Initial residual = 0.00050968, Final residual = 5.31923e-06, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.000478035, Final residual = 3.74021e-06, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.000274494, Final residual = 8.97499e-06, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.00760123, Final residual = 4.27931e-05, No Iterations 5
time step continuity errors : sum local = 4.131e-08, global = -9.20692e-10, cumulative = -4.54905e-06
GAMG:  Solving for p, Initial residual = 0.00159695, Final residual = 8.68188e-07, No Iterations 8
time step continuity errors : sum local = 8.38347e-10, global = 2.37631e-11, cumulative = -4.54902e-06
DILUPBiCG:  Solving for omega, Initial residual = 0.00013723, Final residual = 0.00808422, No Iterations 1001
bounding omega, min: -91365.5 max: 126745 average: 823.181
DILUPBiCG:  Solving for k, Initial residual = 0.000550298, Final residual = 6.53352e-06, No Iterations 2
ExecutionTime = 7531.18 s  ClockTime = 7552 s

Courant Number mean: 0.21883 max: 9.99752
deltaT = 0.000256884
Time = 1.6997

DILUPBiCG:  Solving for Ux, Initial residual = 0.676568, Final residual = 2.13076e-06, No Iterations 4
DILUPBiCG:  Solving for Uy, Initial residual = 0.196484, Final residual = 5.24048e-06, No Iterations 3
DILUPBiCG:  Solving for Uz, Initial residual = 0.341194, Final residual = 7.24752e-06, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.0180928, Final residual = 0.000168046, No Iterations 4
time step continuity errors : sum local = 1.63347e-07, global = 1.0021e-08, cumulative = -4.539e-06
GAMG:  Solving for p, Initial residual = 0.851594, Final residual = 9.52237e-07, No Iterations 22
time step continuity errors : sum local = 1319.89, global = 33.0263, cumulative = 33.0263
DILUPBiCG:  Solving for omega, Initial residual = 1, Final residual = 2404.53, No Iterations 1001
bounding omega, min: -1.30647e+29 max: 1.85931e+28 average: 3.31278e+27
DILUPBiCG:  Solving for k, Initial residual = 0.401606, Final residual = 8.71858e-14, No Iterations 1
ExecutionTime = 7543.81 s  ClockTime = 7564 s

Courant Number mean: 2.44738e+09 max: 3.10915e+15
deltaT = 8.2622e-19
--> FOAM Warning :
    From function Time::operator++()
    in file db/Time/Time.C at line 1055
    Increased the timePrecision from 6 to 7 to distinguish between timeNames at time 1.6997
Time = 1.699702

DILUPBiCG:  Solving for Ux, Initial residual = 0.0894012, Final residual = 7.65914e-06, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.0923353, Final residual = 4.08565e-06, No Iterations 3
DILUPBiCG:  Solving for Uz, Initial residual = 0.0676507, Final residual = 1.5654e-07, No Iterations 5
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.0062001, No Iterations 5
time step continuity errors : sum local = 1.56288e-08, global = 2.77555e-09, cumulative = 33.0263
GAMG:  Solving for p, Initial residual = 0.0547949, Final residual = 8.07853e-07, No Iterations 15
time step continuity errors : sum local = 5.48245e-12, global = 2.31016e-13, cumulative = 33.0263
DILUPBiCG:  Solving for omega, Initial residual = 0.872425, Final residual = 3.91289e-15, No Iterations 1
bounding omega, min: -7.84392e+23 max: 3.3123e+25 average: 4.83913e+20
DILUPBiCG:  Solving for k, Initial residual = 0.59834, Final residual = 1.01994e-06, No Iterations 1
ExecutionTime = 7545.14 s  ClockTime = 7566 s

Courant Number mean: 3.23617e-06 max: 3.33316
deltaT = 9.91465e-19
--> FOAM Warning :
    From function Time::operator++()
    in file db/Time/Time.C at line 1055
    Increased the timePrecision from 7 to 9 to distinguish between timeNames at time 1.6997
Time = 1.69970204

DILUPBiCG:  Solving for Ux, Initial residual = 0.480664, Final residual = 4.31292e-06, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.243898, Final residual = 4.23116e-07, No Iterations 3
DILUPBiCG:  Solving for Uz, Initial residual = 0.205338, Final residual = 6.8143e-09, No Iterations 4
GAMG:  Solving for p, Initial residual = 0.415985, Final residual = 0.00243229, No Iterations 5
time step continuity errors : sum local = 2.13818e-08, global = -3.8758e-09, cumulative = 33.0263
GAMG:  Solving for p, Initial residual = 0.0894586, Final residual = 9.57609e-07, No Iterations 17
time step continuity errors : sum local = 3.86934e-12, global = -1.4289e-13, cumulative = 33.0263
DILUPBiCG:  Solving for omega, Initial residual = 0.999999, Final residual = 4.84813e-16, No Iterations 1
bounding omega, min: -2.76819e+08 max: 1.25494e+20 average: 5.66184e+18
DILUPBiCG:  Solving for k, Initial residual = 0.0253493, Final residual = 4.64098e-06, No Iterations 1
ExecutionTime = 7546.5 s  ClockTime = 7567 s

Courant Number mean: 3.28599e-06 max: 2.21188
deltaT = 1.18976e-18
--> FOAM Warning :
    From function Time::operator++()
    in file db/Time/Time.C at line 1055
    Increased the timePrecision from 9 to 10 to distinguish between timeNames at time 1.6997
Time = 1.699702042

DILUPBiCG:  Solving for Ux, Initial residual = 0.114399, Final residual = 6.89955e-07, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.06611, Final residual = 2.26779e-07, No Iterations 3
DILUPBiCG:  Solving for Uz, Initial residual = 0.0669545, Final residual = 3.6835e-06, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.28107, Final residual = 0.00247868, No Iterations 3
time step continuity errors : sum local = 1.12853e-08, global = -6.20024e-09, cumulative = 33.0263
GAMG:  Solving for p, Initial residual = 0.0393327, Final residual = 7.7993e-07, No Iterations 15
time step continuity errors : sum local = 2.61803e-12, global = -1.80286e-13, cumulative = 33.0263
DILUPBiCG:  Solving for omega, Initial residual = 0.806297, Final residual = 1.26278e-07, No Iterations 1
bounding omega, min: -1642.61 max: 9.39139e+18 average: 3.42263e+18
DILUPBiCG:  Solving for k, Initial residual = 0.0195621, Final residual = 3.98071e-06, No Iterations 1
ExecutionTime = 7547.77 s  ClockTime = 7568 s

Courant Number mean: 3.63278e-06 max: 1.91779
deltaT = 1.42771e-18
--> FOAM Warning :
    From function Time::operator++()
    in file db/Time/Time.C at line 1055
    Increased the timePrecision from 10 to 11 to distinguish between timeNames at time 1.6997
Time = 1.6997020424

DILUPBiCG:  Solving for Ux, Initial residual = 0.0618404, Final residual = 4.24347e-07, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.039055, Final residual = 6.00065e-06, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.0439426, Final residual = 1.02212e-06, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.232859, Final residual = 0.00135369, No Iterations 3
time step continuity errors : sum local = 5.35432e-09, global = -3.08645e-09, cumulative = 33.0263
GAMG:  Solving for p, Initial residual = 0.0231275, Final residual = 8.75852e-07, No Iterations 13
time step continuity errors : sum local = 2.85843e-12, global = -1.94964e-13, cumulative = 33.0263
DILUPBiCG:  Solving for omega, Initial residual = 0.751881, Final residual = 1.51817e-07, No Iterations 1
bounding omega, min: -16.8098 max: 4.45048e+18 average: 2.32139e+18
DILUPBiCG:  Solving for k, Initial residual = 0.0161215, Final residual = 3.89756e-06, No Iterations 1
ExecutionTime = 7548.98 s  ClockTime = 7570 s

Courant Number mean: 4.13771e-06 max: 1.81145
deltaT = 1.71325e-18
--> FOAM Warning :
    From function Time::operator++()
    in file db/Time/Time.C at line 1055
    Increased the timePrecision from 11 to 12 to distinguish between timeNames at time 1.6997
Time = 1.69970204241

DILUPBiCG:  Solving for Ux, Initial residual = 0.0513661, Final residual = 3.59432e-07, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.0339824, Final residual = 6.19584e-06, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.037676, Final residual = 3.88398e-07, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.188385, Final residual = 0.000806453, No Iterations 3
time step continuity errors : sum local = 3.3117e-09, global = -1.69835e-09, cumulative = 33.0263
GAMG:  Solving for p, Initial residual = 0.016901, Final residual = 8.37202e-07, No Iterations 11
time step continuity errors : sum local = 2.97013e-12, global = 6.06414e-14, cumulative = 33.0263
DILUPBiCG:  Solving for omega, Initial residual = 0.712028, Final residual = 1.71042e-07, No Iterations 1
bounding omega, min: -118.702 max: 2.72813e+18 average: 1.67489e+18
DILUPBiCG:  Solving for k, Initial residual = 0.0137967, Final residual = 3.95515e-06, No Iterations 1
ExecutionTime = 7550.13 s  ClockTime = 7571 s

Courant Number mean: 4.7452e-06 max: 1.80144
deltaT = 2.0559e-18
--> FOAM Warning :
    From function Time::operator++()
    in file db/Time/Time.C at line 1055
    Increased the timePrecision from 12 to 13 to distinguish between timeNames at time 1.6997
Time = 1.699702042413

DILUPBiCG:  Solving for Ux, Initial residual = 0.0502811, Final residual = 2.77759e-07, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.0329875, Final residual = 5.51329e-06, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.0364956, Final residual = 2.06106e-07, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.16494, Final residual = 0.000553674, No Iterations 3
time step continuity errors : sum local = 2.48088e-09, global = -8.1257e-10, cumulative = 33.0263
GAMG:  Solving for p, Initial residual = 0.0137505, Final residual = 6.89309e-07, No Iterations 12
time step continuity errors : sum local = 2.68517e-12, global = 2.31621e-13, cumulative = 33.0263
DILUPBiCG:  Solving for omega, Initial residual = 0.682025, Final residual = 1.90232e-07, No Iterations 1
bounding omega, min: -0.0388451 max: 1.86296e+18 average: 1.25543e+18
DILUPBiCG:  Solving for k, Initial residual = 0.0120913, Final residual = 4.33272e-06, No Iterations 1
ExecutionTime = 7551.31 s  ClockTime = 7572 s

I tried to decrease the CFL number to 1, but deltaT is very low, in particular its magnitude is 10^-28.

Regarding the mesh: checkMesh is OK, with a maximum nonOrthogonality of 65.9 degrees and a average nonOrtho of 7 degrees.

Regarding the fvSchemes:

ddtSchemes
{
    default         Euler;
}

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

divSchemes
{
    default         none;

    div(phi,U)      Gauss linearUpwindV grad(U);
    div(phi,k)      Gauss upwind;
    div(phi,omega)  Gauss upwind;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default                     Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p;
}

How can I stabilize the simulation? Do you have any idea?

Thanks,
Roberto

[pbachant]

What does your fvSolution file look like?

[robyTKD]

Hi Pete, thank you for the quick reply. Here is my fvSolution

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.01;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    pFinal
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    "(U|k|omega|epsilon)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-05;
        relTol          0.1;
    }

    "(U|k|omega|epsilon)Final"
    {
        $U;
        tolerance       1e-05;
        relTol              0;
    }
}

PIMPLE
{
    // momentumPredictor                yes;
    nOuterCorrectors            1;
    nCorrectors                 2;
    nNonOrthogonalCorrectors    0;
    pRefCell                    0;
    pRefValue                   0;
}

relaxationFactors
{
    fields
    {
    }
    equations
    {
        "p.*"           0.3;
        "U.*"           0.7;
        "k.*"           0.7;
        "omega.*"     0.7;
    }
}

[alexeym]

Hi,

as you've got diverging solution (primary reason for omega going unbound and catastrophic decrease of time step), you can:

1. Increase number of outer correctors
2. Increase number of non-orthogonal correctors
3. Switch from fixed number of outer correctors to residual control to exit outer corrector loop

[robyTKD]

Thanks.

1. Already done, no luck.
2. Is it necessary although I have low nonOrthogonality?
3. How can I switch to residual control to exit outer correction loop?

Another important thing I just saw in my case: running on 48 processors I met the problem above while running on 56 it ran without problems.

Have you any idea?

Thank you in advance

[alexeym]

Well,

1. Can you tell that increased value?

2. You've said - "...maximum nonOrthogonality of 65.9 degrees..."

3. Change your PIMPLE sub-dictionary from

PIMPLE
{
    // momentumPredictor                yes;
    nOuterCorrectors            1;
    nCorrectors                 2;
    nNonOrthogonalCorrectors    0;
    pRefCell                    0;
    pRefValue                   0;
}

to

PIMPLE
{
    // momentumPredictor                yes;
    nOuterCorrectors            100;
    nCorrectors                 2;
    nNonOrthogonalCorrectors    1;
    pRefCell                    0;
    pRefValue                   0;

    residualControl
    {
        "(p|U|k|omega)"
        {
            tolerance 1e-2;
            relTol 0;
        }
    }
}

4. How many cells you mesh has?

    p
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.01;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

with nCellsInCoarsestLevel 10; the solver can go in a very strange direction. Change it to PCG first, if it works, then change nCellsInCoarsestLevel to sqrt(number of cells in mesh).

5. Are you sure about linearUpwindV? Change it to upwind first, if solution will not diverge, then you can go with higher order schemes.

And finally about 48 and 56 processors: as you do not check convergence of the solution on every time step, there can be anything (different time stepping, more lucky distribution of cells).

Update: Also with 56 processors, you've got smaller subdomains, maybe with these smaller subdomains GAMG behaves better with your nCellsInCoarsestLevel 10 settings.

[robyTKD]

Thank you very much for your advices.

I will perform these tests as soon as possible and I will post the results.

[robyTKD]

Sorry for my late reply.

I performed simulation with your advices and now it works very well. Thank you very much.

Just for clarity:

- Increasing the number of Outer and nonOrthogonal correctors didn't give me any stability;
- Following your advice for GAMG increased the stability, also in parallel!!!
- Using upwind scheme increased the stability!

Thank you again.
Roberto

[alexeym]

Hi,

Originally Posted by robyTKD

- Using upwind scheme increased the stability!

Keep in mind, upwind is a first order scheme, so in some cases with this scheme you can flow can loose certain flow features (for example in certain melting problems first order schemes do not resolve all flow cells).

About increased number of outer corrector iterations. When you use residual control to exit outer corrector loop, you at least can be sure that your solution is converged. Otherwise you can't be sure in the results you get.

[1988]

hi guys
My simulation has the same problem that here is mentioned ,I have read all these posts but I don't know how to solve it .these are my fv scheme and fv solution

ddtSchemes
{
    default         Euler;
}

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

divSchemes
{
    default         none;
    div(phi,U)      Gauss Gamma 0.2;
    div(phi,k)      Gauss Gamma 0.2;
    div(phi,epsilon) Gauss Gamma 0.2;
    div(phi,R)      Gauss Gamma 0.2;
    div(R)          Gauss Gamma 0.2;
    div(phi,nuTilda) Gauss Gamma 0.2;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         none;
    laplacian(nuEff,U) Gauss linear corrected;
    laplacian((1|A(U)),p) Gauss linear corrected;
    laplacian(DkEff,k) Gauss linear corrected;
    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
    laplacian(DREff,R) Gauss linear corrected;
    laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
    interpolate(U)  linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p               ;
}

solvers
{
     p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          1e-05;
    }

    pFinal
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-08;
        relTol          0;
    }

    "(U|k|epsilon)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-05;
        relTol          0.1;
    }

    "(U|k|epsilon)Final"
    {
        $U;
        tolerance       1e-05;
        relTol          0;
    }
}

PIMPLE
{
    nOuterCorrectors 4;
    nCorrectors     2;
    nNonOrthogonalCorrectors 3;
    pRefCell        0;
    pRefValue       0;
residualControl
    {
        U
        {
            tolerance 1e-2;
            relTol 0;
        }

        p
        {
            tolerance 1e-2;
            relTol 0;
        }
    }
}

relaxationFactors
{
    fields
    {
    }
    equations
    {
        "p.*"           0.3;
        "U.*"           0.7;
        "k.*"           0.4;
        "epsilon.*"     0.5;
    }
}

thanks