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 steady-state incompressible solver simpleFoam, with a k-epsilon type (standard or RNG) turbulence model coupled with a WF wall treatment. The Reynolds number is about 10^6, and using the upwind first-order interpolation scheme for the convection terms the numerical convergence is completely reached after 800-1000 iterations (at 1000 iterations all the residuals are below 10^-05, with the exeption of the pressure residual wich is of about 5-6*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.

First of all, thanks for the reply. Secondly, those are the required dicitonaries:

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((1|A(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 1e-06;
relTol 0.01;
}

U
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
}

k
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
}

epsilon
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
}

R
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
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?

Here it is the checkMesh response:


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 (non-closed singly connected)
floor 26808 13741 ok (non-closed singly connected)
inlet 596 334 ok (non-closed singly connected)
outlet 544 306 ok (non-closed singly connected)
side 2208 1194 ok (non-closed singly connected)
top 1646 909 ok (non-closed singly connected)
body 65889 33196 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (-1.26 0 0.1945) (6.27 1.44 1.1945)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (3.81773e-20 1.86644e-19 8.03877e-20) OK.
Max cell openness = 2.59975e-16 OK.
Max aspect ratio = 6.2611 OK.
Minumum face area = 1.57523e-06. Maximum face area = 0.0146675. Face area magnitudes OK.
Min volume = 9.77998e-10. Max volume = 0.000524117. Total volume = 10.7878. Cell volumes OK.
Mesh non-orthogonality Max: 58.2126 average: 16.0731
Non-orthogonality 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.

After about 8 corrections the solver stops to iterate because the pressure residual falls under 10^-06. However, setting the number of correctors to 8 the simulation blows up again after 4-5 iterations (I've tried only with the limitedLinear scheme for the convection terms, just to have an idea if it could be a useful solution).

The conditions are the same posted above, replacing the upwind scheme with the higher-order scheme (limitedLinear/limitedLinearV, vanLeer/vanLeerV, SFCD/SFCDV, etc.)

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 k-eps 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.

Hi Felix, thanks for the reply. You're right: in my prior trials I've changed all the divSchemes at the same time. So, I will test your idea and see what happens. Apart from this, I want to pose you another question: what will be, in your opinion, the impact on the overall solution accuracy of a "mixed" convection scheme treatment (I mean, for instance, a higher order scheme for the velocity flux and a first order scheme for the turbulence quantities fluxes)? Let's say that it will work: will it be significantly more accurate than the "all-first-order" option? Of course, I know that a rigorous answer would require some numerical tests, but I'd just like to know your personal-knowledge-based opinion (and maybe the opinions of the other forum users) about it

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.

Thank you very much for all the answers. Ok, now I'll try first with the velocity convection scheme and then with the other options (underrelaxation parameters and so on). Wait for additional news...

V. K.

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 QUICK-based bounded scheme).

In addition to these considerations, I will add the charts of the residuals on a semi-logarithmic scale.
Waiting fro your comments...

V.

Hi, Vesselin,


your residual histories are looking great (except for the QUICK-case), 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.

Following Jasak (http://powerlab.fsb.hr/ped/kturbo/op...GammaPaper.pdf), this scheme is slightly unbounded and diffusive. Nevertheless some implementations can make it bounded. In FOAM docs (http://www.openfoam.com/docs/user/fvSchemes.php) it is presented as a First/second order, bounded squeme.

From some cases I've run limiterLinear performs better than QUICK as vkrastev showed.

Best.

Hi Felix,
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 under-relaxation 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 (steady-state 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 regular-shaped (and height-controlled) prisms, with a first-cell-height 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 log-layer 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 "real-world" 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

Hi Santiago,

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.