I'm using explicit Mac Cormack scheme to solve a discontinuity problem. When I do not add artificial viscosity, the solution has irregular oscillation. When I add strong artifical viscosity, the oscillation become rather regular - smooth amplitude and regularly over grid node (1-D) alternately. Could you guru tell me the possible mechanism in it please? Is there any simple method to eliminate or lessen osciliation of Mac Cormack scheme? What if I can sacrifice accuracy?

(1). Well, that's the nature of the method. (2). There's an old paper which covers 1-D shock problem with different schemes. You may want to take a look at it. ( G.A.Sod, J. Comp. Physics,27(1) (1978)) (3). I took a couple of methods from the paper and wrote a program in BASIC using PC back in early 80's. (4). The method development in this area (transient compressible flows) has been very active through out 80's. So, you should be able to find many papers on the subject from Journal of Computational Physics (1978 to 1988). (5).I think the MacCormack method is easy to understand and easy to program. As for the oscillation in the solution, I didn't try to fix it. So, I don't have the answer for you. That is how the method " capture" the shock. In other words, the shock with oscillations is just a simulation.

Thanks for your reply. Maybe you can give me more clue. Actually I'm not working on shock, but diffusion. The discontinuity is in concentration of species. Pressure is almost constant and velocity is almost zero - the problem is that the near-zero velocity is taken into account. Then, what I got is oscillation in velocity and pressure, but the discontinuity is in species! Pressure is calculated from equation of state seperately, which may be the real catch too.

(1). It is important to specify your problem clearly in the first place. (2). Well, as you said the velocity is very low, then the Mach number is probably very low. That's just my guess. Most methods suitable for supersonic compressible flows have the same difficulty at low Mach numbers. This is because the flow at low Mach number behaves like incompressible flow. (3). At the incompressible flow limit, the pressure field can be derived from the velocity field. And the method of solution has to be changed. This is why there are pressure-based methods developed for incompressible flows.

Hi there,

it seems that there are about three points that you might want to consider when you change and implement your method of solution. The three problems that I see are the following:

- The discontinuity of course

- The numerical differentiation

- The small velocity and constant pressure

The discontinuity is always a source of oscillations (sometimes refered to as the Gibbs phenomenon) and the use of an artificial viscosity should reduce their amplitude especially far from the discontinuity. This should also smooth out the discontinuity over more grid points (increasing the number of grid points such that the initial discontinuity is spread over more points should also help). Since this does not seem to be the case in your flow, the two next points should be considered.

Some methods of numerical differentiation are such that the even grid points are coupled together, and the odd too, but the even and odd are not coupled together in any manner. THis is equivalent to solving separately for the odd and the even points and therefore getting two different solutions for the even and odd points (for example the odd have a given boundary condition at 1, but the even don't have any at node number 2). This creates the two points oscialltions. To remdediate to this problem one has to change the scheme, or alter it slightly, for example by introducing a kind of weight function (like in finite elements) which will help to couple the even and odd points to a certain degree.

If the velocities are small and the pressure almost constant, one has to solve for the small departure of the velocity and pressure only and not for the whole quantities. Namely, it is better to write:

V=V0+U

P=P0+p , etc...

where V0=0 and P0=constant are the steady state solution (or the known unperturbed analytical solution), and then put P and V in the equations keeping in mind that P0 and V0 are known and constant and solving only for u and p. Yor new variables are p and u. Then P and V are obtained by adding P0 and p, and V0 and u, etc.. Since P is almost constant, the small changes of P (namely p) are lost in the errors when calculating P. So it is better to get rid of the errors on P (P0).

I hope this helps. PG.