I am working on simulating the internal flow in the pipe of fuel injection system. I used Explicit Mc Cormack scheme + FCT in order to eliminate the oscillation near the shock. After appearance the first shock FCT works very well but then after reflection at the end of the pipe I got trouble with a negative values which is physically impossible and very strong oscillation.

Did someone have experience with this problem ?

thanks

arief

Howdy,

Which quantities are going negative?

If its conserved quantites like density, for example, then you definitely have a problem with your implementation, and there is probably a bug. So start looking. Make sure your timestep is OK, and so on..

If its derived quantities such as pressure or temperature then this is not as big a deal. Just use MAX functions to cap the quanties from below. There are some numerical issues with deriving pressure from the equation of state and plugging that back in to your simulation. Especially when the shocks become particularily strong.

Dan W.

It's not un-common to have negative density in the compressible flow calculations. So, you are not alone. The negative density issue is not just related to the shock wave calculation ( strictly speaking, shock is almost a discontinuity ). Just a few weeks ago, when I was using a commercial code ( written by a professor ) to solve a 3-D turbine blade design problem, I had negative density in the low speed inlet region, after around 20,000 iterations. In my case, the non-uniform mesh used in the inlet coupled with the accelerated local time-step was problbly the cause of the problem. In your case, a first step to do is to create a mesh which is consistent with your shock systems ( easier when solving steady-state solution ). The other parameter is the time -step. Here, try to use uniform time-step over the computational domain ( this may take longer to converge ). Between the two factors, the mesh has the greatest impact on the solution behavior. For the transient problem, ideally, you would like to use a mesh system which also change in time so that it is consistent with your solutions. Based on my experience with the McCormack explicit method, it is possible to get sharp shock system when the mesh is consistently defined.

Hi Dan, Thanks for your advice,

You're right, I have problem with a negative pressure in my case. Could you tell me more about MAX function and maybe some papers where I can find about it.

thanks in advance

arief

Arief,

There are not really any papers regarding this issue when using FCT codes. It's more that it's just a thing that people do with FCT codes. It's pretty simple actually, just cap the pressure from below:

do i = 1, n pressure(i) = whatever from your equation of state pressure(i) = MAX(pressure(i),pmin) enddo

where pmin is a minimum pressure you would expect in your simulation. Some people use 0.1*pmin instead as well. Since pressure is not a conserved monotone quantity, doing this is not really a big deal.

This problem arises from solving the momentum and energy equations separately and coupling them with the equation of state. The momentum and energy get out of phase, say by one grid point, around the area of the shock and you end up with negative pressures when you derive them from the EOS. It does'nt really have anything to do with the mesh or time-step as suggested in the posting by John Chien. It also happens on really fine meshes or with whatever time steps you choose. It's really a problem with the way in which you solve the conservation equations.

You can do the same for temperature and other derived quantites as well.

Dan.

John,

It's not uncommon to have negative densites in compressible flow computations if you are using codes that were not really meant for that. If you have a proper implementation of a hyperbolic compressible flow codes, conserved postitive definite quantities should remain that way. It's just the way it is. I suggest that you maybe try using a different code, something based on FCT, PPM, ENO, or some other TVD technique. Maybe look in the books by Hirsch for some guidance.

You seem to be using an implicit solver to get steady state solutions. Arief on the other hand is using a fully explicit technique for a transient flow calculation. It is easy to argue that implicit solvers are not very good for these sorts of things. Especially when you want to resolve quite fast phenomenon, where the advantage of using a large timestep with an implicit solver is not that great.

However, if you are using an implicit solver without any kind of monotonicity preserving algorithms, your advice is perfectly correct. I guess I'm just trying to point out that implicit solvers are not the greatest idea for transient compressible flows.

Dan.

If that what you mean, I have done it, and it didn't work very well. Let me ask you a little bit more detail. for mac cormack scheme I used (abs(u)+a)*DT/DX <=1 as a stability criteria, and then I coupled it with FCT. Refer to Fletcher, FCT approach has more restrictive stability, u*DT/DX<0.5. Do you think I should use the second stability also when I couple Mac Cormack with FCT ? I have tried the second stability for both scheme but I got a very smeared profiles of shock.

thanks,

arief

Arief,

It's odd that you did not have much success with the limiting from below. I've done tons of computations using it, and it worked great for me. Are you saying that doing it doesn't do anything at all for you? What happens? Does the computation just crash or what?

The stability criteria really depends on your particular implementation. This is something you should work out for yourself.

The classical FCT originally developed by Boris and Book is stable for all courant numbers less than 1, but only guarantees *monotone* solutions for courant numbers less than 0.5, that is the proper distinction. Try running a shock tube problem sometime with like 0.75, it will run, but you'll have oscillations all over the place. I typically used courant numbers of 0.4 with most of my work. There are FCT implementations which are monotone for all courant numbers below 1 as well, but I never really tried them out. I guess that's why I'm saying you should work out for yourself where your scheme will produce monotone solutions.

When you say you see smearing using lower courant numbers, how much smearing are we talking about. A shock that is resolved over 3-4 grid points is probably fine. If you are talking like 8-10, this would be a lot. FCT codes should be able to resolve a shock over a few grid points without any significant diffusive effects.

Again, if you are still having problems, you should carefully check your implementation. I know people that have done good work with MacCormack-FCT implementations which worked quite well for them.

You might also try the LCPFCT implementation from NRL, there are full 2D and direction split versions available on the NRL web site.

BTW, Who is "Fletcher"? I've only read the original FCT stuff done by all the people from the Naval Research Lab where it came from, i.e., Boris and his colleauges.

Maybe we should take this chat strictly to email, so as not to choke up the CFD site with more discussion.

Dan W.

The method used in the code I mentioned is explicit method. Since it is a black box, I was only trying to share my experience with others on Internet.

There has been a lot of success in improving the robustness and accuracy of shock calculations during the 1980s.

I recall several papers between 1984 and 1986 on this subject. If my memory is correct check out: J. Comp Physics, 1984 for a pair of articles by Woodward and Colella. They talk about many different methods while presenting their work on PPM. I believe this is still a landmark paper.

I believe Gary Sod also review several methods (c. 1979).

If you are considering changing your approach (as suggested by Dan Williams, and I concur), you might want to review the work of Bram van Leer (1982-1985). His techniques (Flux Vector Splitting) deal with monotonicity quite well, are not limited by Mach number, and are simple and elegant. I believe his work was extended in the late 1980s with TVD techniques at NASA Langley.

I believe the PPM code is still considered the best, but is not generally available. It requires considereable effort to recreate this from the articles in the literature.

Arief

I have used FCT + Mac Cormack explicit and faced the same problem of getting negative values. I over come it by:- (a) Changing order of differencing in the predictor and corrector at alternate time cycles. (b) While applying FCT, Choosing the order of differencing for calculation of fluxes near boundaries so that the resulting flux is always physically consistent, e. g., heat flux always directed to a heat sink.

Ref: AIAA Journal, 34, #9, 1835-1842