
[Sponsors] 
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 kepsilon 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.2234893e05, No Iterations 1 Min/max T:290.60892 746.00013 GAMG: Solving for p_rgh, Initial residual = 0.0037615237, Final residual = 1.1872936e05, No Iterations 3 time step continuity errors : sum local = 4.7313027e07, global = 2.6202936e08, cumulative = 12.049972 Min/max rho:0.47318397 1.2147625 DILUPBiCG: Solving for epsilon, Initial residual = 0.00086839005, Final residual = 1.9941133e05, 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.8090526e05, 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.9433187e05, No Iterations 1 Min/max T:290.65192 746.00013 GAMG: Solving for p_rgh, Initial residual = 0.003663336, Final residual = 1.8532424e05, No Iterations 3 time step continuity errors : sum local = 7.4935637e07, global = 2.1118708e08, cumulative = 5.0903559e+09 Min/max rho:0.47318398 1.214584 DILUPBiCG: Solving for epsilon, Initial residual = 0.00083188985, Final residual = 1.8973471e05, 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/OpenFOAM1.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/OpenFOAM1.7.0/lib/linux64GccDPOpt/libOpenFOAM.so" [7] #1 Foam::error::abort() in "/home/stefano/OpenFOAM/OpenFOAM1.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/OpenFOAM1.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/OpenFOAM1.7.0/lib/linux64GccDPOpt/libbasicThermophysicalModels.so" [7] #4 main in "/home/stefano/OpenFOAM/OpenFOAM1.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/OpenFOAM1.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).  I have attached the /system and /constant files. The case ran well with KwSST model. Thank you in advance. Best Regards, Stefano run6topAirbottomAirStreamlines2Ts.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,
Thank you for your reply. Quote:
Quote:
It returns the following error: Code:
[1] > FOAM FATAL IO ERROR: [1] attempt to read beyond EOF [1] [1] file: /home/stefano/OpenFOAM/stefano1.7.0/run/parallelPlatesrun6A/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. 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; 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.036451e05, 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.7624836e05, No Iterations 1 Min/max T:291.29734 746.00059 GAMG: Solving for p_rgh, Initial residual = 0.0041758552, Final residual = 2.1626866e05, No Iterations 4 time step continuity errors : sum local = 9.3353508e07, global = 8.8156252e08, cumulative = 0.81054556 Min/max rho:0.47318442 1.2120324 DILUPBiCG: Solving for epsilon, Initial residual = 0.00044237679, Final residual = 2.0286949e06, 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:
Quote:
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*(1A(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; } The y+ for this model is up to 900. That is probably the reason of it blew up. I also had to put all underrelaxation 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:
Quote:
Code:
solvers { rho { solver PCG preconditioner DIC; tolerance 1e10; relTol 0; } p_rgh { solver GAMG; tolerance 1e10; relTol 0; smoother GaussSeidel; cacheAgglomeration true; nCellsInCoarsestLevel 20; agglomerator faceAreaPair; mergeLevels 2; nSweeps 2; nPreSweeps 1; nPostSweeps 2; } "(Uhkepsilon)" { solver PBiCG; preconditioner DILU; tolerance 1e09; 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; } My Mesh is attached. I just used blockMesh. Quote:
 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 

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


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
SimpleFoam k and epsilon bounded  nedved  OpenFOAM Running, Solving & CFD  15  December 2, 2016 03:40 
SimpleFoam k and epsilon bounded  nedved  OpenFOAM Running, Solving & CFD  1  November 25, 2008 21:21 
Compressible > incompressible.  Jinwon  Main CFD Forum  6  November 23, 2007 22:07 
Problem with LaunderSharma compressible turbulence model wrong formulation  luca  OpenFOAM Running, Solving & CFD  4  July 3, 2006 06:24 
Compressible vs. Incompressible formulations  Fernando Velasco Hurtado  Main CFD Forum  3  January 7, 2000 17:51 