CFD Online Discussion Forums - Compressible epsilon blows up

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

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

- - Compressible epsilon blows up (https://www.cfd-online.com/Forums/openfoam/81354-compressible-epsilon-blows-up.html)

Compressible epsilon blows up

Hi All,

I'm simulating a conjugate heat transfer across a flat plate in a rectangular channel.
I modified the multiRegionHeater tutorial in OF 1.7.0, and used k-epsilon turbulence model.

The Case geometry setup is shown in the attached PNG file.
Basically, the top and bottom air is turbulent (Re ~ 10^6).
The top air is hot at 746K, while the bottom air is cold at 288K.
The two flow regions are separated by a solid plate.

I ran the case with realizableKE model, and the simulation blew up after 38 iterations.
Here is the last few iterations before it blowing up:

Code:

Time = 36



Solving for fluid region bottomAir

DILUPBiCG:  Solving for Ux, Initial residual = 0.3204164, Final residual = 0.0049385518, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.43623786, Final residual = 0.01306048, No Iterations 1

DILUPBiCG:  Solving for Uz, Initial residual = 0.46554744, Final residual = 0.015560302, No Iterations 1

DILUPBiCG:  Solving for h, Initial residual = 0.14392941, Final residual = 0.010392872, No Iterations 1

Min/max T:119.30898 433.59925

GAMG:  Solving for p_rgh, Initial residual = 0.049386189, Final residual = 0.00042936621, No Iterations 4

time step continuity errors : sum local = 68.178516, global = -12.111985, cumulative = -12.049972

Min/max rho:0.2 2

DILUPBiCG:  Solving for epsilon, Initial residual = 0.0060067685, Final residual = 0.000382735, No Iterations 1

bounding epsilon, min: -810751.63 max: 1.6587192e+11 average: 2266471.6

DILUPBiCG:  Solving for k, Initial residual = 0.60163212, Final residual = 0.020304535, No Iterations 1

bounding k, min: -196.25661 max: 1.3913795e+12 average: 28036037

maxResidual: 0.46554744  convergence criterion: 0



Solving for fluid region topAir

DILUPBiCG:  Solving for Ux, Initial residual = 0.010704973, Final residual = 0.000326993, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.0189083, Final residual = 0.00060774355, No Iterations 1

DILUPBiCG:  Solving for Uz, Initial residual = 0.018775159, Final residual = 0.00059849536, No Iterations 1

DILUPBiCG:  Solving for h, Initial residual = 0.0093179606, Final residual = 9.2234893e-05, No Iterations 1

Min/max T:290.60892 746.00013

GAMG:  Solving for p_rgh, Initial residual = 0.0037615237, Final residual = 1.1872936e-05, No Iterations 3

time step continuity errors : sum local = 4.7313027e-07, global = -2.6202936e-08, cumulative = -12.049972

Min/max rho:0.47318397 1.2147625

DILUPBiCG:  Solving for epsilon, Initial residual = 0.00086839005, Final residual = 1.9941133e-05, No Iterations 1

DILUPBiCG:  Solving for k, Initial residual = 0.010316428, Final residual = 0.00084437086, No Iterations 1

maxResidual: 0.0189083  convergence criterion: 0



Solving for solid region solidWall

DICPCG:  Solving for T, Initial residual = 0.025633671, Final residual = 0.0021160135, No Iterations 1

Min/max T:287.64755 296.48869

maxResidual: 0.025633671  convergence criterion: 0

ExecutionTime = 200.4 s  ClockTime = 201 s



Time = 37





Solving for fluid region bottomAir

DILUPBiCG:  Solving for Ux, Initial residual = 0.50298705, Final residual = 0.012268166, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.46978779, Final residual = 0.013510551, No Iterations 1

DILUPBiCG:  Solving for Uz, Initial residual = 0.45493707, Final residual = 0.011565384, No Iterations 1

DILUPBiCG:  Solving for h, Initial residual = 0.025038369, Final residual = 0.0019166388, No Iterations 1

Min/max T:-1557.9293 422.26051

GAMG:  Solving for p_rgh, Initial residual = 0.19809226, Final residual = 0.0019424248, No Iterations 4

time step continuity errors : sum local = 2.9035403e+10, global = 5.090356e+09, cumulative = 5.0903559e+09

Min/max rho:0.2 2

DILUPBiCG:  Solving for epsilon, Initial residual = 0.021184386, Final residual = 5.8090526e-05, No Iterations 3

bounding epsilon, min: -3578671.8 max: 5.937534e+11 average: 10175646

DILUPBiCG:  Solving for k, Initial residual = 0.99999971, Final residual = 0.050065216, No Iterations 2

maxResidual: 0.50298705  convergence criterion: 0



Solving for fluid region topAir

DILUPBiCG:  Solving for Ux, Initial residual = 0.010419636, Final residual = 0.00031881238, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.018670427, Final residual = 0.00059330634, No Iterations 1

DILUPBiCG:  Solving for Uz, Initial residual = 0.018306335, Final residual = 0.0005741435, No Iterations 1

DILUPBiCG:  Solving for h, Initial residual = 0.0090171776, Final residual = 8.9433187e-05, No Iterations 1

Min/max T:290.65192 746.00013

GAMG:  Solving for p_rgh, Initial residual = 0.003663336, Final residual = 1.8532424e-05, No Iterations 3

time step continuity errors : sum local = 7.4935637e-07, global = -2.1118708e-08, cumulative = 5.0903559e+09

Min/max rho:0.47318398 1.214584

DILUPBiCG:  Solving for epsilon, Initial residual = 0.00083188985, Final residual = 1.8973471e-05, No Iterations 1

DILUPBiCG:  Solving for k, Initial residual = 0.0098047206, Final residual = 0.0008637932, No Iterations 1

maxResidual: 0.018670427  convergence criterion: 0



Solving for solid region solidWall

DICPCG:  Solving for T, Initial residual = 0.025238773, Final residual = 0.0020939892, No Iterations 1

Min/max T:287.7893 295.63824

maxResidual: 0.025238773  convergence criterion: 0

ExecutionTime = 204.98 s  ClockTime = 206 s



Time = 38





Solving for fluid region bottomAir

DILUPBiCG:  Solving for Ux, Initial residual = 0.3018213, Final residual = 0.004181555, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.39316677, Final residual = 0.0066244412, No Iterations 1

DILUPBiCG:  Solving for Uz, Initial residual = 0.4739378, Final residual = 0.014966079, No Iterations 1

DILUPBiCG:  Solving for h, Initial residual = 0.01415212, Final residual = 0.0010154766, No Iterations 1

[7] 

[7] 

[7] --> FOAM FATAL ERROR: 

[7] Maximum number of iterations exceeded

[7] 

[7]     From function specieThermo<thermo>::T(scalar f, scalar T0, scalar (specieThermo<thermo>::*F)(const scalar) const, scalar (specieThermo<thermo>::*dFdT)(const scalar) const) const

[7]     in file /home/stefano/OpenFOAM/OpenFOAM-1.7.0/src/thermophysicalModels/specie/lnInclude/specieThermoI.H at line 67.

[7] 

FOAM parallel run aborting

[7] 

[7] #0  Foam::error::printStack(Foam::Ostream&) in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/libOpenFOAM.so"

[7] #1  Foam::error::abort() in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/libOpenFOAM.so"

[7] #2  Foam::hPsiThermo<Foam::pureMixture<Foam::constTransport<Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> > > > >::calculate() in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/libbasicThermophysicalModels.so"

[7] #3  Foam::hPsiThermo<Foam::pureMixture<Foam::constTransport<Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> > > > >::correct() in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/libbasicThermophysicalModels.so"

[7] #4  main in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/applications/bin/linux64GccDPOpt/chtMultiRegionSimpleFoam"

[7] #5  __libc_start_main in "/lib64/libc.so.6"

[7] #6  Foam::regIOobject::writeObject(Foam::IOstream::streamFormat, Foam::IOstream::versionNumber, Foam::IOstream::compressionType) const in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/applications/bin/linux64GccDPOpt/chtMultiRegionSimpleFoam"

--------------------------------------------------------------------------

MPI_ABORT was invoked on rank 7 in communicator MPI_COMM_WORLD 

with errorcode 1.



NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.

You may or may not see output from other processes, depending on

exactly when Open MPI kills them.

--------------------------------------------------------------------------

--------------------------------------------------------------------------

mpirun has exited due to process rank 7 with PID 5379 on

node barracuda exiting without calling "finalize". This may

have caused other processes in the application to be

terminated by signals sent by mpirun (as reported here).

--------------------------------------------------------------------------

Please if anyone can help me....It is urgent.
I have attached the /system and /constant files.

The case ran well with KwSST model.

Thank you in advance.

Best Regards,
Stefano

Attachment 5083

Attachment 5084

Hi Stefano

did you try to switch to another div scheme for epsilon?
You are using different Schemes in top and bottom air.

My suggestion is, try in the limtedUpwind scheme, that you are using in the bottom area or a limitedLinear scheme. You also can dry to reduce the relaxtionFactors for epsilon to 0.5.

Regards

Hi Friedrich,

Thank you for your reply.

Quote:

Originally Posted by FG_HSRM (Post 280610)

did you try to switch to another div scheme for epsilon?
You are using different Schemes in top and bottom air.

Yes, I changed the div scheme on bottom air because only the bottom air got unstable.

Quote:

Originally Posted by FG_HSRM (Post 280610)

My suggestion is, try in the limtedUpwind scheme, that you are using in the bottom area or a limitedLinear scheme. You also can dry to reduce the relaxtionFactors for epsilon to 0.5.

I have tried the limitedLinear scheme; however, the simulation did not run.
It returns the following error:

Code:

[1] --> FOAM FATAL IO ERROR: 

[1] attempt to read beyond EOF

[1] 

[1] file: /home/stefano/OpenFOAM/stefano-1.7.0/run/parallelPlates-run6A/processor1/../system/bottomAir/fvSchemes::divSchemes::div(phi,U) at line 30.

[1] 

[1]     From function ITstream::read(token& t)

[1]     in file db/IOstreams/Tstreams/ITstream.C at line 83.

Any of the "limited" schemes returned the above error.

I found the following schemes ran:

Code:

divSchemes

{

    default         none;

    div(phi,U)      Gauss linearUpwindV cellLimited Gauss linear 1;

    div(phi,h)      Gauss linear;

    div(phi,k)      Gauss linear;

    div(phi,epsilon) Gauss linear;

    div(phi,R)      Gauss linear;

    div(R)          Gauss linear;

    div((muEff*dev2(grad(U).T()))) Gauss linear;

    div(U,p)        Gauss linearUpwind cellLimited Gauss linear 1;

I have also reduced the epsilon and h relaxation factor to 0.3.

After about 50 iterations, the simulation blew up again.
Epsilon reached -4839598, and k reached -3402.
Temperature and rho on bottom air also reached unrealistic extremes.

The weird thing is that the top air and solid seems okay.

Here is the extract of the last iteration:

Code:

Time = 53



Solving for fluid region bottomAir

DILUPBiCG:  Solving for Ux, Initial residual = 0.19770931, Final residual = 0.0020336792, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.41642211, Final residual = 0.005267386, No Iterations 1

DILUPBiCG:  Solving for Uz, Initial residual = 0.37520706, Final residual = 0.0037334979, No Iterations 1

DILUPBiCG:  Solving for h, Initial residual = 0.1065923, Final residual = 0.0032230233, No Iterations 1

Min/max T:-48.660131 609.49704

GAMG:  Solving for p_rgh, Initial residual = 0.25635642, Final residual = 0.0024364087, No Iterations 4

time step continuity errors : sum local = 1.5766647, global = -0.14254934, cumulative = 0.81054565

Min/max rho:0.2 2

DILUPBiCG:  Solving for epsilon, Initial residual = 0.044754031, Final residual = 0.00022822142, No Iterations 1

bounding epsilon, min: -4893538.9 max: 54046780 average: 16582.176

DILUPBiCG:  Solving for k, Initial residual = 0.13676716, Final residual = 0.0030502085, No Iterations 2

bounding k, min: -3402.294 max: 53095.119 average: 472.50111

maxResidual: 0.41642211  convergence criterion: 0



Solving for fluid region topAir

DILUPBiCG:  Solving for Ux, Initial residual = 0.0061144655, Final residual = 9.036451e-05, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.013482903, Final residual = 0.000134842, No Iterations 1

DILUPBiCG:  Solving for Uz, Initial residual = 0.016297898, Final residual = 0.00014666149, No Iterations 1

DILUPBiCG:  Solving for h, Initial residual = 0.0048919717, Final residual = 3.7624836e-05, No Iterations 1

Min/max T:291.29734 746.00059

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

time step continuity errors : sum local = 9.3353508e-07, global = -8.8156252e-08, cumulative = 0.81054556

Min/max rho:0.47318442 1.2120324

DILUPBiCG:  Solving for epsilon, Initial residual = 0.00044237679, Final residual = 2.0286949e-06, No Iterations 1

DILUPBiCG:  Solving for k, Initial residual = 0.0046153281, Final residual = 0.00041521771, No Iterations 1

maxResidual: 0.016297898  convergence criterion: 0



Solving for solid region solidWall

DICPCG:  Solving for T, Initial residual = 0.020776695, Final residual = 0.0017699078, No Iterations 1

Min/max T:278.99619 298.02374

maxResidual: 0.020776695  convergence criterion: 0

ExecutionTime = 301.91 s  ClockTime = 304 s

Please if anyone can help.
It's getting pretty desperate here.

Best Regards,
Stefano

Hi Stefano
how are your result with upwind?
I know it's just 1st order and I also prefer higher order schemes, but it is stable.

Can you please post you fvSchemes with the limited schemes, that is not working.
I would like to take a look and understand why not.
Maybe somebody with more experience can help you,
until now I have no real new suggestion.

Hi Friedrich,

My apologies for delay in responding to your post.
The upwind result was not satisfactory when compared to Fluent.
That is the reason I tried the 2nd order scheme.

Quote:

Originally Posted by FG_HSRM (Post 280734)

how are your result with upwind?
I know it's just 1st order and I also prefer higher order schemes, but it is stable.

Quote:

Originally Posted by FG_HSRM (Post 280734)

Can you please post you fvSchemes with the limited schemes, that is not working.
I would like to take a look and understand why not.

Unfortunately, I have overwritten my previous fvSchemes.
This is my fvSchemes that works:

Code:

divSchemes

{

    default         none;

    div(phi,U)      Gauss linearUpwindV cellMDLimited Gauss linear 1;

    div(phi,h)      Gauss linearUpwind cellMDLimited Gauss linear 1;

    div(phi,k)      Gauss linearUpwind cellMDLimited Gauss linear 1;

    div(phi,epsilon) Gauss linearUpwind cellMDLimited Gauss linear 1;

    div(phi,R)      Gauss linearUpwind cellMDLimited Gauss linear 1;

    div(R)          Gauss linearUpwind cellMDLimited Gauss linear 1;

    div((muEff*dev2(grad(U).T()))) Gauss linear;

    div(U,p)        Gauss linearUpwind cellMDLimited Gauss linear 1;

}



laplacianSchemes

{

    default         none;

    laplacian(muEff,U) Gauss linear limited 0.333;

    laplacian((rho*(1|A(U))),p_rgh) Gauss linear limited 0.333;

    laplacian(alphaEff,h) Gauss linear limited 0.333;

    laplacian(DkEff,k) Gauss linear limited 0.333;

    laplacian(DepsilonEff,epsilon) Gauss linear limited 0.333;

    laplacian(DREff,R) Gauss linear limited 0.333;

}



interpolationSchemes

{

    default         linear;

}



snGradSchemes

{

    default         limited 0.333;

}

I also found that epsilon is highly sensitive to y+.
The y+ for this model is up to 900. That is probably the reason of it blew up.

I also had to put all under-relaxation factors to 0.1 for it to work.
Even then, I had poor pressure convergence even after 4000 iterations!

With OF, it is hard to distinguish if the problem is caused by improper BC setup, poor mesh or solver setup.
I tried to eliminate the first two factors by using standard BC as everyone recommended and structured grid.

Another thing, for compressible solver with high temperature, is it better to use fixedValue for velocity inlet or flowRateInletVelocity?
The flowRateInletVelocity always gave the wrong inlet velocity when I tried it.

Thank you very much for your help.
I hope we can continue this discussion.

Best Regards,
Stefano

Dear Stefano,

Could you please tell me what you mean by y+?

I am not aware of it.

Thank you

Hari

Dear all,

Could you please tell me what y+ mean ?

I am not aware of it.

Thank you

Hari

Hallo hariya03,

y+ is a dimensionless factor which is used by wall function to compute the velocities at the wall, i.e. the velocities are not obtained form a empirical mathematical model. See Versteeg.

swahono,

you set relaxation factor to 0.1: could you please tell me how your maxCo number is bounded. Moreover, could you post you fvSolution and fvSchemes and a cross section plot of your mesh.

Thank you Mr.Claus,

I hope the Velocity at the wall will be zero. So May I know how it computes the velocity?

Did in "versteeg " the same denotion were used sir?

Thank you for your reply.

Hari.

Quote:

Originally Posted by idrama (Post 284727)

you set relaxation factor to 0.1: could you please tell me how your maxCo number is bounded.

I used chtMultiRegionSimpleFoam (steady RANS solver) for my case above. So there is no URF :)

Quote:

Originally Posted by idrama (Post 284727)

Moreover, could you post you fvSolution and fvSchemes and a cross section plot of your mesh.

I have posted my fvSchemes in the post above. My fvSolution is:

Code:

solvers

{

    rho

    {

        solver          PCG

        preconditioner  DIC;

        tolerance       1e-10;

        relTol          0;

    }



    p_rgh

    {

        solver           GAMG;

        tolerance        1e-10;

        relTol           0;

        smoother         GaussSeidel;

        cacheAgglomeration true;

        nCellsInCoarsestLevel 20;

        agglomerator     faceAreaPair;

        mergeLevels      2;

        nSweeps  2;

        nPreSweeps  1;

        nPostSweeps  2;

    }



    "(U|h|k|epsilon)"

    {

        solver           PBiCG;

        preconditioner   DILU;

        tolerance        1e-09;

        relTol           0;

    }

}



SIMPLE

{

    nNonOrthogonalCorrectors 0;

    pRefCell                0;

    pRefValue               101325;

    rhoMin      rhoMin [1 -3 0 0 0] 0.2;

    rhoMax      rhoMax [1 -3 0 0 0] 2;

}



relaxationFactors

{

    rho         0.3;

    p_rgh       0.1;

    U           0.1;

    h           0.1;

    nuTilda     0.1;

    k           0.1;

    epsilon     0.1;

}

You can increase you relaxation factor to speed up convergence at later iterations.

My Mesh is attached. I just used blockMesh.

Quote:

Originally Posted by hariya03 (Post 284738)

I hope the Velocity at the wall will be zero. So May I know how it computes the velocity?

Did in "versteeg " the same denotion were used sir?
Hari.

- Velocity at the wall is zero (you set it up as BC in 0/U).
- What Hari referred to was the use of "wall function" to compute the U, k and epsilon profiles near the wall (inside the boundary layer). This is based on the friction velocity , u*. Simply put, if your yPlus is small, means your first grid point is well within the boundary layer.

Most wall function works between yPlus of 30 - 100.

In SimpleFoam, pisoFoam, rhoSimpleFoam and rhoPisoFoam you can compute yPlus by using
>> yPlusRAS

This utility does not work with chtMultiRegionFoam without modification to the source code. I used Fluent to checked my yPlus in this case.
If I have time I will try to modify the yPlusRAS to work with chtMultiRegionFoam.

Best regards,
Stefano

Unfortunately, I didn't post you maxCo number, but I would recommended to reduce it a bit. But more dramatically, in you SIMPLE entries the keyword
nCorrectors does not appear. You should write, e.g. before nOrthogonalCorrectors,

nCorrectors 5;

Start with five is fine and then increasing or decreasing it accordingly. In your scheams directors you use limited 0.333 but nOrthogonalCorrestors is here set to zero. Do you have non-Ortho. cells? If yes then set

nOrthogonalCorrectors 3;

as starting point. In the run of the simulation you can adjust them accordingly. For now, that's all. Just post if stuck.