Higher order convection schemes with unstructured grids
Hi everybody, I was just wondering who may help me in selecting an appropriate convection scheme for a 3D aerodynamics simulation...The grid is composed of about 2 milion of cells (1.8 milion are tetrahedrons, the rest of them are pentahendrons because of the presence of a boundary layer around the solid body). The simulation is performed with the steadystate incompressible solver simpleFoam, with a kepsilon type (standard or RNG) turbulence model coupled with a WF wall treatment. The Reynolds number is about 10^6, and using the upwind firstorder interpolation scheme for the convection terms the numerical convergence is completely reached after 8001000 iterations (at 1000 iterations all the residuals are below 10^05, with the exeption of the pressure residual wich is of about 56*10^05).
The problem is that I would like to use something more accurate than the pure upwind scheme, but irrespective of the bounded first/secon order scheme (vanLeer, limitedLinear, SFCD, etc.) I've tried to use, the simulation blows up after a few iterations... Thank you in advance V. K. 
Any replies, please?

Could you post your dictionaries:
system/controlDict system/fvSchemes system/fvSolution 0/* to check BC's and discretization and solver settings? Regards. 
Quote:
controlDict // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application simpleFoam; startFrom latestTime; startTime 0; stopAt endTime; endTime 1000; deltaT 1; writeControl runTime; writeInterval 200; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression uncompressed; timeFormat general; timePrecision 6; runTimeModifiable yes; // ************************************************** *********************** // fvSchemes ddtSchemes { default steadyState; } gradSchemes { default Gauss linear; grad(p) Gauss linear; grad(U) Gauss linear; } divSchemes { default none; div(phi,U) Gauss upwind; div(phi,k) Gauss upwind; div(phi,epsilon) Gauss upwind; div(phi,R) Gauss upwind; div(R) Gauss linear; div((nuEff*dev(grad(U).T()))) Gauss linear; } laplacianSchemes { default none; laplacian(nuEff,U) Gauss linear corrected; laplacian((1A(U)),p) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected; laplacian(DepsilonEff,epsilon) Gauss linear corrected; laplacian(DREff,R) Gauss linear corrected; } interpolationSchemes { default linear; interpolate(U) linear; } snGradSchemes { default corrected; } fluxRequired { default no; p ; } // ************************************************** *********************** // fvSolution // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { p { solver PCG; preconditioner DIC; tolerance 1e06; relTol 0.01; } U { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0.1; } k { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0.1; } epsilon { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0.1; } R { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0.1; } } SIMPLE { nNonOrthogonalCorrectors 3; } relaxationFactors { p 0.3; U 0.7; k 0.7; epsilon 0.7; R 0.7; } // ************************************************** *********************** // 0/U // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 1 0 0 0 0]; internalField uniform (40 0 0); boundaryField { inlet { type fixedValue; value uniform (40 0 0); } outlet { type zeroGradient; } body { type fixedValue; value uniform (0 0 0); } floor { type fixedValue; value uniform (0 0 0); } side { type zeroGradient; } top { type zeroGradient; } symmetry { type symmetryPlane; } } // ************************************************** *********************** // 0/p // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 2 0 0 0 0]; internalField uniform 0; boundaryField { inlet { type zeroGradient; } outlet { type fixedValue; value uniform 0; } body { type zeroGradient; } floor { type zeroGradient; } symmetry { type symmetryPlane; } side { type zeroGradient; } top { type zeroGradient; } } // ************************************************** *********************** // 0/k // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 2 0 0 0 0]; internalField uniform 0.0096; boundaryField { inlet { type fixedValue; value uniform 0.0096; } outlet { type zeroGradient; } body { type kqRWallFunction; value uniform 0.0096; } floor { type kqRWallFunction; value uniform 0.0096; } side { type zeroGradient; } top { type zeroGradient; } symmetry { type symmetryPlane; } } // ************************************************** *********************** // 0/epsilon // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 3 0 0 0 0]; internalField uniform 0.00155; boundaryField { inlet { type fixedValue; value uniform 0.00155; } outlet { type zeroGradient; } body { type epsilonWallFunction; value uniform 0.00155; } floor { type epsilonWallFunction; value uniform 0.00155; } side { type zeroGradient; } top { type zeroGradient; } symmetry { type symmetryPlane; } } // ************************************************** *********************** // 0/nut // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 1 0 0 0 0]; internalField uniform 0; boundaryField { inlet { type calculated; value uniform 0; } outlet { type calculated; value uniform 0; } body { type nutWallFunction; value uniform 0; } floor { type nutWallFunction; value uniform 0; } side { type calculated; value uniform 0; } top { type calculated; value uniform 0; } symmetry { type symmetryPlane; } } // ************************************************** *********************** // 
Skewness might be an issue. Can you post a checkMesh output?

Quote:
Create time Create polyMesh for time = 0 Time = 0 Mesh stats points: 435948 faces: 4149850 internal faces: 4014985 cells: 1991792 boundary patches: 7 point zones: 0 face zones: 0 cell zones: 0 Overall number of cells of each type: hexahedra: 0 prisms: 197667 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 1794125 polyhedra: 0 Checking topology... Boundary definition 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 symmetry 37174 19931 ok (nonclosed singly connected) floor 26808 13741 ok (nonclosed singly connected) inlet 596 334 ok (nonclosed singly connected) outlet 544 306 ok (nonclosed singly connected) side 2208 1194 ok (nonclosed singly connected) top 1646 909 ok (nonclosed singly connected) body 65889 33196 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (1.26 0 0.1945) (6.27 1.44 1.1945) Mesh (nonempty, nonwedge) directions (1 1 1) Mesh (nonempty) directions (1 1 1) Boundary openness (3.81773e20 1.86644e19 8.03877e20) OK. Max cell openness = 2.59975e16 OK. Max aspect ratio = 6.2611 OK. Minumum face area = 1.57523e06. Maximum face area = 0.0146675. Face area magnitudes OK. Min volume = 9.77998e10. Max volume = 0.000524117. Total volume = 10.7878. Cell volumes OK. Mesh nonorthogonality Max: 58.2126 average: 16.0731 Nonorthogonality check OK. Face pyramids OK. Max skewness = 1.02516 OK. Mesh OK. End 
I was about to ask for the same as Sandeep, thx him to remember this point. Skewness seems to be OK. Maybe you can reach nonOrthogonalCorrecions to an extreme, i.e. 20 and try again, if problem is solved then tune this parameter properly (in order to reduce calculation time).
Regards. 
Hi, there,
what are the initial conditions you're using for starting up the simulation with the second order convection schemes? In case you haven't already tried: sometimes it helps taking the (converged) results of the first order simulations as the initial field for further calculations. Greetings. 
Quote:

Quote:
I've tried also to start up the simulation from the last iteration of the upwind run (with all the residuals significantly below 10^05), but the results all almost the same for all the higher order scheme mentioned above... 
Hi all, I have some news:
 I've tried to use the blended differencing scheme, but if the simulation has to start from the timestep "0" it blows up again after at most 20 iterations (even with a blending factor of 0.01, wich if I'm not wrong means an upwind "dominance" of 99%...)  I've also tried to start from the iteration number 600 (for the upwind run with the keps standard model it means substantial numerical convergence for many of the flow parameters): in this case I was able to push the blending factor till 0.15, obtaining a good convergence pattern for the next 100 iterations...it is better than nothing, but to my knowledge the solution will be still very diffusive (am I wrong about it?) However, any suggestions? PS  when I use the blended scheme, the fvSchemes dictionary is the same as posted above, replacing Gauss upwind with Gauss blended <blending factor>) 
Hi,
sorry, I can't really help you with the blended interpolation scheme, didn't even know it existed in OpenFOAM before. But here's another idea: you're replacing all the divSchemes with higher order schemes when switching to second order, right? You might want to try out the following: only change one convection scheme at a time to second order (e.g. div(phi,U) Gauss upwind; to div(phi,U) Gauss SFCDV;) and see if the simulation blows up. Sometimes the convective terms of the turbulent quantities are very sensitive to the respective discretization scheme. Greetings. 
Quote:
Thank you once again V. K. 
Hey, there,
I'm looking forward to your testing results. In case that the convection scheme of velocity flux works in second order while the turbulent quantities still blow up (when using higher order discretization), we have knowledge what might be causing the instabilities when using second order discretization (namely the turbulent quantities). As a further step it might be helpful to experiment with the underrelaxation parameters for k and epsilon then. According to your other question: it usually isn't easy to tell about accuracy when speaking of the order of the discretization schemes since it is influenced by many parameters (such as type of flow, mesh, boundary conditions, solver and so on). What we know for sure is that the numerical error we make with e.g. Gauss linear is of second order, while Gauss upwind produces an error of first oder etc. We only know the order of the error, not the magnitude. In my experience it usually improves the numerical results when using second order discretization methods only on the velocity flux as the velocity field per se is of a very special interest in many CFD applications, but like I said: the effect regarding accuracy is depending on many things. First order discretization on the convection terms of the turbulent quantities is not that bad in my opinion, since the turbulence models are based upon empirical correlations anyway  but nonetheless: it of course is better to try getting the order of the discretization error as high as possible, if the goal is to have accurate results. A general recommendation is very hard to make. Greetings. 
Quote:
V. K. 
5 Attachment(s)
I'm back...
Ok, once again Felix was right: the turbulent quantities seem to be much more sensitive to the convection scheme type than the mean velocity field...assuming that, I've made several runs changing the scheme for the velocity and that's what I have obtained: limitedLinearV 1, GammaV 1, SFCDV : the solution exhibits good convergence patterns, but a bit slower than the run with the upwind scheme. The GammaV 1 and SFCDV schemes seem to behave better than the limitedLinearV 1 one. By the other hand, the QUICKV and UMISTV schemes seem to be quite unstable, and this is not excactly what I've expected, because to my knowledge the QUICK implementation in OpenFOAM should be bounded (and, judging by some papers involving the UMIST scheme, it also should be a QUICKbased bounded scheme). In addition to these considerations, I will add the charts of the residuals on a semilogarithmic scale. Waiting fro your comments... V. 
Hi, Vesselin,
your residual histories are looking great (except for the QUICKcase), I'm glad I could help you out a bit. Now it might be interesting to see, if you can use convection schemes of higher order on the turbulent quantities when highly underrelaxing k and epsilon. Another option would be refining the mesh in (according to turbulence) critical regions like the boundary layer and possible wakes. Regarding your observation when using the QUICK scheme: I am not entirely sure, but I think QUICK  in it's original form  is an unbounded scheme. Maybe someone else can tell us more about that. I don't know if OpenFOAM uses a bounded version of this scheme  I'm gonna have a look at the implementation tomorrow or the day after and tell you more about it. Greetings, Felix. 
Quote:
From some cases I've run limiterLinear performs better than QUICK as vkrastev showed. Best. 
4 Attachment(s)
Quote:
my intention for the next steps is to select the two most promising schemes (which seem to be the Gamma and the SFCD ones) and then try to reduce the underrelaxation factors for k and epsilon. Of course, some additional work on the mesh could always be done, but for now (and in my opinion) this is not a priority, because I think that for this kind of simulation (steadystate RANS with Wall Functions at the wall) the quality of the mesh in question can be considered satisfactory (or, at least, good enough to make some numerical tests as I've been doing so far): about this I can tell you that near the solid body there is a boundary layer with regularshaped (and heightcontrolled) prisms, with a firstcellheight of 2 mm (which, at the Reynolds number I'm using for the runs, corresponds to an yplus varying between 40 and 130, so in good agreement with a loglayer based clasical WF); from the boundary layer, a fully unstructured tetrahedron mesh is built, with a minimum cell dimension of 5 mm (side of the tetra control volume) and with a gradual growing apart from the body and into the wake; as you can see from my previous posts, the quality check from the checkMesh utility seems to be quite good, an so seems the check made with Tgrid (the max Fluent skewness reported is of about 0.79). But maybe I'm spending too much words to describe something that could be much better explain with an image, so I'll post some details of the mesh below... About the QUICK scheme, I also know (both from Jasak's paper and from Leonard original paper of 1979) that in its original formulation It is not bounded (and, theoretically talking, third order accurate), but In the OpenFOAM usr's guide it is reported as a bounded first/second order scheme, so here it is my surprise in founding it much more unstable than other bounded schemes...However, to be honest, I have not investigated the "realworld" implementation of QUICK in OpenFOAM, so as you sade we should check it before doing too much hypothesis... However, thanks for your comments and "see" you soon Regards Vesselin 
Quote:
this is exactly wath I'm talking about: as prof. Jasak said, the QUICK original scheme should be unbounded, but in the OpenFOAM usr's guide it is presented as a bounded first/second order bounded scheme (and, as you have noticed, it seems to be more unstable than other bounded schemes)... Regards V. 
All times are GMT 4. The time now is 17:51. 