HI All,

I am trying to simluate channel flow using Large eddy simulation, I have been facing a new problem, which is quite complex to comprehend.

Basically for the simulation, I use a Finite Volume Method to discretize the NS eqns, the Viscous terms and the Convection terms are discretized using central difference scheme and the time scheme is Crank-Nicolson scheme, so basically all of the time and space discretization is 2nd order. I also use a Cyclic-Tridiagonal solver to solve the Periodic boundaries.

But when I run the simulation, the problem i am encountering is that the residuals fall for some time, and after some ten of thousands of iterations they become to increase...and the flow field has a lot of oscillations and the code blows up ultimately...!!!!

But this problem is not there, when i use a Hybrid scheme instead of CD for Convective terms...I am not sure whether this problem is due to my careless coding r inherent nature of CD scheme....if so how to combat the same ????

Kindly help me with your expertise....

Thanks in advance,

Glen.

i think it is one of the biggets problems of CDS, there are some bounded version available, their application might help.

Hi Zxaar,

Thanks for the reply, can you eloborate more on the bounded CDS and also the reason, why it keeps the solution stable, can you also refer to me some literature in this regard...

thanks,

Glen.

This is what fluent user manual says:


26.2.7 Bounded Central Differencing Scheme The central differencing scheme described in Section 26.2.6 is an ideal choice for LES in view of its meritoriously low numerical diffusion. However, it often leads to unphysical oscillations in the solution fields. In LES, the situation is exacerbated by usually very low subgrid-scale turbulent diffusivity. The bounded central differencing scheme is essentially based on the normalized variable diagram (NVD) approach [ 189] together with convection boundedness criterion (CBC). The bounded central differencing scheme is a composite NVD-scheme that consists of a pure central differencing, a blended scheme of the central differencing and the second-order upwind scheme, and the first-order upwind scheme. It should be noted that the first-order scheme is used only when the CBC is violated.
this is what they refer in [189],

189 B.P. Leonard. The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection. Comp. Methods Appl. Mech. Eng., 88:17-74, 1991.

I really do not have much idea about BCDS, there might be other guys who will be able tell suggest more, best of luck.

i am new to cfd, and didnt reach les yet, but in general the oscillatory problem of cds for the advection term is well known, and what the fix was uds, but it suffers since it is fijrst order, so QUICK, by Leonard 1979, is the best choice in many simulations I know, and it is third order but I doubt it is really quick wen the mesh fines,however.

Your scheme would appear to be a 'classical' implemention used by many without significant problems. The small amount of viscous dissipation is almost always sufficient to control problems with a reasonable central differencing scheme. Bounded schemes may be required for some equations where the transport quantity is required to be bounded on physical/numerical grounds but this is not the case for the momentum equation.

I list a few areas below which might be related the problem you observe (apart from coding errors which are always quite likely for all of us).

Does your subgrid model allow negative viscosities?

Do you have a variable density?

Do you use grid staggering or pressure smoothing?

Do you check for violation of the convective and the diffusive stability limits at each grid point and each iteration?

Are the boundary conditions for your intermediate velocity component (assuming you have one) correctly handled?

What are your boundary conditions for driving the flow?

Where do the oscillations start grow? Next to a boundary?

Hi Taw,

Well I know that QUICK or infact "Deferred correction" is an alternate approach to avoid Central Difference face to face, but the point is QUICK scheme is a high ordered Upwind scheme, which makes it quite unsuitable for LES, as LES requires a scheme which would produce the least viscosity, for the subgrid scale viscosity to govern the flow domain, I am quite skeptic about trying deffered correction as the same is applicable only when the soultion is converged !!!!

Glen.

It should be possible to carry out a stable computation of a plane channel flow with a CD scheme: many people have done that before. With a Quick scheme (or similiar stuff) you just try to hyde errors in your computation, so that doesn't seem to be the correct approach. Some possible errors: bug(s) in your code, too large time steps (try an explicit scheme), the SGS model doesn't dissipate enough energy (what kind of SGS model do you use?) ....

Hi Andy,

Thanks for the list u have provided me, I have tried to answer the same, let me know abt it...

Does your subgrid model allow negative viscosities?

Ans: I use a Standard smagorinsky model, not a dynamic one, so Negative viscosity is avoided.

Do you have a variable density? Ans: It is a Incompressible flow, so a constant density flow.

Do you use grid staggering or pressure smoothing?

Ans: I do use a staggered grid, for Primitive and Non-primitive variables.

Do you check for violation of the convective and the diffusive stability limits at each grid point and each iteration?

Ans: Actually in order to get rid of the problem, I used a Variable time stepping procedure, based on user-defined (me) courant number, the time-step calcualtion will be based on the Maximum velocity inside the flow field and also the Minimum grid size, usually i set the courant number to be equal to 0.1. I had not checked on the Peclet number yet...

Are the boundary conditions for your intermediate velocity component (assuming you have one) correctly handled?

Ans: I am using a simple algorithm, and the mass residual behaves properly (Until the oscillations).

What are your boundary conditions for driving the flow?

Ans: Along the streawise direction, I use a Constant mass flow rate, which is periodic BC.

Where do the oscillations start grow? Next to a boundary?

Ans: As per my observation, it happens some where inside the domain, if it happens near the boundary, should i check on my BC...kindly advice !!!!

Thanks,

Glen.

Hi Tom,

Thank you for ur comments, actually I am sticking onto CD and I am encountering this problem for the past 2 months, and I have checked and re-checked my code humpty number of times, I am not sure why is this happening, and I am not into QUICK scheme as mentioned by you....and I am trying to simulate a plane channel flow.

I do use an implicit scheme with adavptive time stepping, based on a constant courant number set by me, such that the Time-step is calculated based on the Maximum velocity and the minimum grid size in the domain...

I used a Standard smagorinsky model, with a constant Cs of 0.1, infact i tried to increse this to 0.23, but I have been facing the same problem of spurious oscillations...

thanks,

Glen.

If you are using an explicit time advance of the momentum equations you need to check the diffusion number stability limit. On a boundary/shear layer type clustered grid this is often the stability limit and not the courant number.

> I am using a simple algorithm, and the mass residual behaves properly (Until the oscillations).

I do not understand. SIMPLE is a solution procedure for solving to steady-state and not a time accurate scheme. Do you use a relaxation factor anywhere or just the time derivative term?

> Along the streawise direction, I use a Constant mass flow rate, which is periodic BC.

Again, I do not understand. How do you impose a constant mass flow rate if you use periodic conditions? How do you handle the drop in pressure in the streamwise direction? How do the stresses at the wall arise to balance the drop in pressure in the streamwise direction?

CDS for your problem should be enough, forget about UCDS, QUICK the problem is not in the convective term. Nor in the CFL, 0.1 seems pretty good to me. Are the mesh of good quality (cell not too skewed)? Does increasing the grid size stabilize the problem? Change BC to see what happens? Are you sure the code is reliable?, has been tested or is the first try?

Hi Andy,

On the whole I do not use an explicit scheme for the momentum eqns and I use an Semi-Implicit Crank-Nicolson scheme for the Time discretisation. I think this shd answer the my mis-understanding of using "intermediate velocity component", which in my case there is none !!!

With respect to the constant mass flow rate, the same was posted by me abt 2 weeks before, and Archie mailed me some literature on the same which states that maintaining a constant flow rate across the channel (in which case u shd know the bulk velocity),which in turn will be obtained numerically by relating pressure gradient variations to the flow rate departures from the const intial value....this will help u gain a speedier convegence rather than maintaing a constant pressure-drop across the channel.

In your semi-implicit scheme you have an "intermediate" velocity which you obtain after advancing the momentum equations, you then solve a pressure equation and correct the intermediate velocity to obtain the final velocity? The splitting/factorization/projection or whatever you want to call it of the full set of equations to produce the semi-implicit scheme often leads to the intermediate velocity not being a velocity but a combination of velocity and pressure gradient. If this is the case, the boundary conditions imposed during the first step must reflect this to get a consistent scheme. I strongly doubt violation would be strong enough to make the scheme blow up but it is something to look at.

In my experience, the settling down period is not improved greatly by imposing a constant mass flow but the boundary conditions would require you to calculate and maintain the overall force balance at every time step since you are preventing the flow from accelerating. If you do not maintain a force balance at the boundary the turbulent flow within will have to "adjust" to the force imbalance but without accelerating. Blowing up in the way you describe seems quite likely to me.

If you have doubts about how you maintain the force balance then changing to the normal way of driving the flow via a constant pressure gradient body force and allowing the flow to accelerate should allow you to test if this is where the problem lies.

Hi Glen, I had the same problem as you, and the solution was to rebuild the mesh so that there was almost no stretching.

The commutative errors are increased if you increase mesh stretching, and their effect may be of a negative viscosity (if they stretch in the oposite direction of the velocity vector)

I would suggest you to try first a regular mesh, since everything else looks fine (CFL number, etc...).

How much is your stretching towards the wall, for instance?

Good luck

Hi Luiz,

Thanks for the reply, I will try the same without stretching, but the point as u mentioned that viscosity going negative, it that the case even in the standard Smagorinsky model..I have read some lit saying that cases of dynamic model going negative....due to inherent nature...but i am not sure, in what context will viscosity go negative in this case, can you kindly eloborate...

thanks,

Glen.

Hi Glen, It is not that the added (SGS) viscosity goes negative. It is that the effect of grid stretching (due to commutation errors, which happens when the mesh is not uniform) is equivalent of a negative viscosity, eventought you have not included any explicit viscosity in your equations.

There is a paper from Germano in one of ERCOFTAC conferences, I cannot remember if it is the III or the V...

good luck

PS- I am curious, how much is your grid stretching?

Hi there,

Instead of using some disipative schemes that will have a real nehative inmpact on Turbulent flow simulations (check Mittal And Moin's work), then focus on implementing a kinetic energy conserving scheme (see the work of Morinishi, lund et al. Jour. comput. Physics 1998)

Hope this helps

Fred Felten

Hi Luiz,

Thanks for the help, I am no trying to re-run my code without grid stretching, also the grig stretching that I used was an exponential stretching towards the wall, not much though....say about 1/4 th of the grid points lie within 1/10th of the height of the channel near ot he walls...I did it crudly so that I can save on my run time....will tht have an impact ?

thanks,

Glen.

How much stretching ration (from a adjacent cell to another) this would represent?

In my case it completely solved my problem, but in a channel, there is not much velocity from the wall to the center line.

So I think Fred is right: you may end up having to implement a energy conserving scheme.

But I also think it is worth to see what happens when you remove grid stretching. And you can let a case running while you implement the conservative scheme...

Good luck