Introduction

The Muscl-Hancock scheme is a 2nd order numerical method. The higher the order of accuracy, the more convergent is the scheme. For instance, here is a plot of the Godunov 1st order upwind scheme (aka CIR scheme, 1st order monotone, what I have been using), the Lax Friedrichs scheme (1st order monotone), Lax Wendroff scheme (2nd order, non-monotone), and Warming-Beam (2nd order, non-monotone).

As you can see, the 1st order monotone schemes are more 'diffusive' than the higher order schemes, i.e. they clip extreme values of the solution. By the 1250 time steps, this effect is highly exaggerated:

On the other hand, the higher order schemes are 'dispersive'. That is, they suffer from slight sign errors in the position of the wave, with the numerical solution either lagging behind or ahead of the true solution.

The above plots show that 2nd order schemes largely perform better than 1st order schemes in regions of 'smooth flow'. However, it is well known that these same higher order accurate schemes can produce spurious oscillations near large gradients. This is due to the higher order schemes no longer being monotone (i.e. having strictly non-negative coefficients in their conservation update formula). Here is a plot showing the same 4 schemes listed above, but for a step function initial condition:

Note the ugly, spurious oscillations near the large gradients of the step function. Unfortunately, all monotone schemes (those that do not contain artificial oscillations) are 1st order at best. Thus, the challenge is to reduce these oscillations in regions of discontinuous flow while maximizing accuracy in regions of smooth flow. One such class of schemes are known as 'Total Variation Diminishing' (TVD). These schemes implicitly add artificial viscosity, resulting in smoother flows near large gradients. Such a scheme is the Muscl-Hancock scheme.

The Muscl Hancock (MH) Scheme

The goal of this method, like others we have looked at, is to solve the system of equations known as the 1D Euler equations,

\frac{\partial \vec{u}}{\partial t} + \frac{\partial f(\vec{u})}{\partial x} = 0

where u is the vector of conserved variables and f is the flux function. Adhering to the IVBP known as the Riemann Problem, we solve this system of equations using the fully discrete, explicit, conservative formula:

u_i^ {n+1} = u_i ^n + \frac{\triangle t}{\triangle x}[f_{i-1/2} - f_{i+1/2}]

for each element of u, where f is the numerical flux. For the 1st order Godunov scheme studied previously, f was taken to simply be the physical flux of 2 adjacent cells, evaluated using the solution to the Local Riemann Problem at the intercell boundary. We have seen various approximation methods to this Godonov scheme, some which estimate the solution of the Riemann problem itself along these intercell boundaries, and others that instead approximate the numerical flux at these boundaries. Next, we will see how the MH scheme produces a 2nd-order accurate solution to the problem. The MH scheme proceeds as follows, being broken down into 3 distinct steps.

1. 'Data reconstruction' - this is the first step of the MH scheme, and is responsible for providing higher order accuracy by fitting a piece-wise linear function to the initial piece-wise constant data. Note that this steps maintains conservation. Such a piece-wise linear function is of the form:

u_i(x) = u_i^n + (x - x_i) \frac{\triangle _i }{\triangle x}, x \in [0, \triangle x]

where \triangle _i is called a 'slope' for u_i ^n , that is for the elements of \vec{u} in cell i on the nth time level. Using a linear function produces 2 extreme values of each cell, one on the left and on the right. These values are crucial in the MH scheme and their role will be presented shortly. A schematic of the situation is as follows:

Image(MusclReconstruction.png, 35%

At this point, we are left with a higher order accurate code, but not one that is free of spurious oscillations near large gradients. To circumvent this, we need to add a TVD measure.

Algorithm Outline

Results

,

,

,