[subhojitkadiacfd]

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

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

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

[Tobi]

Originally Posted by subhojitkadiacfd

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]

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

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;
}

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

 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

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

          - 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

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

[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:

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).

[subhojitkadiacfd]

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

[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

[subhojitkadiacfd]

Originally Posted by 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:

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?

[Tobi]

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



Including limiters in fvOptions

[subhojitkadiacfd]

Originally Posted by Tobi

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



Including limiters in fvOptions

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

[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.

[subhojitkadiacfd]

Originally Posted by 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.

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

[Tobi]

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

[subhojitkadiacfd]

Many thanks Tobi. Lets see how it goes.

[lpz456]

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