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)

 swahono October 25, 2010 03:34

Compressible epsilon blows up

2 Attachment(s)
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.

Best Regards,
Stefano

Attachment 5083

Attachment 5084

 FG_HSRM October 25, 2010 05:44

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

 swahono October 25, 2010 18:44

Hi Friedrich,

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```
It's getting pretty desperate here.

Best Regards,
Stefano

 FG_HSRM October 26, 2010 03:29

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.
until now I have no real new suggestion.

 swahono October 28, 2010 01:45

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

 hariya03 November 25, 2010 03:22

Reg - Y+

Dear Stefano,

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

I am not aware of it.

Thank you

Hari

 hariya03 November 25, 2010 03:23

Reg - Y+

Dear all,

Could you please tell me what y+ mean ?

I am not aware of it.

Thank you

Hari

 idrama November 25, 2010 03:46

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.

 hariya03 November 25, 2010 05:07

Reg - Y+

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?

Hari.

 swahono November 25, 2010 19:18

1 Attachment(s)
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

 idrama November 26, 2010 06:38

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.

 All times are GMT -4. The time now is 08:51.