Numerical methods

All variables are defined at the same location (collocated). Several different numerical methods are implemented, to allow testing of their accuracy and robustness.

Advection terms $\nabla\cdot\left(\mathbf{b}V_{||}f\right)$

Flux splitting, MinMod limiter

The default method uses a combination of HLL-style flux splitting and MinMod slope limiting. Terms of the form $\nabla\cdot\left(\mathbf{b} f\right)$ are implemented as fluxes through cell boundaries:

\[\nabla\cdot\left(\mathbf{b} V f\right)_i \simeq \frac{1}{J\Delta y} \left[ F_{i+1/2} - F_{i-1/2}\right]\]

where $F$ is the flux. This is calculated by linearly interpolating the velocity to the cell edges

\[V_{i+1/2} = \frac{1}{2}\left(V_{i} + V_{i+1}\right)\]

The field being advected, $f$, is reconstructed from the cell centre values $f_i$ onto cell edges $f^L_i$ and $f^R_i$:

\[f^L_i = f_i - \frac{1}{2}s \qquad f^R_i = f_i + \frac{1}{2}s\]

where the slope $s$ is limited using the MinMod method:

\[\begin{split}s = \left\{\begin{array}{ll} 0 & \textrm{if $\operatorname{sign}(f_{i+1} - f_{i}) \neq \operatorname{sign}(f_{i} - f_{i-1})$} \\ f_{i+1} - f_{i} & \textrm{if $\left|f_{i+1} - f_{i}\right| < \left|f_{i} - f_{i-1}\right|$} \\ f_{i} - f_{i-1} & \textrm{otherwise} \end{array}\right.\end{split}\]

In order to handle waves travelling both left and right, flux splitting handles characteristics moving left differently from characteristics moving right. In general this is problem dependent and computationally expensive, so here we adopt a simple approximation similar to an HLL splitting [2]. We assume that the fastest waves in the system travel with speed $a$ (the sound speed) with respect to the flow, so that there are waves travelling with $V+a$ and $V-a$. If the flow speed is supersonic then these waves are only in one direction, but for subsonic flows there is a flux in both directions. The fluxes between cells are calculated using:

\[\begin{split}F_{i+1/2} = \left\{\begin{array}{ll} f^R_iV_{i+1/2} & \textrm{if $V_{i+1/2} > a$} \\ f^L_{i+1}V_{i+1/2} & \textrm{if $V_{i+1/2} < -a$} \\ f^R_i\frac{1}{2}\left(V_{i+1/2} +a\right) + f^L_{i+1}\frac{1}{2}\left(V_{i+1/2} - a\right) & \textrm{otherwise} \end{array}\right. \label{eq:splitfluxes}\end{split}\]

Hence for subsonic flows the flux becomes $V_{i+1/2}\frac{1}{2}\left(f^R_i + f^L_{i+1}\right) + \frac{a}{2}\left(f^R_i - f^L_{i+1}\right)$, where the second term is a diffusion. When the solution is smooth, $f^R_{i}\simeq f^L_{i+1}$, the numerical method becomes central differencing and the diffusion goes to zero as $\Delta x^2$. Oscillatory solutions introduce dissipation, and the method becomes increasingly upwind as the flow becomes sonic.

Nonlinear fluxes

When advecting quantities which are a nonlinear combination of variables, such as $nV_{||}$, conservation properties can be slightly improved by using the following interpolation [3] [4] [5]:

\[\left(fg\right)^R = \frac{1}{2}\left(f^Rg^C + f^Cg^R\right)\]

where superscript $C$ denotes cell centre, and $R$ right hand side. This method is implemented, using MinMod interpolation for each variable.

Central differencing

Central difference schemes have an advantage over upwind schemes, in that they do not need to take account of wave speeds. The simple central differencing scheme produces large unphysical oscillations, due to the decoupling of odd and even points in collocated schemes, but can (usually) be stabilised by adding dissipation. It is implemented here for comparison with other schemes.

Skew symmetric central differencing

A simple modification to the central differencing scheme improves numerical stability, coupling nearby points [6] [7] The idea is to split the divergence terms into a “skew-symmetric” form

\[\nabla\cdot\left(\mathbf{b} V_{||} f\right) = \frac{1}{2}\left[ \nabla\cdot\left(\mathbf{b} V_{||} f\right) + V_{||}\mathbf{b}\cdot\nabla f + f\nabla\cdot\left(\mathbf{b} V_{||}\right)\right]\]

Each of the terms on the right are then discretised with standard 2nd-order central differences. This method can avoid the need for additional dissipation, or be stabilised with a smaller viscosity than the simple central differencing method.

Artificial viscosity

Artificial viscosity (viscos input) is implemented as a diffusion of momentum in index space, so that the diffusion coefficient varies as $\Delta y^2$.

\[\frac{\partial}{\partial t}\left(nV_{||}\right)_i = \ldots + \nu\left[ \left(V_{i+1} - V_i\right)J_{i+1/2} - \left(V_i - V_{i-1}\right)J_{i-1/2} \right]/J_i\]

where $J$ is the Jacobian, subscript $i$ indicates cell index, and $J_{i+1/2} = \left(J_i + J_{i+1}\right)/2$.