A Lagrangian cell-centered discontinuous Galerkin ...

50 downloads 0 Views 7MB Size Report
total energy are approximated with linear Taylor expansions on a reference element.The ... 1 of 30. American Institute of Aeronautics and Astronautics ...
A Lagrangian cell-centered discontinuous Galerkin hydrodynamic method for 2D Cartesian and RZ axisymmetric coordinates Xiaodong Liu∗, Nathaniel R. Morgan†, and Donald E. Burton‡ Los Alamos National Laboratory, P.O. Box 1663, Los Alamos, NM, USA We present a new Lagrangian discontinuous Galerkin (DG) hydrodynamic method for solving the gas dynamic equations on unstructured hybrid meshes in 2D Cartesian and RZ axisymmetric coordinates. For 2D RZ axisymmetric coordinates, the method is based on a true volume scheme. The physical conservation laws for the mass, momentum and total energy are discretized using the DG approach. The specific volume, velocity and specific total energy are approximated with linear Taylor expansions on a reference element.The nodal velocity, and the corresponding forces, are calculated by solving a multidirectional approximate Riemann problem at the element nodes. An effective limiting strategy is presented that ensures monotonicity of the primitive variables. This new Lagrangian DG hydrodynamic method conserves mass, momentum, and total energy in 2D Cartesian and RZ axisymmetric coordinates. A suite of test problems are presented to demonstrate the robustness and expected second order accuracy of this new method.

Nomenclature J j w an w0 AN Ω Λs (x, y) (Z, R) (X, Y ) (ξ, η) u τ ρ σ

Jacobian matrix det(J) Volume of the current element Surface area normal of current element Volume of the initial element Surface area normal of initial element Reference element volume Surface area normal of reference element Natural 2D Cartesian coordinates Natural RZ axisymmetric coordinates Initial coordinates Reference coordinates Velocity Specific total energy Density Stress tensor

I.

p Pressure S Source term e Specific internal energy t Time M Mass matrix ψ Vector of basis functions Superscript ∗ Riemann solution at an element node k A vector of Taylor expansion coefficients n Temporal value Subscript p An element node c An element corner i An element surface segement g A Gauss quadrature point

Introduction

The Discontinuous Galerkin (DG) approach15, 16, 47 has been widely studied for high-order solutions of the gas dynamics equations in the Eulerian framework.10, 11, 25–27, 56, 59 As in classical finite element methods, accuracy for DG is obtained by means of high-order polynomial expansions within an element (also called a cell), rather than by using wide stencils as in the case of finite volume CCH methods. The DG(P) method ∗ Post-doc

research associate, X-Computational Physics Division. X-Computational Physics Division. Corresponding author: [email protected] ‡ Scientist, X-Computational Physics Division. † Scientist,

1 of 30 American Institute of Aeronautics and Astronautics

with P > 0 can be regarded as the natural extension of finite volume methods to higher order methods. The DG approach has many distinguishing features for Eulerian hydrodynamics.32–34 This motivates us to research and develop a DG method for the Lagrangian hydrodynamics.3, 5, 8, 17, 19, 20, 28–30, 35, 37, 41, 42, 58 Prior Lagrangian DG methods research has focused on developing methods for 2D Cartesian coordinates. Loubère et. al.31 proposed a Lagrangian DG method for 2D Cartesian triangular meshes and used a Bernstein basis. Jia21 developed a 2D Cartesian Lagrangian DG method with a spectral basis expansion in the Lagrangian coordinates (i.e., the original undeformed configuration). Vilar54, 55 developed a Lagrangian DG method for polytopal 2D Cartesian meshes. The unknowns were approximated with Taylor expansions in the Lagrangian coordinates. Li and Jia22 proposed a 2D Cartesian Lagrangian DG method with Taylor expansions in the natural coordinates, which differs from the previous Lagrangian DG methods that used Taylor expansions in the Lagrangian coordinates. Recent works by Liu et. al.24 and Morgan et. al.43 created a new 2D Cartesian Lagrangian DG method that used Taylor expansions in a reference element. These authors explored several ways to calculate the density variation over the element and proposed symmetry preserving limiting techniques for ensuring monotone solutions on problems with strong shocks that also preserve accuracy on smooth flows. Prior Lagrangian DG research efforts did not address RZ axisymmetric coordinates,8, 9, 12, 13, 23, 40, 51–53, 57 which motivates us to explore this topic. In this paper, we build on the method by Liu et. al.24 and propose the first Lagrangian DG hydrodynamic method for solving the gas dynamics equations in 2D Cartesian and RZ axisymmetric coordinates. The physical conservation laws for the mass, momentum and the total energy are discretized using the DG approach. The specific volume, velocity, and specific total energy fields are approximated with a Taylor expansion about the center of mass of a reference element. An explicit TVD second-order Runge-Kutta method is employed for time marching. A limiting strategy is presented that yields accurate results and ensures monotonicity in the specific volume, velocity, and pressure fields. The layout of the paper is structured as follows. The governing equations are described in Section II. The details on the new Lagrangian DG method are presented in Section III. The discrete evolution equations for momentum, total energy and specific volume (the inverse of density) are presented in Sections V, VI and IV respectively. The details of the nodal Riemann solver are discussed in Section VII. The results from a large suite of test problems are reported in Section IX. Concluding remarks are given in Section X.

II.

Governing equations

The analytic Lagrangian equations for specific volume, velocity, and specific total energy evolution in 2D coordinates are given by, dv =∇·u dt

(1)

d(u)T =∇·σ dt

(2)

ρ

ρ

dτ = ∇ · (σ · u) (3) dt where ρ is the density, v = 1/ρ is the specific volume, u is the velocity, σ is the stress, and τ is the specific total energy. The velocity vector and the right hand side of the governing evolution equations are, " # u u= Z , (4) uR ρ

∂ 1 ∂ (uZ ) + α (Rα uR ) ∂Z R ∂R " # T " # ∂ 1 ∂ α (σ ) + (R σ ) 0 α ZZ ZR R ∂R ∇ · σ = ∂Z − α σθθ 1 ∂ ∂ α (R σ ) + (σ ) RR RZ Rα ∂R ∂Z R ∇·u =

∇ · (σ · u) =

(5) T

 ∂ 1 ∂  α (σZZ uZ + σZR uR ) + α R (σRZ uZ + σRR uR ) ∂Z R ∂R 2 of 30

American Institute of Aeronautics and Astronautics

(6)

(7)

When α = 0, x corresponds to Z and y corresponds to R, which recovers the governing equations for 2D Cartesian coordinates24 (refer to Fig. 2 for more details). When α = 1, the governing equations for 2D R-Z coordinates can be recovered (refer to Fig. 1 for additional details). A convention is followed throughout the paper for denoting scalars, vectors and tensors. A slanted character is used to denote a scalar such as density ρ, specific total energy τ , and specific internal energy e. A slanted bold character is used to denote a vector such as velocity u, surface normal vector(n or s), Taylor basis function(ψ), and Riemann forces(F ∗i ). An upright character denotes a tensor such as stress σ, deformation gradient F, and Jacobian matrix J. At times it is useful to derive very general equations for an arbitrary unknown quantity that might be a scalar or a vector; in this case, a blackboard bold character U is used to denote the arbitrary unknown. The fluxes in a generalized evolution equation will be denoted as H , which can be either a vector or tensor. An arbitrary array of scalars or an array of vectors is denoted as Uk .

Figure 1: The reduction of the 3D problem to a 2D plane problem in Cartesian coordinates is shown. A 3D element is formed by stretching the corresponding 2D cell in x-y plane along axis z about δz = 1 for one unit length. Therefore the volume of the 3D element can be obtained by dw = δzdxdy = dxdy.

III.

Discretization

The solution domain is decomposed into non-overlapping discrete control volumes called cells or elements. The initial element volume is w0 , where the superscript 0 denotes the initial time t = t0 . Each element has a constant mass, dm dt = 0. The elements will deform with the flow and the volume at a later time will be w(t) (Fig. 3). Then the element volume w(t) can be calculated by dw = Rα dA with dA = dRdZ. The coordinate system for the initial configuration of the element is (X, Y )T and the mapping from the initial configuration to the current configuration is given by R = R (X, t) . The deformation gradient is F = configuration is

∂R ∂X .

An isoparametric mapping from a reference element Ω to the current

R = R (ξ, t) =

X

bp R p .

p

Rp are the nodal coordinates of the element, bp are the canonical shape functions at the nodes (for the element type), and the reference coordinates are (ξ, η). The Jacobian matrix is J = ∂R ∂ξ . Both the deformation

3 of 30 American Institute of Aeronautics and Astronautics

Figure 2: The reduction of the 3D rotational symmetry problem to a 2D meridian cut plane problem in R-Z coordinates is shown. The meridian cut plane is marked in blue lines. A 3D element is formed by rotating the corresponding 2D cell in the meridian cut plane around the axis Z for one unit radian. Therefore the volume of the 3D element can be obtained by dw = δθRdRdZ = RdRdZ.

gradient and the Jacobian matrix (F and J respectively) will vary in time.

Figure 3: The map from the initial configuration and the map from a reference element are graphically illustrated.

A.

Taylor expansion

Unknown fields (i.e., v, u, and τ ) can be represented with Taylor expansions on the reference element Ω about the center of mass. The Taylor expansion for an unknown U is ∂U ∂U U(ξ) = Ucm + (ξ − ξcm ) + (η − η ) (8) cm ∂ξ cm ∂η cm The subscript cm denotes the center of mass. The Taylor basis functions, for a linear expansion about the center of mass, are h i ψ = 1 ξ − ξcm η − ηcm (9)

4 of 30 American Institute of Aeronautics and Astronautics

The basis functions in the vector ψ are constant in time and are solely a function of the reference coordinates, which are temporally constant, as such, dψ =0 dt A Taylor expansion of an unknown can be succinctly written as,

(10)

U(ξ, t) = ψ (ξ) · Uk (t)

(11)

where 

 Ucm   Uk =  ∂U ∂ξ |cm  ∂U ∂η |cm

(12)

The Taylor expansion coefficients Uk are spatially constant on each element but they vary temporally. Uk is an array of scalars for the specific volume v and specific total energy τ fields. Uk is an array of vectors for the velocity field u. B.

The discontinuous Galerkin approach

The unknown fields (i.e., v, u, and τ ) in this work are polynomial expansions on a reference element so the governing equations must be transformed from the natural coordinates to the reference coordinates. The discontinuous Galerkin discretization will require converting volume integrals (that are over the current element configuration) to be either surface integrals or volume integrals over the reference element. The transformation of a volume integral to the reference element is, Z Z Z α ∇ · H(ξ)dw = ∇ · H(ξ)R dRdZ = ∇ · H(ξ)Rα jdΩ (13) w(t)



A(t)

where j = det(J) and H is either a vector or a tensor varying in space and time. We use the Gauss’s theorem for expressing the volume integrals associated with divergence operators on a vector h (Eq. 14) and a symmetric second-order tensor σ (Eq. 15) as follows, Z Z ∇ · h(ξ)dw = n · h(ξ)Rα da (14) w(t)

Z

∂A(t)

Z

n · σ(ξ)Rα da − α

∇ · σ(ξ)dw = w(t)

∂A(t)

Z σθθ (ξ)e R dA

(15)

A(t)

Here e R = [0 1], and n is a row vector denoting the normal vector. These relations are essential for deriving the Lagrangian discontinuous Galerkin method presented in this paper. The governing evolution equations have the generalized form, dU =∇·H dt where U is an unknown field (scalar or a vector) and H is a flux (vector or a tensor). The evolution equation is multiplied by the Taylor basis functions ψ and then integrated over the current element configuration.   Z dU ψT ρ − ∇ · H dw = 0 dt ρ

w(t)

which creates a system of q equations for q unknown Taylor expansion coefficients. Substituting a Taylor expansion into the above equation gives,

5 of 30 American Institute of Aeronautics and Astronautics

Z

dUk ρψ ψdw · = dt

Z

T

w(t)

ψ T (∇ · H)dw

(16)

w(t)

For the term on the left side of Eq. (16), the basis functions were factored out of the time derivative because dψ/dt = 0; likewise, the time derivative of the basis coefficients dUk /dt was factored out of the integral because they are solely a function of time for each element. The integral on the left side of Eq. (16) is the mass matrix, Z Z M= ρψ T ψdw = ρψ T ψRα jdΩ (17) Ω

w(t)

The evolution of the unknown basis coefficients is then, with integration by part for the term on the right side of Eq. (16),   Z Z ∇ · (ψ1 H) k dU   ∇ψ T · Hdw (18) = M· ∇ · (ψ2 H) dw − dt ∇ · (ψ3 H) w(t) w(t) Using Gauss’s Theorem Eq. (14) and (15), the resulting evolution equations for the unknown basis coefficients, i.e. v k , u k and τ k are, Z  Z  dv k ψ T (n · u ∗ )Rα da − ∇ξ ψ T · jJ−1 · uRα dΩ (19) M· = dt Ω

∂A(t)



du k = dt

Z

ψ T (n · σ∗ )Rα da − α

∂A(t)



dτ k = dt

Z

ψ T e R σθθ dA −



ψ T [n · (σ∗ · u )]Rα da −

 ∇ξ ψ T · jJ−1 · σRα dΩ

(20)



A(t)

Z

Z 

Z 

 ∇ξ ψ T · jJ−1 · (σ · u)Rα dΩ

(21)



∂A(t)

where the first term of the right hand side (RHS) in Eq. (19), (20), and (21) is a Riemann problem on the element surface. The Riemann problem is solved at the nodes of the element following Lagrangian finite volume methods (Section VII). The Riemann problem can be solved on the deformed linear element ∂w(t) by transforming back to the natural coordinates. The volume integrals over the reference element Ω are replaced with a summation over the Gauss quadrature points in Ω. The Lagrangian DG method reduces to a first-order finite volume method when using only the element average in the Taylor expansion. A second-order accurate solution is achieved when using the linear Taylor expansion. C.

Mass conservation

Mass conservation is enforced by guaranteeing that the mass matrix (Eq. 17) is constant in time. The total derivative of the mass matrix is  R R d  T T dM d α α = ρψ ψR jdΩ = ρψ ψR j dΩ dt dt dt Ω Ω  (22) R d d ψ T ψ dt (ρRα j) + ρRα j dt (ψ T ψ) dΩ = Ω

The terms inside the integral are dψ T ψ =0 dt d α (ρR j) = dt

0

Substituting these relations into Eq. (22) gives, dM =0 dt 6 of 30 American Institute of Aeronautics and Astronautics

(23)

D.

Temporal integration

A second-order explicit Runge-Kutta (RK) method is used to evolve the system of equations from time level k n to n + 1. The RK method has two stages and is applied to a system of equations of the form M · dU dt = Rk . The mass matrix M is constant in time (Section C), and Rk is the right hand side of Eq. (18). The time integration, for a vector of unknown Taylor expansion coefficients, is

Uk

n+1

s1 n Uk = Uk + ∆tM−1 · Rk n n s1 1 = 21 Uk + 21 Uk + 2 ∆tM−1 · Rk s1

(24)

The superscript s1 denotes the solution after the first stage. Meanwhile, the nodal position x p is updated using, s1

(Rp )

n+1

n

(Rp ) = (Rp ) + ∆t(u ∗p )n n s1 = 12 (Rp ) + 21 (Rp ) + 12 ∆t(u ∗p )s1

(25)

u ∗p (Section VII) denotes the Riemann velocity at a node.

IV.

Specific volume evolution

In this paper, we evolve a Taylor expansion of the specific volume field v over the reference element Ω. v (ξ, t) = ψ (ξ) · v k (t)

(26)

The density field is calculated at various locations in the element (e.g., Gauss points) by evaluating the Taylor expansion of the specific volume. Using the DG discretization in Eq. (19), the system of the equations for the evolution of the specific volume field is,  ρ¯w 0    0 

   0 0 P P vcm  α α (ρψ2 ψ2 )g Rg jg Ωg (ρψ2 ψ3 )g Rg jg Ωg  ∆  ∂v   · ∆t  ∂ξ |cm  = g∈Ω g∈Ω  P P ∂v (ρψ3 ψ2 )g Rgα jg Ωg (ρψ3 ψ3 )g Rgα jg Ωg ∂η |cm g∈Ω g∈Ω  P   ai Riα n i · u ∗p 0   i∈w(t) P    P (∇ξ ψ2 · jJ−1 · uRα )g Ωg     ψ2i ai Riα n i · u ∗p  =  − g∈Ω   P i∈w(t)  −1 α   P (∇ξ ψ3 · jJ · uR )g Ωg ψ a Rα n · u ∗ 3i i

i

i∈w(t)

i

p

(27)

g∈Ω

ρ¯ denotes the element average density and w is the current element volume. The Riemann velocity at the node is u ∗p (Section VII). The subscript g denotes the Gauss point locations in the reference element, the subscript p denotes the value at an element node, and the subscript i denotes the element surface segments connected to a node. ai Ri n i denotes the surface area normal vector of a segment i.

V.

Momentum evolution

The DG method in this work evolves a Taylor expansion of the velocity field u over the reference element Ω. u (ξ, t) = ψ (ξ) · u k (t) k

(28)

Here, u is an array of vectors. Using the DG discretization in Eq. (20), the system of the equations for the evolution of the velocity field is,

7 of 30 American Institute of Aeronautics and Astronautics

   0 0 P P u cm  α α (ρψ2 ψ2 )g Rg jg Ωg (ρψ2 ψ3 )g Rg jg Ωg  ∆  ∂u   · ∆t  ∂ξ |cm  = g∈Ω g∈Ω  P P ∂u (ρψ3 ψ2 )g Rgα jg Ωg (ρψ3 ψ3 )g Rgα jg Ωg ∂η |cm g∈Ω g∈Ω    P P   ai Riα n i · σ∗c (ψ1 jσθθ e R )Ωg 0  g∈Ω   i∈w(t)   P (∇ ψ · jJ−1 · σRα ) Ω  P   P    ξ 2 g g ψ2i ai Riα n i · σ∗c  (ψ2 jσθθ e R )Ωg  =   − g∈Ω  − α  g∈Ω  i∈w(t)   P −1 α P     P (∇ξ ψ3 · jJ · σR )g Ωg (ψ3 jσθθ e R )Ωg ψ3i ai Riα n i · σ∗c  ρ¯w 0    0 

g∈Ω

g∈Ω

i∈w(t)

(29)

ρ¯ denotes the element average density and w is the current element volume. The Riemann stress at an element node is σ∗c (Section VII). The subscript g denotes Gauss point locations in the reference element, the subscript p denotes the value at an element node, and the subscript i denotes the element surface segments connected to a node. ai Ri n i denotes the surface area normal vector of a segment i. The nomenclatures for the surface segments in the corner c1 connected to a node p is shown in Fig. 4. The contribution of node p to the second term of RHS is expressed as,   2   P P α ∗ ai Riα n i · σ∗c a R n · σ  i=1 i i i c    i∈w(t)     P 2  P α ∗  α ∗ ψ a R n · σ  2 i i c i i (30) = ψ2i ai Ri n i · σc      i∈w(t) i=1   P   2  P ψ3i ai Riα n i · σ∗c α ∗ ψ a R n · σ 3 i i c i i i∈w(t) p

. Here, a1 = 12 ap−p , a2 = 12 app+ , R1α ψ11 = R1α ψ21 =

α α 3Rp +Rp −

6 α α 3Rp +Rp − 6

ψ 1p + ψ 2p +

α α Rp +Rp −

6 α α Rp +Rp − 6

α α Rp +Rp −

α α 3Rp +Rp −

ψ1p− , R2α ψ12 = ψ2p− , R2α ψ22 =

α α 3Rp +Rp +

6 α α 3Rp +Rp + 6 α α 3Rp +Rp +

i=1

ψ 1p + ψ 2p +

α α Rp +Rp +

6 α α Rp +Rp + 6 α α Rp +Rp +

ψ 1p + , ψ 2p + ,

R1α ψ31 = ψ 3p + ψ3p− , R2α ψ32 = ψ 3p + ψ 3p + . 6 6 6 6 The subscript nomenclature of p− and p+ is from Maire et. al.35, 37 and is used to denote the two neighboring surface element nodes; in addition, the subscripts p − p and pp+ denote the corresponding faces (Fig. 4). Since the basis function component ψ1 at any node is always equal to 1.0, then R1α = R1α ψ11 = and

R2α

=

R2α ψ12

=

α α 2Rp +Rp +

3

α α 2Rp +Rp −

3

,

38, 49

; therefore, these parameters satisfy the geometry conservation law.

VI.

Total energy evolution

The DG method in this work evolves a Taylor expansion of the total energy field τ over the reference element Ω. τ (ξ, t) = ψ (ξ) · τ k (t)

(31)

Using the DG discretization in Eq. (21), the system of the equations for the evolution of the total energy is,

8 of 30 American Institute of Aeronautics and Astronautics

Figure 4: A patch of four elements is shown in the natural coordinate system. The element surface area is decomposed into smaller segments, where the surface area normal vector of a segment is ai Riα n i in the natural coordinates. The control surface for the Riemann problem is constructed from all segments connected to the node p and is denoted as i ∈ p. The control surface for the Riemann problem is highlighted in the figure with a dashed-line. The element center of mass is denoted as cm, and the centroid is z.

 ρ¯w 0    0 

   0 τcm  α (ρψ2 ψ3 )g Rg jg Ωg  ∆  ∂τ   · ∆t  ∂ξ |cm  = g∈Ω g∈Ω  P P ∂τ (ρψ3 ψ2 )g Rgα jg Ωg (ρψ3 ψ3 )g Rgα jg Ωg ∂η |cm g∈Ω g∈Ω  P   ai Riα n i · σ∗c · u ∗p 0   i∈w(t)   P (∇ ψ · jJ−1 · σ · uRα ) Ω   P ∗ α ∗   ξ 2 g g · u ψ a R n · σ 2 i i p c i i =  − g∈Ω   P i∈w(t)  −1 α  P (∇ ψ · jJ · σ · uR ) Ω α ∗ ∗ ξ 3 g g ψ3i ai Ri n i · σc · u p g∈Ω P

0 (ρψ2 ψ2 )g Rgα jg Ωg

P

(32)

i∈w(t)

The Riemann velocity at the node is u ∗p (Section VII). Please refer to Section V for other variables.

VII.

Riemann problem

In this work a multidirectional approximate Riemann problem is solved at the nodes of the element following the finite volume method in Burton et. al.5 The first multidirectional approximate Riemann problem for Lagrangian finite volume CCH was proposed by B. Després and C. Mazeran17 for gas dynamic problems. Maire and various authors35, 37 extended the work in17 and proposed a new multidirectional approximate Riemann problem that had improved mesh robustness properties. Burton et. al.5 proposed another robust multidirectional approximate Riemann solution that was suitable for materials with strength. Following the method by Burton et. al.,5 the Riemann force acting on a segment of the element surface is, ˆ c |ai Riα u ∗p − u c F ∗i = ai Riα n i · σ∗c = ai Riα n i · σc + µc |n i · e



(33)

The nomenclature is illustrated in Fig. 4. The acoustic impedance is µc = ρc where c is the sound speed of ˆ is a unit vector in the approximate direction of the shock, eˆ = (u avg the element. e − u c )/||(u avg − u c )|| p p avg where u p is the average of the surrounding corner velocities u c . Momentum conservation at the node requires that the summation of all forces around the node to be equal to zero. X ∗ Fi = 0 (34) i∈p

9 of 30 American Institute of Aeronautics and Astronautics

A single Riemann velocity is found at the node that guarantees momentum conservation. P ˆ c |ai Riα u c − ai Riα n i · σc ) (µc |n i · e i∈p ∗ P up = µc |n i · eˆ c |ai Riα

(35)

c∈p

Total energy is also conserved because,  X

F ∗i

 ∗

· up = 

i∈p

 X

F ∗i 

· u ∗p = 0

(36)

i∈p

Analysis is presented by Liu et. al.24 to show that the new Lagrangian DG method, combined with this Riemann problem, satisfies the second-law of thermodynamics. This approximate Riemann problem creates ˆ c | will apply the dissipation in the direction dissipation that is greater than or equal to zero. The term |n i · e ˆ c , which is an estimate of the shock direction.43 of the unit vector e The new Lagrangian DG hydrodynamic method proposed in this paper conserves mass, momentum, and total energy in 2D Cartesian and RZ axisymmetric coordinates. Furthermore, this Lagrangian DG method can preserve symmetry for 1D flow on an equal angle grid in 2D Cartesian coordinates,37, 41 but it doesn’t preserve symmetry in 2D RZ axisymmetric coordinates.6, 38

VIII.

Limiting

Limiting certain fields toward the element average value (i.e., piecewise constant over the element) is essential for monotone solutions on problems with shocks and discontinuities. A set of design goals will be followed for creating a limiting approach for the Lagrangian DG method. The first goal is to ensure that the pressure and velocity in an element corner is bounded by the neighboring element average values around the node, which ensures that the inputs to the Riemann problem are bounded. The second goal is to ensure that the evolved Taylor expansions (specific volume, momentum, and total energy) are within the permissible bounds defined by the neighboring elements around each node. The third goal is to guarantee second-order accuracy on smooth flows. A limiting approach will be described that satisfies these design goals. In this work, the velocity, specific volume, pressure and specific total energy fields are limited towards a constant field. The limiting is done after the pressure is calculated; in other words, the pressure is obtained using unlimited fields. The limiting on the specific total energy is only used in the temporal integration. All the limiting is done after each time integration step (Eq. 24), and the limited fields are used in the RK time integration. Some fields are only known at the discrete points inside the element (i.e., the pressure, density, and specific internal energy), while other fields are represented by a Taylor expansion (i.e., the specific volume, velocity, and specific total energy). A new limiting approach is presented that works for modal and nodal fields in 2D Cartesian and RZ axisymmatric coordinates. Likewise, limiters are presented for scalar and vector fields that preserve symmetry on equal-angle polar meshes in 2D Cartesian coordinates. A.

Scalar limiter, φ

The limiting of a scalar field in an element is achieved by, U (ξ)

lim

= U + φ U (ξ) − U



(37)

The superscript lim denotes the limited field, U is the element average value for the field, and U (ξ) is the unlimited field. The scalar limiter reduces the field to be equal to the element average with φ = 0; likewise, the field is unlimited when φ = 1. The scalar limiter is found using the Barth-Jesperson limiter,1    max  η(U −U)   if U ξ −U>0 min 1,  p  ξp )−U   U(min   φ= (38) −U) min 1, η(U if U ξ −U