Elimination of hourglass distortions by means of the Lagrangian–Gauss-point-mass method for compressible fluid flows

Accepted Manuscript Elimination of Hourglass Distortions by Means of the Lagrangian-Gauss-PointMass Method for Compressible Fluid Flows Junxia Cheng, Baolin Tian PII: DOI: Reference:

S0045-7930(15)00128-0 http://dx.doi.org/10.1016/j.compfluid.2015.04.013 CAF 2862

To appear in:

Computers & Fluids

Received Date: Revised Date: Accepted Date:

27 February 2014 10 June 2014 12 April 2015

Please cite this article as: Cheng, J., Tian, B., Elimination of Hourglass Distortions by Means of the LagrangianGauss-Point-Mass Method for Compressible Fluid Flows, Computers & Fluids (2015), doi: http://dx.doi.org/ 10.1016/j.compfluid.2015.04.013

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Elimination of Hourglass Distortions by Means of the Lagrangian-Gauss-Point-Mass Method for Compressible Fluid Flows Junxia Cheng and Baolin Tian Institute of Applied Physics and Computational Mathematics, Beijing, China

Abstract Staggered Lagrangian schemes on quadrilaterals or hexahedrons with constant stress have been troubled by hourglass distortions for many years. Hourglass viscosities have been often used to suppress hourglass distortions resulting in the nonconservative total energy, even so computations often halt abnormally because of spurious large grid distortions. It is difficult to discern the spurious distortions from the physical in multidimensional problems, so we want to get the scheme eliminating hourglass distortions. This paper gives a new Lagrangian scheme for compressible flows to eliminates the grid hourglass motions by means of the Lagrangian Gauss-point Mass(LGM) method. In the LGM method, pressures are varied in a zone similar to the subzonal pressure method, but there are several important differences between them. Numerical examples show that the LGM method can eliminate the hourglass motions for many problems, and it has the good property of preserving the conservation of total energy nearly in the two-dimensional plane geometry. Because of the area-weighted scheme, it cannot preserve the total energy but improves it to a great extent. The procedure of the LGM method is simple to be implemented into the existing 2D finite element codes on quadrilaterals, and it can be extended to the 3D finite element codes on hexahedrons without theoretical difficulties. Keywords: Staggered Lagrangian scheme; Hourglass distortions; Finite element method; Quadrilaterals

1. Introduction Lagrangian hydrodynamic algorithms have been widely used for perhaps the longest time of numerical methods employed for the solution of the complex problem of multi-material fluid flows[1], because they can track the material interfaces clearly and avoid the numerical error of the advection term. Lagrangian hydrodynamics algorithms can be classified into the staggered Lagrangian and cell-centered Lagrangian algorithms. Cell-centered Lagrangian algorithms define physical variables and velocities at zone centers, which need to solve the Riemann problem on the cell boundary usually. Staggered Lagrangian algorithms define physical variables of density, internal energy and pressure at zone centers, but velocities are defined at grid points, which need the artificial viscosity usually to resolve the shock. Staggered Lagrangian schemes have the overwhelming advantages of being simple and inexpensive[2], so they have been widely used in applications for many years. However, constant stress staggered Lagrangian algorithms on quadrilaterals have been troubled by hourPreprint submitted to Computers & Fluids

glass motions since their inception. The hourglass mode derives its name from the fact that the deformed element literally resembles an hourglass, which is also called keystone mode or spurious zero-energy mode. Hourglass mode will grow without bound and turn the element inside out soon without hourglass control, shown in Fig. 1. The hourglass modes were first noted in the work by Maenchen and Sack[1, 3], and they used some rotational hourglass viscosity in Tensor code to control the hourglass distortion. Wilkins[6] developed the tensor triangle viscosity to prevent the hourglass pattern grid distortions. Thei triangle viscosity is formulated for each of the triangles in the four zones surrounding the vertex. Hallquist[7] presented several hourglass viscosities including tensor triangle viscosity, Maenchen-Sack rotational viscosity, Flanagan-Belytschko viscosity[8] and Key’s HONDO viscosity in DYNA2D. The choice of hourglass viscosity depends on the problem, moreover incorrect computational results can be given possibly if hourglass viscosity or its coefficient is unsuitable. The hourglass viscosity, which is not physical April 23, 2015

0.6

0.4

0.4

R(cm)

R(cm)

0.6

0.2

0.2

0 0

0.2

0.4

0

0.6

0

Z(cm)

0.2

0.4

0.6

Z(cm)

Figure 1: Comparison the mesh of Sedov problem without hourglass control with that with hourglass control at time t=0.2. left: result without hourglass control, right: result with hourglass control

but numerical[9], is added into the momentum equation but not into the energy equation, so the staggered Lagrangian scheme with the hourglass viscosity can not preserve the conservation of total energy. Except for hourglass viscosities, other methods were developed to suppress hourglass modes of grids in the staggered Lagrangian scheme for compressible fluid flows. An edge-centered tensor artificial viscosity[4] was given to resolve the shock which is used as the replacement of the Richtmyer-von Neumann viscosity. This tensor artificial viscosity works on the edge which is under compression, so it can provide some resistance to the hourglass distortion. Caramana and Shashkov[1] eliminated hourglass-type motions and spurious grid distortions by means of Lagrangian subzonal masses and associated pressures. In the case of two dimensional quadrilaterals, the Lagrangian cell is divided into four subzones with Lagrangian masses. Four subzonal densities are different from the density of the Lagrangian zone center, hence subzonal pressures are different also. Subzonal pressures were used multiplied by some coefficient in applications, which depends on the problem also. The tensor artificial viscosity using a mimetic finite difference algorithm was developed to eliminate hourglass distortions and spurious grid distortions[10] after the subzonal pressure method was developed. Numerical examples showed that the tensor viscosity can reduce the dependence of the solution on the relation of the grid to the flow structure[10]. Barlow[11] extended the compatible Lagrangian algorithm into the compatible finite element Lagrangian scheme. A cell-centered multi-dimensional approximate Riemann solver was used to build a form of artificial viscosity in the formulation of the compatible staggered Lagrangian hydrodynamics scheme[12].

Works of references[1, 10, 11, 12] were carried out on the compatible staggered Lagrangian hydrodynamics formulation introducing the concept of the corner force and the corner mass, which is different from the traditional staggered Lagrangian scheme[3, 13, 16]. Here, we will carry out our work on the traditional staggered Lagrangian scheme which has been used for many years. In order to eliminate the hourglass distortion of the traditional hydrocode, a Lagrangian scheme is given here based on the finite element method for compressible fluid flows. Four Gauss-points of the Lagrangian quadrilateral cell carry Lagrangian masses in the twodimensional geometry, so each Gauss-point has its own pressure. This method is similar to the subzonal pressure method to some extent, but it needs not the coefficient depending on the problem used in the suzonal pressure method. The difference between the LGM method and the subzonal pressure method will be explained in detail in section 2.3. 2. The isoparametric finite element scheme on quadrilaterals and the LGM method 2.1. Preliminaries of the bilinear isoparametric finite element method Firstly, we would like to introduce some definitions for clarity. The spatial(Euler) coordinates are denoted by xi or x, and material coordinates are denoted by Xi or X which are independent on the time t. The standard technique in finite element is to express variables in terms of element coordinates, which are sometimes called as parent coordinates or natural coordinate[9]. Here, ξi denote the element coordinates, which are functions of Xi . xi is equivalent to (r, z) and ξi is equivalent to (ξ, η) in the two-dimensional geometry in this paper. 2

Three types of notations such as xi , x, (r, z) will be used in different conditions for convenience and clarity in the following. We just focus on the 4-node quadrilateral element in this paper. The domain Ω is discretized into a series of quadrilateral elements, and Ωe denotes the element domain. The transformation between spatial coordinates and element coordinates is given by on the 4-node quadrilateral element xi =

4 ∑

NI xiI ,

J is the determinant of Jacobian matrix J, which will be used extensively in the following. More details about this section can be found in the reference book[9]. 2.2. Finite element scheme of fluid dynamics equations The governing equations of Lagrangian hydrodynamics are written as  dρ   + ρ∇ · (v) = 0,    dt dv (7) ρ   dt = −∇p,    ρ de = −p∇ · v, dt

(1)

where v is the velocity vector, ρ is density, e is specific internal energy, p is the pressure, dtd is the material derivative. The integral forms of the mass equation and the momentum equation are written into ∫  d  ρdV = 0,   dt   ∫ Ω ∫ (8)  d     dt ρedV = − p∇ · vdV.

I=1

where 4 nodes of the element are ordered anticlockwise, NI is the shape function, xiI are the nodal spatial coordinates. The displacement field is give as ui =

4 ∑

NI uiI ,

(2)

I=1

Ω

where uiI are the nodal displacements. Because same functions are used for the mapping (1) and the displacement interpolation (2), this element is called an isoparametric element. The bilinear isoparametric shape function NI is defined for the 4-node quadrilaterals as NI =

1 (1 + ξI ξ)(1 + ηI η), 4

Divergence of velocity vector of two-dimensional plane and axisymmetrical problem can be written into the unified form ∇·v=

∂z

∂z

∂η

j

tn are denoted with superscripts n. By means of Petrov Galerkin finite element approach, weak form of the momentum equation is discretized into ∫ ∫ dv δui ρ drdz = − δui ∇pdrdz, (11) dt

(5)

Ω

= 14 ξI (1 + ηηI ), = 14 ηI (1 + ξξI ).

Ω

where δui is the virtual displacement, δui =

where derivatives of NI with respect to ξ and η can be obtained from Eq.(3), ∂NI ∂ξ ∂NI ∂η

(9)

where E j is the ∫ internal energy of element j, ρ j is its density, V j = Ωe rα drdz. In this paper, variables at time

∂η

where J is the Jacobian matrix, which is given by  4  4  ∑ ∂NI ∑  ∂N I rI zI  ∂r ∂z    ∂ξ   ∂ξ ∂ξ   I=1 ∂ξ I=1  , J =  ∂r ∂z  =  4 4   ∑ ∂NI ∑ ∂N ∂η ∂η I   rI zI ∂η ∂η I=1 I=1

∂vr ∂vz vr + +α , ∂r ∂z r

where α = 0 in the two-dimensional plane geometry, and α = 1 in the axisymmetrical geometry where z is the axis of symmetry. Using the time central difference method to the energy conservation equation, discretization of Eq.(8) on element j are written as  n+1  ρnj V nj = ρn+1 = ρ0j V 0j ,  j Vj (10)  n+1/2  E n+1 − E n = −p (V n+1 − V nj ), j j j j

(3)

where ξI = (−1, 1, 1, −1),ηI = (−1, −1, 1, 1). The transformations are one-to-one between the two of spatial coordinates xi , element coordinates ξi and material coordinates Xi . Partial derivatives of the shape functions with respect to the spatial coordinates r, z are obtained from Eq.(4) and Eq.(5),  ∂NI  ( ∂NI )  ∂ξ ∂η   ∂NI   ∂ξ   ∂ξ      −1   ∂r ∂r ∂r =  ∂ξ   ∂NI  = J  ∂NI  , (4) ∂NI ∂η  ∂z

Ω

N ∑

NI δuiI .

I=1

Integrated by parts, Eq.(11) is written into ∫ ∫ dv δui ρ drdz = −(pδui )|∂Ω + p∇(δui )drdz, (12) dt

(6)

Ω

3

Ω

where the first term of right hand is given by the boundary condition. The momentum equation is discretized by the areaweight scheme(Petrov Galerkin) in order to preserve the symmetry property. Eq.(12) is discretized on the inner elemental domain into ∫ ∫ dvk ρNI Nk drdz = p∇NI drdz. (13) dt Ωe

used the bulk viscosity defined in Eq. (17)[7, 17, 18], which has been widely used for many years. { a2 ρ(l∇ · v)2 − a1 ρcl∇ · v, ∇ · v < 0, Q= (17) 0, ∇ · v ≥ 0, where ρ is the density, c is the sound speed, l is the characteristic length, v is the velocity vector, a1 , a2 are the given coefficients which are chosen here a1 = 0.06, a2 = 1.5 as usual[16]. In this paper, this scheme, composed of Eq.(10), Eq.(15), Eq.(16), Eq.(17) and the hourglass viscosity, is called the old scheme in contrast with the new scheme developed in the following context. The finite element scheme of one-point Gauss quadrature is equivalent to the traditional difference schemes of Tensor and HEMP. When the velocity v2I or v1I of nodes on the quadrilateral element is (1, −1, 1, −1)in the two-dimensional plane and axisymmetrical geometry, the divergence of velocity is zero on the Gauss point(ξ = 0, η = 0). In other words, this type of distortions, called hourglass modes here, will not produce any restoring forces and will persist in the computation. Without hourglass control, hourglass distortions will spread over all the elements and turn the element inside out soon. In the past years, hourglass viscosities have been widely used to suppress the hourglass-type distortions in the complex compressible large-distortion problems. Examples of the old scheme in this paper are the results computed with hourglass viscosity. Hourglass viscosities have some coefficients depending on the physical problem, so it is not known whether the grid distortion is physical or not for the multi-dimensional computation. Besides, the conservation of total energy cannot be preserved because one part of energy will be consumed by the hourglass viscosity. Consequently, the Lagrangian schemes without hourglass motions have been pursued for compressible fluid flows for many years.

Ωe

Eq.(13) has the consistent mass matrix. When the mass is lumped, the mass matrix is diagonal and the elemental area mass is distributed evenly to its nodes. When the lumped mass matrix is used, the momentum equation on one element is give by 1 e dvI m = 4 dt

∫

1 −1

∫

1

−1

p∇NI Jdξdη,

(14)

∫ where me is the area mass of the element, me = Ωe drdz, J is the determinant of Jacobian matrix defined in Eq.(4). Eq.(14) is much simpler than Eq.(13) to be solved especially in explicit time discretization scheme, but a lumped mass matrix is more than just a computational simplification, and it also gives better answers for impulsive loads[16]. We denote the time increments ∆tn+1/2 = tn+1 − n t , ∆tn = 12 (∆tn−1/2 + ∆tn+1/2 ). The momentum equation of Eq.(14) are discretized by the time central difference scheme into ) 1 e ( n+1/2 m vI − vn−1/2 = ∆tn I 4

∫

1

−1

∫

1 −1

(p∇NI J)n dξdη.

The central difference method is also used to advance the nodes’ positions in time by xiIn+1 = xiIn + vn+1/2 ∆tn+1/2 . iI

(15)

2.3. The LGM method eliminating hourglass motions of the finite element scheme When an element deforms in an hourglass mode, its density and internal energy will not change shown in Eq.(10), so any forces will not be produced to resist this type of distortion because the pressure depends on the density and the internal energy. When the pressure is a constant in the element, the scheme can not represent the hourglass modes. So, the key cause of hourglass distortions is that the pressure is constant on the element. The subzonal pressure method defined four subzonal pressures on four subzones of the Lagrangian cell, and it did suppress the hourglass-type distortions and spurious vorticities on high-aspect-ratio grids[1]. Inspired by

In the formulation of traditional finite element scheme, the pressure is constant in one element, so onepoint Gauss quadrature is enough to compute the integral momentum equation of Eq.(14) accurately, ) 1 e ( n+1/2 m vI − vn−1/2 = I 4 4∆tn (p∇NI J)n (ξ = 0, η = 0).

(16)

In order to resolve the shock, the artificial viscosity has to be included working with the pressure[17, 18]. So the pressure in the momentum equation and internal energy should be considered as p + Q in the above. We 4

this idea and based on our analysis of the old scheme, we develop a finite element scheme with varied pressures on the element to suppress hourglass distortions for compressible fluid flows. Integral forms of mass equation and momentum equation are carried out by Gauss quadrature. Mass equation of Eq.(10) is written into ∫ ∫1 ∫1 d ρrα drdz = dtd −1 −1 ρrα Jdξdη = dt Ωe (18) nq nq ∑ d ∑ ρq Vq = dtd mq = 0, dt q=1

q=1

Figure 2: Quadrilateral element with Gauss points denoted by △ and subzones. point 5, 6, 7, 8 are the midpoints of the edges, point 9 is the element center[1].

where nq is the number of Gauss points, ωq is the weighting value ( ) at Gauss point q, Vq = α ωq r J ξ = ξq , η = ηq is called the volume of Gauss point q. Here, we assume that the mass of Gauss point mq is Lagrangian, and nq is equal to 4. When the nodes’ coordinates are known, densities of Gauss points can be derived according to the assumption that masses of Gauss points are Lagrangian. Here it should be stressed that masses of Gauss point are different from subzonal masses on the quadrilateral shown in Fig.2. The area of Gauss point is AGI = J(ξ = ξI , η = ηI ). We can get the difference between the area of Gauss point and the corresponding subzonal area: √ 1 (+3 − 2 3)[+(r1 + r3 )z24 − 48 r2 (z14 + z34 ) + r4 (z12 + z32 )] √ 1 = (−3 + 2 3)[−(r2 + r4 )z13 + 48 r1 (z23 + z43 ) − r3 (z41 − z12 )], √ 1 = (−3 + 2 3)[−(r1 + r3 )z24 + 48 r2 (z14 + z34 ) − r4 (z12 + z32 )], √ 1 = (+3 − 2 3)[+(r2 + r4 )z13 − 48 r1 (z23 + z43 ) + r3 (z41 − z12 )],

The semi-discrete equations of Eq.(14) and Eq.(21) are computed by the use of four-point Gauss quadrature,  4 ∑  1 e dviI   m dt = pq Liq ,   4   q=1 (22)  ∫ 4  ∑  d   pq Vq ,   dt ρedV = − Ωe

) I where Liq = ωq ∂N J ξ = ξq , η = ηq , Vq = ωq ∇ · ∂x i ( ) vrα J ξ = ξq , η = ηq . The four Gauss points locate as follows: Point 1: (− √13 , − √13 ),

AG1 − A1598 =

AG2 − A2695

AG3 − A3796

AG4 − A4897

Point 2: (+ √13 , − √13 ), Point 3: (+ √13 , + √13 ),

Point 4: (− √13 , + √13 ), and the weighting value for each integration point is 1. The full discretization of Eq.(22) is written as  4  n−1/2  1 e n+1/2 n ∑  m (v − v ) = ∆t (pq Liq )n ,   iI iI 4   q=1 (23)  4    n+1 n n+1/2 ∑  (pq Vq )n+1/2 .   E j − E j = −∆t

(19)

q=1

In this new scheme, the internal energy is constant on the element, but pressures and densities are varied. The time discretization schemes of the momentum equation and the energy conservation equation are the same with the old scheme, shown in Eq.(23). In order to resolve the shock, artificial viscosity has to be included working with the pressure[17, 18]. In the LGM scheme, the artificial viscosity Q shown in Eq.(17) will be varied on the element just like pressures, where ∇·v, sound speed c should be computed on each Gauss point. Numerical examples using the LGM method are computed in this paper using the identical coefficients a1 = 0.06, a2 = 1.5 with the old scheme.

where r34 = r3 − r4 , z12 = z1 − z2 , z34 = z3 − z4 . From Eq.(19), we can see that the area of Gauss point is not equivalent to the subzonal area except for rectangles. The integral form of the energy conservation equation is ∫ ∫ d α ρer drdz = − p∇ · vrα drdz. (20) dt Ω

q=1

(

Ω

On an arbitrary element domain Ωe , Eq.(20) can be written into ∫ ∫ 1∫ 1 d ρedV = − p∇ · vrα Jdξdη. (21) dt Ωe −1 −1 5

Additionally, in the two-dimensional plane geometry, both the old scheme and the LGM method satisfy the geometrical conservation law[21] naturally when the timecentral discretization scheme is used, shown in Eq. (24). But, both schemes cannot satisfy the geometrical conservation law in the axisymmetrical geometry.

d dt

∫

∫ drdz =

Ωe

∫ ∇ · vdrdz =

Ωe

1 −1

∫

1

−1

ric finite element method, ( ∂NI ) 1 ∂r = 16J ∂NI [( ∂z ) [(1 − ξ)(z41 ) + (1 + ξ)z32 ]ξI (1 + ηI η) (−[(1 − ξ)r41 + (1 + ξ)r32 ]ξI (1 + ηI η) )] −[(1 − η)z21 + (1 + η)z34 ]ηI (1 + ξI ξ) + , [(1 − η)r21 + (1 + η)r34 ]ηI (1 + ξI ξ)

where ri j = ri − r j , zi j = zi − z j . When the mesh deforms in the hourglass pattern, combined with Eq.(25) and Eq.(26), we can get the divergence of velocity on the element on the Gauss point (ξ = 0, η = 0),

∇ · vJdξdη.

(24) The LGM scheme is described fully at the moment. We would like to emphasize that the LGM method differs from the subzonal pressure method in three following aspects: Firstly, the LGM method assumes the Lagrangian Gauss point mass, and the subzonal pressure method assumed the Lagrangian subzonal mass. Masses of Gauss points are different from subzonal masses on quadrilaterals except for some special conditions, shown in Eq.(19). Secondly, the LGM method solves the energy conservation equation using the Gauss quadrature when pressures of Gauss points are given by the equation of state, and the subzonal pressure method integrated the energy conservation equation on zones using zones’ pressures and subzonal pressures working with some coefficient. Thirdly, it is the most important that the LGM method solves the weak form of the momentum equation using the Petrov-Galerkin isoparametric finite element method, but the compatible hydrodynamics scheme with subzonal pressure method[1] solved the integral form of the momentum equation by the use of finite volume method.

∇ · v(ξ = 0, η = 0) = vI · ∇NI (ξ = 0, η = 0) = 0. (27) The mesh is deformed when the nodes move in the hourglass mode but the divergence of velocity is zero on the Gauss point (ξ = 0, η = 0) in the old scheme from Eq.(27). In other words, the volume or the internal energy of element do not change in hourglass mode. Since pressure depends on density and internal energy, the pressure field does not change[14]. Consequently, hourglass distortions without any restoring forces will grow without bound. When the quadrilateral element deforms in the hourglass pattern, for instance (vr1 , vr2 , vr3 , vr4 )= (1, −1, 1, −1), (vz1 , vz2 , vz3 , vz4 ) = (0, 0, 0, 0), we can get the divergences of velocity on the four Gauss-points: /√ /√ √ (−1 3, −1 3) : ∇ · v = (z2 − z4 )/(2 3), /√ /√ √ (+1 3, −1 3) : ∇ · v = (z1 − z3 )/(2 3), /√ /√ √ (28) (+1 3, +1 3) : ∇ · v = (z4 − z2 )/(2 3), /√ /√ √ (−1 3, +1 3) : ∇ · v = (z3 − z1 )/(2 3).

2.4. Analysis of the old scheme and LGM scheme about hourglass modes

It is obvious from Eq.(28) that the divergence of velocity on each Gauss point is nonzero in most cases in the LGM method when the element deforms in this hourglass mode, although the element volume does not change. It is easy to show that divergences of velocity on four Gauss-points are nonzero in most cases in the hourglass mode,(vr1 , vr2 , vr3 , vr4 ) = (0, 0, 0, 0), (vz1 , vz2 , vz3 , vz4 ) = (1, −1, 1, −1). Consequently the density and pressure of Gauss point will change to produce the restoring force to suppress hourglass modes. Therefore, we can say that the LGM scheme eliminates the hourglass modes on the 4-node quadrilateral element.

There are hourglass distortions with the old scheme depicted in section 2.2, but how about the LGM scheme? It is not easy to give the theoretical analysis because pressures are nonlinear functions of nodal displacements, so we will think about the problem conversely. What will happen if the element deforms in hourglass modes? On an arbitrary 4-node quadrilateral element, velocities of nodes in hourglass patterns are (v1I ) = (1, −1, 1, −1),(v2I ) = (1, −1, 1, −1), I = 1, · · · , 4. The divergence of velocity on an arbitrary element is ∇ · v = vI · ∇NI .

(26)

3. Numerical results

(25)

In this section, we will give the numerical results computed by the use of two schemes, the old scheme

where ∇NI can be derived using the 4-node isoparamet6

suppressing the hourglass mode with hourglass viscosity and the new scheme of LGM eliminating the hourglass mode by means of Lagrangian Gauss-point masses and associated pressures. We will compare two schemes in grid distortions and the conservation of total energy on problems involving the ideal gas. We should give the definition of the total energy of the staggered Lagrangian scheme before comparisons of the conservation of total energy. From the above descriptions, it is obvious that the velocity and the internal energy are defined at the different time and space style in the staggered Lagrangian scheme. The velocity is defined at the time tn+1/2 at the node, but the internal energy is defined at time tn at the element. Here, we use the approach proposed by Benson[16] to compute the total energy of two-dimensional system at time tn , Tn =

N ∑ 1 I=1

2

mI (vn−1/2 vn+1/2 + vn−1/2 vn+1/2 )+ rI rI zI zI

K ∑

Figure 3: Results comparison of Sod problem using the LGM scheme and old scheme at t=0.26.

problem, which has the exact solution: the shock front lies at the place of radial radius 0.6, and the density after the shock front is 16. there is isentropic heating ahead of the shock, ρ = 1 + t/R. Cylindrical Noh problem is computed on the Cartesian grid, which is 50×50 zones. Only a quarter of volume is modeled, and two boundaries on the r and z axes are reflective. Results of the LGM scheme and the old scheme are shown in Fig. 4. Densities of all nodes are plotted in Fig. 5, which are the average values of their neighbor cells or Gauss points. Results of LGM scheme are more symmetrical than those of the old scheme. The LGM scheme changes results a little, but result grids by the LGM scheme are less distorted comparing with those by the old scheme. Here, we do not talk about the wall heating of Noh problem, which has been an open problem for Lagrangian schemes since the inception.

mk enk ,

k=1

where N is the number of nodes, K is the number of elements, vn+1/2 is the r-axis velocity of node I at time rI tn+1/2 , vn+1/2 is the z-axis velocity of node I at time zI n+1/2 t , mI is the mass of node I which is the sum of one quarter mass of its neighbor elements, mk is the mass of element k, enk is the specific internal energy of element k at time tn . 3.1. Sod problem

old

new

Sod problem is one-dimensional physically, which has not apparent hourglass distortions without hourglass control. It is used to show that LGM method will not change the results of problems almost where apparent hourglass distortions do not exist. The left state is a high pressure gas characterized by (ρL , pL , vL ) = (1, 1, 0), and the right state is a low pressure gas characterized by (ρR , pR , vR ) = (0.125, 0.1, 0). Two materials are both described by the ideal gas equation of state with γ = 75 . All boundaries are reflective. Results of the two schemes are shown in Fig.3. The difference between them is too small to be discerned from the figure, so the two schemes are nearly identical in Sod problem, which has not apparent hourglass distortions without hourglass control.

0.3

R(cm)

0.2

0.1

0

-0.2

-0.1

0

0.1

0.2

Z(cm)

Figure 4: Grid zoom of cylindrical Noh problem using the LGM scheme and old scheme on Cartesian grid at time t=0.6. For the convenience of comparison, grid of LGM scheme is reversed around r axis. black grid: result of old scheme, red grid: result of LGM scheme

We compute the spherical Noh problem also on the one-dimensional grid, which has 100 radial cells. The LGM scheme preserves the symmetry, moreover it improves the energy conservation to a great extent for the two-dimensional axisymmetrical geometry, although it cannot preserve the energy conservation because of the area-weighted scheme of the momentum equation,

3.2. Noh problem Noh problem[19] was proposed on shock errors involving an ideal gas(γ = 53 ). The initial state of the ideal gas: the internal energy is zero, the density is 1, the gas is given an initial unit inward velocity in the radial direction[16]. Here we compute the cylindrical Noh 7

shown in Fig.7.

1

18 16

0.995

exact old new

14

totener

12

new old

0.99

0.985

Density

10

0.98

8 6

0.975

4

0.97

2

0

0.5

1

1.5

2

2.5

3

3.5

time 0 0.0

0.2

0.4

0.6

Figure 7: Comparison of total energy conservation of spherical Noh problem computed by the LGM scheme and old scheme. black line: result of old scheme, red line: result of LGM scheme

0.8

Radius

3.3. Sedov problem

Figure 5: Results comparison of cylindrical Noh problem using the LGM scheme and old scheme on Cartesian grid at time t=0.6. black circles: result of old scheme, red circles: result of LGM scheme

The Sedov blast wave problem is used to demonstrate the symmetry of the scheme on the asymmetrical grid. It is run in the two-dimensional plane geometry. The initial state consists of a square grid with an edge of length 1.1. The space is divided into 45×45 square zones[20]. The initial velocity is zero and the initial density is unity. The specific internal energy is zero except for the cell containing the origin where it is 409.7. The self-similar solution is: at time t=1, the radius of diverging shock is unity, and the peak density is 6 for the ideal gas with γ = 1.4. Results of the LGM scheme and the old scheme are shown in Fig.8, in which the nodes’ densities are the average values of their neighbor elements or Gauss points. The peak density of the LGM scheme is closer to 6 than that of the old scheme, and the shock place is almost unity. The LGM scheme is much better than the old scheme in the conservation of total energy, shown in Fig.10. The LGM scheme is almost conservative for Sedov problem. There is asymmetry phenomena about mesh near the axis shown in Fig. 9, which is related to the artificial viscosity. So, we will improve the artificial viscosity in the future.

1 new old

totener

0.995

0.99

0.985

0.98

0

0.2

0.4

0.6

0.8

time

Figure 6: Comparison of total energy conservation of cylindrical Noh problem computed by the LGM method and the old scheme. black line: result of old scheme, red grid: result of LGM scheme

3.4. Dukowicz-Meltz problem The Dukowitz-Meltz problem[2] contains the refraction of shock which interacts with an inclined interface 8

1

new old

6

0.995

totener

new old

0.99

4

Density

0.985

0.98

2

0

0.5

1

1.5

2

time

Figure 10: Energy conservation of cylindrical Sedov problem at time t=1 using the LGM scheme and old scheme. black line: results of old scheme, red line: results of LGM scheme.

0 0

1

2

Radius

possessing a density discontinuity. In Dukowicz-Meltz problem, a piston moves from the left to the right at the velocity of 1.48, the density of left material is unity, and the density of the right material is 1.5. Two materials are ideal gases with γ = 1.4, and their initial internal energy per volume are 2.5. The right boundary is at an angle of 60◦ relative to the vertical. The grid of the right material is slanted but uniformly spaced at an angle of 60◦ . The top and the bottom boundaries are reflective. The initial grid is shown in Fig. 11, where the inner boundary distinguishing two materials does not exist in the computation. The grid of the left material is 36×30, and the grid of right material is 40×30.

Figure 8: Density of Sedov problem at time t=1 using the LGM scheme and old scheme, density of all nodes are plotted. black line: results of old scheme, red line: results of LGM scheme.

1

R(cm)

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

Figure 11: Initial grid of Dukowicz-Meltz problem

1.2

Z(cm)

The grid is noticeably distorted in the vicinity of lower part of the interface when the old scheme is used. Almost all the grid tangles disappear when the LGM scheme is used shown in Fig.13, and the internal energy contour has changed slig htly (Fig. 12). The re-

Figure 9: Grids of Sedov problem at time t=1 using the LGM scheme.

9

time = 1.30

time = 1.30

1.5

1.5

1

1

old

old

Level ENERGY 9 8 7 6 5 4 3 2 1

0

-0.5

420 400 380 360 340 320 300 280 260

Level den

0.5

R(cm)

R(cm)

0.5

new

0

-0.5

-1

-1.5

9 8 7 6 5 4 3 2 1

-1

2

3

4

-1.5

5

4.4 4 3.6 3.2 2.8 2.4 2 1.6 1.2

new

2

3

Z(cm)

4

5

Z(cm)

Figure 12: Internal energy contours and pressure contours of Dukowicz-Meltz problem at time t=1.3

time = 1.30

0.4 time = 1.30

1.5

1

old

0.2

old

R(cm)

R(cm)

0.5

0

0

-0.5 -0.2 new

-1

-1.5

new

2

3

4

-0.4

5

Z(cm)

2

2.2

2.4

2.6

2.8

Z(cm)

Figure 13: Grid(left) and zoom of bottom boundary(right) of Dukowicz-Meltz problem at time t=1.3

10

sult shows that the LGM method can suppress the nonphysical distortions, which is more robust than the old scheme.

1

old new

3.5. Saltzman problem

0.9

totener

The Saltzman problem is a well-known problem to evaluate the Lagrangian scheme. This problem is the propagation of a strong one-dimensional shock wave. Many Lagrangian codes have difficulty with this test problem, particularly since deviations from the expected one-dimensional behavior are so easy to detect[2]. In Saltzman problem, a piston moves with unit velocity from the left to the right, sending a shock across the grid that is slightly slanted with the vertical. The z coordinates of nodes is: z(i,j)=(i-1)*0.01+(11-j)*sin(π (i1)/100)*0.01, and the initial mesh is uniformly spaced by 10 in the r direction. The gas is described by an ideal gas equation of state with γ = 53 . The right, top and bottom boundaries are reflective. The computation is carried out in the two-dimensional plane geometry. The exact solution is: at time t=0.925, the shock is at the place of z=0.95. In this problem, almost all the grid distortions are spurious because it is a one-dimensional problem physically, so it is a good example to test the robustness of the scheme.

0.8

0.7

0

0.1

0.2

0.3

0.4

0.5

time

Figure 17: Comparison of the total energy conservation of t=2.0

Results in Fig.16 show that the LGM scheme preserves the symmetry much better than the old scheme. The shock of results by the LGM scheme is almost a circle at time t=0.5, comparing with the shock by the old scheme which has lost the circle shape at time t=0.2. The LGM scheme preserves the conservation of total energy much better than the old scheme shown in Fig.17. The result of the old scheme has lost more than 30% of total energy at time t=0.5, so it is unbelievable. 3.7. Three material instability problem Three material instability problem is often used to assess the robustness of Lagrangian scheme, which involves the physical vortices. The initial condition is shown in Fig. 18. The grids of the new scheme do not tangle at the time 2.5 although some of them are largedeformation, but those of the old scheme tangles at early time, and there are apparent tangles at time 2.5, Fig. 19. Comparing with the old scheme, the new scheme improves the energy conservation, Fig. 20.

Figure 14: Initial grid of Saltzman problem

From Fig. 15, almost all the grid tangles in the result of the old scheme disappear when the LGM scheme is used. The shock place is closer to z=0.95, and the LGM scheme keeps the one-dimensional property much better than the old scheme. This test example shows that the LGM scheme is better to suppress the spurious distortions on the high-aspect-ratio grid comparing with the old scheme.

3

3.6. Diverging blast wave from a point similar to Sedov problem

ρ3 =0.125 p3 =0.1 γ3 =1.6

ρ1 =1

2

R

p1 =1 γ1 =1.5

ρ2 =1

1

Here we use the grid and the ideal gas of the Sedov problem, only that the cell of high internal energy moves from the origin to the center of the region. A diverging blast wave from a point is produced, which is used to test the properties of symmetry and conservation of two schemes here.

p2 =0.1 γ2 =1.4

0

7

6

5

4

3

2

1

0

Z

Figure 18: Initial condition of three material instability problem

11

0.1

time = 0.925

0.1

time = 0.925

old

old

0.05

0.05

10 9 8 7 6 5 4 3 2 1 new

0

-0.05

-0.1 0.92 0.94 0.96 0.98

1

den 20.8 19.6889 18.5778 17.4667 16.3556 15.2444 14.1333 13.0222 11.9111 10.8

R(cm)

R(cm)

Level

0

-0.05

new

-0.1 0.92 0.94 0.96 0.98

1.02 1.04

Z(cm)

1

1.02 1.04

Z(cm)

Figure 15: Density contour and grid of saltzman problem at t=0.925

time = 0.2

time = 0.5

1

1

R(cm)

old

P(GPa)

0.5

42 38 34 30 26 22 18 14 10 6 2

0

new

-0.5

R(cm)

old

0.5

16 14 12 10 8 6 4 2

0

new

-0.5

-1 0

P(GPa)

-1 0.5

1

1.5

0

Z(cm)

0.5

1

Z(cm)

Figure 16: Diverging shocks from a point at the time of t=0.2 and t=0.5

12

1.5

time = 2.5

3 2.5 2 1.5

old

R(cm)

1 0.5 0

-0.5 -1 -1.5

new

-2 -2.5 -3

7

6.5

6

5.5

5

4.5

4

3.5

3

2.5

2

1.5

1

0.5

Z(cm) Figure 19: Grid comparison between the old scheme and the new scheme

13

0

Acknowledgements 1

We thank Shuanghu Wang, Weidong Shen and Yibing Chen for simulating discussions. The authors are supported by NSFC(Grant No.11271051 and No. 11171037, No. 91130021) and CAEP foundation(2012A0202010).

old new

0.998

totener

0.996

References

0.994

[1] E. J. Caramana and M. J. Shashkov, Elimination of Artificial Grid Distortion and Hourglass-type Motions by Means of Lagrangian Subzonal Masses and Pressures, Journal of Computational Physics 142 (1998) 521-561. [2] John K. Dukowicz and Bertrand J. A. Meltz, Vorticity Errors in Multidimensional Lagrangian Codes, Journal of Compuational Physics 99 (1992) 115-134 . [3] G. Maenchen and S. Sack, The Tensor Code, Methods in Computation Physics, Vol. 3, Fundamental methods in Hydrodynamics, Academic Press, 181-210, 1964. [4] W. D. Shultz, Two-dimensional Lagrangian Hydrodynamic Difference Schemes, Methods in Computation Physics, Vol. 3, Fundamental methods in Hydrodynamics, Academic Press, 1-45, 1964. [5] Mikhail Shashkov and Burton Wendroff, A Composite Scheme for Gas Dynamics in Lagrangian Coordinates, Journal of Computational Physics 150 (1999) 502-517. [6] Mark L. Wilkins, Use of Artificial Viscosity in Multidimensional Fluid DYnamics Calculations, Journal of Computational Physics, 36 (1980) 281-303. [7] J. O. Hallquist, Lectures Notes for Bulk and Hourglass Viscosity in Wave Propagation Codes, UCRL-89156, 1983. [8] D. P. Flanagan and T. Belytschko, A Uniform Strain Hexahedron and Quadrilateral with Orothogonal Hourglass Control, Internal. J. Numer. Methos. Eng. 17 (1981) 679-706. [9] Ted Belytschko, Wing Kam Liu, Brian Moran, Nonlinear Finite Elments for Continua and Structures, John WIley & Sons Ltd., 2000. [10] J. C. Campbell and M. J. Shashkov, A Tensor Artificial Viscosity Using a Mimetic Finite Difference Algorithm, Journal of Compuational Physics 172 (2001) 739-765. [11] A. J. Barlow, A Compatible Finite Element Multi-material ALE Hydrodynamic Algorithm, Internation Journal for Numerical Methods in Fluids 56 (2008) 953-964. [12] Raphael Loubere, Pierre-Henri Maire and Pavel Vachal, A Second-order Compatible Staggered Lagrangian Hydrodynamics Scheme Using a Cell-centered Multidimensional Approximate Riemann Solver, Procedia Computer Science 1 (2010) 1931-1939. [13] M. L. Wilkins, Calculation of Elastic-Plastic Flow, Methods in Computation Physics, Vol. 3, Fundamental methods in Hydrodynamics, Academic Press, 211-264, 1964. [14] L. G. Margolin and J. J. Pyun, A Method for Treating Hourglass Patterns, LA-UR-87-439, 1987. [15] L. G. Margolin, D. E. Burton, W. P. Crowley and B. C. Trent, Computers Simulation of Nuclear Weapons Effects, UCRL98438, 1988. [16] D. J. Benson, Computational Methods in Lagrangian and Eulerian Hydrocodes, Computer Methods in Applied Mechanics and Engineering 99 (1992) 235-394. [17] J. von Neumann and R. D. Richtmyer, A Method for the Numerical Calculation of Hydrodynamic Shocks, Journal of Applied Physics 21 (1950) 232-237.

0.992

0

1

2

3

4

time

Figure 20: Energy conservation comparison of the old scheme with the new scheme

4. Concluding remarks This paper gives a method to suppress the hourglass motions on quadrilaterals based on the traditional Lagrangian scheme for compressible fluid flows, which is abbreviated to LGM for convenience. The LGM method has no other coefficients except for the coefficients of the artificial viscosity. The core of the LGM method is the assumption that the mass of Gauss point is Lagrangian. When densities of Gauss points are known, their pressures can be derived by the use of equations of state. The 4-node quadrilateral element has four pressures if four-points Gauss quadrature is used. Numerical examples show that the LGM method eliminates the hourglass-type distortions and spurious distortions on the high-aspect-ratio grid, at the mean time it improves the conservation of total energy to a great extent. The LGM method is almost conservative in the twodimensional plane geometry. In the two-dimensional axisymmetrical geometry, it does improve the conservation of total energy to a great extent, although it is not conservative because the area-weight scheme is used to the momentum equation for symmetry. It is simple to implement the LGM method into the existing 2D finite element hydrocodes on quadrilaterals, moreover it can be extended to the 3D finite element hydrocodes on hexahedrons without theoretical difficulties. The LGM method works well for compressible fluid flows with ideal gas equation of state, but it is unknown whether it works well or not for materials with complex equations of state and elastic-plastic materials. We will work on to make that clear in the future. 14

[18] Rolf Landshoff, A Numerical Method for Treating Fluid Flow in the Presence of Shocks, LA-1930, 1955. [19] W. F. Noh, Errors for Calculations of Strong Shocks Using an Artificial Viscosity and an Artificial Heat Flux, UCRL-53669, 1985. [20] E. J. Caramana, D. E. Burton, M. J. Shashkov and P. P. Whalen, The Construction of Compatible Hydrodynamics Algorithms Utilizing Conservation of Total Energy, Journal of Computaitonal Physics 146 (1998) 227-262. [21] Pierre-Henri Maire, A High-order Cell-centered Lagrangian Scheme for Compressible Fluid Flows in Two-dimensional Cylindrical Geometry, Journal of Computational Physics 228 (2009) 6882-6915.

15

Highlights （1） （2） （3） （4） （5）

Hourglass distortions are analyzed on quadrilaterals by the finite element method. LGM method is given for hourglass distortions assuming Lagrangian Gauss-point mass. LGM method is similar to the subzonal pressure method, but they are different. LGM method can suppress hourglass distortions on quadrilaterals. LGM method improves the total energy conservation to a great extent.