CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM

Bounding epsilon error while running LRR/SSG using interFoam

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 1 Post By Tobi
  • 1 Post By Tobi

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 26, 2021, 08:27
Default Bounding epsilon error while running LRR/SSG using interFoam
  #1
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 19
Rep Power: 2
subhojitkadiacfd is on a distinguished road
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

Last edited by Tobi; March 3, 2021 at 04:45.
subhojitkadiacfd is offline   Reply With Quote

Old   March 2, 2021, 08:02
Default
  #2
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,552
Blog Entries: 6
Rep Power: 46
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
Quote:
Originally Posted by subhojitkadiacfd View Post
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.
subhojitkadiacfd likes this.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   March 2, 2021, 08:54
Default
  #3
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 19
Rep Power: 2
subhojitkadiacfd is on a distinguished road
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

Last edited by Tobi; March 3, 2021 at 04:46. Reason: Tobi: Cleaned the mess
subhojitkadiacfd is offline   Reply With Quote

Old   March 2, 2021, 11:22
Default
  #4
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,552
Blog Entries: 6
Rep Power: 46
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
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).
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   March 2, 2021, 11:57
Default
  #5
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 19
Rep Power: 2
subhojitkadiacfd is on a distinguished road
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
subhojitkadiacfd is offline   Reply With Quote

Old   March 2, 2021, 12:16
Default
  #6
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,552
Blog Entries: 6
Rep Power: 46
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
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
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   March 3, 2021, 04:32
Default
  #7
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 19
Rep Power: 2
subhojitkadiacfd is on a distinguished road
Quote:
Originally Posted by Tobi View Post
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).
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?
subhojitkadiacfd is offline   Reply With Quote

Old   March 3, 2021, 04:47
Default
  #8
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,552
Blog Entries: 6
Rep Power: 46
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
Very important -> https://letmegooglethat.com/?q=fvOptions+limit+velocity



Including limiters in fvOptions
subhojitkadiacfd likes this.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   March 3, 2021, 08:01
Default
  #9
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 19
Rep Power: 2
subhojitkadiacfd is on a distinguished road
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
subhojitkadiacfd is offline   Reply With Quote

Old   March 3, 2021, 08:05
Default
  #10
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,552
Blog Entries: 6
Rep Power: 46
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
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.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   March 3, 2021, 08:45
Default
  #11
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 19
Rep Power: 2
subhojitkadiacfd is on a distinguished road
Quote:
Originally Posted by Tobi View Post
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.
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
subhojitkadiacfd is offline   Reply With Quote

Old   March 3, 2021, 09:23
Default
  #12
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,552
Blog Entries: 6
Rep Power: 46
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
I am not sure if the turbulence models you use are suitable for your case (VOF). No idea actually.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   March 3, 2021, 13:13
Default
  #13
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 19
Rep Power: 2
subhojitkadiacfd is on a distinguished road
Many thanks Tobi. Lets see how it goes.
subhojitkadiacfd is offline   Reply With Quote

Reply

Tags
bounding error, interfoam diverging, multiphase flow, reynolds stress model

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
interfoam parallel running exit code 145 fatlem OpenFOAM Running, Solving & CFD 5 July 27, 2019 07:48
SimpleFoam k and epsilon bounded nedved OpenFOAM Running, Solving & CFD 16 March 4, 2017 08:30
How to write k and epsilon before the abnormal end xiuying OpenFOAM Running, Solving & CFD 8 August 27, 2013 15:33
SimpleFoam k and epsilon bounded nedved OpenFOAM Running, Solving & CFD 1 November 25, 2008 20:21
Bounding epsilon and K with rasInterFoam openfoam_user OpenFOAM Running, Solving & CFD 0 October 23, 2008 08:48


All times are GMT -4. The time now is 17:01.