
[Sponsors] 
Nonphysical pressure field when changing from first to secondorder time advance 

LinkBack  Thread Tools  Search this Thread  Display Modes 
March 8, 2016, 20:20 
Nonphysical pressure field when changing from first to secondorder time advance

#1 
Member
Join Date: Aug 2015
Posts: 37
Rep Power: 11 
I'm running LES of a turbulent jet using ~1.6 million cells and a Courant Number of 0.3. My solver is customized, but all of the nonstandard features are turned off and it should be behaving identically to reactingFoam with all reactions turned off. The solver works fine with secondorder spatial schemes and firstorder (Euler) time advance, but when I try to switch to secondorder (backward) time advance, the pressure field develops nonphysical "beads" of high and low pressure. If I leave the simulation for long enough, many other nonphysical things happen to the pressure field, and the simulation often diverges or stalls (delta t goes to zero).
The evolution of the p field is illustrated here: http://imgur.com/a/fcF2I [note: the colour scale varies between images]. The test case: a turbulent jet with laminar coflow The boundary conditions:
Code:
Mesh stats points: 1617892 faces: 4824873 internal faces: 4797927 cells: 1603800 faces per cell: 6 boundary patches: 5 point zones: 0 face zones: 0 cell zones: 0 Overall number of cells of each type: hexahedra: 1603800 prisms: 0 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 0 polyhedra: 0 Checking topology... Boundary definition OK. Cell to face addressing OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces... Patch Faces Points Surface topology JET 405 424 ok (nonclosed singly connected) PILOT 864 900 ok (nonclosed singly connected) COFLOW 1404 1440 ok (nonclosed singly connected) OUTERRADIUS 21600 21636 ok (nonclosed singly connected) OUTLET 2673 2692 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (0.07172601826 0.07172601826 0) (0.07172601826 0.07172601826 0.36) Mesh (nonempty, nonwedge) directions (1 1 1) Mesh (nonempty) directions (1 1 1) Boundary openness (2.5360889850782e16 2.502408869494276e16 3.295433281148142e16) OK. Max cell openness = 4.165398281849697e16 OK. Max aspect ratio = 22.5683629471026 OK. Minimum face area = 5.055728864375342e08. Maximum face area = 3.708057695313814e05. Face area magnitudes OK. Min volume = 1.872805105800212e11. Max volume = 4.104974063197774e08. Total volume = 0.003402726340535457. Cell volumes OK. Mesh nonorthogonality Max: 22.23980401213051 average: 2.729797707905458 Nonorthogonality check OK. Face pyramids OK. Max skewness = 0.5165905340833714 OK. Coupled point location match (average 0) OK. Mesh OK. Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 2.1.0   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default backward; } gradSchemes { default Gauss linear; } divSchemes { div(phi,rho) Gauss limitedLinear phi 1; default Gauss limitedLinear 1; div(phi,U) Gauss filteredLinear2V 0.2 0; div(phi,Yi_h) Gauss multivariateSelection { CO2 limitedLinear01 1; H2O limitedLinear01 1; NO limitedLinear01 1; N2 limitedLinear01 1; h filteredLinear2 0.2 0; }; div((muEff*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } fluxRequired { default no; p; } // ************************************************************************* // Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 2.1.0   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { "(prho)" { solver PCG; preconditioner DIC; tolerance 1e08; relTol 0; // relTol 0.01; } "(prho)Final" { $p; relTol 0; } "(UhZvarZYikepsilon)" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e06; relTol 0; // relTol 0.01; } "(UhZvarZYikepsilon)Final" { $U; relTol 0; } } PIMPLE { momentumPredictor yes; nOuterCorrectors 50; nCorrectors 3; nNonOrthogonalCorrectors 1; turbOnFinalIterOnly false; residualControl { U { tolerance 1e4; relTol 0; } p { tolerance 1e4; relTol 0; } } } relaxationFactors { fields { p 0.3; pFinal 1; } equations { "(Ukepsilon)" 0.7; "(Ukepsilon)Final" 1; } } // ************************************************************************* // My question: Why is the p field becoming nonphysical, and what can I do to fix it while maintaining secondorder accuracy? I have tried tightening the solver tolerances and adding adding more inner and nonorthogonal correctors, but they do not seem to help. I was thinking that I should add more PIMPLE iterations, but 56 already seems like a lot to me. Should I be reducing the Courant number? I would be happy to do this for a transition period, but if I have to keep Co very low for the remainder of my simulation then I think that this could become prohibitively expensive. 

March 12, 2016, 09:33 

#2 
Senior Member
Artur
Join Date: May 2013
Location: Southampton, UK
Posts: 372
Rep Power: 20 
Hi,
Have you tried the CrankNicolson time scheme? It also gives you an option to blend in a bit of firstorder time which could help with your issue. Also, have you tried running the solver with a single outer loop, i.e. in PISO mode, or is that not an option for your case? I found that the latter sometimes helps with limiting problems with convergence and reduces run times, although may require you to reduce the time step a bit to keep the residuals under control. Furthermore, I take it you are using adjustable time step? If so, maybe try switching it off and fixing delta t to something which you already know will keep the Co<0.5 or whatever criterion you are following. I've found the adjustable time step with backward scheme to yield very similar looking issues with one of the old OF versions when using a multiphase solver. All the best, A 

March 12, 2016, 14:56 

#3 
Member
Join Date: Aug 2015
Posts: 37
Rep Power: 11 
Hi Artur,
Thanks for the reply. These are all new ideas to me, so I will try them out! Why would you expect that switching from PIMPLE to PISO would be more stable? Is the idea to switch from (PIMPLE with bigger time step) to (PISO with smaller time step), with the expectation that the two changes will kind of cancel out, or maybe even lead to improvements in both simulation time and stability? 

March 13, 2016, 06:21 

#4 
Senior Member
Artur
Join Date: May 2013
Location: Southampton, UK
Posts: 372
Rep Power: 20 
Not sure exactly what your solver is doing, but for some of my cases with cavitation I've found that big time steps with pimple may cause trouble in other transport equations, which then pollutes your solution and eventually causes divergence. You said you switched all the extra pieces off but just maybe there is still something there which is causing issues, hence my idea. Could be entirely wrong though


March 16, 2016, 00:22 

#5 
Member
Join Date: Aug 2015
Posts: 37
Rep Power: 11 
Hi Artur,
I tried some of your suggestions. Here are my results so far:


March 16, 2016, 04:45 

#6 
Senior Member
Artur
Join Date: May 2013
Location: Southampton, UK
Posts: 372
Rep Power: 20 
I agree with you in terms of rough cost equivalence of the two algorithms. It's probably a good idea to try PIMPLE with reduced time step just to get a full picture of how the problem behaves in different configurations. If it still persists there are always the standard questions of BCs and whether the mesh is sufficiently dense. Perhaps it'd also be worth to run the case with a builtin pimple/piso solver using the same grid and setup just to verify there isn't anything in the custom one you're using that's causing the issue?


March 16, 2016, 14:11 

#7 
Member
Join Date: Aug 2015
Posts: 37
Rep Power: 11 
Hi Artur,
My newest results:
I can keep reducing the time step, but I feel like Co = 0.19 is all ready pretty low, no? 

March 16, 2016, 15:02 

#8 
Senior Member
Artur
Join Date: May 2013
Location: Southampton, UK
Posts: 372
Rep Power: 20 
Yes it is very low from a flow point of view, although from reaction point of view I have no idea to be honest. I suppose this brings the checklist to the BCs and mesh... After that I'm really out of ideas since your schemes and solvers look OK to me (I normally wouldn't use smoothSolver myself but that should give these strange results I reckon). Did you base your BCs on an existing case or a tutorial?


March 16, 2016, 16:13 

#9 
Member
Join Date: Aug 2015
Posts: 37
Rep Power: 11 
My BCs are based on the "pilotedDiffusionFlame" tutorial case from the "flameletfoam" solver: https://openfoamwiki.net/index.php/E...n/flameletFoam. The geometry in the flameletfoam test case is quite similar to mine, but they run RANS (and thus use a 2D domain with wedge boundaries) whereas I run LES (and thus use a 3D domain with those boundaries removed). Their tutorial files are on github at https://github.com/flameletFoam/flam...ionFlame/ras/0. For reference, my jet Reynolds number (rho * v_bulk * d_jet / mu) is 22,400, which I believe is much larger than the value for the flameletfoam test case.
I made the mesh myself. I don't consider myself a meshing expert, but I assume that it isn't horribly broken given the fact that the checkMesh utility has no complaints about it. I can post the blockMeshDict if that might help. 

March 17, 2016, 03:35 

#10 
Senior Member
Artur
Join Date: May 2013
Location: Southampton, UK
Posts: 372
Rep Power: 20 
The mesh is probably fine in terms of quality, I'm no expert in nonwallbounded flows so not sure what kind of resolution you need to resolve turbulent stresses to a level sufficient for the sgs model, so it's probably worth checking that with the literature unless you're confident in the values you've chosen.
On thing which caught my attention is the proximity of the outlet and the fact that you (the tutorial) use inletOutlet on U and fixedValue on p while from your images one can tell there's quite a lot of turbulence reaching the outlet. I realise this is quite expensive but maybe moving the outlet further and switching the U BC to pressureInletOutletVelocity would help? I found it to work very well with fV on p. Optionally, I think you might also look into using a convective outlet BC on U which afaik is supposed to be design for the exact kind of situation like what you're having now (turbulence reaching edge of the domain), there's some code and more discussion here: http://www.cfdonline.com/Forums/ope...adyflows.html I'm sorry but if none of these help then I'm officially out of ideas A 

April 15, 2016, 20:29 

#11  
Member
Ali Shamooni
Join Date: Oct 2010
Posts: 44
Rep Power: 15 
Quote:
Hi, What is ur diameter, velocity and level of turbulence intensity of the jet? 

March 31, 2020, 04:25 

#12 
New Member
Hongsheng
Join Date: Dec 2017
Posts: 2
Rep Power: 0 
Hi knuckles,
did you solve the nonphysical p fields problem? i recently meet a similar pressure issue when i use the 2nd time scheme in reactingfoam to calculate jet flame, or can you please give me some advice? 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
static vs. total pressure  auf dem feld  FLUENT  17  February 26, 2016 13:04 
Superlinear speedup in OpenFOAM 13  msrinath80  OpenFOAM Running, Solving & CFD  18  March 3, 2015 05:36 
How to write k and epsilon before the abnormal end  xiuying  OpenFOAM Running, Solving & CFD  8  August 27, 2013 15:33 
dynamic Mesh is faster than MRF????  sharonyue  OpenFOAM Running, Solving & CFD  14  August 26, 2013 07:47 
Could anybody help me see this error and give help  liugx212  OpenFOAM Running, Solving & CFD  3  January 4, 2006 18:07 