Viscoelastic flows with conservation laws

TL;DR

This paper introduces a novel symmetric hyperbolic system of conservation laws for viscoelastic Maxwell fluids, extending classical models to include fluids with memory and providing a framework for theoretical analysis and numerical simulation.

Contribution

It presents the first symmetric hyperbolic conservation law system for Maxwell fluid flows, generalizing existing models and enabling new computational approaches.

Findings

Proposes a new hyperbolic PDE system for Maxwell fluids.

Extends the Saint-Venant system to viscoelastic flows.

Includes a finite-volume discretization for numerical simulations.

Abstract

We propose in this work the first symmetric hyperbolic system of conservation laws to describe viscoelastic flows of Maxwell fluids, i.e. fluidswith memory that are characterized by one relaxation-time parameter. Precisely, the system of quasilinear PDEs is detailed for the shallow-water regime, i.e. for hydrostatic incompressible 2D flows with free surface under gravity. It generalizes Saint-Venant system to viscoelastic flows of Maxwell fluids, and encompasses previous works with F. Bouchut. It also generalizes the (thin-layer) elastodynamics of hyperelastic materials to viscous fluids, and to various rheologies between solid and liquid states that can be formulated using our new variable as material parameter.The new viscoelastic flow model has many potential applications, additionally to falling into the theoretical framework of (symmetric hyper-bolic) systems of conservation laws. In computational rheology, it offers a new approach to the High-Weissenberg Number Problem (HWNP). Fortransient geophysical flows, it offers perspectives of thermodynamically-compatible numerical simulations, with a Finite-Volume (FV) discretization say. Besides, one FV discretization of the new continuum model is proposed herein to precise our ideas incl. the physical meaning of the solutions. Perspectives are finally listed after some numerical simulations.

Equations

$${\partial }_{t}H+\mathrm{div}(H\mathbf{U})=0,$$

$${\partial }_t(H\mathbf{U})+\mathrm{div}(H\mathbf{U}\otimes \mathbf{U}+H(P+{\Sigma }_{zz})\mathbf{I}-H{\mathbf{\Sigma }}_h)=0,$$

$$U^x\approx \frac{1}{H}{\int }_0^{H}dz{u}^x{U}^y\approx \frac{1}{H}{\int }_0^{H}dz{u}^y$$

$${\mathbf{\Sigma }}_h={\Sigma }_{xx}{\mathbf{e}}_x\otimes {\mathbf{e}}_x+{\Sigma }_{yy}{\mathbf{e}}_y\otimes {\mathbf{e}}_y+{\Sigma }_{xy}{\mathbf{e}}_x\otimes {\mathbf{e}}_y+{\Sigma }_{yx}{\mathbf{e}}_y\otimes {\mathbf{e}}_x$$

$$\begin{gathered} {\Sigma }_{xx}=2\nu {\partial }_x{U}^x{\Sigma }_{yy}=2\nu {\partial }_y{U}^y{\Sigma }_{xy}=\nu ({\partial }_x{U}^y+{\partial }_y{U}^x)={\Sigma }_{yx} \\ {\Sigma }_{zz}=-\nu ({\partial }_x{U}^x+{\partial }_y{U}^y) \end{gathered}$$

$$H{\int }_0^{H}dz{u}^i{u}^j\approx \left({\int }_0^{H}dz{u}^i\right)\left({\int }_0^{H}dz{u}^j\right)$$

$${\partial }_{t}H+{\partial }_x(H{U}^x)+{\partial }_y(H{U}^y)=0$$

$${\partial }_t(H{U}^x)+{\partial }_x(H{U}^x{U}^x+HP+H{\Sigma }_{zz}-H{\Sigma }_{xx})+{\partial }_y(H{U}^x{U}^y-H{\Sigma }_{xy})=-KH{U}^x$$

$${\partial }_t(H{U}^y)+{\partial }_x(H{U}^x{U}^y-H{\Sigma }_{yx})+{\partial }_y(H{U}^y{U}^y+HP+H{\Sigma }_{zz}-H{\Sigma }_{yy})=-KH{U}^y$$

$${\partial }_t(HE)+\mathrm{div}(HE\mathbf{U}+H(P+{\Sigma }_{zz})\mathbf{U}-H{\mathbf{\Sigma }}_h\cdot \mathbf{U})=-KH|\mathbf{U}{|}^2-HD$$

$${\partial }_{t}e+(\mathbf{U}\cdot \nabla )e+(P+{\Sigma }_{zz})\mathrm{div}\mathbf{U}-{\mathbf{\Sigma }}_h:\nabla \mathbf{U}=-D.$$

$$E=\frac{1}{2}\left(|\mathbf{U}{|}^2+gH\right),$$

$${\partial }_t(HE)+\mathrm{div}(HE\mathbf{U}+H(P+{\Sigma }_{zz})\mathbf{U}-H{\mathbf{\Sigma }}_h\cdot \mathbf{U})\le -KH|\mathbf{U}{|}^2$$

$$D_t{\mathbf{\Sigma }}_h-{\mathbf{L}}_h{\mathbf{\Sigma }}_h-{\mathbf{\Sigma }}_h{\mathbf{L}}_h^T=(2\nu {\mathbf{D}}_h-{\mathbf{\Sigma }}_h)/\lambda$$

$$D_t{\Sigma }_{zz}+2\mathrm{div}\mathbf{U}{\Sigma }_{zz}=(-2\nu \mathrm{div}\mathbf{U}-{\Sigma }_{zz})/\lambda$$

$$\boxed{\begin{array}{rl}& {\partial }_{t}H+\mathrm{div}(H\mathbf{U})=0\\ & {\partial }_t(H\mathbf{U})+\mathrm{div}(H\mathbf{U}\otimes \mathbf{U}+(g{H}^2/2+GH{B}_{zz})\mathbf{I}-GH{\mathbf{B}}_h)=-KH\mathbf{U}\\ & {\partial }_t{\mathbf{B}}_h+\mathbf{U}\cdot \nabla {\mathbf{B}}_h-{\mathbf{L}}_h{\mathbf{B}}_h-{\mathbf{B}}_h{\mathbf{L}}_h^T=(\mathbf{I}-{\mathbf{B}}_h)/\lambda \\ & {\partial }_t{B}_{zz}+\mathbf{U}\cdot \nabla B_{zz}+2B_{zz}\mathrm{div}\mathbf{U}=(1-B_{zz})/\lambda \end{array}}$$

$$H{\mathbf{\Sigma }}_h=|{\mathbf{F}}_h{|}^{-1}{\partial }_{{\mathbf{F}}_h}\left(\frac{G}{2}{\mathbf{F}}_h:{\mathbf{F}}_h\right){\mathbf{F}}_h^T\text{ and }H{\mathbf{\Sigma }}_{zz}=|F_c^z|{\partial }_{F_c^z}\left(\frac{G}{2}|F_c^z{|}^2\right)F_c^z$$

$${\partial }_t\mathbf{F}+(\mathbf{u}\cdot \nabla )\mathbf{F}=\mathbf{L}\mathbf{F}.$$

$${\partial }_t{\mathbf{F}}_h+(\mathbf{u}\cdot \nabla ){\mathbf{F}}_h={\mathbf{L}}_h{\mathbf{F}}_h$$

$$E=(U^2+V^2)/2+E_H+E_{\Sigma }$$

$$e\equiv E_H+E_{\Sigma }:=\frac{g}{2}|{\mathbf{F}}_h{|}^{-1}+\frac{G}{2}({\mathbf{F}}_h:{\mathbf{F}}_h+|{\mathbf{F}}_h{|}^{-2}).$$

$$D_t{\mathbf{A}}_h=({\mathbf{F}}_h^{-1}{\mathbf{F}}_h^{-T}-{\mathbf{A}}_h)/\lambda ,$$

$$D_t{A}_{cc}=(H^{-2}-A_{cc})/\lambda .$$

$$E_H+E_{\Sigma }=\frac{g}{2}H+\frac{G}{2}\left(\mathrm{tr}({\mathbf{B}}_h)-\ln \left(\det {\mathbf{B}}_h\right)+B_{zz}-\ln (B_{zz})\right),$$

$${\mathbf{B}}_h={\mathbf{F}}_h{\mathbf{A}}_h{\mathbf{F}}_h^T{B}_{zz}=H^2{A}_{cc}>0,$$

$$\boxed{\begin{array}{rl}& {\partial }_{t}H+\mathrm{div}(H\mathbf{U})=0\\ & {\partial }_t(H{\mathbf{F}}_h)+\mathrm{div}(H\mathbf{U}\otimes {\mathbf{F}}_h-H{\mathbf{F}}_h\otimes \mathbf{U})=0\\ & {\partial }_t(H\mathbf{U})+\mathrm{div}(H\mathbf{U}\otimes \mathbf{U}+(g\frac{H^2}{2}+G{H}^3{A}_{cc})\mathbf{I}-GH{\mathbf{F}}_h{\mathbf{A}}_h{\mathbf{F}}_h^T)=-KH\mathbf{U}\\ & {\partial }_t(H{\mathbf{A}}_h)+\mathrm{div}(H\mathbf{U}\otimes {\mathbf{A}}_h)=H({\mathbf{F}}_h^{-1}{\mathbf{F}}_h^{-T}-{\mathbf{A}}_h)/\lambda \\ & {\partial }_t(H{A}_{cc})+\mathrm{div}(H\mathbf{U}{A}_{cc})=H(H^{-2}-A_{cc})/\lambda \end{array}}$$

$${\partial }_t(H{\mathbf{A}}_h^{-2})+\mathrm{div}(H\mathbf{U}\otimes {\mathbf{A}}_h^{-2})=H{\mathbf{A}}_h^{-1}(\mathbf{I}-{\mathbf{A}}_h^{-1}{\mathbf{F}}_h^{-1}{\mathbf{F}}_h^{-T}+{\mathbf{F}}_h^{-T}{\mathbf{F}}_h^{-1}{\mathbf{A}}_h^{-1}){\mathbf{A}}_h^{-1})/\lambda ,$$

$${\partial }_t(H{A}_{cc}^{1/4})+\mathrm{div}(H\mathbf{U}{A}_{cc}^{1/4})=H(H^{-2}{A}_{cc}^{-3/4}-A_{cc}^{1/4})/4\lambda ,$$

$$\tilde{E}=(|\mathbf{U}{|}^2+gH)/2+G\left(\mathrm{tr}({\mathbf{F}}_h{\mathbf{A}}_h{\mathbf{F}}_h^T)+H^2{A}_{cc}\right)/2.$$

$$\mathcal{A}:=\{H>0,{\mathbf{A}}_h^{-1}={\mathbf{A}}_h^{-T}>0,A_{cc}^{-1}>0\}.$$

$${\nabla }_{H^{-1},A_{cc}^{-1}}^2{\tilde{E}}_1=\left(\begin{array}{cc}2g{H}^3+6G{H}^4{A}_{cc}& -2G{H}^3{A}_{cc}^{3/4}\\ -2G{H}^3{A}_{cc}^{3/4}& 2G{H}^2{A}_{cc}^{1/2}\end{array}\right).$$

$$\begin{gathered} \mathrm{tr}({\mathbf{H}}_{\theta }{\mathbf{Y}}_{\theta }^{\frac{1}{2}}({\mathbf{H}}_{\theta }{)}^T)>\mathrm{tr}\left({\mathbf{H}}_{\theta }(\theta {\mathbf{Y}}_1^{\frac{1}{2}}+(1-\theta ){\mathbf{Y}}_2^{\frac{1}{2}})({\mathbf{H}}_{\theta }{)}^T\right) \\ =\mathrm{tr}\left(({\theta }^2+\theta (1-\theta ){\mathbf{Y}}_1^{-\frac{1}{4}}{\mathbf{Y}}_2^{\frac{1}{2}}{\mathbf{Y}}_1^{-\frac{1}{4}}){\mathbf{F}}_1^T{\mathbf{F}}_1+((1-\theta {)}^2+\theta (1-\theta ){\mathbf{Y}}_2^{-\frac{1}{4}}{\mathbf{Y}}_1^{\frac{1}{2}}{\mathbf{Y}}_2^{-\frac{1}{4}}){\mathbf{F}}_2^T{\mathbf{F}}_2\right) \\ =\mathrm{tr}\left({\mathbf{F}}_{\theta }^T{\mathbf{F}}_{\theta }+\theta (1-\theta ){\mathbf{D}}_{\theta }^T{\mathbf{D}}_{\theta }\right)\ge \mathrm{tr}({\mathbf{F}}_{\theta }^T{\mathbf{F}}_{\theta }) \end{gathered}$$

$$E_{\Sigma }\equiv \frac{G}{2}\left(\mathrm{tr}({\mathbf{F}}_h{\mathbf{A}}_h{\mathbf{F}}_h^T)+H^2{A}_{cc}-\ln (|{\mathbf{F}}_h{|}^2\det {\mathbf{A}}_h{H}^2{A}_{cc})\right)$$

Full text

Precisely, the system of quasilinear PDEs is detailed for the shallow-water regime, i.e. for hydrostatic incompressible 2D flows with free surface under gravity. It generalizes Saint-Venant system to viscoelastic flows of Maxwell fluids, and encompasses previous works with F. Bouchut. It also generalizes the (thin-layer) elastodynamics of hyperelastic materials to viscous fluids, and to various rheologies between solid and liquid states that can be formulated using our new variable as material parameter.

The new viscoelastic flow model has many potential applications, additionally to falling into the theoretical framework of (symmetric hyperbolic) systems of conservation laws. In computational rheology, it offers a new approach to the High-Weissenberg Number Problem (HWNP). For transient geophysical flows, it offers perspectives of thermodynamically-compatible numerical simulations, with a Finite-Volume (FV) discretization say. Besides, one FV discretization of the new continuum model is proposed herein to precise our ideas incl. the physical meaning of the solutions. Perspectives are finally listed after some numerical simulations.

1 Introduction

Many mathematical models have been proposed for flows of real fluids. Like the celebrated Navier-Stokes equations [41], they mainly account for viscosity as a manifestation of the fluid non-ideality. But other macroscopic manifestations of the fluid microstructure in continuum mechanics can also be accounted for, typically through complex deviatoric stresses [61].

However, most of those models have a diffusion form ; thus they violate the physical principle of a finite speed of propagation for material/energetic perturbations. This violation may not be so important in some applications ; but it seems important to us e.g. in geophysics, when one is interested in the propagation of perturbations though large systems.

Moreover, a diffusion form also often entails unnecessary energy dissipation in numerical simulations ; this prevents one from capturing realistic transient dynamics on long time ranges.

A few hyperbolic systems of quasilinear PDEs have been proposed to model fluid flows with realistic shear-waves propagating at finite-speed. They rely on seminal ideas of Maxwell [45, 46, 47] after Poisson [59] to produce viscous (frictional/dissipative) effects through the relaxation of an elastic deformation (i.e. a source term produces entropy in shear flows).

For instance, [58, 24] proposes a model for slightly compressible 2D flows of viscoelastic fluids, see also [30] for recent numerical simulations. Although the latter model does not preserve mass, it is quite close to a model that we proposed in [10] and where mass conservation can be ensured, see [11]. But the latter model has a non-conservative form (see [11]) with genuinely nonlinear products that seem difficult to interpret for multi-dimensional flows. Another very interesting model has also been proposed recently in [57, 56], which aims at general flows of viscous fluids. Like in new our model, that very inspiring work introduces an additional “state variable” for viscoelastic matter. But the latter very recent proposition still deserves discussing with respect to its performance in applications, to its generalization to various rheologies and to its justification. Here, we propose an alternative similar in spirit, but mathematically different.

In the present work, we propose a multi-dimensional symmetric-hyperbolic system of conservation laws to model flows of (compressible) viscoelastic fluids of Maxwell type, which we believe new and useful for any flow and for various rheologies. The new flow model of Maxwell fluids is detailed for 2D isothermal flows of incompressible fluids with a free-surface and a hydrostatic pressure, see eq. (20) in Section 2. We have detailed our new viscoelastic approach in a low-dimensional case here, i.e. an analog of 2D Euler equations proposed e.g. by Saint-Venant [18] for shallow-water flows, because (i) our ideas are easier to understand theoretically and numerically then than in the most general case, (ii) that framework is useful already to a number of applications in geophysics (see below), and (iii) it is sufficient to see that our model generalizes to 3D flows and more complex fluids.

The new viscoelastic model contains our previous 1D model [9] as a closed subsystem for translation-invariant solutions, and our previous 2D models [10, 11] without a conservative formulation, see eq. (11) below. In particular, it contains the standard shallow-water system of Saint-Venant in the zero elasticity limit $G\rightarrow 0$. To encompass the limitations of our previous models, we have interpreted them as particular “closed” hyperbolic subsystems of the new model. We have generalized the elastodynamics system of hyperelastic materials to Maxwell fluids thanks to a new state variable $\boldsymbol{A}$ to that aim, and our new model thus also contains standard elastodynamics equations when $\boldsymbol{A}$ is uniform in space and time (i.e. at large Deborah, or Weissenberg number $\lambda\gg 1$).

Our new model has a molecular (“microscopic”) justification, formally at least, which is reminiscent of the molecular theory of polymer solutions [54, 23], see Rem. 1. The new state variable is one kind of “order parameter” which accounts for the flow-induced distortion of the microsctructure, relaxes due to thermal agitation at non-zero temperature and then creates viscous effects. This is thermodynamically compatible : one can formulate the second principle with an entropy functional, or with a Helmholtz free-energy that decays following a Clausius-Duhem inequality at non-zero temperature.

Our model copes well with a general formalism desired for the mathematical modelling of continuum mechanics according to the admitted physical principles [49]. As a consequence, our model straightforwardly benefits from the numerous efforts toward a precise mathematical understanding of (solutions to the Cauchy problem for) symmetric-hyperbolic systems of conservation laws [17], and toward accurate numerical simulations of physically-meaningful solutions after discretization e.g. by the Finite-Volume (FV) method [27, 40]. In particular, our model is close to well-known systems: the Euler gas dynamics system (formally equivalent to Saint-Venant shallow-water system when 2D and barotropic) and the elastodynamics system for hyperelastic materials [65].

Of course, our model is also limited by the present state of the mathematical theory (multi-dimensional solutions are not well defined yet, see [15] for instance), and of the FV method. For “entropy-stability”, one would like the free-energy to decay at the discrete level; however it is not globally convex with respect to the whole set of conservative variables (unlike the energy). Moreover, involutions like in the elastodynamics system for hyperelastic materials are important but hardly preserved discretely [65, 37].

However, systems of conservation laws are still much studied theoretically [19, 20, 16] and numerically [21]. So we believe our model has good perspectives. And even if at present, one has to pragmatically resort to heuristics and empiricism in numerical simulations, our model can already have many applications. For instance, it could be useful in environmental hydraulics to model complex fluid flows (turbulent/non-Newtonian) with generalized Saint-Venant equations. See in section 2.1 the modelling issue for real fluids in Saint-Venant 2D framework.

In section 2.2, our new Saint-Venant-Maxwell (SVM) model is introduced starting from the more standard SVUCM model for viscoelastic flows of Upper-Convected Maxwell fluids. We hope that our new viscoelastic approach solves a number of difficulties in computational rheology like the High-Weissenberg Numebr Problem (HWNP) [55], and generalization to compressible models [4] for thermodynamically-compatible mass transfers. In particular, we recall that the standard viscoelastic models usually require strict incompressibility with a non-zero retardation time (a background viscosity) [36, 55]. Now, whereas the standard approach can be useful numerically in computational rheology for small-data solutions, to avoid singularities [62], the trick does not seem to work well at high Weissenberg numbers [36, 55]; moreover, it naturally limits the scope of application.

In Section 3 we propose a Finite-Volume (FV) discretization of our new model, using a 1D Riemann solver of relaxation type (“à la Suliciu”) endowed with some discrete stability properties. Numerical approximations are naturally interesting for quantitative estimations in applications. But they are also necessary qualitatively at present, if one wants to precisely discuss solutions that are physically reasonable.

Note that our new model cannot be straightforwardly treated by the standard FV strategy recalled in section 3.1 for the sake of clarity. Indeed, (i) the mathematical entropy (i.e. the energy here, strictly convex with respect to conservative variables) should be replaced with a Helmholtz free-energy (a priori non convex with respect to conservative variables) to ensure a second-principle in presence of source terms here (i.e. at non-zero temperature), and (ii) our model (without source term) is a convex extension of a system with involutions, and we need the exact preservation of those involutions at the discrete level to ensure stability – which is a well-known difficulty in multi-dimensional numerical discretizations.

In section 3.2, we present a modification of the standard strategy. We suggest to consistently reconstruct at interfaces and in cells those (conservative) variables which do not correspond to fundamental conservation laws in physics (unlike mass, energy, momentum). Additionally, our numerical discretization gives some insight onto the physical meaning of the new system (with a view to selecting physically reasonable solutions, recall).

We draw a conclusion in Section 4 after showing some numerical results.

2 Viscoelastic Saint-Venant equations

Before introducing our viscoelastic Saint-Venant equations for Maxwell fluids, let us recall the usual Saint-Venant system of equations, for 2D flows of simpler fluids, i.e. the standard (nonlinear) shallow-water model for thin-layer (i.e. shallow) free-surface flows governed by a hydrostatic pressure.

2.1 Saint-Venant models for shallow free-surface flows

We consider a Eulerian description in a Galilean frame equipped with a Cartesian system of coordinates $(\boldsymbol{e}_{x},\boldsymbol{e}_{y},\boldsymbol{e}_{z})$, under a constant gravity field $\boldsymbol{g}=-g\boldsymbol{e}_{z}$. The shallow free-surface flows of homogeneous incompressible fluids above a flat impermeable plane $z=0$ can be modelled with a 2D velocity field $\boldsymbol{U}(t,x,y)=U^{x}(t,x,y)\boldsymbol{e}_{x}+U^{y}(t,x,y)\boldsymbol{e}_{y}$ for “infinitesimal fluid columns” under a non-folded free-surface $z=H(t,x,y)\geq 0$ using mass and momentum balance laws:

[TABLE]

where $\boldsymbol{I}$ denotes the identity second-order tensor.

Assuming the pressure hydrostatic $p(t,x,y,z)\approx g(H-z)$ so $P=gH/2\approx\frac{1}{H}\int_{0}^{H}dzp$ in (2), the 2D velocity $\boldsymbol{U}=(U,V)$ is interpreted as a depth-average

[TABLE]

of the horizontal components of a 3D velocity $\boldsymbol{u}=u^{x}\boldsymbol{e}_{x}+u^{y}\boldsymbol{e}_{y}+u^{z}\boldsymbol{e}_{z}$, possibly the mean of a statistically stationary field e.g. in turbulent flows.

The “stress” term111 Note that $P$, $\Sigma_{zz}$ and the components $\Sigma_{h,ij}$ have the dimension of energies per unit mass rather than per unit volume here, unlike standard Cauchy stresses and pressures: one could term them specific forces. But for the sake of simplicity, we will omit the label “specific”, insofar as there is no ambiguity here. $-\Sigma_{zz}\boldsymbol{I}+\boldsymbol{\Sigma}_{h}$ is usually symmetric, i.e. $\Sigma_{yx}=\Sigma_{xy}$ in

[TABLE]

and it accounts, depending on the closure, for:

• •

depth-averaged Cauchy stresses when

[TABLE]

(the case of Newtonian fluids with constant viscosity $\nu>0$ see e.g. [26, 43] for a justification based on asymptotic analysis of the depth-averaged Navier-Stokes equations for Newtonian fluids), or

• •

empirical corrections in the depth-averaged acceleration terms

[TABLE]

when the horizontal components $(u^{x},u^{y})$ of the (time-averaged) velocity do not have a uniform profile in vertical direction222 The effect is called dispersion in hydraulics. , or when the second-order moments of a turbulent velocity 333 They are usually termed turbulent stresses. cannot be neglected in the dynamics of the (depth-averaged) mean velocity field $\boldsymbol{u}$.

New 2D models keep being developped for turbulent flows (see e.g. [63, 25, 11]) as well as for real non-Newtonian fluids (see e.g. [5]).

One approach to the modelling of non-Newtonian stresses is to start with 3D models and next close 2D depth-averaged models using scaling assumptions (like [26, 43] for Newtonian fluids, [50] or our previous work [10] otherwise).

But the 3D viscoelastic flows models still raise many questions (see Introduction), and we believe a direct 2D analysis in the shallow framework could be a useful alternative first step toward better 3D models. Let us recall precisely the modelling issue, which still holds for (1),(2), or its more usual formulation with source terms modelling friction over rugous bottom in real application [31]:

[TABLE]

whose smooth physically-meaningful solutions are expected to satisfy an additional conservation law for a free-energy $E=\frac{1}{2}\left(|\boldsymbol{U}|^{2}+e\right)$ with $D\geq 0$ :

[TABLE]

or equivalently for a specific internal energy $e$:

[TABLE]

In the pure viscous case (4) (Newtonian fluids), it is well-known that (5) holds as equality with $D=2\nu\left(|D(\boldsymbol{U})|^{2}+2|\operatorname{div}\boldsymbol{u}|^{2}\right)\geq 0$ and $e=gH$ for smooth solutions with a Helmholtz free-energy $E$ termed the mechanical energy

[TABLE]

while (5) only holds as inequality ($\leq$) for weak solutions that implicitly model irreversible flows [64].

In the ideal case $\Sigma_{zz}=0,\boldsymbol{\Sigma}_{h}={\boldsymbol{0}}$ for the widely-used inviscid shallow-water model (or Euler isentropic 2D flow model for perfect polytropic gases with $\gamma=2$), the equality (5) also holds for smooth solutions to the symmetric hyperbolic quasilinear system (note Godunov-Mock theorem applies since (7) is a strictly convex function of $(H^{-1},\boldsymbol{U})\in\mathbb{R}_{>0}\times\mathbb{R}^{2}$) see e.g. [27, 1]. Moreover, on requiring the inequality associated with (5) ($D=0$ in the ideal case):

[TABLE]

one can define physically-admissible entropy solutions where irreversible processes dissipate the mechanical energy (7), that are unique in 1D within a translation-invariant solution class (see e.g. [3] when $H>0$ and the shallow-water system is strictly hyperbolic).

None of the two standard cases above model well the propagation of shear stress/strain at finite-speed444 The propagation of information at finite-speed is not only physical. It also allows for a precise computation of fronts, that can be compared with experimental observations. This is the reason why the inviscid model is more often used in practice, e.g. by hydraulic engineers, although it cannot sustain shear motions like the viscous model (i.e. a non-trivial steady state in shear flows). . So let us look for a hyperbolic quasilinear model of Saint-Venant type (a shallow-water model) for 2D viscoelastic flows. Precisely, we look an additional law like (5) that properly defines a notion of viscosity and accounts for vortices in stationary flows.

2.2 Saint-Venant models generalized to Maxwell fluids

Viscoelastic fluids of Maxwell type are characterized by an elasticity modulus $G$ (in stress units) and a finite relaxation time scale555 Typically characterizing the time needed for the stress in a fluid initially at rest to relax to a viscous state (proportional to strain-rate by the viscosity factor $\nu=G\lambda>0$) after suddenly straining the fluid at a fixed maintained rate. $\lambda$ (termed Weissenberg number when non-dimensionalized with a time scale of the flow like $|\boldsymbol{\nabla}\boldsymbol{U}|^{-1}$).

Closure formulas for the non-Newtonian Cauchy stresses of Maxwell fluids can be obtained e.g. following the same depth-averaged analysis as in [26, 43], starting with full 3D models for free-surface flows of Maxwell fluids see e.g. [10]. When the time rate of change for the stress tensor $\boldsymbol{\Sigma}$ is the upper-convected derivative (i.e. we consider Upper-Convected Maxwell fluids), one obtains

[TABLE]

without background viscosity (i.e. with zero retardation time in viscoelastic terminology) in some asymptotic regime, denoting $\boldsymbol{L}_{h}$ the horizontal velocity gradient in (9) which corresponds to the horizontal stress tensor (3).

Viscous stresses with a viscosity $\nu>0$ arise from (9),(10) in (2) when $G=\lambda\nu$, $\lambda\rightarrow 0$. On the contrary, when $\lambda\rightarrow\infty$, the fluid becomes purely elastic and governed by the homogeneous quasilinear system (11) which is similar to a thin-layer approximation of the elastodynamics system governing the deformations of a Hookean hyperelastic continuum (see details below).

With $\boldsymbol{B}_{h}:=\lambda\boldsymbol{\Sigma}_{h}/\nu+\boldsymbol{I}$ and $B_{zz}:=\lambda\Sigma_{zz}/\nu+1$, (1),(2),(9),(10) rewrites as the hyperbolic quasilinear system (11) (see App. B for details and a proof):

[TABLE]

which is similar to another hyperbolic system for viscoelastic flows already known in the literature, see Rem. 3.

But hyperbolicity is only a necessary condition for the definition of a sensible initial-value problem with a quasilinear system (in a very weak sense, see e.g. [1] in the particular case of constant coefficients), it is not sufficient. The initial-value problem can be shown well-posed (on small times $t\in[0,T)$ for smooth initial data) for the symmetric hyperbolic quasilinear systems [35, 42]. For instance, the systems of conservation laws are symmetric hyperbolic when they are endowed with an additional conservation law for a strictly-convex functional termed “entropy” [28, 27].

Although the quasilinear SVUCM system is endowed with an additional conservation law for a convex free-energy functional (see $E$ below in (14)), (1),(2),(11) is obviously not in conservation form. So, the additional conservation law is not useful, and the meaning of weak solutions remains unclear.

One could try to make sense of SVUCM as such with non-conservative products like in [38, 2]. However, note that SVUCM is unlikely to possess “generalized symmetrizers”, see the analysis in [53] for a close system in 2D, while “generalized symmetrizability” seems a minimum requirement (see e.g. [48]) to define a meaningful concept of solution (like dissipative measure-valued solutions satisfying a weak-strong uniqueness principle as in [20]).

Besides, let us also recall that the question how to correctly formulate (multidimensional) equations for flows of non-Newtonian fluids is not settled in general, especially for compressible viscoelastic flows [4], although it has received a number of answers (the close 2D system discussed in [53] see also Rem. 3 is for slightly compressible viscoelastic fluids).

That is why we propose here to consider a modification of (1),(2),(11) that bears the same physical meaning but conforms with the existing theory for symmetric hyperbolic system of conservation laws. Our strategy is to devise an enlarged system of conservation laws that contains (1),(2),(11).

To that aim, observe first that when $\lambda\rightarrow\infty$, $\boldsymbol{B}_{h}$ and $B_{zz}$ in fact identify with the horizontal and vertical components (resp.) of the Cauchy-Green (left) deformation tensor $\boldsymbol{B}=\boldsymbol{F}\boldsymbol{F}^{T}$ in a hyperelastic incompressible homogeneous continuum, where $\boldsymbol{F}=\partial_{a,b,c}(x,y,z)$ is the deformation gradient with respect to a reference configuration in the Cartesian coordinate system $(\boldsymbol{e}_{a},\boldsymbol{e}_{b},\boldsymbol{e}_{c})$.

Precisely, the Cauchy stress terms in the momentum equation rewrite

[TABLE]

when $H\equiv F^{z}_{c}=|\boldsymbol{F}_{h}|^{-1}>0$. So, when $\lambda\rightarrow\infty$, eq. (11) in the SVUCM system is a consequence of the elastodynamics of Hookean incompressible materials with uniform mass density in the shallow-water regime, using

[TABLE]

In fact, the SVCUM system (1),(2),(11) coincides with the elastodynamics of a 2D hyperelastic continuum, with $\boldsymbol{B}_{h}\equiv\boldsymbol{F}_{h}\boldsymbol{F}_{h}^{T}$ computed from $\boldsymbol{F}_{h}$

[TABLE]

when $\lambda\rightarrow\infty$, recall (12) as well as $B_{zz}\equiv|\boldsymbol{F}_{h}|^{-2}$ by incompressibility. The energy conservation law (5) is satisfied by the Helmholtz polyconvex free-energy

[TABLE]

with $D=0$ (recall $H=|\boldsymbol{F}_{h}|^{-1}$; $\lambda\rightarrow\infty$ means no source) and internal energy

[TABLE]

This is consistent with the fact that Cauchy stress term $H\boldsymbol{\Sigma}_{h}-H\boldsymbol{\Sigma}_{zz}\boldsymbol{I}\equiv GH(\boldsymbol{B}_{h}-B_{zz}\boldsymbol{I})$ in the momentum equation actually equals $H(\partial_{\boldsymbol{F}_{h}}e)\boldsymbol{F}_{h}^{T}$.

So, at least when $\lambda\rightarrow\infty$, one can make sense of SVUCM system to model time-evolutions. SVUCM coincides with the 2D elastodynamics of a hyperelastic incompressible materials, which is equipped with a polyconvex energy like (15) and a symmetric hyperbolic conservative formulation (see e.g. [17, 65]): the Cauchy problems are well-posed on small times given smooth initial conditions.

For the general case with $\lambda>0$, we now propose to embed SVUCM into a quasilinear system of conservation laws for the 2D visco-elastodynamics of a hyperelastic incompressible continuum with memory. To that aim, we introduce new state variables666 One may also want to call them “internal” variables. $\boldsymbol{A}_{h}=\boldsymbol{A}_{h}^{T}>0$, $A_{cc}>0$ that account for “viscous” deformations of the microstructure in a coordinate system attached to the reference configuration, which is reminiscent of the distortion metric used in elasto-plasticity [51, 39]. We postulate simple constitutive laws (see Remark 1):

[TABLE]

Then, using in (14) the internal energy of SVUCM (recall [67, 9] e.g.)

[TABLE]

with horizontal (symmetric) and vertical positive strains defined by

[TABLE]

the mass, momentum and energy conservation laws lead to the system:

[TABLE]

which we term Saint-Venant-Maxwell or SVM in short.

Proposition 2.1**.**

The quasilinear system of conservation laws (20) written for $(H,H\boldsymbol{U},H\boldsymbol{F}_{h},HA_{cc}^{1/4},H\boldsymbol{A}_{h}^{-2})$ with ($\boldsymbol{A}_{h}^{-1}$ is $\boldsymbol{A}_{h}^{-2}$ square-root matrix):

[TABLE]

is equipped (on neglecting the source) with a mathematical entropy $H\tilde{E}$ where

[TABLE]

*It is therefore symmetric hyperbolic on the convex admissibility domain *

[TABLE]

Proof.

Since $H\tilde{E}$ obviously satisfies (5) in the case without source, it suffices to show that (20) is symmetric hyperbolic on $\mathcal{A}$ (obviously convex), i.e. that $H\tilde{E}$ is stricly convex with respect to a full set of conserved variables, using Godunov-Mock theorem [27]. Moreover, $H\tilde{E}$ is (strictly) convex in $(H,H\boldsymbol{U},H\boldsymbol{F}_{h},HA_{cc}^{1/4},H\boldsymbol{A}_{h}^{-2})$ if and only if $\tilde{E}$ is (strictly) convex in $(H^{-1},\boldsymbol{U},\boldsymbol{F}_{h},A_{cc}^{1/4},\boldsymbol{A}_{h}^{-2})$ [6], i.e. if $\tilde{E}_{1}=gH+G\left(H^{2}A_{cc}\right)$ and $\tilde{E}_{2}=\operatorname{tr}(\boldsymbol{F}_{h}\boldsymbol{A}_{h}\boldsymbol{F}_{h}^{T})$ are (strictly) convex in $(H^{-1},A_{cc}^{1/4})$ and $(\boldsymbol{F}_{h},\boldsymbol{A}_{h}^{-2})$ respectively, like $|\boldsymbol{U}|^{2}/2$ in $\boldsymbol{U}$. Now, the smooth function $\tilde{E}_{1}$ is strictly convex insofar as its Hessian matrix is strictly positive:

[TABLE]

On the other hand, consider two couples of matrix values $(\boldsymbol{F}_{1},\boldsymbol{Y}_{1}:=\boldsymbol{A}_{1}^{-2})$ and $(\boldsymbol{F}_{2},\boldsymbol{Y}_{2}:=\boldsymbol{A}_{2}^{-2})$ for $(\boldsymbol{F}_{h},\boldsymbol{A}_{h}^{-2})$. for any $\theta\in[0,1]$, using $\boldsymbol{F}_{\theta}=\theta\boldsymbol{F}_{1}+(1-\theta)\boldsymbol{F}_{2}$, $\boldsymbol{Y}_{\theta}=\theta\boldsymbol{Y}_{1}+(1-\theta)\boldsymbol{Y}_{2}$ and $\boldsymbol{H}_{\theta}=\theta\boldsymbol{F}_{1}\boldsymbol{Y}_{1}^{-\frac{1}{4}}+(1-\theta)\boldsymbol{F}_{2}\boldsymbol{Y}_{2}^{-\frac{1}{4}}$, $\boldsymbol{D}_{\theta}\boldsymbol{F}_{1}\boldsymbol{Y}_{1}^{-\frac{1}{4}}\boldsymbol{Y}_{2}^{\frac{1}{4}}-\boldsymbol{F}_{2}\boldsymbol{Y}_{2}^{-\frac{1}{4}}\boldsymbol{Y}_{1}^{\frac{1}{4}}$ it holds

[TABLE]

hence $\operatorname{tr}(\boldsymbol{H}_{\theta}(\boldsymbol{H}_{\theta})^{T})>\operatorname{tr}(\boldsymbol{F}_{\theta}\boldsymbol{Y}_{\theta}^{-\frac{1}{2}}\boldsymbol{F}_{\theta}^{T})$ since $\boldsymbol{Y}_{\theta}^{-\frac{1}{2}}$ is symmetric positive definite, which is the desired result $\theta E_{2}(\boldsymbol{F}_{1},\boldsymbol{Y}_{1})+(1-\theta)E_{2}(\boldsymbol{F}_{2},\boldsymbol{Y}_{2})>E_{2}(\boldsymbol{F}_{\theta},\boldsymbol{Y}_{\theta})$. ∎

Proposition 2.1 allows one to check that the system (20) makes sense as a flow model with any smooth source term, for small times at least. But in fact, the full SVM system (20) can also be shown thermodynamically compatible i.e.:

Corollary 2.1**.**

Given smooth initial conditions, the Cauchy problems for strong solutions to (20) are well-posed. These strong solutions preserve the relation $H\equiv F^{z}_{c}=|\boldsymbol{F}_{h}|^{-1}$. They also satisfy the companion conservation law (5) for the same free-energy (14) as SVUCM where the usual hydrostatic potential term $E_{H}=gH/2$ of Saint-Venant is complemented by the viscoelastic term

[TABLE]

and by thermodynamically compatible source terms, dissipating energy at rate:

[TABLE]

Proof.

Corollary 2.1 is a consequence of Proposition 2.1: the well-posedness of Cauchy problems for strong solutions to symmetric hyperbolic systems is classical, see e.g. [1]. Then, given smooth solutions, one can directly check that $|\boldsymbol{F}_{h}|^{-1}$ follows the same evolution as $H$. The SVUCM system can be retrieved exactly using (19), (16), (17). This is the reason why the same companion conservation law (5) holds for (20) as for SVUCM, where the source term $D$ is thermodynamically compatible – and also an upper bound for $E_{\Sigma}$, thus $E$ ! ∎

Some comments about the structure of SVM are now in order.

First, note that we have not been able to find a full set of conserved variables such that the free-energy $HE$ is convex strictly on the whole admissible domain $\mathcal{A}$ when it is defined as in (14) by $e\equiv E_{H}+E_{\Sigma}$, $E_{H}=\frac{g}{2}H$ and (24) for $E_{\Sigma}$. This is why we use $H\tilde{E}$ rather than $HE$ to show that the SVM system is symmetric hyperbolic and thus a good model, with well-posed Cauchy problems. But we insist on the importance of $HE$ with $E_{\Sigma}$ defined by (24). It allows one to show that the source terms in SVM which formally yield Navier-Stokes equilibrium asymptotically when $\lambda\rightarrow 0$ are thermodynamically compatible. A result like Corrolary 2.1 is necessary, to justify weak solutions satisfying the inequality:

[TABLE]

as second thermodynamics principle complementing (20) for time-evolution with irreversible processes (the non-convexity of the free-energy is physical in presence of a source which modifies the energy landscape at a fixed temperature).

Second, recall that the SVUCM system can be retrieved exactly from (smooth solutions) to SVM using (19), (16), (17) (cf. Corr. (2.1)). Now, the SVUCM system with $\boldsymbol{B}_{h}\equiv\boldsymbol{F}_{h}\boldsymbol{A}_{h}\boldsymbol{F}_{h}^{T}$, $B_{zz}=H^{2}A_{cc}$ is a closed Galilean-invariant quasilinear system which is hyperbolic (strongly, see App. B). The Jacobian matrix (say for 1D waves along $\boldsymbol{e}_{x}$) is diagonalizable. The 7 eigenvalues

[TABLE]

with as many eigenvectors are also 7 eigenvalues of the extension termed SVM. Moreover, the eigenvalue $\lambda_{0}=U$ is of multiplicity $5$ at least for SVM, when using $\boldsymbol{A}_{h}$ as independent variable ($B_{zz}$ is equivalent to $A_{cc}$ as long as $H=|\boldsymbol{F}_{h}|^{-1}$). In fact, being symmetric hyperbolic, our SVM system above is also strongly hyperbolic when “freezing” the Jacobian: it is easily seen that [math] is the remaining (real) eigenvalue of SVM (above, in Eulerian coordinates), with multiplicity 2. So, in general, resonance can occur [33, 7] as well as numerical difficulties like for systems with discontinuous flux, see e.g. [52]. This difficulty was not apparent in our previous work on the (closed) 1D subsystem [9] but it is inline with standard results in multidimensional elastodynamics (the case $\boldsymbol{A}_{h}=\boldsymbol{I}$, $A_{cc}=1$), see e.g. [65, 37] and references therein, and possibly with the numerical difficulties observed close to vacuum $H=0$ in our previous 1D works on SV(UC)M [12].

Last, for relevant applications to a specific context, one may still want to precise the physical interpretation of the new state variable $\boldsymbol{A}$. At present, we think of it as a material parameter that simply follows the flow at “zero-temperature”, and that relaxes to balance mechanical deformations otherwise. For application to polymer suspensions, one may think of it as a mean-field approximation, see Rem. 1 for a “micro-macro” interpretation. More generally, $\boldsymbol{A}$ is some kind of viscous strain measuring the distortion of the volume elements (or equiv. microstructure deformations) in a coordinate system attached to a reference configuration. This is mathematically different from, but similar in spirit to, the tensor state variable $\boldsymbol{A}$ introduced in [57] (governed by the non-conservative equations (9),(10) like the elastic strain measure $\boldsymbol{B}=\boldsymbol{F}\boldsymbol{A}\boldsymbol{F}^{T}$ there). In any case, the fact that $\boldsymbol{A}$ describes material properties of the flows can be entlightened on writing SVM with a Lagrangian description, when a bijective flow map $\phi:(t,x,y)\rightarrow(t,a,b)$ between the current and reference configurations is well-defined with $\boldsymbol{U}\equiv\partial_{t}\phi^{-1}$ and $\boldsymbol{F}\equiv\boldsymbol{\nabla}\phi^{-1}$. Recall

[TABLE]

is the Eulerian description (20) for $(H,F^{i}_{\alpha},U^{i},A_{\alpha\beta},A_{cc})$ ($i,j\in\{x,y\};\alpha,\beta\in\{a,b\}$) in a Cartesian basis where $\boldsymbol{F}_{h}=F^{i}_{\alpha}\boldsymbol{e}_{i}\otimes\boldsymbol{e}_{\alpha}$ and $\boldsymbol{A}_{h}=A_{\alpha\beta}\boldsymbol{e}_{\alpha}\otimes\boldsymbol{e}_{\beta}$ so, denoting $\sigma_{xy}=1=-\sigma_{yx}$, it holds $\boldsymbol{F}_{h}^{-1}=|\boldsymbol{F}_{h}|^{-1}(\sigma_{ij}\sigma_{\alpha\beta}F^{j}_{\beta})\boldsymbol{e}_{\alpha}\otimes\boldsymbol{e}_{i}$ with $|\boldsymbol{F}_{h}|=\sigma_{ij}\sigma_{\alpha\beta}F^{i}_{\alpha}F^{j}_{\beta}$. Then, an equivalent777In the case of sufficiently smooth solutions [66, 65]. Lagrangian description holds using

[TABLE]

if $H|\boldsymbol{F}_{h}|=1$ i.e. the so-called Piola’s identities, which reads :

[TABLE]

while (26) for $E=\frac{1}{2}\left(\sum_{i}|U^{i}|^{2}+gH+GF^{i}_{\alpha}A_{\alpha\beta}F^{i}_{\beta}+GH^{2}A_{cc}-\log(A_{cc}|\boldsymbol{A}_{h}|)\right)$

[TABLE]

(simplified with $H|\boldsymbol{F}_{h}|=1$) also has a conservative Lagrangian equivalent:

[TABLE]

with the same dissipation $D>0$ given in (25). Introducing $\mathcal{G}^{i}_{\alpha}=GF^{i}_{\beta}A_{\alpha\beta}$,

[TABLE]

we obtain a simple reformulation of (the 3 first lines of) (29) and (31) as:

[TABLE]

[TABLE]

which can now be easily compared to the usual Lagrangian formulation of elastodynamics [17, 65] (see Rem. 2): $G\boldsymbol{A}_{h}$, $GA_{cc}$ can be understood as variable anisotropic elastic properties, which induce a viscous behaviour through friction on a time-scale $\lambda\rightarrow 0$ inline with Maxwell ideas [46, 47, 59]. Furthermore, like in standard elastodynamics, the Lagrangian equations above should be useful for variational calculus with SVM (see e.g. [44, Lecture 2]) as well as for numerical approximation (see e.g. [32, 22] and our last section). Note however that, also like in standard elastodynamics, the link between the Eulerian and Lagrangian formulations requires $H|\boldsymbol{F}_{h}|=1$ as well as (28). While this is known to hold for all times if it holds initially at a continuous level (this is called an involution), it is a well-known difficulty at a discrete level [37].

Remark 1**.**

*The system (20) cannot be retrieved from the standard kinetic interpretation in statistical physics, when non-Newtonian stresses $\boldsymbol{\Sigma}_{h}=G(\boldsymbol{B}_{h}-\boldsymbol{I})$ and $\Sigma_{zz}=G(B_{zz}-1)$ are due to Brownian elastic “dumbbells” diluted in a fluid suspension with conformation matrix $\boldsymbol{B}=\mathbb{E}(\boldsymbol{R}\otimes\boldsymbol{R})$, $\boldsymbol{R}(t,\boldsymbol{x})$ being a random end-to-end vector solution to an overdamped Langevin equation [54, 23]. However, in a similar spirit, one could interpret $\boldsymbol{A}_{h}$ and $A_{cc}$ as the mean-field approximations $\mathbb{E}(\boldsymbol{G}_{h}\boldsymbol{G}_{h}^{T})$ and $\mathbb{E}(|G^{c}_{z}|^{2})$ of stochastic processes : *

[TABLE]

[TABLE]

*modelling (relaxation of) the elastic microstructure with thermal fluctuations. *

Remark 2** **(Direct derivation of the

model in Lagrangian description).

The Lagrangian formulation (29) of our model could also be straightforwardly derived as the 2D visco-elastodynamics of a hyperelastic incompressible continuum with internal energy (18) in a coordinate system attached to a reference configuration (a Lagrangian description), see e.g. [65, 17, 60]. Like the Eulerian formulation, it can be seen as a formal thin-layer approximation of the 3D (visco)elastodynamics of an incompressible continuum with memory and a free-surface when the terms containing $F^{z}_{\alpha}$ are negligible: the shallow-water regime.

Remark 3** (Hyperbolic models of viscoelastic flows).**

In introduction, we have already mentionned that we were aware of a few other interesting hyperbolic models of viscoelastic flows for Maxwell fluids. The 2D model for slightly compressible flows in [58, 24], see also [53, 30, 29] for simulations, is similar to (11) (though it does not conserve mass). It is hyperbolic under the same physically-natural conditions $H,B_{zz},\operatorname{tr}\boldsymbol{B}_{h},\operatorname{det}\boldsymbol{B}_{h}\geq 0$. Another hyperbolic model for 3D flows of Maxwell fluids is in [57]. Although similar in spirit, noet that it does not seem to explicitly compare with our new proposition. Last, note that one can derive another hyperbolic closed subsystem of SVM, similar but different from the SVUCM equations (11), on choosing $-\boldsymbol{F}_{h}$ as a conservative variable in (20). The energy functional remains the same however. Without source terms, that model was discovered independently by [63]. It can be interpreted as a building block for Reynolds-avergaged turbulence models in Saint-Venant framework (see SVTM in [11], and references therein).

3 Discretization by a Finite-Volume method

To investigate quantitatively the features of (20) as a model for viscoelastic flows (under gravity in the shallow-water regime), we need to define precisely how to compute solutions. To that aim, we consider discrete solutions to Cauchy problems using a Finite-Volume (FV) method [27, 40].

The standard FV strategy requires one to handle 1D Riemann problems (to define numerical fluxes), but this is not precise enough. Indeed, the discrete FV fields $H,\boldsymbol{F}_{h},\boldsymbol{U},\boldsymbol{A}_{h},A_{cc}$ (piecewise-constants on the polygonal cells $V_{i}$, $i\in\mathbb{N}$ of 2D meshes) should be “stable enough” and satisfy not only the conservation laws (20), but also the preservation of the domain $H,\boldsymbol{A}_{h}=\boldsymbol{A}_{h}^{T},A_{cc}>0$ and the physically-meaningful inequality (26).

In section3.1 we recall how to standardly build stable Riemann-based FV approximations for a quasilinear model

[TABLE]

like (20) when a companion law with a dissipation888 Typically induced by a thermodynamically-compatible source term like $B_{i}(q)=(q_{i}^{\infty}(q)-q_{i})/\lambda_{i}$ where $q^{\infty}(q)$ lies in the convex domain for $q$, and $\lambda_{i}(q)>0$ [14]. $D(q)\equiv-\nabla_{q}S(q)\cdot B(q)\geq 0$

[TABLE]

holds for a true mathematical entropy $S(q)$ convex in the Galilean-invariants of the state $q$ on the whole domain $H,\boldsymbol{A}_{h}=\boldsymbol{A}_{h}^{T},A_{cc}>0$ [27, 13]. Typically:

• •

the Cauchy problems are numerically solved by time-splitting: at each time step, the nonlinear flux terms are computed first by a forward method, while source terms are computed last backward in time,

• •

in the first (forward) fractional step one uses 1D Riemann solutions e.g. with relaxed conservation laws such that the inequality (37) is approximated at the same time as (36).

Such “standard” FV approximations are fully computable, and they converge to the smooth solutions on small times at least [13, 34]. Despite the lack of well-defined global solutions, they usually allow one to numerically explore a hyperbolic system of conservation laws in some useful regimes thanks to a physically-based guarantee of stability999 Even though they do not converge in all cases, and some meaningul solutions are not captured. In particular, we are aware that such FV approximations are likely to remain consistent only for small times, which is a problem to capture e.g. steady states when $\lambda\ll 1$. Asymptotics preserving schemes are however more involved and will be studied later in future works, in a second numerical exploration of our model. .

But in our case, (26) plays the role of (37) with $S$ replaced by $HE$ which is not convex, so it seems we cannot use the standard procedure.

Therefore, after recalling the standard case for the sake of clarity, we present a modification of the standard FV strategy in Section 3.2, which relies on a previous analysis of the Lagrangian reformulation as it is usual for Eulerian systems (36) like SVM that possess contact-discontinuity waves (see e.g. [7]).

3.1 Finite-Volume approach to standard conservation laws

Given a tesselation of $\mathbb{R}^{2}$ using polygonal cells $V_{i}$ ($i\in\mathbb{N}$), consider first $q_{h}(t)=\sum_{i}q_{i}(t)1_{V_{i}}$ a semi-discrete FV approximation of $q$ solution to a Cauchy problem for (36) with $q_{h}(0)\equiv q_{h}^{0}\approx q(0)$ at $t=0\equiv t^{0}$. To define $q_{h}(t)$ at $t>0$, integration from time $t^{n}$ to $t^{n+1}=\sum_{k=0}^{n-1}\tau^{k}$ ($\tau^{k}>0$) for standard systems of conservation laws like (36) is usually splitted into two sub-steps as follows.

First, for each $n\in\mathbb{N}$, one integrates the flux terms forward so $q_{h}^{n+1,-}$ approximates the solution at $t^{n+1,-}$ to the Cauchy problem with $q_{h}(t^{n})$ at $t^{n}$ for

[TABLE]

on $[t^{n},t^{n+1})$. Denoting $\boldsymbol{n}_{i\rightarrow j}$ the unit normal from $V_{i}$ to $V_{j}$ at $\Gamma_{ij}\equiv\overline{V_{i}}\cap\overline{V_{j}}$, domain-preserving 1D Riemann solvers allow one to define admissible FV approximations $q^{n+1,-}_{h}=\sum_{i}q^{n+1,-}_{i}1_{V_{i}}$, i.e. in the domain of $q$, by the formula

[TABLE]

through numerical fluxes $\boldsymbol{F}_{i\rightarrow j}(q_{i}^{n},q_{j}^{n};\boldsymbol{n}_{i\rightarrow j})$ precised below, under a CFL condition on $\tau^{n}$ (see Prop. 3.1). Moreover, one gets a fully admissible $q_{h}(t^{n+1})=\sum_{i}q^{n+1}_{i}1_{V_{i}}$, which is also consistent with the dissipation inequality associated with (37) in addition to (38), when $\boldsymbol{F}_{i\rightarrow j}$ use entropy-consistent Riemann solvers (see Prop. 3.2). Second, source terms are integrated backward (see Prop. 3.3):

[TABLE]

Recall that by Galilean invariance, in (39), standard numerical fluxes read

[TABLE]

with $\boldsymbol{O}_{i,j}^{-1}q$ in a local basis $(\boldsymbol{n}_{i\rightarrow j},\boldsymbol{n}_{i\rightarrow j}^{\perp})$ rather than in $(\boldsymbol{e}_{x},\boldsymbol{e}_{y})$, and with

[TABLE]

defined simply with a 1D Riemann solver $R(\xi,\tilde{q}_{i},\tilde{q}_{j})$, i.e. a well-defined solution to the 1D Riemann problem for (38) with initial condition $\tilde{q}_{i}1_{a<0}+\tilde{q}_{j}1_{a>0}$ on $\mathbb{R}\ni a$, which is a function of $\xi=a/t$, or some conservative approximation (termed simple in the latter case when $R(\cdot,\tilde{q}_{i},\tilde{q}_{j})$ is piecewise constant).

Recall also that information propagates at finite speed in explicit FV approximation $q_{h}^{n+1,-}$, like in (38). Moreover, that speed is consistent when the maximal speed $s(q_{l},q_{r})>0$ of the waves in $R(\cdot,q_{l},q_{r})$ is bounded continuously as a function of $q_{l},q_{r}$. And $q_{h}^{n+1,-}$ is admissible if the Riemann solver $R(\cdot,q_{l},q_{r})$ in (42) preserves the domain of $q$ (i.e. the values assumed by $\xi\rightarrow R(\xi,q_{l},q_{r})$ belong to the admissibility domain for $q$) under CFL condition [7]:

Proposition 3.1**.**

If a numerical flux is given by (41), (42) with a 1D Riemann solver $R(\cdot,q_{l},q_{r})$ that preserves a convex domain for $q$ and has bounded maximal wavespeed $s(q_{l},q_{r})>0$, then the FV approximate solution (39) to (38) also preserves the domain for $q$ under the CFL condition (43)

[TABLE]

Proof.

It suffices to rewrite (39) with (42) as

[TABLE]

i.e. as a convex combination under the CFL condition (43). ∎

A direct consequence of the admissibility of $q_{h}^{n+1,-}$ is that $q_{h}(t^{n+1})$ computed by (40) is also admissible, in the domain of $q$. Moreover it holds

[TABLE]

i.e. the (convex) mathematical entropy necessarily decreases. But this is not enough yet for $q_{h}(t^{n+1})$ to be fully admissible, i.e. to approximate the inequality

[TABLE]

as an admissibility criterion formulating the thermodynamics second principle.

Next, if $R(\cdot,q_{l},q_{r})$ in (41), (42) is also entropy-consistent (with (46) in 1D, see (49) below) then, under a CFL condition more stringent than (43), a discrete version of (46) holds when $k=0=D$ (Prop. 3.2), and when $k\geq 0,D\geq 0$ after backward integration of the sources (Prop. 3.3):

Proposition 3.2**.**

If the flux (41), (42) uses a 1D Riemann solver $R(\cdot,q_{l},q_{r})$ that preserves the domain of $q$ and is entropy-consistent with (46) when $k=0=D$ in the sense that, given admissible state-vectors $q_{l},q_{r}$ and a direction $\boldsymbol{n}$, a discrete entropy-flux vector $\tilde{G}_{\boldsymbol{n}}(q_{l},q_{r})=-\tilde{G}_{\boldsymbol{n}}(q_{r},q_{l})$ satisfies

[TABLE]

then, under the CFL condition (48)

[TABLE]

where $s_{i}^{n}:=\max_{j}s(\boldsymbol{O}_{i,j}^{-1}q_{i}^{n},\boldsymbol{O}_{i,j}^{-1}q_{j}^{n})$, the FV approximation (39) preserves the domain of $q$ and satisfies the following discrete version of (46) (with $D=0$):

[TABLE]

Proof.

We follow the 1D proof in [7] and first rewrite (39) with (42) as

[TABLE]

i.e. as a convex combination that depends only on the Riemann solver under the stringent CFL condition (48). We can now use Jensen inequality with (50), and next (47) with $\int\displaylimits_{0}^{+\infty}S\left(R(\xi,q_{l},q_{r})\right)=-\int\displaylimits_{-\infty}^{0}S\left(R(\xi,q_{r},q_{l})\right)$, to get

[TABLE]

(recall $S$ is a convex function of the rotation-invariants of the state-vector $q$), and we finally obtain (49) with $\sum_{j}|\Gamma_{ij}|\boldsymbol{n}_{i\rightarrow j}={\boldsymbol{0}}$. ∎

The bound on $q_{h}(t^{n+1})$ provides one with more stability101010 The convex domain for $q$ is indeed preserved as a consequence of Prop. 3.1, insofar as (48) is more stringent than (43). than Prop. 3.1. In particular, using (49), Prop. 3.2 provides one with an a priori error estimate for FV approximations of smooth solutions to the conservation laws (38), see e.g. [13]. And backward integration (40) of the source term next provides one with a fully admissible approximation $q_{h}(t^{n+1})$:

Proposition 3.3**.**

For any $\tau^{n}>0$, using (40) with $q_{h}^{n+1,-}$ from Prop. 3.2 satisfying (49) yields $q_{h}^{n+1}$ satisfying the following discrete version of (46):

[TABLE]

Proof.

To (49), add (40) tested against $\boldsymbol{\nabla}_{q}S(q_{i}^{n+1})$ for all cells $V_{i}$ (which is possible since $q_{i}^{n+1}$ is admissible as a convex combination in a convex admissibility domain): (52) results from the convexity of $S$ and the definition of $D$. ∎

For SVM neither Prop. 3.2 nor the simpler consequence (45) of Prop. 3.1 can be straightforwardly used because we are not aware of conservative variables $q$ such that $HE$ is convex on the whole admissible domain, i.e. is a mathematical entropy. But note that $HE$ equals $\tilde{S}=H\tilde{E}-H\log(|\boldsymbol{A}_{h}|A_{cc})$ when $H|\boldsymbol{F}|=1$, and it is convex with respect to $\tilde{q}=(H,H\boldsymbol{F},H\boldsymbol{U},HA_{cc}^{1/4})$, and $\boldsymbol{A}_{h}$ which is simply transported. So the standard FV approximation procedure above can be used on slightly modifying the time-splitting: segregating the time-evolution of $\tilde{q}$ and of $\boldsymbol{A}_{h}$ in the first split step. We propose such a discretization of SVM in the sequel that also takes advantage of the possibility to rewrite the SVM system in Lagrangian coordinates to strike a balance between entropy stability, and numerical accuracy (especially for the contact discontinuities satisfied by $\boldsymbol{A}_{h},A_{cc}$). Following [6, 7], we construct entropy-consistent Riemann solvers for (the $\tilde{q}$ sub-system of) a relaxed SVM system in Eulerian coordinates obtained with the help of a BGK approximation in Lagrangian coordinates. The SVM system in Lagrangian coordinates is studied in Appendix A. In the next Sec. 3.2, we use the results of Appendix A to discretize SVM in Eulerian coordinates.

3.2 Finite-Volume approach

to Saint-Venant-Maxwell

We adapt the framework of Section 3.1 to compute discrete FV fields

[TABLE]

that solve SVM on a space-tesselation (i.e. in Eulerian description, with cells paving the same space at all time steps). It splits time-integration into 2 steps.

First step: Given a FV approximation $q_{i}^{n}$ of $q$ at time $t^{n}$, $n\in\mathbb{N}$, we consider the homogeneous SVM system without source term. We use a transport-projection method [8] based on an approximation of (Eulerian) SVM:

[TABLE]

where $\tilde{H}\approx|\boldsymbol{F}_{h}|^{-1}$, $\tilde{\boldsymbol{U}}\approx\boldsymbol{U}$, $\tilde{\boldsymbol{F}}_{h}\approx\boldsymbol{F}_{h}$ still have to be defined such that not only the involution $\tilde{H}|\boldsymbol{F}_{h}|=1$ of SVM is (approximately) preserved, but also

[TABLE]

Indeed, using $dx_{i}=\tilde{U}^{i}dt+\tilde{F}^{i}_{\alpha}da_{\alpha}$ and equality (54) in the smooth case, one can retrieve from (53) the SVM equations in Lagrangian coordinates $(a_{\alpha})$. Then, the entropy-stability of a flux-splitting FV scheme for SVM in Lagrangian coordinates (see Lemma A.2 in Appendix A.1) can be transferred to a simple Riemann solver in Eulerian coordinates, which also captures well contact discontinuities [6]. Precisely, for the Lagrangian-to-Eulerian mapping, we complement (53) by

[TABLE]

and we now propose to use

[TABLE]

where $E^{j}_{\alpha}\approx HF^{j}_{\alpha}$ is the $(\alpha,j)$ entry of the cofactor matrix of $\boldsymbol{F}_{h}^{-1}$. Then, (53), (55) can be closed using for $E^{j}_{\alpha}$ (see e.g. [65])

[TABLE]

which preserves (54). It is noteworthy that preserving (54) ”discretely” with 1D Riemann solvers using (56) reduces to capturing a contact discontinuity.

Now, in this first time-integration step of SVM by a splitting approach, we can update the FV approximation of $\tilde{q}$ by (39) on the one hand, using $F(q_{i}^{n},q_{j}^{n};\boldsymbol{n}_{i\rightarrow j})$ computed from a fully-admissible 1D Riemann solver for the system (53), (55), (56) with a flux that is discontinuous through contact waves for $\boldsymbol{A}_{h}$. On the other hand, $\boldsymbol{A}_{h}$ can be transported and projected in a segregated sub-step preserving full-admissibility of the solution ($\boldsymbol{A}_{h}$ remains SDP with the upwind scheme, and $E$ decays as a single convex function in $\boldsymbol{A}_{h}$).

Precisely, we compute $F(q_{i}^{n},q_{j}^{n};\boldsymbol{n}_{i\rightarrow j})$ with the following simple 1D Riemann solver in direction $\boldsymbol{e}_{m}=\boldsymbol{n}_{i\rightarrow j}$ 111111In (57), $m$ is fixed: Einstein convention is not used. for (53), (55), (56) (motivated by a flux-splitting for SVM in Lagrangian coordinates, see App. A.1):

[TABLE]

where $\delta\in\{\parallel,\perp\}$ denotes the two components in a Cartesian basis $(\boldsymbol{n}^{\parallel},\boldsymbol{n}^{\perp})$ for the geometric coordinates (Eulerian description), and $(\boldsymbol{e}_{e},\boldsymbol{e}_{f})$ is a basis for the material coordinates (Lagrangian description) yet to be precised.

Our motivation for (57) is the possibility to reformulate it as a mapping to Eulerian coordinates of the fully-admissible 1D Riemann solver constructed in App. A.1 for SVM in Lagrangian coordinates:

Proposition 3.4**.**

If $\boldsymbol{e}_{m}=\boldsymbol{n}^{\parallel}$, $F^{\parallel}_{f}\equiv 0$, $E^{m}_{e}F^{\perp}_{f}\equiv 1$, then (57) also writes

[TABLE]

for any initial condition such that $\mathcal{Z}_{ee}:=\Pi^{\parallel}_{e}F^{\perp}_{f}-\Pi^{\perp}_{e}F^{\parallel}_{f}=\Pi^{\parallel}_{e}(E^{m}_{e})^{-1}$.

Proof.

One obtains (58) from (57) by direct computation and for the equivalence it suffices to see that, under assumptions of Prop. 3.4, it holds:

[TABLE]

∎

In particular, the Lagrangian eigenstructure is preserved by the mapping: (58) shows that the system has 3 linearly degenerate waves with speed

[TABLE]

that are ordered $\lambda_{-}\leq\lambda_{0}\leq\lambda_{+}$ if $E^{m}_{e}\geq 0$. Moreover, the 3-wave solutions to (57) that are initialized at $t^{n}$ with left/right values in $V_{i}$/$V_{j}$ such that

[TABLE]

are (formally) consistent with 1D SVM solutions provided

• i)

choosing $\boldsymbol{e}_{f}$ such that $F^{\parallel}_{f}\equiv 0$ is consistent (i.e. $F^{\parallel}_{f}\approx 0$ in reality), and

• ii)

the following evolution equations hold (in the smooth case for some $c^{2}>0$)

[TABLE]

recalling the analysis in App. A.1 for SVM in Lagrangian coordinates (compare with the equation (78) satisfied by the true $\mathcal{P}^{\delta}_{e}$ in Lagrangian description).

Of course, in practice, given any two neighbour cells $V_{i}$/$V_{j}$ and $\boldsymbol{e}_{m}=\boldsymbol{n}_{i\rightarrow j}=\boldsymbol{n}^{\parallel}$, there is no reason why there should exists one direction $\boldsymbol{e}_{f}$ such that $F^{\parallel}_{f}\equiv 0$ in $V_{i}$/$V_{j}$. However, recall that the meaning of the tensor conservative variable $H\boldsymbol{F}_{h}$ in SVM is not purely physical. For the Riemann solver, one can therefore use a reconstruction $F^{\delta}_{f}\neq H^{-1}(HF^{i}_{\alpha})n^{\delta}_{i}f_{\alpha}$ where $n^{\delta}_{i}$, $f_{\alpha}$ are respectively the coordinates of $\boldsymbol{n}^{\delta}$ and $\boldsymbol{e}_{f}$ in $(\boldsymbol{e}_{x},\boldsymbol{e}_{y})$ and $(\boldsymbol{e}_{a},\boldsymbol{e}_{b})$ provided $F^{\delta}_{f}$ retains the physics: $F^{\delta}_{f}$ are the current coordinates in $(\boldsymbol{e}_{x},\boldsymbol{e}_{y})$ (Eulerian description) of a material vector $\boldsymbol{e}_{f}$ attached to the reference configuration (useful in Lagrangian description). And for similar reasons, one can also “project” $H\boldsymbol{F}_{h}$ at the end of the transport-projection method used in the present first step of our splitting scheme to make cell-values compatible with mass and energy conservation.

Thus, at the beginning of the transport-projection method, we reconstruct $F^{\delta}_{e,f}$ at each interface such that the mass is conserved on each side of the interface, such that the elastic energy lost by enforcing $F^{\parallel}_{f}=0$ is minimal, $F^{\parallel}_{e}\geq 0$ and the error due to enforcing $F^{\parallel}_{f}=0$ is small. For consistency and stability reasons, the initial values $c_{l/r}$ should also be well-chosen in the Riemann solver (57) see our analysis for SVM in Lagrangian coordinates in App. A.1.

Proposition 3.5**.**

*Given $\boldsymbol{e}_{m}=\boldsymbol{n}_{i\rightarrow j}=\boldsymbol{n}^{\parallel}$ at some interface in between two cells $V_{i}/V_{j}$, define $(\boldsymbol{e}_{e},\boldsymbol{e}_{f})$ such that $\boldsymbol{e}_{f}$ is the arithmetic mean of two right eigenvectors of $\boldsymbol{F}|_{V_{i}},\boldsymbol{F}|_{V_{j}}$ with singular value closest to 1 and $F^{\parallel}_{e}\geq 0$.

Next, reconstruct $F^{\parallel}_{e}=(\lambda H)^{-1}$, $(F^{\parallel}_{f},F^{\perp}_{f})=(0,\lambda)$, $F^{\perp}_{e}$ in $V_{i}/V_{j}$ such that*

[TABLE]

where $A_{ee},A_{ef},A_{ff}$ are the coefficients of the tensor $\boldsymbol{A}_{h}$ in the basis $(\boldsymbol{e}_{e},\boldsymbol{e}_{f})$, and $\lambda\geq 0$ is chosen close to 1 in $\left|\lambda^{2}-\frac{A_{ee}\mathcal{E}}{2|\boldsymbol{A}|}\right|\leq\frac{A_{ee}\sqrt{\mathcal{E}^{2}-4|\boldsymbol{A}|/H^{2}}}{2|\boldsymbol{A}|}$, hence

[TABLE]

is a real solution to (61). (We choose $x-\log x-1$ as distance to 1 for $x>0$.)

Then, the 1D Riemann solver (57) initialized with (60) and Lemma A.3 for $c$ is fully-admissible in the sense of Prop. 3.2, for $\tilde{S}$ and any direction $\boldsymbol{n}=\boldsymbol{n}_{i\rightarrow j}$ with flux $\tilde{G}_{\boldsymbol{n}}=(H\tilde{S}\mathcal{V}_{e}+\Pi^{\delta}_{e}U^{\delta})|_{\Gamma_{ij}}$ provided $H=|\boldsymbol{F}_{h}|^{-1}$ in $V_{i}/V_{j}$.

Proof.

First, note that the reconstruction using $(\boldsymbol{e}_{e},\boldsymbol{e}_{f})$ always exists. In particular, $\mathcal{E}^{2}\geq 4|\boldsymbol{A}|/H^{2}$ holds when $H=|\boldsymbol{F}_{h}|^{-1}$. Indeed, the inequality rewrites

[TABLE]

which obviously holds for any symmetric positive matrix $\boldsymbol{F}_{h}^{T}\boldsymbol{A}\boldsymbol{F}_{h}$.

Second, note the solution has exactly the same structure with same intermediate states as in Lagrangian coordinates. Then, admissibilty $H>0$ as well as the discrete entropy inequality (thus, full-admissibility) follow from the full-admissibility of the Riemann solver in Lagrangian coordinates, provided it is well-initialized for the equivalence to hold. ∎

To ensure $H=|\boldsymbol{F}_{h}|^{-1}$ in each cell $V_{i}$, we propose the following “projection” at the end of the transport-projection method. We modify the singular values $\lambda_{a^{\prime}},\lambda_{b^{\prime}}$ in the SVD decomposition of $\boldsymbol{F}_{h}=U^{T}\mathop{\rm diag}(\lambda)V$ such that in each cell, the mass $(\lambda_{a^{\prime}}\lambda_{b^{\prime}})^{-1}=|\boldsymbol{F}_{h}|^{-1}=H$ is conserved and the loss $|\lambda_{a^{\prime}}V_{a^{\prime}\alpha}A_{\alpha\beta}V_{a^{\prime}\beta}+\lambda_{b^{\prime}}V_{b^{\prime}\alpha}A_{\alpha\beta}V_{b^{\prime}\beta}-\operatorname{tr}(\boldsymbol{F}_{h}\boldsymbol{A}\boldsymbol{F}^{T}_{h})|$ of elastic energy is minimized. When there exist more than one solution, we choose the one closest to $\boldsymbol{F}_{h}$ before projection in “energy norm” $\operatorname{tr}(\boldsymbol{F}_{h}\boldsymbol{A}\boldsymbol{F}_{h}^{T})$. Of course, if the transport-projection method for the 2D FV discretization preserved $H=|\boldsymbol{F}_{h}|^{-1}$ we would not need that additional step. But this preservation is a well-known difficulty in the discretization of hyperbolic systems with involutions.

In the other sub-step $\boldsymbol{A}_{h}$ can be updated as a solution to transport equations with a scheme preserving the convex domain of $A_{aa}>0,A_{bb}>0,A_{ab}/\sqrt{A_{aa}A_{bb}}\in(-1,1)$. As a matter of fact, if we use an upwind scheme with $\boldsymbol{U}$ given at interfaces by the Riemann problems of the first sub-step, which is consistent with the decrease of free-energy in each cell thanks to the convexity of $\tilde{S}$ with respect to $A_{aa},A_{bb},A_{ab}/\sqrt{A_{aa}A_{bb}}$, then the second sub-step can be done at the same time as the first one.

Second step: Source terms can be integrated as usual with a backward formula in a second time-splitting step.

4 Conclusion

To conclude, we show some numerical simulations illustrating our new model after FV discretization, we discuss the results and we list some perspectives.

Test case 1

First, recall that our model contains the 1D model SVUCM that was derived in [9]: SVUCM is a closed subsystem of our new model SVM, that admits solutions preserving the translation-invariance of initial conditions. Moreover, our 1D Riemann solver is also similar to the one constructed in [9]. And our face and cell reconstructions of the deformation gradient variable $\boldsymbol{F}_{h}$ a priori preserve the translation-invariance of an initial condition on a 1D-conforming grid: $HF^{x}_{a}=1=F^{y}_{b}$, $F^{x}_{b}=0=F^{y}_{a}$ is preserved by time-evolution on a Cartesian grid $\boldsymbol{e}_{x},\boldsymbol{e}_{y}$.

But our FV discretization of SVM does not use exactly the same variables as the one for SVUCM in [9]: $HA_{aa}$, $HA_{cc}$ is used here, while $HB^{xx}=HF^{x}_{a}A_{aa}F^{x}_{a}\equiv A_{aa}H^{-1}$, $HB^{zz}=HF^{z}_{c}A_{cc}F^{z}_{c}\equiv A_{cc}H^{3}$ was used in [9].

Now, the 2D numerical results obtained here in the translation-invariant case $HF^{x}_{a}=1=F^{y}_{b}$, $F^{x}_{b}=0=F^{y}_{a}$, $A_{ab}=0$, $A_{bb}=1$ with a $2^{7}\times 2^{7}=128\times 128$ Cartesian grid for $(x,y)\in[0,8]^{2}$, $t\in[0,.2)$ and an initial condition

[TABLE]

at rest $U^{x}=0=U^{y}$, $A_{aa}=H^{2}=A_{cc}^{-1}$ with $G=1$, $\lambda=.1$, $g=10$ compare well with the 1D results in [9, section 5.5] (Test Case 1). Without source term, the 1D solution consists exactly in a left-going rarefaction wave, a right-going shock wave, and a contact-discontinuity wave [12]. Here, $\lambda$ is quite large in comparison with $T=.2$ and this is as well-captured in Fig. 1 as in [9]. Note however that the latter translation-invariant 2D solution is by no way the unique, and we have indeed observed that other solutions could be captured depending on the cell reconstruction at the end of the transport-projection method.

Test case 2

Second, we now run the previous test case with the initial condition rotated by $\pi/4$ on the same 2D grid (see view in Fig. 3).

Although the results in Fig. 2 compare with the previous simulations, one now sees that translation invariance is broken, see Fig. 3. We believe it is mainly because of our mass-conforming projection of $\boldsymbol{F}_{h}$ in each cell at the end of the transport-projection method, whose implementation through re-balancing singular values in SVD decomposition does not preserve well the symmetry. One consequence is that the contact discontinuity is smeared.

Test case 3

We next simulate an axisymetric test case (with rotation invariance): the collapse of a 2D column initially at rest. Initially, we choose:

[TABLE]

in polar coordinates $(\boldsymbol{e}_{r},\boldsymbol{e}_{\theta})$, with $\boldsymbol{F}=H^{-1}\boldsymbol{e}_{r}\otimes\boldsymbol{e}_{r}+\boldsymbol{e}_{\theta}\otimes\boldsymbol{e}_{\theta}$.

The results of Fig. 4 are consistent with an axisymetric (rotation-invariant) solution and compares with the unidirectional (translation-invariant) solutions of cases 1 and 2. Note however that the discrete solution loses some symmetry ($\pi/2$-rotation) and the contact discontinuity wiggles, see Fig. 5.

Test case 4

Last, we simulate a standard permanently-sheared flow on a long time range $t\in[0,10]$, expecting a steady viscous behaviour asymptotically. Precisely, we consider a lid-driven cavity on the same 2D grid as in cases 1,2 and 3, with a uniform initial condition $\boldsymbol{F}_{h}=\boldsymbol{I}$ at rest, $G=1$, $g=10$ and $\lambda=.1$. The following conditions are prescribed at boundaries121212 We have not precisely studied the 2D Initial-Boundary Value Problem and we impose steady values at boundary cells for the sake of simplicity. Note: boundary conditions were not meaningful in cases 1,2 and 3 so long as the front is far enough from the boundary. :

[TABLE]

with $\boldsymbol{A}$ at equilibrium (so relaxation source term is zero).

We recall the physical interpretation for the coefficients $F^{i}_{\alpha}$. At time $t$, a basis of “material vectors” $\boldsymbol{e}_{a},\boldsymbol{e}_{b}$ in the reference configuration (Lagrangian description) has been stretched and turned into two “geometrical vectors” $\boldsymbol{F}_{h}\boldsymbol{e}_{a},\boldsymbol{F}_{h}\boldsymbol{e}_{b}$ resp. of the current configuration (Eulerian description). So our condition on the “left” boundary (the cavity lid) actually means that $\boldsymbol{e}_{a}$ (supposedly aligned with $\boldsymbol{e}_{x}$ in the reference configuration) is permanently stretched and turned by shearing.

Clearly, standard inviscid Saint-Venant equations would not sustain a developped vortex in such conditions, while the standard viscous Saint-Venant would immediately develop a developped vortex in the whole cavity after start. On the contrary, our model can capture the transient development of a vortex until reaching a stationary state, see Fig. 6.

With a time-step that remains approximately constant equal to $.007$ (under our $1/2$ CFL condition, recall section 3.1), our FV discretization quickly reaches a nearly stationary state with $H$ (thus, pressure) almost uniform in space except at the boundary, close to which one can observe the stagnation point of a (nearly stationary) vortex, see Fig. 7.

Last, note that in each case we have run our simulations on finer and finer grids: they seem to converge globally in space, despite symmetry losses. However, we are aware that multi-dimensional conservation laws can admit many entropy weak solutions, and a precise study remains to be done.

Finally, let us summarize the main features of our new symmetric hyperbolic system of conservation laws to model viscoelastic flows of (Upper-Convected) Maxwell fluids with shear-waves propagating at finite-speed.

For the 2D gravity flows with a free surface detailed herein, the resulting generalization of Saint-Venant system has not only a theoretically interesting formulation. But moreover, physically-reasonable numerical simulations have shown that our model indeed satisfies our initial goal. This is encouraging from the conceptual viewpoint of physical principles as well as for practical applications e.g. to transient geophysical flows.

However, to fully validate our new model and make it as useful as possible, a number of points need to be tackled.

• •

Although it is satisfying that real (viscous) fluids seem well-modeled by symmetric-hyperbolic conservation laws, it is also well-known that the mathematical theory is far from complete at present. There remain challenging issues to precisely define physically-meaningful solutions at large times, and next simulate them numerically.

Our FV approach certainly needs improving to that aim.

• •

The large-time asymptotic stability of sheared-flows needs to be investigated theoretically and numerically, as boundary-value problems.

• •

In many practical applications, one would also want to use our new model with 3D flows ; adequate choices of pressures should therefore be more carefully studied. Moreover, specific applications would also require one to modify the deviatoric stresses. It remains to investigate the various non-Newtonian possibilities in detail, with care as regards “non-isothermal” flows (recall our microscopic interpretation of the new state variable $\boldsymbol{A}$ as a mean-field approximation of material distortion thermally-agitated).

Appendix A

Riemann solvers based on flux-splitting

The entropy-consistency condition (47) can be achieved with simple 1D Riemann solvers which yield flux-splitting FV schemes with a kinetic interpretation [6]. Here, we consider 1D Riemann solutions of the Lagrangian system (29):

[TABLE]

that are consistent in the sense of (47) with the following entropy inequality:

[TABLE]

when the flux-coefficients $A_{\alpha a}$ are given so the free-energy functional:

[TABLE]

is actually a (strictly convex) mathematical entropy for (62).

A.1 Entropy-consistent solver in Lagrangian coordinates

The entropy-consistent simple 1D Riemann solvers, satisfying (47) with $\boldsymbol{n}=\boldsymbol{e}_{a}$, can be characterized in the flux-splitting case [6, p.643]. Recall: $E$ is a strictly convex entropy for a system $\partial_{t}u+\partial_{a}F(u)=0$ like (62), is equivalent to, the entropy variable $v\equiv\partial_{u}E$ symmetrizes the system and $F=\partial_{v}\psi$. Then [6]:

Lemma A.1**.**

The FV schemes built with flux-splitting $F=F^{+}+F^{-}$ for $\partial_{t}u+\partial_{a}F(u)=0$ is entropy-consistent if $F^{+}=\partial_{v}\psi^{+}$, $F^{-}=\partial_{v}\psi^{-}$ are built from two scalar functions $\psi^{+},-\psi^{-}$ convex in $v$.

Moreover, the entropy-consistency (47) can be fully analyzed after kinetic interpretation like in [6]. Then, considering $A_{cc}$ as flux-coefficient like $A_{\alpha\beta}$ for $u=(H^{-1},F^{x}_{a},F^{y}_{a},U^{x},U^{y})$ solution to (62) here, let us split $F=\partial_{v}\psi$ where

[TABLE]

Let us choose the following ”convex” splitting:

[TABLE]

[TABLE]

with a single parameter $c>0$ (to be fixed later, e.g. for entropy-consistency). One can check that the entropy fluxes $G^{\epsilon}=v\cdot F^{\epsilon}-\psi^{\epsilon}$ associated with $F^{\pm}$, such that $\partial_{u}G^{\epsilon}=\partial_{u}E\cdot\partial_{u}F^{\epsilon}$, recall $\partial_{u}\psi^{\epsilon}=(\partial_{uu}^{2}E)\partial_{v}\psi^{\epsilon}$, do “dissipate” [6], i.e.

[TABLE]

holds true for all $u_{1},u_{2}$. (Note indeed that (65) is equivalent to the convexity of $\epsilon\psi^{\epsilon}$, it can be checked by a direct computation here on noting $G^{\epsilon}=\psi^{\epsilon}$ and $\epsilon G^{\epsilon}=\epsilon\psi^{\epsilon}$ is a convex function of $v=\partial_{u}E$ such that $F^{\epsilon}=\partial_{v}\psi^{\epsilon}$). Then, following [6], to fully check the entropy-consistency condition (47) for the flux-splitting above as numerical flux – indeed a simple Riemann solver –, one can use its kinetic interpretation as a discretized BGK model with relaxation form131313Where $\varepsilon$ assumes its usual meaning for relaxation systems [6]: some time $\varepsilon\rightarrow 0$ characteristic of the (numerical) relaxation and infinitesimal in comparison with the time-step of the numerical scheme. It should not be mixed with $\epsilon\in\{+,-\}$.:

[TABLE]

also endowed with a diagonal $8\times 8$ formulation:

[TABLE]

that clearly shows the system (66) has only linearly degenerate fields.

In particular, one can add $\partial_{t}A_{cc}=0$ without changing the structure of the Lagrangian system (66), as well as

[TABLE]

to check (47) with $\tilde{G}_{\boldsymbol{e}_{b}}(u_{l},u_{r})=(\Pi^{i}_{a}U^{i})|_{a/t=0}$ as in [6]. It amounts to add

[TABLE]

in (67). Then, note that for (47) to hold, it is sufficient that

[TABLE]

holds for all states $u$ in Riemann solutions of the extended system (67)–(68) i.e.

[TABLE]

for all intermediate state $u_{1}$ in the Riemann solution such that $(\hat{e}-|\Pi^{x}_{a}|^{2}/2c^{2}-|\Pi^{y}_{a}|^{2}/2c^{2})\equiv E(u_{2})-\frac{G^{+}(u_{2})-G^{-}(u_{2})}{c}$ with $u_{2}$ the ”outward” neighbour state (chosen in the direction opposite to the interface, which is directly fixed by the boundary condition here in the 3-wave case).

Now, the split-form of the extended system (67)–(68) allows one to check (69) independently for left and right (intermediate) states $u_{1}$ (with $u_{2}$ resp. the left and right boundary condition here), when $A_{cc}$ is a variable solution to a contact disconitnuity wave, and when the coefficient $A_{\alpha\beta}$ is discontinuous in

[TABLE]

Recalling (65) and $c(u_{2}-u_{1})+F(u_{2})-F(u_{1})=0$ then (69) rewrites

[TABLE]

[TABLE]

for boundary states $u_{2}=u_{l/r}$ resp. with neigbour intermediate state $u_{1}=u^{*}_{l/r}$

[TABLE]

which can be checked e.g. using

[TABLE]

[TABLE]

When $u_{1}=u_{2}$, (70) is satisfied, and for $u_{2}=(-\mathcal{P},GA_{aa}F^{x}_{a},GA_{aa}F^{y}_{a},)_{2}\approx u_{1}$

[TABLE]

so on recalling $\partial_{H^{-1}}\mathcal{P}=-(gH^{3}+3GH^{4}A_{cc})$, a sufficient condition for entropy-consistency (of the 3-wave Riemann solver in Lagrange coordinates) reads:

[TABLE]

on an admissible neighbourhood of $u_{1}$ containing $u_{2}$ (such that $H>0$).

Lemma A.2**.**

The 1D Riemann solver (67)–(68) for (the $\tilde{q}$ subsystem of) SVM in Lagrangian coordinates is entropy-consistent in the sense of Prop. 3.2, for the entropy $E$ given by (64) with numerical flux $\tilde{G}_{\boldsymbol{e}_{b}}=(\Pi_{a}^{i}U^{i})|_{a/t=0}$, provided (75) holds on admissible neighbourhoods of the left/right states that resp. contain the left/right intermediate state.

Then, one can look for numerical values of $c$ satisfying (75), e.g. with $c$ solution to $\partial_{t}c=0$ as in [7] to enforce $H_{l/r}^{*}>0$, and it remains to map the solver above into Eulerian coordinates (see Section 3.2). Note however a consistency limitation on the choice of $\boldsymbol{e}_{a}$ to approximate SVM Riemann solutions in material coordinates as above, with 3 waves only. Indeed, the 1D solutions computed with (66) are translation-invariant solutions to a 2D hyperbolic system of conservation laws in Lagrangian coordinates like for instance:

[TABLE]

where $\alpha,\beta\in\{a,b\}$ refers to the axes of one particular Cartesian coordinate system. This relaxation generalizes to

[TABLE]

with the same symmetric-positive coefficient-matrices $\tilde{c}^{2}_{\alpha,\gamma}=c^{2}_{\alpha,\gamma}=c^{2}_{\gamma,\alpha}$ for $(\Pi^{i}_{\alpha})_{\alpha}$, $(Z_{\beta\alpha})_{\alpha}$ in the wave equations resp. for $U^{i}$, $\mathcal{V}_{\beta}$. (The 1D Riemann solver (66) coincides with 1D solutions of (77) along one principal directions). But it seems impossible to find a 2D system that admits the 1D Riemann solver (66) as particular solution, and that is a relaxation of the 2D Lagrangian SVM system (29) for all (smooth) solutions. Indeed, recalling $\partial_{H^{-1}}\mathcal{P}=-(gH^{3}+3GH^{4}A_{cc})$, $\Pi^{i}_{\alpha}$ in (77) above cannot consistently approximate any exact solution $\mathcal{P}^{i}_{\alpha}$ to

[TABLE]

insofar as the term $\partial_{\alpha}U^{j}$ in the last line of (78) is always missing (it vanishes though for 1D solutions along some direction $\beta$ such that $\sigma_{ij}F^{j}_{\beta}=0$ for all $i$). So the 1D Riemann solver above with 3-wave cannot be consistent with any solution of the 2D SVM system, even smooth.

A.2 Parameter initialization

To achieve (75) in a Riemann problem, one can require left/right initial values of the ”relaxation parameter” $c$ as follows

Lemma A.3**.**

Defining for $\delta\in\{\parallel,\perp\}$ and $o=\{l,r\}$ in a Riemann problem

[TABLE]

then conditions (75) are ensured using as initial values (with $o\neq o^{\prime}\{l,r\}$):

[TABLE]

Proof.

The computation is straightforward and similar to Lemma 5.3 in [9]. ∎

Appendix B Hyperbolicity of SVUCM

We investigate the hyperbolicity of SVUCM and variations with a non-zero slip-parameter, the so-called Johnson-Segalman models see e.g. [12, 11].

Proposition B.1**.**

Among all Gordon-Schowalter derivatives with slip-parameter $\zeta\in[0,2]$, only $\zeta=0$ (i.e. the Upper-Convected case) ensures hyperbolicity of (1–2–9–10) under strain-free constraints, namely: $H,B_{zz}>0$ and $\boldsymbol{B}_{h}=\boldsymbol{B}_{h}^{T}>0$.

The proof follows from rotation-invariance, after computing the eigenvalues of the jacobian in a 1D projection of the system (1–2–9–10) like

[TABLE]

similarly to the proof in [24] for a similar system written in stress variables $\boldsymbol{\Sigma}$ when $\zeta=0$ (though without vertical stress and strain components, which allow here mass preservation). Denoting $\Delta=2gh+G\left(2(3-2\zeta)B_{zz}+\zeta B_{yy}-3\zeta B_{xx}\right)=2gh+6GB_{zz}+G\zeta\left(B_{yy}-4B_{zz}-3B_{xx}\right)$, four eigenvalues read

[TABLE]

and are real if, and only if, the following strain-parametrized inequality holds

[TABLE]

Unless $\zeta=0$, the quadratic condition (82) on $\zeta$ is not clearly satisfied for all values $B_{xx},B_{yy},B_{xy},B_{zz}$ of the strain. We therefore consider only (1–2–9–10) when $\zeta=0$ (the SVUCM case), where hyperbolicity is ensured with eigenvalues $u\pm\sqrt{gh+3GB_{zz}+GB_{xx}}$, $u\pm\sqrt{GB_{xx}}$ and $u$ (with multiplicity 3) under the physcially-natural constraints $h\geq 0,B_{zz}\geq 0,B_{xx}\geq 0$.