A space–time smooth artificial viscosity method for nonlinear conservation laws

Abstract

We introduce a new methodology for adding localized, space–time smooth, artificial viscosity to nonlinear systems of conservation laws which propagate shock waves, rarefactions, and contact discontinuities, which we call the C-method. We shall focus our attention on the compressible Euler equations in one space dimension. The novel feature of our approach involves the coupling of a linear scalar reaction–diffusion equation to our system of conservation laws, whose solution $C(x\text{,}t)$ is the coefficient to an additional (and artificial) term added to the flux, which determines the location, localization, and strength of the artificial viscosity. Near shock discontinuities, $C(x\text{,}t)$ is large and localized, and transitions smoothly in space–time to zero away from discontinuities. Our approach is a provably convergent, spacetime-regularized variant of the original idea of Richtmeyer and Von Neumann, and is provided at the level of the PDE, thus allowing a host of numerical discretization schemes to be employed.

We demonstrate the effectiveness of the C-method with three different numerical implementations and apply these to a collection of classical problems: the Sod shock-tube, the Osher–Shu shock-tube, the Woodward–Colella blast wave and the Leblanc shock-tube. First, we use a classical continuous finite-element implementation using second-order discretization in both space and time, FEM-C. Second, we use a simplified WENO scheme within our C-method framework, WENO-C. Third, we use WENO with the Lax–Friedrichs flux together with the C-equation, and call this WENO-LF-C. All three schemes yield higher-order discretization strategies, which provide sharp shock resolution with minimal overshoot and noise, and compare well with higher-order WENO schemes that employ approximate Riemann solvers, outperforming them for the difficult Leblanc shock tube experiment.

Introduction

The initial-value problem for a general nonlinear system of conservation laws can be written as an evolution equation,

$${\partial }_{t}U(x\text{,}t)+\mathrm{div}F(U(x\text{,}t))=0\mathrm{with}U{|}_{t=0}=U_0\text{,}$$

for an m-vector U defined on (D + 1)-dimensional space–time. Such partial differential equations (PDE) are both ubiquitous and fundamental in science and engineering, and include the compressible Euler equations of gas dynamics, the magneto-hydrodynamic (MHD) equations modeling ionized plasma, the elasticity equations of solid mechanics, and numerous related physical systems which possess complicated nonlinear wave interactions.

It is well known that solutions of (1) can develop finite-time shocks, even when the initial data is smooth, in which case, discontinuities of U are propagated according to the so-called Rankine-Hugoniot conditions (see Section 2.1). It is important to develop stable and robust numerical algorithms which can approximate shock-wave solutions. Even in one-space dimension, nonlinear wave interaction such as two shock waves colliding, is a difficult problem when considering accuracy, stability and monotonicity. The challenge is maintaining higher-order accuracy away from the shock while approximating the discontinuity in an order-$\mathrm{\Delta }x$ smooth transition region where $\mathrm{\Delta }x$ denotes the spatial grid size.

As we describe below, a variety of clever discretization schemes have been developed and employed, particularly in one-space dimension, to approximate discontinuous solution profiles in an essentially non-oscillatory (ENO) fashion. These include, but are not limited to, total variation diminishing (TVD) schemes, flux-corrected transport (FCT) schemes, weighted essentially non-oscillatory (WENO) schemes, discontinuous Galerkin methods, artificial diffusion methods, exact and approximate Riemann solvers, and a host of variants and combinations of these techniques.

We develop a robust parabolic-type regularization of (1), which we refer to as the C-method, which couples a modified set of m equations for U with an additional linear scalar reaction–diffusion equation for a new scalar field $C(x\text{,}t)$. Thus, instead of (1) we consider a system of m + 1 equations, which use the solution $C(x\text{,}t)$ as a coefficient in a carefully chosen modification of the flux. As we describe in detail below, the solution $C(x\text{,}t)$ is highly localized in regions of discontinuity, and transitions smoothly (in both x and t) to zero in regions wherein the solution is smooth. Further, as $\mathrm{\Delta }x\to 0$, we recover the original hyperbolic nonlinear system of conservation laws (1).

In the case of 1-D gas dynamics, the construction of non-oscillatory, higher-order, numerical algorithms such as ENO by Harten et al. [1] and Shu and Osher [2], [3]; WENO by Liu et al. [4] and Jiang and Shu [5]; MUSCL by Van Leer [6], Colella [7], and Huynh [8]; or PPM by Colella and Woodward [9] requires carefully chosen reconstruction and numerical flux.

Such numerical methods evolve cell-averaged quantities; to calculate an accurate approximation of the flux at cell-interfaces, these schemes reconstruct kth-order ($k⩾2$) polynomial approximations of the solution (and hence the flux) from the computed cell-averages, and thus provide kth-order accuracy away from discontinuities. See, for example, the convergence plots of Greenough and Rider [10] and Liska and Wendroff [11]. Given a polynomial representation of the solution, a strategy is chosen to compute the most accurate cell-interface flux, and this is achieved by a variety of algorithms. Centered numerical fluxes, such as Lax–Friedrichs, add dissipation as a mechanism to preserve stability and monotonicity. On the other hand, characteristic-type upwinding based upon exact (Godunov) or approximate (Roe, Osher, HLL, HLLC) Riemann solvers, which preserve monotonicity without adding too much dissipation, tend to be rather complex and PDE-specific; moreover, for strong shocks, other techniques may be required to dampen post-shock oscillations or to yield entropy-satisfying approximations (see Quirk [12]). Again, we refer the reader to the papers [10], [11] or Colella and Woodward [13] for a thorough overview, as well as a comparison of the effectiveness of a variety of competitive schemes.

Majda and Osher [14] have shown that any numerical scheme is at best, first-order accurate in the presence of shocks or discontinuities. The use of higher-order numerical schemes is, nevertheless, imperative for the elimination of error-terms in the Taylor expansion (in mesh-size) and the subsequent limiting of truncation error. Moreover, higher-order schemes tend to be less dissipative than there lower-order counterparts, as discussed by Greenough and Rider [10]; therein, a comparison between a 2nd-order PLMDE scheme and a 5th-order WENO scheme demonstrates the improved resolution of intricate fine structure afforded by 5th-order WENO, while simultaneously providing far less clipping of local extrema than PLMDE.

In multi-D, similar tools are required to obtain non-oscillatory numerical schemes, but the multi-dimensional analogues to those described above are generally limited by mesh considerations. For structured grids (such as products of uniform 1-D grids), dimensional splitting is commonly used, decomposing the problem into a sequence of 1-D problems. This technique is quite successful, but stringent mesh requirements prohibits its use on complex domains. Moreover, applications to PDE outside of variants of the Euler equations may be somewhat limited. For further discussion of the limitations of dimensional splitting, we refer the reader to Crandall and Majda [15], and Jiang and Tadmor [16]. For unstructured grids, dimensional splitting is not available and alternative approaches must be employed, necessitated by the lack of multi-D Riemann solvers. WENO schemes on unstructured triangular grids have been developed in Hu and Shu [17], but using simplified methods, which employ reduced characteristic decompositions, can lead to a loss of monotonicity and stability.

Algorithms that explicitly introduce diffusion provide a simple way to stabilize higher-order numerical schemes and subsequently remove non-physical oscillations near shocks. In the mathematical analysis of conservation laws (and in the truncation error of certain discretization schemes), the simplest parabolic-regularization is by the addition of a uniform linear viscosity. Choosing a constant $\beta >0$, which depends upon mesh-size $\mathrm{\Delta }x$ and sometimes velocity or wave-speed, and adding

$$\beta (\mathrm{\Delta }x){\partial }_x^{2}U(x\text{,}t)$$

to the right hand side of (1) provides a uniformly parabolic regularization of the hyperbolic conservation laws, and its discrete implementation smears sharp discontinuities across $O(\mathrm{\Delta }x)$-regions and thus adds stabilization, but unfortunately, at the cost of accuracy. With the addition of uniform linear viscosity, shocks and discontinuities are captured in a non-oscillatory fashion, but the transition region from left to right state, which approximates the discontinuity, tends to grow over time. Moreover, since viscosity is applied uniformly over the entire domain $\mathcal{I}$, the benefits of a higher-order scheme (away from the discontinuity) may be lost, and the accuracy may reduce to merely first-order (at best). In numerical schemes, the use of viscosity should be localized in regions of shock (so as to stabilize the scheme), limited at contact discontinuities (to avoid over-smearing the sharp transition), and very small in smooth regions away from discontinuities. Achieving these requirements allows higher-order approximation of smooth flow and sharp, non-oscillatory, resolution of shocks and discontinuities. Naturally, this necessitates that the amount of added viscosity be a function of the solution.

The pioneering papers of Richtmyer [18], Von Neumann and Richtmyer [19], Lax and Wendroff [20], and Lapidus [21] suggest the introduction of nonlinear artificial viscosity to Eqs. (1) in a form similar to the following expression:

$$\beta (\mathrm{\Delta }x{)}^2{\partial }_x(|{\partial }_{x}u(x\text{,}t)|{\partial }_{x}U(x\text{,}t))\text{.}$$

We refer the reader to the classical papers of Gentry et al. [22] and Harlow and Amsden [23] for an interesting discussion on artificial viscosity. Specifically, Gentry et al. [22] define the nonlinear viscosity of the type (3) to be artificial viscosity, and show that the linear viscosity (2), scaled by the magnitude of local velocity, arises as truncation error (in finite-difference approximations). The latter is responsible for stabilizing the transport of sound waves, while (3) stabilizes the steepening of sound waves.¹

We are primarily concerned with the steepening of sound waves, and shall term artificial viscosity of the type (3) as classical artificial viscosity. Formally, the use of (3) produces the required amount of viscosity near shocks but allows for second-order accuracy in smooth regions. On the other hand, the diffusion coefficient $|{\partial }_{x}u(x\text{,}t)|$ is precisely the quantity which loses regularity (or smoothness) near shock discontinuities. Also, the constant β must be larger than one to control numerical oscillations behind the shock wave, which in turn overly diffuses the waves and produces incorrect wave speeds.

Alternative procedures have been proposed. For streamline upwind Petrov–Galerkin schemes (SUPG), Hughes and Mallet [24] and Shakib et al. [25] use residual-based artificial viscosity. Guermond and Pasquetti [26] present a similar, entropy-residual-based scheme for use in spectral methods. Persson and Peraire [27] develop a method based upon decay of local interpolating polynomials for discontinuous Galerkin schemes. Later, Barter and Darmofal [28] use a reaction–diffusion equation to provide a regularized variant of this approach.

Our approach is similar to [28] in that it uses a reaction–diffusion equation to calculate a smooth distribution of artificial viscosity. Instead of regularizing a DG-based noise-indicator that allows for the growth of viscosity near shocks, we regularize the classical artificial viscosity of [21], using a gradient based approach for this source term. This approach yields both a discretization- independent and PDE-independent methodology which can be generalized to multiple dimensions by regularizing a similar viscosity to that in Löhner et al. [29].

In 1-D, our approach proves to be a simple way of circumventing the need for characteristic or other a priori information of the exact solution to remove oscillations in higher-order schemes. Due to the simple and discretization-independent nature of our method, we expect our methodology to be useful for a wide range of applications.

In Section 2, we introduce the C-method for the compressible Euler equations in one space dimension. We show that the C-method is Galilean invariant and that solutions of the C-method converge to the entropy solutions of the Euler equations in the limit of zero mesh size. We also show the relative smoothness of our new viscosity coefficient with respect to the classical artificial viscosity of Richtmyer and Von Neumann, and we demonstrate the ability of the C-method to remove downstream oscillation in slowly moving shocks.

In Section 3, we give a brief outline of the numerical schemes whose solutions are used in this paper. First, we outline a second-order, continuous Galerkin finite-element method. Second, we outline a simple WENO-based finite-volume scheme which performs upwinding using only the sign of the velocity (no Riemann-solvers or characteristic decompositions in primitive variables). The resulting schemes applied to the C-method are referred to as FEM-C and WENO-C, respectively. Third, we outline the central-finite-difference scheme of Nessyahu and Tadmor (NT), a simple scheme, easily generalizable to multi-D [30]. Like our FEM-C scheme, the NT-scheme is at best, second-order, and does not require specialized techniques for upwinding. Fourth, we outline a Godunov-type characteristic decomposition-based WENO scheme (WENO-G) developed by Rider et al. [31] which utilizes a variant of a Godunov/Riemann-solver as upwinding, providing a very competitive scheme for modeling the collision of very strong shocks.

In Section 4, we consider the classical shock-tube problem of Sod. With the Sod shock problem, we apply our FEM-C scheme and compare with the classical viscosity approach. We then compare the FEM-C scheme with the two standalone methods, NT and WENO-G.

In Section 5, we consider the moderately difficult problem of Osher–Shu, modeling the interaction of a mild shock with an entropy wave. We compare FEM-C to NT and WENO-G in which the differences are more significant than in the Sod-shock comparisons. We show that WENO-C compares well with WENO-G; on the other hand, the simple WENO scheme without the C-method and without the Gudonov-based characteristic solver also does well in modeling the Osher–Shu test case.

In Section 6, we consider the numerically challenging Woodward–Colella blast wave simulation, which models the collision of two strong interacting shock fronts. Though the FEM-C scheme performs better than NT, both second-order schemes are somewhat out-performed by the higher-order WENO-G method (with characteristic solver). On the other hand, WENO-C compares well with WENO-G, having slightly less damped amplitudes with the same shock resolution.

Finally, in Section 7, we consider the Leblanc shock-tube, an extremely difficult test case consisting of a very strong shock. For this problem, we devise two strategies to demonstrate the use of the C-method. In the first strategy, we use our simplified WENO-C scheme with a right-hand side term for the energy equation that relies on a second C-equation which smooths gradients of $E/\rho$. We obtain an excellent approximation of the notoriously difficult contact discontinuity for internal energy, while maintaining an accurate shock speed; simultaneously, we avoid generating large overshoots at the contact discontinuity, which would indeed occur without the use of the C-method. For our second strategy, we show that WENO with the Lax–Friedrichs flux can be significantly improved with the addition of the C-method. We call this algorithm WENO-LF-C, and show that by using just one C-equation (as we have for all of the other test cases), we can sharply resolve the contact discontinuity for the internal energy, with accurate wave speed, and without overshoots.