I am solving a concentration equation with the NS equations. Everthing is fine except that, occationally, the concentration values become negative in the area where the value should be very small. I am looking for a method which can numerically enforce the positive requirment. Any suggestion or reference is highly appreciated.

There are advection-diffusion schemes which enforce the positivity condition on the scalar variables in simulations. I do not have a reference with me right now, but you can search in the Journal of Comp. Phys. index using strings "positive" & "scheme" in the title.

If you have difficulty finding any of these papers, I can try to search for one I had copied long time ago.

Please describe completely about governing equations that you should be solved. Do you involve a convection-diffusion problem or diffusion only ? If you have a problem with convection-diffusion, it is more complex than diffusion only. So, you must take care about boundary condition about linearization or others. Thanks

I am a newcomer to CFD and presentlt I am trying to write the Baldwin -lomax mixing length turbulence model .could somebody give me ideas and possibly an algorithm that I can use to implement the model.

Thanks

John

A simple way is just to solve for the log of the function rather than for the function itself.

I am presently doing this for the density and internal energy in the Full Navier-Stokes Equations. I have changes in density or many orders of magnitude and the best way to avoid sharp gradient, small quantities, negative quantities is to use the log of the function. So I am solving for log(rho) where rho is the density and log(e) where e=3/2NkT is the internal energy of the gas. In my case the equations are actually simpler (they have a simpler form) than using the original functions themselves.

For sure you will NEVER get a Negative value, since log is a positive function.

I hope this can help.

Patrick Godon

After comparing several methods, I find this one more interesting. Would you please show me some detail? I am solving a convection-diffusion equation of a concertration field together with NS equations. The equation looks like: (for simplicity, in 1D form)

dc/dt + u*dc/dx = d/dx(k*dc/dx))

What is the equation of log(c) look like? I am not quite sure whether or not the sequence of log operation and derivitives can be changed. Also, I have to model the diffusivity k, which is not a constant. Does this matter?

Thanks

Hi,

Thank you for your interest in my suggestion.

let f = log(c) ; then c = dexp(f)

In Fortran for example alog and dexp are the functions you need and are in basis e (natural logarithmic) NOT in basis 10. So you can keep the basis e, that's just fine.

Then you differentiate:

dc/dx=c * df/dx

and again

d2c/dx2 = dc/dx*df/dx+c*d2f/dx2 =

= c*(df/dx)**2 + c*(d2f/dx2)

Then you put all this back in the original equation (the same hold for time derivatives):

c*df/dt+u*c*df/dx=c*dk/dx*df/dx+k*c*(df/dx)**2+k*c*(d2f/dx2) ,

where d/dx(k*dc/dx)=dk/dx*dc/dx +k*d2c/dx2

Now you can divide by c (if c is not zero) and get an new equation for f:

df/dt+u*df/dx=dk/dx*df/dx+k*(df/dx)**2+k*(d2f/dx2) ,

and you notice that the equation does not even include f, but only its derivatives! So the equation is even simpler. This might be a little different if you have other terms where c does not appear, because then you will get terms like 1/c that you need to write 1/dexp(f). The difficult part is probably just to pay attention and to replace c correctly everywhere else c=dexp(f).

And the results are wonderfull, you can get extremely small value for c!

Cheers,

Patrick Godon

Thank you very much. Your method is very practical. I will extend it to 3D situation and apply to my computation. Thanks again.

Interesting, but your scheme is non conservative anymore, i.e. a finite volume method discretizing the 'c' mass conservation equation in log form will not guarantee mass to be conserved.. So in assuring positive values of 'c' you lose the conservative property of the scheme..

Or am I mislead?

-bernard

I am using a high order accuracy method and still conserve quantities. But of course depending on the method you are using you can have trouble conserving some quantities. THere are given forms used to conserve momentum for example, where you solve for the momenta rather than velocities; or to conserve energy where you include the kinetic energy , etc.. All this is a matter of taste. For sure if you have a problem of negative pressure and negative density you will first be conerned with that and only after that you will want to worry about concerving the mass. You will need to check for the specific problem that you are solving to see where and how much mass is not concerved, if this can be neglected or if you find a better approach, etc... You might want to introduce rectification terms or geometric terms to help conserve quantities.

(1). You have to tell the system that the concentration value must be on the positive side. Since you did not include this property in the formulation, it is possible to get the negative value. (2). So, getting the negative value is normal. (it is like an organization without a cfd department and trying to run cfd codes, most of the time, they are going to get diverged solutions or wrong solutions. I am using this example, because I have been studying the failure modes of angineering organizations.) So, it is important to have someway to correct the course, that is to have human interaction with the codes. (3). Therefore, if you are looking for non-negative solutions, you must include this constraint in the formulation. It shouldn't be difficult to throw out the negative values in the calculation. As a matter of fact, it is a good idea to set a minimum possible level for the concentration. (4). In most applications, the concentration is normalized to one at the inlet, in this way, you have a cut-off value relative to one to follow. The convection and the diffusion of the concentration should not cause the value to become negative. And if for any reason, it becomes negative, you can use the lower limit to correct the behavior. It is also a good idea to find out "why" the concentration is getting negative values. (5). For reacting flows, it is another story. But the same principle applies. I would say that higher-order schemes tends to produce numerical wiggles in a convection-diffusion system. In this case, you will have to refine the mesh until the wiggles are gone. Otherwise, higher-order scheme will simply create non-physical solutions. (6). This is the basic limitation of the finite-difference or finite-volume methods. Unless you re-invent new physics, it is hard to get the real solutions with finite size mesh. Something has to give.