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.

Hi, Vesselin, are you modelling the Ahmed body? If this were the case or in other aerodynamics cases boundary layer modelation and pressure recovery in the back of the bodies is crucial, surely you know about that. I comment this relative to the value of y+ and the mesh pics, it's necessary to fulfill the requirements of being in the log layer, with a value of y+ ~ 30 if it is possible. This implies to change the size of near wall elements along the wall.
But, in the meanwhile, good luck with the div schemes.

Regards.

Yes, I'm modeling the Ahmed body, starting with the simplest case (symmetrical domain, steady state run, body without stilts, 25° slant angle case). You're right about the mesh dimension near the solid body: I know that it should be slightly refined, in fact when I talk about "satisfactory quality" I'm referring mainly to convergence, numerical stability and physical "acceptability" of the solution. When I'll finish with this preliminary testing phase, for sure I will follow your advice and refine the grid near the walls.
Thanks for your comment

Regards

V.

Yes, sometimes it requires to refine the mesh, but in other parts mesh must be coarsened in order not to have finer nor coarser meshes. In other words y+ mustn't be in laminar sublayer zone nor outer region.


Regards.

Hi, Vesselin,


your intention for the next steps sounds good to me. I hope you'll be successful that way!
And please excuse me if I wasn't perfectly clear regarding the mesh. I was aware of the checkmesh-results and didn't think your mesh is not suitable or something like that, the results are indeed very good and there won't be any problems because of skewness or non-orthogonality, I think. And judging by the pictures you attached everything looks very fine to me.
What I actually meant was that one of the (mesh-wise) critical regions for a blunt body like that might be the wake behind the body. Especially the turbulent kinetic energy can have large gradients in such regions which might require more cells in said region. Remember: higher-order schemes are much more sensitive to cell size (amongst other things), which may cause the instabilities, if the mesh is to coarse in regions of high gradients. But this is just guesswork and it may be helpful to have a look at the contour plots of k and epsilon from the simulations with lower order discretization to get an idea, where it may be neccessary to refine the mesh. It could be an issue, but it doesn't have to be - like you said: your mesh is pretty good already.

Regarding the QUICK scheme: I just had a quick peep into the source files but I didn't came to a proper conclusion, I need to have a deeper look, but I'm too tired right now, sorry.
But there's a clue inside the description:

It could be the downwind contribution which might cause the instabilities, but I can't really tell right now. Maybe this gives you an idea?


Greetings,
Felix.

Hi Felix,
don't worry about the mesh: I posted a more detailed description and the images only for giving you and the other forum readers a better idea of how the mesh itself was built...after this "numerical trial" phase I'll start to compare the results also from a quantitative and physical point of view (there is a very large numerical and experimental literature about the Ahmed body) and I'm quite sure that further mesh refinement and assestment will be needed....
About the QUICK scheme: to my actual knowledge, I can't tell if this kind of downwind contribution could be an instability source...this matter need further investigation, so indeed I also have to go depper in the code (and to enlarge my theorical background)...Hope we'll reach soon some agreement about it
About my tests with the underrelaxation factors: I think that in a few days I will be able to post some results, so wait for me...

Regards

V.

Hi all,

I'm sorry I'm a little late...here there are some results of my numerical tests on the convection scheme for the turbulent quantities. As I said before, assuming that k and epsilon are much more sensitive to the convection scheme type than the velocity field, I've tried to "impose" a higher order scheme for the turbulent quantities by "playing" with the under-relaxation factors. My considerations about this matter are the following:

-Having in mind that my tests are focused (for now) only on the realizable k-epsilon model running in steady-state configuration (simpleFoam solver) and mainly on the Gamma differencing scheme, trying to "force" the convergence of the turbulent quantities (once a higher order scheme for convection, such as the Gamma one, is imposed) by lowering a lot the under-relaxation factors seems to be not a good idea...First of all, as you can see from the attachments, the convergence patterns are not satisfactory because the residuals of k and epsilon tend to remain quite unstable and the residuals of the other quantities remain too high (compared with my previous runs). Secondly, the results are also definetly not physical (the tke distribution in the wake of the body as well as the velocity fields are far from what one should expect by considering previous experimental and numerical works on the Ahmed model).

-Of course, my testing campaign is incomplete (unfortunately I have not so much time neither computational facilities), because I've not covered all the range of combinations between the under-relaxation factors and differencing schemes and I've kept the mesh and the eddy viscosity turbulence model as "fixed parameters", but maybe one can deduce a kind of "trend" in the turbulent quantities behaviour for similar cases and meshes.

About the images: three of them are referred to relaxation factor values identical for k and epsilon (0.1, 0.3, 0.5), the fourth is a combination of two different values (0.1 for epsilon, 0.5 for k). In all of the trials, relaxation factors for U and p are kept as their "standard" values for the SIMPLE algorithm (0.3 for p, 0.7 for U)

V.

Hi all, I'm back again...
It seems that OpenFOAM (at least the 1.6 version) is really, really, sensitive to the convection schemes when someone tries to use a tetrahedron-based unstructured 3D mesh...If you have a look to all the posts above, at least I came to the conclusion that a good compromise between accuracy and convergence of the (steady) solution should be a NVD/TVD scheme for the velocities (such as GammaV 1 or limitedLinearV 1) and the good old first order upwind scheme for the turbulent quantities (k and epsilon in my case)...well, I was wrong, because when I use a practically identical mesh, but with very small differences in the cell distribution (for example, the first cell height of the prism-shaped boundary layer is 1.5 mm instead of 2 mm), the same settings in the fvSchemes tend to diverge slowly, until the solution blows up...So, after, all, I'm still asking for some advice if I'm doing something wrong with the overall settings...In particular I want to know:

1) If with such settings the very instable behaviour of the code (simpleFoam solver) is, by your experience, quite "normal"...
2) Can I change something apart from the divSchemes which could improve the stability of the solution?
3) Now I'm starting to do some tests with the linerUpwind scheme for the velocity: can someone give me some more information abput this scheme? I've searched in the forum and had a look in the user's Guide, but still I cannot understand properly how this scheme is formulated...Am I wrong in saying that it's a second order upwind scheme?

Down below there are my fvSchemes and fvSolution dictionaries. For the BC's you can have a look on the previous posts

Thank you very much

V.

fvSchemes:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
default steadyState;
}

gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}

divSchemes
{
default none;
div(phi,U) Gauss GammaV 1;
div(phi,k) Gauss upwind;
div(phi,epsilon) 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;
}

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;
}


}

SIMPLE
{
nNonOrthogonalCorrectors 3;
}

relaxationFactors
{
p 0.3;
U 0.7;
k 0.7;
epsilon 0.7;
}


// ************************************************** *********************** //

Vesselin,

I am really glad that you and felix had this discussion here. I am also going through a similar phase with my research in modeling tip vortices.

One thing which I noticed with the divSchemes ( as you guys have already pointed out ) is that , the turbulence parameters are very sensitive to the scheme. I am using openfoam 2.0.1 and I get converged results ( less than 10^-5) when I used Gauss linearUpwind. Any other higher order turbulence model for all of the divergence terms resulted in a blowing up of the solution.

Changing the scheme just for the velocity term to a higher order ( QUICK, QUICKV) atleast keeps the run going but the residuals never converge to decent values( hanging around 10^-3).

From all this, I tend to feel that the not just the div schemes but may be any other schemes(especially the interpolation schemes) might have an effect on the stability and the convergence of the solution. May be people who have had experience can explain better.

Kalyan

I am also facing the problem of nonconvergence .Can you give me some advice.The case files are here:http://www.cfd-online.com/Forums/ope...t351846.Thanks in advance.

As far as i remember, QUICK requires hexa mesh, didnt it? This experience is from Fluent, but it would explain behaviour of your QUICK residuals :)