# Compressible epsilon blows up

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

 October 25, 2010, 03:34 Compressible epsilon blows up #1 Member   Stefano Wahono Join Date: Aug 2010 Location: Melbourne, Australia Posts: 42 Rep Power: 8 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::T(scalar f, scalar T0, scalar (specieThermo::*F)(const scalar) const, scalar (specieThermo::*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 > > > >::calculate() in "/home/stefano/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/libbasicThermophysicalModels.so" [7] #3 Foam::hPsiThermo > > > >::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 run6-topAir-bottomAir-Streamlines2Ts.jpg parallelPlates.tar.gz

 October 25, 2010, 05:44 #2 Member   Join Date: Dec 2009 Posts: 36 Rep Power: 8 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

October 25, 2010, 18:44
#3
Member

Stefano Wahono
Join Date: Aug 2010
Location: Melbourne, Australia
Posts: 42
Rep Power: 8
Hi Friedrich,

Quote:
 Originally Posted by FG_HSRM 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 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]     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(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

 October 26, 2010, 03:29 #4 Member   Join Date: Dec 2009 Posts: 36 Rep Power: 8 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.

October 28, 2010, 01:45
#5
Member

Stefano Wahono
Join Date: Aug 2010
Location: Melbourne, Australia
Posts: 42
Rep Power: 8
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 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 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(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;
}

{
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

 November 25, 2010, 03:22 Reg - Y+ #6 New Member   Hari Krishnan Join Date: Mar 2009 Location: chennai, Tamil nadu, india Posts: 27 Rep Power: 9 Dear Stefano, Could you please tell me what you mean by y+? I am not aware of it. Thank you Hari

 November 25, 2010, 03:23 Reg - Y+ #7 New Member   Hari Krishnan Join Date: Mar 2009 Location: chennai, Tamil nadu, india Posts: 27 Rep Power: 9 Dear all, Could you please tell me what y+ mean ? I am not aware of it. Thank you Hari

 November 25, 2010, 03:46 #8 Senior Member   Claus Meister Join Date: Aug 2009 Location: Wiesbaden, Germany Posts: 241 Rep Power: 10 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.

 November 25, 2010, 05:07 Reg - Y+ #9 New Member   Hari Krishnan Join Date: Mar 2009 Location: chennai, Tamil nadu, india Posts: 27 Rep Power: 9 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.

November 25, 2010, 19:18
#10
Member

Stefano Wahono
Join Date: Aug 2010
Location: Melbourne, Australia
Posts: 42
Rep Power: 8
Quote:
 Originally Posted by idrama 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 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 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
Attached Images
 mesh.jpg (93.4 KB, 32 views)

 November 26, 2010, 06:38 #11 Senior Member   Claus Meister Join Date: Aug 2009 Location: Wiesbaden, Germany Posts: 241 Rep Power: 10 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.

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post nedved OpenFOAM Running, Solving & CFD 15 December 2, 2016 03:40 nedved OpenFOAM Running, Solving & CFD 1 November 25, 2008 21:21 Jinwon Main CFD Forum 6 November 23, 2007 22:07 luca OpenFOAM Running, Solving & CFD 4 July 3, 2006 06:24 Fernando Velasco Hurtado Main CFD Forum 3 January 7, 2000 17:51

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