
[Sponsors] 
SimpleFoam high timestep continuity error on a simple geometry 

LinkBack  Thread Tools  Search this Thread  Display Modes 
October 24, 2014, 12:10 
SimpleFoam high timestep continuity error on a simple geometry

#1 
New Member
Chris Stiapis
Join Date: Mar 2014
Posts: 19
Rep Power: 12 
Dear al,
My name is Christos and I'm currently trying to simulate a distribution system for my PhD. The distributions system has an inlet and 4 outlets. And I'm using simpleFoam. The inlet is at the bottom and the outlets on the right. I meshed the Geometry with the use of Salome and I got this mesh (on a branch) which is not bad on my opinion The problem is that when I;m using simpleFoam with kEpsilon the timecontinuity explodes. When I'm using realizablekE timestep error increases but but slowly. When I'm using kOmega again the timestep explodes. I'm attaching the 0 folder and the system folder so if anyone wants may have a look. Mesh orthogonality 50 and y+37. The velocity boundary conditions boundaryField { inlet_extruded { type fixedValue; value uniform (0 0 0.44211); } outlet_extruded { type inletOutlet; inletValue uniform (0 0 0); value uniform (0 0 0); } walls_extruded { type fixedValue; value uniform (0 0 0); } frontAndBack { type empty; } } // And the pressure boundaryField { inlet_extruded { type zeroGradient; } outlet_extruded { type fixedValue; value uniform 44500; } walls_extruded { type zeroGradient; } frontAndBack { type empty; } } // The fvSolution solvers { p { solver GAMG; tolerance 1e08; relTol 0.00005; smoother GaussSeidel; cacheAgglomeration true; nCellsInCoarsestLevel 20; agglomerator faceAreaPair; mergeLevels 1; } U { solver smoothSolver; smoother GaussSeidel; nSweeps 2; tolerance 1e07; relTol 0.00005; } k { solver smoothSolver; smoother GaussSeidel; nSweeps 2; tolerance 1e07; relTol 0.01; } epsilon { solver smoothSolver; smoother GaussSeidel; nSweeps 2; tolerance 1e07; relTol 0.1; } } SIMPLE { nNonOrthogonalCorrectors 1; pRefCell 0; pRefValue 0; } relaxationFactors { fields { p 0.3; } equations { U 0.5; k 0.5; epsilon 0.5; omega 0.5; } } I really don't uderstand what is happening ???????? I get GAMG: Solving for p, Initial residual = 0.276687, Final residual = 8.57431e06, No Iterations 14 smoothSolver: Solving for Uy, Initial residual = 0.00637261, Final residual = 5.2594e08, No Iterations 8 smoothSolver: Solving for Uz, Initial residual = 0.0119267, Final residual = 6.9414e08, No Iterations 8 GAMG: Solving for p, Initial residual = 0.0189897, Final residual = 4.52346e07, No Iterations 15 time step continuity errors : sum local = 0.255507, global = 0.00466411, cumulative = 6.20251 smoothSolver: Solving for epsilon, Initial residual = 0.0277328, Final residual = 0.00100213, No Iterations 2 bounding epsilon, min: 1.56859e+12 max: 3.84711e+13 average: 3.64094e+10 GAMG: Solving for p, Initial residual = 0.0189897, Final residual = 4.52346e07, No Iterations 15 time step continuity errors : sum local = 0.255507, global = 0.00466411, cumulative = 6.20251 smoothSolver: Solving for k, Initial residual = 3.51264e12, Final residual = 3.51264e12, No Iterations 0 ExecutionTime = 2561.74 s ClockTime = 3104 s 

October 24, 2014, 15:29 

#2 
Senior Member

Hi,
and what are ICs and BCs for k, epsilon, and omega? Can you also post fvSchemes? Also, looking at your geometry I can't grasp why do you need triangular mesh. 

October 27, 2014, 07:03 

#3 
New Member
Chris Stiapis
Join Date: Mar 2014
Posts: 19
Rep Power: 12 
Alexey thanks for replying,
I'm using triangular mesh cause this geometry and mesh is generated by an automated routine in Salome, where i use the salome functionality to batch generate geometry and meshes. I tried and experimented with many schemes and now I'm using that one that is considered very stable although diffusive. I took a look in the bibliography and it states that in case of multiple outlets only a velocity boundary at the inlet needs to be specified and pressure at the outlet. The velocity U boundaryField { inlet_extruded { type fixedValue; value uniform (0 0 0.44); } outlet_extruded { type zeroGradient; } walls_extruded { type fixedValue; value uniform (0 0 0); } frontAndBack { type empty; } And the pressure P dimensions [ 0 2 2 0 0 0 0 ]; internalField uniform 44900; boundaryField { inlet_extruded { type zeroGradient; } outlet_extruded { type fixedValue; value uniform 44900; } walls_extruded { type zeroGradient; } frontAndBack { type empty; } } Also the fVSchemes ddtSchemes { default steadyState; } gradSchemes { default cellLimited leastSquares 1.0; } divSchemes { default none; div(phi,U) bounded Gauss linearUpwind cellLimited Gauss linear 1; div(phi,k) bounded Gauss upwind; div(phi,epsilon) bounded Gauss upwind; div(phi,R) bounded Gauss upwind; div(R) Gauss linear; div(phi,nuTilda) bounded Gauss upwind; div((nuEff*dev(T(grad(U))))) Gauss linear; div(phi,omega) bounded Gauss upwind; } laplacianSchemes { default Gauss linear limited 0.5; laplacian(1,p) Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default limited 0.333; } fluxRequired { default no; p ; } 

October 27, 2014, 08:10 

#4 
Senior Member

According to your log file, epsilon equation is a source of the problems:
Code:
time step continuity errors : sum local = 0.255507, global = 0.00466411, cumulative = 6.20251 smoothSolver: Solving for epsilon, Initial residual = 0.0277328, Final residual = 0.00100213, No Iterations 2 bounding epsilon, min: 1.56859e+12 max: 3.84711e+13 average: 3.64094e+10 Also if you know pressure are the outlet, why not use pressureDirectedInletOutletVelocity BC for velocity? And are you sure all your outlets will have the same pressure? 

October 27, 2014, 08:50 

#5 
New Member
Chris Stiapis
Join Date: Mar 2014
Posts: 19
Rep Power: 12 
Dear Alexey, this distribution system is in the bottom of a vesssel(not pressured) so the pressure on the outlet is the hydrostatic pressure and as long as the outlets share the same height i suppose that the pressure will be equal.
I decided to neglect the turbulence Effects till i get my case converged in laminar. I did as u suggested and I'm using the pressureDirectedInletOutletVelocity condition for velocity and the the total pressure. I read somewhere in this forum that this coupling is stable. I'm getting again a explosion in timestep continuity ????? pressureDirectedInletOutletVelocity smoothSolver: Solving for Uy, Initial residual = 0.243622, Final residual = 1.00739e08, No Iterations 13 smoothSolver: Solving for Uz, Initial residual = 0.271367, Final residual = 1.00003e08, No Iterations 13 GAMG: Solving for p, Initial residual = 0.410951, Final residual = 4.08045e06, No Iterations 413 time step continuity errors : sum local = 2.12258, global = 0.445773, cumulative = 0.296711. The velocity boundary now is dimensions [ 0 1 1 0 0 0 0 ]; internalField uniform (0 0 0); boundaryField { inlet_extruded { type pressureDirectedInletOutletVelocity; inletDirection uniform (0 0 1); value uniform (0 0 0.44); } outlet_extruded { type zeroGradient; } walls_extruded { type fixedValue; value uniform (0 0 0); } frontAndBack { type empty; } } // and the pressure dimensions [ 0 2 2 0 0 0 0 ]; internalField uniform 44900; boundaryField { inlet_extruded { type totalPressure; U U; phi phi; rho none; psi none; gamma 1.4; p0 uniform 244900; } outlet_extruded { type fixedValue; value uniform 44900; } walls_extruded { type zeroGradient; } frontAndBack { type empty; } } // ????? 

October 27, 2014, 09:04 

#6 
Senior Member

Well,
If your case is laminar, why even bother about turbulence model? If your case is turbulent, it won't converge without turbulence model. We can continue to play in a game called "guess real physical conditions in your system" (to suggest right boundary conditions), or you can try to describe in more details what you're trying to model. Right now your BCs don't make sense. 

October 27, 2014, 10:48 

#7 
New Member
Chris Stiapis
Join Date: Mar 2014
Posts: 19
Rep Power: 12 
Thank you Alexey,
My case is turbulent (RE>50000). What i was saying earlier was to minimize the mesh sensitivity by minimizing the velocity so it will be laminar in order to test the boundary conditions without worrying about the mesh sensitivity. Now, my case. There is a distribution system in the bottom of a vessel. Therefore the pressure on the outlet is known from the hydrostatic pressure. The pressure on the inlet is known, cause there is a pump. The velocity at the inlet is known because of the flowrate of the pump. Can I ask you, what do u think is wrong on the previously described boundary conditions??? My epsilon BCS internalField uniform 0.00000001; boundaryField { inlet_extruded { type fixedValue; value $internalField; } outlet_extruded { type zeroGradient; } walls_extruded { type epsilonWallFunction; value $internalField; } } frontAndBack { type empty; } And My K boundary conditions internalField uniform 0.0000001; boundaryField { inlet_extruded { type fixedValue; value uniform 0.00750571558131; } outlet_extruded { type zeroGradient; } walls_extruded { type kqRWallFunction; value $internalField; } frontAndBack { type empty; } } although at some point converges due to the very loose tolerances 

October 28, 2014, 03:49 

#8 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
Chris, you use the wrong syntax for the underrelaxation factors. Accordingly, OpenFOAM doesn't recognize them. You need to put
Code:
relaxationFactors { fields { "(p)" 0.3; } equations { "(U)" 0.7; "(k)" 0.8; "(epsilon)" 0.8; } } If you post code, hit the "Go advanced" button below the text box and use the "#" sign, to wrap "code" tags. This is much easier to read. Also some tabs will help to understand the code. Try this and post again, if it worked or not.
__________________
The skeleton ran out of shampoo in the shower. 

October 28, 2014, 04:33 

#9 
New Member
Chris Stiapis
Join Date: Mar 2014
Posts: 19
Rep Power: 12 
Thank you guys for replying, really appreciated.
I changed it and it doesn't seem to make a difference. I think that in OpenFoam 2.3.0 the previous systax is allowed as well as the one suggested. No I'm doing the simulation and it seems to converge very slowly, #Time = 159 smoothSolver: Solving for Uy, Initial residual = 0.000998984, Final residual = 5.55782e10, No Iterations 14 smoothSolver: Solving for Uz, Initial residual = 0.000555796, Final residual = 7.87531e10, No Iterations 13 GAMG: Solving for p, Initial residual = 0.434494, Final residual = 9.3188e09, No Iterations 1000 time step continuity errors : sum local = 9.16547e08, global = 2.3667e10, cumulative = 2.67703e09 smoothSolver: Solving for epsilon, Initial residual = 0.000423026, Final residual = 4.34165e10, No Iterations 9 smoothSolver: Solving for k, Initial residual = 8.5116e10, Final residual = 8.5116e10, No Iterations 0 ExecutionTime = 44.03 s ClockTime = 47 s Time = 160 smoothSolver: Solving for Uy, Initial residual = 0.00092316, Final residual = 4.05633e10, No Iterations 14 smoothSolver: Solving for Uz, Initial residual = 0.000517306, Final residual = 7.15796e10, No Iterations 13 GAMG: Solving for p, Initial residual = 0.434494, Final residual = 9.3188e09, No Iterations 1000 time step continuity errors : sum local = 9.16547e08, global = 2.3667e10, cumulative = 2.67703e09 smoothSolver: Solving for epsilon, Initial residual = 0.000423026, Final residual = 4.34165e10, No Iterations 9 GAMG: Solving for p, Initial residual = 0.434494, Final residual = 9.3188e09, No Iterations 1000 time step continuity errors : sum local = 9.16547e08, global = 2.3667e10, cumulative = 2.67703e09 smoothSolver: Solving for k, Initial residual = 8.5116e10, Final residual = 8.5116e10, No Iterations 0 ExecutionTime = 45.38 s ClockTime = 48 s # I beievee that this convergence is due to the very stict tolerances that I have imposed in fvSolution #p { solver GAMG; tolerance 6e10; relTol 0.00000001; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; cacheAgglomeration on; agglomerator faceAreaPair; nCellsInCoarsestLevel 10; mergeLevels 1; } "(UkepsilonomegaRnuTilda)" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e09; relTol 0.0000001; }# I really don't know what is going wrong. I neverthought that this case will get so complicated . The turbulence mdoel I'm using is realizableKE. 

October 28, 2014, 04:45 

#10 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
"relTol" normally doesn't need to be lower than 0.1 or 0.01 in simple. The nonlinear error of the outer iterations is much larger anyway. In any case, 0.00000001 is insane
I think for SIMPLE, you should try to have just a few iterations in all equations, such as maximum 5. Also, you don't actually need the nonorthogonal correction, as you mesh isn't that bad. Hrvoje Jasak stated somewhere in this forum, that you should only use them, if your mesh orthogonal quality is larger than 70°. In your fvSchemes:
__________________
The skeleton ran out of shampoo in the shower. 

October 28, 2014, 05:09 

#11 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
Edit: Ok, you are right, my fault.
__________________
The skeleton ran out of shampoo in the shower. 

October 28, 2014, 10:01 

#12 
Senior Member

Hi,
Just as PoC I've made similar case and it converges quite quickly. My BCs (I've decided not to mess with turbulence, so it's just p and U): inlet: U > fixedValue p > zeroGradient outlet: U > zeroGradient (as I don't expect backflow, otherwise something more fancy can be used) p > fixedValue (as you've said) walls: U > nonslip p > zeroGradient You can find a case attached to the message. Concerning your turbulence BCs  where did you get these constants for k and epsilon? RNG? Experimental values? Estimations from turbulence intensity and hydraulic radius? 

October 28, 2014, 10:05 

#13 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
I can highly recommend to use Alexey's settings. These are the most basic / common settings for "p" and "u" and thus (let's hope) the most stable ones. I use them for all my pipe simulations, also in Fluent. I would not try to use any others, just because some person in the forum told you so, except you really know what you are doing.
Regarding my upper post, I just found this: http://www.openfoam.org/mantisbt/view.php?id=1410#c3246 A statement of henry: Also when running the case I improved the schemes, in particular it is a VERY bad idea to limit all gradients: the pressure gradient should NEVER be limited.
__________________
The skeleton ran out of shampoo in the shower. 

October 28, 2014, 11:56 

#14 
New Member
Chris Stiapis
Join Date: Mar 2014
Posts: 19
Rep Power: 12 
For the inlet
turbulent intensity I=0.16Re−1/8. The turbulent length scale can be estimated as k=32(UI)2, ϵ=cμ3/4k3/2l−1.andl=0.07L, with L a characteristic length. For internal flows this may take the value of the inlet duct (or pipe) width (or diameter) or the hydraulic diameter. This is what I'm using in order to calculate the l=0.07L And the k file Code:
internalField uniform 0.075; boundaryField { inlet_extruded { type fixedValue; value uniform 0.075; intensity 0.160000; Code:
internalField uniform 0.02; boundaryField { inlet_extruded { type turbulentIntensityKineticEnergyInlet; intensity 0.16; // 16% turbulent intensity value $internalField; } boundaryField { inlet_extruded { type turbulentMixingLengthDissipationRateInlet; mixingLength 0.02; // 0.02m  half pipe diameter value $internalField; } The weird thing is that it converges although the timestep continuity is high??? Code:
Time = 3072 smoothSolver: Solving for Uy, Initial residual = 0.000147551, Final residual = 8.79075e07, No Iterations 7 smoothSolver: Solving for Uz, Initial residual = 0.00014121, Final residual = 7.46354e07, No Iterations 7 GAMG: Solving for p, Initial residual = 0.00998099, Final residual = 9.77039e05, No Iterations 7 time step continuity errors : sum local = 0.000324409, global = 4.5041e06, cumulative = 0.0117465 smoothSolver: Solving for epsilon, Initial residual = 3.33489e05, Final residual = 1.33563e07, No Iterations 4 smoothSolver: Solving for k, Initial residual = 8.77387e10, Final residual = 8.77387e10, No Iterations 0 ExecutionTime = 584.71 s ClockTime = 586 s SIMPLE solution converged in 3072 iterations 

October 28, 2014, 12:13 

#15 
New Member
Chris Stiapis
Join Date: Mar 2014
Posts: 19
Rep Power: 12 
The problem that I'm suspecting now,as alexey was saying, maybe there is an issue with the IC.
I'm a starter now in CFD so I'm a little bit lost 

October 28, 2014, 14:30 

#16 
Senior Member

Well, if ICs affect your steady state solution, I guess there's something wrong with the solution
If you take a look at src/finiteVolume/cfdTools/incompressible/continuityErrs.H: Code:
{ volScalarField contErr(fvc::div(phi)); scalar sumLocalContErr = runTime.deltaTValue()* mag(contErr)().weightedAverage(mesh.V()).value(); scalar globalContErr = runTime.deltaTValue()* contErr.weightedAverage(mesh.V()).value(); cumulativeContErr += globalContErr; Info<< "time step continuity errors : sum local = " << sumLocalContErr << ", global = " << globalContErr << ", cumulative = " << cumulativeContErr << endl; } 

October 29, 2014, 02:39 

#17 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
1) Did you set the schemes as I suggested in post # 10?
2) You can increase the convergence criteria? Now it says it converged just because all initial residuals are below 1e3. Set them to 1e12 and see how low they can get.
__________________
The skeleton ran out of shampoo in the shower. 

October 29, 2014, 07:57 

#18 
New Member
Chris Stiapis
Join Date: Mar 2014
Posts: 19
Rep Power: 12 
I did as u suggested, that's why took so long to answer,many iterations,
Although, the continuity has been greatly reduced the cumulative error remains high Code:
Time = 13146 smoothSolver: Solving for Uz, Initial residual = 3.10449e06, Final residual = 2.44747e09, No Iterations 8 GAMG: Solving for p, Initial residual = 8.14567e05, Final residual = 7.23971e08, No Iterations 10 time step continuity errors : sum local = 1.39972e07, global = 4.68047e09, cumulative = 0.000204088 smoothSolver: Solving for Uy, Initial residual = 2.03736e06, Final residual = 1.46767e09, No Iterations 8 smoothSolver: Solving for Uz, Initial residual = 3.13649e06, Final residual = 2.42425e09, No Iterations 8 smoothSolver: Solving for epsilon, Initial residual = 5.74286e07, Final residual = 2.65494e10, No Iterations 6 smoothSolver: Solving for Uz, Initial residual = 3.09894e06, Final residual = 2.41098e09, No Iterations 8 smoothSolver: Solving for Uy, Initial residual = 2.05149e06, Final residual = 1.42138e09, No Iterations 8 smoothSolver: Solving for Uz, Initial residual = 3.06732e06, Final residual = 2.39964e09, No Iterations 8 smoothSolver: Solving for epsilon, Initial residual = 5.66788e07, Final residual = 2.62298e10, No Iterations 6 smoothSolver: Solving for Uy, Initial residual = 2.02899e06, Final residual = 1.45363e09, No Iterations 8 smoothSolver: Solving for Uy, Initial residual = 2.091e06, Final residual = 1.47467e09, No Iterations 8 smoothSolver: Solving for k, Initial residual = 5.72958e07, Final residual = 4.05752e10, No Iterations 6 smoothSolver: Solving for Uz, Initial residual = 3.08764e06, Final residual = 2.3985e09, No Iterations 8 

October 29, 2014, 08:18 

#19 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
1) Chris, to get the inlet profiles, I usually run a case with just a short pipe / channel with periodic boundary conditions with the shape of the inlet of my second case. Then, I use "mapFields" to map the outlet of my first case (U, k, omega, epsilon, nu) to the inlet of my second case. Now, I have physically meaningful inlet conditions. See this post for a sketch:
http://www.cfdonline.com/Forums/openfoampreprocessing/141196mapfieldsquestion.html#post508830 2) For the other problem: I have never seen, that the cumulative error remains high, while all others decrease. What pressure and velocity inlet and outlet b.c. do you use right now? 3) If it converges, you can also start to change convective schemes to "linearUpwind". If you get numerical problems you can set a limiter, but just for the upwind gradient scheme: Code:
gradSchemes { default Gauss linear; grad(linUpwind) faceLimited Gauss linear 1.0; } divSchemes { div((nuEff*dev(T(grad(U))))) Gauss linear; div(phi,U) Gauss linearUpwind grad(linUpwind); div(phi,k) Gauss linearUpwind grad(linUpwind); div(phi,epsilon) Gauss linearUpwind grad(linUpwind); }
__________________
The skeleton ran out of shampoo in the shower. 

October 29, 2014, 10:33 

#20 
New Member
Chris Stiapis
Join Date: Mar 2014
Posts: 19
Rep Power: 12 
My velocity and Pressure
U Code:
FoamFile { version 2.0; format ascii; class volVectorField; object U; } dimensions [ 0 1 1 0 0 0 0 ]; internalField uniform (0 0 0); boundaryField { inlet_extruded { type fixedValue; value uniform (0 0 0.44211); } outlet_extruded { type zeroGradient; } walls_extruded { type fixedValue; value uniform (0 0 0); } frontAndBack { type empty; } } // ************************************** Code:
dimensions [0 2 2 0 0 0 0]; internalField uniform 49000; boundaryField { inlet_extruded { type zeroGradient; } outlet_extruded { type fixedValue; value uniform 49000; } walls_extruded { type zeroGradient; } frontAndBack { type empty; } } This procedure with the mapfields looks promising although i have never used this function before. 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Rapidly decreasing deltaT for interDyMFoam  chrisb2244  OpenFOAM Running, Solving & CFD  3  July 1, 2014 16:40 
AMI interDyMFoam for mixer nu problem  danny123  OpenFOAM Programming & Development  8  September 6, 2013 02:34 
plot over time  fferroni  OpenFOAM PostProcessing  7  June 8, 2012 07:56 
Full pipe 3D using icoFoam  cyberbrain  OpenFOAM  4  March 16, 2011 09:20 
Convergence moving mesh  lr103476  OpenFOAM Running, Solving & CFD  30  November 19, 2007 14:09 