A Runge Kutta Discontinuous Galerkin Method for Lagrangian Compressible Euler Equations in Two-Dimensions

. This paper presents a new Lagrangian type scheme for solving the Euler equations of compressible gas dynamics. In this new scheme the system of equations is discretized by Runge-Kutta Discontinuous Galerkin (RKDG) method, and the mesh moves with the ﬂuid ﬂow. The scheme is conservative for the mass, momentum and total energy and maintains second-order accuracy. The scheme avoids solving the geometrical part and has free parameters. Results of some numerical tests are presented to demonstrate the accuracy and the non-oscillatory property of the scheme


Introduction
In determining a numerical method for the multi-dimensional fluid flow, there are two typical choices. One is the Lagrangian framework in which the mesh is embedded in the fluid and moves with it. And the other is the Eulerian framework in which the mesh is treated as a fixed reference frame when the fluid moves. More generally, the grid points may be moved in some arbitrarily specified way that it is called the Arbitrary Lagrangian-Eulerian method (ALE). Most ALE algorithms consist of three steps. One is the Lagrangian step in which the solution and the grid are updated. The second step is a rezoning mesh in which the grid nodes are moved to a more optimal position and the third is remapping values in which the Lagrangian solution values are transferred to the new grid. In the numerical simulations of multi-material compressible fluid flows, the difficulty is how to handle the moving medium interface. Since the Lagrangian method can calculate fluid interface clearly, we present a new Lagrangian type scheme for solving the Euler equations of compressible gas dynamics in this paper.
An important point of constructing a Lagrangian discrete scheme is to decide where to locate degrees of freedom. It is generally separated into two kinds. One is the staggered scheme in which the velocity is defined at the nodes while the other variables are located inside the cells. This type of scheme was first introduced by von Neumann and Richtmyerfor in [2] for one-dimensional flows. The two-dimensional scheme was extended by Wilkins in [23] based on an internal energy formulation. The scheme was not conservative and admitted numerical spurious modes. In spite of these drawbacks, the scheme has been widely used for many years. Of course, many improvements have been made in order to solve the previous problems. Caramana and Shashkov constructed a staggered scheme which ensures the conservation of total energy in [1]. In [3], based on a mimetic finite difference scheme, Campbell and Shashkov improved the discretization of artificial viscosity so that the staggered Lagrangian scheme is an accurate and robust method. The other is the cell-centered scheme in which all variables are defined inside the cells. Here we concentrate our interests in the cell-centered schemes. This is because the staggered Lagrangian schemes use different control volumes for primary variables, it is very difficult to construct coherent high-order schemes for all these variables. On the contrary, the cell-centered Lagrangian schemes use only one control volume for primary variables, thus it is possible to construct coherent high-order schemes for all these variables.
In general, there are two kinds of the cell-centered Lagrangian schemes. The first is that the mesh moves with the flow velocity (e.g. [4][5][6][7][8]). In [4], Cheng and Shu developed a class of Lagrangian cell-centered schemes on quadrilateral meshes. Their schemes are based on the finite volume method and achieve a higher order accuracy by using the high-order ENO reconstruction. The schemes are conservative for the density, momentum and total energy, are essentially non-oscillatory, have free parameters, and can maintain formal high order accuracy both in space and time. But the ENO reconstruction requires the information from the surrounding cells. Therefore, this method has less compactness for high order schemes. In [6], Maire et al. developed a new Lagrangian cell-centered scheme. In their scheme the vertex velocities and numerical fluxes through the cell interfaces are not computed independently as usual but in a consistent manner. The scheme feature is the introduction of four pressures on each edge, two for each node on each side of the edge. In the limit of a one-dimensional flow computed by their two-dimensional solver, the scheme recovers the classical Godunov approximate Riemann solver. This scheme is first order accurate.
The second is that the mesh is fixed in the Lagrangian space (e.g. [9][10][11][12]19]). Després and Mazeran proposed a new and canonical way of writing the equations of gas dynamics in a fully Lagrangian form in [19]. They showed that the physical part is symmetrizable and that the weak hyperbolicity is due to shear contact discontinuities. Based on this formulation, a new conservative and entropy consistent Lagrangian scheme in two-dimension was obtained by using the finite volume method. The node velocity is computed in a coherent manner with the face fluxes. However, it appears that the scheme leads to a nodal velocity depending on the cell aspect ratio in the case of one-dimensional flows. In [9], Hui et al. introduced a unified coordinate system. In this system the flow variables are considered to be functions of time and of some permanent identification of pseudo-particles which move with velocity hU, where U is the velocity of fluid particles, and h is a parameter. This turns out to be a unified description, ranging from Eulerian when h = 0 to Lagrangian when h = 1, and the freedom in choosing h makes it possible to avoid the disadvantages of excessive diffusion across slip lines in the Eulerian description and of severe grid deformation in the Lagrangian description. But this method has to choose a proper h in every test case. In [10], Jia and Zhang solved a fully Lagrangian form of the gas dynamics equations which contains the physical part and the geometrical part. They discretized this system in the Lagrangian space by discontinuous Galerkin method using a spectral basis. The feature is that the vertex velocities and numerical fluxes are computed by using the method in [6]. In [12], Zhao and Yu also used discontinuous Galerkin method to solve the gas dynamics equations which contains the physical part and the geometrical part, and they used the Lax-Friedriches flux in their DG scheme. The schemes in [10] and [12] can maintain high-order accuracy both in space and time and has free parameter, but their system has eight equations and it is weakly hyperbolic system.
On the basis of previous studies, in this paper we give a new Lagrangian scheme for solving the compressible Euler equations. In this new scheme the system of equations is discretized by the RKDG finite element method. The mesh moves with the fluid flow and the scheme uses the Lax-Friedrichs and the HLLC numerical flux. In order to control spurious numerical oscillation near the discontinuity, we use a slope limiter in the scheme. The new scheme is conservative for the mass, momentum and total energy. The scheme avoids solving the geometrical part, and has no parameters to be tuned for individual test cases. Several two-dimensional numerical examples show that our scheme can achieve uniformly second-order accuracy on moving and distorted Lagrangian meshes. Moreover, it is essentially non-oscillatory.
The paper is organized as follows. In Section 2 we recall the gas dynamics equations written in the Lagrangian form, and give their weak form. In Section 3 we use the RKDG method to solve this weak form. In Section 4 we validate our new scheme with several test cases which demonstrate the non-oscillatory property and the accuracy. Finally, we give concluding remarks in Section 5.

Derivation of the Euler equations in the Lagrangian formalism
In this section we recall briefly how compressible Euler equations are written in the Lagrangian formalism which follow the approach developed in [6], and give its weak form.
Let D be a region filled with an inviscid fluid. Considering a fluid particle, its initial position is the point M, and at time t >0 it is moving through the point m. The coordinates of point M are denoted by (X,Y) and are named Lagrangian coordinates. The coordinates of point m are denoted by (x,y) and are named Eulerian coordinates. The Eulerian coordinates are obtained from the trajectory equations where V = (u,v) is the fluid velocity. If the velocity field is smooth enough, there exists a unique smooth solution (x(t,X,Y),y(t,X,Y)). The map is defined as where (x,y) is the unique solution of (2.1). With fixed time t > 0, this map advances each fluid particle from its position at time t = 0 to its position at time t.
If Ω denotes a region in D, then χ t (Ω) = W is the volume Ω moving with the fluid. We assume that for each t > 0, χ t is invertible. The Jacobian of the map is Since J(X,Y,0)=1, and χ t is invertible, for each t>0, we have J(X,Y,t)>0. We know that the time differentiation of (2.3) is The system of compressible Euler equations (see [20]) is where ρ is the density, p is the press, and E is the specific total energy. We denote the specific internal energy by ε = E− 1 2 V 2 . The thermodynamic closure of (2.7) is given by the equation of state p ≡ p(ρ,ε). Now, by using (2.6) for the conservative variables ϕ = ρ,ρV,ρE and substituting them into (2.7), we get the compressible Euler equations in Lagrangian form as We notice that (2.8) is only a semi-Lagrangian formula since the gradient and divergence operate on variables which depend on Eulerian coordinates. In order to write (2.8) in a full Lagrangian way, one has to express the gradient and divergence operators with Lagrangian coordinates using the map χ t .
In order to use the discontinuous Galerkin (DG) method to perform the space discretization in the next section, we give the weak form of system (2.8). That is, both sides of Eq. (2.8) are multiplied by the function ψ(X,Y) ∈ L 2 (Ω) and integrated on the Lagrangian domain Ω: We set ψ·χ −1 t (x,y) = ϕ(x,y). By knowing that χ t (Ω) = W and J ·dΩ = dW, we can get where V =(u,v) is the velocity. We call that (2.10) is the weak form of (2.8). We notice that d dt ψ(X,Y) = 0, Eq. (2.10) are obtained by performing a mapping transformation on (2.9). We set ϕ = 1 and use Green's formula, then (2.10) becomes where n is the unit outward normal vector of ∂W. This system (2.11) represents the conservation of mass, momentum and energy so that the system (2.10) also satisfies the conservation of mass, momentum and energy.

RKDG method 3.1 Spatial discretization
In this section, we present a spatial discretization of (2.10) by using the DG method. Note that W t=0 = Ω, Ω is discretized into M× N computational cells. As time progresses, we use the velocity at the vertex to get the new fluid region W and the new parti- We define a finite-element space consisting of piecewise polynomials where P k (W ij ) denotes a set of polynomials of degree up to k defined on the cell W ij . We use the DG method to solve (2.10) as follows: find the approximate function Then (3.1) can be rewritten as By using Green's formula we get where n=(n x ,n y ) is the unit outward normal vector of ∂W ij ,(F, G)·n is the numerical flux, which is a single-valued function defined at the cell interfaces. Here we set H = (F, G)·n.
In this paper, we use the following two typical numerical fluxes: (1) The L-F (Lax-Friedrichs) flux where U ± h are the values of U h inside and outside the cell W ij , and α ij is taken as an upper bound for the eigenvalues of the Jacobian. Here in the Lagrangian formulation, We refer to [14] for the details of the HLLC flux. This flux has been successfully used to compute the compressible Euler equations in the Lagrangian framework [4].

The basis functions
Let the Taylor basis be a set of basis functions of V k h . For the P 1 case, assuming that u h is a component of U h , we use the following expression for the approximate solution where the basis functions are Here (x c ,y c ) is the centroid of the cell W ij . ∆x = 0.5(x max −x min ), ∆y = 0.5(y max −y min ), and x max , x min , y max and y min are the maximum and minimum coordinates in the cell W ij in x-, and y-direction, respectively. Because for all test functions 3) denotes the mass matrix, and L is a discrete spatial operator. Notice that P k (W ij ) and B lij (l = 1,2,3) are defined in the cell W ij . As time progresses, W ij is updated, so P k (W ij ) and B lij (l = 1,2,3) should be updated at every time step. The detail is described in Section 3.4.

The determination of the vertex velocity
In the Lagrangian formula, the grid moves with the fluid velocity, thus we would need to know the velocity at the vertex to move the grid. Since the velocity is a deduced quantity, we would need to obtain it from the conserved variables. In the following we describe how to determine the vertex velocity in our scheme, which follows the approach in [4].
Considering a vertex (i,j) shared by four edges, we give these four edges a serial number k = 1,2,3,4, and define the direction of each edge to be the direction of the incremental index i or j, for example the direction of the edge with two endpoints (i−1,j) and We then split the left and right velocities into normal and tangential components along the edge k. Let (n k x ,n k y ) be the clockwise unit normal of the edge k, w k− t and w k+ t be their tangential components, and w k− n and w k+ n be their normal components. Then we define the tangential velocity of the vertex (i,j) along the edge k As to the normal velocity, we define it as the Roe average of the normal velocities from its two sides Thus by the formulas (3.6) and (3.7), we can get four x-velocities and y-velocities at the vertex (i,j) which have the following form Finally, the velocity at the vertex (i,j) is obtained as follows

Time discretization
The time marching for the semi-discrete scheme (3.5) is implemented by a class of Runge-Kutta (RK) type methods [21]. Since the mesh changes with the time advancing in the Lagrangian simulation, the velocity at the position of each vertex and the size of each cell needs to be updated at each RK stage. In this section the variables with the superscripts n and n+1 represent the values of the corresponding variables at the nth and n+1th time steps, respectively. For the first-order spatial discretization, the corresponding time discretization uses the first-order forward Euler step. Now the basis function is B ij = 1. So where S ij is the area of cell W ij . We know the value of each variable at the nth time step. By using the vertex velocity we get Therefore, the first-order time discretization of (3.5) is M n+1 ij ·u n+1 hij = M n ij ·u n hij +∆t n L(U n h ). (3.10) For the second-order spatial discretization, the corresponding time discretization uses the second-order RK method. We know the value of each variable at the nth time step. The time marching is as follows.
Step 1: Using the method in Section 3.3, we get the vertex velocity (u (1) ij ,v (1) ij ), then Through the four vertices (x (1) i,j ,y (1) and ∆y (1) . Then Step 2: Using variables obtained from Step 1 and the method in Section 3.3, we get the vertex velocity (u n ij ,v n ij ), then Through the four vertices ( , ∆x n+1 and ∆y n+1 . Then The time step ∆t n is chosen as follows where ∆l n ij is the shortest edge length of the cell W ij , and c n ij is the sound speed within this cell. The Courant number λ in the following tests is set to be 0.1 unless otherwise stated.

Boundary conditions
Usually boundary conditions are expected to provide numerical fluxes at the boundaries. We add some fictitious cells adjacent to the boundary real cells, and give the values of variables in the fictitious cells according to the boundary conditions, then we can compute the boundary fluxes as done for the interior cells. In this section, we prescribe fictitious data values in the fictitious cells, and the variables with the superscripts of f and r represent the values of the corresponding variables in fictitious cells and the real cells, respectively. Here we only consider two types of boundaries: prescribed normal velocity and transmissive boundaries.

Prescribed normal velocity.
Let w * n be the value of the prescribed normal velocity along the boundary edge. Let (w f n ,w f t ),(w r n ,w r t ) be the normal and tangential velocity along the boundary edge for the fictitious cell and the real cell, respectively. Then we set Thus we can get the x-velocity and the y-velocity: where (n x ,n y ) is the unit normal vector of the fictitious cell. We also set where ρ f and ρ r are the density, e f and e r are the specific internal energy. Then we can get the specific total energy Thus we get the conserved variables (ρ f ,ρ f u f x ,ρ f u f y ,ρ f E f ). Transmissive boundaries. For a transmissive boundary we give the boundary condition as follows:

The local slope limiting
In the case of piecewise-constant approximations, the artificial viscosity is introduced by the numerical flux which is enough to keep the stability. But for high-order piecewise polynomial approximations, the influence of the numerical fluxes can not be enough to guarantee the absence of spurious oscillations of the scheme. To enhance the stability of the scheme and eliminate possible spurious oscillations in the approximate solution, we introduce a local slope limiting developed in [21]. The slope limiting is performed on u ij2 and u ij3 in (3.4). By using the differences of the means. For a scalar equation, u ij2 would be limited bỹ (3.12) where the functionm is the TVB corrected minmod function defined bȳ m(a 1 ,a 2 ,··· ,a n ) = a 1 , if |a 1 | ≤ H ·∆x 2 , m(a 1 ,a 2 ,··· ,a n ), else, (3.13) with the minmod function m(a 1 ,a 2 ,··· ,a n ) = s·min|a i |, For an estimate of the TVB constant H in (3.13), see [15]. Similarly, u ij3 is limited bỹ (3.14) just with a change of ∆x to ∆y in (3.13).
For the system, when we obtain the physical variables and the moving control volume W n+1 at the time t = t n+1 , we perform the limiting in the local characteristic variables at this fixed time t = t n+1 and this fixed fluid region W n+1 (see [21]). Here we just give a simple description.
The Euler equations (2.7) can be rewritten as To limit the vector U ij2 in the element W n+1 ij , we proceed as follows: (1) Find the matrix R and its inverse R −1 which diagonalize the Jacobian evaluated at the mean in the element W n+1 ij in the x-direction, where Λ is a diagonal matrix containing the eigenvalues of the Jacobian.
(2) Transform all quantities needed for the limiting, i.e., the three vectors U ij2 ,Ū i+1,j −Ū i,j , andŪ i,j −Ū i−1,j , to the characteristic fields. This is achieved by left-multiplying these three vectors by R −1 .
(3) Apply the scalar limiter (3.12) to each of the components of the transformed vectors.
(4). The result is transformed back to the original space by left multiplying R on the left.

Numerical results
It is much more difficult to simulate a 2D problem than to simulate a 1D one in the Lagrangian framework, mainly because of the mesh distortion in multi-dimensions. In this section, we present several tests in order to validate our numerical scheme. We use the L-F flux and the HLLC flux in our scheme. Their results are similar, so we mainly show the results obtained by using L-F flux.

Accuracy test
First we compute two problems to check the accuracy of our schemes.
The domain is taken as [0,1]×[0,1], and periodic boundary conditions in both directions are used. The exact solution is The errors and the convergence order of our first, second-order Lagrangian schemes with L-F flux at t = 0.1 are listed in Tables 1 and 2, respectively. Through the convergence results, we can see the desired first and second-order accuracy.  Example 4.2 (The two-dimensional periodic vortex problem). The periodic vortex problem is described as follows: The mean flow is ρ=1, p=1 and (u,v)=(1,1) (diagonal flow). We add to this mean flow an isentropic vortex perturbation in (u,v) and the temperature T = p ρ , no perturbation in the entropy S = p ρ γ , that is  where (x,ȳ) = (x−5,y−5), r 2 =x 2 +ȳ 2 , and the vortex strength is ε = 5.
The domain is taken as [0,10]×[0,10], and periodic boundary conditions in both directions are used. It can be readily verified that the Euler equations with the above initial conditions admit an exact solution which is convected with speed (1,1) in the diagonal direction. The errors and the convergence order of our first, second-order Lagrangian schemes with L-F flux at t = 1.0 are listed in Tables 3 and 4, respectively. The initial mesh and the mesh at t = 1.0 are displayed in Fig. 1. Through the convergence results and Fig. 1, we can see the desired first and second-order accuracy on moving and distorted Lagrangian meshes.

Non-oscillatory tests
We compute five problems to check the non-oscillatory property of our scheme. Example 4.3 (Two-dimensional Sod's shock tube problem [12]). This problem can assess the reasonableness and non-oscillatory property of our schemes. Its initial condition is as follows The domain is taken as [0,1]×[0,1]. The problem can be interpreted as the two-dimension form of one-dimensional Sod's shock tube problem. Let nx and ny be the number of cells in the x-and y-directions. We set nx = 200 and ny = 10. In Fig. 2, we present the results by using the L-F flux at t = 0.2, y = 0.5, and we plot the piece-wise straight line connecting the two ends of each cell in x-direction for the density and pressure results of the second scheme. From the results we can see that the second-order scheme satisfy the non-oscillatory property. Example 4.4 (The Saltzman problem [18]). This is a well known difficult test to validate the robustness of a Lagrangian scheme when the mesh is not aligned with the fluid flow. The problem consists of a rectangular box whose left end is a piston. The piston moves into the box with a constant velocity of 1.0. The one-dimensional symmetry is broken by the initial mesh composed of 100×10 non-uniform cells defined as follows: x(i, j) = (i−1)∆x+(11− j)sin(0.01π(i−1))∆y, y(i, j) = (j−1)∆y, where ∆x = ∆y = 0.01. The initial mesh is displayed in Fig. 3. Note that the initial mesh is deliberately distorted to set it as a more demanding test case. The working fluid is described by an ideal gas with γ = 5/3. The initial conditions involve a stationary gas with a unity density and an internal energy of 10  wave hits the face x = 1. We can observe that the numerical result of our second-order scheme is well except for the region near the up and bottom wall boundaries where the results are affected by the boundary conditions. Example 4.5 (The Sedov blast wave problem [12]). We consider the propagation of a high intensity cylindrical shock wave generated by a strong explosion; see, for instance, [12]. The initial density is unity and the initial velocity is zero. The internal energy is zero except in the regional center where it has a value of 1. The domain is taken as [−1.21,1.21]×[−1.21,1.21], and the mesh is of 121×121 cells. All the boundary conditions are wall conditions. The analytical solution gives a shock at radius unity at time unity with a peak density of 6. Due to the symmetry of the problem, we only give the grid and the density contours of the first quadrant at top and middle of Fig. 5, respectively. But we display the result of the density in all the cells with respect to the cell center radius at the bottom of Fig. 5. In Fig. 6 we give the density profiles of the x axis of our secondorder scheme and the same order scheme in [12]. From the direct comparison with the exact solution, we can clearly see that our second-order scheme is more precise than the scheme in [12].
Example 4.6 (Explosion problem [17]). The domain is taken as [0,2]×[0,2]. Its initial condition is as follows All the boundary conditions are wall conditions. We remark that in initializing the explosion problems, we modify the initial data on quadrilateral cells cutting the initial discontinuity, by assigning modified area-weighted values to the appropriate cells at the initial time t = 0. This procedure avoids forming small amplitude waves created at early times by the staircase data. Fig. 7 shows the results of our scheme by using the L-F flux with 200×200 initially uniform cells at t = 0.25. The up pictures in Fig. 9 show the results of the WAF (Weighted Average Flux) scheme developed in [17] with 200×200. In order to obtain the convergence of our scheme, we compute this problem by using the WAF scheme with 1000×1000 cells as the reference "exact" solution. Fig. 8 shows a comparison between the reference "exact" solution and our solutions with the three grids 100×100, 200×200 and 300×300 along the radial line which towards the origin (1,1). We can observe the convergence of the numerical solutions toward the reference "exact" solution.  Figure 6: The density profiles along the x axis of our second-order scheme and that of the same order scheme in [12]. All the boundary conditions are wall conditions. We modify the initial data on quadrilateral cells cutting the initial discontinuity same as the method in Example 4.6. This problem is created by reversing the initial data of (4.1). Now shock focusing takes place  as part of the solution. Fig. 10 shows the results of our scheme by using the L-F flux with 200×200 initially uniform cells at t = 0.25. The down pictures in Fig. 9 show the results of the WAF(Weighted Average Flux) scheme developed in [17] with 1000×1000 initially uniform cells at t = 0.25 as the reference "exact" solution. Fig. 11 shows a comparison between the reference "exact" solution and our solutions with the three grids 100×100, 200×200 and 300×300 along the radial line which towards the origin (1,1). We can ob-serve the convergence of the numerical solutions toward the reference "exact" solution, and the numerical result of the second-order scheme is obviously much better than that of the first-order scheme.

Concluding remarks
In this paper we have proposed a new Lagrangian type scheme to solve the compressible Euler equations. In this new scheme the system of equations is discretized by the RKDG method, and the mesh moves with the fluid flow. The Lax-Friedrichs numerical flux and the HLLC numerical flux are used in our scheme. A slope limiter is used to control spurious numerical oscillations near discontinuity. This new scheme is conservative for the mass, momentum and total energy. It can maintain second-order accuracy, avoids solving the geometrical part, and has free parameter for individual test cases. Several two-dimensional numerical examples have been presented to demonstrate the good performance of the scheme.