1 Introduction

We are chiefly concerned with solving the two-dimensional compressible gas dynamics equation written in the Lagrangian form. In this paper,we intent to construct a high-order Lagrangian cell-centered conservative gas dynamics scheme on unstructured meshes. Before we formulate our algorithm,we briefly review the Lagrangian gas dynamics algorithm.

A Lagrange fluid mechanics algorithm is featured by the cell with the fluid motion,which means that the cell vertices advance with a computed velocity,and the cell face is determined by the vertex location. This guarantees that there is no mass flux passing through the boundary of the Lagrangian moving cell. The Lagrangian fluid dynamics algorithm can ferret out contact discontinuities in the multi-material flow. However,in the Lagrangian framework,the numerical fluxes of the physical conservation laws must be compatible with the vertex velocity computation. The algorithm to meet this requirement is to use a staggered discretization in which the position,velocity,and kinetic energy are centered at the points,while the density, pressure,and internal energy are within the cells. The dissipation of kinetic energy into internal energy through shock waves is guaranteed by an artificial viscosity term. Based on the pioneering work of von Neumann and Richtmyer^[1] and Wilkins^[2],many Lagrangian staggered discretization schemes were developed,which improved the accuracy and the robustness of the Lagrangian staggered discretization scheme^[3,4,5],and a compatible staggered discretization deduced total energy conservation scheme in a rigorous manner was constructed^[6,7]. Another choice to the previous discretizations is to establish a Lagrangian cell-centered scheme according to the Godunov method^{[8,9,10,11,12]}. Unlike staggered discretizations,cell-centered schemes have the characteristics of conservation,which have no artificial viscosity term,and permit the direct implementation of the conservation remapping method in the context of the arbitrary Lagrangian Eulerian (ALE) algorithm. The major supposition of all of these algorithms is that the cell masses do not vary with time. Caramana et al.^[6] extended the idea of a Lagrangian zone mass to define subzonal Lagrangian masses,and eliminated the artificial cell distortion and hourglass type motion generated in the presence of irregular meshes by the correct utilization of subzonal Lagrangian masses and associated densities and pressures. Ge^[13] constructed a scheme with the piecewise constant pressure of the cell. In the following,we describe a high-order Lagrangian cell-centered conservative gas dynamics scheme on general polygonal meshes. With the high-order piecewise pressure of the cell constructed by means of the monotonic upwind scheme for conservation laws (MUSCL)^[12],the high-order subcell forces are constructed. The vertex velocities are evaluated by means of momentum conservation. A high-order Lagrangian cell-centered conservative gas dynamics scheme on unstructured meshes is performed by means of the Taylor expansion of fluxes in time centered. The remainder of this paper is structured as follows. Previous work is revisited in Section 2. A high-order Lagrangian cell-centered conservative gas dynamics scheme on unstructured meshes is detailed in Section 3. Extensive numerical experiments are reported in Section 4. Conclusions are given in Section 5. 2 Notation and previous work 2.1 Lagrangian hydrodynamics

We discretize the following equations of the Lagrangian hydrodynamics:

where d/ dt expresses the time derivative. V (t) expresses the control volume,and S(t) expresses its boundary. ρ,U = (u,v)^T,P,and E are the density,velocity,pressure,and specific total energy of the fluid,respectively. N expresses the unit outward normal vector to the moving boundary S(t). Equations (1),(3),and (4) indicate the conservation of mass,momentum,and total energy. Equation (2) is also named as the geometric conservation law (GCL). It can also be expressed as the local kinematic equation where X expresses the coordinates defining the control volume surface at time t > 0,and x expresses the coordinates at time t = 0. Then,X = X(x,t) can be implicitly determined by the local kinematic equation,which is also named as the trajectory equation. The equation of state is where ε = E − 1/ 2||U||². 2.2 Notation and assumptions

Let us consider the physical domain V (0) that is initially filled with the fluid. It is mapped by a set of polygonal cells Ω_c (c = 1,2,· · · ,I) without voids or superpositions,where I is the total number of the cells. Each cell is designated by an index c,and is expressed by Ω_c. Each vertex of the cell Ω_c is designated by an index p. We express by p(c) the counter clockwise ordered list of the points of the cell. The vertices of the cell Ω_c satisfy periodicity,i.e.,pp(c)+1 = p1. We express by c(p) the counter clockwise ordered list of the set of the cell around the vertex p. c(p) satisfies periodicity,i.e.,Ω_c(p)+1 = Ω1. p− and p⁺ express the previous and next points with respect to p in the ordered list of the cell vertices (see Fig. 1).

2.3 Problem statement

Maire^[14] discretized the momentum flux by introducing two pressures and at each node p of the cell Ω_c (see Fig. 1). Using these nodal pressures,Maire^[14] constructed a high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes.

where ac,n denotes the local isentropic sound speed,and Γc denotes a material-dependent parameter. The unit corner vector Npc,n is denoted by M_pc,n is a 2 × 2 matrix expressed by

Set

Then,

Maire^[14] claimed that if (9) was satisfied,then the momentum and total energy could be conserved. However,we find that (9) is incorrect,because Maire^[14] introduced two pressures at each node p of the cell Ω_c while the pressure discontinuous line exists pqc inside the cell Ω_c (see Fig. 1). Similarly,the pressure discontinuous lines pqc (c = 1,2,· · · ,c(p)) exist inside each cell Ω_c (c = 1,2,· · · ,c(p)). Due to the presence of these pressure discontinuous lines pqc (c = 1,2,· · · ,c(p)),(9) is incorrect. 3 Conservative spatial approximation 3.1 Compatible discretization of GCL

According to Fig. 2,(2) over the cell Ω_c can be rewritten as

where ∂Ω_c is the boundary of Ω_c,N is the unit outward normal to ∂Ω_c,and s is the length element on ∂Ω_c. Introduce the discrete divergence operator over cell Ω_c as follows: The volume of the cell Ω_c,V_c,is a function of the coordinate Xp of the point p for p ∈ p(c). We compute this volume by performing the triangular decomposition of the cell Ω_c displayed in Fig. 2. where the unit vector e_Z is expressed by e_Z = e_X×e_Y . (e_X,e_Y ) is the orthonormal basis,where e_X,e_Y ∈ R². The time differentiation of this equation is Combining the previous results and this definition,we obtain 3.2 Computation of momentum flux

We next define what is meant by a piecewise constant pressure of a polygonal cell Ω_c. Although these expositions are given according to a quadrilateral cell in two dimension,all of the expositions are used to any type of cells in any dimension. In Fig. 3,we show a quadrilateral cell Ω_c,and pick out one of its defining vertices p. We label the side midpoints of a given cell Ω_c as b₁,b₂,b₃,and b₄,respectively. Employing the cell vertices p,p⁺,pk,and p−,the side midpoints b₁,b₂,b₃,and b₄,and the center point c of the cell Ω_c,we can partition a quadrilateral Ω_c into four quadrilateral subcells cb₁pb₂,cb₂p⁺b₃,cb₃pkb₄,and cb₄p−b₁. Considering the quadrilateral subcells cb₁pb₂,cb₂p⁺b₃,cb₃pkb₄,and cb₄p−b₁,we can express the pressure by (p = 1,2,· · · ,p(c)). Then,the pressure of the cell Ω_c can be expressed as piecewise constant pressures (p = 1,2,· · · ,p(c)) in the quadrilateral subcells cb₁pb₂,cb₂p⁺b₃,cb₃p_kb₄,and cb₄p−b₁. We know the density ρc,0 and the coordinates of the vertices (Xp,0,Yp,0) (p = 1,2,· · · ,p(c)) of each cell Ω_c at the initial time. We calculate the quadrilateral subcell volumes (p = 1,2,· · · ,p(c)) at the initial time. We calculate the quadrilateral subcell masses (p = 1,2,· · · ,p(c)) from the initial density ρc,0 and the initial quadrilateral subcell volumes (p = 1,2,· · · ,p(c)). Then,we know the fluid variables ρc,n,ac,n,Pc,n and the geometrical characteristics Xp,n,Yp,n (p = 1,2,· · · ,p(c)). We calculate the current time subcell volumes(p = 1,2,· · · ,p(c)),the current time subcell densities and the current time subcell pressures:

Defining the high-order piecewise pressure of the cell H^c _p,n (p = 1,2,· · · ,p(c)) as follows^[12]: where ac,n is the isentropic sound speed of the cell Ω_c,and Δt is the time step. The forces created by the high-order piecewise pressure can be calculated along the boundaries of these subcells. We define a high-order subcell force related to the point p and the cell Ω_c by The length d and the unit outward normal D are related to the edge cb₁ of the quadrilateral subcell cb₁pb₂,and the length d and the unit outward normal D are related to the edge cb₂ of the quadrilateral subcell cb₂p⁺b₃. Then,the momentum equation can be written as To examine momentum conservation,we establish the global balance of momentum without considering the boundary conditions. The sum of the momentum equation (19) over all the cells is Switching the sum over cells and the sum over vertices in the right-hand side,one can obtain The condition of momentum conservation is Since the sum of all forces exerting on one point p is equal to zero,the subcell forces satisfy the condition (22). 3.3 Computation of energy flux 3.4 Construction of solver at vertice velocity U_p,n

We can introduce only one pressure at the point p. This pressure is created by the half Riemann problem defined in the direction of the unit corner vector N_pc,n.

where p = 1,2,· · · ,p(c),and Zc is the acoustic impedance of the cell Ω_c. The subcell force corresponding to the nodal pressure Πpc for each cell Ω_c that surrounds the point p (see Fig. 4) reads where is a 2 × 2 symmetric positive definite matrix expressed by

This generic vertex p is not on the boundary of the domain so that it is surrounded by the cells Ω_c (c = 1,2,· · · ,c(p)) (see Fig. 4). Since the sum of all forces exerting on one point p is equal to zero,we get

Using conditions (15),(17),and (22),we can establish the linear system satisfied by the components (u^*_p,n,v^*_p,n) of the velocity Up,n of an internal vertex (see Fig. 5). Substitute (30) in (31),and set

Then,the system satisfied by the point velocity Up,n can be written as

The vector S is defined by

where

The right-hand sides of (34) are the components of the vector S. The determinant of (34), Δ = AB − C²,is always positive. We have shown that the system (34) can determine the vertex velocity (u^*_p,n,v^*_p,n). 3.5 Second-order time discretization

The current paper designs a version of second-order time integration based on the Taylor expansion of the fluxes in time centered. All the physical and geometrical variables are supposed to be known at the beginning of time tn. We express them using the subscript n. Their updated values at time t_n+1 = tn + Δt,where Δt is the current time step,are obtained by means of

3.6 Time derivatives of subcell force and nodal velocity

Using the chain rule derivative,we can obtain the time differentiation of (18) as follows:

We calculate the time derivative of the physical variables by utilizing the gas dynamics equations in the non-conservative form: Using the chain rule of derivative,we can obtain the time differentiation of (32) as follows: The corner normal and the corner matrix are expressed as functions of the normal vectors to the two edges of the cell c intersecting at the node p.

where Z^c is the acoustic impedance. The time derivative of the corner matrix is calculated by the following chain rule:

As the matrix is symmetric and positive definite,(38) gives a unique time derivative for the nodal velocity. 3.7 Momentum conservation

Let us write the global balance of momentum without considering the boundary conditions. The sum of the momentum equation (35) over all the cells obtains

Switching the sum over cells and the sum over nodes in the right-hand side of (41),one gains The condition of momentum conservation is Since the sum of all forces intersecting at one point is equal to zero,the condition (43) is satisfied,and the momentum is conserved. 3.8 Boundary conditions

In the Lagrangian formalism,we must consider two types of boundary conditions. One is that the pressure is imposed on the boundary,and the other is that the normal component of the velocity is imposed on the boundary. Set p on the boundary. p is arounded by c(p) cells contained in the domain. There are c(p)+1 edges intersecting at p. They are numbered counter clockwise (see Fig. 6). The edge [M₁,p] of the cell Ω1 is on the boundary. The edge [M_c(p)+1,p] of the cell Ω_c(p) is on the boundary. and are the outward normals to the two edges [M₁,p] and [M_c(p)+1,p] intersecting at p.

Case 1 A prescribed pressure

Π^*₁ and Π^*_c(p)+1 express the pressures that are imposed on the edges [M₁,p] and [M_c(p)+1,p] (see Fig. 6). We establish a balance around the vertex p by means of the boundary conditions and conservation relation. Then,we can obtain

After substituting (44) in (45),we can write the system satisfied by the point velocity U_p,n as follows: The vector is defined by The right-hand sides of (47) are the components of the vector . The system (47) can determine the velocity (u^*_p,n,v^*_p,n).

The time differentiation of (47) gains

and are the components of the vector defined by

The linear system (48) can provide the time derivative of the vertex velocity ( du^*_p,n/ dt,dv^*_p,n/ dt).

Case 2 A prescribed normal velocity

Let W₁ and W_c(p)+1 be the values of the prescribed normal velocities on the boundary edges [M₁,p] and [M_c(p)+1,p] intersecting at p. We discern the following two cases.

and are not colinear.

The value of the vertex velocity Up,n is computed by the boundary conditions. The components (u^*_p,n,v^*_p,n) of the vertex velocity Up,n are the solution of the following linear system:

The time differentiation of the system (49) obtains This linear system (50) can determine the time derivative of the vertex velocity ( du^*_p,n/ dt, dv^*_p,n/ dt). Since the normals are not colinear,the coefficient determinant of the system (50) is not zero.

(ii) and are colinear.

We cannot give directly the vertex velocity U_p,n. We need to know the balance around the vertex p which considers the boundary conditions. We gain

After substituting (51) in (52),we can obtain the following equation:

The system satisfied by the vertex velocity Up,n is written as

The vector is defined by

where Π^* is a mean pressure imposed on the external side of the edges [M1,p] and [Mc(p)+1,p] intersecting at p. We do not know what is the value of the pressure. However,we can build an additional equation by the boundary condition

We obtain

Here, and are the components of the vector . We set

The determinant of the liner system (54) is Δ = −AH² + 2CGH − BG². Using the fact that AB − C² > 0,we can obtain Δ < 0. The system (54) can solve the vertex velocity (u^*_p,n,v^*_p,n). The time differentiation of (54) obtains

Here,d/ dt and d/ dt are the components of the vector d/ dt defined by

The system (55) can determine the time derivative of the vertex velocity ( du^*_p/ dt,dv^*_p/ dt). 3.9 Time step

The time step Δt is evaluated by the following criterion. At time tn,for each cell Ω_c,we express by λc,n the minimal value of the distance between two vertexes of the cell Ω_c. We define

where is the maximal value of the local isentropic sound speeds of the quadrilateral subcells of the cell Ω_c. 3.10 Description of algorithm

(i) Initialization

We know the density ρ_c,0 and the coordinates of the vertices (X_p,0,Y_p,0) (p = 1,2,· · · ,p(c)) of each cell Ω_c. We calculate the quadrilateral subcell masses (p = 1,2,· · · ,p(c)) from the initial density ρ_c,0 and the coordinates of the vertices (X_p,0,Y_p,0) of Ω_c (p = 1,2,· · · ,p(c), c = 1,2,· · · ,I). At the time t = t_n,we know the fluid variables U_c,n,E_c,n,ρ_c,n,a_c,n,P_c,n and the geometrical characteristics (X_p,n,Y_p,n) (p = 1,2,· · · ,p(c)) in each cell Ω_c. We calculate the quadrilateral subcell volumes (p = 1,2,· · · ,p(c)) and the densities of quadrilateral subcells (p = 1,2,· · · ,p(c)). Using the pressure P_c,n and the densities of the quadrilateral subcells (p = 1,2,· · · ,p(c)),we can get the isentropic sound speed of the quadrilateral subcells (p = 1,2,· · · ,p(c)). Using (15),(17),(18),and (36),we can calculate the high-order subcell force and the time differentiation of the high-order subcell force d/ dt .

(ii) Nodal solver

For each internal vertex,we calculate the velocity (u^*_p,n,v^*_p,n) and the time differentiation of the velocity ( du^*_p,n/ dt,dv^*_p,n/ dt) by solving the linear system (34) and (38). For each boundary vertex,we calculate the velocity (u^*_p,n,v^*_p,n) and the time differentiation of the velocity ( du^*_p,n/ dt, dv^*_p,n/ dt) by solving the linear system (47),(48),(49),(50),(54),and (55).

(iii) Time step limitations

We calculate Δt by solving (56).

(iv) Update of the geometrical quantities and the physical variables

We calculate (X_p,n+1,Y_p,n+1),τ_c,n+1,U_c,n+1,and E_c,n+1 from (35).

(v) Equation of state

The internal energy ε_c,n+1 is given by ε_c,n+1 = E_c,n+1 − ||U_c,n+1||²/2. Then,we can gain the pressure P_c,n+1 and the isentropic sound speed ac,n+1 from the equation of state. 4 Numerical results 4.1 Shock refraction problem

The computation domain is made up of two adjacent domains. The left border of the first domain is vertical,and the right boundary of the first domain is slanted at an angle of 600 relative to the vertical. The second domain is slanted but uniformly spaced at an angle of 60◦. A reflection boundary condition is imposed on the upper,lower,and right boundaries of the computation domain. The first domain is composed of 30 × 40 cells,and has an initial unit density. The second domain generates unstructured quadrilateral meshes,and has an initial density of 1.5. For the material,we utilize the equation of state of ideal gas (γ = 5/3). The left boundary of the first domain moves at the normal velocity u^* p = 1 from left toward right. The problem runs to the time t = 1.3 just before the shock starts to run off the second domain. We show results at this time. In Fig. 7,the initial mesh of a shock refraction problem is shown. In Figs. 8 and 9,the cells at t = 1.3 for the scheme of Maire^[14] and a high-order cell-centered conservative scheme are shown. In Figs. 10 and 11,the contour maps of the density at t = 1.3 are shown for the scheme of Maire^[14] and a high-order cell-centered conservative scheme. The importance of this problem is that we have both physical vorticity and feigned cell distortion^[3,10]. The cell in Fig. 8 is apparently distorted,especially near the lower part of the interface where the cell tends to tangling. The cell in Fig. 9 has not feigned cell distortion and hourglass type motion. The contour map of the density in Fig. 10 does not present physical vorticity. The contour map of the density in Fig. 11 presents physical vorticity. These results do not concede one to discriminate between the two calculations,but the high-order cell-centered conservative scheme is clearly the first rate from the viewpoint of robustness of the computed physical vorticity and the quality of the Lagrangian cell.

4.2 Saltzman’s shock tube

The computation domain is the quadrilateral (x,y) ∈ [0,1] × [0,0.1]. We use two types of initial cells. The first initial cell is gained by converting regular 50 × 5 cells with the mapping: X_tr = x + (0.1 − y) sin(xπ),Y_tr = y. Then,we generate unstructured quadrilateral meshes, unstructured hexagonal meshes,unstructured quadrilateral meshes,and unstructured triangular meshes in the triangle A1A2A3,the quadrilateral A2A4A5A3,the quadrilateral A4A6A7A5, and the triangle A7A6A8 of the domain. Here,the coordinates of A₁,A₂,A₃,A₄,A₅,A₆,A₇, and A₈ are (0.24,0.1),(0.3,0),(0.3,0.1),(0.33,0),(0.34,0.1),(0.35,0),(0.35,0.1),and (0.45,0), respectively. The second initial mesh is gained by converting regular 100×10 cells with the map ping: X_tr = x+(0.1−y) sin(xπ),Y_tr = y. Then,we generate unstructured quadrilateral meshes, unstructured hexagonal meshes,unstructured quadrilateral meshes,and unstructured triangular meshes in the triangle A₁A₂A₃,the quadrilateral A₂A₄A₅A₃,the quadrilateral A₄A₆A₇A₅, and the triangle A₇A₆A₈ of the domain. Here,the coordinates of A₁,A₂,A₃,A₄,A₅,A₆,A₇, and A₈ are (0.24,0.1),(0.3,0),(0.3,0.1),(0.33,0),(0.34,0.1),(0.35,0),(0.35,0.1),and (0.45,0), respectively. For the material,we utilize the equation of state of ideal gas (γ = 5/3). The initial condition is (ρc,0,Pc,0,Uc,0) = (1,0,0). The boundary condition is that there is a normal velocity u^* p = 1 (the inflow velocity) at x = 0. All the other boundaries are imposed reflection boundary conditions. In Fig. 12,the first initial mesh is shown. In Fig. 13,the second initial mesh is shown. The analytical solution is a shock wave moving at the speed D = 4/3. The correct answer for the density is 4.0 in the singly shocked region,10.0 in the doubly shocked region,and 20.0 in the triply shocked region. The locations of the shock front are 0.966 667 at t = 0.8 and 0.96 and t = 0.93,respectively. In Figs. 14 and 15,the first cell and a contour map of the density at t = 0.8 are shown for the high-order Larangian cell-centered conservative scheme. In Figs. 16 and 17,the second cell and a contour map of the density at t = 0.8 are shown for the high-order Larangian cell-centered conservative scheme. In Figs. 18 and 19,the first cell and a contour map of the density at t = 0.93 are shown for the high-order Larangian cell-centered conservative scheme. In Figs. 20 and 21,the second cell and a contour map of the density at t = 0.93 are shown for the high-order Larangian cell-centered conservative scheme. The main point of the result comparison is that a high-order cell-centered conservative scheme can remove the artificial cell distortion and hourglass type motion,and the results agree well with the precise solution. Finally,we record in Table 1 the mean value of the location of the shock front at t = 0.8. We use two meshes: the first cell and the second cell. The mean value of the location of the shock front at t = 0.8 with the first cell is 0.967 298 for the high-order Larangian cell-centered conservative scheme. The mean value of the location of the shock front at t = 0.8 with the second cell is 0.966 878 for the high-order Larangian cell-centered conservative scheme. The experimental order of the convergence is 1.58 for the high-order Larangian cell-centered conservative scheme.

5 Conclusions

We have constructed a high-order Lagrangian cell-centered conservative gas dynamics scheme on polygonal cells. The main features of the algorithm are that the high-order subcell forces are constructed by the piecewise constant pressures of the cell. The high-order Lagrangian cell-centered conservative gas dynamics scheme is constructed by means of the Taylor expansion of the fluxes in time-centered. The numerical results demonstrate that the high-order Lagrangian cell-centered conservative gas dynamics scheme is a genuinely two-dimensional highorder scheme.