CFD Online Discussion Forums - Bounding epsilon error while running LRR/SSG using interFoam

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)

- - Bounding epsilon error while running LRR/SSG using interFoam (https://www.cfd-online.com/Forums/openfoam/234208-bounding-epsilon-error-while-running-lrr-ssg-using-interfoam.html)

Bounding epsilon error while running LRR/SSG using interFoam

Hello to every one,

I am trying to simulate supercritical multiphase flow (interFoam) using different RAS models to get the best solution. The simulation is running with k-e, RNG, k-wSST etc., but for LRR and SSG models, I am getting bounding epsilon message from the very beginning of the simulation. The epsilon blows up and the simulation stops after 1-2 seconds of simulation time.

Further, in LRR, SSG

Code:

div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; was not accepted and

div(((rho*nu)*dev2(T(grad(U))))) Gauss linear; was accepted.

The error message appears just after a few iterations

Code:

ExecutionTime = 89.92 s  ClockTime = 98 s



Courant Number mean: 0.298636 max: 0.895395

Interface Courant Number mean: 0.00270084 max: 0.895395

deltaT = 0.000471525

Time = 0.08411459



smoothSolver:  Solving for alpha.water, Initial residual = 0.00299815, Final residual = 1.53669e-09, No Iterations 3

Phase-1 volume fraction = 0.97951  Min(alpha.water) = -2.46878e-07  Max(alpha.water) = 1

MULES: Correcting alpha.water

MULES: Correcting alpha.water

Phase-1 volume fraction = 0.97951  Min(alpha.water) = -0.000233734  Max(alpha.water) = 1.00009

GAMG:  Solving for p_rgh, Initial residual = 0.00769725, Final residual = 3.59405e-05, No Iterations 2

GAMG:  Solving for p_rgh, Initial residual = 6.21739e-05, Final residual = 2.90274e-07, No Iterations 4

time step continuity errors : sum local = 6.90028e-08, global = 3.90763e-08, cumulative = 4.71314e-06

GAMG:  Solving for p_rgh, Initial residual = 6.95841e-05, Final residual = 2.37694e-07, No Iterations 3

GAMG:  Solving for p_rgh, Initial residual = 8.27069e-07, Final residual = 9.28369e-09, No Iterations 3

time step continuity errors : sum local = 2.20722e-09, global = -9.53455e-10, cumulative = 4.71219e-06

GAMG:  Solving for p_rgh, Initial residual = 1.6805e-05, Final residual = 4.85751e-08, No Iterations 3

GAMG:  Solving for p_rgh, Initial residual = 1.42639e-07, Final residual = 6.06823e-09, No Iterations 2

time step continuity errors : sum local = 1.44347e-09, global = -8.06204e-10, cumulative = 4.71138e-06

DILUPBiCG:  Solving for epsilon, Initial residual = 0.999478, Final residual = 3.1044e-09, No Iterations 3

bounding epsilon, min: -1.43351e-08 max: 0.00811 average: 4.5564e-08

DILUPBiCG:  Solving for Rxx, Initial residual = 0.967164, Final residual = 5.10025e-10, No Iterations 1

DILUPBiCG:  Solving for Rxy, Initial residual = 0.778008, Final residual = 1.6466e-09, No Iterations 1

DILUPBiCG:  Solving for Ryy, Initial residual = 0.98587, Final residual = 4.91203e-10, No Iterations 1

DILUPBiCG:  Solving for Rzz, Initial residual = 0.989121, Final residual = 5.00174e-10, No Iterations 1

ExecutionTime = 90.44 s  ClockTime = 99 s

Your suggestions, comments would be highly appreciated.

Subhojit

Quote:

Originally Posted by subhojitkadiacfd (Post 797521)

Further, in LRR, SSG

div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; was not accepted and
div(((rho*nu)*dev2(T(grad(U))))) Gauss linear; was accepted.

This is based on the fact that the stress tensor (which is split into deviatoric part --> tau and hydrostatic part --> pressure) looks like that and the implementation in the models are accordingly. However, I am not familiar with these two models but commonly I would expect to have the effective viscosity here (but maybe it is calculated differently here). One needs to check how the turbulent contribution is applied to the equations. Nevertheless, in RANS modelling approaches we commonly increase the mixing by increasing the molecular viscosity according to the turbulent viscosity.

Maybe a bug.

To the other topic you mention. If you use the LLR turbulence model, are you under-relaxing the REqn? You can make it at least diagonal dominant by setting the matrix relaxation to 1 (same for epsilon). Same is valid for SSG.

In addition, you did not provide any error message or other useful hints. fvSolution and the schemes you are using are of interest too.

Many thanks for your valuable comment Tobi. I have used relaxation of 1 for the above-mentioned terms. The fvSolution and fvSchemes are:

Code:

ddtSchemes

{

    default         Euler;

}



gradSchemes

{

    default         faceLimited leastSquares 0.5;

    grad(U)            faceLimited leastSquares 0.5;

    //grad(U)         cellLimited leastSquares 0.5;

    //grad(k)         cellLimited leastSquares 0.5;

    //grad(epsilon)   cellLimited leastSquares 0.5;

    //default         Gauss linear;

}



divSchemes

{

    default                none;

    div(rhoPhi,U)          Gauss linearUpwind grad(U);//Gauss LUST grad(U);

    div(phi,alpha)      Gauss vanLeer;

    div(phirb,alpha)     Gauss linear;

    div(phi,k)          Gauss upwind;

    div(phi,epsilon)     Gauss upwind;

    div(phi,omega)         Gauss upwind; //for kOmega, kOmegaSST

    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;

    div(phi,R)          Gauss upwind; //For LRR, SSG

    div(R)              Gauss linear; //For LRR, SSG

    div((rho*R))          Gauss linear; //For LRR, SSG

    div(((rho*nu)*dev2(T(grad(U))))) Gauss linear; //For LRR, SSG

    div(nonlinearStress)    Gauss linear; //For LRR, SSG

}



laplacianSchemes

{

    default Gauss linear limited corrected 0.33;

    //default         Gauss linear corrected;

}



interpolationSchemes

{

    default         linear;

}



snGradSchemes

{

    default limited corrected 0.33;

    //default         corrected;

}



wallDist // for LRR, SSG

{

    method            meshWave;

}



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

Code:

 solvers // water channel, dam break tutorials

{

    "alpha.water.*"

    {

        nAlphaCorr      2;

        nAlphaSubCycles 1;

        cAlpha          1;



         MULESCorr       yes;

        nLimiterIter    3;

        

        solver          smoothSolver;

        smoother        symGaussSeidel;

        tolerance       1e-8;

        relTol          0;

    }



    "pcorr.*"

    {

        solver          PCG;

        preconditioner

        {

            preconditioner  GAMG;

            tolerance       1e-05;

            relTol          0;

            smoother        GaussSeidel;

        }

        tolerance       1e-05;

        relTol          0;

        maxIter         50;

    }



    p_rgh

    {

        solver          GAMG;

        tolerance       1e-08;

        relTol          0.01;

        smoother        GaussSeidel;

    }



    p_rghFinal

    {

        $p_rgh;

        relTol          0;

    }



    "(U|k|epsilon|omega|R)" //add R for LRR and SSG, omega for kOmega, kOmegaSST

    {

        solver          PBiCG;

        preconditioner  DILU;

        tolerance       1e-06;

        relTol          0;

    }



    "(U|k|epsilon|omega|R)Final" //add R for LRR and SSG, omega for kOmega, kOmegaSST

    {

        $U;

        tolerance       1e-08;

        relTol          0;

    }

}



PIMPLE

{

    momentumPredictor off; // off for multiphase

    nCorrectors     3;

    nNonOrthogonalCorrectors 1;

    

    //residualControl

    //{

        U        1e-07;

        p_rgh    1e-06;

    //}

}



relaxationFactors

{

    equations

    {

        ".*"                            1;

        "(U|k|epsilon|omega|R)"            1;    

    }

}

After some time the epsilon blows up and the simulation stops

Code:

ExecutionTime = 804.44 s  ClockTime = 827 s



Courant Number mean: 1.97399e+30 max: 1.25118e+35

Interface Courant Number mean: 1.34134e+09 max: 1.12614e+13

deltaT = 2.07961e-39

Time = 1.719277



smoothSolver:  Solving for alpha.water, Initial residual = 1.14868e-11, Final residual = 1.14868e-11, No Iterations 0

Phase-1 volume fraction = 0.741719  Min(alpha.water) = -1.10889e-05  Max(alpha.water) = 1.00002

MULES: Correcting alpha.water

MULES: Correcting alpha.water

Phase-1 volume fraction = 0.741719  Min(alpha.water) = -1.10889e-05  Max(alpha.water) = 1.00002

GAMG:  Solving for p_rgh, Initial residual = 1, Final residual = 0.0050873, No Iterations 6

GAMG:  Solving for p_rgh, Initial residual = 0.00245794, Final residual = 1.94538e-05, No Iterations 5

time step continuity errors : sum local = 7.6869e-10, global = 5.68414e-10, cumulative = 3.94798e+30

GAMG:  Solving for p_rgh, Initial residual = 0.00318051, Final residual = 1.26821e-05, No Iterations 4

GAMG:  Solving for p_rgh, Initial residual = 7.93767e-05, Final residual = 6.21374e-07, No Iterations 4

time step continuity errors : sum local = 2.66209e-11, global = -1.98525e-11, cumulative = 3.94798e+30

GAMG:  Solving for p_rgh, Initial residual = 0.000329261, Final residual = 2.90286e-06, No Iterations 5

GAMG:  Solving for p_rgh, Initial residual = 2.09426e-05, Final residual = 7.86557e-09, No Iterations 8

time step continuity errors : sum local = 3.3573e-13, global = -2.25506e-13, cumulative = 3.94798e+30

DILUPBiCG:  Solving for epsilon, Initial residual = 0.999135, Final residual = 1.41005e-14, No Iterations 1

bounding epsilon, min: -3.97239e+38 max: 3.03117e+52 average: 2.01431e+48

DILUPBiCG:  Solving for Rxx, Initial residual = 0.979011, Final residual = 14.1644, No Iterations 1000

DILUPBiCG:  Solving for Rxy, Initial residual = 0.000176047, Final residual = 0.0894335, No Iterations 1000

DILUPBiCG:  Solving for Ryy, Initial residual = 0.999409, Final residual = 9.48781e-09, No Iterations 5

DILUPBiCG:  Solving for Rzz, Initial residual = 1, Final residual = 136864, No Iterations 1000

ExecutionTime = 813.45 s  ClockTime = 836 s



Courant Number mean: 3.99863e-05 max: 0.212594

Interface Courant Number mean: 7.22945e-09 max: 1.22022e-05

deltaT = 2.49554e-39

Time = 1.719277



smoothSolver:  Solving for alpha.water, Initial residual = 6.81458e-10, Final residual = 6.81458e-10, No Iterations 0

Phase-1 volume fraction = 0.741719  Min(alpha.water) = -1.10889e-05  Max(alpha.water) = 1.00002

MULES: Correcting alpha.water

MULES: Correcting alpha.water

Phase-1 volume fraction = 0.741719  Min(alpha.water) = -1.10889e-05  Max(alpha.water) = 1.00002

GAMG:  Solving for p_rgh, Initial residual = 0.26929, Final residual = 0.00255839, No Iterations 10

GAMG:  Solving for p_rgh, Initial residual = 3.55552e-73, Final residual = 3.55552e-73, No Iterations 0

time step continuity errors : sum local = 2.61007e-07, global = -1.32084e-07, cumulative = 3.94798e+30

GAMG:  Solving for p_rgh, Initial residual = 4.16704e-73, Final residual = 4.16704e-73, No Iterations 0

GAMG:  Solving for p_rgh, Initial residual = 4.16704e-73, Final residual = 4.16704e-73, No Iterations 0

time step continuity errors : sum local = 3.05898e-07, global = -1.08856e-07, cumulative = 3.94798e+30

GAMG:  Solving for p_rgh, Initial residual = 4.22678e-73, Final residual = 4.22678e-73, No Iterations 0

GAMG:  Solving for p_rgh, Initial residual = 4.22678e-73, Final residual = 4.22678e-73, No Iterations 0

time step continuity errors : sum local = 3.10284e-07, global = -1.06426e-07, cumulative = 3.94798e+30

DILUPBiCG:  Solving for epsilon, Initial residual = 5.59504e-21, Final residual = 5.59504e-21, No Iterations 0

DILUPBiCG:  Solving for Rxx, Initial residual = 1, Final residual = 1, No Iterations 1000

DILUPBiCG:  Solving for Rxy, Initial residual = 0.999986, Final residual = 0.999986, No Iterations 1000

DILUPBiCG:  Solving for Ryy, Initial residual = 1, Final residual = 1, No Iterations 1000

DILUPBiCG:  Solving for Rzz, Initial residual = 1, Final residual = 1, No Iterations 1000

ExecutionTime = 825.17 s  ClockTime = 848 s

-------------------------------------------------------
Primary job terminated normally, but 1 process returned
a non-zero exit code.. Per user-direction, the job has been aborted.
-------------------------------------------------------
--------------------------------------------------------------------------
mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

Process name: [[21610,1],0]
Exit code: 145

In ReynoldsStress.c file at line 243

Code:

- fvc::div(this->alpha_*rho*this->nu()*dev2(T(fvc::grad(U))))

and in bound.c file, the code of bounding a scalar field is present

Code:

const scalar minVsf = min(vsf).value();



    if (minVsf < lowerBound.value())

    {

        Info<< "bounding " << vsf.name()

            << ", min: " << minVsf

            << " max: " << max(vsf).value()

            << " average: " << gAverage(vsf.primitiveField())

            << endl;



        vsf.primitiveFieldRef() = max

        (

            max

            (

                vsf.primitiveField(),

                fvc::average(max(vsf, lowerBound))().primitiveField()

              * pos0(-vsf.primitiveField())

            ),

            lowerBound.value()

        );



        vsf.boundaryFieldRef() = max(vsf.boundaryField(), lowerBound.value());

    }



    return vsf;

}

So, what I understood is that I am getting -ve value of alpha and epsilon at some places which is not physical.

Thanking you.

Subhojit

First of all, please clean your posts. Its just a mess. We have the code tag to include dictionaries, code and others in order to keep the format.

Probably your problem starts somewhere before:

Code:

Courant Number mean: 1.97399e+30 max: 1.25118e+35

Interface Courant Number mean: 1.34134e+09 max: 1.12614e+13

deltaT = 2.07961e-39

Time = 1.719277

Anything went wrong here previously. Check out your states in paraview, e.g., at 1.70 s, 1,71 s 1,715 s to check where the problem starts and to analyze why your simulation crashes. As you can see, the Co Number kills everything. Probably somewhere is a high velocity (which can be limited by the fvOptions to make your simulation more stable).

Thanks. I tried several fvSchemes to limit the alpha between 0 and 1, but still some -ve values are there which is possibly generating the -ve epsilon. The problem is starting with bounding epsilon and after some time getting a very high value of epsilon (both +ve and -ve).

Subhojit

Bounding problems occur from time to time and therefore, we have the bounding functions that change the values to physical ones. The alpha field can have small negative values (this is a numerical topic) but I had such values in my calculations too. However, I never used the turbulence models you are using and hence, no idea what can cause the problem. In addition, we (the community) do not have insight into your case so well ... a bit tricky :)

Quote:

Originally Posted by Tobi (Post 797642)

First of all, please clean your posts. Its just a mess. We have the code tag to include dictionaries, code and others in order to keep the format.

Probably your problem starts somewhere before:

Code:

Courant Number mean: 1.97399e+30 max: 1.25118e+35

Interface Courant Number mean: 1.34134e+09 max: 1.12614e+13

deltaT = 2.07961e-39

Time = 1.719277

I have not used fvOptions file (in constant directory) in my simulations. Could you please give some example/reference of a fvOptions file to limit the velocity?

Quote:

Originally Posted by Tobi (Post 797733)

Very important -> https://letmegooglethat.com/?q=fvOptions+limit+velocity

https://www.cfd-online.com/Forums/op...fvoptions.html

Thanks for the suggestion. I am looking forward to use it. Do you feel I need to include the source term of epsilon in fvOptions too?

Subhojit

For epsilon you don't need anything and definitely no source term ;). I guess you meant if you need some settings in fvOptions for epsilon. Well, the limitVelocity fvOptions will try to limit the velocity if you exceed the maximum number you specified. However, that does not mean that your simulation will run fine/correct/without errors. Its just a try to stabilize the simulation.

However, something is going wrong here as you can calculate a few time steps. So I expect that anything causes the solver to crash. Limiting the velocity might help but will not resolve the error.

Quote:

Originally Posted by Tobi (Post 797757)

However, something is going wrong here as you can calculate a few time steps. So I expect that anything causes the solver to crash. Limiting the velocity might help but will not resolve the error.

Many thanks for the suggestion.
What I noticed is that at most of the cells in the domain, the epsilon is close to the minimum value of 2.2e-16, and the calculated k (from R) is close to the minimum value of 3.3e-16. And at some cells close to the interface, the epsilon is increasing rapidly.
The maximum velocity is below 10 m/s which is okay for my case.

Subhojit

I am not sure if the turbulence models you use are suitable for your case (VOF). No idea actually.

Many thanks Tobi. Lets see how it goes.

Hi, have you solved your problem? Recently， I meet the same problem，do you have some experience to share with me?