the value of k and Epsilon is constantly increasing

saidc. · February 29, 2020, 10:04

Hi, i have a kEpsilon case. I created a geometry (like the attachment) and meshed with salome. Then i did import the mesh with "ideasUnvToFoam to a "pitzDaily" file which is the steady-state turbulent simpleFoam solver. I calculated the k and Epsilon values. When i run the case, the value of k and Epsilon is constantly increasing. What should i do, any suggestion?

I used;
Salome 9.3. (NETGEN 1D-2D-3D, 25k nodes)
OpenFoam-7

Hdraulic diameter= 5e-4
Reynolds= 25000
nu= 2e-7

controlDict

Code:

application     simpleFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         200;

deltaT          1e-4;

writeControl    adjustableRunTime;

writeInterval   1e-4;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

functions
{
    #includeFunc streamlines
}

Epsilon

Code:

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

internalField   uniform 1078.1;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 1078.1;
    }
    outlet
    {
        type            zeroGradient;
    }
    upperWall
    {
        type            epsilonWallFunction;
        value           uniform 1078.1;
    }
    lowerWall
    {
        type            epsilonWallFunction;
        value           uniform 1078.1;
    }
    frontAndBack
    {
        type            empty;
    }
}

k

Code:

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

internalField   uniform 0.375;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 0.375;
    }
    outlet
    {
        type            zeroGradient;
    }
    upperWall
    {
        type            kqRWallFunction;
        value           uniform 0.375;
    }
    lowerWall
    {
        type            kqRWallFunction;
        value           uniform 0.375;
    }
    frontAndBack
    {
        type            empty;
    }
}

U

Code:

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (-10 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    upperWall
    {
        type            noSlip;
    }

    lowerWall
    {
        type            noSlip;
    }

    frontAndBack
    {
        type            empty;
    }
}

p

Code:

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

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    upperWall
    {
        type            zeroGradient;
    }

    lowerWall
    {
        type            zeroGradient;
    }

    frontAndBack
    {
        type            empty;
    }
}

ERROR

Code:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
Build  : 7-1ff648926f77
Exec   : simpleFoam
Date   : Feb 29 2020
Time   : 17:34:15
Host   : "said"
PID    : 20416
I/O    : uncollated
Case   : /home/said/OpenFOAM/openfoam7/tutorials/denemeler/pitzDeneme1
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


SIMPLE: Convergence criteria found
        p: tolerance 0.01
        U: tolerance 0.001
        "(k|epsilon|omega|f|v2)": tolerance 0.001

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type RAS
Selecting RAS turbulence model kEpsilon
RAS
{
    RASModel        kEpsilon;
    turbulence      on;
    printCoeffs     on;
    Cmu             0.09;
    C1              1.44;
    C2              1.92;
    C3              0;
    sigmak          1;
    sigmaEps        1.3;
}

No MRF models present

No finite volume options present

Starting time loop

streamLine streamlines:
    automatic track length specified through number of sub cycles : 5

Time = 0.0001

smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 0.0350657, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 1, Final residual = 0.0327314, No Iterations 1
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.0728737, No Iterations 8
time step continuity errors : sum local = 0.000345373, global = -1.05665e-05, cumulative = -1.05665e-05
smoothSolver:  Solving for epsilon, Initial residual = 0.999999, Final residual = 0.0693174, No Iterations 5
bounding epsilon, min: -3.3678e+08 max: 1.2544e+10 average: 7.12552e+08
smoothSolver:  Solving for k, Initial residual = 1, Final residual = 7.96532e-07, No Iterations 1
bounding k, min: -24.3664 max: 702.772 average: 0.386658
ExecutionTime = 34.84 s  ClockTime = 35 s

streamLine streamlines write:
    seeded 0 particles
    Tracks:0
    Total samples:0
Time = 0.0002

smoothSolver:  Solving for Ux, Initial residual = 0.784733, Final residual = 0.0692963, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 0.815212, Final residual = 0.0663261, No Iterations 1
GAMG:  Solving for p, Initial residual = 0.982169, Final residual = 0.088101, No Iterations 1
time step continuity errors : sum local = 2.74335, global = 0.000341709, cumulative = 0.000331142
smoothSolver:  Solving for epsilon, Initial residual = 0.0125579, Final residual = 0.000768383, No Iterations 1
bounding epsilon, min: -1.67538e+07 max: 3.7362e+09 average: 3.5674e+07
smoothSolver:  Solving for k, Initial residual = 0.748306, Final residual = 0.0705025, No Iterations 20
bounding k, min: -12361.6 max: 88308.7 average: 0.682547
ExecutionTime = 54.37 s  ClockTime = 55 s

streamLine streamlines write:
    seeded 0 particles
    Tracks:0
    Total samples:0
Time = 0.0003

smoothSolver:  Solving for Ux, Initial residual = 0.809226, Final residual = 0.0099766, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 0.731181, Final residual = 0.0124558, No Iterations 1
GAMG:  Solving for p, Initial residual = 0.031464, Final residual = 0.00108104, No Iterations 2
time step continuity errors : sum local = 8.80402e+11, global = -4.86576e+08, cumulative = -4.86576e+08
smoothSolver:  Solving for epsilon, Initial residual = 1, Final residual = 0.0887976, No Iterations 3
bounding epsilon, min: -7.93897e+22 max: 6.5463e+23 average: 2.40067e+19
smoothSolver:  Solving for k, Initial residual = 0.999996, Final residual = 0.0472814, No Iterations 2
bounding k, min: -3.61155e+17 max: 2.11732e+18 average: 1.10185e+14
ExecutionTime = 70.58 s  ClockTime = 71 s

streamLine streamlines write:
    seeded 0 particles
    Tracks:0
    Total samples:0
Time = 0.0004

smoothSolver:  Solving for Ux, Initial residual = 0.830135, Final residual = 0.0220289, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.839851, Final residual = 0.0284726, No Iterations 2
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.0870518, No Iterations 4
time step continuity errors : sum local = 6.79143e+35, global = -1.49171e+34, cumulative = -1.49171e+34
smoothSolver:  Solving for epsilon, Initial residual = 1, Final residual = 0.0892089, No Iterations 2
bounding epsilon, min: -6.69723e+81 max: 1.97989e+98 average: 4.39393e+92
smoothSolver:  Solving for k, Initial residual = 1, Final residual = 1.35535e-13, No Iterations 1
bounding k, min: -4.75429e+59 max: 2.35257e+60 average: 6.087e+54
ExecutionTime = 88.71 s  ClockTime = 89 s

streamLine streamlines write:
    seeded 0 particles
    Tracks:0
    Total samples:0
Time = 0.0005

smoothSolver:  Solving for Ux, Initial residual = 0.99729, Final residual = 0.0409, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 0.997922, Final residual = 0.0452263, No Iterations 2
#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::GAMGSolver::scale(Foam::Field<double>&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, Foam::Field<double> const&, unsigned char) const at ??:?
#4  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
#5  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#6  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
#7  Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/opt/openfoam7/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
#8  Foam::fvMatrix<double>::solve() in "/opt/openfoam7/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
#9  ? in "/opt/openfoam7/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
#10  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11  ? in "/opt/openfoam7/platforms/linux64GccDPInt32Opt/bin/simpleFoam"

aow · February 29, 2020, 16:52

Hi Said,

Looking at your values for

and $\varepsilon$ , it seems like they might be a bit off. Taking $\ell_\mathrm{turb}=0.07D_\mathrm{Hydraulic}$ , $I=0.16Re^{-1/8}$ , $k=\frac{3}{2}(UI)^2$ , and $\varepsilon=C_\mu^{3/4}k^{3/2}/\ell_\mathrm{turb}$ , where $C_\mu=0.09$ , you get $\ell_\mathrm{turb}=0.000035~\mathrm{m}$ ,

, $k=0.3054~\mathrm{m^2/s^2}$ , $\varepsilon=792.3~\mathrm{m^2/s^3}$ based on the Wikipedia TKE webpage.

However, these relatively small differences are probably not what's causing your simulation to diverge. I would try one or both of two things:
1) reduce the relaxation parameters in the fvSolution file for your case, 2) gradually increase the velocity/turbulence instead of starting from rest and immediately trying to jump to extremely high values for a very turbulent flow condition. The best way to do this is probably to ramp up the velocity using the timeVaryingMappedFixedValue boundary condition for your inlet boundary velocity. Reading the header text at the beginning of the source code header file (see https://github.com/OpenFOAM/OpenFOAM...FvPatchField.H) helps to explain what the setup for this BC is like and searching through the forums you'll find references to it as well.

Also, have you used the checkMesh utility to check if there are any poor-quality cells in your mesh? If not, I would do this first to verify that there are not any severely non-orthogonal or skewed cells.

saidc. · March 1, 2020, 08:22

Dear aow,

Thanks for your reply. I did what you said but there is nothing change. And also i change my geometry, mesh, deltaT, Reynolds number but the results still same, k and epsilon is constantly increasing. I did try a lot of things but my case still doesn't working. Any suggestion, what should i do?

aow · March 1, 2020, 18:21

Sorry that didn't help at all. Did you use blockMesh and snappyHexMesh or are you importing the mesh entirely from Salome? I've never used Salome so I'm not familiar with what sort of results it produces. Anyway, if you did use blockMesh/snappyHexMesh, could you post your dictionary files for them?

Also, have you tried playing around with the time step? Maybe 1e-4 is initially too large. You could try starting with 1e-5 and 1e-6 instead to see if there is any improvement. Your tolerances for p, U, k, etc... might be a bit large too so maybe try reducing them a bit as well. If you could show your fvSolution and fvSchemes files that might help identify the problem too.

saidc. · March 2, 2020, 00:32

Quote:

Originally Posted by aow

Sorry that didn't help at all. Did you use blockMesh and snappyHexMesh or are you importing the mesh entirely from Salome? I've never used Salome so I'm not familiar with what sort of results it produces. Anyway, if you did use blockMesh/snappyHexMesh, could you post your dictionary files for them?

Also, have you tried playing around with the time step? Maybe 1e-4 is initially too large. You could try starting with 1e-5 and 1e-6 instead to see if there is any improvement. Your tolerances for p, U, k, etc... might be a bit large too so maybe try reducing them a bit as well. If you could show your fvSolution and fvSchemes files that might help identify the problem too.

Dear aow,

Thanks again. I did import the mesh with salome, I didn't use blockMesh or snappyHexMesh. Do you think this could be a problem? And I tried playing around with deltaT but it diverging same time (like 0.001 sec). Here my other case files:

fvSolution

Code:

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-6;
        relTol          1e-3; //0.1
        smoother        GaussSeidel;
    }

    "(U|k|epsilon|omega|f|v2)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          1e-3; //0.1
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent      yes;

    residualControl
    {
        p               1e-2;
        U               1e-3;
        "(k|epsilon|omega|f|v2)" 1e-3;
    }
}

relaxationFactors
{
    equations
    {
        U               0.8; // 0.9 is more stable but 0.95 more convergent
        ".*"            0.8; // 0.9 is more stable but 0.95 more convergent
    }
}

fvSchemes

Code:

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);
    div(phi,k)      bounded Gauss limitedLinear 1;
    div(phi,epsilon) bounded Gauss limitedLinear 1;
    div(phi,omega)  bounded Gauss limitedLinear 1;
    div(phi,v2)     bounded Gauss limitedLinear 1;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
    div(nonlinearStress) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method meshWave;
}

aow · April 23, 2020, 22:18

Hi Saidc,

Sorry for the slow reply... I've been busy with work and have not paid attention to CFD Online lately.

Hmm, well there doesn't seem to be anything obviously wrong with your solution and schemes files to me. What's a bit weird to me is that you only have solution options specified for p and "(U|k|epsilon|omega|f|v2)", but I'm not used to running steady state solutions in OpenFOAM and after examining some of the steady state tutorials it looks like you're not doing anything necessarily wrong there. Another thing that I noticed is that it seems like your residual control values are rather large at 1e-2 and 1e-3. It looks like you copied the pitzDaily tutorial based on comparing the scripts you've posted to the tutorial files. Perhaps reduce the residual control values and try again to force the simulation to go through more iterations and reach a more converged state. For example, the incompressible/simpleFoam/airFoil2D tutorial uses 1e-5 as the residual control value for all field variables being solved for in the simulation.

Also, have you tried recording and plotting solution residuals for the variables being solved for in your simulation? You should verify they're actually converging to a steady value in each time step if it's an unsteady model, or, that it reaches a single converged value in a steady state model.

Otherwise, if your solution and schemes setup are okay and so are your initial/boundary conditions, then the problem most likely lies with your mesh. Is it possible for you to post an image of it or at least of the domain? I'm wondering if you have too coarse of a mesh somewhere or high skew and/or non-orthogonality of some cells. As I mentioned before, you can check this using the checkMesh utility; however, even if all checks pass this is not a fool-proof approach of checking your mesh. You may need to refine it in regions of your model, especially where large gradients in pressure, velocity, etc... are expected to occur, but be careful how you do this. If you use snappyHexMesh, it will refine very abruptly with large jumps in cell size. It is better to use uniform grading like blockMesh can do. If you're using Salome, you may need to figure out a way to do this to get smooth cell size variation when refining the mesh.

aow · April 23, 2020, 22:35

One more thing I thought I should mention is that your k and epsilon values have already hit "bounding" values on your first solver iteration at time = 0.0001 s, which tells me something is fundamentally wrong with your initial setup. Maybe your outlet boundary conditions for k and epsilon are causing part of this problem? Perhaps you should try changing them to the following example, which should handle backflow better than the zeroGradient BC if that is happening in your simulation. The other BCs for k and epsilon look okay to me, it's just that your initial values are very large compared to what I've seen doing water wave impact simulations and it might be part of the problem. I would verify that they are indeed the values you should be using once more to be sure.

Code:

outlet
{
    type            inletOutlet;
    inletValue      $internalField;
    value           $internalField;
}

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
SimpleFoam k and epsilon bounded	nedved	OpenFOAM Running, Solving & CFD	16	March 4, 2017 08:30
simpleFoam error - "Floating point exception"	mbcx4jc2	OpenFOAM Running, Solving & CFD	12	August 4, 2015 02:20
Calculation of k and epsilon freezes	Nigirim	OpenFOAM Running, Solving & CFD	1	November 14, 2012 07:52
epsilon and K blowing up.	sivakumar	OpenFOAM Running, Solving & CFD	1	October 25, 2012 04:50
SimpleFoam k and epsilon bounded	nedved	OpenFOAM Running, Solving & CFD	1	November 25, 2008 20:21

February 29, 2020, 16:52		#2
aow Member Andrew O. Winter Join Date: Aug 2015 Location: Seattle, WA, USA Posts: 78 Rep Power: 10	Hi Said, Looking at your values for and $\varepsilon$ , it seems like they might be a bit off. Taking $\ell_\mathrm{turb}=0.07D_\mathrm{Hydraulic}$ , $I=0.16Re^{-1/8}$ , $k=\frac{3}{2}(UI)^2$ , and $\varepsilon=C_\mu^{3/4}k^{3/2}/\ell_\mathrm{turb}$ , where $C_\mu=0.09$ , you get $\ell_\mathrm{turb}=0.000035~\mathrm{m}$ , , $k=0.3054~\mathrm{m^2/s^2}$ , $\varepsilon=792.3~\mathrm{m^2/s^3}$ based on the Wikipedia TKE webpage. However, these relatively small differences are probably not what's causing your simulation to diverge. I would try one or both of two things: 1) reduce the relaxation parameters in the fvSolution file for your case, 2) gradually increase the velocity/turbulence instead of starting from rest and immediately trying to jump to extremely high values for a very turbulent flow condition. The best way to do this is probably to ramp up the velocity using the timeVaryingMappedFixedValue boundary condition for your inlet boundary velocity. Reading the header text at the beginning of the source code header file (see https://github.com/OpenFOAM/OpenFOAM...FvPatchField.H) helps to explain what the setup for this BC is like and searching through the forums you'll find references to it as well. Also, have you used the checkMesh utility to check if there are any poor-quality cells in your mesh? If not, I would do this first to verify that there are not any severely non-orthogonal or skewed cells. saidc. and shuji_type97 like this. Last edited by aow; February 29, 2020 at 22:04.

March 1, 2020, 08:22		#3
saidc. Member saidc Join Date: Feb 2020 Location: nowhere Posts: 62 Rep Power: 6	Dear aow, Thanks for your reply. I did what you said but there is nothing change. And also i change my geometry, mesh, deltaT, Reynolds number but the results still same, k and epsilon is constantly increasing. I did try a lot of things but my case still doesn't working. Any suggestion, what should i do?

March 1, 2020, 18:21		#4
aow Member Andrew O. Winter Join Date: Aug 2015 Location: Seattle, WA, USA Posts: 78 Rep Power: 10	Sorry that didn't help at all. Did you use blockMesh and snappyHexMesh or are you importing the mesh entirely from Salome? I've never used Salome so I'm not familiar with what sort of results it produces. Anyway, if you did use blockMesh/snappyHexMesh, could you post your dictionary files for them? Also, have you tried playing around with the time step? Maybe 1e-4 is initially too large. You could try starting with 1e-5 and 1e-6 instead to see if there is any improvement. Your tolerances for p, U, k, etc... might be a bit large too so maybe try reducing them a bit as well. If you could show your fvSolution and fvSchemes files that might help identify the problem too. saidc. likes this.

April 23, 2020, 22:18		#6
aow Member Andrew O. Winter Join Date: Aug 2015 Location: Seattle, WA, USA Posts: 78 Rep Power: 10	Hi Saidc, Sorry for the slow reply... I've been busy with work and have not paid attention to CFD Online lately. Hmm, well there doesn't seem to be anything obviously wrong with your solution and schemes files to me. What's a bit weird to me is that you only have solution options specified for p and "(U\|k\|epsilon\|omega\|f\|v2)", but I'm not used to running steady state solutions in OpenFOAM and after examining some of the steady state tutorials it looks like you're not doing anything necessarily wrong there. Another thing that I noticed is that it seems like your residual control values are rather large at 1e-2 and 1e-3. It looks like you copied the pitzDaily tutorial based on comparing the scripts you've posted to the tutorial files. Perhaps reduce the residual control values and try again to force the simulation to go through more iterations and reach a more converged state. For example, the incompressible/simpleFoam/airFoil2D tutorial uses 1e-5 as the residual control value for all field variables being solved for in the simulation. Also, have you tried recording and plotting solution residuals for the variables being solved for in your simulation? You should verify they're actually converging to a steady value in each time step if it's an unsteady model, or, that it reaches a single converged value in a steady state model. Otherwise, if your solution and schemes setup are okay and so are your initial/boundary conditions, then the problem most likely lies with your mesh. Is it possible for you to post an image of it or at least of the domain? I'm wondering if you have too coarse of a mesh somewhere or high skew and/or non-orthogonality of some cells. As I mentioned before, you can check this using the checkMesh utility; however, even if all checks pass this is not a fool-proof approach of checking your mesh. You may need to refine it in regions of your model, especially where large gradients in pressure, velocity, etc... are expected to occur, but be careful how you do this. If you use snappyHexMesh, it will refine very abruptly with large jumps in cell size. It is better to use uniform grading like blockMesh can do. If you're using Salome, you may need to figure out a way to do this to get smooth cell size variation when refining the mesh. saidc. likes this.

April 23, 2020, 22:35		#7
aow Member Andrew O. Winter Join Date: Aug 2015 Location: Seattle, WA, USA Posts: 78 Rep Power: 10	One more thing I thought I should mention is that your k and epsilon values have already hit "bounding" values on your first solver iteration at time = 0.0001 s, which tells me something is fundamentally wrong with your initial setup. Maybe your outlet boundary conditions for k and epsilon are causing part of this problem? Perhaps you should try changing them to the following example, which should handle backflow better than the zeroGradient BC if that is happening in your simulation. The other BCs for k and epsilon look okay to me, it's just that your initial values are very large compared to what I've seen doing water wave impact simulations and it might be part of the problem. I would verify that they are indeed the values you should be using once more to be sure. Code: outlet { type inletOutlet; inletValue $internalField; value $internalField; } MaySea and saidc. like this.