Matter flow method for alleviating checkerboard oscillations in triangular mesh SGH Lagrangian simulation
Li Zhao, Bo Xiao, Ganghua Wang, Haibo Zhao, Jinsong Bai, Chunsheng Feng, Shi Shu
JJournal of Computational Physics (2020)
Contents lists available at ScienceDirect
Journal of Computational Physics
Matter flow method for alleviating checkerboard oscillations in triangular meshSGH Lagrangian simulation
Li Zhao a , Bo Xiao b, ∗ , Ganghua Wang b , Haibo Zhao c , Jinsong Bai b , Chunsheng Feng a , Shi Shu a a School of Mathematics and Computational Science, Xiangtan University, Xiangtan 411105, China b Institute of Fluid Physics, CAEP, Mianyang 621999, China c Department of Mechanics and Aerospace Engineering, Southern University of Science and Technology, Shenzhen 518055, China
A R T I C L E I N F O
Article history : Keywords:
Compressible hydrodynamicsLagrangian methodcheckerboard oscillationMatter flow methodParallel computation
A B S T R A C TWhen the SGH Lagrangian based on triangle mesh is used to simulate com-pressible hydrodynamics, because of the sti ff ness of triangular mesh, the prob-lem of physical quantity cell-to-cell spatial oscillation (also called ”checker-board oscillation”) is easy to occur. A matter flow method is proposed toalleviate the oscillation of physical quantities caused by triangular sti ff ness.The basic idea of this method is to attribute the sti ff ness of triangle to the factthat the edges of triangle mesh can not do bending motion, and to compensatethe e ff ect of triangle edge bending motion by means of matter flow. Threee ff ects are considered in our matter flow method: (1) transport of the mass,momentum and energy carried by the moving matter; (2) the work done onthe element, since the flow of matter changes the specific volume of the gridelement; (3) the e ff ect of matter flow on the strain rate in the element. Numer-ical experiments show that the proposed matter flow method can e ff ectivelyalleviate the spatial oscillation of physical quantities.c (cid:13)
1. Introduction
The motion of compressible multi-material large deformation fluid is a common hydrodynamic process in thefields of high energy density such as detonation, inertial confinement fusion, superhigh velocity collision, astro-physics. It is also a di ffi cult point in hydrodynamic numerical simulation. Currently, among the main technicalschemes to simulate the motion of compressible multi-material large deformation fluid are Euler method [1–17] andArbitrary Lagrange-Euler (ALE) method [18–24]. The ALE method here usually refers to the ALE method that al-lows the interface mesh to move across di ff erent matters, also known as MMALE (Multi-Material ALE). Becauseof the cross-matter motion of the grid in the ALE method, there will be mixed matter grid elements similar to theEulerian method, and the mixing can cause the dispersion of the material interface. In order to control the interface ∗ e-mail: [email protected]; a r X i v : . [ phy s i c s . c o m p - ph ] A ug Li Zhao et al. / Journal of Computational Physics (2020) dispersion, it is necessary to introduce interface reconstruction (such as VOF [2, 12, 20] and MOF [19, 20, 22, 23])into Euler method and ALE method.At present, Lagrangian method [25–36], despite its advantage in edge-tracking for multi-material fluids, has notbecome the mainstream method of compressible multi-material large deformation hydrodynamic simulation. Its mainreasons include mesh distortion, physical quantity oscillation (mainly in two-dimensional triangular mesh and three-dimensional tetrahedral mesh), and it is not easy to deal with the interface topology changes caused by materialcollision or fracture. The oscillation of physical quantities is the focus of this paper.In Lagrangian hydrodynamic simulation, when triangular or tetrahedral meshes are used, it is easy to appear thecell-to-cell oscillation distribution of physical quantities between mesh elements, that is, the problem of checkerboardoscillation of physical quantities. The reason can be attributed to the mesh sti ff ness of triangular or tetrahedral meshes.For the nonphysical checkerboard oscillation problem, some research work [35, 36] has been done at present. In 2012,G.Scovazzi [35] took the lead in discussing the use of ”Flux” to alleviate sti ff ness. Starting from variable multiscaleanaysis (VMS), scovazzi makes a linear approximation of the mesh motion on a finer scale and transforms it intoflux on the edge of a large-scale mesh. After a series of approximations, Scovazzi discarded many complex terms,and finally retained an energy di ff usion term proportional to the pressure gradient in flux. Scovazzi’s ”Flux” methodcan alleviate the oscillation in some typical shock wave problems. However, this method also causes some other nonphysical e ff ects, such as the density increases instead of decreasing at the wall heating. This may be due to the oversimplification of Scovazzi’s ”Flux” term. In 2015, N. R. Morgan [36] also discussed a method of using ”Flux” toalleviate the sti ff ness. The author thinks that in the Point-Centered Lagrangian hydrodynamics (PCH) discretization,sti ff ness originates from the volume error, and then proposes a matter flow term to correct the error. The so-calledmatter flow means that matter is allowed to be transported from one control body to another. In addition to mass, theenergy and momentum carried by matter will be transported along with the matter.Similar to the method in the literature [35, 36], this paper also constructs a ”Flux” method to alleviate the physicalquantity oscillation caused by triangle sti ff ness in two-dimensional SGH Lagrange simulation. The core idea of thismethod is to attribute the sti ff ness of triangular mesh to the fact that there is no proper bending on the edge of the meshelement, and a ”Flux” is constructed to replace the edge bending e ff ect. ”Flux” method in this paper considers a varietyof e ff ects, including the mass, momentum, and energy transport caused by the transport of matter between elements(similar to that done in the Morgan), the energy transport between elements due to the work done by the specificvolume change of the elements (this is a bit like that done in the Scovazzi), and the influence of material transport onthe strain rate of the element. To facilitate the application of matter flow method in parallel hydrodynamic program,this paper also designs a parallel implementation scheme of matter flow SGH Lagrange method based on OpenMP[37].The following chapters are arranged as follows: Section 2 introduces the compressible hydrodynamic equationsand discrete format. Section 3 presents a matter compensation flow method. Section 4 discusses the parallel imple-mentation algorithm based on shared memory. Section 5 gives numerical examples and analysis. Section 6 summa-rizes and discusses the research work of this paper.
2. Compressible hydrodynamic equations and discrete format
Consider the following two dimensional compressible hydrodynamic equations in this article. Equations (1) - (4)are mass equation, momentum equation, internal energy equation and equation of state,respectively.1 ρ d ρ dt + ∇ · u = , (1) ρ d u dt = −∇ ( p + q ) , (2) ρ dedt = − ( p + q ) ∇ · u , (3) p = p ( ρ, e ) . (4)where ρ , e , u are the density, specific internal energy and velocity of the fluid, p is the pressure and q is the artificialviscosity. The specific forms of di ff erential operators ∇· and ∇ are respectively ∂∂ x + ∂∂ y and ( ∂∂ x , ∂∂ y ) T , and ddt is theLagrangian time derivative. i Zhao et al. / Journal of Computational Physics (2020) 3
In the example of this paper, the equation of state (4) is taken as the ideal gas equation of state, which is expressedas follows: p ( ρ, e ) = ( γ − ρ e (5)where γ is the gas adiabatic index. In computational fluid dynamics, according to locations (in elements or on grid points) on which physical quan-tities in the Lagrangian method are defined, it can usually be divided into Staggered-Grid Hydrodynamics (SGH)method, Cell-Centered Hydrodynamics (CCH) method and Point-Centered Hydrodynamics (PCH) method. In theSGH method, the pressure, density and internal energy are defined in the center of the element, and the velocity andkinetic energy are defined on the nodes. In this paper, the SGH Lagrangian finite volume method is used to discretizethe control equations, and the physical quantities defined in cells are treated as piece-wise constant. The discretecontrol volume is shown in the Fig. 1. (a) triangular element c (b) node p Fig. 1. Control volume diagram
For the SGH Lagrangian finite volume method, the discrete form of compressible hydrodynamics equations (1)-(3) are expressed as follows: m c = constant , (6) m p = (cid:88) c ∈ T ( p ) m c , (7) u n + p = u np + (cid:80) c ∈ T ( p ) f npc m p ∆ t , (8) x n + p = x np + u np ∆ t + (cid:80) c ∈ T ( p ) f npc m p ( ∆ t ) , (9) E n + c = E nc − (cid:88) p ∈ P ( c ) f npc · ( x n + p − x np ) . (10)where m c represents the mass of the element c , m p denotes the mass of the node p , u np and u n + p respectively representthe velocity of the node p at the moments t n and t n + , x np and x n + p respectively represent the position of the nodes p at t n and t n + moments, E nc and E n + c respectively indicate the internal energy of element c at t n and t n + , T ( p ) representsall the element sets containing nodes p , P ( c ) represents all the node sets in element c , f npc represents the force of theelement c on the node p , its expression is as follows: f pc = ( p c + q c ) − y pc − y pc x pc − x pc , where p c , q c denotes the pressure and viscous in the element c . The expression of q c is as follows: q c = − c visc ρ c ˙ v , (11) Li Zhao et al. / Journal of Computational Physics (2020) where c visc is the viscosity coe ffi cient, ρ c is the density in the grid element c , and ˙ v is the relative change rate ofvolume. The viscosity coe ffi cient is selected as follows: c visc = k max (cid:8) − vh , v s h (cid:9) . (12)where h is the maximum value of the three sides of the triangular element, and v s is the sound velocity and k is anadjustable factor, which is taken as 2.0 in this paper. To retain time stability, the time step ∆ t in the evolution of Lagrange hydrodynamics needs to satisfy the stabilitycondition. In this paper, the time step is selected according to the following conditions. ∆ t = min c ∈ C (cid:110) ∆ t cs , ∆ t cv , ∆ t cf v , ∆ t cf a (cid:111) . where C is the set of grid elements, ∆ t cs , ∆ t cv , ∆ t cf v , ∆ t cf a for each element are as follows.(1) Time step determined by sound velocity: ∆ t cs = C sa f e · h min v s . where C sa f e = .
05 is an adjustable safe factor, which is taken to be 0.05 all through the following, h min is theminimum height of the triangle, v s is the sound velocity of the triangle.(2) Time step determined by viscosity: ∆ t cv = C sa f e · h min c visc . where c visc is the viscosity coe ffi cient of the triangle.(3) The time step determined by the velocity of matter flow in a triangle: ∆ t cf v = C sa f e · min (cid:110) h | u f low | , h | u f low | , h | u f low | (cid:111) . where h , h , h is the three heights of the triangle, u f low , u f low , u f low are the matter flow velocities on the threesides of the triangle (see section 3).(4) The time step determined by the acceleration of triangular matter flow: ∆ t cf a = C sa f e · min (cid:110) (cid:115) h | a f low | , (cid:115) h | a f low | , (cid:115) h | a f low | (cid:111) . where a f low , a f low , a f low is the acceleration of matter flow on the three sides of the triangle (see section 3).
3. Matter compensation flow method
When the SGH Lagrangian method based on section 2.2 is used to simulate the motion of compressible hydro-dynamics, it is easy to appear the phenomenon of physical quantity cell-to-cell oscillation caused by the sti ff ness oftriangular mesh. Fig. 2 gives an intuitive description of the physical quantity spatial oscillation caused by the sti ff nessof the triangular mesh: suppose Fig. 2(a) the quadrilateral mesh abcd be filled with fluid, nodes b and c are fixed,nodes a and d move in the direction of the arrow. Under these conditions, fluid density in the quadrilateral abcd willdecrease as the mesh area increases. While in Fig. 2(b), a quadrilateral grid abcd is divided into four triangular grids.The rest of the conditions do not change. With the movement of nodes a and d , the length of the edge ad decreases,and the area of the triangle ade becomes smaller. Finally, the density and pressure in the triangle ade increase sig-nificantly higher than that in the adjacent triangle. As a result, the cell-to-cell oscillation phenomenon of physicalquantities appears, as shown in Fig. 2(c).Inspired by the above analysis, we attribute the sti ff ness of a triangle to the fact that the edges of the trianglecannot do bending motion. Taking Fig. 2(c) as an example, if the edge ae and de of the triangle ade can do bending i Zhao et al. / Journal of Computational Physics (2020) 5(a) Quadrilateral computing grid (b) Triangular computing grid (c) Non-physical oscillation phenomena
Fig. 2. Oscillation of physical quantities in triangular meshes [32] motion, the edge ae and de will bend outward with the increase of pressure in the triangle ade , which will compensatefor the decrease of the area of the curved triangle ade . Thus, the oscillation of physical quantity is alleviated (Thisalso tells us that in principle, if we adopt a Lagrangian method [18] which allows the grid to bend, the oscillation ofphysical quantities can be alleviated). From this point of view, this paper proposes a method of matter compensationflow to approximate the e ff ect of triangular side bending motion, so as to alleviate checkerboard oscillation. The basicidea of this method is shown in Fig. 3. In the Fig. 3(a), let triangle bd f pressure be greater than triangle ab f . If thesides of a triangle can bend, under pressure di ff erential, the edge b f will become a curved bg f . Because in the usualLagrangian simulation, the mesh is actually not allowed to bend, so we can consider using the material compensationflow between cells to replace the e ff ect of the edge bending motion (Fig. 3(a) the area of the shadow part determinesthe amount of matter flow). Because the curved edge b f g is not easy to obtain, in order to calculate conveniently, theshadow part is approximated to a triangle, as shown in Fig. 3(b). The triangle can be considered to be formed by themidpoint h of the edge b f moving to node g under the action of pressure di ff erence. (a) Triangle mesh edge bending (b) Approximate the curved area with a triangle. Fig. 3. Triangle mesh edge bending diagram
There are three e ff ects related to the matter compensation flow. The first e ff ect is that the mass, momentum andenergy carried by the matter are transferred from the grid cell to the adjacent grid cell in the process of flow. Thesecond e ff ect is that the specific volume of the grid element is changed due to the ”squeezing in” and ”extrusion” ofthe matter from the grid cell, which will produce work e ff ect on the original matter. The third e ff ect is that the volumestrain rate of the grid element is also a ff ected by the matter flow, which leads to the change of the viscous stress of thegrid element, which will eventually a ff ect the evolution of the internal energy. Based on the above three e ff ects, thesteps of the matter compensation flow method are designed as follows (Fig. 4 gives an illustration and defines someof the symbols for the description of the steps). Step 1
Using the accelerations of node B and node D to calculate the component of the average acceleration in thenormal direction of the midpoint I of the edge BD :¯ a I = s n · a + a . (13)where s n is the out of unit normal vector of the edge BD in the element K , a is the acceleration of node B , a is the acceleration of node D . Li Zhao et al. / Journal of Computational Physics (2020)(a) ∆ x ≥ ∆ x < Fig. 4. Illustration of the matter compensation flow method.
Step 2
The acceleration of node I is also calculated by the pressure di ff erence between triangular element K and K Nb : a I = p K − p K Nb m I , (14)where p K and p K Nb are the pressures of elements K and K Nb respectively, and the mass of node I is: m I = m K + m K Nb . m K is the mass of element K and m K Nb is the mass of K Nb . Step 3
The di ff erence between a I and ¯ a I determines the acceleration of matter flow: a f low = a I − ¯ a I , (15) a f low = a f low s n . where the direction of a f low is the normal direction of the edge BD . Step 4
T hematter f lowaccelerationa f low determines the imaginary movement of the node I : ∆ x = u nf low ∆ t + a nf low ( ∆ t ) , u n + f low = ( u nf low + a nf low ∆ t )(1 − C diss ) . (16)where u nf low is the size of the matter flow velocity of the t n moment, a nf low is the amount of matter flow acceler-ation size at the t n moment, ∆ t is the time step, and C diss represents the artificial dissipation factor of the flowvelocity, which is proportional to the viscous coe ffi cient: C diss = max (cid:110) S K c Kvisc , S K Nb c K Nb visc (cid:111) here, c Kvisc and S K are the viscosity coe ffi cient and area of the element K , and c K Nb visc and S K Nb are the viscositycoe ffi cient and area of the element K Nb , respectively. Remark 1.
The direction of ∆ x and u n + f low is the normal direction of edge BD. Step 5
Mass compensation, according to conservation of mass: m (cid:48) K = m K − ∆ Mm (cid:48) K Nb = m K Nb + ∆ M (17)where ∆ M is the mass carried by the matter flow, and its calculation formula is ∆ M = (cid:40) ∆ xL ρ K , ∆ x ≥ ∆ xL ρ K Nb , ∆ x < i Zhao et al. / Journal of Computational Physics (2020) 7 L is the length of edge BD , and the ρ K and ρ K Nb are the density of K and K Nb , respectively.According to the definition of node mass, the mass of nodes A and C is modified as m (cid:48) A = m A − ∆ mm (cid:48) C = m C + ∆ m (19)where ∆ m : = ∆ M , the mass of nodes B and D does not change. Step 6
The change of internal energy should not only consider the ∆ E carried by the matter flow, but also the extrawork W extra caused by the change of element volume. E (cid:48) K = E K − ∆ E − W extra E (cid:48) K Nb = E K Nb + ∆ E + W extra (20)where E K and E K Nb are the internal energies of element K and K Nb respectively, and the calculation formula of ∆ E is ∆ E = ∆ Mm K E K , ∆ M ≥ ∆ Mm KNb E K Nb , ∆ M < W extra formula (See Remark 2 for details) is W extra = (cid:104) ( p K + q K ) ∆ M ρ K + ( p K Nb + q K Nb ) ∆ M ρ K Nb (cid:105) . (22)here q K is the viscous force of element K and q K Nb is the viscous force of element K Nb Remark 2.
The extra work done by the movement of matter to cause changes in element volume is: W extra = ( p + q ) ∆ v , where p and q are the pressure and viscous force of the element. ∆ v is the volume change of the element, andits calculation formula is: ∆ v = ∆ Mm v = ∆ M ρ , So the extra work can be written as: W extra = ( p + q ) ∆ M ρ , To maintain conservation of energy, the extra work is averaged over the element K and K Nb , and then writtenas: W extra =
12 [( p K + q K ) ∆ M ρ K + ( p K Nb + q K Nb ) ∆ M ρ K Nb ] . Step 7
The change of volume relative rate caused by matter flow is as follows: ∆ ˙ v = ∆ Mm K ∆ t , ∆ M ≥ ∆ Mm KNb ∆ t , ∆ M < ff ects the value of viscous by (11), and then a ff ects the change of internalenergy and other physical quantities. Step 8
Velocity compensation: u (cid:48) A = p (cid:48) A m (cid:48) A u (cid:48) C = p (cid:48) C m (cid:48) C (24)where p (cid:48) A = p A − ∆ p , p (cid:48) C = p C + ∆ p , p A and p C are the momentum of nodes A and C , respectively. ∆ p satisfiesthe following optimization problems: ∆ p = arg min p ∈ R (cid:8) | p − λ | (cid:9) a (cid:104) p , p (cid:105) + (cid:104) b , p (cid:105) + c = . (25) Li Zhao et al. / Journal of Computational Physics (2020) here λ = ∆ m u B + u D , a = m (cid:48) A + m (cid:48) C > b = m (cid:48) C p C − m (cid:48) A p A , c = ( m (cid:48) A − m A ) (cid:104) p A , p A (cid:105) + ( m (cid:48) C − m C ) (cid:104) p C , p C (cid:105) . The derivation of formula (25) is given below . Take Fig. 4(a) as an example (Fig. 4(b) has the same result).We suppose that node A transports the momentum of ∆ p (to be solved) to node C , then the conservation ofmomentum is used p (cid:48) A = p A − ∆ p , p (cid:48) C = p C + ∆ p . (26)According to conservation of kinetic energy12 m A (cid:104) u A , u A (cid:105) + m C (cid:104) u C , u C (cid:105) = m (cid:48) A (cid:104) u (cid:48) A , u (cid:48) A (cid:105) + m (cid:48) C (cid:104) u (cid:48) C , u (cid:48) C (cid:105) (27)Since p = m u , (27) can be written as1 m A (cid:104) p A , p A (cid:105) + m C (cid:104) p C , p C (cid:105) = m (cid:48) A (cid:104) p (cid:48) A , p (cid:48) A (cid:105) + m (cid:48) C (cid:104) p (cid:48) C , p (cid:48) C (cid:105) (28)We substitute (26) into (28), and we get1 m A (cid:104) p A , p A (cid:105) + m C (cid:104) p C , p C (cid:105) = m (cid:48) A (cid:104) p A − ∆ p , p A − ∆ p (cid:105) + m (cid:48) C (cid:104) p C + ∆ p , p C + ∆ p (cid:105) (29)We can get from (29) a (cid:104) ∆ p , ∆ p (cid:105) + (cid:104) b , ∆ p (cid:105) + c = a = m (cid:48) A + m (cid:48) C > b = m (cid:48) C p C − m (cid:48) A p A , c = ( m (cid:48) A − m A ) (cid:104) p A , p A (cid:105) + ( m (cid:48) C − m C ) (cid:104) p C , p C (cid:105) .The di ff erence between the momentum carried by the matter flow and the momentum I the midpoint of the edge BD reaches a minimum, i.e. min (cid:8) | ∆ p − λ | (cid:9) , in the ∆ p satisfying (30).Thus, we can get the optimization problem (25). Lemma 1.
The binary quadratic function f ( p ) = a (cid:104) p , p (cid:105) + (cid:104) b , p (cid:105) + c (a > ), p ∈ R , if ∃ p , satisfies f ( p ) ≤ ,then f ( p ) = must have a real solution. Theorem 2.
The optimal solution exists in optimization problem (25) , and the expression of the optimal solutionis: ∆ p = rd λ + (1 − rd ) p ∗ , d (cid:44) p ∗ + r n , d = where p ∗ = − a b , r = (cid:115) (cid:104) b , b (cid:105) a − ca, d = (cid:112) (cid:104) λ − p ∗ , λ − p ∗ (cid:105) , n = (1 , T . P roof . Firstly, the existence of solutions is proved.
We assume that f ( p ) = a (cid:104) p , p (cid:105) + (cid:104) b , p (cid:105) + c (32)where a = m (cid:48) A + m (cid:48) C > b = m (cid:48) C p C − m (cid:48) A p A , c = ( m (cid:48) A − m A ) (cid:104) p A , p A (cid:105) + ( m (cid:48) C − m C ) (cid:104) p C , p C (cid:105) .Let p = ∆ mm A p A , have f ( p ) = a (cid:104) p , p (cid:105) + (cid:104) b , p (cid:105) + c = a ∆ m m A (cid:104) p A , p A (cid:105) + ∆ mm A (cid:104) b , p A (cid:105) + c = ( 1 m (cid:48) A + m (cid:48) C ) ∆ m m A (cid:104) p A , p A (cid:105) + ∆ mm A (cid:104) m (cid:48) C p C − m (cid:48) A p A , p A (cid:105) + ( 1 m (cid:48) A − m A ) (cid:104) p A , p A (cid:105) + ( 1 m (cid:48) C − m C ) (cid:104) p C , p C (cid:105) = (cid:104) ( 1 m (cid:48) A + m (cid:48) C ) ∆ m m A − ∆ mm A m (cid:48) A + m (cid:48) A − m A (cid:105) (cid:104) p A , p A (cid:105) + ∆ mm A m (cid:48) C (cid:104) p A , p C (cid:105) + ( 1 m (cid:48) C − m C ) (cid:104) p C , p C (cid:105) (33) i Zhao et al. / Journal of Computational Physics (2020) 9
We substitute (19) into (33), and we get f ( p ) = − m C ∆ mm A ( m C + ∆ m ) (cid:104) p A , p A (cid:105) + ∆ mm A ( m C + ∆ m ) (cid:104) p A , p C (cid:105) − ∆ mm C ( m C + ∆ m ) (cid:104) p C , p C (cid:105) = − ∆ mm A m C ( m C + ∆ m ) (cid:16) m C (cid:104) p A , p A (cid:105) − m A m C (cid:104) p A , p C (cid:105) + m A (cid:104) p C , p C (cid:105) (cid:17) = − ∆ mm A m C ( m C + ∆ m ) (cid:16) m C p A − m A p C (cid:17) Since ∆ m > m A > m C >
0, then f ( p ) ≤
0. According to Lemma 1, we know that there must be a realsolution in formula (34). a (cid:104) p , p (cid:105) + (cid:104) b , p (cid:105) + c = Then the expression of the optimal solution ∆ p is derived. For the convenience of derivation, we set p = ( x , y ) T , b = ( b , b ) T in (34), have a ( x + y ) + b x + b y + c = x + b a ) + ( y + b a ) = b + b a − ca (36)Because there is a real solution, then b + b a − ca ≥
0. We introduce the notation r : = (cid:113) b + b a − ca , p ∗ : = − a b = ( − b a , − b a ) T , d : = | λ − p ∗ | = (cid:112) (cid:104) λ − p ∗ , λ − p ∗ (cid:105) .(1) If r =
0, it is easy to know that ∆ p = p ∗ .(2) If r >
0, then the trajectory of the (36) formula represents a circle with the center of the circle as the p ∗ ,radius as the r . The geometric meaning of the solution of the optimization problem (25) indicates that thedistance from the point on the circle to the λ is the minimum. (a) d = < d < r (c) d = r (d) d > r Fig. 5. Diagram of a circle
1) If d = p ∗ = λ , see Fig. 5(a)), then any point on the circle is the solutionof optimization problem (25). For simplicity, let’s take ∆ p = p ∗ + ( r , T .2) If 0 < d < r (i.e. λ is in the circle, see Fig. 5(b)), then ∆ p = rd λ + (1 − rd ) p ∗ .3) If d = r (i.e. λ is on a circle, see Fig. 5(c)), then ∆ p = λ .4) If d > r (i.e. λ is outside the circle, see Fig. 5(d)), then ∆ p = rd λ + (1 − rd ) p ∗ .In summary, the optimal solution of optimization problem (25) is as follows: ∆ p = rd λ + (1 − rd ) p ∗ , d (cid:44) p ∗ + r n , d = r = (cid:115) (cid:104) b , b (cid:105) a − ca , p ∗ = − a b , d = (cid:112) (cid:104) λ − p ∗ , λ − p ∗ (cid:105) , n = (1 , T . (cid:3) The above eight steps are the complete operation process of the matter compensation flow method. et al. / Journal of Computational Physics (2020)
4. Parallel implementation of matter compensation flow based on OpenMP
OpenMP is a thread level parallel Application Programming Interface (API) based on shared memory. It is com-posed of a set of compilation guidance, run-time routines and environment variables. It has the advantages of simpleprogramming, portability and expansibility, and is widely used in the field of scientific computing. In this paper, aParallel Matter Flow Lagrangian (P-MFL) algorithm is designed for SGH Lagrangian simulation with matter flowbased on OpenMP. The flow chart of the algorithm is shown in Fig. 6.
BeginRead DataInitializationMesh Partition ··· NT − NT ··· NT − NT ··· NT − NT ··· NT − NT ··· NT − NT ··· NT − NT t ≤ T EndNo Yes
Fig. 6. P-MFL algorithm flow chart
Tab. 1 introduces the functions of the six modules in Fig. 6.
Table 1. Function description of the six core modules in the P-MFL
Core modules Function P-DetermineDeltT Parallel computing time step ∆ t P-DynamicEvolve Parallel evolution of a time step P-SetAllDependentVariables Parallel setting of all dependent variables on the element P-CalculateMatterFlowAcc Parallel calculation of matter flow acceleration on three edges of triangular element P-MatterFlowEvolve Parallel evolution of matter compensation flow P-CalculateVertexForce Parallel calculation of node forcesIn the Mesh Partition link in the P-MFL algorithm flow chart, this paper uses the graph partition toolbox METIS[38]to decompose the grid T into NT sub-grids T i ( i = , · · · , NT ). Fig. 7 shows the mesh partition diagram for NT = i Zhao et al. / Journal of Computational Physics (2020) 11 satisfies the minimum principle, so the parallel partition often has good parallel performance. Let’s suppose that thecell index set in T is C and the node index set is N , and the index set in sub grid T i is C i and the node index set is N i , then (1) C = NT (cid:83) i = C i , N = NT (cid:83) i = N i ; (2) C i ∩ C j = φ , N i ∩ N j = φ , ∀ i (cid:44) j , i , j = , , · · · , NT . Fig. 7. Mesh partition ”Fork-Join” is the standard parallel mode of OpenMP, as shown in Fig. 8(a). The code is divided into serial regionand parallel region. The serial region is executed by the main thread. When executing to the parallel region, the slavethread is forked by the system. In the parallel region, the parallel task is completed by the main thread and the slavethread. After the calculation of the parallel region is finished, all threads will join together again. The derived slavethread will exit or block, no longer work, and control the flow return to the main thread and proceed to the next task. (a) ”Fork-Join” mode (b) Diagram of interface point p Fig. 8. ”Fork-Join” mode and Diagram of interface point p In parallel region, thread i ( i = , · · · , NT ) is responsible for computing tasks in subgrid T i , that is, only cellindex set C i and node index set N i are considered, and the calculation task of each thread is about 1 / NT of serialtask. There will be data competition for shared memory cells at the subgrid boundary that may lead to inaccuratecalculation results, such as computing node force at interface p point. Let’s consider two subgrids T i (green area) and T j (yellow area), as shown in Fig. 8(b). Assuming that thread i is calculating the force f Ap of element A on node p ,and thread j is also calculating the force f Fp of element F on node p , and they read the data f p of the shared memoryunit at the same time, then the updated node force will either take the value f (cid:48) p = f p + f Ap in thread i or f (cid:48)(cid:48) p = f p + f Fp in thread j , while the correct result should be f (cid:48)(cid:48) p = f (cid:48) p + f Fp . To ensure the correctness of the calculation results andthe security of the data, techniques such as ”atomic operation” and ”critical region” are used. Because the subgridinterface is short and the cell aggregation is strong in the Mesh Partition step, the probability of the above situation isvery small (the larger the scale, the smaller the possibility), which almost does not a ff ect the parallel e ffi ciency of theprogram. et al. / Journal of Computational Physics (2020)
5. Example and analysis
This section examines the previous matter flow method based on Saltzman Piston Problem[39, 40], Noh Implo-sion Problem[41] and Sedov Explosion Problem[42]. All three examples contain a highly transient shock, and theorientation of the wavefront is inconsistent with that of the grid. The conventional Lagrangian method is easy toappear the physical quantity oscillation problem caused by the sti ff ness of the triangular mesh. Example 5.1.
Consider model question (1) - (4) , domain Ω = [0 , × [0 , . , simulation time t ∈ [0 , . , gas adiabatic γ = . Initial condition: initial density is 1, and the pressure is 0. Boundary condition: the left boundary adoptspiston boundary condition (i.e. the boundary moves to the right at constant unit velocity), and the right and upperand lower boundaries adopt solid wall boundary condition (i.e. normal velocity or displacement is 0). Two grid types of type I and type II, as shown in Fig. 9(a) and Fig. 9(b) respectively, are used to simulate theSaltzman piston problem. (a) Mesh T Sa , I (b) Mesh T Sa , II Fig. 9. The type I grid with T Sa , I of 100 ×
10 and the type II grid with T Sa , II of 100 × Saltzman piston problem with the initial mesh T S a , I and T S a , II is simulated with both the regular SGH Lagrangianmethod (labeled as ”no-matterflow”) and the method implementing the matter flow (labeled as ”matterflow”) in orderto test the e ff ectiveness of the proposed matter flow method.Fig. 10 and Fig. 11 show the grid diagram, density and pressure contour diagram, and Fig. 12 shows the scatterplots of density, pressure and velocity (x), at T S a , I and T S a , II with or without matter flow, at t = .
5. The cell-to-celloscillation is clear in the regular Lagrangian simulation. When the matter flow method is implemented, the pressureoscillation is flattened to a nearly perfect state, which is an expected result for the method, while the oscillation indensity or velocity distributions with and without the matter flow method are with about the same size. Even thoughthe density and velocity distributions do not gain remarkable improvement as the pressure, we can say the result isimproved as a whole.To explore the dependence of the e ff ect of the matter flow method on grid size, we have also carried out simulationsof Saltzman problem using the matter flow method based on refined meshes T S a , Ii and T S a , IIi , i = , T S a , I and T S a , II , as listed in Tab. 2. The simulation results are shown inFig. 13. As the refined mesh corresponds to a smaller viscos, the simulation results are closer to the ideal solutionfor refined mesh, as expected. The oscillation, on the other hand however, does not shrink with the mesh refinement.We interpret this phenomenon as follows: the e ff ect of the matter flow method depends on the smoothness of thephysical quantities distributions along cells, the smoother the better; when the mesh is refined, the viscous alsobecomes smaller, this leads to sharper distributions of physical quantities along space for shock wave problem, andif transferring from space to cells, it leads to physical quantity distribution smoothness approximately independent ofthe mesh size, and so the oscillation is also approximately independent of the mesh size. So, this phenomenon can beviewed as a nature of the matter flow method combined with the SGH Lagrangian scheme. This analysis also point a i Zhao et al. / Journal of Computational Physics (2020) 13 way of how to obtain result with the oscillation better alleviated than that in Fig. 13: to construct more refined meshand, at the same time, increase the viscous coe ffi cient. (a) Mesh, no-matterflow, T Sa , I (b) Mesh, matterflow, T Sa , I (c) Density, no-matterflow, T Sa , I (d) Density, matterflow, T Sa , I (e) Pressure, no-matterflow, T Sa , I (f) Pressure, matterflow, T Sa , I Fig. 10. T Sa , I mesh diagram, density and pressure contour diagram at t = . . (a) Mesh, no-matterflow, T Sa , II (b) Mesh, matterflow, T Sa , II (c) Density, no-matterflow, T Sa , II (d) Density, matterflow, T Sa , II (e) Pressure, no-matterflow, T Sa , II (f) Pressure, matterflow, T Sa , II Fig. 11. T Sa , II mesh diagram, density and pressure contour diagram at t = . . et al. / Journal of Computational Physics (2020)(a) Density, T Sa , I (b) Pressure, T Sa , I (c) Velocity (x), T Sa , I (d) Density, T Sa , II (e) Pressure, T Sa , II (f) Velocity (x), T Sa , II Fig. 12. The scatter plot of density, pressure and velocity (x) of T Sa , I and T Sa , II with or without matter flow at t = . . (a) Density, T Sa , Ii ( i = , ,
3) (b) Pressure, T Sa , Ii ( i = , ,
3) (c) Velocity (x), T Sa , Ii ( i = , , T Sa , IIi ( i = , ,
3) (e) Pressure, T Sa , IIi ( i = , ,
3) (f) Velocity (x), T Sa , IIi ( i = , , Fig. 13. The density, pressure and velocity x scatter diagrams of T Sa , Ii and T Sa , IIi , i = , , at the time of t = . . i Zhao et al. / Journal of Computational Physics (2020) 15
Table 2. Type I and type II grid size information
Type Notation Mesh resolution T S a , I T S a , I T S a , I T S a , II T S a , II T S a , II Example 5.2.
Consider the model problem (1) - (4) , where domain Ω = [0 , . × [0 , . , simulation time t ∈ [0 , . (General reference t = . , but the grid distortion in our calculation is too severe, we can not calculate this time, sowe take t = . ), gas adiabatic γ = . Initial condition: initial density is 1, specific internal energy is 0, and velocityis 1 (the direction points to the origin, i.e. lower left corner (0, 0) position). Boundary condition: the left and lowerboundaries of the domain adopts the solid wall boundary condition, and the right and upper boundaries adopt freesurfaces conditions. The domain Ω is divided by a uniform grid of 40 ×
40, 80 ×
80 and 160 ×
160 (see Fig. 14), whichare recorded as T No , T No and T No , respectively. Remark 3.
Radius represents the radius, i.e.
Radius = (cid:112) x + y , and the radial velocity represents the velocity alongthe radius direction, which is (cid:113) u x + u y and points to the origin. (a) T No (b) T No (c) T No Fig. 14. Three kinds of initial meshes
Firstly, the influence of matter flow on the experimental results is explored for the initial mesh T No . Fig. 15shows the grid and density, pressure contour diagram at t = .
4. Fig. 16 shows the density, pressure, radial velocityscatter diagram at t = .
4. The big di ff erence between the numerical solution and the ideal solution in Fig. 16 is due tothe large artificial viscosity corresponding to the coarse mesh. As can be seen from the diagram in Fig. 15 and Fig. 16,the introduction of the matter flow method can greatly alleviate the physical quantity oscillation in SGH Lagrangesimulation. As a side e ff ect, matter flow method also reduces the distortion of the mesh in the simulation of Nohproblems. Then, the influences of di ff erent grid sizes on the results are explored for three groups of initial meshes T No , T No and T No . Fig. 17 shows the corresponding scatter diagram of density, pressure and velocity at t = . et al. / Journal of Computational Physics (2020)(a) Local mesh, no-matterflow (b) Density, no-matterflow (c) Pressure, no-matterflow(d) Local mesh, matterflow (e) Density, matterflow (f) Pressure, matterflow
Fig. 15. Mesh, density, pressure contour diagram of T No at t = . . (a) Density, matterflow (b) Pressure, matterflow (c) Velocity(r), matterflow Fig. 16. The density, pressure and radial velocity scatter diagram of T No at t = . . (a) Density, matterflow (b) Pressure, matterflow (c) Velocity(r), matterflow Fig. 17. The density, pressure and velocity scatter diagram of di ff erent grid sizes at t = . . i Zhao et al. / Journal of Computational Physics (2020) 17
Example 5.3.
Consider the model problem (1) - (4) , where domain Ω = [0 , × [0 , and simulation time t ∈ [0 , ,gas adiabatic γ = . . We use × , × , × , × , × , × uniform mesh for thedomain Ω . Tab. 3 shows the corresponding notation, number of grid cells and number of nodes, and Fig. 18 showsthe first three sets of grids. Initial condition: initial density is 1, velocity is , square domain in the lower left corner [0 , h i ] × [0 , h i ] (consists of a pair of triangles). The internal energy density is E / h i ( i = , , · · · , for the pair ofcells at the down-left corner, where E = . , h i = . × − ( i + , i = , , · · · , , h = − , and is 0 for the restcells. Boundary condition: the left and lower boundary of the domain adopts the solid wall boundary conditions, Theright and upper bounds adopt free surfaces conditions. Table 3. Six kinds of initial grid information of Sedov
Notation Mesh resolution Number of elements Number of nodes T S e ×
40 3200 1681 T S e ×
80 12800 6561 T S e ×
160 51200 25921 T S e ×
320 204800 103041 T S e ×
640 819200 410881 T S e × (a) T Se (b) T Se (c) T Se Fig. 18. Sedov the first three sets of initial meshes of the problem.
Firstly, the e ff ect of material flow on numerical simulation results is discussed for initial mesh T S e . Fig. 19shows the contour diagram of the grid, density and pressure when the T S e grid does not use matter flow and usesmatter flow at t =
1. Fig. 20 shows the scatter diagram of density, pressure and radial velocity at t =
1. It can be seenfrom the figure that the introduction of the matter flow method can greatly alleviate the physical oscillation in the SGHLagrangian simulation. As a side e ff ect, the matter flow method also reduces the mesh distortion in the simulation ofthe Sedov problem. Secondly, the e ff ects of di ff erent grid sizes on the numerical solutions are investigated for the initial meshes T S e , T S e , T S e . Fig. 21 shows the scatter diagram of the density, pressure, radial velocity of the three initial meshsizes at t =
1, respectively. The conclusions of these simulations are similar to that of Saltzman and Noh problems.
Finally, we take the Sedov problem as an example to explore the parallel performance of the SGH La-grangian program after adding the matter flow (the parallel test results of the other two examples are basicallyconsistent with the examples). (1) To validate the correctness of the parallel algorithm, considering the initial grid T S e , the influence of di ff erentthread number NT = , , , ,
16 on the numerical solution is studied.Fig. 22 shows the curves of total mass, total energy and total momentum of di ff erent threads with time. Fig. 23 showsthe scatter plot of density, pressure and radial velocity of di ff erent threads at t = et al. / Journal of Computational Physics (2020)(a) Mesh, no-matterflow (b) Density, no-matterflow (c) Pressure, no-matterflow(d) Mesh, matterflow (e) Density, matterflow (f) Pressure, matterflow
Fig. 19. The grid, density and pressure contour diagram of T Se at t = . (a) Density (b) Pressure (c) Velocity (r) Fig. 20. The scatter diagram of density, pressure and radial velocity of T Se at t = . (a) Density (b) Pressure (c) Velocity (r) Fig. 21. The density, pressure and radial velocity scatter diagram of T Se , T Se and T Se at t = . i Zhao et al. / Journal of Computational Physics (2020) 19(a) Total mass (b) Total energy (c) Total momentum(x) (d) Total momentum(y)
Fig. 22. Total mass, total energy, total momentum variation curves of di ff erent threads. It can be seen from the diagram that the totalmass and total energy are conserved, and the experimental results of multithreading are completely consistent, and the total momentumincreases with the increase of simulation time, and the total momentum change curve of multithreading completely overlaps. (a) Density (b) Pressure (c) Radial velocity Fig. 23. The scatter plot of density, pressure and radial velocity of di ff erent threads at t = . It can be found that the density, pressure andradial velocity calculated by di ff erent threads are exactly the same. Figs. 22-23 shows that the total mass, total energy and total momentum change curve of multithreading completelycoincide with the change curve of single thread. Density, pressure and radial velocity scatter plot of multithreading at t = scalability of parallel programs, three groups of large-scale grids T S e , T S e , T S e , fixed iterationtimes of 1000 times, the number of threads are NT = , , , , ff erent sizes in di ff erent thread numbers, re-spectively. Table 4. P-MFL CPU-Time (s) NT T S e T S e T S e +
02 1.42E +
03 3.35E +
032 1.91E +
02 7.53E +
02 1.75E +
034 1.01E +
02 3.92E +
02 9.11E +
028 5.74E +
01 2.21E +
02 5.19E + +
01 1.48E +
02 3.52E + Table 5. P-MFL Speedup NT T S e T S e T S e \ \ \ Remark 4.
The parallel speedup S n = T T n , where the T is the time of execution of one processor, T n the time ofexecution of n processor.Tab. 4 shows that when the grid amount is fixed and the number of threads increases, P-MFL CPU wall timegradually decreased. When the number of threads is fixed and the grid size increases, P-MFL CPU wall increasesalmost linearly with grid size. Tab. 5 shows that, when the grid is fixed, P-MFL speedup increases with the increaseof the number of threads. When the number of threads is fixed, the speedup increases gradually with the grid amountincrease. Especially, If the grid is T S e (the number of cells is 2 million) and the number of threads is 16, the speedupreached 9.51. et al. / Journal of Computational Physics (2020)
To sum up, the P-MFL parallel algorithm based on OpenMP is correct and has good parallel scalability.
6. Summary and discussion
In this paper, aiming at the checkerboard oscillation problem of triangular mesh SGH Lagrangian hydrodynamicsimulation, a matter flow method is designed to alleviate the physical quantity oscillation, and parallelization iscarried out. The matter flow method is similar to that of Scovazzi[35] and Molgan[36] —- by introducing somephysical quantity transport terms between elements. However, compared with these two methods, we think that themethod in this paper takes into account the physical quantities that need to be transported more comprehensively.Three kinds of e ff ects are considered in the matter flow method. Firstly, the mass, energy and momentum transportcaused by matter transport. Secondly, energy transport caused by work due to the change of element density. Finally,the e ff ect of matter flow on strain rate in the element. In contrast, Scovazzi’s method only considers energy transport,while Molgan’s method only considers the first kind of e ff ect. The e ff ectiveness of the proposed method is verified bynumerical experiments.Although the matter flow method in this paper has achieved some good results, there are still many limitations.Firstly, it is only suitable for scalar viscosity, and how to extend it to tensor viscosity needs further study. Secondly,when one hopes to simulate multi-material fluids, the problem of how to implement matter flow between cells withdi ff erent materials needs to be solved. Finally, more work is needed in the parallel implementation algorithm of matterflow algorithm, such as designing parallel method based on MPI. Acknowledgments
This work was supported by the National Natural Science Foundation of China (NSFC Project number 11971414).The authors would like to thank Long Xie and Shuchao Duan for useful discussions.
References [1] J. W. Banks, D. W. Schwendeman, A. K. Kapila, et al, A high-resolution godunov method for compressible multi-material flow on overlappinggrids, Journal of Computational Physics 223 (2007) 262–297.[2] D. De Niem, E. Kuhrt, U. Motschmann, A volume-of-fluid method for simulation of compressible axisymmetric multi-material flow, Com-puter Physics Communications 176 (2007) 170–190.[3] J. G. Zheng, T. S. Lee, S. H. Winoto, Numerical simulation of richtmyer-meshkov instability driven by imploding shocks, Mathematics andComputers in Simulation 79 (2008) 749–762.[4] S. K. Sambasivan, H. S. Udaykumar, Sharp interface simulations with local mesh refinement for multi-material dynamics in strongly shockedflows, Computers & Fluids 39 (2010) 1456–1479.[5] H. W. Zheng, C. Shu, Y. T. Chew, N. Qin, A solution adaptive simulation of compressible multi-fluid flows with general equation of state,International Journal for Numerical Methods in Fluids 67 (2011) 616–637.[6] P. Movahed, E. Johnsen, A solution-adaptive method for e ffi cient compressible multifluid simulations, with application to the richtmyer-meshkov instability, Journal of Computational Physics 239 (2013) 166–186.[7] H. Chen, W. B. Zhu, X. B. Zhang, et al, Application of real ghost fluid method to simulation of compressible multi-fluid flows, Explosion &Shock Waves 33 (2013) 29–37.[8] W. Zhang, W. Ye, J. Wu, et al, Hydrodynamic instabilities of laser indirect-drive inertial-confinement-fusion implosion, Sci. Sin.-Phys.Mech. Astron. 44 (2014) 1–23.[9] A. Main, C. Farhat, A second-order time-accurate implicit finite volume method with exact two-phase riemann problems for compressiblemulti-phase fluid and fluid-structure problems, Journal of Computational Physics 258 (2014) 613–633.[10] A. Kapahi, C. Hsiao, G. L. Chahine, A multi-material flow solver for high speed compressible flows, Computers & Fluids 115 (2015) 25–45.[11] S. Diot, M. M. Francois, E. D. Dendy, A higher-order unsplit 2d direct eulerian finite volume method for two-material compressible flowsbased on the mood paradigms, International Journal for Numerical Methods in Fluids 76 (2014) 1064–1087.[12] C. D. Sijoy, S. Chaturvedi, An eulerian multi-material scheme for elastic-plastic impact and penetration problems involving large materialdeformations, European Journal of Mechanics B-fluids 53 (2015) 85–100.[13] D. Pavlidis, J. L. M. A. Gomes, Z. Xie, et al, Compressive advection and multi-component methods for interface-capturing, InternationalJournal for Numerical Methods in Fluids 80 (2016) 256–282.[14] Z. He, B. Tian, Y. Zhang, et al, Characteristic-based and interface-sharpening algorithm for high-order simulations of immiscible compress-ible multi-material flows, Journal of Computational Physics 333 (2017) 247–268.[15] L. Wang, G. M. D. Currao, F. Han, et al, An immersed boundary method for fluid-structure interaction with compressible multiphase flows,Journal of Computational Physics 346 (2017) 131–151.[16] M. L. Wong, S. K. Lele, High-order localized dissipation weighted compact nonlinear scheme for shock- and interface-capturing in com-pressible flows, Journal of Computational Physics 339 (2017) 179–209.[17] C. Liu, C. Hu, Adaptive thinc-gfm for compressible multi-medium flows, Journal of Computational Physics 342 (2017) 43–65.i Zhao et al. / Journal of Computational Physics (2020) 21[18] A. Barlow, P. Maire, W. J. Rider, et al, Arbitrary lagrangian-eulerian methods for modeling high-speed compressible multimaterial flows,Journal of Computational Physics 322 (2016) 603–665.[19] H. R. Anbarlooei, K. Mazaheri, Moment of fluid interface reconstruction method in multi-material arbitrary lagrangian eulerian (mmale)algorithms, Computer Methods in Applied Mechanics and Engineering 198 (2009) 3782–3794.[20] M. Kucharik, R. V. Garimella, S. P. Schofield, et al, A comparative study of interface reconstruction methods for multi-material ale simula-tions, Journal of Computational Physics 229 (2010) 2432–2452.[21] B. Tian, W. Shen, S. Jiang, et al, A global arbitrary lagrangian-eulerian method for stratified richtmyer-meshkov instability, Computers &Fluids 46 (2011) 113–121.[22] S. Galera, J. Breil, P. Maire, A 2d unstructured multi-material cell-centered arbitrary lagrangian-eulerian (ccale) scheme using mof interfacereconstruction, Computers & Fluids 46 (2011) 237–244.[23] Z. Jia, J. Liu, S. Zhang, An e ff ective integration of methods for second-order three-dimensional multi-material ale method on unstructuredhexahedral meshes using mof interface reconstruction, Journal of Computational Physics 236 (2013) 513–562.[24] Q. H. Zeng, Mmale numerical simulation for multi-material large deformation fluid flows, Journal of Physics Conference 510 (2014) 012047.[25] J. Cheng, C. Shu, Positivity-preserving lagrangian scheme for multi-material compressible flow, Journal of Computational Physics 257(2014) 143–168.[26] M. Wicke, D. Ritchie, B. M. Klingner, et al, Dynamic local remeshing for elastoplastic simulation, Acm Transactions on Graphics 29 (2010)49.[27] R. L. Wang, Z. Lin, L. Wei, Reconnection-based lagrangian-local remeshing method for large deformations, Journal of ComputationalPhysics 28 (2011) 501–506.[28] R. L. Wang, Z. Lin, W. Z. Wen, et al, Development and application of adaptive multi-media lagrangian fluid dynamics software lad2d,Computer Aided Engineering 23 (2014) 1–7.[29] J. Liu, A second-order changing-connectivity ale scheme and its application to fsi with large convection of fluids and near contact ofstructures, Journal of Computational Physics 304 (2016) 380–423.[30] J. Cheng, C. Shu, Second order symmetry-preserving conservative lagrangian scheme for compressible euler equations in two-dimensionalcylindrical coordinates, Journal of Computational Physics 272 (2014) 245–265.[31] G. Georges, J. Breil, P. Maire, A 3d gcl compatible cell-centered lagrangian scheme for solving gas dynamics equations, Journal ofComputational Physics 305 (2016) 921–941.[32] H. B. Zhao, Two dimensional Lagrangian simulation of large deformation motion of multi-material compressible fluid, Master’s thesis, ChinaAcademy Of Engineering Physics, 2018.[33] H. B. Zhao, B. Xiao, J. S. Bai, et al, Simulation of two-dimensional multi-material compressible flows using lagrangian methods, Journal ofHigh Pressure Physics 144(04) (2018) 47–59.[34] J. Waltz, N. R. Morgan, T. R. Canfield, et al, A nodal godunov method for lagrangian shock hydrodynamics on unstructured tetrahedral grids,International Journal for Numerical Methods in Fluids 76 (2014) 129–146.[35] G. Scovazzi, Lagrangian shock hydrodynamics on tetrahedral meshes: A stable and accurate variational multiscale approach, Journal ofComputational Physics 231 (2012) 8029–8069.[36] N. R. Morgan, J. Waltz, D. E. Burton, et al, A godunov-like point-centered essentially lagrangian hydrodynamic approach, Journal ofComputational Physics 281 (2015) 614–652.[37] B. Chapman, G. Jost, R. Pas, Using openmp:portable shared memory parallel programming, Journal of Computer Science & Technology 10,no. 3 (2010).[38] G. Karypis, V. Kumar, A software package for partitioning unstructured graphs, partitioning meshes, and computing fill-reducing orderingsof sparse matrices, Landolt B¨ornstein Group III Condensed Matter (1998) 372–374.[39] J. K. Dukowicz, B. J. A. Meltz, Vorticity errors in multidimensional lagrangian codes, Journal of Computational Physics 99 (1992) 115–134.[40] L. G. Margolin, A centered artificial viscosity for cells with large aspect ratio, Nasa Sti //