B Després

E Labourasse
email: emmanuel.labourasse@cea.fr

Angular Momentum preserving cell-centered Lagrangian and Eulerian schemes on arbitrary grids

Keywords: Compressible fluid dynamics, cell-centered Lagrangian and Eulerian schemes, general grids, angular momentum conservation, conservation laws

We address the conservation of angular momentum for cell-centered discretization of compressible fluid dynamics on general grids. We concentrate on the Lagrangian step which is also sufficient for Eulerian discretization using Lagrange+Remap. Starting from the conservative equation of the angular momentum, we show that a standard Riemann solver (a nodal one in our case) can easily be extended to update the new variable. This new variable allows to reconstruct all solid displacements in a cell, and is analogous to a partial Discontinuous Galerkin (DG) discretization. We detail the coupling with a secondorder Muscl extension. All numerical tests show the important enhancement of accuracy for rotation problems, and the reduction of mesh imprint for implosion problems. The generalization to axi-symmetric case is detailed.

Introduction

This work intends to contribute to a long lasting CFD debate which is the enhancement of the accuracy of compressible fluid solvers for vortical flow. In our case we concentrate more on cell-centered Lagrangian compressible schemes on moving grids. But as demonstrated by the numerical results, the proposed approach is also valid for Eulerian calculations on a fixed grid using a Lagrange+Remap procedure.

A seminal and inspirational work is the one of Dukowicz and Meltz [START_REF] Dukowicz | Vorticity errors in multidimensional lagrangian codes Original[END_REF] where the authors analyze the spurious vorticity errors of the Lagrangian Caveat scheme [START_REF] Adessio | Caveat: a computer code for fluid dynamics problems with large distortion and internal slip[END_REF] and develop a procedure to correct these errors. It is in our mind representative of situations where vorticity of the numerical flow is seen as a potential source of problems, that must be controlled. Such kind of procedure has also been developped in [START_REF] Caramana | Curl-q": A vorticity damping artificial viscosity for essentially irrotational Lagrangian hydrodynamics calculations[END_REF] for staggered schemes (the curl-Q pseudo-viscosity). A general review of vorticity in Finite Volume schemes is in [START_REF] Roe | Preserving vorticity in finite-volume schemes, Finite volumes for complex applications II[END_REF]. See also [START_REF] Toth | Divergence-and curl-preserving prolongation and restriction formulas[END_REF].

In our case we consider that the situation of cell-centered Lagrangian schemes has somewhat changed since the Dukowicz-Meltz contribution. Cell-centered Lagrangian are now becoming a mature ensemble of techniques, owing to the preservation of the GCL (geometric conservation law) and the compatibility with the entropy principle which is rendered possible by the use of nodal-based Riemann solvers instead of the standard edge-based solvers used in [START_REF] Dukowicz | Vorticity errors in multidimensional lagrangian codes Original[END_REF]. It started in [START_REF] Després | Symmetrization of Lagrangian gas dynamics and Lagrangian solvers[END_REF][START_REF] Després | Lagrangian gas dynamics in 2D and lagrangian systems[END_REF], and was developed in [START_REF] Maire | A cell-centered Lagrangian scheme for 2D compressible flow problems[END_REF]. The numerical formalism of the general multi-D Glace scheme is developed in [START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF]. The Eucclhyd formalism on 2D grids was later proposed in [START_REF] Maire | A unified sub-cell force-based discretization for cell-centered Lagrangian hydrodynamics on polygonal grids[END_REF]. The difference between the nodal Riemann solvers of Glace and Eucclhyd is quite small, similar in a sense to different quadrature formulas in the theory of Finite Element Methods for elliptic equations. We quote [START_REF] Cheng | Positivity-preserving Lagrangian scheme for multi-material compressible flow[END_REF][START_REF] Boscheri | A direct arbitrary-Lagrangian Eulerian ADER-WENO finite volume scheme on unstructured tetrahedral meshes for conservative and nonconservative hyperbolic systems in 3D[END_REF][START_REF] Toth | Divergence-and curl-preserving prolongation and restriction formulas[END_REF] for recent related works for compressible Lagrangian fluid dynamics. It has also been used to develop artificial viscosities for staggered schemes [START_REF] Burbeau-Augoula | A Node-Centered Artificial Viscosity Method for Two-Dimensional Lagrangian Hydrodynamics Calculations on a Staggered Grid[END_REF][START_REF] Loubère | 3D staggered Lagrangian hydrodynamics scheme with cell-centered Riemann solver-based artificial viscosity[END_REF]. Extension to elastoplastic solvers has been recently performed in a series of papers: it started with [START_REF] Kluth | Discretization of hyperelasticity on unstructured mesh with a cell-centered Lagrangian scheme[END_REF] where the form of the nodal elastic Riemann solver is defined in the context of very general hyper-elastic models, then later extended in [START_REF] Maire | A nominally secondorder cell-centered Lagrangian scheme for simulating elastic-plastic flows on twodimensional unstructured grids[END_REF][START_REF] Burton | A cell-centered Lagrangian Godunov-like method for solid dynamics[END_REF] for simplified hypoelastic models. Definition of a stabilization procedure named subzonal entropy is proposed in [START_REF] Després | Stabilization of cell-centered compressible Lagrangian methods using subzonal entropy[END_REF]. A proof of weak consistency is given in [START_REF] Després | Weak consistency of the cell centered lagrangian glace scheme on general meshes in any dimension[END_REF]. Most of these methods have a wide domain of efficiency in terms of stability and accuracy, the main reason being the compatibility with the GCL and with the entropy principle. Moreover cell-centered Lagrangian schemes are naturally adapted to remapping procedure which means that any stability issue of the mesh can be addressed using ALE (Arbitrary Lagrange Euler), still maintaining the conservation properties and the accuracy for shock calculations, see for example [START_REF] Galera | A two-dimensional unstructured cell-centered multi-material ALE scheme using VOF interface reconstruction[END_REF][START_REF] Pino | Metric-based mesh adaptation for 2D Lagrangian compressible flows[END_REF][START_REF] Boscheri | A direct arbitrary-Lagrangian Eulerian ADER-WENO finite volume scheme on unstructured tetrahedral meshes for conservative and nonconservative hyperbolic systems in 3D[END_REF]. All these Lagrangian schemes can be run in a purely Eulerian mode by using Remap at every time steps. Based on these advances we consider that the time is less to consider that vorticity is a spuriosity that must be controlled or eliminated (another drawback being of course that physical vorticity might be treated like spurious vorticity), but more to enhance the accuracy of flow with strong vortical parts.

Our analysis starts from a well known physical principle which is that angular momentum

w = u ∧ x (1)
is solution of a conservation law

∂ t (ρw) + ∇ • (ρu ⊗ u ∧ x) + curl(px) = 0. (2)
At the analytical level this law is redundant with the inertial momentum equation

∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = 0. (3)
But at the numerical level, many basic simulations show without any doubt that angular momentum is far to be preserved by standard cell-centered flow solvers on general grids (situation is less severe on Cartesian fixed grids). This is why we concentrate in this work on the development and analysis of a general numerical method for the preservation of the angular momentum variable [START_REF] Adessio | Caveat: a computer code for fluid dynamics problems with large distortion and internal slip[END_REF]. We will show that angular momentum discretization can be understood as a special DG discretization of (3), which seems to us a new result with respect to the literature [START_REF] Loubère | Lagrangian discontinuous Galerkin-type method on unstructured meshes to solve hydrodynamics problems[END_REF][START_REF] Shu | A brief survey on discontinuous Galerkin methods in computational fluid dynamics[END_REF][START_REF] Vilar | Cell-centered discontinuous Galerkin discretization for two-dimensional Lagrangian hydrodynamics[END_REF]. We notice also that the vorticity is easily recovered from the knowledge of angular momentum and of the inertial momentum. Let x 0 be a given point in the domain and consider w 0 (x) = w(x)-u(x)∧x 0 = u(x)∧(xx 0 ). Then ∆w 0 (x) = ∆u(x) ∧ ∇ (xx 0 ) -2∇ ∧ u(x). Therefore one has the identity ∇ ∧ u(x 0 ) = -1 2 ∆w 0 (x) which shows that the vorticity ω = ∇ ∧ u can be computed once one knows the inertial momentum and the angular momentum. Conservation of angular momentum is also an important question for many different fundamental problems in fluid dynamics. We just give few examples. A first problem is fluid simulations of the atmosphere around the earth. Indeed it is known that angular momentum of the atmosphere may interact with the angular momentum of the planet itself, an early work on this topic is to be found in [START_REF] Oort | Angular Momentum in the Atmosphere-Ocean-Solid Earth System[END_REF]. Quoting indeed a recent PhD thesis [START_REF] Driscoll | The Earth's Atmospheric Angular Momentum Budget and its Representation in Reanalysis Observation Datasets and Climate Models[END_REF], the Angular Momentum budget represents a beautiful example of how the atmosphere, oceans and solid earth interact. In this context an accurate computation of angular momentum is necessary to simulate such systems with fluid flow solvers. For this example we are not aware of any use of standard Finite Volume CFD schemes. A completely different physical problem is rotation of MHD flows in Tokamaks for which angular preservation is clearly fundamental issue. It is addressed in the context of MHD solvers, either full MHD or reduced MHD, a general review is to found in [START_REF] Jardin | Review of implicit methods for the magnetohydrodynamic description of magnetically confined plasmas[END_REF]. We notice that Finite Volume techniques are rarely used in the Tokamaks community. On the other hand Godunov solvers are widely used for astrophysical flows, and angular momentum is a key feature for an accurate numerical treatment of the rotation of stars and planets: many works are devoted to this issue on Cartesian fixed grids and we quote only on few of them such as [START_REF] Mignone | A conservative orbital advection scheme for simulations of magnetized shear flows with the PLUTO code[END_REF][START_REF] Skinner | The Athena Astrophysical Magnetohydrodynamics Code in Cylindrical Geometry[END_REF]. It has also a big impact for the chemical reactions into the combustion chamber of engines, in which the intake valve is usually placed to give the mixture a pronounced swirl [START_REF] Amsden | KIVA-2: A computer program for chemically reactive flows with sprays[END_REF]. The initial stage of turbulent flows is also clearly dominated by the strong vortices inside the flow. In the context of this work, we will show that preservation of the angular momentum enhances a lot the accuracy of implosion calculations near the focusing point, and that it minimizes the mesh imprint for such problems. A simple proof will be given that explains this fact. We stop here the list of such examples, but it is clear that vortical flows and related problems challenge the quality of flow solvers on arbitrary grids in many areas of applied science and computational fluid dynamics.

The plan of the works is as follows. The basis of our method is detailed in section 2, where we propose to add a local degree of freedom to respect the preservation of angular momentum. The structure of the new scheme is detailed using the Glace formalism. We also explain how the new scheme can be recast as a special DG method. Next in section 3, we analyze the stability with the entropy principle. Section 4 is devoted to some key implementation details, in particular how to design an angular momentum scheme compatible with the Muscl techniques which are in our case essential to obtain a stable second-order. Implementation of the method in axi-symmetric formulation is addressed in section 5. We then turn to dedicated numerical examples in section 6 and conclude.

Conservation laws and cell-centered discretization

We first briefly discuss the conservation laws associated with Euler equations, and in particular the one related to angular momentum. Then we recall the semi-discrete form characteristic of cell-centered schemes and propose an extension of these schemes to take into account the angular momentum in semi-discrete form.

Integral form of conservation laws

A conventional integral form of the Euler equations is

               d t Ω(t) 1 = Ω(t) div (u) , volume conservation, d t Ω(t) ρu = - Ω(t)
∇p, inertial momentum conservation,

d t Ω(t) ρe = - Ω(t) div (pu) , total energy conservation, (4)
where

Ω(t) ⊂ R d is a bounded open domain in dimension d ∈ 1, 2, 3
. This domain moves with the fluid velocity u, and ρ is the density of the fluid, p is the pressure, e = ε + u 2 /2 is the total energy, and ε is connected to ρ and p by the equations of states (eos). The

d-measure of Ω(t) is |Ω(t)| = Ω(t) 1.
Without loss of generality, the volume equation of system (4), also called GCL for Geometric Conservation Law, can be replaced by the mass conservation equation

d t Ω(t) ρ = 0,
where

Ω(t)
ρ is the mass.

(

) 5
This equation is trivially satisfied by the semi-discrete Lagrangian. Our concern in this work is the conservation of the angular momentum w = u ∧ x. This conservative equation that can be written for this quantity is less commonly mentioned in the literature, probably because it is redundant with the conservation of inertial momentum. Using the inertial momentum equation of the system (4) and the Reynolds theorem yield the following integral formulation of angular momentum equation

d t Ω(t) ρw = - Ω(t) curl (px) . ( 6
)
The equation is here written in integral form. The equivalent Eulerian conservation law is [START_REF] Aguirre | A vertex centred Finite Volume Jameson-Schmidt-Turkel (JST) algorithm for a mixed conservation formulation in solid dynamics[END_REF]. It can be checked that usual cell-centered schemes do not preserve this law. This is particularly evident on the simulations reported in [START_REF] Kucharik | Enhancement of Lagrangian slide lines as a combined force and velocity boundary condition[END_REF][START_REF] Clair | A one-mesh method for the cell-centered discretization of sliding[END_REF].

Semi-discrete cell-centered equations

Before going further, we remind the general form of cell-centered schemes which can be found for instance in the reference [START_REF] Kluth | Discretization of hyperelasticity on unstructured mesh with a cell-centered Lagrangian scheme[END_REF]. This formulation gives a general simple framework for the presentation of both Glace and Eucclhyd schemes. For sake of simplicity, we present the intermediate stages of construction of the scheme into the Glace formalism which will serve as a paradigm for the presentation. Most of the ideas are immediate to generalize to many Lagrangian and Eulerian flow solvers.

We consider a decomposition of the computational domain Ω into control volumes V j (or cells) indexed by j. We assume that the boundaries (or faces E l , indexed by l) of these control volume are entirely defined by a set of vertex (or nodes) which are indexed by r. The numerical simulations in 2D use standard elements with straight edges and two vertices by edge. The decomposition of Ω is assumed to fulfill the following conditions: the full domain is covered and there is no overlap between cells

∪ j V j = Ω, and dim (V i ∩ V j ) < d if i = j.
For a given edge indexed by r, the set of neighboring cells is denoted a C(r). Similarly, the set of nodes that belong to the j-th cell is denoted N (j), and the set of faces that belong to the j-th cell could be denoted as E(j). In addition, we denote x r (resp. u r ) the position (resp. the velocity) of the r-th node within the mesh. These are vectors of dimension d. In the text, we will often refer to nodal position and nodal velocity.

Let's introduce the corner vectors C jr defined as the gradient of the volume |V j | with respect to the nodal positions x r :

C jr = ∇ xr |V j |.
Using this notation, the semi-discretization of the system (4) is the following:

             m j d t τ j (t) = r∈N (j) C jr • u r , m j d t u j (t) = - r∈N (j) C jr p jr , m j d t e j (t) = - r∈N (j) C jr • u r p jr . ( 7
)
where m j , u j and e j are respectively the Lagrangian mass, the mean velocity and the mean specific total energy of the cell j, and d t ϕ ≡ (∂ t + u • ∇) • ϕ is the Lagrangian derivative of a generic function ϕ. The nodal unknowns u r (the velocity of the r-th node), and p jr (the nodal pressure) are computed thanks to a nodal Riemann solver for which many possibility exist nowadays. The one we consider consists in solving the following set of linear equations

   p jr -p j + ρ j c j (u r -u j ) • n jr = 0, j∈C(r) C jr p jr = 0, (8)
where n jr = C jr |C jr | is a normalized direction. In the remaining of the paper, we will implicitly use the notation that j ≡ j∈C(r) and r ≡ r∈N (j) . The nodal solver (8) can be rewritten as

A r u r = j A jr u j + C jr p j , F jr = C jr p j + A jr (u j -u r ) . (9)
In this system, the Glace and Eucclhyd schemes differ ultimately only by the definition of the R d×d matrices A jr and of the vector F jr , which shows this formalism embodies both the Glace and the Eucclhyd schemes, and many similar methods in the same class. We consider in the following the Glace scheme [START_REF] Burton | A cell-centered Lagrangian Godunov-like method for solid dynamics[END_REF] for which

A jr = ρ j c j C jr ⊗ C jr |C jr | , A r = j A jr and F jr = C jr p jr . (10)
Under very reasonable conditions on the mesh, A r is a symmetric positive matrix, and A jr is a symmetric non-negative matrix of rank 1 for Glace. A jr is a a symmetric positive matrix for Eucclhyd.

The semi-discrete system [START_REF] Burbeau-Augoula | A Node-Centered Artificial Viscosity Method for Two-Dimensional Lagrangian Hydrodynamics Calculations on a Staggered Grid[END_REF] recasts

             m j d t τ j = r C jr • u r , m j d t u j = - r F jr , m j d t e j = - r F jr • u r . (11)
Two important properties are the conservativity of the scheme and the correct entropy law. Concerning the conservativity, it is easy to check (see also [START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF]). The increase of the entropy can be shown as follows. One has

         m j T j d t S j = m j d t ε j + p j m j d t τ j , = m j d t (e j -|u j | 2 /2) + p j m j d t τ j , = m j d t (e j -|u j | 2 /2) + p j m j d t τ j , = m j d t e j -m j u j • d t u j + p j m j d t τ j ,
where S j and T j account for the entropy and the temperature of the cell j. Replacing the time derivatives and the spatial derivatives by there respective discrete expressions from [START_REF] Burbeau-Augoula | A Node-Centered Artificial Viscosity Method for Two-Dimensional Lagrangian Hydrodynamics Calculations on a Staggered Grid[END_REF], we obtain

m j T j d t S j = - r F jr • u r + u j • r F jr + p j r C jr • u r , (12)
= r (u r -u j ) • A jr (u r -u j ) ≥ 0, (13)
using the expression of F jr provided in [START_REF] Caramana | Curl-q": A vorticity damping artificial viscosity for essentially irrotational Lagrangian hydrodynamics calculations[END_REF]. Since the matrices A jr are symmetric non negative, the entropy increases. However, another property trivially true at the continuous level is not guaranteed at the discrete level: it is the conservation of angular momentum. Let us consider the discrete mean angular momentum for the cell j: w j = 1/|V j | j x ∧ u. The angular momentum balance is

m j d t w j = - r F jr ∧ x j + m j d t x j ∧ u j , (14)
where x j = 1/|V j | j x denotes the center of mass of the cell. A summation of ( 14) over all cells yields j m j d t w j =j r

F jr ∧ x j + j m j d t x j ∧ u j . ( 15
)
Remark 2.1 (First term in the right hand side of equation ( 15)). Unlike the situation for staggered schemes, the flux p jr in F jr = C jr p jr depends on both the cell and the node. It can not get out of the sum, even by reversing the summation over j and r.

Remark 2.2 (Second term in the right hand side of equation ( 15)). For cell-centered schemes one has a priori that d t x j = u j since u j is just a cell-wise mean value of the velocity.

These two remarks imply that cell-centered schemes are not natively conservative in angular momentum.

Cell-centered discretization of angular momentum

Let us now consider that w j is no more dependent of u j . Then, as observed for instance in [START_REF] Burton | A cell-centered Lagrangian Godunov-like method for solid dynamics[END_REF], a semi-discrete form of (6) can easily be derived

m j d t w j = - r∈N (j) F jr ∧ x r . (16)
The only difference between equations ( 14) and ( 16), is that in [START_REF] Després | Symmetrization of Lagrangian gas dynamics and Lagrangian solvers[END_REF] the vector position of the node x r replaces the vector position of the cell-center x j in the fluxes, and that the term m j d t x j ∧ u j has been omitted. However, the angular momentum balance [START_REF] Després | Stabilization of cell-centered compressible Lagrangian methods using subzonal entropy[END_REF] now recasts into

j m j d t w j = - j r F jr ∧ x r , (17)
= r x r ∧ j F jr , (18)
= 0, (19)
since the Riemann solver enforces j F jr = 0. This way, we have a conservation law on the angular momentum w, but it has no more connexion with the velocity field u.

To tide over this difficulty, we propose to add a degree of freedom to the algorithm, in enriching the velocity field

v j (x) = a j + b j ∧ x. (20)
where a j and b j are constant per cell. That is the enriched velocity field is locally made of all rigid body displacements. Obvioulsy, if b j = 0, ∀j, that is only translations are considered, we recover the usual constant by cell approximation of the first-order Glace and Eucclhyd schemes. Additional interpretation in terms of partial Discontinuous Galerkin extension will be provided in next section. The corresponding mean velocity is defined by

u j := 1 |V j | j a j + b j ∧ x, = a j + b j ∧ x j . (21)
We can also deduce w j

w j := 1 |V j | j (a j + b j ∧ x) ∧ x, = (a j + b j ∧ x j ) ∧ x j + 1 |V j | j (b j ∧ (x -x j )) ∧ (x -x j ). (22)
Equations ( 21) and ( 22) assume implicitly that x j is the center of mass.

In the following we will denote 0 < H j = H t j ∈ R d×d the positive symmetric matrix defined by

∀b ∈ R d , (H j b, b) := 1 |V j | j |b ∧ (x -x j )| 2 . ( 23
)
One has that |H j | = O(h 2 ) where h accounts for the characteristic length of the mesh, and then the second term in the right-hand-side of the previous equation is a second-order correction of the angular momentum field. Then equation ( 22) recasts

w j = u j ∧ x j -H j b j . (24)
The formula [START_REF] Galera | A two-dimensional unstructured cell-centered multi-material ALE scheme using VOF interface reconstruction[END_REF] can be used to redefine the mean total energy as

e j := ε j + 1 |V j | j |v j (x)| 2 2 , = ε j + |u j | 2 /2 + H j b j • b j /2. ( 25
)
Once again H j b j • b j is a second-order in space correction to the total energy e j .

Proposition 2.3. Since H j is non-singular, there is a one to one correspondence between the physical unknowns u j , w j and the couple a j , b j

b j = H -1 j (u j ∧ x j -w j ), a j = u j -b j ∧ x j . ( 26
)
So a natural extension of the scheme ( 7) is

                     m j d t τ j = r C jr • u r , m j d t u j = - r F jr , m j d t w j = - r F jr ∧ x r . m j d t e j = - r F jr • u r . ( 27
)
This system is naturally closed using equations [START_REF] Loubère | The force/work differencing of exceptional points in the discrete, compatible formulation of Lagrangian hydrodynamics[END_REF] for the reconstructed field, equation [START_REF] Kucharik | Enhancement of Lagrangian slide lines as a combined force and velocity boundary condition[END_REF] to recompute the internal energy, and the equation of state for the pressure. It remains to compute the new fluxes as in equation ( 9). Let us emphasis, we have to take into account the enrichment of the field v j (x) into the equation ( 9). Then the nodal acoustic Riemann solver becomes

A r u r = j A jr v j (x r ) + C jr p j , F jr = C jr p j + A jr (v j (x r ) -u r ) . (28)
Comparing this new solver with the previous one ( 9), we observe that the only difference is that the constant by cell velocity u j has been replaced by the reconstructed value at the node v j (x r ), which is a usual procedure for high-order extension of finite volume Godunov schemes.

Interpretation in terms of Discontinuous Galerkin (DG)

The previous point of view somehow treats inertial and angular momentum as distinct conservation laws. The understanding in terms of DG [START_REF] Loubère | Lagrangian discontinuous Galerkin-type method on unstructured meshes to solve hydrodynamics problems[END_REF][START_REF] Shu | A brief survey on discontinuous Galerkin methods in computational fluid dynamics[END_REF][START_REF] Vilar | Cell-centered discontinuous Galerkin discretization for two-dimensional Lagrangian hydrodynamics[END_REF] methods that we present now is quite different.

Concentrating on the velocity variable, we note the vectorial set of polynomial of degree n in dimension d (that is x = (x 1 , . . . , x d ) ∈ R d ) as

P n = Span (in R d ) x j 1 1 . . . x j d d where 0 ≤ j 1 + • • • + j p ≤ n , so that u ∈ P n is equivalent to the expansion u = 0≤j 1 +•••+j d ≤n a j 1 ,...,j d x j 1 1 . . . x j d d , a j 1 ,...,j d ∈ R d .
The dimension of P n increases with n: dim(P 0 ) = d, dim(P 1 ) = d(d + 1) = d 2 + d and so on. Let us define a new space

Q = Span {a + Bx} ( 29
)
where a ∈ P 0 is arbitrary and the vectorial product is noted with the multiplication by an antisymmetric matrix

B t = -B which is arbitrary is a space of dimension d(d-1) 2 . Therefore dim(Q) = d + d(d-1) 2 = 1
2 dim(P 1 ). Moreover one has the embeddings

P 0 ⊂ Q ⊂ P 1 .
That is Q is an intermediate space between P 0 and P

1 . In dimension d = 2 dim(P 0 ) = 2 ≤ dim(Q) = 3 ≤ dim(P 0 ) = 6. In dimension d = 3, the dimensions are dim(P 0 ) = 3 ≤ dim(Q) = 6 ≤ dim(P 0 ) = 12.
With these notations, the standard Finite Volume framework or assumption recasts as u j ∈ P 0 for all j while the new method recasts as u j ∈ Q for all j.

In DG methods, the space of polynomials is also used for trial functions. We detail only the consequence for the Eulerian impulse equation (3) since it concentrates all the structure and allows simpler notations. Let us write this equation under the form ∂ t (ρu)+ ∇ • M = 0 where M is a symmetric matrix M = M t = ρu ⊗ u + pI d . One integrates in a cell Θ, against the trial function u

d t Θ ρu • u + Θ (∇ • M ) • u = 0, ∀ u ∈ Q.
The usual integration by parts of the second integral yields

d t Θ ρu • u - Θ M : ∇ u = ∂Θ (M u, n) dσ, ∀ u ∈ Q. ( 30
)
Proposition 2.4. The interior integral vanishes for all test functions in Q.

Indeed ∇ u = B is an antisymmetric matrix. Since M is symmetric, the contraction M : B = 0 vanishes. The proof is ended.

Using then any kind of discrete fluxes, one gets a scheme which is ultimately conservative for all test functions in Q. Choosing a trial function with B the null matrix in [START_REF] Maire | A cell-centered Lagrangian scheme for 2D compressible flow problems[END_REF], gives the conservation of inertial momentum. Choosing a trial function with a = 0 gives the conservation of angular momentum. This DG method with the intermediate space Q is just an alternative numerical method to construct a scheme that preserves the angular momentum.

Remark 2.5. It is also reasonable to postulate that any DG methods such as [START_REF] Loubère | Lagrangian discontinuous Galerkin-type method on unstructured meshes to solve hydrodynamics problems[END_REF][START_REF] Shu | A brief survey on discontinuous Galerkin methods in computational fluid dynamics[END_REF][START_REF] Vilar | Cell-centered discontinuous Galerkin discretization for two-dimensional Lagrangian hydrodynamics[END_REF] based on P n (n ≥ 1) applied to the impulse equation [START_REF] Amsden | KIVA-2: A computer program for chemically reactive flows with sprays[END_REF]

yields two family of conservation laws: one for the inertial impulse variable and the other one for the angular momentum variable.

This observation shows that the equation of angular momentum has a very special status. On the one hand it is a conservation law so there is great numerical benefit to use Finite Volume techniques (which are basically P 0 techniques) to discretize this equation.

On the other hand it can be understood as a partial increase of order with the space Q.

In this work we rely more on the first interpretation. The fact that angular momentum corresponds to partial increase of order will appear also evident in section 4.1.

Properties of the new scheme (27-28)

The system [START_REF] Loubère | 3D staggered Lagrangian hydrodynamics scheme with cell-centered Riemann solver-based artificial viscosity[END_REF] shares with the seminal scheme [START_REF] Després | Lagrangian gas dynamics in 2D and lagrangian systems[END_REF] the conservation of mass (and volume), inertial momentum and total energy. Proposition 3.1. Additionally, the semi-discrete scheme [START_REF] Loubère | 3D staggered Lagrangian hydrodynamics scheme with cell-centered Riemann solver-based artificial viscosity[END_REF] preserves total angular momentum j m j d t w j = 0.

(

) 31
The proof of ( 31) is a consequence of [START_REF] Després | Lagrangian gas dynamics in 2D and lagrangian systems[END_REF]. The same property holds of course for usual time discretizations (as in subsection 4.2). We now pay attention to the entropy balance, which is known to be strongly related to the physical relevance of the method and to its numerical stability.

Proposition 3.2. The entropy balance writes

m j T j d t S j = Q j + R j . ( 32
)
where Q j ≥ 0 and R j is a residual term given in [START_REF] Mignone | A conservative orbital advection scheme for simulations of magnetized shear flows with the PLUTO code[END_REF] that will be discussed later on.

The algebra is as follows. One has from the chain rule and the fundamental principle of thermodynamics

                                       m j T j d t S j = m j d t ε j + p j m j d t τ j = m j d t (e j -|u j | 2 /2 -H j b j • b j /2 + p j m j d t τ j = m j d t e j -m j u j • d t u j -m j d t H j b j • H -1 j H j b j /2 + p j m j d t τ j = m j d t e j -m j u j • d t u j -m j b j • d t (H j b j ) -m j (H j b j ) • d t (H -1 j )H j b j /2 + p j m j d t τ j = m j d t e j -m j u j • d t u j -m j b j • d t (H j b j ) + m j b j • d t (H j )b j /2 + p j m j d t τ j = m j d t e j -m j u j • d t u j + m j b j • d t (w j -u j ∧ x j ) +m j b j • d t (H j )b j /2 + p j m j d t τ j = m j d t e j -m j u j • d t u j + m j b j • d t w j -m j b j • d t u j ∧ x j + p j m j d t τ j + R j = m j d t e j -m j a j • d t u j + m j b j • d t w j + p j m j d t τ j + R j .
where the residual R j is

R j = -m j b j • u j ∧ d t x j + m j b j • (d t (H j )b j /2) . (33)
Let us now recast the first part of the right-hand-size of the entropy balance, using the definition of v j (x r ), and the fact that r C jr = 0 (translational invariance of the volume) and r C jr ∧ x r = 0 (rotational invariance of the volume). One has that

                 Q j = m j d t e j -m j a j • d t u j + m j b j • d t w j + p j m j d t τ j , = - r F jr • u r + a j • r F jr -b j • r F jr ∧ x r + p j r C jr • u r , = - r F jr • u r + r F jr • v j (x r ) + p j r C jr • u r , = r (u r -v j (x r ) • A jr (u r -v j (x r )),
where we have used the equation ( 28). The proof is ended.

Remark 3.3. Using the language of thermodynamics, the identity

m j T j d t S j = m j d t e j -m j a j • d t u j + m j b j • d t w j + p j m j d t τ j + R j
shows that -a and b are the intensive variables for the extensive variables u and w.

We now perform a dimensional analysis of both Q j and R j . The first term Q j is obviously non-negative. As explained in [START_REF] Després | Weak consistency of the cell centered lagrangian glace scheme on general meshes in any dimension[END_REF], the dimension of A jr is O(h d-1 ), where we recall that h is the characteristic length of the mesh and d the dimension of the problem. The mass of the cell is m j = O(h d ). For smooth flows, we can reasonably assume u rv j (x r ) = O(h). Therefore Q j = O(h d+1 ), and is non-negative. This term is the generalization in our context of the usual entropy production of Lagrangian schemes. The second term is

R j = -m j b j • u j ∧ d t x j S j + m j b j • d t (H j )b j /2 T j . ( 34
)
It may be unsigned, and can strongly pollute the entropy production term. For example the entropy production becomes negative if Q j +R j < 0. Therefore, it must be controlled. Such a control is equivalent to get a control of the magnitude of b j , u j ∧d t x j and d t H j . To ease the discussion, we assume that b j is bounded: this is reasonable from a theoretical point of view, and is achieved in all our implementations. Then, since H j = O(h 2 ), and d t (H j ) = O(h 2 ), the dimensional analysis tells us that T j = O(h d+2 ). Therefore, we can expect T j << Q j for sufficiently fine grids. Situation is different for S j . Indeed, nothing enforces the equality of u j and d t x j , contrary to nodal quantities for which u r = d t x r . Nevertheless, any physically admissible numerical simulation is such that u j ≈ d t x j . To continue the analysis, we assume that u j ∧ d t x j = O(h). Then, S j = O(h d+1 ) is of the same order of magnitude than Q j . It shows that S j can be the source of a non-physical entropy production, and one can expect stability issues in such cases. Let us discuss how to have a control on S j .

One first possibility to address this issue is to enforce d t x j = u j at the numerical level: a first-order in time numerical discretization could be x n+1 j = x n j + ∆tu n j . Let us emphasis that in this case, x j is no more coincident with the center of mass of the cell, which probably introduces other type of numerical errors in view of equations ( 21) and [START_REF] Hoch | Entropic and Galilean invariant second-order extension for cell-centered Lagrangian schemes[END_REF]. Another option could be to redefine the total energy ( 25) as

e j := ε j + |u j | 2 2 , (35)
since the canceled term is just a second-order contribution to the total energy. In this case, the term R j = 0 vanishes and Q j becomes

Q j = r (u r -u j ) • A jr (u r -v j (x r )). ( 36
)
The drawback is that, once again, Q j is no more signed. In summary, the stability analysis of the new scheme cannot rely only on the entropy production. On the other hand the new variable that provides the preservation of the angular momentum can be interpreted as a kind of second-order discretization of the velocity field. This is why the stability of the scheme is addressed in combination with a usual second-order Muscl extension (and limitation) of the basic scheme. As usual for compressible fluid dynamics solvers, some links are possible between limitations techniques and the entropy increase of the scheme. In practice satisfactory calculations are both non oscillatory (so with convenient limitation techniques) and stable.

Numerical implementation

The numerical implementation has be found challenging because of the previously underlined stability issues, which cannot be addressed only with theoretical tools. We describe in the following the numerical strategy we have employed to obtain an stable and efficient second-order discrete scheme. It has also its own interest since it illustrates the strong numerical interaction of the new degree of freedom with Muscl reconstruction techniques. The second-order Muscl-type approach we use for the scheme without angular momentum conservation is described for instance in [START_REF] Hoch | A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on Local Convex Hull Preservation[END_REF][START_REF] Hoch | Entropic and Galilean invariant second-order extension for cell-centered Lagrangian schemes[END_REF].

Second-order extension and limitation

As stressed previously, the reconstruction performed on the velocity field can also be view as a partial second-order extension. Indeed, combining equations ( 20) and [START_REF] Hoch | A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on Local Convex Hull Preservation[END_REF] gives

v j (x) = u j + B j (x -x j ), (37)
where

B j = -B t j is the d × d anti-symmetric matrix defined by ∀x, B j x = b j ∧ x. (38)
Our objective here is to provide a second-order extension for the new scheme, as close as possible to this one, in order to keep desirable properties as Galilean invariance. The unlimited reconstructed velocity proposed in [START_REF] Hoch | A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on Local Convex Hull Preservation[END_REF][START_REF] Hoch | Entropic and Galilean invariant second-order extension for cell-centered Lagrangian schemes[END_REF] reads

v j (x) = u j + ∇u C j (x -x j ), (39)
where the superscript C in ∇u C j refers to the Classical definition of the gradient. For the new scheme, we assume that the anti-symmetric part B j of the velocity gradient is implicitly imposed by the conservation of angular momentum. Then, we compute the velocity gradient as usual with a least-squares procedure, and modify it the following way

∇u j = D j + B j , D j = 1 2 ∇u C j + ∇u C j t . ( 40
)
It means that the antisymmetric part of the velocity tensor is deduced from the new angular momentum conservation law, while the symmetric part D j remains computed thanks to an usual Mucl-type algorithm. Finally, it leads to the reconstructed velocity

v j (x) = u j + ∇u j (x -x j ). (41)
Once the velocity gradient ∇u j is computed, we have the choice for the limitation strategy. Currently, we use the VIP limiter [START_REF] Luttwak | Slope limiting for vectors: A novel vector limiting algorithm[END_REF] based on the strategy described in [START_REF] Hoch | A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on Local Convex Hull Preservation[END_REF][START_REF] Hoch | Entropic and Galilean invariant second-order extension for cell-centered Lagrangian schemes[END_REF]. Let us emphasis, that we have experienced another strategy, consisting in first limiting B j , and then the symmetric part of the gradient, in order to privilege the rotational part of the reconstruction, but it leads to less satisfactory results.

In the following, the formula ( 37) is referred to as a first-order reconstruction of v j (x), while formula ( 41) is referred as a second-order reconstruction. In the case of secondorder reconstruction, we also use a standard second-order Muscl reconstruction of the pressure, with Barth-Jespersen [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF] limiter

p j (x) = p j + ∇p j • (x -x j ).
(42)

Time discretization and generic work flow

The time integration is performed using an Euler scheme (first-order) or a Runge-Kutta scheme (second-order). For one time-step of the Euler scheme, or one step of the Runge-Kutta algorithm, mean quantities ρ j , u j , w j , p j and e j are initially available in each cell. We observe that an important issue occurs, when using the formula [START_REF] Kucharik | Enhancement of Lagrangian slide lines as a combined force and velocity boundary condition[END_REF] for updating specific internal energy. Since

ε j = e j -|u j | 2 /2 -H j b j • b j /2. ( 43
)
it is obvious the positivity of ε j requires a careful control of the amplitude of b j . In this case, it is found necessary to update ε j just after a convenient evaluation and eventual limitation of b j . In our test problems, we observe that the second-order limitation proposed in subsection 4.1, is sufficient to ensure the positivity of ε j , provided the limitation is performed just before the update of ε j (if not, technics described in [START_REF] Hoch | Entropic and Galilean invariant second-order extension for cell-centered Lagrangian schemes[END_REF] can be applied).

A generic work flow is then the following (we describe here the first-order in time algorithm, second-order Runge-Kutta extension is straightforward). The updated states ρ * j , u * j , w * j , p * j and e * j are obtained from a previous state ρ # j , u # j , w # j , p # j and e # j with the following steps 1. From the mesh, compute the geometrical features

C * jr , x * j , x * r , H * j , 2. Compute b * j from (26): b * j = H * ,-1 j (u * j ∧ x * j -w * j ), (44)
from which the antisymmetric matrix B * j is deduced (refer to equation ( 38)).

3. Limit the reconstructed velocity as described in subsection 4.1, and deduce B # j and eventually D # j (refer to equations [START_REF] Sod | A survey of finite difference methods for systems of nonlinear conservation laws[END_REF] or ( 40)). 4. From equation ( 41), compute

v # j (x * r ) = u * j + ∇u # j (x * r -x * j ), (45)
5. if needed reconstruct a second-order nodal pressure using equation (42):

p # j (x r ) = p * j + ∇p # j • (x * r -x * j ). (46)
6. Solve the system [START_REF] Loubère | Lagrangian discontinuous Galerkin-type method on unstructured meshes to solve hydrodynamics problems[END_REF] with the reconstructed quantities:

A * r u # r = j A * jr v # j (x * r ) + C * jr p # j (x * r ), F # jr = C * jr p # j (x * r ) + A * jr v # j (x * r ) -u # r . (47)
7. Compute u # j , w # j and e # j using equation ( 27):

                 u # j = u * j - ∆t * m j r F # jr , w # j = w * j - ∆t * m j r F # jr ∧ x * r . e # j = e * j - ∆t * m j r F # jr • u # r . (48)
8. Compute the new internal energy using equation (43):

ε # j = e # j -|u # j | 2 /2 -H * j b # j • b # j /2. ( 49
)
9. Move the mesh vertex at the velocity u # r . Compute the new density using the conservation of mass.

10. Update the cell pressure and sound-speed using the equations of state. 11. Compute the new time step.

Axi-symmetric formulation

The axi-symmetric formulation is deduced from the 3D version in averaging along the homogeneous direction (azimuth). Consequently, the velocity field (u r , u z , 0) T in the cylindrical frame, writes in the Cartesian frame (x, y, z), u = (cos θ u r , sin θ u r , u z ) T . The expression for angular momentum is then w = (sin θ(zu rru z ), cos θ(ru zzu r ), 0) T . Injecting this expression into the formula [START_REF] Boscheri | A direct arbitrary-Lagrangian Eulerian ADER-WENO finite volume scheme on unstructured tetrahedral meshes for conservative and nonconservative hyperbolic systems in 3D[END_REF], and recasting it into a continuous formulation for convenience, one gets ρd t (sin θw) = z∂ y psin θr∂ z p, ρd t (cos θw) = z∂ x pcos θr∂ z p,

where w is zu rru z . The two equations of the system (50) are not independent. We can combine them to eliminate θ in multiplying the fist line by sin θ and the second line by cos θ and summing the two resulting equations. Combining with the identities sin θd t (sin θw) + cos θd t (cos θw) = d t w and sin θ∂ y + cos θ∂ x = ∂ r , it yields

ρd t w = z∂ r p -r∂ z p. (51)
This equation admits two natural divergent formulations.The first one is

ρd t w = ∂ r (pz) -∂ z (pr). (52)
By multiplying the equation (51) by r, the second one is

ρrd t w = ∂ r (prz) -∂ z (pr 2 ) -zp. (53)
The equivalent integral form of ( 52) is

d t S(t) ρwdrdz = S(t) (∂ r (pz) -∂ z (pr)) drdz, (54)
where S(t) accounts for the axi-symmetric volume Ω(t). The equivalent integral form of ( 53) is

d t S(t) ρwrdrdz = S(t) ∂ r (pzr) -∂ z (pr 2 ) -zp drdz. (55)
The advantage of the first formulation is that the right-hand-side is in divergent form. We call it the area-weighted formulation, by analogy with the inertial momentum axisymmetric formulation. The advantage of the second formulation is that the Lagrangian mass naturally appears (ρrdrdz). We call it volume-weighted formulation.

We employ the following semi-discrete form to approach the first formulation

ρ j (t)S j (t)d t w j = - r F jr ∧ x r , (56)
where S j (t) is the planar surface of the cell j, and F jr are computed using the planar Riemann solver (including the definition of C jr = ∇ r S j ). We employ the following semidiscrete form to approach the second formulation

m j d t w j = - r (F jr -C jr p j ) ∧ x r , (57)
where F jr are computed using the axi Riemann solver (including the definition of C jr = ∇ r V j ). Refer for instance to [START_REF] Whalen | Algebraic limitations on two dimensionnal hydrodynamics simulations[END_REF] for more details on these two axi-symmetric formulations. We consider that the choice of the angular momentum formulation between (56) and (57) is to be made accordingly to the similar choice for the inertial momentum formulation.

Numerical results

In this section, we present the result of some test problems which are computed with either the basic scheme, meaning without angular momentum variable so not conservative in angular momentum, or the modified scheme, meaning with the angular momentum variable and so conservative in angular momentum. We also consider first-order space discretization and second-order space discretization. It makes four combinations, which are all compared in the rotating shock problem. The tests are performed in dimension two: it simplifies the implementation since the angular momentum variable and the matrix H j are just scalars; moreover reference solutions are much more easy to define in dimension two. We also detail a test problem in 2D axi-symmetric formulation which means a 3D physical problem. The order in time is two thanks to second-order Runge Kutta in time. The CFL is the usual one. We do not report on the extra-cost of the new scheme since it is weak in 2D: indeed we solve for 5 unknowns instead of 4 and the Riemann solver is the same; only quadrature formulas such as [START_REF] Jardin | Review of implicit methods for the magnetohydrodynamic description of magnetically confined plasmas[END_REF] slightly increase the numerical cost.

Cylindrical Sod shock tube

We first propose an irrotational sanity check based on a classical Sod shock problem [START_REF] Sod | A survey of finite difference methods for systems of nonlinear conservation laws[END_REF] in cylindrical configuration. The goal is to measure the intrinsic stability or instability of the new scheme.

The heavy (ρ = 1) compressed (p = 1) fluid is in the outer shell (0.5 < R = √ r 2 + z 2 < 1), and the light (ρ = 0.125) expanded (p = 0.1) fluid is in the inner shell (0 < R < 0.5), so that the flow is convergent. We perform the calculation on an equal angle and equal layers zoned polar mesh with 20 slices and 100 layers visible on figure 1. The calculation is run until t = 0.2.

Our aim is not to argue about the quality of the solution in term of average quantities, what has been widely discussed in previous papers, but to measure the amplitude of b j . Results are displayed on figures 1 and 2. The color map of b shows a maximum value of ||b|| ∞ ≈ 5 × 10 -10 . This value seems reasonable in view of equation [START_REF] Kluth | Discretization of hyperelasticity on unstructured mesh with a cell-centered Lagrangian scheme[END_REF], since |H| ≈ ×10 -7 for this test problem. The density plot shows that the angular momentum conservation has no impact on the quality of the solution for this case.

Rotating shell (2D planar Lagrangian)

This test problem consists in the solid body rotation of a cylinder. The initial configuration is depicted on figure 3. The red and blue part of the mesh have no physical meaning, and just help for the visualization of the solution. We run the problem on a 400 slices ×40 layers mesh until t = 2π, corresponding to one revolution of the cylinder. We initialize the pressure in order to achieve the dynamical equilibrium, and impose zero normal velocity boundary conditions on the external and internal boundaries. Consequently, the problem is theoretically steady, in the sense that initial and final solutions are strictly identical, since we achieve a whole revolution. We use the second-order scheme unlimited scheme, with and without angular momentum conservation.

Without angular momentum conservation, the calculation stop around t = 2.56. It is due to a strong deformation of the mesh at the vicinities of the boundary conditions, as depicted on figure 4. The mean interface between red and blue domain, and consequently mean angular velocity, are almost correct, but the mesh is widely stretched close to boundaries. It was also what we observed in [START_REF] Clair | A one-mesh method for the cell-centered discretization of sliding[END_REF]. See also [START_REF] Kucharik | Enhancement of Lagrangian slide lines as a combined force and velocity boundary condition[END_REF] where a similar behavior is reported on a sliding ring.

The results of the same problem but with angular momentum conservation are displayed on figure 5. The problem runs now until the end, and the mesh suffers almost no modification. In right part of figure 5, we have superposed the blue domain at the end of calculation on the red domain at the beginning of the calculation, and zoom at the vicinity of the interface. Very small differences, almost invisible, can nevertheless be observed at the vicinity of the boundaries.

On figures 6, we focus on the value of the angular momentum w. The reference solution is a constant angular momentum in any cell during the whole calculation. The figure 6 shows little discrepancies of the w field for the calculation without angular momentum conservation, while the w field is identical to the analytical solution for the calculation with angular momentum conservation. It is related to a perfect preservation of the integral angular momentum as shown on figure 7, while calculation without conservation shows a linear decay of this integral.

Taylor Vortex (2D planar Eulerian)

The isentropic Taylor-Vortex test problem is a preliminary and natural candidate to investigate the effect of angular momentum conservation on the accuracy of a flow dominated by vorticity. In particular we will see that it confirms that the original decomposition proposed in section 4.1 preserves the order of the method (second-order in our case). It also validates the coupling with a natural remapping procedure.

The problem is run on a purely Eulerian mode, that is we use a standard La-grange+Remap approach where the Remap stage is a projection on the grid of all conserved quantities. Our Lagrangian and Remap steps use a second-order reconstruction with limiters. To update w during the remap step, we simply project it the same way than u. The initial conditions are

ρ(x) = 1 - (γ -1)β 2 8γπ 2 e 1-r 2 1 γ-1 , p(x) = ρ(x) γ , u(x) = β 2π e 1-r 2 2
where β = 5 corresponds to the strength of the vortex, and r = ||x|| is the distance to the center. We run this problem on four Cartesian meshes 25×25 cells, 50×50 cells, 100×100 cells and 200×200 cells, until t = 2π. For this problem the exterior boundary condition has an unfortunate influence on the numerical diagnostics since the mass, impulse, energy and angular momentum are modified at the boundary cells. Indeed the domain of computation exchanges information with the exterior domain.

The errors ||ρρ analytic || 2 with and without angular momentum conservation (amc) are displayed on the table 1. As expected, the convergence order is almost the same for the two calculations which validates the general approach and the original second-order decomposition proposed in this work. We see that the conservation of angular momentum lowers the error by a factor ≈ 1. We observe that the order of our Muscl procedure is not degraded by the introduction of the angular momentum variable. The results are even slightly better.

Implosion test problems

We show that the numerical accuracy of implosion test problems on skewed meshes such as the ones illustrated in figure 8 takes great advantage of the preservation of angular momentum. To our knowledge this is the first time such a rigorous explanation of the enhanced accuracy of angular momentum preserving schemes for implosion problems is proposed.

We validate the concept on dedicated test problems in cylindrical and spherical geometries. Our test problems are in the same vein as the Gurderley problem for which we refer to [START_REF] Loubère | The force/work differencing of exceptional points in the discrete, compatible formulation of Lagrangian hydrodynamics[END_REF]. Unfortunately we found out that the initial data of the Guderley problem are not fully documented in the open literature. Moreover the Guderley problem does not have an explicit solution: an additional code is needed to generate the reference solution. This is why we propose hereafter our version of the problem where we give all initial values.

Symmetry analysis

We first focus on the symmetry properties for a focusing shock. For practical problems meaning for problems where the mesh is designed before knowing exactly the focusing point as in figure 8, one cannot guaranty that the mesh will be aligned with the shock (polar mesh for a cylindrical or spherical implosion). It may induce mesh imprints in the prediction of center of the implosion, and more generally on the position of the shock with respect to this center. We now establish a theoretical link between inertial and angular momentum conservation and center of implosion M.

Let I = Ω ρu, be the integrated inertial momentum over the whole domain, and W = Ω ρu ∧ x, be the integrated angular momentum over the whole domain. Assume M corresponds to the center of implosion, so: ∃ t f such that ∀x, one has u∧(x-M) = 0. Consequently, we necessarily have I∧M = W. This is a linear system where the unknown is M. Since the matrix is singular, M cannot be fully recovered: however this equation provides one (respectively two) scalar equation(s) in 2D (respectively 3D) from which M can be partially determined. The remaining equation can be for instance provided by symmetry considerations, such has the ones considered in our test problems.

It implies that if sources terms (for instance due to boundary conditions) are negligible, and when a symmetry consideration is provided, the center of implosion M can be determined from initial conditions. For numerical solutions, if the scheme used for the simulation is conservative in inertial momentum and angular momentum, then the same proof can be performed. In this case, the numerical scheme will predicts the correct value of M.

Cylindrical problem (2D planar Lagrangian)

To verify the above assertion, we consider a flying shell test problem on a mesh not aligned with the shock. The configuration is depicted on figure 8.

The blue domain corresponds to a heavy perfect gas (ρ = 1, p = 1) with an initial incoming radial velocity u r = -30. The shell is completely centered with respect to M, so that initially u ∧ (x -M) = 0 for all x. The internal radius of the shell is 6.456, while outer radius is 7.5. It implodes on a perfect light gas (ρ = 10 -2 , p = 1) at rest, corresponding to the red part of the mesh. As depicted on the figure, the corresponding mesh is not aligned with the flow, and the gap between M, and the center of the O-grid is 1.2912. The boundary condition is symmetry on the axis y = 0 and prescribed pressure p = 10 on the outer surface of the shell, so that the shell is initially accelerated. We use a quite coarse mesh, to exhibit clearly the gain obtained by the conservation of the angular momentum. 20 slices are used for angular discretization. We use 10 layers in the shell and 50 in the inner material. We perform this calculation in planar configuration. The shock focuses approximately at t = 0.135. We run the calculation until t = 0.16 and look at the shape of the rebounded shock wave.

For first-order calculations, we depict on figure 9 the pressure maps corresponding with and without the conservation of angular momentum. These figures show the huge enhancement in term of symmetry of the rarefaction wave obtain by the conservation of angular momentum. We have plot on figure 10, the pressure in all the cells at t = 0.16 for the two calculations. It confirms the important gain in symmetry.

For second-order calculations, we depict on figure 11 the pressure maps corresponding with and without the conservation of angular momentum. The gain obtained by the conservation of angular momentum is less visible, even if one observes a small shift of the shock front with respect to the black circle for the calculation without angular momentum conservation. We assume that the global quality of the shock position is already widely improved by the second-order extension. However the 1D plots, on figure 12, shows much less spreading in the smooth part with angular momentum conservation than without.

Spherical problem (Lagrangian 2D axisymmetric)

Spherical implosion problem are naturally stiffer than cylindrical ones. In consequence, the initial condition of the spherical test problem must be adapted from the cylindrical one. We consider here a simpler configuration, a sphere uniformly filled by a perfect gas of initial density (ρ = 0.01) and pressure (p = 1) as depicted on the right part of figure 8. The total energy in the domain is much less than for the cylindrical problem. On the other hand a higher pressure 10 5 is applied on the external surface of the sphere imploding the gas. We keep the same mesh than for the cylindrical test problem. We run the problem using the volume-weighted formulation (refer to annexe 5 for more details). The shock focuses approximately at t = 5 × 10 -3 .

For first-order calculations, we depict on figure 14 the pressure maps corresponding with and without the conservation of angular momentum. The shock positions are almost the same with and without angular momentum conservation, which is the effect of the preservation of angular momentum. One can argue a slight amelioration on the axis with preservation of the angular momentum: in any case this is less visible compared to the first-order cylindrical calculations reported in figure 9. It was somewhat expected, since our modification of the scheme has no impact on 1D solutions and spherical implosion problems are more dominated by 1D implosion than cylindrical problems are. The figure 15 shows some enhancement in term of symmetry of the rarefaction wave obtained by the conservation of angular momentum. The pressure is plot on figure 15 in all the cells at t = 4.6 × 10 -3 for the two calculations. We observe strong symmetry discrepancies into the tail of the profile, due our choice of axi-symmetric formulation (which does not respect symmetry) on this coarse mesh.

Second-order calculations are shown on figure 16 the pressure maps corresponding to with and without the conservation of angular momentum. For this axi-symmetric test problem the gain obtained by the conservation of angular momentum is here visible. The general quality of the second-order calculation without preservation of the angular momentum is similar to the one of a first-order calculation. One the other hand the second-order calculation with preservation of the angular momentum seems to be much more symmetric. The 1D radial cuts, on figure 17, strengthen this observation. The pressure in all the cells at t = 4.6 × 10 -3 for the two calculations, shows less spreading with angular momentum conservation than without.

Rotating shock tube (2D planar Lagrangian)

Last test problem, a rotating Sod shock tube, is completely dedicated to this study. The initial and boundary conditions are described on figure 18.

The two shells rotate initially with the same angular velocity, and are in equilibrium, except at the contact discontinuity, where pressures and densities are discontinuous. The velocity is continuous. It is a 1D test problem, in the sense that all quantities depend only on the radius. Consequently, the pressure gradient and the position are always aligned and m j d t w j = V j ∇p ∧ x = 0 all along the calculation. The radial interfaces are fictitious and, as for test problem 6.2, only for visualization purposes. The circles correspond to the inner surface (R int , zero normal velocity boundary condition), contact discontinuity (R shock ), and outer surface (R ext , zero normal velocity boundary condition).

The problem is run until t = 0.1, corresponding to a quarter of a revolution, on a 90 slices ×50 layers polar grid. We use the first-order and second-order extensions in space and time, and a little amount of subzone entropy [START_REF] Després | Stabilization of cell-centered compressible Lagrangian methods using subzonal entropy[END_REF] for stability purpose, and we compare the results to a reference calculation with the second-order extension but on a 180 slices ×100 layers grid (angular momentum conservation has very few visible effect with this resolution). We display on figure 19 the results obtained for the reference calculation, in term of Lagrangian deformation of the initial calculation domain. Figure 20 corresponds to a zoom on the bottom part of the figure 19.

The calculation domain is heavily deformed. This is expected since w j = u j ∧ x j has to be preserved in each cell. Then, the more the cells are brought nearer to the center, the more their velocities increase. We state that the higher deformation corresponds to the contact discontinuity, what was expected, since the shock is initialized at this spot. Moreover, because of zero normal velocity boundary conditions, theoretical final positions of boundary nodes are known. Since these vertex have constant radius during the calculation, they should exactly have done a quarter of a revolution at final time. In particular, we focus on vertex M int , M shock and M ext , which initial location are displayed The result obtained in using the first-order scheme without angular momentum on the coarse mesh, are displayed on figures 22 and 21 (zoom). We observe large discrepancies with respect to the reference calculation, corresponding to the white line for the interfaces, and the green circles for the vertex locations.

In enforcing the angular momentum calculation (see figures 23 and 24), we find here an major enhancement on the quality of the results with respect to figures 22 and 21. Indeed M ext and the interface position are now correct, even if M int is still quite far to the exact position. It confirms nevertheless that angular momentum conservation is an essential feature for this problem. Now on figure 25, we have depicted the results for the second-order scheme without angular momentum conservation. Figure 26 corresponds to a zoom on the bottom of the figure 25. The accuracy is similar to the one for the first-order scheme with angular momentum conservation, the latter giving actually slightly better results. Locations of M int and M ext are still very far from the theoretical values.

Performing a second-order calculation, but with angular momentum conservation, improves the results a lot (see figures 27 and 28) which are almost perfect even on this coarse grid. It shows that many of the numerical discrepancies of the previous calculations are related to the non conservation of angular momentum. For this configuration, we also display the mesh on figure 29. Cells are highly streched, in particular at the viscinity of the contact discontinuity.

Since fluid are perfect gas, we can calculate explicitely an entropy. We plot on figure 30 the entropy ratii profiles as a function of the radius. First statement is that entropy ratio is never below 1, what means that entropy increases everywhere, even for the schemes with angular momentum conservation. It demonstrates that we correctly control the difficulty discussed in section 3 concerning the entropy production. As expected, entropy production is much higher for the first-order schemes than for the second-order ones, in particular at the shock front. For the first-order scheme, the conservation of angular momentum drastically decreases the entropy production. The peak corresponding to the shock position is misplaced for the first-order scheme without angular momentum conservation, what is coherent with the color maps. For second-order schemes, there is much lesser difference between the curves, except close to the inner boundary.

Another interesting feature of the use of the angular momentum is that the overall scheme seems to be much less dependent on any kind of additional stabilization procedure that one might think to add. In our tests we use often a small amount of the additional stabilization procedure called subzone entropy techniques, as described in [START_REF] Després | Stabilization of cell-centered compressible Lagrangian methods using subzonal entropy[END_REF]. Canceling the subzone entropy everywhere, except on the internal boundary (where it is really needed, to avoid the simulation crash), we obtain for the original scheme not conservative in angular momentum the result depicted on figures 32 and 31. This is a far better result than the one depicted on figures 25 and 26, even if a spurious behavior remains at the vicinity of the internal boundary, where subzone entropy is non-zero. So it shows that the stabilization procedure based on subzone entropy increases a lot the overall dissipation of the scheme for this highly dominated by rotation test problem. On the contrary subzone entropy has almost no impact on the result for the new scheme. This is why we consider that conservation of angular momentum is a valuable enhancement of this scheme, since it makes the result independent of the stabilization procedure.

Finally it seems to us the conservation of angular momentum has no impact on the converged result, since both second-order schemes give the same result on the fine 180

Conclusion

We detailed in this work a way to preserve angular momentum for cell-centered Lagrangian solvers on general grids. The method is simple since it introduces few local additional degrees of freedom and does not change the basic Riemann solver. This is why this method is relatively easy and cheap to implement in an existing code, as we have done either in a purely Lagrangian mode or in a Eulerian mode with a simple La-grange+Remap approach. The method corresponds also to a partial DG-like increase of order of the velocity, and the compatibility with second-order Muscl reconstruction is easy. Even if we concentrated on Lagrangian schemes, the basic principles of the method can be used for many compressible flow solvers. In this direction, further researches could be devoted to enhanced remapping-transport algorithms for the angular momentum, a topic that was only marginally developed in this work.

Numerical tests show a strong enhancement of the numerical accuracy for pure rotation test problems, for implosion problems and finally for rotation test problem with shock. Two interesting features is that the method can be much more accurate than usual second-order Muscl extensions and that it strongly reduces numerical discrepancy at sliding boundaries. Our explanation is clearly that conservation laws, and angular momentum satisfies such a conservation law, are fundamental equations for compressible flows that require conservative discretization. It seems to us the method has strong potential for any locally, or globally, rotation dominated flows.

Figure 1 :

 1 Figure 1: Cylindrical Sod shock tube. Mesh and map of b ≈ 5 × 10 -10 at t = 0.2.

Figure 2 :

 2 Figure 2: Cylindrical Sod shock tube. Scatter plot of the density versus radius for calculations with (blue symbols) and without (red symbols) angular momentum conservation. The results are the same.

Figure 3 :

 3 Figure 3: Initial mesh for the rotating shell problem. The blue and red parts are defined on purpose of numerical visualization.

Figure 4 :

 4 Figure 4: Calculation without angular momentum conservation. Black line corresponds to theoretical interface between red and blue domain at t = 2.56. Left: mesh for the rotating shell test problem at t = 2.56; the calculated solution is globally correct. Right: zoom; at boundaries one sees the huge error between the exact solution and the calculated solution. Waiting a little longer, the calculation crashes due to extreme deformation of the mesh.

Figure 5 :

 5 Figure 5: Calculation with angular momentum conservation. Left: mesh for the rotating shell test problem at t = 2 × π. Right: zoom, blue domain is at the end of the calculation, and red at the beginning of the calculation. The mesh is globally and locally correct.

Figure 6 :

 6 Figure 6: Rotating shell. Plot of the angular momentum at t = π/2: lower half shell corresponds to the analytic solution; upper left quarter shell corresponds to calculation without angular momentum conservation (no amc); upper right quarter corresponds to calculation with angular momentum conservation (amc). The important boundary error without angular momentum conservation is visible at the top-left interior boundary.

Figure 7 :

 7 Figure 7: Rotating shell. Right: time evolution of the integrated angular momentum for calculation without angular momentum conservation (blue curve), and calculation with angular momentum conservation (red curve).As expected the integrated angular momentum is constant for the calculation with the new scheme, and is dissipated for the calculation for the standard cell-centered.

Figure 8 :

 8 Figure 8: Example of a initial skewed mesh, which is the one used for the flying shell problem. We also plot the initial data of the cylindrical problem detailed in section 6.4.2. The center of implosion is at M.

Figure 9 :

 9 Figure 9: Pressure maps for the cylindrical flying shell problem at t = 0.16. First-order calculations. Black circle is centered on M and corresponds to the reference solution. Right and left part of the map are troncated to focus on the center of the configuration. Irregular aspect is due to the poor resolution. Left: without angular momentum conservation. Right: with angular momentum conservation.

Figure 10 :

 10 Figure 10: Pressure plot for the cylindrical flying shell problem at t = 0.16 (blue symbols). First-order calculations. Black dots correspond to a reference solution on a very fine mesh. Left: without angular momentum conservation. Right: with angular momentum conservation.

Figure 11 :

 11 Figure 11: Pressure maps for the cylindrical flying shell problem at t = 0.16. Second-order calculations. Black circle is centered on M. Right and left part of the map are troncated to focus on the center of the configuration. Irregular aspect is due to the poor resolution. Left: without angular momentum conservation. Right: with angular momentum conservation.

Figure 12 :

 12 Figure 12: Pressure versus radius plot for the spherical flying shell problem at t = 0.16 (blue symbols).Second-order calculations. Black dots correspond to a reference solution on a very fine mesh. Left: without angular momentum conservation. Right: with angular momentum conservation. The dispersion is more important without preservation of angular momentum.

Figure 13 :

 13 Figure 13: The spherical implosion problem.

Figure 14 :

 14 Figure 14: Pressure maps for the spherical implosion problem at t = 4.6 × 10 -3 for the first-order schemes. White circle is centered on M (green circle) and corresponds to the reference solution on a very fine mesh. Zoom at the center of the configuration. Irregular aspect is due to the poor resolution. Left: without angular momentum conservation. Right: with angular momentum conservation.

Figure 15 :

 15 Figure 15: Pressure plot for the spherical implosion test problem at t = 4.6 × 10 -3 , for the first-order schemes. Left: without angular momentum conservation. Right: with angular momentum conservation.

Figure 16 :

 16 Figure 16: Pressure maps for the spherical implosion problem at t = 4.6 × 10 -3 for the second-order schemes. Black circle is centered on M and corresponds to the second-order solution. Zoom at the center of the configuration. Irregular aspect is due to the poor resolution. Left: without angular momentum conservation. Right: with angular momentum conservation.

Figure 17 :

 17 Figure 17: Pressure versus radius plot for the spherical implosion problem at t = 4.6 × 10 -3 for the second-order schemes. Left: without angular momentum conservation. Right: with angular momentum conservation.

Figure 18 :

 18 Figure 18: Rotating shock tube problem configuration.

Figure 19 :

 19 Figure 19: Rotating shock tube problem. Shape of the calculation domain at t = 0.1 for the reference calculation (two times finer in each dimension). Green circles correspond to the final positions of the vertex: initial positions are depicted in figure 18.

Figure 20 :

 20 Figure 20: Reference calculation for the rotating shock tube problem. Zoom on the bottom part of figure 19.

Figure 21 :

 21 Figure 21: Rotating shock tube problem. First-order scheme, without angular momentum conservation. Zoom on the bottom part of figure 22. The red circle corresponding to M int is out of the figure.

Figure 22 :

 22 Figure 22: Rotating shock tube problem. Shape of the calculation domain at t = 0.1 for the first-order scheme, without angular momentum conservation. Red circles correspond to the final positions of the vertex, which initial position was depicted on figure 18. Green circles correspond to the position of the same vertex, but for the reference calculation. White line is the interface corresponding to the reference calculation.

Figure 23 :

 23 Figure 23: Rotating shock tube problem. Shape of the calculation domain at t = 0.1 for the firstorder scheme, with angular momentum conservation. Red circles correspond to the final positions of the vertex, which initial position was depicted on figure 18. Green circles correspond to the position of the same vertex, but for the reference calculation. White line is the interface corresponding to the reference calculation.

Figure 24 :

 24 Figure 24: Rotating shock tube problem. Zoom on the bottom part of figure 23. The red circle corresponding to M int is out of the figure.

Figure 25 :

 25 Figure 25: Rotating shock tube problem. Shape of the calculation domain at t = 0.1 for the secondorder calculation without angular momentum conservation. White line corresponds to the profile of the calculation domain for a reference calculation. Red circles correspond to the final positions of the vertex, which initial position was depicted on figure 18. Green circles correspond to final positions of the same vertex, but for the reference calculation.

Figure 26 :

 26 Figure 26: Rotating shock tube problem. Zoom on the bottom part of figure 25. The red circle corresponding to M int is out of the figure.

slices

  ×100 layers grid.

Figure 27 :

 27 Figure 27: Rotating shock tube problem. Shape of the calculation domain at t = 0.1 for the secondorder calculation with angular momentum conservation. White line corresponds to the profile of the calculation domain for a reference calculation. Red circles correspond to the final positions of the vertex, which initial position was depicted on figure 18. Green circles correspond to final positions of the same vertex, but for the reference calculation. When green and red circles are superposed, we omit the green circle.

Figure 28 :

 28 Figure 28: Rotating shock tube problem. Zoom on the bottom part of figure 27.

Figure 29 :

 29 Figure 29: Rotating shock tube problem. Mesh at t = 0.1 for the second-order calculation with angular momentum conservation.

Figure 30 :

 30 Figure 30: Rotating shock tube problem. Ratio of the entropy at final time over initial entropy as a function of the radius. Left: first-order schemes. Right: Second-order schemes. Blue curves: not conservative for the angular momentum. Red curves: conservative for the angular momentum.

Figure 31 :

 31 Figure 31: Zoom on the bottom part of figure 32.

Figure 32 :

 32 Figure 32: Rotating shock tube problem. Shape of the calculation domain at t = 0.1 for the calculation without angular momentum conservation, and subzone entropy only for the cells at the vicinity of the internal boundary. White line corresponds to the profile of the calculation domain for a reference calculation. Red circles correspond to the final positions of the vertex, which initial position was depicted on figure 18. Green circles correspond to final positions of the same vertex, but for the reference calculation. When green and red circles are superposed, we omit the green circle.

Table 1 :

 1 7. Convergence table for the Taylor Vortex test problem.

		with amc		no amc
	number of cells	L 2 error	order	L 2 error	order
	25 × 25 50 × 50 100 × 100 200 × 200	1.09 × 10 -3 2.86 × 10 -4 1.93 4.98 × 10 -4 1.42 1.33 × 10 -3 8.31 × 10 -5 1.78 1.83 1.4 × 10 -4 2.17 × 10 -5 1.94 3.73 × 10 -5 1.91