Algebraic Flux Correction II

Abstract

Flux limiting for hyperbolic systems requires a careful generalization of the design principles and algorithms introduced in the context of scalar conservation laws. In this chapter, we develop FCT-like algebraic flux correction schemes for the Euler equations of gas dynamics. In particular, we discuss the construction of artificial viscosity operators, the choice of variables to be limited, and the transformation of antidiffusive fluxes. An a posteriori control mechanism is implemented to make the limiter failsafe. The numerical treatment of initial and boundary conditions is discussed in some detail. The initialization is performed using an FCT-constrained L ² projection. The characteristic boundary conditions are imposed in a weak sense, and an approximate Riemann solver is used to evaluate the fluxes on the boundary. We also present an unconditionally stable semi-implicit time-stepping scheme and an iterative solver for the fully discrete problem. The results of a numerical study indicate that the nonlinearity and non-differentiability of the flux limiter do not inhibit steady state convergence even in the case of strongly varying Mach numbers. Moreover, the convergence rates improve as the pseudo-time step is increased.

Appendix

In this appendix, we derive the artificial diffusion operator for the piecewise-linear Galerkin approximation to the one-dimensional Euler equations

$$\frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}=0.$$

(118)

In the 1D case, we have

$$U=\left [ \begin{array}{c}\rho\\ \rho v\\ \rho E\end{array} \right ],\qquad F=\left [\begin{array}{c}\rho v\\\rho v^2+p\\ \rho Hv\end{array} \right ].$$

(119)

The differentiation of F by the chain rule yields the equivalent quasi-linear form

$$\frac{\partial U}{\partial t}+A\frac{\partial U}{\partial x}=0, $$

(120)

where $A=\frac{\partial F}{\partial U}$ is the Jacobian matrix. It is easy to verify that

$$A=\left [ \begin{array}{c@{\quad}c@{\quad}c}0 & 1 & 0\\ \frac{1}{2}(\gamma-3)v^2 & (3-\gamma)v & \gamma-1\\[0.1cm]\frac{1}{2}(\gamma-1)v^3-vH\quad& H-(\gamma-1)v^2 &\gamma v\end{array} \right ]. $$

(121)

The eigenvalues and right/left eigenvectors of A satisfy the system of equations

$$A\mathbf{r}_k=\lambda_k\mathbf{r}_k,\quad \mathbf{l}_kA=\lambda_k\mathbf{l}_k,\quad k=1,2,3$$

(122)

which can be written in matrix form as AR=RΛ and R ⁻¹ A=ΛR ⁻¹. Thus,

$$A=R\varLambda R^{-1},\quad\varLambda=\mbox{diag}\{v-c,v,v+c\}$$

(123)

in accordance with (9). The matrices of eigenvalues and eigenvectors are given by

(124)

(125)

and

$$R^{-1}=\left [ \begin{array}{c@{\quad}c@{\quad}c}\frac{1}{2} (b_1+\frac{v}{c} ) &\frac{1}{2} (-b_2v-\frac{1}{c} ) & \frac{1}{2} b_2\\[0.1cm]1-b_1 & b_2v & -b_2\\[0.1cm]\frac{1}{2} (b_1-\frac{v}{c} ) &\frac{1}{2} (-b_2v+\frac{1}{c} ) & \frac{1}{2} b_2\end{array} \right ]=\left [\begin{array}{c} \mathbf{l}_1\\\mathbf{l}_2\\\mathbf{l}_3\end{array} \right ],$$

(126)

where

$$b_1=b_2\frac{v^2}{2},\qquad b_2=\frac{\gamma-1}{c^2}.$$

On a uniform mesh of linear finite elements, the coefficients of the lumped mass matrix M _L and of the discrete gradient operator C are given by

$$m_i=\varDelta x,\qquad c_{ij}=\left \{ \begin{array}{l@{\quad}l}1/2, &j=i+1,\\-1/2, &j=i-1.\end{array} \right .$$

(127)

The lumped-mass Galerkin approximation is equivalent to the central difference scheme which can be written in the generic conservative form

$$\frac{\mathrm {d}\mbox {\textsc {u}}_i}{\mathrm {d}t}+\frac{\mbox {\textsc {f}}_{i+1/2}-\mbox {\textsc {f}}_{i-1/2}}{\varDelta x}=0,$$

(128)

where

$$\mbox {\textsc {f}}_{i+1/2}=\frac{\mbox {\textsc {f}}_i+\mbox {\textsc {f}}_{i+1}}{2}.$$

The numerical flux of the low-order scheme with d _i+1/2 defined by (42) is

$$\mbox {\textsc {f}}_{i+1/2}=\frac{\mbox {\textsc {f}}_i+\mbox {\textsc {f}}_{i+1}}{2}- \frac{1}{2}|\mbox {\textsc {a}}_{i+1/2}|(\mbox {\textsc {u}}_{i+1}-\mbox {\textsc {u}}_i),$$

(129)

where a _i+1/2 is the 1D Roe matrix. The so-defined approximation is known as Roe’s approximate Riemann solver [54]. A detailed description of this first-order scheme can be found in many textbooks on gas dynamics [24, 37, 64]. Roe’s method fails to recognize expansion waves and, therefore, may give rise to entropy-violating solutions (rarefaction shocks) in the neighborhood of sonic points. Hence, some additional numerical diffusion may need to be applied in regions where one of the characteristic speeds approaches zero [20, 21]. This trick is called an entropy fix.

The use of scalar dissipation (46) leads to a Rusanov-like low-order scheme with

$$\mbox {\textsc {f}}_{i+1/2}=\frac{\mbox {\textsc {f}}_i+\mbox {\textsc {f}}_{i+1}}{2}-\frac {a_{i+1/2}}{2}(\mbox {\textsc {u}}_{i+1}-\mbox {\textsc {u}}_i), $$

(130)

where a _i+1/2 denotes the fastest characteristic speed. Zalesak [73] defines it as

$$a_{i+1/2}=\frac{|v_i|+|v_{i+1}|}{2}+\frac{c_i+c_{i+1}}{2}.$$

For reasons explained in [5], our definition of the Rusanov flux (130) is based on

$$a_{i+1/2}:=\max\{|v_{i}|+c_{i},|v_{i+1}|+c_{i+1}\}.$$

This formula yields a very robust and efficient low-order method for FCT [33].