# Capturing nonclassical shocks in nonlinear elastodynamic with a conservative finite volume scheme

Nina Aguillon

► To cite this version:

Nina Aguillon. Capturing nonclassical shocks in nonlinear elastodynamic with a conservative finite volume scheme. Interfaces and Free Boundaries: Mathematical Analysis, Computation and Applications, 2016, 18 (2), pp.137-159. 10.4171/IFB/360. hal-01116582

HAL Id: hal-01116582

https://inria.hal.science/hal-01116582

Submitted on 13 Feb 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Capturing nonclassical shocks in nonlinear elastodynamic with a conservative finite volume scheme

Nina Aguillon*

February 13, 2015

## Abstract

For a model of nonlinear elastodynamics, we construct a finite volume scheme which is able to capture nonclassical shocks (also called undercompressive shocks). Those shocks verify an entropy inequality but are not admissible in the sense of Liu. They verify a kinetic relation which describes the jump, and keeps an information on the equilibrium between a vanishing dispersion and a vanishing diffusion. The scheme presented here is by construction exact when the initial data is an isolated nonclassical shock. In general, it does not introduce any diffusion near shocks, and hence nonclassical solutions are correctly approximated. The method is fully conservative and does not use any shock-tracking mesh. This approach is tested and validated on several test cases. In particular, as the nonclassical shocks are not diffused at all, it is possible to obtain large time asymptotics.

**Key words and phrases:** nonclassical shocks, undercompressive shocks, kinetic relation, finite volume schemes

**2010 Mathematics Subject Classification:** 65M08, 76M12, 35L65, 35L67

## Introduction

Hyperbolic systems of conservation laws often arise as limits of systems of partial differential equations including small scales effects, like diffusion and dispersion, when the parameters driving the small scales effects vanish. For example, the compressible Euler equations are the (formal) limit of the compressible Navier-Stokes equations, when the diffusion parameter tends to zero. However, the solutions of the hyperbolic system do not inherit all the properties of the augmented system, and it particular they are not always smooth. More importantly, in the class of weak solutions, the Cauchy problem for the hyperbolic system

<footnote>*nina.aguillon@math.u-psud.fr, Université Aix Marseille. This work has been carried out in the framework of the Labex Archimide (ANR-11-LABX-0033) and of the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the "Investissements d'Avenir" French Government programme managed by the French National Research Agency (ANR). The authors also thanks the Université Paris Sud where this work began.</footnote>
may admit infinitely many solutions, plenty of them being non relevant toward the augmented system. Thus it is necessary to add some criterions that select the solutions that are indeed limits, as the parameters vanish, of the solutions of the augmented diffusive-dispersive system. One of this criterion yields to the selection of classical solutions, which roughly speaking corresponds to the case where diffusion is dominant. But when two small scale effects are of comparable strength (say diffusion and dispersion) nonclassical solutions appear. They contain shocks that do not verify the classical criterion but do arise as limits of the diffusive-dispersive system.

Let us sketch a general framework for nonclassical solutions of hyperbolic systems. Consider the Cauchy problem for a system of conservation laws

$$ \left\{ \begin{array}{l} \partial_t U(t,x) + \partial_x f(U(t,x)) = 0, \\ U(t=0,x) = U^0(x), \end{array} \right. \qquad (1) $$

where $t$ is positive, $x$ belongs to $\mathbb{R}$, $U$ is the vector regrouping the $N \ge 1$ unknowns belonging to an open convex set $\Omega$ of $\mathbb{R}^N$, and $f = \mathbb{R}^N \to \mathbb{R}^N$ is the flux function, which we assume to be smooth. We suppose that System (1) is strictly hyperbolic, i.e. that the Jacobian matrix of $f$ is diagonalizable with simple eigenvalues $\lambda_1(U) < \dots < \lambda_n(U)$. We denote by $r_i(U)$ a right eigenvector associated to the eigenvalue $\lambda_i(U)$.

Solutions of hyperbolic systems are in general not smooth, even when the initial data is smooth. Discontinuities appear in finite time and it is necessary to consider weak solutions. Doing so, uniqueness is lost, and an additional criterion has to be imposed to select a unique solution. For example, the vector of unknown $U$ is asked to verify, in addition to the initial System (1), a so called entropy inequality

$$ \partial_t U(U) + \partial_x W(U) \le 0. \qquad (2) $$

The entropy $U$ and the entropy flux $W$ usually come from considerations on an augmented version of the system under study, where small physical terms like diffusion and dispersion are not neglected. Such an augmented system has additional properties, for example the total energy is conserved. Inequality (2) states that at the limit, this energy is dissipated. In general, from the augmented system only one entropy inequality is deduced, thus the solutions of (1) are not expected to verify other entropy inequalities than (2). When the system is genuinely nonlinear, i.e. when

$$ \forall i \in \{1, \dots, n\}, \quad \forall U \in \Omega, \quad \nabla \lambda_i(U) \cdot r_i(U) \neq 0, \qquad (3) $$

and when the entropy $U$ is strictly convex, the entropy inequality (2) contains enough information to select only one weak solution of (1). In the sequel we suppose that the entropy $U$ is strictly convex. In the general case where (3) does not hold, the addition of a single entropy inequality is not sufficient to select a unique weak solution.

The question of uniqueness of the weak solution is strongly related to the choice of a criterion to decide whether a shock is admissible or not. Two states $U_L$ and $U_R$ of $\Omega$ are linked by a shock if there exists a real $s$, called the speed of the shock, so that the Rankine-Hugoniot relations

$$ s(U_L - U_R) = f(U_L) - f(U_R) \qquad (4) $$
hold. The shock is entropy satisfying if it verifies (2) in a weak sense, i.e. if

$$s(U(U_L) - U(U_R)) - (W(U_L) - W(U_R)) \le 0. \qquad (5)$$

When the system is genuinely nonlinear (i.e. when (3) holds), (5) is equivalent to the Lax inequalities [Lax57]

$$\exists i \in \{1, \dots, n\} : \lambda_i(U_R) < s < \lambda_i(U_L), s < \lambda_{i+1}(U_R) \text{ and } s > \lambda_{i-1}(U_L), \qquad (6)$$

which themselves are an extension to systems ($N > 1$) of Oleinik's criterion [Ole57].

In [Liu74], Liu generalized this criterion to the non genuinely nonlinear case where $\nabla \lambda_i \cdot r_i$ may vanish. Liu proved that his criterion selects a unique selfsimilar solution to the Riemann problem (1-2). We recall that a Riemann problem is the case of a piecewise constant initial data

$$U^0(x) = U_L \mathbf{1}_{x<0} + U_R \mathbf{1}_{x>0}, \quad (U_L, U_R) \in \Omega^2.$$

For systems that are not genuinely nonlinear, the Liu criterion is stronger than (5) (whereas it is equivalent for genuinely nonlinear systems). An admissible shock in the sense Liu always verifies (5), but some shocks may verify (5) without satisfying the Liu criterion. These shocks are called nonclassical. Moreover, unlike in the genuinely nonlinear case, the addition of the single entropy inequality (2) is not enough to regain uniqueness of the weak solution, because too many shocks are allowed. In other words, the inequality entropy (2) does not contain enough information on the small scale effects in the augmented system. Uniqueness can be regain by strongly constraining the nonclassical shocks to verify an algebraic relation of the form

$$U_L = \Phi^b(U_R).$$

The function $\Phi^b$ is called a kinetic relation.

It has been proved that the addition of such a kinetic relation selects a unique weak solution to the Riemann problem for several models, including models of phase transitions in solids in [AK91], [LT01] and [LT02] or in liquids [LT00], [HL00], [SY95] and a model of magnetohydrodynamics [HL00]. The link between kinetic relations and augmented diffusive dispersive equations is explored in [JMS95], [LeF02] and [HL97] in the scalar case and in [Sle89] for gas dynamics with a Van der Waals pressure law. The literature on the subject is so wide that we only cited a few key references. The interested reader will find a more comprehensive bibliography on this topic in [LeF02].

The numerical approximation of nonclassical solutions is challenging, because they are by nature very sensitive to the equilibrium between the small scale effects. The stake is to preserve this equilibrium at the numerical level. Roughly speaking, usual finite volume schemes introduce a numerical viscosity which destroys the equilibrium between the small scales effects. The diffusion becomes dominant and the schemes converge toward the classical solution, even though they are based on nonclassical Riemann solvers.

To overcome this difficulty, two types of schemes have been proposed. On the one hand, it is possible to use high order schemes consistent with the augmented system (see for example [HL98], [LMR02], [LR00] and [KR10]). As outlined in [HL98], it is necessary
to discretize the flux of the hyperbolic part with a high enough order, once again to keep the exact balance between the small scales. On the other hand, some schemes relies on the tracking of nonclassical shocks (see for example [CL03], [CG08], [BCLL08], [Per11], [CCER12] and [CDMG14]). In that case the kinetic relation is taken into account the scheme, typically through the use of an exact nonclassical Riemann solver.

The conservative finite volume scheme presented here belongs to this last category. It is built to be exact when the initial data is a nonclassical shock. In general, it does not introduce any numerical diffusion near nonclassical shocks (a particular class of nonclassical solutions) and turns out to capture correctly nonclassical solutions. This scheme extends to hyperbolic systems the discontinuous reconstruction scheme introduced in the scalar case in [BCLL08] and recently used in [CDMG14].

In the first section, we present the hyperbolic system admitting nonclassical solution for which we construct the scheme, namely the model of nonlinear elasticity studied in [LT01] and [HL00]. The second section is devoted to the construction of the scheme itself. We present a way to select cells in which a nonclassical shock is reconstructed, explain how this shock is reconstructed and give the numerical flux. Eventually, we prove that the scheme is exact when the initial data is a nonclassical shock. In the third and last section, we propose several test cases involving nonclassical shocks. It appears that the proposed scheme capture nonclassical shocks very sharply. Let us outline that the scheme is fully conservative (unlike Glimm type schemes [CL03]) and uses a fixed grid, which allows the computation of solutions containing several interacting nonclassical shocks.

**Acknowledgment** The author warmly thanks Frédéric Lagoutière for his support and advice.

# 1 The Riemann problem for a nonlinear elasticity model

This paper is devoted to the numerical approximation of the solutions of the system of conservation laws

$$ \left\{ \begin{array}{l} \partial_t v - \partial_x \sigma(w) = 0, \\ \partial_t w - \partial_x v = 0, \\ v(t=0, x) = v^0(x), \\ w(t=0, x) = w^0(x), \end{array} \right. \qquad (7) $$

where the stress $\sigma$ is twice differentiable, and verifies

$$ w\sigma''(w) > 0, \sigma'(w) > 0 \quad \text{and} \quad \lim_{|w| \to +\infty} \sigma'(w) = +\infty. \qquad (8) $$

We are interested in weak solutions of (7) that are morally limits when $\epsilon$ tends to $0^+$ of the augmented system

$$ \left\{ \begin{array}{l} \partial_t v^\epsilon - \partial_x \sigma(w^\epsilon) = \epsilon \partial_{xx} v + \alpha \epsilon^2 \partial_{xxx} w, \\ \partial_t w - \partial_x v = 0. \end{array} \right. \qquad (9) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Notations used for solving the Riemann problem. $\Phi^{-\natural}(w)$ is the first Liu shock reachable from $w$. The entropy dissipation between $w$ and $\Phi_{\infty}^{\text{b}}$ is zero. The states in red can be linked to $w$ by a Liu shock, the blue ones by a nonclassical shock.</figcaption></figure>

The parameter $\alpha$ is positive. From (9) we deduce the following entropy inequality for System (7) (see [LeF02]):

$$ \partial_t \left( \frac{v^2}{2} + \int_0^w \sigma(z) dz \right) + \partial_x (-v\sigma(w)) \le 0. \quad (10) $$

In this section we briefly recall the results of Thanh and LeFloch [LT01] on the Riemann problem for System (7), i.e. the case where

$$ \begin{cases} v^0(x) = v_L \mathbf{1}_{x<0} + v_R \mathbf{1}_{x>0}, \\ w^0(x) = w_L \mathbf{1}_{x<0} + w_R \mathbf{1}_{x>0}. \end{cases} \quad (11) $$

We adopt the notation of this paper. It is recalled on Figure 1.

If the left state ($v_L, w_L$) and the right state ($v_R, w_R$) are linked by a 1-shock, its speed is equal to $-s(w_L, w_R)$, where

$$ s(w_L, w_R) = \sqrt{\frac{\sigma(w_L) - \sigma(w_R)}{w_L - w_R}}, \quad (12) $$

and $v_R$ is given by

$$ v_R = v_L + s(w_L, w_R)(w_R - w_L) := \mathcal{H}_1(w_R, v_L, w_L). \quad (13) $$
If they are linked by a 2-shocks, its speed is +s(w_L, w_R), and

$$v_L = v_R - s(w_L, w_R)(w_L - w_R) := \mathcal{H}_2(w_L, v_R, w_R).$$

Graphically, the speed of a shock s(w_L, w_R) corresponds to the slope of the segment joining (w_L, σ(w_L)) to (w_R, σ(w_R)).

Let us now focus on when a 1-shock is admissible or not, either in the sense of the entropy dissipation (10) or in the sense of Liu. The results for 2-shocks are the same, with the roles of w_L and w_R reversed in what follows. A careful study of Equation (10) applied to the shock yields that it verifies the entropy inequality if and only if

$$w_L^2 \le w_L w_R \quad \text{or} \quad w_R w_L \le \Phi_\infty^b(w_L) w_L,$$

where Φ_∞^b(w_L) is the real having the opposite sign than w_L for which the entropy dissipation (5) vanishes. The function Φ_∞^b is its own inverse.

For System (7), a shock verifies the Liu criterion if and only if for all w between w_L and w_R,

$$\frac{\sigma(w) - \sigma(w_R)}{w - w_R} \ge \frac{\sigma(w_L) - \sigma(w_R)}{w_L - w_R}.$$

Geometrically, it means that the segment joining (w_L, σ(w_L)) to (w_R, σ(w_R)) is above the graph of σ if w_R > w_L, and below if w_R < w_L, see Figure 1. It follows that there exists an invertible function Φ^‡ such that a shock verifies the Liu criterion if and only if

$$w_L^2 \le w_L w_R \quad \text{or} \quad w_R w_L \le \Phi^{\ddagger}(w_L) w_L.$$

The fact that System (7) is not genuinely nonlinear yields

$$w_L \Phi_\infty^{\ddagger}(w_L) < w_L \Phi_\infty^b(w_L).$$

Thus, if w_R lies in between Φ_∞^b(w_L) and Φ_∞^‡(w_L), the shock between (v_L, w_L) and (v_R, w_R) dissipates the entropy (10) but is not admissible in the sense of Liu. Such shocks are called *nonclassical shocks*.

Many Riemann problems for System (7-10) admit an infinity of solutions (see [LT01], Section 3). Uniqueness can be obtained by imposing a *kinetic relation* that strongly constrain the nonclassical shocks by imposing

$$\begin{cases} w_L = \Phi^{b,1}(w_R) & \text{for a 1-nonclassical shock,} \\ w_R = \Phi^{b,2}(w_L) & \text{for a 2-nonclassical shock,} \end{cases} \qquad (14)$$

where Φ^{b,1} and Φ^{b,2} both verify

$$w \Phi_\infty^b(w) \le w \Phi^{b,i}(w) \le w \Phi^i(w)$$

but are not necessarily equal. We now state the main result of [LT01].

**Theorem 1.1.** [Thanh, LeFloch] If the kinetic functions Φ^{b,1} and Φ^{b,2} are monotone decreasing, the Riemann Problem (7-10-14) has a unique selfsimilar solution.
The exact Riemann solver associated to Theorem 1.1 is the foundation of the scheme presented in the next section. Moreover, we will use the inner structure of the Riemann solution. Thus for the sake of completeness, we describe below the forward 1-wave and the backward 2-wave.

**Proposition 1.2.** [Thanh, LeFloch] Let $(v_L, w_L)$ be a fixed left state, and $(v_M, w_M)$ be a state linked to $(v_L, w_L)$ by a 1-wave. Then the nature of this 1-wave is the following:

* if $w_L w_M > w_L^2$, it is a classical shock;

* if $0 \le w_L w_M \le w_L^2$, it is a rarefaction wave;

* if $w_L \varphi^{-b,1}(w_L) < w_L w_M \le 0$, it first contains a 1-rarefaction linking $(v_L, w_L)$ to the state $(\mathcal{H}_1(\varphi^{b,1}(w), v_M, w_M), \varphi^{b,1}(w_M))$, which is itself linked to $(v_M, w_M)$ by a 1-nonclassical shock;

* if $w_L w_M < w_L \varphi^{-b,1}(w_L)$, two cases arise:

- if $-s(w_L, \varphi^{b,1}(w_M)) < -s(\varphi^{b,1}(w_M), w_M)$, it first contains a 1-classical shock linking $(v_L, w_L)$ to the state $(\mathcal{H}_1(\varphi^{b,1}(w), v_M, w_M), \varphi^{b,1}(w_M))$, which is itself linked to $(v_M, w_M)$ by a 1-nonclassical shock;

- otherwise, it is just a 1-classical shock (in which $w$ changes sign).

**Proposition 1.3.** [Thanh, LeFloch] Let $(v_R, w_R)$ be a fixed right state, and $(v_M, w_M)$ be a state such that $(v_M, w_M)$ and $(v_R, w_R)$ are linked by a 2-wave. Then the nature of this 2-wave is the following:

* if $w_R w_M > w_R^2$, it is a classical shock;

* if $0 \le w_R w_M \le w_R^2$, it is a rarefaction wave;

* if $w_R \varphi^{-b,2}(w_R) < w_R w_M \le 0$, it first contains a 2-nonclassical shock linking $(v_M, w_M)$ to the state $(\mathcal{H}_2(\varphi^{b,2}(w), v_M, w_M), \varphi^{b,1}(w_M))$, which is itself linked to $(v_R, w_R)$ by a 2-rarefaction wave;

* if $w_R w_M < w_R \varphi^{-b,2}(w_R)$, two cases arise:

- if $s(w_M, \varphi^{b,2}(w_M)) < s(\varphi^{b,2}(w_M), w_R)$, it first contains a 2-nonclassical shock linking $(v_M, w_M)$ to the state $(\mathcal{H}_2(\varphi^{b,2}(w), v_M, w_M), \varphi^{b,1}(w_M))$, which is itself linked to $(v_R, w_R)$ by a 2-classical shock;

- otherwise, it is just a 2-classical shock (in which $w$ changes sign).

The distinction between the last two cases in the last point of the enumerations of Proposition 1.2 and 1.3 can be reexpressed geometrically, as illustrated on Figure 2.

**Proposition 1.4.** The exists a function $\varphi^\sharp$ such that for $i \in \{1,2\}$, for all real $w$ and for all real $w_M$ such that $ww_M < w \varphi^{-b,i}(w)$, the quantity $|s(w, \varphi^{b,i}(w_M))|$ is larger than $|s(\varphi^{b,1}(w_M)), w_M|$ if and only if $w \varphi^{b,i}(w_M)$ is larger than $w \varphi^\sharp(w, w_M)$. In particular if $0 \ge ww_M \ge w_M \varphi^\sharp(w, w_M)$, the solution is a classical shock (in which $w$ changes sign).
<figure><img src="image_3.png" /><figcaption>Figure 2: The solution contains a 2-nonclassical shock (in blue) followed by 2-classical shock (in red) if their speeds are correctly ordered, i.e if the slope of the blue segment is smaller than the red one (right). This depends on the relative position of $\Phi^{b,2}(w_M)$ and $\Phi^{\#}(w_M, w_L)$.</figcaption></figure>

## 2 The scheme

We now describe the scheme. The time discretization is denoted by $t^0 = 0 < t^1 < \dots < t^n < \dots$, and at each time $t^n$, the space is discretized with cells having all the same size $\Delta x$. We denote by $x_{j-1/2}^n$ their extremities and by $x_j^n$ there centers. The centers of the cells will move at each time step but every cell will always be of size $\Delta x$. The speed of the mesh is defined by $V_{\text{mesh}}^n = \frac{x_j^{n+1}-x_j^n}{t^{n+1}-t^n}$. Integrating the first line of (7) over the set

$$ \{(t, x), t^n \le t < t^{n+1}, x_{j-1/2}^n + V_{\text{mesh}}^n t < x < x_{j+1/2}^n + V_{\text{mesh}}^n t\} $$

we obtain the integral formula

$$ \int_{x_{j-1/2}^n}^{x_{j+1/2}^n} v(t^{n+1}, x)dx = \int_{x_{j-1/2}^n}^{x_{j+1/2}^n} v(t^n, x)dx - \left[ \int_{t^n}^{t^{n+1}} (-\sigma(w) - V_{\text{mesh}}^n v)(s, x_{j+1/2}^n + V_{\text{mesh}}^n(s-t^n))ds - \int_{t^n}^{t^{n+1}} (-\sigma(w) - V_{\text{mesh}}^n v)(s, x_{j-1/2}^n + V_{\text{mesh}}^n(s-t^n))ds \right]. $$

The formula is identical for the variable $w$ with the flux $(-v - V_{\text{mesh}}^n w)$ instead of $(-\sigma(w) - V_{\text{mesh}}^n v)$. The principle of finite volume methods is based on this integral formula. A finite
volume scheme writes

$$ \begin{cases} v_j^{n+1} = v_j^n - \frac{t^{n+1}-t^n}{\Delta x} (f_{j+1/2}^{n,v} - f_{j-1/2}^{n,v}), \\ w_j^{n+1} = w_j^n - \frac{t^{n+1}-t^n}{\Delta x} (f_{j+1/2}^{n,w} - f_{j-1/2}^{n,w}), \end{cases} \quad (15) $$

where $v_j^n$ plays the role of the mean value

$$ \frac{1}{\Delta x} \int_{x_{j-1/2}^n}^{x_{j+1/2}^n} v(t^n, x) dx $$

and $f_{j+1/2}^{n,v}$ plays the role of the flux on the edge $x = x_{j-1/2} + V_{\text{mesh}}^n(s - t^n)$

$$ \frac{1}{t^{n+1} - t^n} \int_{t^n}^{t^{n+1}} (-\sigma(w) - V_{\text{mesh}}^n v)(s, x_{j+1/2}^n + V_{\text{mesh}}^n(s - t^n)) ds $$

(and similarly for $w_j^n$ and $f_{j+1/2}^{n,w}$). The scheme is initialized with the exact average of the initial data

$$ \forall j \in \mathbb{Z}, \quad v_j^0 = \frac{1}{\Delta x} \int_{x_{j-1/2}^0}^{x_{j+1/2}^0} v^0(x) dx \quad \text{and} \quad w_j^0 = \frac{1}{\Delta x} \int_{x_{j-1/2}^0}^{x_{j+1/2}^0} w^0(x) dx. \quad (16) $$

A particular finite volume scheme is characterized by a particular choice of a formula expressing the numerical fluxes $f_{j+1/2}^{n,v}$ and $f_{j+1/2}^{n,w}$ in terms of the $(v_k^n)_{k \in \mathbb{Z}}$ and $(w_k^n)_{k \in \mathbb{Z}}$. Our choice is described in the next sections.

Recall that our aim is to derive a scheme which is exact when the initial data is a nonclassical shock. In that case, the initial sampling (16) introduces a small amount of dissipation (unless if by miracle the shock initially falls on an interface). The numerical initial data contains an intermediate value that does not correspond to any pointwise value of $(v^0, w^0)$. More generally, at the end of every time step, finite volume schemes contains a $L^1$-projection on the mesh. Details of the solution are lost in that step, and intermediate values are created in the shocks profile. The main idea of the discontinuous reconstruction scheme is to rebuilt, from the mean values $(v_j^n, w_j^n)_{j \in \mathbb{Z}}$, an initial data that contains more details. This idea is also the foundation of well-known other schemes, like the MUSCL scheme [vL97] and its central version [NT90]. Our reconstruction consists in adding nonclassical shocks in some cells in which we are able to detect that a nonclassical shock where lying inside the cell before the $L^1$-projection. Thus the reconstruction will be very precise near nonclassical shocks, and not on the smooth parts of the solution, taking the opposite strategy of [vL97] and [NT90]

Following this idea, our scheme can be decomposed in three elementary steps.

* First, we detect the special cells which are associated with nonclassical shocks. This detection step is described in Section 2.1.

* Then, we reconstruct nonclassical shocks inside those particular cells, in the sense that we replace the mean values $(v_j^n, w_j^n)$ by piecewise constant functions having the form of a nonclassical shock. This procedure is explained in Section 2.2.
• Eventually, the numerical fluxes are computed by letting the reconstructed nonclassical shocks evolve during the time step $t^{n+1} - t^n$. The use of a moving mesh makes this computation easy, see Section 2.3.

In the last section, we extend the scheme to reconstruct classical shocks as well.

## 2.1 Detection of nonclassical shocks

The key idea of the scheme is to see, whenever it is possible, each mean value $(v_j^n, w_j^n)$ produced by the finite volumes scheme as the average of some nonclassical shock located somewhere inside the cell. Of course, where $(v_{j-1}^n, w_{j-1}^n)$ and $(v_{j+1}^n, w_{j+1}^n)$ are linked by a nonclassical shock, it is the one that should be reconstructed. In general, $(v_{j-1}^n, w_{j-1}^n)$ and $(v_{j+1}^n, w_{j+1}^n)$ are not linked by a single nonclassical shock, but by the full pattern of two waves, each of them likely to contain a succession of a classical shock or a rarefaction followed by a nonclassical shock (see Propositions 1.2 and 1.3). In that case, we reconstruct one of the nonclassical shock appearing in the solution of the Riemann problem. It is chosen thanks to the following Lemma.

**Proposition 2.1.** If the states $(v_L, w_L)$ and $(v_R, w_R)$ are linked by a 1-nonclassical shock, then

$$ w_L w_R < 0, \quad (w_L - w_R)(v_L - v_R) > 0 \quad \text{and} \quad w_L w_R \ge w_R \varphi^{\#}(w_L, w_R). \quad (17) $$

If they are linked by a 2-nonclassical shock, then

$$ w_L w_R < 0 \quad (w_L - w_R)(v_L - v_R) < 0 \quad \text{and} \quad w_L w_R \ge w_L \varphi^{\#}(w_L, w_R). $$

*Proof.* The condition $w_L w_R < 0$ always holds through a nonclassical shock. The second one is a straightforward consequence of the Rankine-Hugoniot relations (4). In the case of nonlinear elastodynamics (7), they write

$$ \begin{cases} (-1)^i s(w_L, w_R)(v_L - v_R) = \sigma(w_R) - \sigma(w_L); \\ (-1)^i s(w_L, w_R)(w_L - w_R) = v_R - v_L. \end{cases} \quad (18) $$

with $i = 1$ for 1-shocks and $i = 2$ for 2-shocks. The positive real $s(w_L, w_R)$ is defined in Formula (12). Eventually, the last condition is a reminder of Proposition 1.4. □

**Definition 2.2.** For all integer $j$ and all positive integer $n$, we denote by $(v_{j,\star}^n, w_{j,\star}^n)$ the intermediate state appearing in the Riemann problem with left state $(v_{j-1}^n, w_{j-1}^n)$ and right state $(v_{j+1}^n, w_{j+1}^n)$. The left and right desired reconstructed states, denoted respectively by $(\bar{v}_{j,L}^n, \bar{w}_{j,L}^n)$ and $(\bar{v}_{j,R}^n, \bar{w}_{j,R}^n)$, are defined as follows.

• If $w_{j-1}^n w_{j+1}^n < 0$, $(w_{j-1}^n - w_{j+1}^n)(v_{j-1}^n - v_{j+1}^n) > 0$, $w_{j-1}^n w_{j+1}^n \ge w_{j+1}^n \varphi^{\#}(w_{j-1}^n, w_{j+1}^n)$ and $w_{j-1}^n w_{j,\star}^n < 0$, the solution contains a 1-shock in which $w$ changes sign.

– If this shock is nonclassical, it is preceded by a classical wave and we set

$$ \left\{ \begin{array}{l} (\bar{v}_{j,L}^n, \bar{w}_{j,L}^n) = (\mathcal{H}_1(\Phi^{\flat,1}(w_{j,\star}^n), v_{j,\star}^n, w_{j,\star}^n), \Phi^{\flat,1}(w_{j,\star}^n)), \\ (\bar{v}_{j,R}^n, \bar{w}_{j,R}^n) = (v_{j,\star}^n, w_{j,\star}^n). \end{array} \right. $$
- If this shock is classical, it is the only contribution in the 1-wave and we set

$$ \begin{cases} (\bar{v}_{j,L}^n, \bar{w}_{j,L}^n) = (v_{j-1}^n, w_{j-1}^n), \\ (\bar{v}_{j,R}^n, \bar{w}_{j,R}^n) = (v_{j,\star}^n, w_{j,\star}^n). \end{cases} $$

* If $w_{j-1}^n w_{j+1}^n < 0$, $(w_{j-1}^n - w_{j+1}^n)(v_{j-1}^n - v_{j+1}^n) < 0$, $w_{j-1}^n w_{j+1}^n \ge w_{j-1}^n \varphi^\#(w_{j-1}^n, w_{j+1}^n)$ and $w_{j+1}^n w_{j,\star}^n < 0$, the solution contains a 2-shock in which w changes sign.

- If this shock is nonclassical, it is followed by a classical wave and we set

$$ \begin{cases} (\bar{v}_{j,L}^n, \bar{w}_{j,L}^n) = (v_{j,\star}^n, w_{j,\star}^n), \\ (\bar{v}_{j,R}^n, \bar{w}_{j,R}^n) = (\mathcal{H}_2(\Phi^{b,2}(w_{j,\star}^n), v_{j,\star}^n, w_{j,\star}^n), \Phi^{b,2}(w_{j,\star}^n)). \end{cases} $$

- If this shock is classical, it is the only contribution in the 2-wave and we set

$$ \begin{cases} (\bar{v}_{j,L}^n, \bar{w}_{j,L}^n) = (v_{j,\star}^n, w_{j,\star}^n), \\ (\bar{v}_{j,R}^n, \bar{w}_{j,R}^n) = (v_{j+1}^n, w_{j+1}^n). \end{cases} $$

* In the other cases, we do not detect any relevant nonclassical shock and we set

$$ \begin{cases} (\bar{v}_{j,L}^n, \bar{w}_{j,L}^n) = (v_j^n, w_j^n), \\ (\bar{v}_{j,R}^n, \bar{w}_{j,R}^n) = (v_j^n, w_j^n). \end{cases} $$

*Remark 2.3.* If there is a nonclassical shock in the Riemann problem with left state $(v_{j-1}^n, w_{j-1}^n)$ and right state $(v_{j+1}^n, w_{j+1}^n)$, then the first two three conditions holds by Proposition 2.1. Note that they are numerically very easy to check. The last condition insures that the Riemann problem indeed contains the expected shock; we verify it only on the cells that passed the first three tests.

## 2.2 Reconstruction

Once a nonclassical shock has been detected, it is placed inside its cell by conservation of the variables $v$ and $w$. We reverse the averaging step of finite volume schemes by replacing the mean value $v_j^n$ by the piecewise constant function

$$ v_{\text{rec}}^n(x) = \bar{v}_{j,L}^n \mathbf{1}_{x<x_{j-1/2}+d_j^{n,v}} + \bar{v}_{j,R}^n \mathbf{1}_{x>x_{j-1/2}+d_j^{n,v}} $$

with

$$ d_j^{n,v} = \Delta x \frac{v_j^n - \bar{v}_{j,R}^n}{\bar{v}_{j,L}^n - \bar{v}_{j,R}^n}. \qquad (19) $$

If $d_j^{n,v}$ belongs to $(0, \Delta x)$ we have

$$ \int_{x_{j-1/2}^n}^{x_{j+1/2}^n} v_{\text{rec}}^n(x) dx = \Delta x v_j^n $$
and no mass is loss when replacing the mean value by $v_{\text{rec}}^n$ inside the j-th cell.
Reasoning similarly for the w variable, we replace $w_j^n$ by

$$w_{\text{rec}}^n(x) = \bar{w}_{j,L}^n \mathbf{1}_{x<x_{j-1/2}+d_j^{n,w}} + \bar{w}_{j,R}^n \mathbf{1}_{x>x_{j-1/2}+d_j^{n,w}},$$

with

$$d_j^{n,w} = \Delta x \frac{w_j^n - \bar{w}_{j,R}^n}{\bar{w}_{j,L}^n - \bar{w}_{j,R}^n}. \qquad (20)$$

We now state one trivial but crucial property of this reconstruction procedure.

**Proposition 2.4.** If $(v_{j-1}^n, w_{j-1}^n)$ and $(v_{j+1}^n, w_{j+1}^n)$ are linked by a single nonclassical shock, and if it exists $\alpha$ in $(0,1)$ such that

$$(v_j^n, w_j^n) = \alpha(v_{j-1}^n, w_{j-1}^n) + (1-\alpha)(v_{j+1}^n, w_{j+1}^n),$$

then

$$(\bar{v}_{j,L}^n, \bar{w}_{j,L}^n) = (v_{j-1}^n, w_{j-1}^n), \quad (\bar{v}_{j,R}^n, \bar{w}_{j,R}^n) = (v_{j+1}^n, w_{j+1}^n) \quad \text{and} \quad d_j^{n,v} = d_j^{n,w} = \alpha \Delta x.$$

In particular if the initial data is a nonclassical shock, we have $v_{\text{rec}}^0 = v^0$ and $w_{\text{rec}}^0 = w^0$: the numerical diffusion introduced by the initial sampling (16) is cancelled by the reconstruction.

In general the two distances $d_j^{n,v}$ and $d_j^{n,w}$ are different, and it is possible that at least one of these distances does not belong to $(0, \Delta x)$. In that case we cancel the reconstruction, considering that seeing $(v_j^n, w_j^n)$ as the mean value of the detected nonclassical shock is not relevant.

**Definition 2.5.** The left and right reconstructed states are defined by:

$$(v_{j,L}^n, w_{j,L}^n) = \begin{cases} (\bar{v}_{j,L}^n, \bar{w}_{j,L}^n) & \text{if } d_j^{n,v} \in (0, \Delta x) \text{ and } d_j^{n,w} \in (0, \Delta x), \\ (v_j^n, w_j^n) & \text{otherwise,} \end{cases} \qquad (21)$$

$$(v_{j,R}^n, w_{j,R}^n) = \begin{cases} (\bar{v}_{j,R}^n, \bar{w}_{j,R}^n) & \text{if } d_j^{n,v} \in (0, \Delta x) \text{ and } d_j^{n,w} \in (0, \Delta x), \\ (v_j^n, w_j^n) & \text{otherwise.} \end{cases} \qquad (22)$$

According to the Rankine-Hugoniot relations (4), the nonclassical shock reconstructed in cell j has velocity

$$s_j^n = \frac{v_{j,R}^n - v_{j,L}^n}{w_{j,L}^n - w_{j,R}^n}$$

and we set arbitrary $s_j^n$ to $\sqrt{\sigma'(w_j)}$ when no reconstruction is performed (which reads $w_{j,L}^n = w_{j,R}^n = w_j^n$).
## 2.3 Advection of the reconstructed discontinuities

The fluxes are computed by letting the reconstructed nonclassical shocks evolve during the time step and by computing exactly what goes through the interfaces. However, two discontinuities reconstructed in adjacent cells can interact and the waves resulting from the interaction can meet the line $x = x_{j+1/2}^n$ within the time step. It follows that if we want to use a fixed grid ($x_j^{n+1} = x_j^n$, or equivalently $V_{\text{mesh}}^n = 0$), the flux along the interface $x = x_{j+1/2}^n$ cannot be computed without resolving the wave interaction, which is obviously extremely costly.

This can been avoiding by using a moving mesh. Let us recall that by moving mesh, we mean that the centers of the cells are moving from time to time, but that their size remains constant equals to $\Delta x$. Thus the numerical difficulties of handling cells with different and varying widths is avoided. The mesh speed is chosen such that it is larger than the maximum of the waves speed:

$$|V_{\text{mesh}}^n| > V_{\text{waves}}^n = \max_{j \in \mathbb{Z}} \sqrt{\sigma'(w_j^n)}, \qquad (23)$$

and the time step such that a wave cannot cross more than an entire cell during the time step:

$$t^{n+1} - t^n \le \frac{\Delta x}{|V_{\text{mesh}}^n| + V_{\text{waves}}^n}. \qquad (24)$$

These two hypothesis insure that any wave created at time $t^n$ inside the j-th cell can only cross the right interface $x = x_{j+1/2}^n + V_{\text{mesh}}^n t$ if $V_{\text{mesh}}^n < 0$ and the left interface $x = x_{j-1/2}^n + V_{\text{mesh}}^n t$ if $V_{\text{mesh}}^n > 0$. It also imply that if two waves interact inside the j-cell during the time step, the waves resulting from their interaction will not have time to catch up the left or the right interfaces. This is depicted on Figure 3. Thus, computing the flux along the j+1/2-th interface $x = x_{j+1/2}^n + V_{\text{mesh}}^n t$ boils down to compute the time at which the nonclassical shock reconstructed in cell j crosses the interface if $V_{\text{mesh}}^n < 0$ (in cell j+1 otherwise). On Figure 3, we also see that if no reconstruction is performed in cells j and j+1, the flux is trivial to compute and is given by a simple left or right decentering (depending on the sign of $V_{\text{mesh}}^n$), which exactly corresponds to the staggered Lax-Friedrichs flux.

We now compute the flux in details, when $V_{\text{mesh}}^n$ is negative and verifies (23), like in Figure 3. Then the nonclassical shock reconstructed in cell j may only cross the right interface $x = x_{j+1/2}^n + V_{\text{mesh}}^n(t - t^n)$. The discontinuity is not located at the same place for the variables v and w, therefore we compute two different crossing times

$$T_{j+1/2}^{w,n} = \frac{\Delta x - d_{j-1}^{n,w}}{s_{j-1}^n - V_{\text{mesh}}^n} \quad \text{and} \quad T_{j+1/2}^{w,n} = \frac{\Delta x - d_{j-1}^{n,w}}{s_{j-1}^n - V_{\text{mesh}}^n} \quad \text{if } V_{\text{mesh}}^n < 0.$$

Once again, Condition (23) insures that the flux passing through the j+1/2 interface only comes from waves arriving from the cell j. It can be clearly seen on Figure 3 that it is
<figure><img src="image_4.png" /><figcaption>Figure 3: Computing the flux through the dashed interfaces can be done without solving waves interactions (interface $j + 1/2$), and is trivial when no reconstruction is performed (cell $j - 1$)</figcaption></figure>

piecewise constant, given by

$$ \left\{ \begin{aligned} f_{j+1/2}^{n,v} &= (-\sigma(w_{j,R}^n) - V_{\text{mesh}}^n v_{j,R}^n) \min(\Delta t^n, T_{j+1/2}^{n,v}) \\ &\quad + (-\sigma(w_{j,L}^n) - V_{\text{mesh}}^n v_{j,L}^n)(\Delta t^n - \min(\Delta t^n, T_{j+1/2}^{n,v})) & \quad \text{if } V_{\text{mesh}}^n < 0 \quad (25) \\ f_{j+1/2}^{n,w} &= (-v_{j,R}^n - V_{\text{mesh}}^n w_{j,R}^n) \min(\Delta t^n, T_{j+1/2}^{n,w}) \\ &\quad + (-v_{j,L}^n - V_{\text{mesh}}^n w_{j,L}^n)(\Delta t^n - \min(\Delta t^n, T_{j+1/2}^{n,w})) \end{aligned} \right. $$

where we denote by $\Delta t^n$ the n-th time step $t^{n+1} - t^n$.

When $V_{\text{mesh}}^n$ is negative, the reasoning is the same but the shock reconstructed in the $j$-th cell now crosses the left interface $x = x_{j-1/2}^n + V_{\text{mesh}}^n(t - t^n)$. The crossing times are given by

$$ T_{j-1/2}^{w,n} = \frac{d_j^{n,w}}{V_{\text{mesh}}^n - s_j^n} \quad \text{and} \quad T_{j-1/2}^{w,n} = \frac{d_j^{n,w}}{V_{\text{mesh}}^n - s_j^n} \quad \text{if } V_{\text{mesh}}^n > 0. $$

and the fluxes are

$$ \left\{ \begin{aligned} f_{j-1/2}^{n,v} &= (-\sigma(w_{j,L}^n) - V_{\text{mesh}}^n v_{j,L}^n) \min(\Delta t^n, T_{j-1/2}^{n,v}) \\ &\quad + (-\sigma(w_{j,R}^n) - V_{\text{mesh}}^n v_{j,R}^n)(\Delta t^n - \min(\Delta t^n, T_{j-1/2}^{n,v})) & \quad \text{if } V_{\text{mesh}}^n > 0 \quad (26) \\ f_{j-1/2}^{n,w} &= (-v_{j,L}^n - V_{\text{mesh}}^n w_{j,L}^n) \min(\Delta t^n, T_{j-1/2}^{n,w}) \\ &\quad + (-v_{j,R}^n - V_{\text{mesh}}^n w_{j,R}^n)(\Delta t^n - \min(\Delta t^n, T_{j-1/2}^{n,w})) \end{aligned} \right. $$

*Remark 2.6.* When no reconstruction is performed, the fluxes (26) and (25) coincide with the staggered Lax-Friedrichs fluxes, i.e., to a simple left or right decentering depending on the sign of $V_{\text{mesh}}^n$. The values assigned to the distances $d_j^{n,v}$ and $d_j^{n,w}$ and to the speed $s_j^n$ (and thus to the crossing times $T_{j+1/2}^{n,v}$ and $T_{j+1/2}^{n,w}$) are of no importance.
Remark 2.7. The idea of using a moving mesh to simplify the computation of the fluxes goes back to the Lax-Friedrichs scheme. The Rusanov scheme follows the same idea by using a local mesh speed $V_{\text{mesh}}^{j,n}$ to make the scheme less diffusive. Higher order extensions based on piecewise polynomial reconstructions are possible, see for example [NT90].

## 2.4 Detection of classical shocks

It is also possible to detect and reconstruct shocks in which $w$ does not change sign. Those shocks are always classical. Their detection is based on the following proposition.

**Proposition 2.8.** If the states $(v_L, w_L)$ and $(v_R, w_R)$ are linked by a 1-shock in which $w$ does not change sign, then either

$$w_R < w_L \le 0 \quad \text{and} \quad v_R < v_L$$

or

$$0 \le w_L < w_R \quad \text{and} \quad v_R > v_L.$$

If they are linked by a 2-shock in which $w$ does not change sign, then either

$$w_L < w_R \le 0 \quad \text{and} \quad v_R < v_L$$

or

$$0 \le w_R < w_L \quad \text{and} \quad v_R > v_L.$$

There is no conflict with the detection of nonclassical shocks of Proposition 2.1, and we can straightforwardly extend Definitions 2.2 and 2.5 to take into account those shocks.

# 3 Exact approximation of isolated nonclassical shocks

The aim of this section is to prove that the scheme described above is exact when the initial data is an isolated nonclassical shock, i.e. (11) with left state $(v_L, w_L)$ and right state $(v_R, w_R)$ verifying the Rankine-Hugoniot relation (18) and constrained by the kinetic relation (14). We recall that in that case the exact solution is

$$\begin{cases} v^{\text{exa}}(t, x) = v_L \mathbf{1}_{x < (-1)^i s(w_L, w_R)t} + v_R \mathbf{1}_{x > (-1)^i s(w_L, w_R)t} \\ w^{\text{exa}}(t, x) = w_L \mathbf{1}_{x < (-1)^i s(w_L, w_R)t} + w_R \mathbf{1}_{x > (-1)^i s(w_L, w_R)t} \end{cases}$$

with $i = 1$ for a 1-shock and $i = 2$ for a 2-shock.

**Theorem 3.1.** The scheme described in the previous section is exact when the initial data is a single nonclassical shock. In other words, the numerical solution is the $L^1$-projection on the mesh of the exact solution at time $t^n$: for all $n \ge 0$ and for all $j \in \mathbb{Z}$,

$$v_j^n = \frac{1}{\Delta x} \int_{x_{j-1/2}^n}^{x_{j+1/2}^n} v^{\text{exa}}(t^n, x) \, dx \quad \text{and} \quad w_j^n = \frac{1}{\Delta x} \int_{x_{j-1/2}^n}^{x_{j+1/2}^n} w^{\text{exa}}(t^n, x) \, dx.$$
<figure><img src="image_5.png" /><figcaption>Figure 4: Structure of w when a 1-phase transition is detected in cell −1. In that case, the Riemann solution (on the right) contains three shocks.</figcaption></figure>

*Proof.* We prove the result for a 1-nonclassical shock such that $w_L < 0 < w_R$, the other cases being exactly similar. With the initial sampling (16), the property holds true at time $t^0$. Suppose that is holds true at some time $t^n$. At this time, we renumbered the cells in a way that the discontinuity lies inside the cell numbered 0, and we denote by $\delta$ its distance to $x_{-1/2}^n$. We have

$$v_j^n = \begin{cases} v_L & \text{if } j < 0, \\ \frac{\delta}{\Delta x} v_L + \frac{\Delta x - \delta}{\Delta x} v_R & \text{if } j = 0, \\ v_R & \text{if } j > 1, \end{cases} \quad \text{and} \quad w_j^n = \begin{cases} w_L & \text{if } j < 0, \\ \frac{\delta}{\Delta x} w_L + \frac{\Delta x - \delta}{\Delta x} w_R & \text{if } j = 0, \\ w_R & \text{if } j > 1. \end{cases}$$

Note that as $w_L < w_R$ and $v_R > v_L$, we have $w_{-1}^n < w_0^n < w_1^n$ and $v_{-1}^n < v_0^n < v_1^n$. The detection step of Proposition 2.1 detects a 1-nonclassical shock in cell 0. By Proposition 2.4, the nonclassical shock is reconstructed in that cell, and it is be placed exactly at the right position in the reconstruction step:

$$(v_{0,L}^n, w_{0,L}^n) = (v_L, w_L), \quad (v_{0,R}^n, w_{j,R}^n) = (v_R, w_R) \quad \text{and} \quad d_j^{0,v} = d_j^{0,w} = \delta.$$

If no reconstruction is performed in the other cells, our scheme gives the correct solution at time $t^{n+1}$ by construction. This is easy to check, but a little tedious, and we refer the reader to [Agu14] for a detailed proof.

Suppose that $w_0^n > 0$. Then by Proposition 2.8, no classical shock is detected in cell 1, and by Proposition 2.1, a 1-nonclassical shock is detected in cell −1. Let us focus on the solution Riemann problem between the state $(v_L, w_L)$ in cell −2 and the state $(v_0^n, w_0^n)$ in cell 0. The notation are recalled on Figure 4. If $\bar{w}_{-1,L}^n \ge w_L$, both $\bar{w}_{-1,L}^n$ and $\bar{w}_{-1,R}^n$ are larger than $w_L$. Thus it is be impossible to have $d_{-1}^{n,\bar{w}} \in (0, \Delta x)$ and no reconstruction occurs in cell −1. Suppose now that $w_L$ is larger than $\bar{w}_{-1,L}^n$. Then $w_L$ and $\bar{w}_{-1,L}^n$ are linked by a 1-classical shock, and $\bar{w}_{-1,L}^n$ is linked to $\bar{w}_{-1,R}^n$ by a 1-nonclassical shock. It is impossible to have $d_{-1}^{n,\bar{w}}$ in $(0, 1)$. Indeed, the kinetic function $\phi^{b,1}$ decreases, thus

$$\bar{w}_{-1,R}^n = \phi^{-b,1}(\bar{w}_{-1,L}^n) > \phi^{-b,1}(w_L) = w_R$$

and there is a classical shock between $(\bar{v}_{-1,R}^n, \bar{w}_{-1,R}^n)$ and $(v_R, w_R)$. The states $(\bar{v}_{-1,R}^n, \bar{w}_{-1,R}^n)$ and $(v_R, w_R)$ are both on the 1-wave curve $W_1^F(v_L, w_L)$. Theorem 4.5 of [LT01] states that
along this curve, $v$ is an increasing function of $w$. Therefore $\bar{v}_{-1,R}^n$ is larger than $v_R$. On the other hand, the Rankine–Hugoniot relations applied to the 2-shock yields that $\bar{v}_{-1,R}^n$ is smaller than $v_{-1}^n$, which contradicts the fact that $v_{-1}^n < v_R$.

The case where a 1-nonclassical shock is detected in cell 1 is simpler. We have $w_0^n < 0$ and nothing is detected on cell $-1$. Moreover, $\bar{w}_{1,L}^n \le \bar{w}_{1,R}^n$, thus it is not possible to reconstruct $w$ in a conservative manner if $\bar{w}_{1,R}^n \le w_R$. On the other hand, if $\bar{w}_{1,R}^n > w_R$, the solution contains a 2-classical shock and (4) yields that $\bar{v}_{1,R}^n \le v_R$. The Rankine–Hugoniot applied to the 1-nonclassical shock implies that $\bar{v}_{1,L}^n \le \bar{v}_{1,R}^n$. Thus it is impossible that $d_1^{n,v}$ belongs to $(0, \Delta x)$. $\square$

*Remark 3.2.* It is crucial to require that both $v$ and $w$ are reconstructed in a conservative manner. In the case where the nonclassical shock is detected on cell 1, it is clear that the proof collapses if we authorize $d_1^{n,v}$ outside $(0, \Delta x)$. Moreover if we authorized $d_1^{n,v}$ outside $(0, \Delta x)$, the solution of the Riemann problem between $(v_1^n, w_1^n)$ and $(v_R, w_R)$ might contain a 1-classical shock, a 1-nonclassical shock and 2-rarefaction wave, in which case

$$ \bar{v}_{1,L}^n \le v_R \le \bar{v}_{1,R}^n $$

and a reconstruction is performed in cell 1.

# 4 Numerical Simulations

For all the numerical simulations presented below, the stress function is $\sigma(w) = w^3 + mw$ with $m \ge 0$. In that case,

$$ \Phi^{\flat}(w) = -\frac{1}{2}w, \quad \Phi^{\flat}_{\infty}(w) = -w \quad \text{and} \quad \Phi^{\flat}(w, w') = -w - w'. $$

For the kinetic functions we take $\Phi^{\flat,1}(w) = \Phi^{\flat,2}(w) = -\beta w$, where $\beta$ belongs to $[-1/2, 1]$. The case $\beta = 1/2$ corresponds to the classical solutions; the choice $\beta = 1$ corresponds to the case where the entropy dissipation (10) is zero across nonclassical choice. It does not fall in the theory of [LeF02], but the Riemann problem can be solved (see [LT01]) and it is possible to explore that case numerically.

## Test 1: Isolated nonclassical shock

This test case illustrates Theorem 3.1. The initial data is the Riemann problem:

$$ \begin{cases} v^0(x) = -10 * \mathbf{1}_{x<0} + 110 * \mathbf{1}_{x \ge 0}, \\ w^0(x) = -6 * \mathbf{1}_{x<0} + 9 * \mathbf{1}_{x \ge 0}, \end{cases} $$

With $m = 1$ and $\beta = 2/3$, the solution is an isolated nonclassical 1-shock. On the top of Figure 5, we plot the exact nonclassical solution and the solution given by the reconstruction scheme. As expected they are exactly the same. On the bottom of Figure 5, we plot the solution given by the Godunov scheme based on an exact *nonclassical* Riemann solver. It
does not capture the nonclassical solution but the classical one, which in that case in a rarefaction followed by a shock. This is a general phenomenon: usual finite volume schemes are not able to capture nonclassical solutions. Thus in the sequel we used the Glimm scheme to compare our scheme with. The CFL number is set to 0.45, the final time is $T = 0.038$ and the space interval $[-0.5, 0.5]$ is discretized with 200 cells.

<figure><img src="image_6.png" /><figcaption>Figure 5: The nonclassical Godunov scheme is unable to capture nonclassical shocks. On the contrary, the reconstruction scheme captures isolated nonclassical shocks exactly.</figcaption></figure>

## Interlude: The Glimm scheme

In the sequel we use the random sampling method of Glimm [Gli65] to compute reference solution. The Glimm scheme is built as follow. Let $(r^n)_{n \in \mathbb{N}}$ be a sequence of i.i.d. random variables uniformly distributed over $[0, \Delta x]$. Suppose that at time $t^n$ a piecewise constant approximation of the solution is given, and denote by $(t, x) \mapsto U_{\text{exa}}^n(t, x)$ the exact solution at time $t \ge 0$ for that initial data. The numerical solution of Glimm at time $t^{n+1} = t^n + \Delta t^n$ is given by

$$ \forall j \in \mathbb{Z}, \quad U_j^{n+1} = U_{\text{exa}}^n(\Delta t^n, x_{j-1/2}^n + r^n). $$
It $\Delta t^n$ is small enough, i.e. if it is smaller than the maximum of the waves speed appearing in the Riemann problems multiplied by the cell size $\Delta x$, $U_{\text{exa}}^n$ is just the juxtaposition of Riemann solutions at each interface.

The main feature of the Glimm scheme is that it does not introduce any numerical diffusion, in the sense that the shock profiles are not smeared out at all. This is why this scheme has been used in [CL03] to approximate nonclassical solutions. The two drawbacks are that the scheme is not conservative and that the exact computation of the solution is costly.

## Test 2: A Riemann problem with two nonclassical shocks

The initial data is now

$$ \begin{cases} v^0(x) = 6 * \mathbf{1}_{x<0} - 10 * \mathbf{1}_{x \ge 0}, \\ w^0(x) = 1 * \mathbf{1}_{x<0} + 2 * \mathbf{1}_{x \ge 0}. \end{cases} $$

The exact solution consists in a 1-classical shock, a 1-nonclassical shock, a 2-nonclassical shock and a 2-rarefaction wave. On Figure 6, we plot the exact solution and the numerical solution given by the reconstruction schemes. We compare the reconstruction where we detect nonclassical shocks only (cf Proposition 2.1), referred to in the legends as RecNC, and the reconstruction scheme were the classical shocks are also detected (cf Proposition 2.8), referred to in the legends as RecNC+C. Both of them capture very sharply the nonclassical shocks; the 1-classical shock is much more diffused when it is not reconstructed, and both scheme behaves in the same way in the rarefaction wave. The CFL number is 0.45, the space interval $[-1, 1]$ contains 200 cells.

## Test 3: Perturbation of a classical shock

This test case is taken from [CL03]. We still have $\beta = 2/3$ but now $m = 2$. The initial data is

$$ \begin{cases} v^0(x) = 1 * \mathbf{1}_{x<0} - 11 * \mathbf{1}_{x \ge 0}, \\ w^0(x) = (1 + \epsilon) * \mathbf{1}_{x<0} - 3 * \mathbf{1}_{x \ge 0}. \end{cases} $$

with $\epsilon \in \{0, 0.05, 0.1\}$ When $\epsilon = 0$, the solution is a 1-classical shock. When $\epsilon > 0$, the classical shock is split into a classical shock followed by a nonclassical shock. Their speeds are different but close to each other. We plot the solutions at time $T = 0.4$ on Figure 7, using 600 cells per unit interval and a CFL number of 0.45. The solutions are very well approached when $\epsilon > 0$. When $\epsilon = 0$, a spike appears. This spike corresponds to the state linked to $(v_R, w_R)$ by a nonclassical shock (on Figure 7 we see that it has the exact same height as the nonclassical shock appearing for $\epsilon > 0$). Remark that $\epsilon = 0$ is exactly the limit between the two cases in the last point of the enumeration of Proposition 1.2, for which the speed of the nonclassical and the classical shock coincide (see also Figure 2).

Numerically the mechanism is the following. After one iteration in time, an intermediate value $(v_i, w_i) = (\alpha v_L + (1 - \alpha)v_R, \alpha w_L + (1 - \alpha)w_R)$ is created. At the second iteration in time, the classical shock (in which $w$ changes sign), is perfectly reconstructed in the corresponding cell. However, a 1-nonclassical shock is also detected by Proposition 2.1 in
<figure><img src="image_7.png" /><figcaption>Figure 6: A Riemann problem with two nonclassical shocks at time 0.15.</figcaption></figure>

the cell just before, and this time the reconstruction succeeds. Thus the scheme is not exact in that case. Note that the proof of Theorem 3.1 uses the kinetic relation and thus collapses when the reconstructed shock does not verify the kinetic relation (and is hence classical).

This phenomena does not prevent the scheme from converging. The numerical solution has the same shape than in the case $\epsilon > 0$: a classical shock followed by a nonclassical shock, but this time they have the same speed, thus they remain at the same position. The spike is only two cells wide when the classical shocks are reconstructed, and hence not diffused. If they are not reconstructed, the spurious classical shock is diffused.

In conclusion, this test case shows a limitation of scheme, namely its incapability to approach exactly classical shock in which $w$ changes sign, but also show its ability to capture finely nonclassical solutions, the gap between the shocks being very thin here. In particular, the apparition of the intermediate state is immediate, while it is linked to the ratio of the width of the gap at time $T$ and of the cell size $\Delta x$ when using a random sampling based method like the Glimm scheme.

## Test 4: Apparition of nonclassical waves from a smooth initial data

We now focus on the case of a smooth initial data

$$ \begin{cases} v^0(x) = 3 \sin(2\pi x), \\ w^0(x) = 1 + 3 \cos(8\pi x), \end{cases} $$
<figure><img src="image_8.png" /><figcaption>Figure 7: Perturbations of an isolated classical shock, with, from top to bottom, $\epsilon = 0$, $\epsilon = 0.05$ and $\epsilon = 0.1$.</figcaption></figure>
with $\beta = 0.95$ and $m = 1$ and periodic boundary condition. Nonclassical shocks appear around time $t = 0.011$ and then propagate in the solution. On Figures 8, 9 and 10 we compare the solution given by the reconstruction schemes (with or without reconstruction of classical shocks) at time $t = 0.015$, $t = 0.06$ and $t = 1$. The space interval contains 1000 cells per unit interval and CFL number is set to 0.45. The reference solution is given by the Glimm scheme with 8000 cells. We see that the reconstruction schemes capture accurately the nonclassical shocks. The result are poorer in smooth areas, because in those regions the reconstruction schemes behave as the Lax-Friedrichs scheme which is quite diffusive. This can be improved by using another scheme on interfaces where no reconstruction is performed. The use of a moving mesh encourages to chose a central scheme like the Nessyahu and Tadmor scheme [NT90]. This scheme is second order accurate in smooth regions of the solutions. It is both easy to implement and fast, as it does not use any information on the Riemann problems.

<figure><img src="image_9.png" /><figcaption>Figure 8: Solution of test 4 at time $t = 0.015$.</figcaption></figure>

## Test 5: Long time simulation in the maximal dissipative case

We reproduce here the test case of [CL03] for which $\beta = 1$ and the final time is large. The case $\beta = 1$ corresponds to the case where the entropy dissipation is zero through nonclassical shocks. It is a limit case in the theoretical framework of [LeF02]. As the Riemann solver is perfectly defined for $\beta = 1$, it can be explored numerically.
<figure><img src="image_10.png" /><figcaption>Figure 9: Solution of test 4 at time $t = 0.06$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Solution of test 4 at time $t = 0.1$.</figcaption></figure>
The initial data is 1-periodic with

$$ (v^0(x), w^0(x)) = \begin{cases} (0.3, 0.4) & \text{if } 0 \le x < 0.3, \\ (0.15, -0.2) & \text{if } 0.3 \le x < 0.3 + 2/3, \\ (0.1, 0.4) & \text{if } 0.3 + 2/3 \le x < 1. \end{cases} $$

Its mean value is null. The parameter $m$ is fixed to 0.05. On Figure 11, we plot the solution at time $t = 40$ with 2000 and 8000 points per unit interval. We can see that at that time $w$ changes sign three times, instead of two in the initial solution, and does not converge to zero. The position of the nonclassical shocks are different with the two schemes. Let

<figure><img src="image_12.png" /><figcaption>Figure 11: Solution of test 5 at time $t = 40$ (the shape of $v$ is similar)</figcaption></figure>

us recall that the Glimm scheme is based on a random sampling of the solution, thus two realizations of the same test give different results. On Figure 12, we plot the histogram of the first nonclassical shock position at time $t = 20$ for 100 independent realizations of the Glimm scheme with 2000 cells (bottom) and the comparison of the two schemes at the same time (with 8000 points). Moreover, the structure with 2 nonclassical shocks very close to each other only appears in 31 out of those 100 realizations. It indicates that the width of this structure is small (and probably smaller than the one appearing in the realization of [CL03]).

## References

[Agu14] Nina Aguillon. A nondissipative reconstruction scheme for the compressible euler equations. Preprint, hal-00967484, 2014.
<figure><img src="image_13.png" /><figcaption>Figure 12: Solution of test 5 at time $t = 40$ (the shape of $v$ is similar)</figcaption></figure>

[AK91] Rohan Abeyaratne and James K. Knowles. Kinetic relations and the propagation of phase boundaries in solids. *Arch. Rational Mech. Anal.*, 114(2):119–154, 1991.

[BCLL08] Benjamin Boutin, Christophe Chalons, Frédéric Lagoutière, and Philippe G. LeFloch. Convergent and conservative schemes for nonclassical solutions based on kinetic relations. *I. Interfaces Free Bound.*, 10(3):399–421, 2008.

[CCER12] Ch. Chalons, F. Coquel, P. Engel, and Ch. Rohde. Fast relaxation solvers for hyperbolic-elliptic phase transition problems. *SIAM J. Sci. Comput.*, 34(3):A1753–A1776, 2012.

[CDMG14] Christophe Chalons, Maria Laura Delle Monache, and Paola Goatin. A conservative scheme for non-classical solutions to a strongly coupled pde-ode problem. *Preprint*, 2014.

[CG08] Christophe Chalons and Paola Goatin. Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling. *Interfaces Free Bound.*, 10(2):197–221, 2008.

[CL03] C. Chalons and P. G. LeFloch. Computing undercompressive waves with the random choice scheme. Nonclassical shock waves. *Interfaces Free Bound.*, 5(2):129–158, 2003.

[Gli65] James Glimm. Solutions in the large for nonlinear hyperbolic systems of equations. *Comm. Pure Appl. Math.*, 18:697–715, 1965.
[HL97] Brian T. Hayes and Philippe G. LeFloch. Non-classical shocks and kinetic relations: scalar conservation laws. *Arch. Rational Mech. Anal.*, 139(1):1–56, 1997.

[HL98] Brian T. Hayes and Philippe G. Lefloch. Nonclassical shocks and kinetic relations: finite difference schemes. *SIAM J. Numer. Anal.*, 35(6):2169–2194 (electronic), 1998.

[HL00] Brian T. Hayes and Philippe G. Lefloch. Nonclassical shocks and kinetic relations: strictly hyperbolic systems. *SIAM J. Math. Anal.*, 31(5), 2000.

[JMS95] Doug Jacobs, Bill McKinney, and Michael Shearer. Travelling wave solutions of the modified Korteweg-de Vries-Burgers equation. *J. Differential Equations*, 116(2):448–467, 1995.

[KR10] Frederike Kissling and Christian Rohde. The computation of nonclassical shock waves with a heterogeneous multiscale method. *Netw. Heterog. Media*, 5(3):661–674, 2010.

[Lax57] Peter D. Lax. Hyperbolic systems of conservation laws. II. *Commun. Pure Appl. Math.*, 10:537–566, 1957.

[LeF02] Philippe G. LeFloch. *Hyperbolic systems of conservation laws*. Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, 2002. The theory of classical and nonclassical shock waves.

[Liu74] Tai-Ping Liu. The Riemann problem for general 2 × 2 conservation laws. *Trans. Am. Math. Soc.*, 199:89–112, 1974.

[LMR02] P. G. Lefloch, J. M. Mercier, and C. Rohde. Fully discrete, entropy conservative schemes of arbitrary order. *SIAM J. Numer. Anal.*, 40(5):1968–1992 (electronic), 2002.

[LR00] Philippe G. Lefloch and Christian Rohde. High-order schemes, entropy inequalities, and nonclassical shocks. *SIAM J. Numer. Anal.*, 37(6):2023–2060 (electronic), 2000.

[LT00] Philippe G. LeFloch and Mai Duc Thanh. Nonclassical Riemann solvers and kinetic relations. III. A nonconvex hyperbolic model for van der Waals fluids. *Electron. J. Differential Equations*, pages No. 72, 19 pp. (electronic), 2000.

[LT01] Philippe G. LeFloch and Mai Duc Thanh. Nonclassical Riemann solvers and kinetic relations. I. A nonconvex hyperbolic model of phase transitions. *Z. angew. Math. Phys.*, 52, 2001.

[LT02] Philippe G. LeFloch and Mai Duc Thanh. Non-classical Riemann solvers and kinetic relations. II. An hyperbolic-elliptic model of phase-transition dynamics. *Proc. Roy. Soc. Edinburgh Sect. A*, 132(1):181–219, 2002.
[NT90] Haim Nessyahu and Eitan Tadmor. Nonoscillatory central differencing for hyperbolic conservation laws. *J. Comput. Phys.*, 87(2):408–463, 1990.

[Ole57] O. A. Olešnik. On the uniqueness of the generalized solution of the Cauchy problem for a non-linear system of equations occurring in mechanics. *Uspehi Mat. Nauk (N.S.)*, 12(6(78)):169–176, 1957.

[Per11] Vincent Perrier. A conservative method for the simulation of the isothermal Euler system with the van-der-Waals equation of state. *J. Sci. Comput.*, 48(1-3):296–303, 2011.

[Sle89] M. Slemrod. A limiting “viscosity” approach to the Riemann problem for materials exhibiting change of phase. *Arch. Rational Mech. Anal.*, 105(4):327–365, 1989.

[SY95] Michael Shearer and Yadong Yang. The Riemann problem for a system of conservation laws of mixed type with a cubic nonlinearity. *Proc. Roy. Soc. Edinburgh Sect. A*, 125(4):675–699, 1995.

[vL97] Bram van Leer. Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method. *J. Comput. Phys.*, 135(2):227–248, 1997.