Hi all:

I am developing a code to solve the pure scalar advection quation with sharp discontinuity. The 1st order upwind schemes are known to smear the interface. On the other hand, the higher order schemes (CDS, QUICK) will suffer from the oscillation problem. It seems that the bounded higher-order schemes are the possible solution for this. However, just like the CDS or QUICK, the resulting coefficients discretised by these bounded HR schemes will not have the same signs and hence violate the rule 3 of Patankar(1980). Consequently, I strongly doubt the stability and convergence of these (bounded) HR schemes for the pure advection equation. It seems that the deferred correction is the only means to cure this. Is it possible to implement these HR convection shcemes directly without the use of deferred correction technique ? Can anyone explain or provide a cure for this problem?

Thanks in advance!

If you wish to follow the well trodden route there are many thousands of papers, several books and tens of conferences on the subject of higher order schemes for advection. A book is probably the best source since it should present many of the methods in a consistent form.

Do not confuse stability and convergence with boundedness. A scheme can both converge and be stable with "negative" coefficients (e.g. QUICK).

Consider a scalar field with one part at value 0.0 and the other part at value 1.0 separated by your discontinuity. If a "negative" coefficient is aligned with a scalar value of 1.0 and all the other coefficients lie on a scalar value of 0.0 the "cell" value must be less than zero (unbounded). In order to prevent this you must do something with the information in the scalar field. Either by evaluation, noticing the violation and doing something about it (deferred correction?) or by using the scalar gradients to manipulate the coefficients (non-linear?) which will degrade to first order in the worst case. However, first order need not mean smeared!

As usual, I would suggest a good cure lies in understanding the physics of what you are simulating and adapting the numerical behaviour of your scheme to that physics. An obvious thing to do would be to add into the simulation an additional unknown (the location of the discontinuity) and matching information describing its behaviour. This, of course, side steps the problem you posed by replacing it with another one. However, this one can be solved to give excellent results at least in some circumstances where you can introduce sound physics to describe the behaviour of the discontinuity.

B.P. Leonard (originator of QUICK/QUICKEST) has published a tremendous amount of research on this very topic. Remember that QUICK was devised for essentially steady-state flows and QUICKEST for unsteady advection. The difference between the two being the characteristic traceback, i.e. knowledge of the advection physics (x-t), that is utilized in QUICKEST.

A good starting point is:

Leonard, B.P., M.K. MacVean, and A.P. Lock, "Positivity-Preserving Numerical Schemes for Multidimensional Advection", NASA Tech. Memo. 106055., ICOMP-93-05.

It is interesting to note that most of the work on highly accurate multi-dimensional advection schemes by Leonard, Roache, vanLeer, etc. is based on an explicit approach. This is contrary to the implicit approach you propose since your worried about the matrix coefficients. Obviously implicit schemes are also valid for pure advection problems, but I would recommend a more physically based approach (also recommended by Andy). For example, reconstruct the discontinuity and geometrically compute the volume fluxes at the control volume faces rather than algebraic expressions for the fluxes with flux limiter (which are inherently first-order for contact discontinuities).

(1). Discontinuity means that the variable and the derivatives of the variable are not unique at a point in space. For example, the temperature on one side of the point can be different from the temperature on the other side of the same point. (2). When you have a discontinuity in the flow field, you basically have a problem with two separate solution regions, that is there is an interface. (3). So, there are two problems to consider when studying the flow with discontinuity. One is the motion of the interface location. The other is the behavior of the interface. (4). In the two fluids case, the discontinuity will remain discontinuous all the time. The interface is an absolute boundary. (5). Other types of discontinuity, such as shock waves, mixing layers, are not true interface problems. Shock wave can be formed by compression waves which are not discontinuities. Mixing layers can also diffuse out as they are being convected downstream. (6). Therefore, a complete understanding of the problem is important before one starts modeling the terms in the governing equations. Because, it may requires different sets of equations at the interface location. The interface may not be included as part of interior points. (7). For shock waves and mixing layers, there exists underlying continuous processes . As a result, a macro-model magnified many many times in scale can exist in the approximate solutions. This has been the " capturing models" concept in deriving the continuous solution of flow with discontinuity. (8). I think, the model should reproduce the physics of the problem one is trying to solve. Otherwise, artificial combination of different schemes can only create more artificial results, which has no real physical meaning at all. (9). Many years ago, I was very surprised that a young engineer said that he was able to obtain separated flow solution by running an inviscid code developed at a national lab. He also had velocity vector plots behind a cylinder to prove it. (this was naturally caused by the heavy use of artificial viscosity terms in the inviscid solution, and probably incorrect use of the boundary conditions.)