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 (u=u(x,t)), and f is the flux function. Adhering to the IVBP known as the Riemann Problem (i.e. constant left and right data states separated by an initial discontinuity), 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 (ul and ur). These values are crucial in the MH scheme and their role will be presented shortly. Their values are given by

u_i ^l = u_i(0) = u_i ^n - \frac{1}{2} \triangle_ i

and

u_i ^r = u_i(\triangle x) = u_i ^n + \frac{1}{2} \triangle_ i

A schematic of the situation is as follows:

Using these end values, one has the new Riemann problem:

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

u(x,0) = \{ ^{u_i^r,~~~~ x <0}_{u_{i+1}^l, ~~x>0}

where 0 is the local origin, i.e. any given intercell boundary. We will discuss meaningful ways to set the extreme values, ul and ur, shortly. For now we will just call the data-reconstruction step that which sets the extreme values using the above formula. Given that the expressions for ul and ur involve the slope \triangle _i, we define

\triangle_i = \frac{1}{2}(1+\omega)\triangle u_{i-\frac{1}{2}} + \frac{1}{2}(1-\omega)\triangle u_{i+\frac{1}{2}}

where w is a free parameter in the real interval [-1,1].

\triangle u_{i-\frac{1}{2}} \equiv u_i^n - u_{i-1}^n

\triangle u_{i+\frac{1}{2}} \equiv u_{i+1}^n - u_{i}^n

2. Boundary evolution - this step 'evolves' the left and right values of u (ul and ur). It is these values of ul and ur that are used in the Riemann problem outlined in step 1. The formulas used to evolve each boundary extrapolated value are:

\bar{u}_i^l = u_i^l + \frac{\triangle t}{2 \triangle x} [f(u_i ^l) - f(u_i^r)]

\bar{u}_i^r = u_i^r + \frac{\triangle t}{2 \triangle x} [f(u_i ^l) - f(u_i^r)]

A schematic of the situation is depicted here for reference:

3. Solve the new Riemann problem

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

u(x,0) = \{ ^{\bar{u}_i^r,~~~~ x <0}_{\bar{u}_{i+1}^l, ~~x>0}

using the evolved extrapolated boundary values. The solution of this problem is then used in the Godunov flux, which is used in updating cell i. Recall the Godunov flux is simply,

f_{i+1/2} = f(u_{i+1/2}(x/t=0))

i.e. the numerical flux at the intercell boundary equals the physical flux evaluated at the similarity solution of the LRP (local Riemann Problem) at that boundary.

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. There are some different ways of doing this, i.e. limited slopes, flux-limiters, and slope limiters. Using one of these methods, you would replace the slope function in steps 1 and 2 with a TVD version of the slope. This leads to a scheme that is both higher order accurate as well as minimizes spurious oscillations.

To summarize, the MH scheme provides a means to solve for the numerical flux of the conservative update formula for cell i . This flux is just the regular Godunov flux, but now evaluated using the solution of the LRP consisting of left and right states that fit a piece-wise linear function. This linear function is chosen to both a) increase the order of accuracy, and b) satisfy TVD conditions. In the end we are left with a numerical solution that is both less dispersive in regions of smooth flow and less oscillatory near sharp gradients.

Results

As you can see here, the MH scheme + Superbee slope limiter + Roe solver produce highly non-diffusive results. However, there still are some small scale spurious oscillations near the boundary of the square wave.

,

,

,