A three-dimensional finite element arbitrary Lagrangian–Eulerian method for shock hydrodynamics on unstructured grids

Computers & Fluids 92 (2014) 172–187

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

A three-dimensional finite element arbitrary Lagrangian–Eulerian method for shock hydrodynamics on unstructured grids q J. Waltz a,⇑, N.R. Morgan a, T.R. Canfield b, M.R.J. Charest a,1, L.D. Risinger c, J.G. Wohlbier c a

Computational Physics Division, Los Alamos National Laboratory, Los Alamos, NM, United States Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, United States c Computational and Computer Sciences Division, Los Alamos National Laboratory, Los Alamos, NM, United States b

a r t i c l e

i n f o

Article history: Received 15 July 2013 Received in revised form 13 November 2013 Accepted 21 December 2013 Available online 29 December 2013 Keywords: Arbitrary Lagrangian–Eulerian Shock hydrodynamics Finite element method

a b s t r a c t We present a three-dimensional (3D) finite element (FE) arbitrary Lagrangian–Eulerian (ALE) method for shock hydrodynamics on unstructured grids. The method is based on an FE Eulerian Godunov scheme for linear tetrahedra that has been extended to include mesh motion in an unsplit, flux-conservative formulation. The proposed method eliminates the splitting errors present in traditional Lagrange-plus-remap methods that occur during the remap phase. Unlike typical unsplit approaches, the mesh velocity is not determined by boundary motion but is instead based on the local fluid velocity. Smoothing operations are then applied to the mesh velocity to avoid mesh tangling. This approach allows the mesh to follow the fluid motion in a robust manner and leverage one of the primary advantages of Lagrangian schemes for shock hydrodynamics, namely that the resolution follows the flow. An approximate Riemann solver is used to calculate fluxes in the co-moving frame of the mesh. Results for a number of standard test problems are presented for 3D meshes of up to 107 tetrahedra. Global convergence rates of 0:8– 1:0 are observed for shock dominated flows and 1:9 for smooth flows. We also demonstrate that the method satisfies the discrete geometric conservation law to truncation error, conserves total energy to machine precision, and preserves symmetry. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Since its inception, the arbitrary Lagrangian–Eulerian (ALE) method [1] has been applied to a wide range of physics and engineering problems. The method attempts to combine the advantages of Eulerian and Lagrangian algorithms to accurately capture large distortions of the continuum while maintaining clear representations of free surfaces and/or material/fluid interfaces. ALE methods have proved especially attractive for shock hydrodynamics problems, which involve the calculation of material response to shocks and related phenomena. The traditional ALE approach for shock hydrodynamics problems is commonly referred to as the Lagrange-plus-remap method. In this formulation, a Lagrange algorithm is used to calculate the stress terms, update the solution variables in the hydrodynamic equations, and move the mesh points at the local fluid velocity (e.g. [2–8]). This Lagrange step is followed by a ‘‘remap’’ step in

q This report was supported by the U.S. Department of Energy through the LANL/ LDRD Program. LANL report no. LA-UR-13-25216. ⇑ Corresponding author. Address. MS B259, Los Alamos National Laboratory, Los Alamos, NM 87545, United States. E-mail address: [email protected] (J. Waltz). 1 Nicholas C. Metropolis Postdoctoral Fellow.

0045-7930/$ - see front matter Ó 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2013.12.021

which the mesh coordinates are modified in some fashion to improve mesh quality [9–12]. Alternatively, a new mesh consisting of higher quality cells is generated [13]. In either case, the solution is transferred from the post-Lagrange mesh to the post-remap mesh using volume or surface integration, i.e. using mesh-mesh intersections or advection. The primary advantage of the Lagrange-plus-remap method is that it allows accurate tracking of discontinuities and material interfaces. The method is also naturally adaptive because the mesh follows flow features, which is of particular benefit for applications that involve significant deformation and/or displacement. Because of these benefits, the method has enjoyed widespread popularity and forms the basis of numerous production codes used to model shock hydrodynamics phenomena [14–18]. More recent works include those of Galera et al. [19], Berndt et al. [20], Jia et al. [21], and Lopez Ortega and Scovazzi [22]. A different class of ALE methods has been applied to various engineering applications such as store separation [23,24], turbomachinery analysis [25], fluid–structure interaction [26,27], and aeroelasticity [28–30]. This class of methods uses an unsplit formulation in which advection due to mesh motion is solved simultaneously with advection due to fluid motion. The mesh velocity is typically specified at a subset of the boundary nodes, e.g. from structural deformation or at the surface of an object moving under

173

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

the influence of gravity and aerodynamic drag. The mesh velocity in the interior of the domain is calculated by smoothing the mesh velocity away from the boundary. Typical smoothing approaches include Laplacian-type algorithms [31], which operate directly on the mesh velocity; and force-based methods [32,33], which calculate equilibrium mesh coordinates and then determine the mesh velocity a posteriori from the mesh displacements. The aforementioned quality-based methods used in Lagrange-plus-remap methods can be applied in a similar manner [34]. Unsplit ALE formulations offer a number of significant advantages over the traditional Lagrange-plus-remap approach: 1. Temporal errors due to operator splitting of the stress and advection terms are eliminated. These errors are either firstor second-order depending on whether Lie or Strang splitting is used [35]. 2. Computational cost is reduced because advection terms are trivially calculated simultaneously with the stress terms; in typical Lagrange-plus-remap methods, e.g. those described in [13–18,21], these terms are calculated via independent passes over the mesh. 3. The advection and stress terms are discretized in a consistent manner, which eliminates the possibility that the two sets of terms will converge at different rates and/or have different truncation errors that lead to e.g. conservation errors, as has been observed with inconsistent discretizations of the momentum and energy equations in Lagrange formulations [5,18].2 4. As will be shown, the discrete system of equations is written in a flux-conservative form. Per the Lax–Wendroff theorem, this ensures that if the numerical method converges, then it will converge to a weak solution and any shocks will be in the correct location [36]. In contrast many Lagrange-plus-remap methods do not conserve total energy by construction and require corrections to do so. 5. The formulation naturally incorporates approximate Riemann solvers and eliminates the need for artificial viscosity treatments [37,38]. This allows direct application of Godunov concepts to Lagrangian and ALE methods, which has been a topic of some recent two-dimensional (2D) work [39,40]. While unsplit ALE formulations have generally been applied to engineering problems with moving boundaries, a small number of researchers have applied this class of methods to shock hydrodynamics. Luo et al. [41] presented an unsplit ALE formulation for the solution of 2D, two-material compressible flows. These authors made use of a node-centered Godunov-like finite volume (FV) spatial discretization on unstructured triangular meshes. The material interface was treated in a Lagrangian manner, with the interface velocity calculated from the approximate solution to a two-material Riemann problem. Although this approach was effective for the class of problems under consideration, the assumption that material interfaces can be treated in a pure Lagrangian manner will hold only when the interfaces remain intact and do not experience significant vorticity and/or shear. The method also was demonstrated for 2D cases only. Kershaw et al. [42] presented an unsplit ALE formulation for the solution of three-dimensional (3D) hydrodynamics problems on unstructured mixed-element meshes. A Godunov-like method was used to calculate fluxes across element boundaries. Their approach was subsequently applied to problems involving radiation-hydrodynamics by Shestakov et al. [43]. However, because 2

None of the implementations referenced above employ consistent discretizations for the Lagrange and advection terms: the Lagrange step typically assumes piecewise constant values within a mesh cell but advection is typically calculated with a piecewise linear reconstruction.

this approach was based on a discontinuous cell-centered finite element (FE) method, the calculation of a unique nodal velocity from the discontinuous fluid velocity was not straightforward and the pure Lagrange limit could not be recovered. A novel 3D unsplit ALE method is proposed herein to model general shock hydrodynamics problems. The proposed method makes use of an edge-based FE discretization for unstructured meshes composed of tetrahedral cells [44]. Approximate Riemann solvers are used to calculate fluxes within the co-moving frame of the mesh. Since the proposed approach is based on a continuous node-centered FE method, the mesh velocity can be set exactly to the fluid velocity with no ambiguity as to how the mesh points should be moved. Smoothing of the mesh velocity is then applied to prevent mesh tangling while simultaneously allowing the mesh to follow the flow as much as possible. Furthermore, because the scheme is purely node-centered, the various issues historically associated with Lagrangian shock hydrodynamics calculations on triangular and tetrahedral meshes are avoided entirely (see for example [8] for a discussion of some of these issues). Compared to the Lagrange-plus-remap method, this mesh velocity treatment has the additional advantage that mesh motion will only occur in (and possibly adjacent to) regions where the fluid velocity is non-zero. This behavior helps preserve mesh features while still adapting the mesh to the flow. In many shock hydrodynamics problems, the initial mesh size can vary drastically and anisotropically throughout the computational domain if there are large changes in material density or other disparate physical scales. When modeling an Inertial Confinement Fusion (ICF) target, for example, the capsule material may require a much finer radial resolution than the fuel, yet both may have the same angular resolution. Similar situations arise in the solution of viscous wallbounded flows where high aspect ratio cells are needed near solid boundaries to accurately represent boundary layers. Mesh-qualitybased approaches that attempt to equalize cell volumes, aspect ratios, interior angles, and/or edge lengths can destroy any such anisotropy in the mesh, degrading the quality of the solution. Additional and often ad hoc controls are required [16] to overcome this problem in traditional ALE formulations. These controls also can be resolution dependent, with the consequence that a given set of inputs may not allow reliable convergence under mesh refinement. In contrast the method described herein converges robustly under mesh refinement at the expected order of accuracy. The method also is shown to satisfy the geometric conservation law (GCL) to truncation error. This paper presents the basic features of the unsplit ALE formulation, with mesh motion based on the fluid velocity, as implemented in the 3D unstructured mesh code CHICOMA [45]. In what follows we describe the FE formulation and its application to the ALE form of the hydrodynamic equations. The time integration scheme and the procedure used to calculate the mesh velocity are subsequently presented. A number of test problems are used to verify and demonstrate the algorithm. Multiple material treatments, adaptive mesh refinement (AMR) [46], and high-order extensions [47] will be addressed in future works. Parallelization has been performed with OpenMP using the techniques described in [44,48,49], but is not specifically addressed in this paper. Table 1 Mesh sizes and the corresponding L1 norms of the error between the exact and predicted density for the GCL problem. Mesh

Points

Tetrahedra

Ave. dx

L1 Error

0

7294

34,181

0:840

3:13  108

1

51,794

273,448

0:414

6:36  109

2

389,139

2,187,584

0:206

1:18  109

3

2,014,277

17,500,672

0:103

2:38  1010

174

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

defined by wi , which is related to the Eulerian (lab) reference frame via

  @ @ @ ¼ þ wi @t w @t @xi

ð4Þ

and the Lagrangian (fluid) reference frame via

  @ D @ ¼  ðui  wi Þ @t w Dt @xi

ð5Þ

where

Fig. 1. Initial (top) and final (bottom) mesh velocity distribution for the GCL test problem on the coarsest mesh.

-7

D @ @ ¼ þ ui Dt @t @xi

ð6Þ

is the Lagrangian time derivative. The Eulerian limit is achieved when wi ¼ 0 and the Lagrangian limit is achieved when wi ¼ ui . In what follows we use the condensed notation

  1 @ VU @F j þ ¼0 V @t @xj

Density m=2

where

0

-7.5

q

ð7Þ

1

B C U ¼ @ qui A qE -8

is the vector of conserved unknowns and

Log (L1)

0

qðuj  wj Þ

1

B C F j ¼ @ qui ðuj  wj Þ þ dij p A qEðuj  wj Þ þ quj

-8.5

ð9Þ

is the flux; we also drop the subscript w from the time derivative for brevity of notation. The above system of conservation laws is closed by an equation of state of the form p ¼ pðq; eÞ, where e is the specific internal energy. Although the numerical method described herein is designed to support arbitrary analytic or tabular equations of state, the test problems are restricted to the ideal gas law p ¼ qeðc  1Þ, where c is the ratio of specific heats.

-9

-9.5

-10 -1

ð8Þ

3. Numerical method -0.9

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0

3.1. Spatial discretization

Log (dx) Fig. 2. L1 norms of the error between the exact and predicted density as a function of mesh resolution for the GCL test problem, including reference second-order slope.

The system (7) is discretized with an edge-based FE method for linear tetrahedra [44], where the unknowns are the conserved solution quantities located at the nodes of the mesh. The solution at any point xi within a tetrahedron Xh is approximated as

2. Governing equations

Uðxi Þ ¼

X v N ðxi ÞU v

ð10Þ

v 2Xh

The equations of interest for this work are the compressible hydrodynamic (Euler) equations. Specifically we seek to solve the flux-conservative ALE form of the hydrodynamic equations, written as

  1 @ðqVÞ @ þ ½qðui  wi Þ ¼ 0 V @t w @xi    @p 1 @ðqui VÞ @  þ qui ðuj  wj Þ þ ¼ 0 V @t @xj @xi w   1 @ðqEVÞ @ @ ðpui Þ þ ½qEðui  wi Þ þ ¼0 V @t @xi @xi w

ð1Þ ð2Þ ð3Þ

where, in the usual notation, q is the density, ui is the fluid velocity, wi is the mesh velocity, E is the specific total energy, p is the pressure, V is the volume of a fluid element, and t is time. The subscript w on the time derivatives is used to indicate that the equations have been transformed to the co-moving frame

Fig. 3. Intermediate (top) and final (bottom) surface mesh and density contours for the shock tube problem on the coarsest mesh.

175

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

-2

Density m=1

-2.2

Log (L1)

-2.4

-2.6

Fig. 4. Intermediate (top) and final (bottom) surface mesh and mesh velocity contours for the shock tube problem on the coarsest mesh.

-2.8 Exact Mesh 0 Mesh 1 Mesh 2 Mesh 3

1

-3

0.8

-3.2 -1

-0.9

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0

Density

Log (dx) Fig. 6. L1 norm of the error between the exact and predicted density as a function of mesh resolution for the shock tube test problem, including reference first-order slope.

0.6

0.4

0.0064

0.0062 0.2

0.006

0.0058 0 0

20

40

60

80

100

0.0056

Fig. 5. Predicted and exact density distributions along the shock tube with increasing mesh resolution.

L1

x 0.0054

0.0052 Table 2 Mesh sizes and the corresponding L1 norms of the error between the exact and predicted density for the shock tube problem. Mesh

Points

Tetrahedra

Ave. dx

L1 Error

0

7294

34,181

0:840

7:72  103

1

51,794

273,448

0:414

4:39  103

2

389,139

2,187,584

0:206

2:46  103

3

2,014,277

17,500,672

0:103

1:49  103

0.005

0.0048

0.0046

0.0044 0.01

0.1

1

10

Log (e_w) where v is a nodal index and U v are the nodal solution values. N v ðxi Þ is a linear FE basis function for node v on element Xh and evaluated at location xi . The semi-discrete version of (7) for node v results from a Galerkin approximation and is written for an interior point3 as: 3 Points on the boundary require additional terms that are detailed in A; those terms have been omitted in this section for brevity.

Fig. 7. L1 norm of the error between the exact and predicted density as a function of w for the shock tube problem.

 v dVU v X X vw d VU ¼ þ Dj F vj þ F w þ Djv w F jv w ¼ 0 j dt dt v w2v v w2v  v X vw d VU ! ¼ Dj F vj w dt v w2v

ð11Þ

176

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

Table 3 Mesh sizes, initial source values, and the corresponding L1 norms of the error between the exact and predicted density for the Sedov problem. Mesh

Points

Tetrahedra

Ave. dx

eS

ps 3

Density L1 3

0

11,304

59,164

0:053

7:25  10

4:86  10

9:17  102

1

65,958

362,363

0:029

5:42  104

3:63  104

5:43  102

2

505,724

2,898,904

0:014

3:46  105

2:32  105

2:38  102

where v and w are nodal indices, vw is the edge connecting nodes v and w, and the sum is over edges containing node v. A lumped-mass approximation is applied to the first term in (7), so that the factor of 1=V on the time derivative in (7) cancels. The quantity Djv w is calculated from the FE basis functions as

Dvj w ¼

 Z

1 X @Nw @Nv w Nv  N dXh 2 X 2v w Xh @xj @xj

-1 Density m=1

-1.1

ð12Þ

-1.2

h

where the sum is over elements connected to edge vw. On each element, Djv w is effectively the area-normal to the dual face separating nodes v and w. Furthermore, because Djv w is antisymmetric, the scheme is conservative to machine precision by construction. A full derivation of the edge-based formulation and an explanation of its meaning is given in Appendix A.

Log (L1)

-1.3

-1.4

-1.5

3.2. Approximate Riemann problem -1.6

The numerical flux, F vj w in (11), is evaluated at the midpoint of each edge via the solution to a local Riemann problem. Both Rusanov [50] and HLLC [51,52] approximate Riemann solvers have been implemented [45,53]. The extension of the approximate Riemann solvers to the ALE formulation requires two modifications. The first is the addition of the mesh velocity to the flux in (9). The second is a modification of the wave speeds to account for the moving reference frame.

-1.7

-1.8 -1.9

-1.8

-1.7

-1.6

-1.5

-1.4

-1.3

-1.2

Log (dx) Fig. 9. L1 norm of the error between the exact and predicted density as a function of mesh resolution for the Sedov problem, including reference first-order slope.

4

Exact Mesh 0 Mesh 1 Mesh 2

In the Rusanov scheme, the numerical flux increment on an edge is written as

3.5

    Div w F iv w ¼ Div w F vi þ F w þ jDiv w jkv w U w  U v i

3

ð13Þ

where the wave speed k is calculated as

Density

2.5

kv w ¼ maxðkv ; kw Þ

ð14Þ

with

2

  Div w þ cv kv ¼ uvi  wvi jDv w j

ð15Þ

i

1.5

cv is the soundspeed at point v. For the ALE HLLC scheme, the left- and right-state wavespeeds are given by

1

sv ¼ minðqv  cv ; qw  cw Þ

0.5

ð16Þ

and 0 0

0.2

0.4

0.6

0.8

1

1.2

x Fig. 8. Predicted and exact density distributions for the Sedov problem with increasing mesh resolution.

sw ¼ maxðqv þ cv ; qw þ cw Þ

ð17Þ

where

  qv ¼ ni uvi  wvi ;

  w qw ¼ ni uw i  wi

ð18Þ

177

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

4 Exact Mesh 2 ALE Mesh 2 EUL

3.5

3

Density

2.5

2

1.5

1

0.5

0 0

0.2

0.4

0.6

0.8

1

1.2

x Fig. 10. Comparison of ALE and Eulerian density distribution on Mesh 2 the Sedov problem. Fig. 12. Pressure contours for ALE (left) and Eulerian (right) calculations of the Sedov problem at intermediate and final times on Mesh 2.

0.0617893 Mesh 0 Mesh 1 Mesh 2

0.0617892 0.0617892 0.0617892

E

0.0617892 0.0617892 0.0617892 0.0617892 0.0617892 0.0617892

Fig. 13. Pressure contours for ALE (left) and Eulerian (right) calculations of the Sedov problem at final time near the shock on Mesh 2.

0.0617892 0

0.2

0.4

0.6

0.8

1

t Fig. 11. Total energy as a function of time for each calculation of the Sedov problem.

The intermediate left- and right-state conserved unknowns are given by

0

ðsv  qv Þqv

1

1 B v C U v ; ¼ v @ ðs  qv Þqv uvi þ ðp  pv Þni A s  sm v ðsv  qv Þqv E þ p sm  pv qv and ni is the outward unit normal on a dual mesh face defined as ni ¼ Div w =jDiv w j. Given the left and right wave speeds, the contact wave speed sm and pressure p are given by

sm ¼

q

q ðs  q Þ  qv qv ðsv  qv Þ þ pv  pw qw ðsw  qw Þ  qv ðsv  qv Þ

w w

w

w

and

0 1 ðsw  qw Þqw 1 B w C  w ¼ w @ ðs  qw Þqw uw i þ ðp  p Þni A s  sm ðsw  qw Þqw Ew þ p sm  pw qw

ð19Þ

U w;

ð20Þ

These values are used with the unmodified left- and right-state unknowns to compute the flux increment as

and

p ¼ qv ðqv  sv Þðqv  sm Þ þ pv

ð21Þ

ð22Þ

178

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

4.5

Eulerian

ALE

4

4

3.5

3.5

3.5

3

3

3

2.5

2.5

2.5

2

Density

4

Density

Density

4.5

4.5

Lagrangian

2

2

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

0

0 0

0.2 0.4 0.6 0.8 1

1.2 1.4 1.6 1.8

0 0

0.2 0.4 0.6 0.8 1

Radius

1.2 1.4 1.6 1.8

0

0.2 0.4 0.6 0.8 1

Radius

1.2 1.4 1.6 1.8

Radius

Fig. 14. Scatter plots of density as a function of radius for pseudo-2D Sedov problem: Lagrangian (left), ALE (middle), and Eulerian (right) variants.

8 ni F vi > > > > < n F v þ sv U v ;  U v  i i Div w F iv w ¼ 2jDiv w j   > ni F w þ sw U w;  U w > i > > : ni F w i

sv > 0 sv < 0 < sm sm < 0 < sw

Diminishing limiter function; in this work we have used Van Leer and minmod functions.

ð23Þ 3.4. Time integration

sw < 0

where the factor of 2 is needed to account for the factor of 1=2 in the definition of Div w .

X vw  v ;k  v ;0 VU ¼ VU  ak Dt Dj F jv w

3.3. Solution reconstruction



  1 1 ð1  kÞ/ðrv Þd1 þ ð1  kÞ/ v d2 4 r 

  1 1 ~ w ¼ uw  ð1 þ kÞ/ðrw Þd3 þ ð1  kÞ/ d2 u 4 rw

ð24Þ ð25Þ

with

@uv  d2 @xi w v d2 ¼ u  u @uw  d2 d3 ¼ 2xiv w @xi

d1 ¼ 2xiv w

ð26Þ

d2 ; d1

rw ¼

d2 d3

ð30Þ

where 1 < k < m is the current stage, m is the total number of stages, Dt is the time step, and

ak ¼

1 1þmk

ð28Þ

ð31Þ

m ¼ 1 yields the explicit Euler time marching scheme and m ¼ 2 is equivalent to the classical second-order Runge–Kutta method. The system is evolved from time step n to n þ 1 by taking U v ;0 ¼ U v ;n and U v ;nþ1 ¼ U v ;m . For linear problems the order of accuracy is m. For nonlinear problems the scheme is limited to second-order but, as shown by Waltz et al. [45], the use of m > 2 does lead to further error reductions. Three-stage integration was therefore used for the test results presented herein. The time step on an edge is calculated as

ð27Þ

Dt v w ¼

cN jxiv w j kv w

ð32Þ

where cN 6 1 is the Courant number. The global time step is the minimum of the Dt v w . An additional time step control based on the volume change at node v is calculated as

v where xiv w ¼ xw i  xi is the vector that defines the edge and

rv ¼

!v ;k1

v w2v

With either solver, the left- and right-state unknowns are computed from a limited MUSCL reconstruction on primitive fluid variables [54]. The Riemann solver is then applied to the reconstructed states. Given a field variable u with nodal values uv and uw , the reconstruction on edge vw is written as

~ v ¼ uv þ u

The time integration method is an explicit multi-stage scheme of the form

ð29Þ

The structured grid analogy would be to replace d1 with uv  uv 1 and d3 with uwþ1  uw . The scheme therefore corresponds to a piecewise linear reconstruction for k ¼ 1 and a piecewise parabolic reconstruction for k ¼ 1=3. The function / in the reconstruction formula can be any symmetric, second-order Total Variation

Dtv ;n ¼ cV Dt n1

 min V v ;n1 ; V v ;n jV v ;n1  V v ;n j

ð33Þ

where cV is a constant typically set to 0:01  0:05. Robustness of the method is enhanced by allowing a time step to restart with a smaller value if any of the following occur:

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

179

Fig. 15. Mesh and density contours at t ¼ 0:5 (left) and t ¼ 1:0 (right) for pseudo-2D Sedov problem: Lagrangian (top), ALE (middle), and Eulerian (bottom) variants.

Table 4 Mesh sizes and the corresponding L1 norms of the error between the exact and predicted velocity magnitude for the Taylor–Green vortex. Mesh

Points

Tetrahedra

Ave. dx

Velocity L1

0.011

4:43  103

0

24,530

75,912

1

148,378

607,296

0.0055

1:20  10

2

997,679

4,858,368

0.0028

3:35  104

the time step restart is activated very infrequently, but does allow calculations to remain stable that might otherwise not. 3.5. Mesh motion

3

The mesh coordinates are integrated in time in sequence with the solution vector as

xvi ;k ¼ xvi ;0 þ ak Dtwiv ;k1  Negative density or total energy.  Negative lumped mass, which indicates tangling.  A change of 10% or more in the Courant time step over a cycle. The time step restart is implemented by simply abandoning the intermediate multi-stage results and resetting k ¼ 1 in (30). The geometry i.e. Dvi w and M vL , also must be reset, or alternatively recalculated from xni . As a consequence of these additional controls, all simulations can be run with cN ¼ 1. If any errors are detected the time step will restart with a smaller Dt. Dt is then allowed to increase until it reaches the maximum allowed value. In practice,

ð34Þ

The edge terms Dvi w and the volume (lumped mass) are then recalculated from the xvi ;k at each stage of the time integration. Because the goal of the algorithm is to leverage the benefits of a Lagrangian formulation, the mesh velocity is initially set to the local fluid velocity, i.e. wki ¼ uki . However, this choice can lead to mesh tangling if the mesh undergoes large deformations or if the flow develops significant vorticity or shear. Therefore some modification of wki is necessary. In the present research we use a simple Laplacian-type smoothing approach. The system

180

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

where c1 and c2 are constants and -i is the fluid vorticity vector. This form allows regions of high vorticity to be treated in a more Eulerian-like manner, and regions of low vorticity to be treated in a more Lagrangian-like manner. In numerical experiments, w ¼ 0:1—1:0; c1 ¼ 0:25—0:5, and c2 ¼ 0:5—1:0 provided reliable and robust performance for a variety of problems. In general, flows with more vorticity require a smaller value of c1 and a larger value of c2 , and vice versa. wi is subject to the same boundary conditions as ui , although where a Dirichlet condition is applied to ui , one can choose either wi ¼ 0 for an inflow/outflow condition or wi ¼ ui for a free surface. The use of the tolerance w rather than a fixed number of iterations has two advantages. First, the required number of iterations will naturally be a function of the smoothness (or lack thereof) of the fluid velocity, which can vary significantly during unsteady problems. Second, the level of convergence will be the same regardless of mesh size, which is not true with a fixed number of iterations. Note that the limiting case w ! 0 will only correspond to the Eulerian limit if the fluid velocity is zero on the boundary of the computational domain. Otherwise a linear mesh velocity distribution is obtained as w ! 0.

-2 Velocity m=2

-2.2

-2.4

Log (L1)

-2.6

-2.8

-3

-3.2

-3.4

3.6. Full timestep -3.6 -2.7

-2.6

-2.5

-2.4

-2.3

-2.2

-2.1

-2

Log (dx) Fig. 16. L1 velocity error as a function of mesh resolution for the Taylor–Green problem, including reference second-order slope.

lr2 wki ¼ 0

ð35Þ

is solved iteratively to some tolerance w using a preconditioned Conjugate Gradient method [56]. The diffusivity l is given by



l ¼ c1 max 0; 1  c2

jj-i jj jj-i jj1

 ð36Þ

A summary of the algorithm for a single time step is now given, assuming m ¼ 1 in (31) for brevity and that a time step has been calculated from (32) and (33): 1. Compute the mesh velocity wi from Eq. (35). 2. Compute the FE gradient of primitive fluid variables q; ui , and e as well as p and c via Eq. (A.20). 3. On each edge: (a) Compute reconstructed solution values of q; ui ; e; p, and c at the edge mid-point from Eqs. (24) and (25).

Fig. 17. Mesh and velocity contours for the Taylor–Green problem at initial, intermediate, and final times on Mesh 0.

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

181

7. Compute M L ; Dvi w , and remaining geometric quantities at n þ 1. 8. Divide the unknowns VU at n þ 1 by the volume (lumped mass) at n þ 1 to obtain updated conserved solution values. 9. Update p and c at n þ 1 via an equation of state evaluation. In the case m > 1 the above procedure is simply repeated as necessary. Also note that for consistency in the Lagrangian limit, the mesh velocity must be recalculated at each stage of the multi-stage scheme. 4. Geometric conservation law 4.1. Statement of the geometric conservation law The geometric conservation law (GCL) governs the geometric parameters of numerical schemes that are designed for solving unsteady problems on moving meshes [29,30,55]. It states that a uniform flow field must be preserved under arbitrary mesh motion, and is a necessary and sufficient condition for a numerical scheme to preserve nonlinear stability [30]. Violation of this law can produce spurious oscillations and even unbounded behavior. Mathematically, the GCL implies that the discrete scheme must satisfy the relation v X vw dV ¼ Dj wvj þ ww j dt v w2v

Fig. 18. Close-up of Taylor–Green mesh at initial and final times illustrating the degree of mesh motion.

ð37Þ

at each point in the mesh, i.e. the deformation of each control volume must be calculated in a manner that is consistent with the motion of the mesh points. This relation results from the substitution of a uniform solution into (11), and the mesh velocity can be specified arbitrarily as long as (37) is satisfied. Detailed discussions of the GCL can be found in [29,30,55]. The primary concern for this work is whether the above relation is satisfied at the discrete level. 4.2. Numerical verification

Fig. 19. Dimensions and initial conditions for the triple point problem.

The GCL can be verified numerically by initializing a uniform field over a mesh, specifying a distribution for wi , evolving the system (11) for some period of time, and calculating the errors relative to the exact uniform solution. A verification study of this type was performed with four increasingly finer meshes on a cylindrical geometry with constant unit density and an initial prescribed mesh velocity in the z-direction given by 2

(b) Compute a numerical flux increment Dvj w F vj w using the reconstructed solution values and the Rusanov or HLLC approximate Riemann solver. 4. Sum Dvj w F vj w over edges connected to each point to form the right-hand-side of Eq. (11). 5. Update the unknowns VU to n þ 1 via Eq. (30). 6. Advance the mesh coordinates to n þ 1 via Eq. (34).

wz ¼ sin

 pz 100

ð38Þ

The problem was evolved for 20 time units, and the L1 norm of the density error between the computed and exact solution was calculated. The L1 norm is defined as

P kek1 ¼

v ^v  uv j v ¼1;N M L ju P v v ¼1;N ML

Fig. 20. Density contours at intermediate and final times for the triple point problem (quadratic scale).

ð39Þ

182

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

Fig. 21. Mesh velocity contours at intermediate and final times for the triple point problem.

Fig. 22. Density contours and surface mesh in the vicinity of the vortex at intermediate and final times for the triple point problem (quadratic scale).

where N is the number of points, MvL is the lumped mass of mesh ^ v and uv are the exact and computed solutions at mesh point v, and u point v, respectively (the above definition of the error norm is also used in the studies that follow). The study was performed using four meshes of increasing resolution. Table 1 summarizes the mesh data and the corresponding numerical errors in the calculated density. The initial (t ¼ 0) and final (t ¼ 20) mesh velocity distributions obtained on the coarsest mesh are shown in Fig. 1. The movement of the mesh over time is clearly visible. The log of the L1 error as a function of the log of the average edge length is shown in Fig. 2. The average convergence rate is 2:3, which is close to the expected second-order for smooth fields. The errors in the GCL result from the application of the discrete gradient operator in (37) to the nonlinear distribution given in (38). As will be shown subsequently, the errors are Oð105 —106 Þ times smaller than the errors measured in actual flow solutions on the same series of meshes. These results confirm that the GCL is satisfied to the truncation error of the numerical scheme. 5. Test problems 5.1. Shock tube The shock tube problem (see e.g. [57]) was calculated using the cylindrical mesh described previously, with c ¼ 1:4 and the following initial conditions:

ui ðxi Þ ¼ 0:0

ð40Þ

eðxi Þ ¼ 2:5

ð41Þ

qðxi Þ ¼ pðxi Þ ¼

1:0

z < 50:0

0:1 z P 50:0

1:0 z < 50:0 0:1 z P 50:0

ð42Þ

ð43Þ

where z is the distance along the length of the tube. Here, the length of the tube is 100 units. The prescribed initial data result in a rightmoving contact and shock with a left moving expansion fan/release

wave. All simulations used the Rusanov scheme with w ¼ 0:1; c1 ¼ 1:0, and c2 ¼ 0:0. Figs. 3 and 4 depict the surface mesh superimposed upon contours of density and mesh velocity, respectively. These results are depicted for two different time levels (t ¼ 10:25 and t ¼ 20) on the coarsest mesh. The initially uniform mesh gets stretched towards z ¼ 100 as the mesh tracks the faster moving contact surface and shock wave. This is illustrated by the higher mesh velocity located in the vicinity of these two features. To assess the spatial accuracy of the proposed numerical scheme, the numerical solution was extracted along a line through the center of the tube and compared with the exact solution for the one-dimensional (1D) Euler equations. The effect of mesh resolution on the predicted 1D profiles for density is presented in Fig. 5, which illustrates the solution along the tube at t ¼ 20:0 for each mesh. The exact solution at t ¼ 20:0 is also provided in the figure. As expected, the predictions improve with mesh resolution. No oscillations are observed near the discontinuities in the predictions and the shocks are accurately located. This illustrates the robustness of the proposed algorithm and its ability to accurately track shock fronts. The computed L1 norms of the error between the exact solution and the predictions for density at t ¼ 20:0 are summarized in Table 2 for each mesh. The log of the L1 error as a function of the average log of the average initial edge length is shown in Fig. 6. These results confirm the expected first-order decrease in error with mesh size that is displayed by limited second-order schemes. Here, the average slope is 0:8. To test the sensitivity of the results to the velocity smoothing algorithm, the L1 norm of the error in predicted density was examined on Mesh 1 for varying values of w . In this study, w was varied between 103 (effectively pure Eulerian) and 10 (effectively pure Lagrangian). The results are illustrated in Fig. 7. The average L1 error is 5:24  103 , with a maximum value of 6:38  103 at w ¼ 10 and a minimum value of 4:39  103 at w ¼ 0:1. This variation in the error of 20% results from the movement of mesh points into or out of regions that cannot be resolved by the spatial discretization. At large w , the shock and contact discontinuity are well resolved, but mesh points have been moved out of the expansion fan region. For small w , the resolution throughout the expansion

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

fan remains relatively fine but fewer mesh points are located near the shock and contact discontinuity. As such, larger amounts of numerical dissipation are observed near the contact and shock. The value of w at which L1 is minimized represents the optimal combination of these two competing errors.

5.2. Sedov blast wave – 3D variant Calculations of the Sedov problem [58] were performed with the Rusanov scheme and a source energy defined to produce a shock location of 1 cm at a time of 1 ls. The problem was modeled as an octant of a sphere with three meshes of increasing resolution. The following initial conditions were used with c ¼ 5=3; w ¼ 0:25, c1 ¼ 0:5, and c2 ¼ 0:0:

ui ðxi Þ ¼ 0:0

ð44Þ

qðxi Þ ¼ 1:0

ð45Þ

( eðxi Þ ¼ ( pðxi Þ ¼

1:0  104

xi – 0:0

es

xi ¼ 0:0

0:67  104

xi – 0:0

ps

xi ¼ 0:0

ð46Þ

ð47Þ

es and ps are initial values that depend on the size of the control volume at the center of the sphere. As the mesh size is decreased, these values must be increased to maintain a constant total energy source. Table 3 summarizes the source values along with the mesh data. Fig. 8 illustrates the density predicted using the ALE formulation along the x-axis for each mesh at t ¼ 1:0. A sharp discontinuity in density is observed at the location of the shock. The solution obtained using the coarsest mesh under-predicts the strength of the shock by a factor of approximately 1:8. However, the predicted solution rapidly approaches the exact solution as the mesh resolution is increased. As in the previous section, the spatial accuracy was assessed by extracting the numerical solution along the x-axis at t ¼ 1:0 and comparing with the exact solution for the 1D Euler equations. The L1 norm of the error between the predicted and exact solution for density was computed for the three meshes. The results for each mesh are included in Table 3. Fig. 9 shows the log of the L1 error as a function of the log of the average edge length in the initial mesh. Similar to the results obtained for the shock tube problem, first-order convergence was observed. The slope of the line in Fig. 9 was approximately equal to 1:0. The volume-integrated total energy over the domain as a function of time for each ALE calculation is shown in Fig. 11. Due to mesh faceting, the total mesh volume and hence total energy increases slightly with resolution. However, the total energy at each resolution is constant with time, indicating strict conservation at the discrete level. Note that the vertical range in the plot is of order 106 . Fig. 10 compares the ALE and Eulerian density profiles on Mesh 2. At this resolution, the peak density is increased by 30% in the ALE calculation relative to the Eulerian calculation. The pressure contours for the two calculations are compared in Figs. 12 and 13: the former shows the full pressure contours at intermediate and final times, and the latter shows a close-up near the shock at the final time. The smaller edge lengths near the shock and larger edge lengths near the origin are clearly visible in the ALE calculation. The movement of mesh with the shock results in a sharper profile and reduced numerical diffusion.

183

5.3. Sedov blast wave – 2D variant Pseudo-2D variants of the Sedov problem also were calculated to test symmetry preservation and to further demonstrate the impact of mesh motion on the solution. The 2D variants used a mesh of 48  48  2 structured hexahedra that were subdivided into approximately 5:3  104 quasi-uniform tetrahedra. Because the mesh has finite thickness, the source strength is different from the analytic source strength for the pure 2D case, which gives rise to the peak density greater than 4:0 observed in the Lagrangian and ALE calculations. Therefore a rigorous comparison to 2D analytic results is not possible and these studies focus on qualitative aspects only. Fig. 14 shows scatter plots of density as a function of radius for Lagrangian (wi ¼ ui ), ALE (w ¼ 0:05; c1 ¼ 1:0; c2 ¼ 0:0), and Eulerian (wi ¼ 0) variants at t ¼ 1:0. As expected, the peak density increases as the mesh velocity approaches the fluid velocity due to decreased mesh size at the shock front. All three variants demonstrate excellent symmetry preservation. Fig. 15 shows the density and computational mesh for each variant at t ¼ 0:5 and t ¼ 1:0. The ALE scheme behaves as expected with more resolution at the shock than the pure Eulerian calculation but less mesh motion in the central region than the pure Lagrangian calculation. 5.4. Taylor–Green vortex Although the intended application of the method presented herein is shock hydrodynamics, 3D simulations of the 2D steadystate Taylor–Green vortex were used to verify the accuracy of the method in the absence of discontinuities. Descriptions of the problem can be found in [7,45]. A computational domain of size 1  1  0:01 was used at three different resolutions and the problem was run to a final time of 0:75 time units using the Rusanov scheme and w ¼ 0:2; c1 ¼ 0:5, and c2 ¼ 0:5. The mesh data are summarized in Table 4. L1 error norms for the velocity relative to the initial conditions were calculated for each mesh along a line through y ¼ 0:25; z ¼ 0:005 and parallel to the x-axis. The L1 errors are included in Table 4, and the log of the L1 error as a function of the log of the average edge length in the initial mesh is shown in Fig. 16. The average convergence rate is 1:9, which is close to the expected value of 2 for smooth problems. Fig. 17 illustrates the instantaneous solution, i.e. computational mesh and fluid velocity, at four different times obtained using the coarse mesh. As time elapses, the mesh becomes distorted due to the vortical motion of the flow. This is especially apparent near the corners of the domain, one of which is shown in detail in Fig. 18 at t ¼ 0 and t ¼ 0:75. However the method exhibits robust convergence despite the mesh distortion. 5.5. Shock-triple point interaction This single material variant of the test problem proposed in [19] involves three regions initialized such that a shock propagates parallel to a contact discontinuity, which in turn generates baroclynic vorticity. In the single material variant the density and pressure vary between the three regions but each has the same ratio of specific heats. The geometry and initial conditions are shown in Fig. 19. Although this test problem does not have an exact solution, it represents a challenge for mesh motion algorithms due to multiple interacting shocks, vorticity development, and shear layers. It is therefore included here primarily as a test of robustness for simulations that must be run closer to the Eulerian limit. The 3D mesh consisted of approximately 2  106 tetrahedra and was one unit thick in the z-direction. The time-dependent density

184

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

0:8–1:0 were observed for shock dominated flows and 1:9 for smooth flows. These rates are consistent with expectations for the method. The method also was demonstrated to satisfy the discrete geometric conservation law to second-order accuracy, with truncation errors that are five to six orders of magnitude smaller than errors observed in flow solutions on the same meshes. The robustness of the method was further demonstrated on a shocktriple point interaction problem. Future work will focus on extensions to multiple materials, adaptive mesh refinement (AMR) [46], and high-order reconstruction techniques [47]. Appendix A. Edge-based finite element discretization

Fig. A.23. A two-dimensional geometric illustration of the edge-based FE formulation. The shaded area represents the dual-mesh control volume (i.e. lumped mass) for node v, and the vector Div w represents the area-normal of the dual mesh face separating nodes v and w within a given triangle (the full expression for Div w in (A.18) is the sum of the individual values for all triangles containing the edge). The loop represents the sum over edges surrounding node v, which is effectively the surface integral over ML .

and mesh velocity distributions are shown in Figs. 20 and 21 for intermediate and final times using the HLLC scheme with w ¼ 0:2; c1 ¼ 0:25, and c2 ¼ 0:5. The effect of c2 – 0 can be seen in Fig. 21 as localized regions of smaller wi . Fig. 22 shows a close-up of the density and the surface mesh in the vicinity of the vortex for intermediate and final times. Mesh quality is largely maintained despite the large amount of vorticity. 6. Conclusions A novel 3D FE ALE method for shock hydrodynamics on unstructured grids has been presented. The method combines an unsplit solution of the flux-conservative ALE hydrodynamic equations with an edge-based FE formulation and a mesh motion algorithm that bases the mesh velocity on the local fluid velocity. The mesh velocity is then smoothed to mitigate the possibility of mesh tangling. Fluxes are computed in the co-moving frame of the mesh from a Rusanov or HLLC approximate Riemann solver. This approach realizes the advantages of an unsplit ALE method while also allowing the mesh to follow the flow – which is an important feature of pure Lagrangian schemes. As a result, the proposed algorithm combines the advantages of Eulerian and Lagrangian approaches. Numerical diffusion at shocks is reduced compared to a pure Eulerian method and mesh tangling effects are reduced compared to a pure Lagrangian method. The proposed formulation offers several advantages over traditional Lagrange-plus-remap ALE approaches. Because the mesh velocity is determined based on the fluid velocity, it preserves anisotropic features of the initial mesh. It also eliminates the splitting errors introduced by the remap step. This reduces the overall truncation error of the scheme, and enables the use of spatial and temporal discretizations with orders of accuracy greater than two. The node-centered formulation, in contrast to staggered or cellcentered formulations, can compute accurate solutions on tetrahedral mesh cells. Lastly, the use of approximate Riemann solvers to compute fluxes eliminates the need for artificial viscosity treatments and allows direct application of Godunov concepts to Lagrangian and ALE formulations. Accuracy and convergence measurements were given for both shock-dominated and smooth flows. Global convergence rates of

In this appendix we present the details of the edge-based FE discretization. The derivation of the method largely follows that presented in [44], but additional discussion has been included to better illustrate the meaning of the method within the context of the numerical solution of hyperbolic conservation laws. The basic problem under consideration is the calculation of discrete spatial derivatives. Assuming a tessellation of the problem domain X into linear tetrahedral elements Xh , a field variable /ðxi Þ stored at the nodes v of the mesh, and linear basis functions N v ðxi Þ defined on each tetrahedral element, the gradient of /ðxi Þ at the nodes of the mesh is given by a Galerkin approximation as

XXZ Xh 2v w2Xh

Nv Nw



Xh

@/ @xi

w dXh ¼

XXZ Xh 2v w2Xh

Xh

Nv

@Nw w / dXh @xi

ðA:1Þ

where the outer summation is over tetrahedra Xh connected to point v, the inner summation is over points w forming Xh (which includes point v), and we have dropped the notational dependence of / and Nv on xi for brevity. The nodal values can be pulled out of the integrand which allows the above expression to be written in the more familiar matrix form,

Mv w



@/ @xi

w

¼ Giv w /w

ðA:2Þ

where M v w is the consistent mass matrix and Giv w (a vector quantity) represents a discrete gradient operator. These matrices are defined respectively as

XXZ

Mv w ¼

Nv Nw dXh

ðA:3Þ

Xh

Xh 2v w2Xh

and

Giv w ¼

XXZ Xh 2v w2Xh

Nv

Xh

@Nw dXh @xi

ðA:4Þ

In what follows we use the lumped mass-matrix M vL defined as

MvL ¼ M v w Iw

ðA:5Þ

where Iw is a unit column vector. The discrete equation therefore becomes

ML

@/ @xi

v

¼ Gvi w /w

ðA:6Þ

Note that the superscript v has been removed from ML to avoid the appearance of an implied summation on the left-hand-side of this expression. The goal of the derivation that follows is to express Gvi w in terms of summations over edges rather than summations over elements. This has two benefits: first, indirect addressing is reduced because an edge contains two nodes whereas a tetrahedra contains four; second, being a one-dimensional quantity an edge forms a natural setting for Godunov methods and approximate Riemann solvers.

185

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

An antisymmetric form of Gvi w also is desirable so that the discrete operator will be exactly conservative at the discrete level when applied to a hyperbolic conservation law (i.e. the method will be mimetic in the sense that the discrete operator mimics the properties of the continuous operator). The first step of the derivation is to separate the terms w ¼ v from the summation in (A.4):

XXZ

Nv Xh

Xh 2v w2Xh

XXZ

@Nw w @Nw w / dXh ¼ ð1  dv w Þ Nv / dXh @xi @xi Xh 2v w2Xh Xh XZ @Nv v þ Nv / dXh ðA:7Þ @xi X 2v Xh h

where dv w is the Kronecker delta. We next define the identity

Z

Z 1 @Nw w N / dXh þ Nv / dXh 2 Xh @xi @xi Xh Z Z 1 @Nw w 1 @Nv w w ¼ Nv / dXh  N / dXh 2 Xh 2 Xh @xi @xi Z 1 þ Nv Nw /w ni dCh 2 Ch  Z

1 @Nw @Nv w w ¼ Nv  N / dXh 2 Xh @xi @xi Z 1 þ Nv Nw /w ni dCh ðA:8Þ 2 Ch

w

Z

1 N / dXh ¼ 2 @xi Xh v @N

w

v @N

w

w

where integration by parts has been used on the second term in the first step, C ¼ @ X; ni is the outward normal on C, and Ch is a triangular face on C. In the case v ¼ w the first term vanishes so that

Z

Nv

Xh

@Nv v 1 / dXh ¼ 2 @xi

Z

Nv Nv /v ni dCh

ðA:9Þ

Ch

Substitution of these expressions into (A.7) gives X X Z v @N w @Nv w  w 1 N / dXh ¼ ð1  dv w Þ N  N / dXh 2 @xi @xi @xi Xh Xh 2v w2Xh Xh Z XX 1 þ ð1  dv w Þ N v N w /w ni dCh 2 Ch 2v w2Ch Ch Z 1X þ Nv Nv /v ni dCh ðA:10Þ 2 C 2v Ch

XXZ Xh 2v w2Xh

v @N

w

w

h

The last term can be simplified by noting that because a tetrahedra is a closed volume,

X @Nv @xi

v 2Xh

¼ 0 ! ð1  dv w Þ

X @Nw @Nv ¼ @x @xi i w2X

ðA:13Þ

h

We may therefore write

XXZ

@Nw v / dXh @xi Xh 2v w2Xh Xh Z XZ @N v v 1 ¼ Nv / dXh ¼ Nv Nv /v ni dCh 2 Ch @xi X 2v Xh

ð1  dv w Þ

Nv

where we have made use of (A.9) in the second step. Eq. (A.12) can now be written as

XXZ

@Nw w / dXh @xi Xh 2v w2Xh Xh X X Z v @Nw @Nv w  w 1 ¼ ð1  dv w Þ N  N ð/ þ /v ÞdXh 2 @x @x i i X h Xh 2v w2Xh XXZ 1 vw þ ð1  d Þ Nv Nw ð/w þ /v Þni dCh 2 Ch 2v w2Ch Ch XZ þ N v N v /v ni dCh ðA:15Þ Nv

Ch 2v

Ch

Focusing on the first term in (A.15), the summation represents an outer sum over all tetrahedra Xh connected to point v, and an inner sum over all points w – v 2 Xh . This summation is equivalent to an outer sum over pairs of points vw (each of which defines an edge vw that connects nodes v and w), and an inner sum over all Xh which contain the edge. The dual summation can therefore be rewritten as

XXZ

@Nw w / dXh @xi Xh 2v w2Xh Xh

 Z 1X X @Nw @Nv w ¼ Nv  N ð/w þ /v ÞdXh 2 v w2v X 2v w Xh @xi @xi h Z X X XZ 1 þ Nv Nw ð/w þ /v Þni dCh þ Nv Nv /v ni dCh 2 v w2v C 2v w Ch C 2v Ch Nv

h

h

The next step in the procedure is to add and subtract the quantity

XXZ

@Nw v / dXh @xi Xh 2v w2Xh Xh X X 1 Z v @Nw @Nv w  v ¼ ð1  dv w Þ N  N / dXh 2 Xh @xi @xi Xh 2v w2Xh XX1Z þ ð1  dv w Þ Nv Nw /v ni dCh 2 Ch X 2v w2X

ð1  dv w Þ

h

Xh 2v w2Xh

Xh

XXZ

ðA:11Þ

@Nw w w / / dXh @xi Xh 2v w2Xh Xh X vw w X vw w ¼ Di ð/ þ /v Þ þ Bi ð/ þ /v Þ þ Bvi /v v w2v

h

Nv rNw /w dXh

 1 @Nw @Nv w ¼ ð1  dv w Þ Nv  N ð/w þ /v ÞdXh 2 @xi @xi Xh 2v w2Xh Xh Z X X 1 þ ð1  dv w Þ Nv Nw ð/w þ /v Þni dCh 2 Ch 2v w2Ch Ch Z XXZ 1X @Nw v þ Nv Nv /v ni dCh  ð1  dv w Þ Nv / dXh 2 C 2v Ch @xi X 2v w2X Xh h

ðA:16Þ where a similar reinterpretation has been applied to the second term. Because the integrals do not depend on the nodal values, we can write the above as

Nv

from (A.10), using the right-hand-side in (A.11) for the additive part and the left-hand-side for the subtracted part. Combining terms gives

XXZ

XXZ

h

ðA:14Þ

h

h

ðA:12Þ

Nv

ðA:17Þ

v w2v

with

 Z

1 X @Nw @Nv w Nv  N dXh 2 X 2v w Xh @xi @xi h Z XZ 1 X Biv w ¼ Nv Nw ni dCh ; Bvi ¼ Nv Nv ni dCh 2 C 2v w Ch C 2v Ch

Div w ¼

h

ðA:18Þ ðA:19Þ

h

With implicit summations the expression for the discrete gradient can now be written as

ML

@/ @xi

v

¼ Div w ð/w þ /v Þ þ Bvi w ð/w þ /v Þ þ Bvi /v

ðA:20Þ

where the first term applies to all points in the mesh and the last two terms apply only to points on the boundary. Note that although the lumped mass matrix is used in the left-hand-side, the full

186

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187

consistent mass matrix also can be used with no impact on the v right-hand-side. Also note the identities Dw ¼ Div w and i

2

X vw X vw Di þ 2 Bi þ Bvi ¼ 0 v w2v

ðA:21Þ

v w2v

where the first follows from (A.18) and the second follows from the requirement that the gradient of a constant field (i.e. /w ¼ /v 8 w) be identically zero. Expressions similar to (A.20) can be written for other discrete operators. For example the discrete divergence of a vector field fi is given by

ML



@fi @xi

v

¼ Div w ðfiw þ fiv Þ þ Biv w ðfiw þ fiv Þ þ Bvi fiv

ðA:22Þ

Within the context of FV schemes, we can interpret the edgebased FE method as analogous to a point-centered FV method for which the solution is assumed to vary linearly within the control volume. The value of Dvi w on each cell represents the inward area normal of the dual mesh face separating nodes v and w. This is illustrated graphically in Fig. A.23, which shows a point v surrounded by a collection of triangles. For a hyperbolic conservation law (or system of conservation laws), the FV analogy applies directly and the flux between points v and w can be calculated from the solution to an approximate Riemann problem in the direction of Div w . Effectively this amounts to treating each edge as a local one-dimensional problem embedded in a three-dimensional space. Furthermore because Dvi w is antisymmetric, the method will be exactly conservative by construction. References [1] Hirt CW, Amsden AA, Cook JL. An arbitrary Lagrangian-Eulerian computing method for all flow speeds. J Comput Phys 1974;14:227–53. [2] Von Neumann J, Richtmyer RD. A Method for the numerical calculation of hydrodynamic shocks. J Appl Phys 1950;21:232. [3] Wilkins ML. Calculations of elastic-plastic flow. In: Methods in computational physics. New York: Academic Press; 1964. [4] Schulz WD. Two-dimensional Lagrangian hydrodynamic difference equations. In: Methods in computational physics. New York: Academic Press; 1964. [5] Caramana EJ, Burton DE, Shashkov MJ, Whalen PP. The construction of compatible hydrodynamics algorithms utilizing conservation of total energy. J Comput Phys 1998;146:227–62. [6] Maire PH, Abgrall R, Breil J, Ovadia J. A cell-centered Lagrangian scheme for two-dimensional compressible flow problems. SIAM J Sci Comput 2007;29:1781. [7] Dobrev VA, Ellis TE, Kolev Tz, Rieben RN. Curvilinear finite elements for Lagrangian hydrodynamics. Int J Numer Methods Fluids 2010;65:1295–310. [8] Scovazzi G. Lagrangian shock hydrodynamics on tetrahedral meshes: a stable and accurate variational multiscale approach. J Comput Phys 2012;231:8029–69. [9] Hansen GA, Douglass RW, Zardecki A. Mesh enhancement. London: Imperial College Press; 2005. [10] Knupp P, Margolin LG, Shashkov M. Reference Jacobian optimization-based rezone strategies for Arbitrary Lagrangian Eulerian methods. J Comput Phys 2002;176:93–128. [11] Knupp P. Algebraic mesh quality metrics for unstructured initial meshes. Finite Elem Anal Des 2003;39:217–41. [12] Berndt M, Moulton DJ, Hansen G. Efficient nonlinear solvers for Laplace– Beltrami smoothing of three-dimensional unstructured grids. Comput Math Appl 2008;55:2791–806. [13] Loubère R, Maire PH, Shashkov M, Breil J, Galera S. ReALE: a reconnectionbased arbitrary-Lagrangian–Eulerian method. J Comput Phys 2010;229:4724–61. [14] Owen JM, Toreja A, Miller D. The KULL ALE package. LLNL report no. UCRLPRES-404387, Livermore National Laboratory, Livermore, CA; 2008. [15] Rieben R, Jun BI, Managan R, Pudliner B. The ALE process in ARES. LLNL report no. UCRL-PRES-404459, Livermore National Laboratory, Livermore, CA; 2008. [16] Waltz J. ALE in the FLAG family of codes. LANL report no. LA-UR-08-04478, Los Alamos National Laboratory, Los Alamos, NM; 2008. [17] Peery JS, Carroll DE. Multi-material ALE methods in unstructured grids. Comput Methods Appl Mech Eng 2000;187:591–619. [18] Barlow AJ. A compatible finite element multi-material ALE hydrodynamics algorithm. Int J Numer Methods Fluids 2008;56:953. [19] Galera S, Maire PH, Breil J. A two-dimensional unstructured cell-centered multi-material ALE scheme using VOF interface reconstruction. J Comput Phys 2010;229:5755–87.

[20] Berndt M, Breil J, Galera S, Kucharik M, Maire P-H, Shashkov M. Two-step hybrid conservative remapping for multimaterial arbitrary Lagrangian– Eulerian methods. J Comput Phys 2011;230:6664–87. [21] Jia Z, Liu J, Zhang S. An effective integration of methods for second-order threedimensional multi-material ALE method on unstructured hexahedral meshes using MOF interface reconstruction. J Comput Phys 2013;236:513–62. [22] Lopez Ortega A, Scovazzi G. A geometrically-conservative, synchronized, fluxcorrected remap for arbitrary Lagrangian–Eulerian computations with nodal finite elements. J Comput Phys 2011;230:6709–41. [23] Formaggia L, Peraire J, Morgan K. Simulation of store separation using the finite element method. Appl Math Model 1988;12:175–81. [24] Murman SM, Aftosmis MJ, Berger MJ. Simulations of 6-DOF motion with a Cartesian method. In: 41st AIAA aerospace sciences meeting, AIAA-2003-1246, Reno, NV; 2003. [25] Nkonga B, Guillard H. Godunov type method on non-structured meshes for three-dimensional moving boundary problem. Comput Methods Appl Mech Eng 1994;113:183–204. [26] Takashi N. ALE finite element computations of fluid–structure interaction problems. Comput Methods Appl Mech Eng 1994;112:291–308. [27] Nguyen V-T. An arbitrary Lagrangian–Eulerian discontinuous Galerkin method for simulations of flows over variable geometries. J Fluids Struct 2010;26:312–29. [28] Rausch RD, Batina JT, Yang HTY. Three-dimensional time-marching aeroelastic analyses using an unstructured-grid Euler method. AIAA J 1993;31:1626–33. [29] Guillard H, Farhat C. On the significance of the geometric conservation law for flow computations on moving meshes. Comput Methods Appl Mech Eng 2000;190:1467–82. [30] Farhat C, Geuzine P, Grandmont C. The discrete geometric conservation law and the nonlinear stability of ALE schemes for the solution of flow problems on moving grids. J Comput Phys 2001;174:669–94. [31] Löhner R, Yang C. Improved ALE mesh velocities for moving bodies. Commun Numer Methods Eng 1996;12:599–608. [32] Murayama M, Nakahashi K, Matsushima K. Unstructured dynamic mesh for large movement and deformation. In: 40th AIAA aerospace sciences meeting, AIAA-2002-0122, Reno, NV; 2002. [33] Zeng D, Ethier CR. A semi-torsional spring analogy model for updating unstructured meshes in 3D moving domains. Finite Elem Anal Des 2005;41:1118–39. [34] Lin P, Baker TJ, Martinelli L, Jameson A. Two-dimensional implicit timedependent calculations on adaptive unstructured meshes with time evolving boundaries. Int J Numer Methods Fluids 2005;50:199–218. [35] Waltz J. Operator splitting and time accuracy in Lagrange plus remap solution methods. J Comput Phys 2013;253:57–67. [36] LeVeque RJ. Numerical methods for conservation laws. Berlin: Bierkhäuser Verlag; 1992. [37] Kolev Tz, Rieben RN. A tensor artificial viscosity using a finite element approach. J Comput Phys 2009;22:8336–66. [38] Lipnikov K, Shashkov M. A framework for developing a mimetic tensor artificial viscosity for Lagrangian hydrocodes on arbitrary polygonal meshes. J Comput Phys 2010;20:7911–41. [39] Barlow AJ, Roe PL. A cell centered Lagrangian Godunov scheme for shock hydrodynamics. Comput Fluids 2011;46:133–6. [40] Burton DE, Carney TC, Morgan NR, Sambasivan SK, Shashkov MJ. A cellcentered Lagrangian Godunov-like method for solid dynamics. Comput Fluids 2013;83:33–47. [41] Luo H, Baum JD, Löhner R. On the computation of multi-material flows using ALE formulation. J Comput Phys 2004;194:304–28. [42] Kershaw DS, Prasad MK, Shaw MJ, Milovich JL. 3D unstructured mesh ALE hydrodynamics with the upwind discontinuous finite element method. Comput Methods Appl Mech Eng 1998;158:81–116. [43] Shestakov AI, Milovich JL, Prasad MK. Combining cell- and point-centered methods in 3D unstructured-grid radiation-hydrodynamic codes. J Comput Phys 2001;170:81–111. [44] Waltz J. Microfluidics simulation using adaptive unstructured grids. Int J Numer Methods Fluids 2004;46:939–60. [45] Waltz J, Canfield TR, Morgan NR, Risinger LD, Wohlbier JG. Verification of a three-dimensional unstructured finite element method using analytic and manufactured solutions. Comput Fluids 2013;81:57–67. [46] Waltz J. Parallel adaptive refinement for unsteady flow calculations on 3D unstructured grids. Int J Numer Methods Fluids 2004;46:37–57. [47] Charest MRJ, Groth CPT, Gauthier PQ. High-order CENO finite-volume scheme for low-speed viscous flows on three-dimensional unstructured mesh. In: Proceedings of the 7th international conference on computational fluid dynamics, The Island of Hawaii, HI; 2012. [48] Waltz J. Performance of a three-dimensional unstructured mesh compressible flow solver on NVIDIA Fermi-class graphics processing unit hardware. Int J Numer Methods Fluids 2013;72:259–68. [49] Risinger LD, Waltz J, Canfield TR, Charest MRJ, Morgan NR, Wohlbier JG. Threading of indirect gather-scatter algorithms for 3D unstructured flow solvers. Comput Methods Appl Mech Eng 2013. submitted for publication. [50] Davis SF. Simplified second-order Godunov-type methods. SIAM J Sci Stat Comput 1988;3:445–73. [51] F Toro E, Sprice M, Speares W. Restoration of the contact surface in the HLL Riemann solver. Shock Waves 1994;4:25–34. [52] Batten P, Clarke N, Lambert C, Causon DM. On the choice of wavespeeds for the HLLC Riemann solver. SIAM J Sci Comput 1997;18:1553–70.

J. Waltz et al. / Computers & Fluids 92 (2014) 172–187 [53] Waltz J, Canfield TR, Morgan NR, Risinger LD, Wohlbier JG. Manufactured solutions for the three-dimensional Euler equations with relevance to Inertial Confinement Fusion. J Comput Phys 2013. submitted for publication. [54] Hirsch C. Numerical computation of internal and external flows, vols. 1 and 2. New York: Wiley; 1988. [55] Zhang H, Regiio M, Trepanier JY, Camarero R. Discrete form of the GCL for moving meshes and its implementation in CFD schemes. Comput Fluids 1993;22:9–23.

187

[56] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in C. 2nd ed. New York: Cambridge University Press; 1992. [57] Kamm JR, Brock JS, Brandon ST, Cotrell DL, Johnson B, Knupp P, et al. Enhanced verification test suite for physics simulation codes. LANL report no. LA-14379, Los Alamos National Laboratory, Los Alamos, NM; 2008. [58] Sedov LI. Similarity and dimensional methods in mechanics. New York: Academic Press; 1959.