Flux-conservative Hermite methods for simulation of nonlinear conservation laws

TL;DR

This paper introduces a novel class of Hermite methods for nonlinear conservation laws that maintain high-order accuracy and enhanced stability, incorporating entropy viscosity to effectively capture shocks.

Contribution

The paper presents a new Hermite method class that improves stability and shock capturing in nonlinear conservation law simulations.

Findings

Achieves high-order accuracy for smooth solutions.

Enhanced stability properties over existing methods.

Effective shock capturing using entropy viscosity.

Abstract

A new class of Hermite methods for solving nonlinear conservation laws is presented. While preserving the high order spatial accuracy for smooth solutions in the existing Hermite methods, the new methods come with better stability properties. Artificial viscosity in the form of the entropy viscosity method is added to capture shocks.

Tables

Table 1. Table 1 : Comparison of Costs in original and flux-conservative Hermite methods, n K subscript 𝑛 𝐾 n_{K} is the number of stages in Runge-Kutta method, d 𝑑 d is the spatial dimension.

	flux computation	interpolation
original	$n_{K}d{\left(2m+2\right)}^{2d}$	$2^{2^{d-1}}{\left(2m+2\right)}^{d+1}$
flux-conservative	$n_{K}d{\left(2m+2\right)}^{2d}$	$\left(1+n_{K}d\right)2^{2^{d-1}}{\left(2m+2\right)}^{d+1}$

Table 2. Table 2 : Convergence study of smooth solution to Burgers’ equation

$h_x$	$\pi /2$	$\pi /4$	$\pi /8$	$\pi /16$	$\pi /32$
$l_{\infty }$-err $m=1$	2.30(-01)	5.85(-02)	1.09(-02)	1.42(-03)	1.80(-04)
Rate		1.97	2.43	2.94	2.98
$l_{\infty }$-err $m=2$	4.85(-02)	5.47(-03)	2.19(-04)	7.25(-06)	1.97(-07)
Rate		3.15	4.64	4.92	5.20
$l_{\infty }$-err $m=3$	1.09(-02)	6.59(-04)	7.71(-06)	4.70(-08)	2.73(-10)
Rate		4.04	6.42	7.36	7.43

Table 3. Table 3 : Convergence study on Euler’s equation with stationary shock.

$h_x$	1.25(-02)	6.25(-03)	3.13(-03)	1.56(-03)	7.81(-04)	3.91(-04)
$L_1$ error	3.44(-03)	2.29(-03)	1.40(-03)	7.96(-04)	4.46(-04)	2.49(-04)
Rate		0.59	0.71	0.81	0.83	0.84

Table 4. Table 4 : Parameters for examples in 1D Euler’s equations.

Problem	CFL	$m$	${\alpha }_{EV}$	${\alpha }_{max}$
Lax	0.2	3	0.5	0.08
Sod	0.15	3	0.2	0.08
Shu-Osher	0.15	3	0.01	0.05
Stationary shock	0.2	3	10	0.3

Table 5. Table 5 : Parameters for examples in 2D Euler’s equations.

Problem	$CFL$	$m$	${\alpha }_{EV}$	${\alpha }_{max}$
Explosion/implosion	0.2	3	0.1	0.2
Vortex-shock interaction 1	0.2	3	0.01	0.04
Vortex-shock interaction 2	0.2	3	0.05	0.07
Jet	0.2	3	0.03	0.2

Equations

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

$$\frac{d}{dt}{\int }_a^{b}udx=f(u(x(a),t))-f(u(x(b),t)).$$

$${\int }_a^b\frac{\partial u_h^j}{\partial t}dx$$

$$=\sum _{j=1}^{N_x}{\int }_{x_{j-1}}^{x_j}\frac{\partial }{\partial x}(-f_h^j)dx$$

$$=\sum _{j=1}^{N_x}\left(f_h^j{|}_{x=x_{j-1}}-f_h^j{|}_{x=x_j}\right)$$

$$=f_h^j(u(a))-f_h^j(u(b)).$$

$${\mathcal{G}}^p=\left\{x_j=x_L+j{h}_x\right\},j=0,\dots ,N_x,$$

$${\mathcal{G}}^d=\left\{x_{j+1/2}=x_L+\left(j+\frac{1}{2}\right)h_x\right\},j=0,\dots ,N_x-1,$$

$$u^p(x,t_n)=\sum _{k=0}^m{c}_k(t_n){\left(x-x_j\right)}^k,x\in {\mathcal{I}}_j^d=[x_{j-1/2},x_{j+1/2}],$$

$$u^d(x,t_n)=\sum _{k=0}^{2m+1}{c}_k(t_n){\left(x-x_{j-1/2}\right)}^k,x\in {\mathcal{I}}_j^p=[x_{j-1},x_j],$$

$$\frac{{\partial }^k}{\partial x^k}{u}^d(x_i,t_n)=\frac{d^k}{d{x}^k}{u}^p(x_i,t_n),k=0,\dots ,m,i=j-1,j.$$

$$(u_h{)}_t+(f_h{)}_x$$

$$(u_h{)}_{tx}+(f_h{)}_{xx}$$

$$(u_h{)}_{txx}+(f_h{)}_{xxx}$$

$$⋮$$

$$f_h=\sum _{k=0}^{2m+1}{b}_k(t)(x-x_{j-1/2}{)}^k.$$

$$\left(\begin{array}{c}{c}_0^{\prime }(t)\\ c_1^{\prime }(t)\\ ⋮\\ c_{2m}^{\prime }(t)\\ c_{2m+1}^{\prime }(t)\end{array}\right)=\left(\begin{array}{c}{b}_1(t)\\ 2b_2(t)\\ ⋮\\ (2m+1)b_{2m+1}(t)\\ 0\end{array}\right).$$

$$u^p(x,t_n)=\sum _{k=0}^{2m+1}{c}_k(t_n){\left(x-x_j\right)}^k,x\in {\mathcal{I}}_j^d=[x_{j-1/2},x_{j+1/2}].$$

$$F_1^p$$

$$F_2^p$$

$$F_3^p$$

$$F_4^p$$

$$w^d(x,t_n)=\sum _{k=0}^{2m+1}{d}_k(t_n){\left(x-x_{j-1/2}\right)}^k,x\in {\mathcal{I}}_j^p=[x_{j-1},x_j],$$

$$\frac{{\partial }^k}{\partial x^k}{w}^d(x_i,t_n)=\frac{d^k}{d{x}^k}{w}^p(x_i,t_n),k=0,\dots ,m,i=j-1,j.$$

$$c_k(t_n+\Delta \hat{t})=c_k(t_n)+\frac{\Delta \hat{t}}{6}(k+1)(b_{k+1}^{(1)}(t_n)+2b_{k+1}^{(2)}(t_n)+2b_{k+1}^{(3)}(t_n)+b_{k+1}^{(4)}(t_n)).$$

$$u^d(x,t_n+\Delta \hat{t})=\sum _{k=0}^{2m+1}{c}_k(t_n+\Delta \hat{t}){\left(x-x_{j-1/2}\right)}^k,x\in {\mathcal{I}}_j^p=[x_{j-1},x_j].$$

$$c_k(t_n+\Delta \hat{t})-c_k(t_n)+\Delta \hat{t}\sum _{i=1}^s{\alpha }_i(k+1)b_{k+1}^{(i)}(t_n)=0,k=0,\dots ,2m+1,$$

$$\sum _j{\int }_{x_{j-1}}^{x_j}\frac{{\partial }^k}{\partial x^k}{u}^d(t_n+\Delta \hat{t},x)dx=\sum _j{\int }_{x_{j-1}}^{x_j}\frac{{\partial }^k}{\partial x^k}{u}^d(t_n,x)dx,k=0,\dots ,2m+2-s.$$

$$\frac{u^d(t_{n+1/2})-u^d(t_n)}{\Delta \hat{t}}=\Delta \hat{t}\sum _i{\alpha }_i\frac{d{F}_i^d}{dx},$$

$$u_t+(f(u){)}_x+(g(u){)}_y=0,$$

$${\mathcal{G}}^p=\{(x_i,y_j)\}=(x_L+i{h}_x,y_B+j{h}_y),i=0,\dots ,N_x,j=0,\dots ,N_y,$$

$$\begin{gathered} {\mathcal{G}}^d=\left\{(x_{i+1/2},y_{j+1/2})=\left(x_L+\left(i+\frac{1}{2}\right)h_x,y_B+\left(j+\frac{1}{2}\right)h_y\right)\right\}, \\ i=0,\dots ,N_x-1,j=0,\dots ,N_y-1, \end{gathered}$$

$$\frac{u^d(t_n+\Delta \hat{t})-u^d(t_n)}{\Delta \hat{t}}=\frac{K_1^d+2K_2^d+2K_3^d+K_4^d}{6}.$$

$$K_1^d$$

$$K_2^d$$

Full text

Keywords— Hermite, conservative, high order method, the entropy viscosity method

1 Introduction

Conservation laws model physical systems arising in traffic flows, aircraft design, weather forecast, and fluid dynamics. Numerical methods for conservation laws ideally conserve quantities like mass or energy, and accurately capture various physical components of the solutions, from small smooth scales to shock waves. The presence of both smooth waves and shock waves, for example in shock-turbulence interaction creates a challenge to the simulation of nonlinear conservation laws.

Shock waves typically appear in solutions to nonlinear conservation laws, and are characterized by very thin regions where the solution changes rapidly. Approximation of shocks and small waves is challenging as the small and large scales need to be solved simultaneously. Historically, low order finite volume and finite difference methods equipped with flux/slope limiters have been used to handle shocks, see for example the textbooks [16, 15]. The drawback with low order methods is that they cannot accurately propagate small scales over long distances and as a result, today the research focus has gravitated towards high order accurate methods with shock capturing capability.

Among high order methods, the weighted essentially non-oscillatory (WENO) method, [23, 22, 17], has proven to be a method popular among practitioners. WENO methods are still relatively dissipative which may be a drawback for turbulent simulations, [18]. Also, discontinuous Galerkin methods combined with shock capturing, [5, 14], or selectively added artificial viscosity, [20, 12, 26], have received significant interest. The latter approach traces back to the artificial viscosity method by Neumann and Richtmyer, [19] and the popular streamline diffusion method, [3, 11], which computes the viscosity based on the residual of the PDE.

In this work, we advocate the combination of a high order method and selectively added artificial viscosity. Specifically, we show how the entropy viscosity by Guermond and Pasquetti, [7], can be implemented in our new flux-conservative Hermite methods.

First introduced by Goodrich, Hagstrom, and Lorenz in [6] for hyperbolic initial-boundary value problems, Hermite methods use the solution and its first $m$ derivatives in each coordinate to construct an approximate solution to the PDE. The formulation by Goodrich et al. computes the flux at the cell centers, which for nonlinear problems leads to discontinuous fluxes at cell interfaces. This discontinuity results in lost of conservation.

To address the lack of conservation, we develop a new class of Hermite methods, which share the basic features with the method in [6], such as interpolation and time evolution, but differs in the computation of the flux function. In the flux-conservative Hermite methods, the numerical flux is computed at cell edges and then interpolated to cell center for time evolution, hence by construction, the numerical flux is continuous at cell interface. Additionally, for nonlinear problems, it is more efficient to use one step methods than the Taylor series approach in [6], see [9, 10]. Here we will use the standard Runge-Kutta method to evolve in time.

The rest of the paper is organized as follows. In Section 2, we derive conservation laws and discrete conservation. Then, in Section 3, we describe Hermite methods as first introduced by Goodrich et al. [6], followed by the description of the flux-conservative Hermite methods in Section 3.2. For shock capturing capability, we implement the entropy viscosity method, which is explained in Section 4. In Section 5, we present the results of simulation on Euler’s equations, see [13] for results on Burgers’ equation.

2 Conservation Laws

A scalar conservation law in one space dimension takes the form

[TABLE]

where $u(x,t)$ is the state variable at location $x$ and time $t$ and $f(u)$ is the flux, or the rate of flow, of the state variable $u$.

The derivation of conservation laws comes from the observation that at any given time $t$, the rate of change of the total quantity of the state variable $u$ over some interval $[a,b]$ must be equal to the net flux $f(u)$ into the interval through the endpoints. Mathematically, this can be expressed as

[TABLE]

When approximating the solution to scalar conservation laws given by equation (1), the PDE is typically discretized on a grid consisting of $N_{x}$ cells where $x_{0}=a$ and $x_{N_{x}}=b$. It is desirable that the numerical method satisfies discrete conservation. If the reconstructed solution $u_{h}^{j}$ and flux $f_{h}^{j}$ on any cell $\mathcal{I}_{j}=[x_{j-1},x_{j}]$ satisfy the condition $f_{h}^{j}|_{x=x_{j}}=f_{h}^{j+1}|_{x=x_{j}}$, $j=1,\ldots,N_{x}-1$, we immediately find

[TABLE]

The property that $f_{h}^{j}|_{x=x_{j}}=f_{h}^{j+1}|_{x=x_{j}}$ does not hold for the original Hermite methods, and our goal here is to design a new Hermite method with this property. Before describing our new method, we briefly describe the original method.

3 Hermite Methods

A Hermite method of order $2m+1$ approximates the solution to a PDE by an element-wise polynomial that has continuous spatial derivatives up to order $m$ at the element’s interfaces. In Hermite methods, the degrees of freedom are function and spatial derivative values, or equivalently the coefficients of the Taylor polynomial at the cell center of each element. The evolution of the degrees of freedom on each element can be performed locally.

3.1 Hermite Method in One Dimension

Consider again equation (1) on the domain $D=[x_{L},x_{R}]$. Let $\mathcal{G}^{p}$ and $\mathcal{G}^{d}$ be the primal grid and the dual grid, defined as

[TABLE]

[TABLE]

where $h_{x}=(x_{R}-x_{L})/N_{x}$ is the distance between two adjacent nodes. Let $u^{p}$ and $u^{d}$ be the approximations to the solution on the primal and dual grids, respectively.

At the initial time $t_{n}=t_{0}+n\Delta t$, we assume that the approximation $u^{p}$ on the primal grid, the global piecewise polynomial

[TABLE]

is known. Starting from time $t=t_{n}$ on the primal grid $\mathcal{G}^{p}$, we evolve the solution one full time step by:

Reconstruction by Hermite interpolation: We construct $u^{d}$, the global Hermite interpolation polynomial of degree $(2m+1)$, on the dual grid. That is,

[TABLE]

where the coefficients $c_{k}(t_{n})$ are uniquely determined by the interpolation conditions

[TABLE]

Time evolution: By rewriting equation (1) as $u_{t}=-f(u)_{x}$, it is obvious that in order to evolve $u^{d}$, we need to compute a polynomial approximation $f^{d}$ to the flux function $f(u)$. We offer two ways to obtain $f^{d}$:

• •

Modal approach: Perform polynomial operations (addition, substraction, multiplication, and division) on $u^{d}$ and truncate the resulting polynomial to suitable degree.

• •

Pseudospectral approach: Compute a local polynomial $f_{h}^{*}$ interpolating $f(u^{d})$ on a quadrature grid $\mathcal{G}^{ps}$ inside $\mathcal{I}^{p}_{j},j=1,\ldots,N_{x}$, and transform $f_{h}^{*}$ to a Taylor polynomial $f^{d}$.

We differentiate $f^{d}$ in polynomial sense to get an approximation to the derivative of the flux function $f(u)_{x}$. We usually use the modal approach unless this option is not applicable. Next, we derive an ODE to evolve $u^{d}$, or equivalently the coefficients of the Hermite interpolant ${\mathbf{c}(t)=(c_{0}(t),\ldots,c_{2m+1}(t))^{T}}$, by insisting that the numerical solution $u^{d}$ satisfy equation (1) along with derivatives of (1), at the cell centers $x=x_{j-1/2},j=1,\ldots,N_{x}$. The resulting system of ODE for $c_{k},k=0,\ldots,2m+1$, can then be evolved independently on each $\mathcal{I}^{p}_{j}$ with any ODE solver. The reconstruction step provides the initial data, $\mathbf{c}(t_{n})$.

Repeat on dual grid: At the end of the half time step, we have evolved the degree $2m+1$ polynomial $u^{d}$ from time $t=t_{n}$ to $t=t_{n+1/2}$. Before repeating the above process, we truncate $u^{d}(x,t_{n+1/2})$ by removing the coefficients $c_{k}$, $k=m+1,\ldots,2m+1$. We then repeat the above process (including the truncation) to obtain $u^{p}$ at time $t=t_{n+1}$, see Figure 1 for illustration.

3.1.1 Example: Burgers’ Equation

To illustrate the specifics of the time evolution, we consider Burgers’ equation, with $f(u)=u^{2}/2$, approximated by $f_{h}=u_{h}^{2}/2$, where $u_{h}$ represents the degree $(2m+1)$ interpolant on either of the grids. We can write

[TABLE]

where

[TABLE]

The coefficients $b_{k}$ are determined by truncated polynomial multiplication, that is $b_{k}=\frac{1}{2}\sum_{l=0}^{k}c_{l}c_{k-l}$. Insisting that the approximate solution $u_{h}$ satisfy equation (9) at the cell centers $x=x_{j-1/2}$, we obtain

[TABLE]

While equation (10) is valid for any cell, the initial data for each cell are different from one cell to another. For a more detailed explanation and open source implementations, see [9, 10, 2].

3.2 Flux-Conservative Hermite Methods

The numerical flux $f_{h}$ obtained by the approach described above, is discontinuous at cell interfaces when the flux function $f(u)$ is nonlinear. Numerically, the discontinuity in the flux induces numerical stiffness. As a result, the time step often needs to be taken very small. To remedy this, we propose new flux-conservative Hermite methods that impose flux continuity by computing the numerical flux at cell edges, and then interpolate the numerical flux to cell center.

To illustrate the difference between the original and flux-conservative Hermite schemes, we plot the numerical flux $f_{h}=u_{h}^{2}/2$ with $m=3$ and for $N_{x}=3$ cells in Figure 2. The numerical flux obtained using the original Hermite method, shown as the blue curve, has discontinuities at cell interfaces. On the other hand, the numerical flux obtained by the flux-conservative Hermite method, shown as the black curve, is continuous everywhere.

3.2.1 The Construction of the New Method

The goal of our construction is to globally conserve the integral of $u_{h}$ and its $m$ first derivatives over one half time step, with $\Delta\hat{t}=\Delta t/2$.

In the flux-conservative methods, we assume that the solution on the primal grid at time $t_{n}$ is given by

[TABLE]

Note that the degree of this polynomial is different than in the original method. We assume that the time stepping will be performed by an explicit one step method requiring stage values. The evolution will be carried out at the cell center where the stage values will be the derivative of the Hermite interpolant of the flux. As this interpolant is $m$ times differentiable at the edges, it will result in a conservative discretization. Now, the time evolution of the approximate solution entails the following steps.

Computation of the stage fluxes at the cell edges: For simplicity, assume that we use the classic fourth order Runge-Kutta, then to construct the Hermite interpolants for the four stages we first compute

[TABLE]

Note that inside the argument of $f$, we keep all the coefficients of the polynomials up to degree $2m+1$, while the nonlinearity $f$ itself, which typically is a higher degree polynomial, is truncated to degree $2m+1$.

Reconstruction of solution and fluxes by Hermite interpolation: Next, we construct $u^{d}$ and $F_{i}^{d}$, $i=1,\ldots,4$, the global Hermite interpolation polynomials of degree $(2m+1)$ of the solution and the flux, respectively. Let $w^{d}$ represent $u^{d}$ or $F_{i}^{d}$ and $w^{p}$ represent $u^{p}$ or $F_{i}^{p}$. Then,

[TABLE]

where the coefficients $c_{k}(t_{n})$ are uniquely determined by the interpolation conditions (at cell edges)

[TABLE]

Time evolution: Let the coefficients of $u^{d}$ be $c_{k}$ and the coefficients of $F_{i}^{d}$ be $b_{k}^{(i)}$. Then, again assuming RK4, we have that for $k=0,\ldots,2m+1$,

[TABLE]

The updated solution on the dual grid is thus

[TABLE]

Repeat on dual grid: At the end of the half time step, we have evolved the degree $2m+1$ polynomial $u^{d}$. We then repeat the above process to obtain $u^{p}$ at time $t=t_{n+1}$.

We note that unlike the original method, the number of degrees of freedom that we keep is twice as large, representing an increase in memory requirement. The number of floating point operations, however, to leading order, is the same as for the original method (see the complexity analysis below).

3.2.2 Conservation

We now consider the conservation properties of the above scheme. Due to the fact that the $s$ fluxes used in the stages have $m$ continuous derivatives we immediately find that for periodic problems the following conservation statements hold.

Theorem: Assume we use the flux-conservative Hermite method to evolve $u_{t}+f(u)_{x}=0$ with periodic boundary conditions. Further assume that $u^{d}(t,x)$ is the periodic degree $2m+1$ Hermite interpolating polynomial and that $F^{d}_{i}$, $i=1,\ldots,s$ are the periodic degree $2m+1$ polynomials Hermite interpolating the fluxes. Further, let the coefficients $c_{k}(t)$ of $u^{d}$ on a cell be evolved from time $t=t_{n}$ to $t=t_{n}+\Delta\hat{t}$ by the one step method

[TABLE]

where $b_{k}^{(i)}$ are the coefficients of $F^{d}_{i}$. Then, the following conservation statement holds.

[TABLE]

Proof: From the RK time stepping method for conservation laws (1) as

[TABLE]

together with the continuity of $m$ first derivatives of each of the $F^{d}_{i}$, the result follows immediately from the update formula. Note also that $F^{d}_{i+1}$ is one order less accurate than $F^{d}_{i}$ due to flux differentiation during stage $i$.

To summarize, in the original Hermite scheme, computation of numerical fluxes is performed at the cell center using the interpolated solution. The flux-conservative Hermite scheme requires a computation (and storage) of numerical fluxes at the cell edges and the interpolation of those fluxes to the cell center. Refer to Figure 3 for an illustrative comparison between the schemes.

3.3 The Flux-Conservative Hermite Method in Two Dimensions

Now, let us consider a conservation law

[TABLE]

on the domain $D=[x_{L},x_{R}]\times[y_{B},y_{T}]$. Let $\mathcal{G}^{p}$ and $\mathcal{G}^{d}$ be the primal and dual grids, defined as

[TABLE]

[TABLE]

where $h_{x}=(x_{R}-x_{L})/N_{x}$ and $h_{y}=(y_{T}-y_{B})/N_{y}$ are the distances between two adjacent nodes in $x$ and $y$ directions respectively.

The extension of the flux-conservative method from one dimension is straightforward. Writing equation (16) as $u_{t}=-f(u)_{x}-g(u)_{y}$ and letting $u^{d}(t)$ represent the two-dimensional tensor product Hermite interpolant of the data on the primal grid we can write the RK4 time stepping of $u^{d}(t_{n})$ to time $t=t_{n}+\Delta\hat{t}$ as

[TABLE]

The left hand side of equation (19) is an approximation to $u_{t}$ and the right hand side is an approximation to stage values of $-(f(u)_{x}+g(u)_{y})$. Similar to the one dimensional case we have

[TABLE]

Here, for example, $F_{i}^{d}$ is the degree $2m+1$ tensor product polynomial that interpolates $f(u^{p}+\gamma_{i}\Delta\hat{t}F^{p}_{i})$ and its $m$ first derivatives in $x$ and $y$ at the four adjacent primal grid-points.

3.4 Comparison of Computational Costs

The time evolution portion of the Hermite methods are performed by a one step method with $n_{K}$ stages, involving computation of the flux function, interpolation of the solution and, in the flux-conservative method, the fluxes, and differentiation of fluxes. For the purpose of this comparison, we assume Burgers’ flux function $f(u)=u^{2}/2$ in 1D or $f(u)=g(u)=u^{2}/2$ in 2D. Each 1-dimensional Hermite interpolation is equivalent to a multiplication by a $(2m+2)$ by $(2m+2)$ matrix. If we use the recipe above, each 2-dimensional Hermite interpolation corresponds to ${2\times(2m+2)}$ one-dimensional interpolations. The factor $(2m+2)$ comes from the the fact that the $y$ dimension brings in $(2m+2)$ interpolations in 1D and the multiplicative factor $2$ comes from the fact that we interpolate in $y$ direction on both the left and right edges of the cell. In 3D, we have 4 interpolations in the $z$ direction, each brings in $(2m+2)$ times interpolations in 2D, and so on. We summarize the cost of the method, corresponding to the number of multiplications involved, in Table 1. The number of interpolation in the flux-conservative Hermite method is $n_{K}d\;2^{2^{d-1}}$ more than the original Hermite method. We note that the flux-conservative scheme can be simplified to just two flux interpolations by adding up the $F$’s and $G$’s, but in this paper, we interpolate each flux separately. There is also an additional cost of differentiation at cell corners, but it is negligible compared to the cost of interpolation.

4 The Entropy Viscosity Method

Given a PDE on the form (1), there exists an entropy function $E(u)$ and its corresponding entropy flux function $F(u)=\int E^{\prime}(u)f^{\prime}(u)\;du$ such that the entropy residual satisfies

[TABLE]

This inequality can be used to select the physically correct solution to (1) or (16). The direction of the inequality can vary from one problem to another, but the inequality becomes strict only at shocks. In essence, the entropy viscosity (EV) method uses the fact that the residual approaches a Dirac delta function centered at shocks to construct a nonlinear dissipation. The resulting dissipation is small away from shocks and just sufficient amount at a shock. The details of EV for conservation laws are described in detail in [8] but we summarize its most important features here.

Consider the conservation law whose right hand side has been replaced by a viscous term, ${u_{t}+\nabla\cdot{f}(u)=\nabla\cdot(\nu\nabla u)}$, with $\nu=\nu_{h}(x,t)$ given by

[TABLE]

where $\nu_{EV}$ is the entropy-based viscosity and $\nu_{\max}$ is a viscosity whose size depends on the largest eigenvalue of the flux function $f(u)$, representing the maximum wave speed. The discretized entropy-based viscosity $\nu_{E}$ is then given by

[TABLE]

where $\beta$ is a positive scalar, $\alpha_{EV}$ is a user defined constant and $C_{1}(u_{h})$ is some PDE-specific normalization.

At shocks, the entropy residual approaches Dirac delta function, so we instead use

[TABLE]

where $\alpha_{\max}$ is another user defined constant, $C_{2}$ serves as the maximum wave speed and $V_{x}$ is some neighborhood of $x$. The $V_{x}$ neighborhood can either be “local”, i.e. containing only a few cells around $x$, or “global”, i.e. $V_{x}=D$, where $D$ is the domain of the PDE. In this work, we use global $V_{x}$.

In recent papers, the parameter $\beta$ is chosen to be 2, but we found that this may prevent convergence for moving shocks, see [13] where we also argue that $\beta=1$ is a more appropriate choice. In essence our argument is as follows. As the entropy residual approaches a Dirac delta distribution, a consistent discretization of the residual with a single shock must satisfy

[TABLE]

where $h_{j}=x_{j+1}-x_{j}$. Thus, we expect $(r_{EV})_{j}\sim h_{j}^{-1}$ near the shock. When $\beta=2$, the two terms $\nu_{EV}$ and $\nu_{\max}$ are both $\mathcal{O}(h)$. Since the parameters are tuned on a coarse grid and the terms $C_{1}$ and $C_{2}$ in (22)-(23) also changes with the grid size, the selection of the minimum of $\nu_{EV}$ and $\nu_{\max}$ does not necessarily “converge” as the grid gets refined. If instead we choose $\beta=1$, then $\nu_{EV}=\mathcal{O}(1)$ while $\nu_{\rm max}=\mathcal{O}(h_{j})$, and the particular choice of ${\alpha}_{EV}$ is thus irrelevant in the limit $h_{j}\to 0$ as the selection mechanism will eventually select $\nu_{\rm max}$ at the shocks.

While the explicit formula for $C_{1}$ and $C_{2}$ varies from one PDE to another, the core of the entropy viscosity method remains the same. The size of the entropy residual gives us a sense of relative distance to the shock, which is then used to take the following actions: near a shock, EV uses sufficiently large dissipation, $\nu_{h}=\nu_{\max}$, and away from a shock, EV uses entropy-based dissipation, $\nu_{h}=\nu_{EV}$.

5 Numerical Examples

We start by confirming that the rates of convergence, $2m+1$, (in space) are the same for the new flux-conservative method as for the original method.

5.1 Convergence for a Smooth Solution

We solve Burgers’ equation on the domain $x\in[-\pi,\pi]$ and impose periodic boundary conditions. The initial data is $u(x,0)=-\sin(x)+0.3$ and we evolve the solution until time $t=0.4$. The timestep is chosen as $\Delta t={\rm CFL}\,\,h_{x}/$ $\max_{x}|u(x,0)|$, with ${\rm CFL}=0.1$.

In Table 2, we display the maximum error at the final time computed against a reference solution computed using $m=5$ and $h_{x}=\pi/64$. As can be seen the rates of convergence appear to approach the predicted rate $2m+1$.

We next present a sequence of experiments displaying the performance of the Hermite-Runge-Kutta-4-Entropy-Viscosity method for Euler’s equations (with artificial viscosity).

5.2 Euler’s equations in One Dimension

We consider Euler’s equations which represent conservation of mass, momentum, and energy,

[TABLE]

Here, $\rho$ is the density, $\rho u$ is the momentum, $u$ is the velocity, and $E$ is the energy. Furthermore, we assume an ideal gas with the equation of state

[TABLE]

where $\gamma=1.4$ is the adiabatic index and $p$ is the pressure.

To regularize equation (25), we add a viscous term $(\nu(\rho,\rho u,E)_{x})_{x}$, where the coefficient $\nu$ is obtained using the entropy viscosity method. Thus, the viscous Eulers’ equations can be written as

[TABLE]

We note that an alternative to this simple viscosity would be to use the full Navier-Stokes equations.

5.2.1 Entropy Viscosity (EV) method for 1D Euler’s equations

The discretized viscosity coefficient $\nu=\nu_{h}$ is given in terms of primitive variables $\rho,p,\text{and }u$,

[TABLE]

where $T_{h}=p_{h}/\rho_{h}$ is the temperature, $h_{x}$ is the grid size and

[TABLE]

is the entropy residual for the entropy function $S_{h}(p_{h},\rho_{h})=\frac{\rho_{h}}{\gamma-1}\log\left(\frac{p_{h}}{\rho_{h}^{\gamma}}\right)$ and its corresponding entropy flux $(uS)_{h}$.

5.2.2 An Improved Entropy Viscosity

The entropy viscosity method discretizes the entropy residual using the numerical solution. In theory, the entropy residual is large at shocks, and zero at contact discontinuities (where no artificial viscosity is needed). However, our experience is that the discretization of the entropy equation may also trigger the maximum viscosity at contact discontinuities. To the left in Figure 4, we see a space-time diagram of the entropy residual for Sod’s problem in logarithm scale. Note that a relatively large amount of residual is produced at the contact discontinuity.

To eliminate this undesired behavior along the contact discontinuity, we use the fact that the velocity of the fluid, $u$, is a Riemann invariant along the second characteristic field. Since $u$ is smooth at the contact discontinuity but not at a shock, the product of $\Delta u_{j}$ and $(r_{E})_{j}$ is small at contact discontinuities but still large at shocks. We incorporate the term $\Delta u$ into the improved entropy viscosity

[TABLE]

To this end we take $\nu_{E}$ and $\nu_{max}$ as a piecewise constant function on each cell. Thus, we compute the discretized density $\rho_{h}$, velocity $u_{h}$, temperature $T_{h}$, entropy function $S_{h}$ and entropy flux $uS_{h}$ at cell center in pointwise manner. Now, to get the entropy residual given in (31), we compute temporal and spatial derivatives using finite differences. Using the notation $S_{h}=S_{j}^{n}$ to denote the approximate flux function $S$ at $x=x_{j}$, $t=t_{n}$, we discretize the term $\partial_{t}S_{j}^{n}$ using second order Backward Difference formula

[TABLE]

Similarly, the term $\partial_{x}(uS)_{j}^{n}$ is approximated by the centered finite difference

[TABLE]

5.3 Experiments in One Dimension with Euler’s Equations

We now present results obtained using our Hermite-RK4-EV method for a stationary shock, the Lax, the Sod, and the Shu-Osher problem. For experiments where we use more than one resolution, the EV parameters are tuned on the coarsest grid. In the 1D Euler’s equations simulations, the timestep is chosen as $\Delta t={\rm CFL}\,h_{x}/\max_{x}|(u\pm c)(x,0)|$, where $c=\sqrt{\gamma p/\rho}$ is the speed of sound, with ${\rm CFL}$ values given in Table 3.

5.3.1 Stationary shock problem

By solving the Riemann problem, we decide the states corresponding to a stationary shock. The goal of this experiment is to investigate the stability and accuracy of EV in the presence of shocks. Since small oscillations coming from shocks could potentially pollute the “smooth” regions, this is also a test for how well EV removes numerical artifacts. The computational domain is $D=[-0.5,0.5]$ with the stationary shock given by

[TABLE]

The boundary condition are imposed by setting the solution at the boundary so that it coincides with the solution at initial time. We perform a grid refinement study and report the errors in the density in Table 3. We also present the ratio between successive errors. It appears that the rate of convergence for the $L_{1}$ error is approaching 7/8.

5.3.2 Lax’s and Sod’s Shock Tube Problems

Lax’s and Sod’s problems come from physical experiments in which a gas tube is separated by a membrane into two sections. The gas in each section is uniform in the $y$ and $z$ direction, so the problem is modeled as a 1-dimensional shock tube. The gas in the left section is kept at a different state than the gas in the right section. At time $t=0$, the membrane is punctured. In the problem setup, the Euler’s equations are solved on the domain $D=[-0.5,0.5]$ with initial data

[TABLE]

for Lax, and

[TABLE]

for Sod. For both problems, we impose fixed boundary condition so that the solution on the boundary is the same as at the initial time. The solution is computed up to time $t=0.16$ for Lax’s problem and time $t=0.1644$ for Sod’s problem.

The solution to Riemann problems such as Lax’s and Sod’s shock tubes contains 3 waves propagating from the discontinuity at the initial time. The second wave is a contact discontinuity, where the discontinuity is translated over time. The first and third waves are nonlinear, and can take either rarefaction waves or shock waves.

The results for density $\rho$, velocity $u$ and pressure $p$ are plotted against the exact solution in Figure 6. In each plot, we use $N_{x}=100$ elements. The entropy viscosity parameters used are given in Table 4. In both problems, the shocks are resolved better than the contact discontinuities. Although the shock strength is only of medium size for both problems, some experts considered these tough test cases for non-characteristic-based high order schemes [23].

5.3.3 The Shu-Osher Problem

The Shu-Osher problem poses difficulties for numerical methods due to the sinusoidal interacting with the shock. Here the domain is $D=[-5,5]$ with initial data

[TABLE]

and with fixed boundary condition so that the solution on the boundary coincides with the solution at the initial time. The solution is computed up to time $t=1.8$ and compared against a computed solution on a much finer grid. We use $N_{x}=80$ to obtain the numerical solution in Figure 7, where we interpolate the solution on to a finer grid. The “exact” solution is computed on a grid with $N_{x}=1280$. Note that even if we use a coarse grid, we can still get roughly the shape of the solution, especially away from the shock. However, when smooth waves are present (see blue oscillatory line to the left of shock) and too close to the shock, these waves get damped.

5.4 Euler’s equations in Two Dimensions

The two dimensional viscous Euler equations are given by

[TABLE]

Here, $\rho$ is the density, $\rho u$ and $\rho v$ are the momentum, $u$ and $v$ are the velocity in $x$ and $y$ direction respectively and $E$ is the energy. Furthermore, we assume an ideal gas with equation of state

[TABLE]

where $\gamma=1.4$ is the adiabatic index and $p$ is the pressure. For all experiments below, the timestep is chosen as ${\Delta t={\rm CFL}\,\,h_{x}/\max_{x}|(u\pm c)(x,0)|}$, where $c=\sqrt{\gamma p/\rho}$ is the speed of sound, with ${\rm CFL}$ values given in Table 5.

5.5 Entropy Viscosity Method for Euler’s Equations in Two Dimensions

The entropy viscosity is identical to the 1D version given in (28), with the exception that it takes the velocity in both directions into account.

[TABLE]

where $T_{h}=p_{h}/\rho_{h}$ is the temperature, $h=\min(h_{x},h_{y})$ is the grid size and

[TABLE]

To discretize the entropy residual $r_{EV}$, we again use BDF for the time derivative and centered finite differences for the spatial derivatives,

[TABLE]

On the domain $[x_{L},x_{R}]\times[y_{B},y_{T}]$, we use the subscript $jk$ to indicate that the variable attached is evaluated at $x=x_{L}+jh_{x}$ and $y=y_{B}+kh_{y}$.

5.5.1 Explosion/Implosion Test

First we solve a radially symmetric Riemann problem from Toro [24]. The computational domain is $D=[-1,1]\times[-1,1]$, and the initial data corresponding to an expanding wave is

[TABLE]

For an imploding wave, the initial data is,

[TABLE]

The boundary conditions are imposed by setting the solution on the boundary so that it stays the same as the solution at the initial time. The simulation is performed until time $t=0.25$, before any waves reach the boundary of the domain. We plot the 2D solution in Figure (8). In Figure (9), we present a cross section of the density at time $t=0.25$ with $N_{x}=N_{y}=100$ elements against computed “exact” solution obtained with $N_{x}=N_{y}=400$ elements.

5.5.2 Shock Vortex Interaction

Next we consider the interaction of a shock and a vortex. In general shock-vortex interactions can produce small scales in the form of acoustic waves, and other interesting wave patterns. It has received substantial interest in the literature, see for example [21, 25, 4, 5].

In this experiment, a strong stationary shock with Mach number $2/\sqrt{1.4}\approx 1.69$ collides with a weak vortex with a Mach number $6/2\pi\approx 0.81$. The computational domain is $D=[-9,3]\times[-4,4]$ and the initial data is

[TABLE]

where

[TABLE]

and $\beta=6$.

As the vortex passes through the shock, the shock is distorted and the vortex is compressed into an elliptical shape. This phenomena is due to the fact that the vortex is relatively weak compared to the shock. The results are consistent with [25]. In Figure 10, we compare snapshots of the density Schlieren using two different sets of entropy viscosity parameters, see Table 5. Although the schlierens are plotted on the same color scale, notice that the structures are more pronounced in the pictures on the right column.

5.5.3 Fluid Flow in Jet

As a final example we simulate a planar Mach $2/\sqrt{1.4}\approx 1.69$ jet. The domain $D=[-15,55]\times[-17.5,17.5]$ is discretized using of $500\times 250$ cells. The initial data is

[TABLE]

We model the jet nozzle by a simple momentum forcing over a $1\times 1$ patch at the left edge of the computational. The jet is started impulsively causing a relatively strong compression to be generated. This wave sharpens up to a shock wave that is handled by the entropy viscosity as it is propagated from the nozzle and out into a damping absorbing layer of super-grid type, see [1].

In Figures (11)-(13) we display snapshots of the vorticity, dilatation and density Schlieren. We note that the viscosity we use here is purely for the regularization of shocks so there is no reason to believe that the flow that we compute resembles reality. Nevertheless, the example illustrates the new methods ability to handle rapidly started flows. Also, it is likely that the particular form of the artificial viscosity does not effect the robustness of the method.

6 Conclusion

In conclusion we have demonstrated that flux-conservative Hermite methods are suitable for solving nonlinear conservation laws, especially in the presence of shocks. The new methods still converges at a rate of $(2m+1)$ for smooth problems.

The adaptation of the entropy viscosity method to Hermite methods successfully suppresses oscillations near shocks, but we find that our current implementation is quite dissipative when solving the Shu-Osher problem. For contact waves we proposed a modification to the entropy viscosity method which eliminates a large amount of the spurious viscosity at contact discontinuities.