A Finite-Volume discretization of viscoelastic Saint-Venant equations for FENE-P fluids

TL;DR

This paper extends a finite-volume discretization method for viscoelastic shallow flow equations to FENE-P fluids, demonstrating practical numerical stability and applicability despite some unresolved theoretical stability guarantees.

Contribution

It introduces an extension of the finite-volume scheme to FENE-P fluids, broadening the applicability of the method to more complex viscoelastic models.

Findings

Numerical simulations showed smooth solutions in practical parameter ranges.

The method preserves physical free-energy dissipation properties.

Stability conditions are not fully guaranteed theoretically but are observed numerically.

Abstract

Saint-Venant equations can be generalized to account for a viscoelastic rheology in shallow flows. A Finite-Volume discretization for the 1D Saint-Venant system generalized to Upper-Convected Maxwell (UCM) fluids was proposed in [Bouchut \& Boyaval, 2013], which preserved a physically-natural stability property (i.e. free-energy dissipation) of the full system. It invoked a relaxation scheme of Suliciu type for the numerical computation of approximate solution to Riemann problems. Here, the approach is extended to the 1D Saint-Venant system generalized to the finitely-extensible nonlinear elastic fluids of Peterlin (FENE-P). We are currently not able to ensure all stability conditions a priori, but numerical simulations went smoothly in a practically useful range of parameters.

Figures

Equations

$${\partial }_{t}h+{\partial }_x(hu)$$

$${\partial }_t(hu)+{\partial }_x\left(h{u}^2+g{h}^2/2+hN\right)$$

$$\lambda \left({\partial }_t{\sigma }_{xx}+u{\partial }_x{\sigma }_{xx}+2(\zeta -1){\sigma }_{xx}{\partial }_{x}u\right)$$

$$\lambda \left({\partial }_t{\sigma }_{zz}+u{\partial }_x{\sigma }_{zz}+2(1-\zeta ){\sigma }_{zz}{\partial }_{x}u\right)$$

$$F=h\left(\frac{u^2}{2}+\frac{gh}{2}-\frac{G}{2(1-\zeta )}\left(\ell \log \left((\ell -({\sigma }_{xx}+{\sigma }_{zz}))/(\ell -2)\right)+\log ({\sigma }_{xx}{\sigma }_{zz})\right)\right)$$

$$-\frac{Gh}{2(1-\zeta)\lambda}\Bigg{(}\sigma_{xx}^{-1}\left(1-\frac{\sigma_{xx}}{1-(\sigma_{zz}+\sigma_{xx})/\ell}\right)^{2}+\sigma_{zz}^{-1}\left(1-\frac{\sigma_{zz}}{1-(\sigma_{zz}+\sigma_{xx})/\ell}\right)^{2}\Bigg{)}\\ =:D\geq\partial_{t}F+\partial_{x}\left(u(F+P)\right)$$

$${\mathcal{U}}^{\ell }:=\{0<h,0<{\sigma }_{xx},0<{\sigma }_{zz},{\sigma }_{xx}+{\sigma }_{zz}<\ell \}.$$

$${\partial }_{t}h+{\partial }_x(hu)$$

$${\partial }_t(hu)+{\partial }_x\left(h{u}^2+\frac{g{h}^2}{2}+hN\right)$$

$${\partial }_t(h\Pi )+{\partial }_x(hu\Pi )$$

$${\partial }_t(h\Sigma )+{\partial }_x(hu\Sigma )$$

$${\partial }_{t}q+A(q){\partial }_{x}q=S(q),$$

$$q_i^{n+1}\approx \frac{1}{\Delta x_i}{\int }_{x_{i-1/2}}^{x_{i+1/2}}q(t,x)dx,i\in \mathbb{Z},t\in (t^n,t^{n+1}]$$

$$q_i^{n+1/2}=q_i^n-\frac{\Delta t^n}{\Delta x_i}\left(F_l(q_i^n,q_{i+1}^n)-F_r(q_{i-1}^n,q_i^n)\right),$$

$$q_i^{n+1}=q_i^{n+1/2}+\Delta t^{n}S(q_i^{n+1}).$$

$$q_i^{n+1}=q_i^n-\frac{\Delta t^n}{\Delta x_i}\left(F_l(q_i^n,q_{i+1}^n)-F_r(q_{i-1}^n,q_i^n)\right)+\Delta t^{n}S(q_i^{n+1})$$

$$F(q_i^{n+1/2})-F(q_i^n)+\frac{\Delta t^n}{\Delta x_i}\left(G(q_i^n,q_{i+1}^n)-G(q_{i-1}^n,q_i^n))\right)\le 0$$

$$F(q_i^{n+1})-F(q_i^n)+\frac{\Delta t^n}{\Delta x_i}\left(G(q_i^n,q_{i+1}^n)-G(q_{i-1}^n,q_i^n))\right)\le \Delta t^{n}D(q_i^{n+1})$$

$$F(q_i^{n+1})-F(q_i^{n+1/2})\le \Delta t^{n}D(q_i^{n+1})\le 0.$$

$$\lambda \left({\sigma }_{xx,i}^{n+1}-{\sigma }_{xx,i}^n\right)/\Delta t^n$$

$$\lambda \left({\sigma }_{zz,i}^{n+1}-{\sigma }_{zz,i}^n\right)/\Delta t^n$$

$${\nabla }_{({\sigma }_{xx},{\sigma }_{zz})}F{|}_{h,u}\cdot S=D$$

$$\begin{array}[]{l}F_{l}(q_{l},q_{r})=F_{0}(q_{l})-\int_{-\infty}^{0}\Bigl{(}R(\xi,q_{l},q_{r})-q_{l}\Bigr{)}d\xi,\\ F_{r}(q_{l},q_{r})=F_{0}(q_{r})+\int_{0}^{\infty}\Bigl{(}R(\xi,q_{l},q_{r})-q_{r}\Bigr{)}d\xi.\end{array}$$

$$R(\xi ,q_l,q_r)=L\mathcal{R}\left(\xi ,{\mathcal{Q}}_l,{\mathcal{Q}}_r\right)$$

$$\{\begin{array}{rl}{\partial }_{t}h+{\partial }_x(hu)& =0\\ {\partial }_t(hu)+{\partial }_x(h{u}^2+\pi )& =0\\ {\partial }_t({\sigma }_{xx}{h}^{2(1-\zeta )})+u{\partial }_x({\sigma }_{xx}{h}^{2(1-\zeta )})& =0\\ {\partial }_t({\sigma }_{zz}{h}^{2(\zeta -1)})+u{\partial }_x({\sigma }_{zz}{h}^{2(\zeta -1)})& =0\\ {\partial }_t(h\pi )+{\partial }_x(hu\pi +u{c}^2)& =0\\ {\partial }_t\left(h(u^2/2+\hat{e})\right)+{\partial }_x\left(hu(u^2/2+\hat{e})+u\pi \right)& =0\\ {\partial }_{t}c+u{\partial }_{x}c& =0\end{array}$$

$${\mathcal{Q}}_o=\left(h_o,(hu{)}_o,h_o^{1-2\zeta }(h{\sigma }_{xx}{)}_o,h_o^{2\zeta -3}(h{\sigma }_{zz}{)}_o,h_{o}P(q_o),(hu{)}_o^2/2h_o+e(q_o),c_o\right)$$

$$\Delta t^n\le \frac{1}{2}\underset{i\in \mathbb{Z}}{\inf }\frac{1}{\Delta x_i}\max \left(u_i^n-c_l(q_i^n,q_{i+1}^n)/h_i^n,u_i^n+c_r(q_i^n,q_{i+1}^n)/h_{i+1}^n\right).$$

$$q_{\xi }:=L\mathcal{R}\left(\xi ,{\mathcal{Q}}_l,{\mathcal{Q}}_r\right)\in {\mathcal{U}}^{\ell }\text{ and }{h}_{\xi }^2{\partial }_h{|}_{h^{2-2\zeta }{\sigma }_{xx},h^{2\zeta -2}{\sigma }_{zz}}P(q_{\xi })\le c_{\xi }^2,\forall \xi \in \mathbb{R}$$

$$\frac{1}{h_l^{\ast }}$$

$$\frac{1}{h_r^{\ast }}$$

$$(h_o^{\ast }{)}^{2(1-\zeta )}(h_o{)}^{2(\zeta -1)}{\sigma }_{zz,o}+(h_o^{\ast }{)}^{2(\zeta -1)}(h_o{)}^{2(1-\zeta )}{\sigma }_{xx,o}<\ell$$

$$w_o{A}_o+w_o^{-1}{B}_o<1\iff 2A_o{w}_o\in \left(1-\sqrt{1-4A_o{B}_o},1+\sqrt{1-4A_o{B}_o}\right)\subset {\mathbb{R}}_{>0}$$

$$(w_o^{+}{)}^{\frac{1}{2(1-\zeta )}}:={\left((1+\sqrt{1-4A_o{B}_o})/(2A_o)\right)}^{\frac{1}{2(1-\zeta )}}\ge \text{αo−1αo}\ge 1/y_o\ge h_o^{\ast }/h_o.$$

$$(w_o^{-}{)}^{\frac{1}{2(1-\zeta )}}\le {\left(1+a_o{h}_o{W}_o/c_o\right)}^{-1}\le h_o^{\ast }/h_o.$$

Full text

MSC (2010): 65M08, 65N08, 35Q30

1 Introduction

Saint-Venant equations standardly model shallow free-surface gravity flows and can be generalized to account for the viscoelastic rheology of non-Newtonian fluids bouchut-boyaval-2015 , Upper-Convected Maxwell (UCM) fluids in particular bouchut-boyaval-2013 . Here, we consider a generalized Saint-Venant (gSV) system for finitely-extensible nonlinear elastic fluids with Peterlin closure (FENE-P fluids) in Cartesian coordinates

[TABLE]

for 1D $\vec{e}_{y}$-translation invariant flow along $\vec{e}_{x}$ under a uniform gravity field $-g\vec{e}_{z}$ with

• •

mean flow depth $h(t,x)>0$ (in case of a non-rugous flat bottom),

• •

mean flow velocity $u(t,x)$ (for uniform cross sections), and

• •

a normal-stress difference $N=G(\sigma_{zz}-\sigma_{xx})/(1-(\sigma_{zz}+\sigma_{xx})/\ell)$ given by conformation variables $\sigma_{zz},\sigma_{xx}>0$ constrained by $0<\sigma_{zz}+\sigma_{xx}<\ell$, a relaxation time $\lambda\geq 0$ and an elasticity modulus $G>0$.

Note that (1-2-3-4) formally reduces to the standard viscous Saint-Venant system with viscosity $\nu\equiv 2\lambda G\geq 0$ when $\ell\rightarrow\infty$, $\lambda\rightarrow 0$ and $G\lambda<\infty$. Moreover we have used the quite general Gordon-Schowalter derivatives with slip parameter $\zeta\in[0,\frac{1}{2}]$ constrained by the hyperbolicity of the system (1-2-3-4). (This follows after an easy computation similar to 2016arXiv161108491B .)

In this work, we discuss a Finite-Volume method to solve (numerically) the Cauchy problem for the nonlinear hyperbolic 1D system (1-2-3-4). Standardly, we need to consider weak solutions (in fact, to (6-7-8-9), see below) plus admissibility constraints that are physically-meaningful dissipation rules formalizing the thermodynamics second principle close to an equilibrium dafermos-2000 . Here, we consider the inequality associated with the companion conservation law for the free-energy

[TABLE]

that is, on denoting the impulse by $P=gh^{2}/2+hN$,

[TABLE]

where the left-hand-side is obviously non-positive on the admissibility domain

[TABLE]

Note that we do not consider the vacuum state $h=0$ as admissible here, see 2016arXiv161108491B .

2 Finite-Volume discretization of FENE-P/Saint-Venant

Piecewise-constant approximate solutions to the Cauchy problem on $(t,x)\in[0,T)\times\mathbb{R}$ for the gSV system can be defined by a Finite-Volume (FV) method. With a view to preserving $\mathcal{U}^{\ell}$ and the dissipation (5) after discretization by a FV method, we choose $q=(h,hu,h\sigma_{xx},h\sigma_{zz})$ as discretization variable. Indeed, the free-energy functional $F$ is convex on the convex domain $\mathcal{U}^{\ell}\ni q$ (this follows after an easy computation from (bouchut-2004, , Lemma 1.3)) while it is not convex in the variable $\left(h,hu,h\Pi,h\Sigma\right)$ whatever smooth invertible functions $\varpi,\varsigma$ are used for the reformulation of gSV

[TABLE]

with $\Pi=\varpi(\sigma_{xx}h^{2(1-\zeta)})$, $\Sigma=\varsigma(\sigma_{zz}h^{2(\zeta-1)})$ (computations are similar to (bouchut-boyaval-2013, , Appendix)). In the sequel, we therefore discretize a quasilinear system with source

[TABLE]

which we recall is not ambiguous here (for those discontinuous solutions built using a Riemann solver, at least) thanks to the dissipation rule (5), see lefloch-2002 ; berthon-coquel-lefloch-2011 ; 2016arXiv161108491B .

2.1 Splitting-in-time

In cell $(x_{i-1/2},x_{i+1/2})$, $i\in\mathbb{Z}$, with volume $\Delta x_{i}=x_{i+1/2}-x_{i-1/2}>0$ and center $x_{i}={(x_{i-1/2}+x_{i+1/2})}/2$, we approximate $q$ solution to (10) on $\mathbb{R}_{\geq 0}\times\mathbb{R}\ni(t,x)$ by

[TABLE]

on a time grid $0=t^{0}<t^{1}<\ldots<t^{n}<t^{n+1}<\ldots<t^{N}=T$ where $\Delta t^{n}=|t^{n+1}-t^{n}|$ will be chosen small enough compared with $\Delta x=\sup_{i\in\mathbb{Z}}\Delta x_{i}<\infty$ to ensure stability.

More precisly, having in mind the numerical approximation of a (well-posed) Cauchy problem for (10) on $\mathbb{R}_{\geq 0}\times\mathbb{R}$ with initial condition $q(t\rightarrow 0^{+})=q^{0}\in L^{\infty}(\mathbb{R})$, and therefore starting from approximations $q_{i}^{0}\approx\frac{1}{\Delta x_{i}}\int_{x_{i-1/2}}^{x_{i+1/2}}q^{0}(x)dx$, $i\in\mathbb{Z}$, we shall define the cell values $q_{i}^{n}$ in two steps for each $n=1,\ldots,N$:

(i) an approximate solution to the homogeneous gSV system (i.e. without the source term $S$) on $[t^{n},t^{n+1})$ is first computed by an explicit three-point scheme

[TABLE]

(ii) an approximate solution to the full gSV system on $(t^{n},t^{n+1}]$ is next computed as

[TABLE]

Then, the scheme is consistent with weak solutions of (1–2) equiv. (6–7)

[TABLE]

provided the two first flux components for the conservative part $(h,hu)$ of the variable $q$ (actually solutions to conservation laws) are conservative $F_{l,h}=F_{r,h}:=F_{h}$, $F_{l,hu}=F_{r,hu}:=F_{hu}$ and consistent $F_{h}(q,q)=hu|_{q}$, $F_{hu}(q,q)=(hu^{2}+gh^{2}/2+hN)|_{q}$ as usual, and with the conservative interpretation (8–9) of (3–4) insofar as we next define $F_{l}$ and $F_{r}$ using a simple approximate Riemann solver harten-lax-vanleer-1983 for (6–7–8–9).

Moreover, with a view to preserving $\mathcal{U}^{\ell}$ and a discrete version of (5)

[TABLE]

for a numerical free-energy flux function consistent with $G(q,q)=u(F+P)|_{q}$ in (5), in the sequel, we shall discuss the relaxation technique introduced by Suliciu as simple Riemann solver in step (i), because it proved satisfying for other close systems bouchut-2003 ; bouchut-2004 ; bouchut-boyaval-2013 equipped with an “entropy” convex in the discretization variable like $F$ here. In the end, for the full scheme (13), a consistent free-energy dissipation

[TABLE]

then holds insofar $h_{i}^{n+1/2}=h_{i}^{n+1}$, $u_{i}^{n+1/2}=u_{i}^{n+1}$ and the convexity of $F$ imply

[TABLE]

Proof

On noting $h_{i}^{n+1/2}=h_{i}^{n+1}$, $u_{i}^{n+1/2}=u_{i}^{n+1}$ it suffices to show that

[TABLE]

imply (16). Now, this is obvious, on noting the convexity of $F|_{h,u}$ in $(\sigma_{xx},\sigma_{zz})$ and

[TABLE]

since $\boldsymbol{\nabla}_{(\sigma_{xx}h^{2(1-\zeta)},\sigma_{zz}h^{2(\zeta-1)})}F\cdot(h^{2(\zeta-1)}S_{h\sigma_{xx}},h^{2(1-\zeta)}S_{h\sigma_{zz}})=D$ by design.

2.2 Suliciu relaxation of the Riemann problem without source

For all time ranges $t\in[t^{n},t^{n+1})$, $n=0\ldots N-1$, let us now define at each interface $x_{i+\frac{1}{2}}$, $i\in\mathbb{Z}$, between cells $i$ and $i+1$ the numerical flux functions $F_{l}$ and $F_{r}$

[TABLE]

invoking an approximate solution $R\left((x-x_{i+1/2})/(t-t^{n}),q_{i}^{n},q_{i+1}^{n}\right)$ to the Riemann problem for (10) with initial condition $q_{i}^{n}1_{x<0}+1_{x>0}q_{i+1}^{n}$ at $t=t^{n}$, and any $F_{0}$.

In this work, we propose as approximate solution that given by Suliciu relaxation

[TABLE]

i.e. the projection (operator $L$) onto $q\equiv(h,hu,h\sigma_{xx},h\sigma_{zz})$ of the exact solution $\mathcal{R}\left(\xi,\mathcal{Q}_{l},\mathcal{Q}_{r}\right)$ of the Riemann problem for the system with relaxed pressure

[TABLE]

and initial condition given by ($o=l,r$)

[TABLE]

where $c_{o}(q_{l},q_{r})$ are chosen so as to ensure stability, that is the dissipation rule (14) here (see below). Note that (19) is a hyperbolic system which fully decomposes into linearly degenerate eigenfields, so $\mathcal{R}$ has an analytic expression (see formulas in bouchut-2004 ; bouchut-boyaval-2013 ). Note also: the Riemann solver $R$ is consistent under the CFL condition

[TABLE]

It remains to specify a choice of functions $c_{l},c_{r}$ preserving $\mathcal{U}^{\ell}$ and ensuring (14).

Although it is not clear whether our construction allows one to approximate solutions on any time ranges $t\in[0,T)$, since the series $\sum_{n}\Delta t^{n}$ may be bounded uniformly for all space-grid choice ($\sup_{i}|u_{i}^{n}|$ may grow unboundedly as $n\rightarrow\infty$), specifying such $c_{l},c_{r}$ fully defines a computable scheme. In particular, (15) then implies that (12) at step (ii) always has at least one solution $q_{i}^{n+1}\in\mathcal{U}^{\ell}$ for any $\Delta t^{n}$ fixed at step (i). (This can be shown using Brouwer fixed-point theorem like in barrett-boyaval-2011 .)

Note however a difficulty here for FENE-P fluids with $c_{l},c_{r}$. Suliciu relaxation approach (19) was retained at step (i) because the solver often allows one to preserve invariant domains like $\mathcal{U}^{\ell}$ and a dissipation rule (14) through well-chosen $c_{l},c_{r}$, see e.g. bouchut-2003 ; bouchut-2004 ; bouchut-boyaval-2013 . Indeed, on noting the exact Riemann solution to (19), to get (14) on choosing G(q_{l},q_{r})=u\Bigl{(}h\bigl{(}\frac{u^{2}}{2}+\hat{e}\bigr{)}+\pi\Bigr{)}|_{\mathcal{R}(0,q_{l},q_{r})}, it is enough that $\forall q_{l},q_{r}\in\mathcal{U}^{\ell}$

[TABLE]

using $c_{\xi}=c_{l}(q_{l},q_{r})$ if $\xi<=u^{*}$ and $c_{\xi}=c_{r}(q_{l},q_{r})$ if $\xi>u^{*}$ with $u^{*}:=\frac{c_{l}u_{l}+\pi_{l}+c_{r}u_{r}-\pi_{r}}{c_{l}+c_{r}}$.

One can easily propose $c_{l},c_{r}$ satisfying the first condition in (22), i.e.

[TABLE]

as usual for Saint-Venant systems, plus the admissibility conditions ($o=l/r$)

[TABLE]

for any $\sigma_{zz,o},\sigma_{xx,o}>0$ satisfying $\sigma_{zz,o}+\sigma_{xx,o}<\ell$ (FENE-P fluids, see below). But the second condition is usually treated for $\phi_{o}:h\rightarrow h\sqrt{\partial_{h}|_{h^{2-2\zeta}_{o}\sigma_{xx,o},h^{2\zeta-2}_{o}\sigma_{zz,o}}P}$ monotone. Unfortunately, a lengthy (but easy) computation shows that the latter is not monotone here, so the standard method to choose $c_{l},c_{r}$ a priori does not apply.

2.3 Choice of relaxation parameter

Let us treat the first part of (22) as usual and define $c_{o}=\max(h_{o}\sqrt{\partial_{h}P(q_{o})}:=h_{o}a_{o},\tilde{c}_{o})$, $o=l/r$ such that the functions $\tilde{c}_{o}(q_{l},q_{r})$ ensure (23–24) and (25).

First, let us inspect (23–24) classically following (bouchut-klingenberg-waagan-2010, , section3.3). Denoting $a_{l}Y_{l}=(u_{l}-u_{r})_{+}+\frac{(\pi_{r}-\pi_{l})_{+}}{h_{l}a_{l}+h_{r}a_{r}}\geq 0$, $a_{r}Y_{r}=(u_{l}-u_{r})_{+}+\frac{(\pi_{l}-\pi_{r})_{+}}{h_{l}a_{l}+h_{r}a_{r}}\geq 0$ so $\frac{1}{h_{o}^{*}}\geq\frac{1-{h_{o}a_{o}Y_{o}}/{c_{o}}}{h_{o}}$, it then holds $(h_{o}^{*})^{-1}\geq(h_{o})^{-1}y_{o}>0$ with $y_{o}:=1-\frac{Y_{o}}{1+\alpha_{o}Y_{o}}\in(\frac{\alpha_{o}-1}{\alpha_{o}},1]$ provided one chooses $\tilde{c}_{o}>0$ such that $c_{o}\geq h_{o}a_{o}(1+\alpha_{o}Y_{o})$ for $\alpha_{o}>1$, which yields $h_{o}^{*}\in(0,h_{o}/y_{o}]$ thus (23–24) in particular.

On the other hand, let us now inspect (25), which rewrites with $h_{o}^{*}>0$

[TABLE]

with $w_{o}=(h_{o}^{*}/h_{o})^{2(1-\zeta)}$, $A_{o}=\sigma_{zz,o}/\ell$, $B_{o}=\sigma_{xx,o}/\ell$ positive such that $A_{o}+B_{o}<1$ (hence $A_{o}B_{o}\leq A_{o}(1-A_{o})\leq\frac{1}{4}$) and $2(1-\zeta)\in[1,2]$. The upper-bound in (26) is satisfied with $\alpha_{o}=(w_{o}^{+})^{\frac{1}{2(1-\zeta)}}/((w_{o}^{+})^{\frac{1}{2(1-\zeta)}}-1)>1$, on noting

[TABLE]

It remains to ensure the lower bound in (26). Obviously, $w_{o}^{-}:=\frac{1-\sqrt{1-4A_{o}B_{o}}}{2A_{o}}<1$ so one only needs to inspect the case $h_{o}^{*}\leq h_{o}$. Now, with $a_{l}W_{l}=(u_{r}-u_{l})_{+}+\frac{(\pi_{l}-\pi_{r})_{+}}{h_{l}a_{l}+h_{r}a_{r}}\geq 0$, $a_{r}W_{r}=(u_{r}-u_{l})_{+}+\frac{(\pi_{r}-\pi_{l})_{+}}{h_{l}a_{l}+h_{r}a_{r}}\geq 0$, if $c_{o}\geq h_{o}a_{o}W_{o}({(w_{o}^{-})^{-\frac{1}{2(1-\zeta)}}}-1)^{-1}$ then holds

[TABLE]

In the end, we claim the following choices

[TABLE]

satisfy simultaneously (23–24) and (25) in a compatible way with $a_{o}=\sqrt{\partial_{h}P(q_{o})}$, $\alpha_{o}=\max(2,(w_{o}^{+})^{\frac{1}{2(1-\zeta)}}/((w_{o}^{+})^{\frac{1}{2(1-\zeta)}}-1))$, $\beta_{o}=(w_{o}^{-})^{\frac{1}{2(1-\zeta)}}/(1-(w_{o}^{-})^{\frac{1}{2(1-\zeta)}}))$, $w_{o}^{-}=\frac{\ell-\sqrt{\ell-4\sigma_{zz,o}\sigma_{xx,o}}}{2\sigma_{zz,o}}$, $w_{o}^{+}=\frac{\ell+\sqrt{\ell-4\sigma_{zz,o}\sigma_{xx,o}}}{2\sigma_{zz,o}}$, for $o=l/r$. Moreover, note that we have chosen $\alpha_{o}$ such that all subcharacteristic conditions (22) are satisfied in the $\ell\rightarrow\infty$ limit, hence also the free-energy dissipation (15). Indeed, $\phi_{o}$ is monotone in the $\ell\rightarrow\infty$ limit and one can then apply the standard method to choose $c_{l},c_{r}$ bouchut-boyaval-2013 .

3 Numerical illustrations

We numerically approximate on $t\in[0,.1]$ the solution to a Riemann problem with

[TABLE]

as initial condition when $g=10$, $\zeta=0$, $G=.1$, $\lambda=.1$. In Fig. 1, we show the initial condition and the result at $t=.1$ when $\Delta x=2^{-8}$ for $\ell=10,100,1000$. Note the influence of the parameter $\ell$ on the stretch $\sigma_{xx}+\sigma_{zz}$. On computing numerically the free-energy dissipation with the choice of relaxation parameter above, we have never observed the wrong sign, while the time-step did not go to zero.