A discontinuous Galerkin method based on a hierarchical orthogonal basis for Lagrangian hydrodynamics on curvilinear grids

Xiaodong Liu [email protected] Nathaniel R. Morgan [email protected] Evan J. Lieberman [email protected] Donald E. Burton [email protected] Remcom Inc; State College, PA, USA X-Computational Physics Division; Los Alamos National Laboratory; P.O. Box 1663, Los Alamos, NM, USA

Abstract

We present a new high-order accurate Lagrangian discontinuous Galerkin (DG) hydrodynamic method to simulate material dynamics (for e.g., gasses, fluids, and solids) with up to fourth-order accuracy on cubic meshes. The variables, such as specific volume, velocity, specific total energy, and deformation gradient fields within a cell, are represented with a polynomial constructed from a novel hierarchical orthogonal basis about the center of mass, which decouples the moments of the solution because the mass matrix is diagonal. The discontinuity in the polynomials at the cell boundary is addressed by solving a multi-directional Riemann problem at the vertices of the cell and a 1D Riemann problem at additional non-vertex quadrature points along the edges so that the surface integral is exact for the polynomial order. The uniqueness lies in that the vertices of the curvilinear grid work as the quadrature points for the surface integral of DG methods. To ensure robust mesh motion, the pressure for the Riemann problem accounts for the difference between the density variation over the cell and a density field from subcell mesh stabilization (SMS). The accuracy and robustness of the new high-order accurate Lagrangian DG hydrodynamic method is demonstrated by simulating a diverse suite of challenging test problems covering gas and solid dynamic problems on curvilinear grids.

keywords:

Lagrangian gas and solid dynamics , discontinuous Galerkin , fourth-order accurate , cubic curvilinear cells , a hierarchical orthogonal basis

1 Introduction

Lagrangian hydrodynamic methods solve the governing physics equations for material dynamics on a mesh that moves with the flow and historically are based on linear cells (also called elements or zones) that have edges defined by a polynomial of degree 1. One of the first Lagrangian hydrodynamic methods was the finite volume (FV) staggered-grid hydrodynamic (SGH) method by von Neumann and Richtmyer [1]. A range of FV SGH methods [2, 3, 4, 5, 6] have been proposed since the work of von Neumann and Richtmyer. The FV Lagrangian cell-centered hydrodynamic (CCH) approach was proposed by Godunov [7, 8] and has garnered significant interest following the work of Després and Mazeran [9] who developed a nodal Riemann solver that greatly improved the robustness of CCH methods over using a Riemann solver at the cell faces [10]. A range of robust nodal Riemann solvers have been proposed for use in Lagrangian CCH methods [11, 12, 13, 14]. Discontinuous Galerkin (DG) [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45] hydrodynamic methods have been developed for Lagrangian motion with linear meshes in Cartesian coordinates [46, 47, 48, 49, 50, 51, 52] and 2D R-Z coordinates [53, 54]. The DG method evolves a polynomial of order Pn forward in time, which is denoted as DG(Pn), whereas FV CCH methods fit the neighboring cells to build a polynomial. Lagrangian hydrodynamic methods with linear cells deliver up to second-order accuracy on multi-dimensional problems.

Recently, some authors have proposed third-order accurate 2D Lagrangian hydrodynamic methods for quadratic cells, which have edges defined by a polynomial of degree 2. Cheng and Shu [55] developed a weighted essentially non-oscillatory (WENO) FV CCH method for quadratic cells. Vilar et al. [48] proposed a third-order accurate Lagrangian DG(P2) method for solving the 2D gas dynamics on polygonal quadratic cells using a hierarchical slope limiter for shock capturing [56, 57]. Spurious mesh motion can still arise on some strong shock problems with quadratic cells. Cheng and Shu addressed the spurious mesh motion by controlling the curvature of the cell surface. In a context of the second-order FV and DG method with a Barth slope limiter [58], Morgan et al. [59] proposed a velocity filter to dissipate spurious mesh motion by modifying the velocity reconstructions. Liu et al. [60, 61] developed a novel subcell mesh stabilization (SMS) scheme for ensuring robust mesh motion with a third-order Lagrangian DG method on quadratic cells. [60] limited the second-order DG(P1) solution truncated from the high-order polynomial near the shock and [61] proposed symmetry-preserving Hermite WENO for shock capturing. Later, Liu et al. [62] also explored the Lagrangian DG(P1) method on quadratic triangles [63] with SMS for improved accuracy on large deformation flows. Lieberman et al.[64] extended the third-order accurate Lagrangian DG method with SMS [60] to support hyper and hypo-elastic plastic models for simulating solid dynamic problems using quadratic cells and to simulate multi-phase high-explosives detonation physics [65].

The objective in this paper is to create a new high-order accurate Lagrangian DG hydrodynamic method to simulate material dynamics (for gases, fluids, and solids) up to fourth-order accuracy with cubic cells, which have edges defined by a polynomial of degree 3. A cell with edges of more degrees of freedom is essential to accurately simulate large deformation problems because it is less likely to tangle compared with quadratic- or linear-order cells. To-date, there are no fourth-order accurate compact stencil Lagrangian hydrodynamic methods with cubic cells. High-order Lagrangian methods that use a compact numerical stencil are desirable for achieving better scaling and performance on modern computer architectures. The SMS method was first proposed by Liu et al. [60] to ensure accurate and robust mesh motion for quadratic cells, and in this work, we extend SMS to cubic cells in the context of this new high-order accurate Lagrangian DG hydrodynamic method. The polynomials in this work are constructed using a novel hierarchical orthogonal basis about the center of mass, which differs from prior Lagrangian DG research efforts that used Taylor basis functions. Using a Taylor basis will create a mass matrix that must be inverted at the start of a calculation, and couples together the equations for the high-order terms, which is problematic when limiting the second-order (DG(P1)) solution truncated from the high-order polynomial for preserving monotonicity near a shock. We show that the hierarchical orthogonal basis combined with SMS is essential for robust mesh motion with cubic cells on strong shock problems using the aforementioned limiting strategy. The orthogonal basis is also constant in time for Lagrangian motion, so they only need to be calculated at the start of a simulation, which is key for a computationally efficient scheme. An explicit multi-step Runge-Kutta (RK) method is employed for time-marching. The new high-order accurate Lagrangian DG hydrodynamic method conserves mass, momentum, and total energy, produces robust solutions on strong shock problems and yields high-order accurate solutions on smooth flows using curvilinear grids.

The layout of the paper is as follows. The governing equations are described in Section 2. The constitutive models are described in Section 2.1. Section 3 discusses the new hierarchical orthogonal basis, the Lagrangian DG method, the evaluation of the right hand side ( $rhs$ ) of the discretized equations, the SMS scheme, and the limiting strategy for shock capturing. The results from a large suite of test problems are reported in Section 4 covering solid and gas dynamics. Concluding remarks are given in Section 5.

2 Governing equations

In Lagrangian hydrodynamics, the differential equations governing the evolution of the specific volume ( $v$ ), velocity (u), specific total energy ( $\tau$ ), and deformation gradient F are given by,

\rho\frac{{dv}}{dt}=\nabla\cdot{\textbf{{u}}},

(1)

\rho\frac{d{\textbf{{u}}}}{dt}=\nabla\cdot{{\upsigma}},

(2)

\rho\frac{d\tau}{dt}=\nabla\cdot({\upsigma}{\textbf{{u}}}),

(3)

\frac{d\text{F}}{dt}=\nabla_{\textbf{{X}}}{\textbf{{u}}},

(4)

here the subscript X in Eq. 4 denotes the gradient with respect to the initial coordinates X instead of the current physical coordinates x, and

\nabla_{\textbf{{X}}}{\textbf{{u}}}=\begin{bmatrix}\frac{\partial u}{\partial X}&\frac{\partial v}{\partial X}\\ \frac{\partial u}{\partial Y}&\frac{\partial v}{\partial Y}\end{bmatrix}.

(5)

The Cauchy stress tensor, pressure, specific internal energy, and specific kinetic energy are denoted as ${\upsigma}$ , $p$ , $e$ , and $k$ respectively. The specific total energy is $\tau=e+k$ . For gas dynamics, the stress is given by ${\upsigma}=-p\text{I}$ , where the pressure is calculated using an equation of state (EOS) for the material, $p=\text{EOS}(\rho,e)$ ; for solid dynamics, the stress is given by ${\upsigma}=-p\text{I}+{\upsigma}^{\prime}$ , where either a hypo- or hyper-elastic-plastic model is used to calculate the deviatoric stress ${\upsigma}^{\prime}$ . The derivations in this work are valid for both gas and solid dynamics, where gas dynamics can be regarded as a special case. The time derivatives are total derivatives that move with the flow. The rate of change of the position is,

\frac{d{\textbf{{x}}}}{dt}={{\textbf{{u}}}}.

(6)

The equations presented in this paper will adhere to the following conventions. A slanted character is used to denote a scalar such as density $\rho$ and specific total energy $\tau$ . A slanted bold character is used to denote a vector such as velocity u and Riemann forces ${\textbf{{F}}}_{i}^{*}$ . An upright character denotes a tensor such as stress $\upsigma$ , deformation gradient F, and Jacobian matrix J. A blackboard bold character $\mathbb{U}$ is used to denote an arbitrary unknown. The fluxes in a generalized evolution equation will be denoted as $\mathbb{H}$ , which can be either a vector or tensor. Since the following sections involve both the inner product of two vectors and the matrix product, in order to describe these operations in the paper more easily and clearly, unless stated, the vectors (e.g., u and n) are column vectors. Therefore, ${\textbf{{n}}}\cdot{\textbf{{u}}}$ denotes the inner product and ${\upsigma}{\textbf{{u}}}$ (or ${\upsigma}{\textbf{{n}}}$ ) means the matrix product. In this work we follow the convention of matrix manipulation if the multiplication or product involves one matrix, and the dot product is used to represent the inner product of two vectors.

2.1 Constitutive models

We will begin with describing the EOS’s for calculating the pressure with gases and solids followed by a discussion on two approaches to calculate the deviatoric stress with a solid material. The gases are modeled using a gamma-law EOS model given by

p=\rho e(\gamma-1),

(7)

where $\gamma$ is the adiabatic index (ratio of specific heats). The pressure in solids is modeled using a Grüneisen EOS model that is given by

$\displaystyle P_{H}(\eta)$	$\displaystyle=\frac{\rho^{0}c_{0}^{2}\eta}{(1-s\eta)^{2}},$	(8)
$\displaystyle E_{H}(\eta)$	$\displaystyle=\frac{\eta P_{H}}{2\rho^{0}},$
$\displaystyle p(\eta,e)$	$\displaystyle=P_{H}+\Gamma\rho(e-E_{H}),$

Here $\eta=1-\rho^{0}/\rho$ , $\rho^{0}$ and $\rho$ are the initial and current density; $e$ is the specific internal energy; $P_{H}$ and $E_{H}$ are the Hugoniot pressure and energy; $\Gamma$ is the Grüneisen coefficient; and $c_{0}$ and $s$ are parameters that relate the shock speed and the particle velocity according to $u_{s}=c_{0}+su_{p}$ .

We consider strength formulations based on infinitesimal strain and also finite strain. We will start by describing the infinitesimal strain formulation followed by the finite strain formulation. The infinitesimal strain is $\upvarepsilon=sym({\text{F}}-{\text{I}})$ and the deviatoric infinitesimal strain is $\upvarepsilon^{\prime}=\upvarepsilon-\frac{1}{3}\text{tr}(\upvarepsilon)$ . The deviatoric stress is

{\upsigma}^{\prime}=2G{\upvarepsilon}^{\prime}_{e},

(9)

where $G$ is the material shear modulus and ${\upvarepsilon}^{\prime}_{e}$ is the deviatoric elastic strain considering the effect of the yield criterion, that is based on a $J_{2}$ radial return method. Assuming no new plastic deformation has occurred, a trial deviatoric stress $\hat{\upsigma}^{\prime}$ can be calculated by $\hat{\upsigma}^{\prime}=2G({\upvarepsilon}^{\prime}-{\upvarepsilon}^{\prime}_{p})$ . The equivalent form of the Cauchy stress tensor is defined as,

{\sigma}_{eq}=\sqrt{3J_{2}}=\sqrt{\frac{3}{2}{\upsigma}^{\prime}:{\upsigma}^{\prime}}

(10)

The material yield stress is $Y$ . If ${\sigma}_{eq}<Y$ , the assumption that no new plastic occurs holds. This means the deviatoric plastic strain ${\upvarepsilon}^{\prime}_{p}$ is kept the same as the one of the previous time step. Therefore,

\begin{split}&{\upsigma}^{\prime}=\hat{\upsigma}^{\prime},\\ &{\upvarepsilon}^{\prime}_{e}={\upvarepsilon}^{\prime}-{\upvarepsilon}^{\prime}_{p}.\end{split}

(11)

If ${\sigma}_{eq}>Y$ , then the $J_{2}$ radial return method is imposed to capture the plastic deformation.

\begin{split}&{\upsigma}^{\prime}=\hat{\upsigma}^{\prime}\frac{Y}{{\sigma}_{eq}},\\ &{\upvarepsilon}^{\prime}_{p}={\upvarepsilon}^{\prime}-\frac{{\upsigma}^{\prime}}{2G}.\end{split}

(12)

We will now shift the discussion to focus on the finite strain approach.

The finite strain model uses the Green-Lagrange strain ${\text{E}}=\frac{1}{2}({\text{C}}-{\text{I}})$ , where the right Cauchy-Green strain ${\text{C}}={\text{F}}^{T}{\text{F}}$ . Accordingly, the above strain definitions can also be extended for the elastic (i.e., ${\text{C}}_{e}$ and ${\text{E}}_{e}$ ) or deviatoric strain (i.e., ${\text{C}}^{\prime}$ and ${\text{E}}^{\prime}$ ) if we use the elastic or deviatoric part (i.e., ${\text{F}}_{e}$ or ${\text{F}}^{\prime}$ ) of the deformation gradient. The deviatoric elastic Green-Lagrange strain ${\text{E}}^{\prime}_{e}$ is calculated using the deviatoric elastic deformation gradient ${\text{F}}^{\prime}_{e}=(\text{det}{\text{F}}_{e})^{-\frac{1}{3}}{\text{F}}_{e}$ . The elastic deformation gradient is calculated from the multiplicative decomposition of the deformation gradient ${\text{F}}_{e}={\text{F}}{\text{F}}_{p}^{-1}$ . The second Piola-Kirchoff (PK2) stress $\tilde{\text{S}}$ is calculated by

\tilde{\text{S}}=2G{\text{E}}^{\prime}_{e}.

(13)

The deviatoric PK2 stress in the intermediate configuration is found by

\tilde{\text{S}}^{\prime}=\tilde{\text{S}}-\frac{1}{3}\left(\tilde{\text{S}}:({\text{F}}_{e}^{T}{\text{F}}_{e})\right)\left({\text{F}}_{e}^{T}{\text{F}}_{e}\right)^{-1}.

(14)

The deviatoric Cauchy stress in the current configuration is found from the deviatoric PK2 stress using,

{\upsigma}^{\prime}=\frac{1}{\text{det}{\text{F}}}{\text{F}}_{e}\cdot\tilde{\text{S}}^{\prime}\cdot{\text{F}}_{e}^{T}.

(15)

Further details on the finite strain elastic-plastic model with the DG(P2) method can be found in Lieberman et al. [64].

3 Discretization

The computational domain is decomposed into non-overlapping cells (also called elements or zones) that connect together to create a mesh. The initial cell volume is $w^{0}$ , where the superscript $0$ denotes the initial time $t=t^{0}$ . Each cell has a constant mass, $\frac{dm}{dt}=0$ . The cells will deform with the flow and the volume at a later time will be $w(t)$ (Fig. 1). The deformation gradient F and Jacobain matrix J is,

{\text{F}}=\frac{\partial{\textbf{{x}}}}{\partial{\textbf{{X}}}}=\begin{bmatrix}\frac{\partial x}{\partial X}&\frac{\partial x}{\partial Y}\\ \frac{\partial y}{\partial X}&\frac{\partial y}{\partial Y}\end{bmatrix},\quad{\text{J}}=\frac{\partial{\textbf{{x}}}}{\partial{\boldsymbol{\xi}}}=\begin{bmatrix}\frac{\partial x}{\partial\xi}&\frac{\partial x}{\partial\eta}\\ \frac{\partial y}{\partial\xi}&\frac{\partial y}{\partial\eta}\end{bmatrix}.

(16)

Here ${\textbf{{X}}},\textit{i.e.},(X,Y)$ and ${\boldsymbol{\xi}},\textit{i.e.},(\xi,\eta)$ denote the coordinate system for the initial and reference configuration. Both F and J vary in time.

Refer to caption — Figure 1: The maps are illustrated for a cubic quadrilateral cell. The cell map from the reference coordinates to the physical coordinates is created using the cubic serendipity Lagrange interpolation basis functions. A fourth-order accurate Lagrangian hydrodynamic method requires a cubic cell.

3.1 Basis functions for DG methods

In DG methods, numerical polynomial solutions ${\mathbb{U}}_{h}$ in each cell can be represented using a range of basis functions including the Taylor basis [15, 16], standard Lagrange finite element basis [66], or the Legendre-Dubiner basis [67, 27]. Taylor basis functions have been used in many Lagrangian DG hydrodynamic methods; however, in this work we create a scheme based on new hierarchical orthogonal basis functions[34]. The polynomial in a cell for $v$ , u, and $\tau$ can be written as

{\mathbb{U}}_{h}({\boldsymbol{\xi}},t)=\sum\limits_{k}{\phi}_{k}\left({\boldsymbol{\xi}}\right)\cdot{\mathbb{U}}_{k}\left(t\right)\quad\quad k=1,...,N(P)

(17)

where $\mathbb{U}_{h}$ denotes the respective field, ${\phi}_{k}$ denotes the hierarchical orthogonal basis functions, ${\mathbb{U}}_{k}$ are the unknown coefficients and $N(P)$ denotes the total number of terms in the numerical solution expansions for the polynomial with degree $P$ ( $N(P)$ is equal to 10 for DG(P3) in 2D). The basis functions ${\phi}_{k}$ are calculated by applying the Gram-Schmidt orthogonalization to the Taylor-basis functions ${\psi}_{k}$ described in the following equation.

\begin{array}[]{l}\phi_{1}=1,\\ \phi_{2}=\psi_{2}-\frac{\langle\psi_{2},\phi_{1}\rangle}{\langle\phi_{1},\phi_{1}\rangle}\phi_{1},\\ \phi_{3}=\psi_{3}-\frac{\langle\psi_{3},\phi_{2}\rangle}{\langle\phi_{2},\phi_{2}\rangle}\phi_{2}-\frac{\langle\psi_{3},\phi_{1}\rangle}{\langle\phi_{1},\phi_{1}\rangle}\phi_{1},\\ \phi_{4}=\psi_{4}-\frac{\langle\psi_{4},\phi_{3}\rangle}{\langle\phi_{3},\phi_{3}\rangle}\phi_{3}-\frac{\langle\psi_{4},\phi_{2}\rangle}{\langle\phi_{2},\phi_{2}\rangle}\phi_{2}-\frac{\langle\psi_{4},\phi_{1}\rangle}{\langle\phi_{1},\phi_{1}\rangle}\phi_{1},\\ \vdots\\ \phi_{n}=\psi_{n}-\frac{\langle\psi_{n},\phi_{n-1}\rangle}{\langle\phi_{n-1},\phi_{n-1}\rangle}\phi_{n-1}-\dots-\frac{\langle\psi_{n},\phi_{0}\rangle}{\langle\phi_{0},\phi_{0}\rangle}\phi_{0},\end{array}

(18)

where the inner product of two functions $g(\boldsymbol{\xi})$ and $h(\boldsymbol{\xi})$ is defined as,

\langle g,h\rangle=\int\limits_{\Omega_{h}}\rho ghj^{0}d\Omega.

(19)

The Taylor basis functions up to P3 for $v_{h}$ , $\textbf{{u}}_{h}$ , and $\tau_{h}$ are about the center of mass ( $cm$ ) and defined as

\begin{array}[]{l}\psi_{1}=1,\\ \psi_{2}=\xi-\xi_{cm},\\ \psi_{3}=\eta-\eta_{cm},\\ \psi_{4}=(\xi-\xi_{cm})^{2},\\ \psi_{5}=(\eta-\eta_{cm})^{2},\\ \psi_{6}=(\xi-\xi_{cm})(\eta-\eta_{cm}),\\ \psi_{7}=(\xi-\eta_{cm})^{3},\\ \psi_{8}=(\eta-\eta_{cm})^{3},\\ \psi_{9}=(\xi-\xi_{cm})^{2}(\eta-\eta_{cm}),\\ \psi_{10}=(\xi-\xi_{cm})(\eta-\eta_{cm})^{2}.\end{array}

(20)

Here, the subscript $cm$ denotes the center of mass, that is defined by ${\boldsymbol{\xi}_{cm}}=\frac{1}{m}{\int\limits_{w_{h}}\rho{\boldsymbol{\xi}}dw}$ , $m$ denotes the mass of the cell, namely $m={\int\limits_{w_{h}}\rho dw}$ . In A, we show that the coefficients(i.e., $\frac{\langle\psi_{2},\phi_{1}\rangle}{\langle\phi_{1},\phi_{1}\rangle}$ , $\frac{\langle\psi_{3},\phi_{2}\rangle}{\langle\phi_{2},\phi_{2}\rangle}$ , $\frac{\langle\psi_{3},\phi_{1}\rangle}{\langle\phi_{1},\phi_{1}\rangle}$ , etc.) are temporally constant, meaning this set of coefficients are only calculated once at the beginning of the calculation. Temporally constant basis functions, $\frac{d{\phi}_{k}}{dt}=0$ , are an essential ingredient in the presented Lagrangian hydrodynamic method.

The discussion will now address the polynomial for the deformation gradient. The basis functions for $\text{F}_{h}$ are about the center of the reference cell $\boldsymbol{\xi}_{c}=\frac{1}{v}{\int\limits_{w^{0}}{\boldsymbol{\xi}}dw}$ because it is a per volume quantity. Likewise, the Taylor basis functions are,

\begin{array}[]{l}\hat{\psi}_{1}=1,\\ \hat{\psi}_{2}=\xi-\xi_{c},\\ \hat{\psi}_{3}=\eta-\eta_{c},\\ \hat{\psi}_{4}=(\xi-\xi_{c})^{2},\\ \hat{\psi}_{5}=(\eta-\eta_{c})^{2},\\ \hat{\psi}_{6}=(\xi-\xi_{c})(\eta-\eta_{c}),\\ \hat{\psi}_{7}=(\xi-\eta_{c})^{3},\\ \hat{\psi}_{8}=(\eta-\eta_{c})^{3},\\ \hat{\psi}_{9}=(\xi-\xi_{c})^{2}(\eta-\eta_{c}),\\ \hat{\psi}_{10}=(\xi-\xi_{c})(\eta-\eta_{c})^{2}.\end{array}

(21)

The orthogonal basis function for $\text{F}_{h}$ are denoted as $\hat{\boldsymbol{\phi}}$ and are calculated using Eq. 19 without the density. The basis functions are constant in the calculation. If the initial density is constant in the cell, then the Taylor basis functions are identical for the per mass quantities and per volume quantities, namely ${\boldsymbol{\psi}}=\hat{\boldsymbol{\psi}}$ .

3.2 Discontinuous Galerkin method

In this section, we present a DG discretization of the governing evolution equations. The DG approach creates a system of equations for evolving the polynomial coefficients for $v$ , $\textbf{{u}}_{h}$ , $\tau_{h}$ , and $\text{F}_{h}$ . The evolution equations for specific volume, velocity, and specific total energy have the general form of

\rho\frac{d\mathbb{U}_{h}}{dt}=\nabla\cdot\mathbb{H}

where $\mathbb{U}_{h}$ is the respective field and $\mathbb{H}$ is the associated flux. This equation is multiplied by a test function and integrated over the cell in the physical coordinates,

\int\limits_{w_{h}}\phi_{m}\left(\rho\frac{d\mathbb{U}_{h}}{dt}-\nabla\cdot\mathbb{H}\right)dw=0,\quad m=1,...,N(P)

Substituting the polynomial expansion for the field gives,

\sum\limits_{n}\int\limits_{w_{h}}\rho{\phi_{m}\phi_{n}}dw\frac{d\mathbb{U}_{n}}{dt}=\int\limits_{\Omega_{h}}{\phi_{m}}(\nabla\cdot\mathbb{H})dw,\quad\quad m=1,...,N(P)

The first term is the mass matrix for the cell, ${M}_{mn}=\sum\limits_{k}\int\limits_{\Omega_{h}}\rho{\phi_{m}\phi_{n}}dw$ , and it is constant in time. Since the basis functions are orthogonal, the mass matrix is diagonal. Using integration by parts for the $rhs$ , the DG equations become,

{\text{M}}_{mm}\frac{d\mathbb{U}_{m}}{dt}=\int\limits_{w(t)}\nabla\cdot({\phi}_{m}\mathbb{H})dw-\int\limits_{w(t)}\mathbb{H}\cdot(\nabla{\phi}_{m})dw,\qquad m=1,...,N(P),

where $\nabla{\phi}_{m}=[\frac{\partial{\phi}_{m}}{\partial x}\quad\frac{\partial{\phi}_{m}}{\partial y}]^{T}$ . The DG equations require solving a Riemann problem on the surface of the reference cell to resolve the discontinuity inside the polynomials. The flux from the Riemann solver is denoted by a superscript $*$ . Using the Riemann solution and transforming the volume integral terms to be in the reference coordinate system, the resulting evolution equations for the unknown basis coefficients for specific volume $v_{n}$ , velocity ${\textbf{{u}}}_{n}$ , and specific total energy $\tau_{n}$ are,

{\text{M}}_{mm}\frac{d{v}_{n}}{dt}=\oint\limits_{\partial w(t)}{\phi}_{m}({\textbf{{n}}}\cdot{\textbf{{u}}}^{*})da-\int\limits_{\mathit{\Omega}}{\textbf{{u}}}\cdot\left(({\text{J}}^{-1})^{T}\nabla_{\xi}{\phi}_{m}\right)jd{\mathit{\Omega}},

(22)

{\text{M}}_{mm}\frac{d{\textbf{{u}}}_{n}}{dt}=\oint\limits_{\partial w(t)}{\phi}_{m}({\upsigma}^{*}{\textbf{{n}}})da-\int\limits_{\mathit{\Omega}}{\upsigma}\left(({\text{J}}^{-1})^{T}\nabla_{\xi}{\phi}_{m}\right)jd{\mathit{\Omega}},

(23)

{\text{M}}_{mm}\frac{d{\tau}_{n}}{dt}=\oint\limits_{\partial w(t)}{\phi}_{m}{\textbf{{n}}}\cdot({\upsigma}^{*}\textbf{{u}}^{*})da-\int\limits_{\mathit{\Omega}}({\upsigma}\textbf{{u}})\cdot\left(({\text{J}}^{-1})^{T}\nabla_{\xi}{\phi}_{m}\right)jd{\mathit{\Omega}}.

(24)

For the $rhs$ , the first (i.e., the surface integral) and second (i.e., the volume integral) terms are evaluated by Gauss quadrature formulas. This will be discussed in detail in Section 3.4. All the volume integrals can be calculated on the reference cell by

dw=jd\Omega,

(25)

where $j$ is the determinant of the Jacobian matrix J.

The discussion will now focus on the DG approach applied to the deformation gradient. Multiplying the evolution equation by the test function and integrating over the initial cell gives,

\int\limits_{w^{0}_{h}}\hat{\phi}_{m}\left(\frac{d{\text{F}}_{h}}{dt}-\nabla_{\textbf{{X}}}{\textbf{{u}}}\right)dw^{0}=0,\qquad m=1,...,N(P).

By introducing the volume matrix for the cell, ${V}^{0}_{mn}=\int\limits_{\Omega_{h}}{\hat{\phi}_{m}\hat{\phi}_{n}}dw^{0}$ , the previous system of equations reduces to

{V}^{0}_{mm}\frac{d\text{F}_{m}}{dt}=\oint\limits_{\partial w_{h}^{0}}\hat{\phi}_{m}(\textbf{{n}}\cdot\textbf{{u}}^{*})dA-\int\limits_{\mathit{\Omega}}\textbf{{u}}\cdot(({\text{J}}^{-1})^{T}{\nabla_{\xi}\,\hat{\phi}_{m}})j^{0}d{\mathit{\Omega}},\quad\quad m=1,...,N(P)

(26)

where ${V}^{0}_{mm}$ is the $m^{th}$ row and $m^{th}$ column of the mass matrix. There is no matrix to invert. The details on the derivation of the DG evolution equations with Taylor basis functions for $v$ , $\textbf{{u}}_{h}$ , $\tau_{h}$ , and $\text{F}_{h}$ can be found in papers [60, 64, 65].

3.3 Temporal integration

In this work, explicit strong-stability-preserving (SSP) Runge-Kutta (RK) schemes [68, 69, 70, 15] are used to evolve the semi-discrete equations for each polynomial coefficient from time level $n$ to $n+1$ . We consider the three-stage third-order accurate TVD RK method (denoted as SSPRK(3,3) in the paper) and the five-stage fourth-order accurate SSP RK method [70] (denoted as SSPRK(5,4) in the paper). The SSPRK(3,3) time integration has been applied to the third-order accurate Lagrangian DG(P2) method and is discussed by a range of papers, see [48, 60, 64]. The equation to evolve a polynomial coefficient ${\mathbb{U}}_{k}$ has the general form of

\frac{d{\mathbb{U}}_{k}}{dt}=\frac{1}{{\text{A}}_{kk}}\mathbb{R}_{k}

(27)

where ${\text{A}}_{kk}=\text{M}_{kk}$ when evolving the coefficients for $v$ , u, and $\tau$ ; likewise, $\text{A}_{kk}=\text{V}_{kk}^{0}$ when evolving the coefficients for F. RK methods for advancing ${\mathbb{U}}_{k}$ forward in time from time level $n$ to $n+1$ can be written as

\begin{array}[]{lll}{\mathbb{U}}^{(0)}_{k}&=&{\mathbb{U}}^{n}_{k}\\ {\mathbb{U}}^{(i)}_{k}&=&\sum\limits_{j=0}^{i-1}\alpha^{i,j}{\mathbb{U}}^{(j)}_{k}+\frac{\Delta t}{{\text{A}}_{kk}}\sum\limits_{j=0}^{i-1}\beta^{i,j}{\mathbb{R}_{k}^{(j)}},\quad i=1,...,s\\ {\mathbb{U}}^{n+1}_{k}&=&{\mathbb{U}}^{(s)}_{k}.\end{array}

(28)

The superscripts $0$ , $i$ , and $j$ denote the solution at a RK step, where $s$ is the number of stages. Table 1 presents the coefficients $\alpha^{i,j}$ and $\beta^{i,j}$ . The mesh position is temporally advanced using the same RK method as used to advance the polynomial coefficients.

\begin{array}[]{lll}{\textbf{{x}}}^{(0)}_{V}&=&{\textbf{{x}}}^{n}_{V}\\ {\textbf{{x}}}^{(i)}_{V}&=&\sum\limits_{j=0}^{i-1}\alpha^{i,j}{\textbf{{x}}}^{(j)}_{V}+\Delta t\sum\limits_{j=0}^{i-1}\beta^{i,j}{\textbf{{u}}}_{V}^{*(j)},\\ {\textbf{{x}}}^{n+1}_{V}&=&{\textbf{{x}}}^{(s)}_{V}.\end{array}

(29)

where ${\textbf{{x}}}_{V}$ is the position of a vertex and ${\textbf{{u}}}_{V}^{*}$ is the velocity of the vertex that comes from the Riemann solver.

In the context of DG(P3) methods, we found no obvious advantage for the SSPRK(5,4) over SSPRK(3,3) in terms of the numerical accuracy on test cases. For instance, the accuracy of the calculations of the Taylor-Green vortex test problem (Section 4.1) was the same with both time integration methods for the mesh resolutions and time step sizes considered. As such, robust, accurate, and computationally efficient solutions are possible using SSPRK(3,3) with this DG(P3) method on cubic cells.

Table 1: The coefficients for two RK time integration methods are shown. SSPRK(3,3) has three steps and is third-order accurate. SSPRK(5,4) has five steps and is fourth-order accurate.

SSPRK(3,3), Stages s=3
$\alpha^{i,j}$	1
	3/4	1/4
	1/3	0	2/3
$\beta^{i,j}$	1
	0	1/4
	0	0	2/3
SSPRK(5,4), Stages s=5
$\alpha^{i,j}$	1
	0.444370493651235	0.555629506348765
	0.620101851488403	0	0.379898148511597
	0.178079954393132	0	0	0.821920045606868
	0	0	0.517231671970585	0.096059710526147	0.386708617503269
$\beta^{i,j}$	0.391752226571890
	0	0.368410593050371
	0	0	0.251891774271694
	0	0	0	0.544974750228521
	0	0	0	0.063692468666290	0.226007483236906

3.4 Quadrature rules

Gauss quadrature formulas are used to evaluate the surface and volume integrals in the DG equations. The quadrature rules are chosen to exactly integrate the integrals, where the number of quadrature points depends on both the polynomial of degree $P$ and the order of the cell. On a linear cell, the number of quadrature points is selected to integrate exactly polynomials of order $2P+1$ on the reference cell. For quadratic and cubic meshes, the number of the quadrature points increases by 1 and 2 respectively in theory. The number of quadrature points used to solve the surface and volume integrals is shown in Table 2.

For Lagrangian hydrodynamics, it is essential to use a nodal Riemann solver to calculate a vertex velocity ${\textbf{{u}}}_{V}^{*}$ that is then used to spatially move the vertex (i.e., for the vertices at the corner and along the edge). Therefore, Gauss-Lobatto quadrature rules are used to calculate the surface integrals as they will spatially overlap some or all of the vertices on the cell surface. For Lagrangian hydrodynamics, the mesh is updated every time step. So another issue is how to guarantee the vertices always overlap the quadrature points especially for the edge vertices. Due to a favorable property in Lagrangian hydrodynamics,

\frac{d\boldsymbol{\xi}}{dt}=0,

(30)

so the reference coordinates ${\boldsymbol{\xi}}$ for any quadrature point is the same as the initial one. In other words, as long as the initial vertex overlaps the Gauss-Lobatto quadrature point, then the vertex can always be used as the quadrature point (Table 3). Table B1 in the Appendix presents the Gauss-Lobatto quadrature point locations and values. The volume integrals are calculated using the Gauss-Legendre quadrature rules because they are more accurate than the Gauss-Lobatto rules using the same number of quadrature points.

For higher-order DG methods, it is necessary to solve 1D Riemann problems at multiple locations along the edge in addition to that at the vertices. As such, the reference coordinate locations for the surface vertices will coincide with some of the Gauss-Lobatto quadrature point locations, see Table 3. For cubic meshes, the number of required Lobatto quadrature points to exactly integrate the $rhs$ varies, meaning that the initial cubic meshes are slightly different for DG methods with different polynomial degrees. The details on the Riemann solvers will be discussed in the following subsections.

Table 2: The information for the Gauss quadrature point number is shown in the context of DG methods. Gauss-Lobatto quadrature points are used for the surface integral because they spatially overlap the vertices where a nodal Riemann solver is used. The Gauss-Legendre quadrature points are adopted for the volume integral due to their high accuracy. For the surface integral, the number in the parentheses is not used since we want to use the Riemann information at the vertices.

DG(P1) DG(P2) DG(P3) Surface integral Linear edge 3 4 5 Quadratic edge 5 (4) 5 6 Cubic edge 6 (5) 6 7 Volume integral Linear cell 2 3 4 Quadratic cell 3 4 5 Cubic cell 4 5 6

Table 3: The position in the reference coordinates for each surface vertex position is shown for differing orders. The reference coordinates are over interval

[-1,1]

. The vertices overlap all or only some of the surface quadrature points. For a cubic edge, different quadrature rules are used for DG methods with differing polynomial orders, which leads to using different initial cubic meshes.

① ② ③ ④ Linear edge -1 1 Quadratic edge -1 0 1 Cubic edge 6-point quadrature rule -1 $-\sqrt{\frac{1}{21}\left(7-2\sqrt{7}\right)}$ $\sqrt{\frac{1}{21}\left(7-2\sqrt{7}\right)}$ 1 7-point quadrature rule -1 $-\sqrt{\frac{5}{11}-\frac{2}{11}\sqrt{\frac{5}{3}}}$ $\sqrt{\frac{5}{11}-\frac{2}{11}\sqrt{\frac{5}{3}}}$ 1

3.4.1 Riemann solvers

The surface integral requires solving Riemann problems at the quadrature points along the cell surface. In this work, a multi-directional approximate Riemann solver (MARS) is used at the vertices of the cell and a 1D Riemann solver is used at additional locations along the cell edge. MARS are widely used in the Lagrangian finite volume CCH where the first one was proposed by B. Després and C. Mazeran [9] for gas dynamics problems. Maire and various authors [11, 12] extended the work in [9] and proposed a new MARS that had improved mesh robustness properties. Burton et al. [13] proposed another robust MARS that was suitable for materials with strength. We use the MARS by Morgan et al. [59] and the discussion that follows will briefly describe that MARS.

The force acting on a surface segment connected to the vertex is given by

{\textbf{{F}}}_{i}^{*}=a_{i}{\upsigma}_{i}^{*}{\textbf{{n}}}_{i}=a_{i}{\upsigma}_{c}{\textbf{{n}}}_{i}+\mu_{i}a_{i}\left({{\textbf{{u}}}_{V}^{*}-{\textbf{{u}}}_{c}}\right).

(31)

The subscript $c$ denotes the value in the cell at the vertex, the subscript $i$ denotes the surface segment, the subscript $V$ is a vertex value, and the superscript $*$ denotes a variable that comes from solving a Riemann problem. The shock impedance [71] is

\mu_{i}=\rho_{z}(c_{z}+\frac{\gamma+1}{2}|({{\textbf{{u}}}_{V}^{*}-{\textbf{{u}}}_{c}})\cdot{\textbf{{n}}}_{i}|),

(32)

where the subscript $z$ denotes the mass average value. Momentum conservation at the vertex $V$ requires that the summation of all forces around the vertex to be equal to zero. A single Riemann velocity is found at the vertex that guarantees momentum conservation.

{\textbf{{u}}}_{V}^{*}=\frac{{\sum\limits_{i\in V}{\left({\mu_{i}a_{i}{\textbf{{u}}}_{c}-a_{i}{\boldsymbol{\upsigma}_{c}}}{\textbf{{n}}}_{i}\right)}}}{{\sum\limits_{i\in V}{\mu_{i}a_{i}}}}.

(33)

Analysis is presented in [50] to show that the Lagrangian DG method, combined with this Riemann solver, satisfies the second-law of thermodynamics. We also use the subcell surfaces (Section 3.5) with the MARS so that there is consistency in solving for the Riemann velocity (and associated forces) between the corner and edge vertices (Fig. 3).

The accuracy of the surface integral for the DG(P3) method is increased by using more quadrature points along the edge. The goal is to increase the quadrature order without adding more kinematic degrees of freedom to the cell. The resulting Riemann solutions correspond to the Gauss-Lobatto quadrature rule where the integration weights and locations in the reference coordinate system are provided in Table B1. The 6-point quadrature rule is exact for integrating a polynomial of order $9$ , and is substantially more accurate than using the 4 points that coincide with the vertices along the cell edge. The discussion that follows will focus on the 1D Riemann solutions at the non-vertex quadrature points that are illustrated in Fig. 2.

For the 1D Riemann problem, there will be two forces acting on a surface segment of the cell.

\begin{split}{\textbf{{F}}}_{1}^{*}=a_{1}{\upsigma}_{1}^{*}{\textbf{{n}}}_{1}=a_{1}{\upsigma}_{c1}{\textbf{{n}}}_{1}+\mu_{1}a_{1}\left({{\textbf{{u}}}_{G}^{*}-{\textbf{{u}}}_{c1}}\right)\\ {\textbf{{F}}}_{2}^{*}=a_{2}{\upsigma}_{2}^{*}{\textbf{{n}}}_{2}=a_{2}{\upsigma}_{c2}{\textbf{{n}}}_{2}+\mu_{2}a_{2}\left({{\textbf{{u}}}_{G}^{*}-{\textbf{{u}}}_{c2}}\right)\end{split}

(34)

Here, the subscripts $1$ and $2$ denote the respective inputs from the two cells that share this surface segment, and ${\textbf{{u}}}_{G}^{*}$ denotes the Riemann velocity at the non-vertex Gauss quadrature point. The surface weighted area normals for the segment are equal and opposite, so $a_{1}{\textbf{{n}}}_{1}=-a_{2}{\textbf{{n}}}_{2}$ . A single Riemann velocity ${\textbf{{u}}}_{G}^{*}$ can be calculated by enforcing momentum conservation $\sum\limits_{i}{\textbf{{F}}}_{i}^{*}=0$ . For a non-vertex Gauss quadrature point, the calculated Riemann velocity ${\textbf{{u}}}_{G}^{*}$ is not necessarily equal to the velocity of the edge ${\textbf{{u}}_{G}}$ . Along an edge, the velocity ${\textbf{{u}}_{G}}$ at the non-vertex Gauss quadrature point is calculated using Lagrangian interpolation where the coefficients are the Riemann velocities at the vertices of this edge ${\textbf{{u}}}_{V}^{*}$ . The fluxes in the specific volume evolution equation at the non-vertex quadrature points are calculated using the edge velocity ${\textbf{{u}}_{G}}$ . The fluxes in the specific total energy evolution equation at the non-vertex quadrature points are calculated by multiplying the forces (Eq. 34) by the edge velocity ${\textbf{{u}}_{G}}$ ,

\begin{split}{\textbf{{F}}}_{1}^{*}\cdot{\textbf{{u}}_{G}}=a_{1}{\upsigma}_{1}^{*}{\textbf{{n}}}_{1}\cdot{\textbf{{u}}_{G}}\\ {\textbf{{F}}}_{2}^{*}\cdot{\textbf{{u}}_{G}}=a_{2}{\upsigma}_{2}^{*}{\textbf{{n}}}_{2}\cdot{\textbf{{u}}_{G}}\end{split}

(35)

These total energy fluxes guarantee the conservation of total energy.

3.5 Subcell Mesh Stabilization (SMS)

In this work, the SMS scheme [60] is extended to work with cubic cells. The goal of SMS is to remove spurious vorticity that can arise on flows with strong shocks and to obtain a more consistent Riemann solution between the corner and edge vertices. For a corner vertex $V$ shown in Fig. 3, there are eight segments contributing to the Riemann problem. This motivates us to use more segments to construct the Riemann problem for an edge vertex by creating subcells that surround the edge vertex. A cubic cell is decomposed into nine subcells (Fig. 3). The subcell surfaces are used in the MARS at the edge vertices. The subcell decomposition also enables us to introduce corrective pressures in the MARS to remove spurious vorticity. The concept in SMS is to calculate a density for each subcell and then correct the pressure inputs to the MARS at every exterior vertex if the subcell density deviates from the modal density field for the cell. The subscells move in a Lagrangian manner so the average density of a subcell is

\bar{\rho}_{s}=\frac{m_{s}}{w_{s}(t)}=\frac{\int\limits_{\Omega_{s}}\rho({\boldsymbol{\xi}},t^{0})j_{s}^{0}d\mathit{\Omega}}{\int\limits_{\Omega_{s}}j_{s}d\mathit{\Omega}}.

(36)

The subscript $s$ denotes the variables for the subcell, the determinant of the Jacobian matrix for the subcell is denoted as $j_{s}$ and the superscript $0$ denotes the initial configuration. Inconsistencies between the cell density field and the subcell density may arise. The average density in the subcell using the specific volume field $v_{h}$ that is evolved with the DG approach is

\bar{\rho}_{vs}=\frac{\int\limits_{w_{s}(t)}\frac{1}{v_{h}}dw}{w_{s}(t)}=\frac{\int\limits_{\Omega_{s}}\frac{1}{v_{h}}j_{s}d{\Omega}}{w_{s}(t)},

(37)

SMS adds a correction term to the density used to calculate the pressure used in the MARS,

\hat{\rho}_{c}=\rho_{c}+\chi\left(\bar{\rho}_{s}-\bar{\rho}_{vs}\right),

(38)

where $\hat{\rho}_{c}$ is the corrected density for all exterior corners of the subcell (i.e., in the corner at an edge vertex and in the cell corner), $\bar{\rho}_{s}$ is the average density of a subcell based on mass conservation (Eq. 36), $\bar{\rho}_{vs}$ is the average density of the subcell based on the specific volume $v_{h}$ distribution for the cell, and $\chi$ is a user settable coefficient in the range of zero to one. If the cell deforms in a manner consistent with the specific volume field then the density correction term is zero, and higher-order accuracy can be achieved with the Lagrangian DG method, i.e., fourth-order accuracy with DG(P3).

3.5.1 Surface integral

This section focuses on the evaluation of the surface integral in Eq. 22, 23, and 24. Taking Eq. 22 as an example, for the 7-point Lobatto quadrature rule on cubic meshes,

\begin{array}[]{l}\int_{V_{3}V_{4}}({\textbf{{n}}}\cdot{\textbf{{u}}}^{*})da=(a_{1}{\textbf{{n}}}_{1})_{G_{1}}\cdot{\textbf{{u}}}_{G_{1}}+(a_{1}{\textbf{{n}}}_{1})_{G_{2}}\cdot{\textbf{{u}}}_{G_{2}}+\sum_{i=1}^{2}(a_{i}{\textbf{{n}}}_{i})_{G_{3}^{r}}\cdot{\textbf{{u}}}_{G_{3}}+\sum_{i=1}^{2}(a_{i}{\textbf{{n}}}_{i})_{G_{3}^{l}}\cdot{\textbf{{u}}}_{G_{3}}\\ +(a_{1}{\textbf{{n}}}_{1})_{G_{4}}\cdot{\textbf{{u}}}_{G_{4}}+\sum_{i=1}^{2}(a_{i}{\textbf{{n}}}_{i})_{G_{5}^{r}}\cdot{\textbf{{u}}}_{G_{5}}+\sum_{i=1}^{2}(a_{i}{\textbf{{n}}}_{i})_{G_{5}^{l}}\cdot{\textbf{{u}}}_{G_{5}}+(a_{1}{\textbf{{n}}}_{1})_{G_{6}}\cdot{\textbf{{u}}}_{G_{6}}+(a_{2}{\textbf{{n}}}_{2})_{G_{7}}\cdot{\textbf{{u}}}_{G_{7}}.\\ \end{array}

(39)

Here, four quadrature points will overlap the vertices, the velocity of the edge ${\textbf{{u}}}_{G}$ is equal to the Riemann velocity at the vertices ${\textbf{{u}}}^{*}_{v}$ , while for the non-vertex quadrature points, the velocity of the edge ${\textbf{{u}}}_{G}$ has been specified in Section 3.4.1. Taking the area normal vector $(a_{2}{\textbf{{n}}}_{2})_{G_{7}}$ as an example, the subscript $G_{7}$ denotes the corner associated with the quadrature point $G_{7}$ . And for the edge vertex $G_{3}$ and $G_{5}$ , the superscript $l$ and $r$ denote the left and right corner associated with the edge vertex $G_{3}$ and $G_{5}$ . The weights used in this 7-point Lobatto quadrature rule have been included in the weighted area $a_{i}$ . The surface area normal vector $s{\textbf{{n}}}$ along the surface $V_{3}V_{4}$ can be expressed as,

\begin{array}[]{l}\int_{V_{3}V_{4}}d{\textbf{{s}}}=\sum_{g}wt_{g}(s{\textbf{{n}}})_{g},\\ \end{array}

(40)

where $s_{g}$ represent the determinant of the Jacobian matrix at the corresponding quadrature point along the surface $V_{3}V_{4}$ . Combining Eq. 40 with 39, we can conclude that

\begin{array}[]{l}(a_{1}{\textbf{{n}}}_{1})_{G_{1}}=wt_{G_{1}}(s{\textbf{{n}}})_{G_{1}},\\ (a_{1}{\textbf{{n}}}_{1})_{G_{2}}=wt_{G_{2}}(s{\textbf{{n}}})_{G_{2}},\\ (a_{2}{\textbf{{n}}}_{2})_{G_{3}^{l}}=\frac{1}{2}wt_{G_{3}}(s{\textbf{{n}}})_{G_{3}^{l}},\\ (a_{1}{\textbf{{n}}}_{1})_{G_{3}^{r}}=\frac{1}{2}wt_{G_{3}}(s{\textbf{{n}}})_{G_{3}^{r}},\\ (a_{1}{\textbf{{n}}}_{1})_{G_{4}}=wt_{G_{4}}(s{\textbf{{n}}})_{G_{4}},\\ (a_{2}{\textbf{{n}}}_{2})_{G_{5}^{l}}=\frac{1}{2}wt_{G_{5}}(s{\textbf{{n}}})_{G_{5}^{l}},\\ (a_{1}{\textbf{{n}}}_{1})_{G_{5}^{r}}=\frac{1}{2}wt_{G_{5}}(s{\textbf{{n}}})_{G_{5}^{r}},\\ (a_{1}{\textbf{{n}}}_{1})_{G_{6}}=wt_{G_{6}}(s{\textbf{{n}}})_{G_{6}},\\ (a_{2}{\textbf{{n}}}_{2})_{G_{7}}=wt_{G_{7}}(s{\textbf{{n}}})_{G_{7}}.\\ \end{array}

(41)

In addition, at the edge vertex $V_{9}$ and $V_{10}$ , we need to pay special attention to segment 2 of the right corner and segment 1 of the left corner in Fig. 2. Likewise,

\begin{array}[]{l}(a_{1}{\textbf{{n}}}_{1})_{V_{9}^{l}}=-(a_{2}{\textbf{{n}}}_{2})_{V_{9}^{r}}=\frac{1}{2}wt_{G_{3}}(s{\textbf{{n}}})_{9t},\\ (a_{1}{\textbf{{n}}}_{1})_{V_{10}^{l}}=-(a_{2}{\textbf{{n}}}_{2})_{V_{10}^{r}}=\frac{1}{2}wt_{G_{5}}(s{\textbf{{n}}})_{10t},\\ \end{array}

(42)

Here, $s_{9t}$ means the determinant of the Jacobian matrix at the vertex $V_{9}$ from the inside surface $V_{6}V_{9}$ . Please refer to [60] for the choice of the weights (i.e., $\frac{1}{2}wt_{G_{3}}$ and $\frac{1}{2}wt_{G_{5}}$ ).

3.6 Limiting

The modal fields and the pressure field (a nodal field) are limited toward the cell average value (i.e., piece-wise constant over the cell) near shocks and discontinuities, which is an essential ingredient for robust solutions. Therefore, the limiting strategy is to deliver high-order (e.g., fourth-order) accuracy in smooth regions of a calculation and to reduce the accuracy towards first-order near discontinuities. A favorable feature in Lagrangian hydrodynamics is that the mesh clusters near the shock, providing natural $h$ refinement for the shock region. Simultaneously, the discretization order is reduced to DG(P1) to guarantee the monotonicity. A merit of using new orthogonal basis functions is that each moment in the Pn polynomial are decoupled from all other terms, allowing flexibility in limiting of each coefficient. Reducing the Pm polynomial space to the Pn one ( $m\geq n$ ) can be done by truncating the higher-order terms that are higher than the Pn polynomial. We found this aforementioned truncation in the context of the Taylor-series polynomials that are higher than P2 generated extremely poor results because the moments are coupled through the mass matrix. Using the orthogonal basis functions was essential for obtaining robust solutions on strong shock problems.

To ensure high-order solutions away from a shock, we use a troubled-cell detector [72, 73, 74, 75, 76, 77, 78] to identify the cells near a discontinuity. In the present work, the troubled-cell detector from [74] is adopted; this detector was successfully used in [60] with a Lagrangian DG(P2) method.

If a cell is tagged as a troubled-cell, then the quadratic and cubic terms are set to zero, and the first-order derivatives in the modal fields ( $v_{h}$ , $\textbf{{u}}_{h}$ , and $\tau_{h}$ ) are limited to ensure the reconstructions at the vertices are bounded by the vertex-neighboring cell-averages just like Lagrangian FV CCH methods [11, 12, 13, 14]. The first-order derivatives in the velocity polynomial reconstructions are limited in either the principle strain direction [79, 13] or the local flow direction [6] to ensure the preservation of symmetry on equal-angle polar meshes. For a troubled-cell, the pressure value in each corner was limited to be within the bounds of the cell-averages of all adjacent cells to the vertex. Additional details on the limiting of the first derivatives for scalar and vector fields with the Lagrangian DG(P1) and DG(P2) methods are provided in [50, 60] respectively.

4 Test problems

In this section, a diverse suite of test problems covering gas and solid dynamics are calculated to demonstrate the accuracy and robustness of this Lagrangian DG hydrodynamic method with a hierarchical orthogonal basis on curvilinear meshes. Both smooth flows and shock driven flows are calculated. We will now provide a brief introduction of each test problem.

The Taylor-Green vortex test case (Section 4.1) and the Gresho vortex test case (Section 4.2) are smooth, vortical flow problems. We use these test cases to quantify the convergence rate (i.e, order of accuracy) of the new Lagrangian DG method using the $L_{2}$ error norm. We will briefly describe the error analysis used on these test problems. Taking pressure $p$ as an example, the $L_{2}$ error norm is defined by,

||p-p_{e}||_{L_{2}}=\sqrt{\sum\limits_{i=1}^{Num}\int_{w}(p-p_{e})^{2}d{w}},\\

(43)

where $p_{e}$ denotes the exact pressure, and $ncell$ is the number of cells in the problem. For all convergence studies, we use $l1$ , $l2$ , $l3$ and $l4$ to denote the various levels of mesh resolution where $l1$ is the coarsest resolution.

The Sedov blast wave test case (Section 4.3) is used to demonstrate the accuracy and stability of the Lagrangian DG method with an exceptionally strong shock. We then present results for an elastic vibration of a beryllium plate test (Section 4.4) problem and a planar Taylor Anvil impact test (Section 4.5) to show the ability of this new scheme to simulate challenging solid dynamic problems.

Curvilinear box meshes are used for all calculations, unless stated otherwise. A few linear mesh results are included for comparison purposes. For a box grid, $40\times 20$ denotes the number of cells along the $x$ and $y$ direction. The initial box grid has been specified in Section 3.4.

4.1 2D Taylor-Green vortex problem

The Taylor-Green vortex problem is a shockless, smooth flow problem [80, 81, 47, 82], which is of particular interest because it is a vortical flow problem with an analytical solution. The material is a gamma-law gas with $\gamma=7/5$ . The initial flow field is

\begin{array}[]{l l l}\rho^{0}&=&1,\\ u_{x}^{0}&=&sin(\pi x)cos(\pi y),\\ u_{y}^{0}&=&-cos(\pi x)sin(\pi y),\\ p^{0}&=&\frac{1}{4}[cos(2\pi x)+cos(2\pi y)]+1,\end{array}

where $u_{x}$ and $u_{y}$ denote the $x$ and $y$ components of velocity, the superscript $0$ represents the initial state. The specific total energy of the field is

\tau=\frac{p}{\rho(\gamma-1)}+\frac{1}{2}(u_{x}^{2}+u_{y}^{2}).\\

An energy source term is included in the total energy equation to maintain a steady-state solution for the compressible inviscid case,

s_{\tau}=\frac{\pi}{4(\gamma-1)}[cos(3\pi x)cos(\pi y)-cos(3\pi y)cos(\pi x)].\\

The computational domain for this test case is defined by $(x,y)=[0,1]\times[0,1]$ . All calculations use a set of uniformly refined quadrilateral meshes, namely $5\times 5$ , $10\times 10$ , $20\times 20$ , and $40\times 40$ respectively.

The meshes and pressure fields at $t=0.75$ using DG(P3) with cubic cells are shown in Fig. 4. DG(P3) delivers very accurate flow details for the pressure field even with a very coarse mesh resolution ( $5\times 5$ ). It can also be observed that the mesh is able to deform with the flow.

Using the analytical solution, the Taylor-Green vortex problem is used to assess the order of accuracy of this Lagrangian DG method with curvilinear meshes. The $L_{2}$ norm of the error at $t=0.1$ and $t=0.4$ is calculated for the density ( $\rho$ ), velocity component ( $u_{x}$ ), pressure ( $p$ ) and total energy field ( $\tau$ ). The associated results are shown in Table 4. The DG(P2) and DG(P3) methods can yield nearly perfect third-order and fourth-order accuracy at $t=0.1$ for all the variables, respectively. With time marching, the mesh has greater deformation and the numerical error accrues. Accordingly, at $t=0.4$ , the convergence rate (i.e., order of accuracy) for all the variables is slightly less than the designed order of accuracy for both DG(P2) and DG(P3). Fig. 5 presents the numerical errors for all the variables at different times vs. the computational costs (i.e., number of degrees of freedom, i.e., nDoF) using DG(P2) and DG(P3). This example demonstrates that higher-order methods can be computationally more efficient to achieve a given level of accuracy compared with lower-order methods.

Finally, we also calculated this problem with DG(P2) using linear and cubic meshes respectively. The numerical error for all the variables are shown in Table 5. DG(P2) with linear meshes delivers second-order accuracy at most for all the variables while DG(P2) with cubic meshes can keep third-order accuracy. In addition, in the context of the DG(P2) method, the numerical error and convergence rate on the cubic meshes is close to that on quadratic meshes. This demonstrates that achieving the $(n+1)^{th}$ order of accuracy requires using quadrilaterals meshes with edges defined by a polynomial of degree $n$ with a DG(Pn) method.

Table 4: Numerical error and convergence rate are shown for the Taylor-Green vortex problem at different times using DG(P2) and DG(P3) on curvilinear meshes. At the early time

t=0.1

, for all the variables (i.e., the density

\rho

, velocity component

u_{x}

, pressure

p

and total energy

\tau

), DG(P2) and DG(P3) can deliver approximately the designed third- and fourth-order accuracy, respectively. While at a later time

t=0.4

, for both DG(P2) and DG(P3), all the variables converge slightly less than the designed order of accuracy due to the acquirement of numerical errors with time marching.

Mesh $\rho$ $u_{x}$ $p$ $\tau$ $L_{2}$ error order $L_{2}$ error order $L_{2}$ error order $L_{2}$ error order DG(P2), t=0.1 $5\times 5$ 1.5995e-3 2.7796e-3 3.2363e-3 6.7452e-3 $10\times 10$ 2.3499e-4 2.77 3.6141e-4 2.95 4.3700e-4 2.90 8.9561e-4 2.92 $20\times 20$ 3.0535e-5 2.96 4.4599e-5 3.03 5.4800e-5 3.00 1.1491e-4 2.97 $40\times 40$ 4.0233e-6 2.93 5.2717e-6 3.09 6.6834e-6 3.04 1.4758e-5 2.97 DG(P2), t=0.4 $5\times 5$ 7.6054e-3 8.2759e-3 1.2875e-2 2.7769e-2 $10\times 10$ 1.7285e-3 2.14 1.7996e-3 2.21 3.2562e-3 1.99 5.4817e-3 2.35 $20\times 20$ 2.7914e-4 2.63 3.1000e-4 2.55 4.8153e-4 2.77 7.7198e-4 2.84 $40\times 40$ 5.0151e-5 2.48 4.6679e-5 2.74 7.2271e-5 2.75 1.3387e-4 2.54 DG(P3), t=0.1 $5\times 5$ 1.3801e-4 2.6790e-4 3.7398e-4 8.3408e-4 $10\times 10$ 1.0177e-5 3.77 1.7242e-5 3.97 2.5624e-5 3.88 5.5205e-5 3.93 $20\times 20$ 7.3830e-7 3.80 1.1475e-6 3.93 1.7080e-6 3.92 3.5667e-6 3.97 $40\times 40$ 5.7127e-8 3.71 8.3370e-7 3.80 1.1058e-7 3.96 2.3785e-7 3.92 DG(P3), t=0.4 $5\times 5$ 1.3688e-3 2.5387e-3 4.4911e-3 7.9309e-3 $10\times 10$ 1.6847e-4 3.03 2.8333e-4 3.17 4.3188e-4 3.39 7.3365e-4 3.45 $20\times 20$ 1.8453e-5 3.20 2.6109e-5 3.45 4.6647e-5 3.22 7.6164e-5 3.28 $40\times 40$ 1.9066e-6 3.28 2.6566e-6 3.31 4.2518e-6 3.47 7.0218e-6 3.45

Table 5: Numerical error and convergence rate are shown for the Taylor-Green vortex problem at

t=0.1

using DG(P2) on both linear and cubic meshes. For all the variables (i.e., the density

\rho

, velocity component

u_{x}

, pressure

p

and total energy

\tau

), DG(P2) delivers second-order accuracy at most on linear meshes. DG(P2) is able to achieve third-order accuracy for all the variables on cubic meshes, and the difference from DG(P2) on quadratic meshes is very small.

Mesh $\rho$ $u_{x}$ $p$ $\tau$ $L_{2}$ error order $L_{2}$ error order $L_{2}$ error order $L_{2}$ error order DG(P2), linear, t=0.1 $5\times 5$ 5.4851e-3 9.0759e-3 7.5490e-3 1.1949e-2 $10\times 10$ 1.0247e-3 2.43 2.2022e-3 2.05 1.1707e-3 2.70 2.7164e-3 2.14 $20\times 20$ 2.7111e-4 1.93 5.4123e-4 2.03 2.3301e-4 2.34 8.0907e-4 1.75 $40\times 40$ 1.1079e-4 1.30 1.3557e-4 2.00 5.5829e-5 2.07 3.1164e-4 1.38 DG(P2), cubic, t=0.1 $5\times 5$ 1.7066e-3 2.8481e-3 3.2726e-3 6.8003e-3 $10\times 10$ 2.4786e-4 2.79 3.8802e-4 2.89 4.4973e-4 2.87 8.8884e-4 2.95 $20\times 20$ 3.2821e-5 2.93 2.3186e-5 2.88 5.9167e-5 2.94 1.1354e-4 2.98 $40\times 40$ 4.3246e-6 2.93 7.7855e-6 2.78 7.6519e-6 2.96 1.4637e-5 2.97

4.2 2D Gresho vortex

The Gresho vortex problem is a smooth vortical flow problem. The velocity in the angular direction $u_{\theta}$ is solely a function of the radius $r$ . The pressure gradient balances the centrifugal forces. The material is a gamma-law gas with $\gamma=7/5$ . The initial flow field is

\begin{array}[]{l l l}\rho^{0}&=&1,\\ u_{r}^{0}&=&0,\\ u_{\theta}^{0}&=&\left\{\begin{array}[]{lll}5r&\textup{if}&0\leq r<0.2,\\ 2-5r&\textup{if}&0.2\leq r<0.4,\\ 0&\textup{if}&0.4\leq r.\\ \end{array}\right.\\ p^{0}&=&\left\{\begin{array}[]{lll}5+\frac{25}{2}r^{2}&\textup{if}&0\leq r<0.2,\\ 9-4ln(0.2+)5+\frac{25}{2}r^{2}-20r+4lnr&\textup{if}&0.2\leq r<0.4,\\ 3+4ln2&\textup{if}&0.4\leq r.\\ \end{array}\right.\\ \end{array}

Here, $u_{r}$ and $u_{\theta}$ denote the radial and tangential components of the velocity. Unlike the Taylor-Green vortex case, the source term is not needed to maintain a steady-state solution for this compressible inviscid case. The density field is uniform and has a value of 1. The Gresho vortex is calculated to a time of 0.62. The computational domain for this test case is defined by $(x,y)=[-0.5,0.5]\times[-0.5,0.5]$ . A set of uniformly refined meshes is used to perform the grid convergence study. In our calculations, the inner and outer boundary are set to be stationary. The mesh resolutions are $16\times 16$ , $32\times 32$ and $64\times 64$ respectively.

The deformed cubic meshes and pressure field at time $t=0.4$ and $t=0.62$ are shown in Fig. 6 with DG(P3) using two mesh resolutions. Indeed the meshes are deforming in a robust and curvilinear manner. In addition, Fig. 7 presents the scatter plots for variables of interest, i.e., density ( $\rho$ ), velocity magnitude ( $v$ ) and pressure ( $p$ ), with DG(P2) and DG(P3). With mesh refinement, the numerical results agree better with the analytical solution for both DG(P2) and DG(P3). Obviously, DG(P3) provides much better agreement with the analytical solution. The time evolution until $t=0.62$ demonstrates the robustness of the Lagrangian DG hydrodynamic method with a hierarchical orthogonal basis.

With the analytical solution, the variant of Gresho vortex problem is used to assess the order of accuracy of this new method. We present the comparison for the numerical error and order of accuracy between DG(P2) and DG(P3). The associated results are shown in Table 6. For this test problem, it has been shown that the velocity is only $C_{1}$ continuous at $r=0.2$ and $r=0.4$ , indicating that numerical methods can deliver at most second-order accuracy. The results in Table 6 show that at $t=0$ both DG(P2) and DG(P3) deliver approximately second-order accuracy for all the variables except pressure, which is slightly higher. Likewise, at time $t=0.4$ , the order of accuracy for both DG(P2) and DG(P3) is still second-order at most, which is expected. Also Fig. 8 compares the efficiency between DG(P2) and DG(P3). The numerical error for variables vs. the computational costs (i.e., nDoF) are presented. It again shows that the higher-order DG(P3) method is more efficient since it can deliver superior accuracy at the same computational cost as the DG(P2) method, even on a problem that is second-order accurate at most.

Table 6: Numerical error and convergence rate are shown for the 2D Gresho vortex problem at different times using the DG method on curvilinear meshes. Even for the initial

L_{2}

projection, DG(P2) and DG(P3) deliver second-order of accuracy at most for velocity component

u_{x}

and pressure

p

. For a later time

t=0.4

, all the variables converge at approximately second-order for both DG(P2) and DG(P3). This is due to the fact that the exact solutions (i.e., velocity and total energy) are only

C^{1}

continuous.

Mesh $\rho$ $u_{x}$ $p$ $\tau$ $L_{2}$ error order $L_{2}$ error order $L_{2}$ error order $L_{2}$ error order DG(P2), t=0 $16\times 16$ 4.8991e-16 3.7666e-3 7.1518e-4 4.8080e-3 $32\times 32$ 4.6949e-16 - 1.2268e-3 1.62 1.2125e-4 2.57 1.3856e-3 1.80 $64\times 64$ 4.6486e-16 - 4.3483e-4 1.50 2.4526e-5 2.31 5.0456e-4 1.46 DG(P2), t=0.4 $16\times 16$ 3.9356e-3 1.5438e-2 1.8418e-2 3.5438e-2 $32\times 32$ 1.2165e-3 1.70 5.2000e-3 1.58 3.8883e-3 2.25 1.0854e-2 1.71 $64\times 64$ 3.8063e-4 1.68 1.8206e-3 1.52 8.8336e-4 2.14 3.6134e-3 1.59 DG(P3), t=0 $16\times 16$ 5.7498e-16 2.1138e-3 2.9084e-4 2.2668e-3 $32\times 32$ 5.6856e-16 - 8.2622e-4 1.37 5.5992e-5 2.39 9.6197e-4 1.24 $64\times 64$ 5.8615e-16 - 2.9127e-4 1.51 9.0702e-6 2.64 3.4154e-4 1.50 DG(P3), t=0.4 $16\times 16$ 1.6372e-3 6.7794e-3 3.9795e-3 1.4093e-2 $32\times 32$ 4.8989e-4 1.75 2.3066e-3 1.56 1.0115e-3 1.98 4.5845e-3 1.63 $64\times 64$ 1.1083e-4 2.15 8.8086e-4 1.39 2.0478e-4 2.31 1.3614e-3 1.76

4.3 Sedov blast problem

The Sedov test problem [83, 84] is an outward traveling blast wave in a gamma-law gas that is initiated by an energy source at the origin. In this work, $\gamma=7/5$ . The initial conditions are given by $(\rho^{0},u_{x}^{0},u_{y}^{0},p^{0})=(1.0,0,0,10^{-6})$ and there is an energy source at the origin. The pressure in the cell containing the origin is given by $p_{o}=(\gamma-1){\rho_{o}}\frac{E_{o}}{w_{o}}$ , where $w_{o}$ denotes the volume of the cell containing the origin and $E_{o}$ is the total amount of released internal energy. The energy source is selected to be equal to $0.244816$ so that the shock front is located at a radius equal to 1 at $t=1$ . The computational domain is defined in 2D Cartesian coordinates $[0,1.2]\times[0,1.2]$ . The mesh resolution is $30\times 30$ and initially uniform. Symmetry boundary conditions are imposed on the left and bottom boundaries. The outer boundaries are fixed. The limiting is performed in the local flow direction.

The final meshes and density contours are shown in Fig. 9 using this new Lagrangian DG hydrodynamic method. The meshes move in a stable and curvilinear manner. There is no spurious mesh motion [59] due to using this new hierarchical orthogonal basis and SMS method. The scatter plots of density at the final time are shown in Fig. 10. The results are in good agreement with the exact solution. From the close-up in Fig. 10b, the density scatter of DG(P3) is a little less than DG(P2). It is necessary to reiterate that the hierarchical orthogonal basis is crucial in the context of limiting for shock problems.

4.4 Elastic vibration of a Beryllium plate

Thin beams or plates have been the subject of many test problems, see for example [85]. The problem considered in this paper is a thick plate with free boundaries. The problem has been previously analyzed using PAGOSA [86], a three-dimensional Eulerian hydrodynamic code. Here the problem is solved in a Lagrangian framework as described in [80]. This is a two dimensional test problem comprised of an elastic beryllium plate with no supports or constraints. The computational domain is defined by $(x,y)\in[-3$ $cm,3$ $cm]\times[-0.5$ $cm,0.5$ $cm]$ . The centerline of the plate initially coincides with the axis $x$ . There is no analytic solution for a thick plate as specified for this problem. However, for a thin unconstrained plate, the solution for the first flexural mode is given by [80],

y(x,t)=Asin(\omega t)\big{(}g_{1}[sinh(\Omega(x+3))+sin(\Omega(x+3))]-g_{2}[cosh(\Omega(x+3))+cos(\Omega(x+3))]\big{)},

where $\Omega=0.7883401241$ $cm^{-1}$ , $\omega=0.2359739922$ ${\mu s}^{-1}$ , $A=0.004336850425$ $cm$ , $g_{1}=56.63685154$ ${\mu s}$ and $g_{2}=57.64552048$ ${\mu s}$ . Then the plate is prescribed with an initial velocity distribution as given by

\begin{array}[]{lll}u(x,t=0)&=&0,\\ v(x,t=0)&=&A\omega\big{(}g_{1}[sinh(\Omega(x+3))+sin(\Omega(x+3))]-g_{2}[cosh(\Omega(x+3))+cos(\Omega(x+3))]\big{)}.\end{array}

The equation of state is given by the Grüneisen model characterized by the following parameters: $\rho_{0}=1.845$ $kg/cm^{3}$ , $c_{0}=1.287$ $cm/{\mu s}$ , $\Gamma_{0}=1.11$ and $s=1.124$ . In addition, the shear modulus $G=1.1519$ $Mbar$ and the yield strength $\sigma_{0}=1$ $Mbar$ . The first flexural mode has a period of $2\pi/\omega=26.6266$ $\mu s$ . All the boundaries are set to have free boundary conditions. A set of uniformly refined meshes (i.e., $10\times 4$ , $20\times 8$ and $40\times 16$ ) are calculated. This case is run to $t=100$ $\mu s$ , demonstrating the robustness and accuracy of the DG(P3) method with cubic meshes. Dissipation errors will artificially dampen the oscillations so the goal is to maintain the peak amplitudes.

Fig. 11 presents the mesh and pressure contours at different times respectively. The legend scales are set differently for each time. The meshes and pressure contours at $t=6.65$ $\mu s$ (one quarter period ) and $t=33.25$ $\mu s$ (three quarters period) (the difference is one period approximately), look very similar qualitatively due to the low numerical dissipation of fourth-order accurate DG(P3) method. In addition, the vertical velocity of the central point for the plate is recorded for assessing the grid convergence that is shown in Fig. 12. There is no visible difference about the evolution of the vertical velocity profile on uniformly refined meshes, showing the accurate solution with DG(P3). Due to its low dissipation, DG(P3) also maintains a uniform amplitude for the vertical velocity. And it can be observed that the period based on the vertical velocity is approximately equal to the analytical one.

4.5 Taylor bar impact on a wall

This problem consists of a 2D aluminum bar impacting a rigid wall. The kinetic energy in the rod is entirely converted into internal energy through plastic dissipation. This type of test was initially proposed by Taylor [87] with a cylindrical rod. In this paper, we considered the 2D planar bar impact test that was simulated in [88]. The equation of state is given by the Grüneisen model characterized by the following parameters: $\rho_{0}=2.785$ $kg/cm^{3}$ , $c_{0}=0.533$ $cm/{\mu s}$ , $\Gamma_{0}=1.5111$ and $s=1.338$ . In addition, the shear modulus $G=0.276$ $Mbar$ and the yield strength $\sigma_{0}=0.003$ $Mbar$ . The computational domain is the initial projectile $(x,y)\in[0$ $cm,5$ $cm]\times[0$ $cm,1$ $cm]$ . A set of uniformly refined meshes (i.e., $10\times 4$ , $20\times 8$ and $40\times 16$ ) are used to perform a grid convergence study. This test case is run to $t=50$ $\mu s$ . The initial velocity is set to $\textit{v}=(-0.015$ $cm/{\mu s},0)$ . There exists no exact solution for this problem, and this case is used to show the robustness of the method and for accuracy bench-marking.

Fig. 13 presents the meshes and density fields. As shown, the mesh is moving in a curvilinear manner where DG(P2) results are shown in Fig. 13a and DG(P3) results are shown in 13b. The wave structure is captured in more details by the DG(P3) method compared with the DG(P2) method.

In this test problem, the kinetic energy is entirely converted into internal energy through the plastic dissipation shown in Fig. 14a and 14b. The temporal evolution of the length of the Taylor bar is also shown in Fig. 14c and 14d. Taking the solution on the finest mesh as a reference, the convergence study shows that the DG(P3) method on cubic meshes converges quickly to the reference solution (i.e., the evolution of energy balance and length). Furthermore, favorable results can be obtained using a very coarse mesh resolution (i.e., $10\times 4$ ).

5 Conclusions

The motivation for the work was to create a robust Lagrangian DG hydrodynamic method with a hierarchical orthogonal basis for material dynamics that is up to fourth-order accurate on cubic cells, i.e., each edge is defined by a cubic polynomial. The subcell mesh stabilization (SMS) method was extended to work with cubic cells. The SMS method is an important part of the new scheme to ensure robust solutions on problems with strong shocks. Another key part to ensuring robust solutions on strong shock problems is using hierarchical orthogonal basis functions for the polynomials of higher degrees of freedom (P3 or higher). We found that the limiting of the polynomials constructed from Taylor basis functions with the DG(P3) method gave very poor results on strong shock problems like the Sedov blast wave. The poor results were attributed to the moments in the mass matrix being coupled. This motivated us to create a new DG method based on a hierarchical orthogonal basis about the center of mass, resulting in a diagonal mass matrix (i.e., no matrix inversion is required) that are constant in time, which ensures accurate and robust solutions on strong shock problems with limiting. According to the operational cost of matrix multiplication, this formulation is also computationally more efficient compared with the Lagrangian DG method based on Taylor-basis functions.

A suite of test cases covering gas and solid dynamics were simulated. Likewise, vortical flows and shock driven flows were simulated to demonstrate the robustness of the proposed new high-order Lagrangian DG hydrodynamic method. The designed order of accuracy (up to fourth order) was demonstrated on the Taylor-Green vortex. The ability to simulate large deformation problems was demonstrated on the Taylor-Green and Gresho vortex problems using cubic cells. In addition, with Taylor-Green vortex, DG(Pn) can deliver the $(n+1)^{th}$ order of accuracy on quadrilaterals meshes with edges defined by a polynomial of degree $n$ or higher in Lagrangian hydrodynamics. The Sedov blast wave problem showed that the new method is robust and accurate on very strong shocks. The results from the Taylor anvil impact test and the vibrating beam problem demonstrate the method can accurately simulate solid dynamics. For the Taylor anvil test case, we also showed that the fourth-order accurate method was superior to the third-order accurate method.

This work fills an existing algorithmic gap by delivering a computationally efficient high-order hydrodynamic method to robustly and accurately evolve cubic meshes. Using cubic meshes is key to improving the fidelity of engineering and physics simulations, because a cubic cell can reproduce the as-built part dimensions in a simulation so that the exact volume, shape, mass, and density can be used. With linear meshes, the part contours do not reproduce the correct volume, so to match the mass, the part contours must be artificially shifted. Mesh refinement with linear meshes will only recover the exact part contours in the limit of a zero mesh size. As such, the impact of this work covers accurate simulations of large deformation flows and shock driven material dynamics to enabling the correct initial conditions to be used in simulations with manageable resolutions.

6 Acknowledgments

We gratefully acknowledge the support of the NNSA through the Laboratory Directed Research and Development (LDRD) program at Los Alamos National Laboratory. The Los Alamos unlimited release number is LA-UR-21-21718. Los Alamos National Laboratory is managed by Triad National Security, LLC for the U.S. Department of Energy’s NNSA .

Appendix A Study of the hierarchical basis functions $\boldsymbol{\phi}$

For this section, we aim at proving the material derivatives for the basis functions $\boldsymbol{\phi}$ are equal to 0, consisting of the following three steps. From the previous work [60],

\frac{d}{dt}(\rho j)=0.

(44)

\begin{split}\frac{d{\boldsymbol{\xi}}}{dt}=0,\frac{d{\boldsymbol{\xi}}_{cm}}{dt}=0,\end{split}

(45)

where, ${\boldsymbol{\xi}}=(\xi,\eta)$ and ${\boldsymbol{\xi}}_{cm}=(\xi_{cm},\eta_{cm})$ . Based on Eq. 44 and 45, it is possible to get

\frac{d\boldsymbol{\psi}}{dt}=0.

(46)

Let’s assume

A_{mn}=-\frac{\langle\psi_{m},\phi_{n}\rangle}{\langle\phi_{n},\phi_{n}\rangle},\quad\quad m=2,...,N(P),\quad 1\leq n<m.\\

(47)

For DG(P3),

\begin{array}[]{l}\begin{bmatrix}\phi_{1}\\ \phi_{2}\\ \phi_{3}\\ \phi_{4}\\ \phi_{5}\\ \phi_{6}\\ \phi_{7}\\ \phi_{8}\\ \phi_{9}\\ \phi_{10}\\ \end{bmatrix}=\begin{bmatrix}1&\\ A_{21}&1\\ A_{31}&A_{32}&1\\ A_{41}&A_{42}&A_{43}&1\\ A_{51}&A_{52}&A_{53}&A_{54}&1\\ A_{61}&A_{62}&A_{63}&A_{64}&A_{65}&1\\ A_{71}&A_{72}&A_{73}&A_{74}&A_{75}&A_{76}&1\\ A_{81}&A_{82}&A_{83}&A_{84}&A_{85}&A_{86}&A_{87}&1\\ A_{91}&A_{92}&A_{93}&A_{94}&A_{95}&A_{96}&A_{97}&A_{98}&1\\ A_{X1}&A_{X2}&A_{X3}&A_{X4}&A_{X5}&A_{X6}&A_{X7}&A_{X8}&A_{X9}&1\\ \end{bmatrix}\begin{bmatrix}\psi_{1}\\ \psi_{2}\\ \psi_{3}\\ \psi_{4}\\ \psi_{5}\\ \psi_{6}\\ \psi_{7}\\ \psi_{8}\\ \psi_{9}\\ \psi_{10}\\ \end{bmatrix}\end{array}

(48)

Therefore,

\frac{d\phi_{1}}{dt}=0.\\

(49)

Then it is possible to derive

\begin{split}\frac{d}{dt}\left({\int\rho\psi_{2}\phi_{1}dw}\right)=\int\frac{d}{dt}(\rho\psi_{2}\phi_{1}j)d\Omega=\int\Big{[}\frac{d}{dt}(\rho j)\psi_{2}\phi_{1}+\frac{d}{dt}(\psi_{2}\phi_{1})\rho j\Big{]}d\Omega=0,\\ \frac{d}{dt}\left({\int\rho\phi_{1}\phi_{1}dw}\right)=\int\frac{d}{dt}(\rho\phi_{1}\phi_{1}j)d\Omega=\int\Big{[}\frac{d}{dt}(\rho j)\phi_{1}\phi_{1}+\frac{d}{dt}(\phi_{1}\phi_{1})\rho j\Big{]}d\Omega=0.\\ \end{split}

(50)

Therefore,

\frac{dA_{21}}{dt}=\frac{d}{dt}\left(-\frac{\int\rho\psi_{2}\phi_{1}dw}{\int\rho\phi_{1}\phi_{1}dw}\right)=0.\\

(51)

As such,

\frac{d\phi_{2}}{dt}=\frac{d}{dt}(\psi_{2}+A_{21}\psi_{1})=0.\\

(52)

Likewise, it is also possible to derive,

\begin{split}&\frac{d}{dt}\left({\int\rho\psi_{3}\phi_{1}dw}\right)=\int\frac{d}{dt}(\rho\psi_{3}\phi_{1}j)d\Omega=\int\Big{[}\frac{d}{dt}(\rho j)\psi_{3}\phi_{1}+\frac{d}{dt}(\psi_{3}\phi_{1})\rho j\Big{]}d\Omega=0,\\ &\frac{d}{dt}\left({\int\rho\phi_{1}\phi_{1}dw}\right)=0.\\ \end{split}

(53)

and

\begin{split}&\frac{d}{dt}\left({\int\rho\psi_{3}\phi_{2}dw}\right)=\int\frac{d}{dt}(\rho\psi_{3}\phi_{2}j)d\Omega=\int\Big{[}\frac{d}{dt}(\rho j)\psi_{3}\phi_{2}+\frac{d}{dt}(\psi_{3}\phi_{2})\rho j\Big{]}d\Omega=0,\\ &\frac{d}{dt}\left({\int\rho\phi_{2}\phi_{2}dw}\right)=0.\\ \end{split}

(54)

Therefore,

\frac{dA_{31}}{dt}=\frac{d}{dt}\left(-\frac{\int\rho\psi_{3}\phi_{1}dw}{\int\rho\phi_{1}\phi_{1}dw}\right)=0,\\ \frac{dA_{32}}{dt}=\frac{d}{dt}\left(-\frac{\int\rho\psi_{3}\phi_{2}dw}{\int\rho\phi_{2}\phi_{2}dw}\right)=0.\\

(55)

Using a similar way, we can get,

\begin{split}\frac{d\boldsymbol{\phi}}{dt}=0,\\ \frac{dA_{mn}}{dt}=0.\\ \end{split}

(56)

It is easy to get the orthogonal basis is still hierarchical as long as the initial basis is ordered in a hierarchical manner.

{\mathbb{U}}_{h}({\boldsymbol{\xi}},t)=\underbrace{\underbrace{\underbrace{\underbrace{{\mathbb{U}}_{1}}_{P0}+{\mathbb{U}}_{2}\phi_{2}+{\mathbb{U}}_{3}\phi_{3}}_{P1}+{\mathbb{U}}_{4}\phi_{4}+{\mathbb{U}}_{5}\phi_{5}+{\mathbb{U}}_{6}\phi_{6}}_{P2}+{\mathbb{U}}_{7}\phi_{7}+{\mathbb{U}}_{8}\phi_{8}+{\mathbb{U}}_{9}\phi_{9}+{\mathbb{U}}_{10}\phi_{10}}_{P3}

(57)

But we need to mention that, although the hierarchy is the same as the original Taylor bases, the unknown solution is not the same as the original unknown individually. For example, ${\mathbb{U}}_{3}$ is the linear combination of $\frac{\partial\mathbb{U}}{\partial\xi}$ and $\frac{\partial\mathbb{U}}{\partial\eta}$ .

Appendix B Gauss-Lobatto quadrature rules

Some Gauss-Lobatto quadrature rules are tabulated for the interval $[-1,1]$ in Table B1.

Table B1: The Gauss-Lobatto point positions and weights for the Gauss-Lobatto quadrature rule.

\xi

and

wt

denote the position and weight of the quadrature points.

$G$ 1 2 3 4 5 6 7 2 points $wt$ 1.0 1.0 $\xi$ -1.0 1.0 3 points $wt$ $\frac{1}{3}$ $\frac{4}{3}$ $\frac{1}{3}$ $\xi$ -1.0 0 1.0 4 points $wt$ $\frac{1}{6}$ $\frac{5}{6}$ $\frac{5}{6}$ $\frac{1}{6}$ $\xi$ -1.0 $-\frac{\sqrt{5}}{5}$ $\frac{\sqrt{5}}{5}$ 1.0 5 points $wt$ $\frac{1}{10}$ $\frac{49}{90}$ $\frac{32}{45}$ $\frac{49}{90}$ $\frac{1}{10}$ $\xi$ -1.0 $\frac{-\sqrt{21}}{7}$ 0 $\frac{-\sqrt{21}}{7}$ 1.0 6 points $wt$ $\frac{1}{15}$ $\frac{\left(14-\sqrt{7}\right)}{30}$ $\frac{\left(14+\sqrt{7}\right)}{30}$ $\frac{\left(14+\sqrt{7}\right)}{30}$ $\frac{\left(14-\sqrt{7}\right)}{30}$ $\frac{1}{15}$ $\xi$ $-1$ $-\sqrt{\frac{1}{21}\left(7+2\sqrt{7}\right)}$ $-\sqrt{\frac{1}{21}\left(7-2\sqrt{7}\right)}$ $\sqrt{\frac{1}{21}\left(7-2\sqrt{7}\right)}$ $\sqrt{\frac{1}{21}\left(7+2\sqrt{7}\right)}$ 1 7 points $wt$ $\frac{1}{21}$ $\frac{124-7\sqrt{15}}{350}$ $\frac{124+7\sqrt{15}}{350}$ $\frac{256}{525}$ $\frac{124+7\sqrt{15}}{350}$ $\frac{124-7\sqrt{15}}{350}$ $\frac{1}{21}$ $\xi$ $-1$ $-\sqrt{\frac{5}{11}+\frac{2}{11}\sqrt{\frac{5}{3}}}$ $-\sqrt{\frac{5}{11}-\frac{2}{11}\sqrt{\frac{5}{3}}}$ 0 $\sqrt{\frac{5}{11}-\frac{2}{11}\sqrt{\frac{5}{3}}}$ $\sqrt{\frac{5}{11}+\frac{2}{11}\sqrt{\frac{5}{3}}}$ 1

References

[1] J von Neumann and R Richtmyer. A method for the calculation of hydrodynamics shocks. Journal of Applied Physics, 21:232–237, 1950.
[2] M. Wilkins. Use of artificial viscosity in multidimensional shock wave problems. Journal of Computational Physics, 36:281–303, 1980.
[3] D. Burton. Multidimensional discretization of conservation laws for unstructured polyhedral grids. Technical Report UCRL-JC-118306, Lawrence Livermore National Laboratory, 1994.
[4] E. Caramana, D. Burton, M. Shashkov, and P. Whalen. The construction of compatible hydrodynamic algorithms utilizing conservation of total energy. Journal of Applied Physics, 146:227–262, 1998.
[5] N. Morgan, K. Lipnikov, D. Burton, and M. Kenamond. A Lagrangian staggered grid Godunov-like approach for hydrodynamics. Journal of Computational Physics, 259:568–597, 2014.
[6] P-H. Maire, R. Loubère, and P. Vachal. Staggered Lagrangian discretization based on cell-centered Riemann solver and associated hydrodynamics scheme. Communications in Computational Physics, 10:940–978, 2011.
[7] S.K. Godunov, A. Zabrodine, M. Ivanov, A. Kraiko, and G. Prokopov. Résolution numréque des problèmes multidimensionnels de la dynamique des gaz. Mir, 1979.
[8] S. Godunov. Reminiscences about difference schemes. Journal of Computational Physics, 153:6–25, 1999.
[9] B. Després and C. Mazeran. Lagrangian gas dynamics in two dimensions and Lagrangian systems. Arch. Rational Mech. Anal., 178:327–372, 2005.
[10] F. Addessio, J. Baumgardner, J. Dukowicz, N. Johnson, B. Kashiwa, R. Rauenzahn, and C. Zemach. CAVEAT: A computer code for fluid dynamics problems with large distortion and internal slip. Technical Report LA-10613-MS-REV.1, Los Alamos National Laboratory, 1990.
[11] P-H. Maire, R. Abgrall, J. Breil, and J. Ovadia. A cell-centered Lagrangian scheme for two-dimensional compressible flow problems. SIAM Journal Scientific Computing, 29:1781–1824, 2007.
[12] P-H. Maire. A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured mesh. Journal Computational Physics, 228:2391–2425, 2009.
[13] D. Burton, T. Carney, N. Morgan, S. Sambasivan, and M. Shashkov. A cell centered Lagrangian Godunov-like method of solid dynamics. Computers $\&$ Fluids, 83:33–47, 2013.
[14] N. Morgan, M. Kenamond, D. Burton, T. Carney, and D. Ingraham. An approach for treating contact surfaces in Lagrangian cell-centered hydrodynamics. Journal of Computational Physics, 250:527–554, 2013.
[15] B. Cockburn and C-W Shu. The Runge-Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems. Journal of Computational Physics, 141:199–224, 1998.
[16] H. Luo, J. Baum, and R. Löhner. A discontinuous Galerkin method based on a Taylor basis for the compressible flows on arbitrary grids. Journal of Computational Physics, 227:8875–8893, 2008.
[17] H. Luo, J. Baum, and R. Löhner. A Hermite WENO-based limiter for discontinuous Galerkin method on unstructured grids. Journal of Computational Physics, 225:686–713, 2007.
[18] H. Luo, J. Baum, and R. Löhner. A p-multigrid discontinuous Galerkin method for the Euler equations on unstructured grids. Journal of Computational Physics, 211:767–783, 2006.
[19] Y. Xia, X. Liu, H. Luo, and R. Nourgaliev. A third-order implicit discontinuous Galerkin method based on a Hermite WENO reconstruction for time-accurate solution of the compressible Navier-Stokes equations. International Journal for Numerical Methods in Fluids, 79:416–435, 2015.
[20] X. Liu, L. Xuan, Y. Xia, and H. Luo. A reconstructed discontinuous Galerkin method for the compressible Navier-Stokes equations on three-dimensional hybrid grids. Computers $\&$ Fluids, 152:271–230, 2017.
[21] X. Liu, Y. Xia, H. Luo, and L. Xuan. A comparative study of Rosenbrock-type and implicit Runge-Kutta time integration for discontinuous Galerkin method for unsteady 3D compressible Navier-Stokes equations. Communications in Computational Physics, 20:1016–1044, 2016.
[22] C. Wang, H. Luo, and A. Kashi. A reconstructed discontinuous Galerkin methods for compressible flows in ALE formulation. 2018 AIAA Aerospace Sciences Meeting, AIAA-2018-0596, Kissimmee, Florida, 2018.
[23] A. Corrigan, A. Kercher, and D. Kessler. A moving discontinuous Galerkin finite element method for flows with interfaces. International Journal for Numerical Methods in Fluids, 89:362–406, 2019.
[24] A. Kercher, A. Corrigan, and D. Kessler. The moving discontinuous Galerkin finite element method with interface condition enforcement for compressible viscous flows. International Journal for Numerical Methods in Fluids, 2020. https://doi.org/10.1002/fld.4939.
[25] J. Cheng, X. Liu, T. Liu, and H. Luo. A parallel, high-order direct discontinuous Galerkin method for the Navier-Stokes equations on 3D hybrid grids. Communications in Computational Physics, 21:1231–1257, 2017.
[26] J. Cheng, Z. Du, X. Lei, Y. Wang, and J. Li. A two-stage fourth-order discontinuous Galerkin method based on the GRP solver for the compressible euler equations. Computers $\&$ Fluids, 181:248–258, 2019.
[27] W. Li, H. Luo, A. Pandare, and J. Bakosi. A p-adaptive discontinuous Galerkin method for compressible flows using Charm++. 2020 AIAA Aerospace Sciences Meeting, AIAA-2020-1565, Orlando, FL, 2020.
[28] J. Cheng, T. Liu, and H. Luo. A hybrid reconstructed discontinuous Galerkin method for compressible flows on arbitrary grids. Computers $\&$ Fluids, 139:68–79, 2016.
[29] X. Liu, Y. Xia, J. Cheng, and H. Luo. Development and assessment of a reconstructed discontinuous Galerkin method for the compressible turbulent flows on hybrid grids. 54st AIAA Aerospace Sciences Meeting, AIAA-2016-1359, San Diego, California, 2016.
[30] J. Cheng and T. Liu. A variational reconstructed discontinuous Galerkin method for the steady-state compressible flows on unstructured grids. Journal of Computational Physics, 380:65–87, 2019.
[31] J. Cheng, F. Zhang, and T. Liu. A high order compact least-squares reconstructed discontinuous Galerkin method for the steady-state compressible flows on hybrid grids. Journal of Computational Physics, 362:95–113, 2018.
[32] P.T. Greene, S.P. Schofield, and R. Nourgaliev. Dynamic mesh adaptation for front evolution using discontinuous Galerkin based weighted condition number relaxation. Journal of Computational Physics, 335:664–687, 2017.
[33] R. Nourgaliev, H. Luo, B. Weston, A. Anderson, S.Schofield, T. Dunn, and J.-P. Delplanque. Fully-implicit orthogonal reconstructed Discontinuous Galerkin method for fluid dynamics with phase change. Journal of Computational Physics, 305:964–996, 2016.
[34] F. Bassi, L. Botti, A. Colombo, D.A. Di Pietro, and P.Tesini. On the flexibility of agglomeration based physical space discontinuous Galerkin discretizations. Journal of Computational Physics, 231:45–65, 2012.
[35] J.S. Park and C. Kim. Higher-order multi-dimensional limiting strategy for discontinuous Galerkin methods in compressible inviscid and viscous flows. Computers $\&$ Fluids, 96:377–396, 2014.
[36] J. Lou, X. Liu, H. Luo, and H. Nishikawa. Reconstructed discontinuous Galerkin methods for hperbolic diffusion equations on unstructured grids. Communications in Computational Physics, 25:1302–1327, 2019.
[37] J. Lou, L. Li, X. Liu, H. Luo, and H. Nishikawa. Reconstructed discontinuous Galerkin methods based on first-order hyperbolic system for advection-diffusion equations. 23rd AIAA Computational Fluid Dynamics Conference, AIAA 2017-3445, Denver, CO, 2017.
[38] J. Lou, L. Li, X. Liu, H. Luo, and H. Nishikawa. First-order hyperbolic system based reconstructed discontinuous Galerkin methods for nonlinear diffusion equations on unstructured grids. 2018 AIAA Aerospace Sciences Meeting, AIAA 2018-2094, Kissimmee, FL, 2018.
[39] J. Lou, L. Li, H. Luo, and H. Nishikawa. Explicit hyperbolic reconstructed discontinuous Galerkin methods for time-dependent problems. 2018 Fluid Dynamics Conference, AIAA 2018-4270, Atlanta, GA, 2018.
[40] L. Li, X. Liu, J. Lou, H. Luo, H. Nishikawa, and Y. Ren. A compact high order finite volume method based on variational reconstruction for compressible flows on arbitrary grids. 23rd AIAA Computational Fluid Dynamics Conference, AIAA 2017-3097, Denver, CO, 2017.
[41] L. Li, X. Liu, J. Lou, H. Luo, H. Nishikawa, and Y. Ren. A discontinuous Galerkin method based on variational reconstruction for compressible flows on arbitrary grids. 2018 AIAA Aerospace Sciences Meeting, AIAA 2018-0831, Kissimmee, FL, 2018.
[42] L. Li, J. Lou, H. Luo, and H. Nishikawa. A new formulation of hyperbolic Navier-Stokes solver based on finite volume method on arbitrary grids. 2018 Fluid Dynamics Conference, AIAA 2018-4160, Atlanta, GA, 2018.
[43] L. Li, J. Lou, H. Luo, and H. Nishikawa. High-order Hyperbolic Navier-Stokes reconstructed discontinuous Galerkin method. AIAA Scitech 2019 Forum, AIAA 2019-1150, San Diego, CA, 2019.
[44] L. Li, J. Lou, H. Luo, and H. Nishikawa. High-order Hyperbolic Navier-Stokes reconstructed discontinuous Galerkin method. AIAA Aviation 2019 Forum, AIAA 2019-3060, Dallas, TX, 2019.
[45] L. Li, J. Lou, H. Luo, and H. Nishikawa. Reconstructed discontinuous Galerkin methods for compressible flows based on a new hyperbolic Navier-Stokes system. Journal of Computational Physics, 427:110058, 2021.
[46] Z. Jia and S. Zhang. A new high-order discontinuous Galerkin spectral finite element method for Lagrangian gas dynamics in two-dimensions. Journal of Computational Physics, 230:2496–2522, 2011.
[47] F. Vilar. Cell-centered discontinuous Galerkin discretization for two-dimensional Lagrangian hydrodynamics. Computers $\&$ Fluids, 64:64–73, 2012.
[48] F. Vilar, P-H. Maire, and R. Abgrall. A discontinuous Galerkin discretization for solving the two-dimensional gas dynamics equations written under total Lagrangian formulation on general unstructured grids. Journal of Computational Physics, 276:188–234, 2014.
[49] Z. Li, X. Yu, and Z. Jia. The cell-centered discontinuous Galerkin method for Lagrangian compressible Euler equations in two-dimensions. Computers $\&$ Fluids, 96:152–164, 2014.
[50] X. Liu, N. Morgan, and D. Burton. A Lagrangian discontinuous Galerkin hydrodynamic method. Computers $\&$ Fluids, 163:68–85, 2018.
[51] E. Lieberman, N. Morgan, D. Luscher, and D. Burton. A higher-order Lagrangian discontinuous Galerkin hydrodynamic method for elastic-plastic flows. Computers $\&$ Mathematics with Applications, 78:318–334, 2019.
[52] C. Wang, H. Luo, and M. Shashkov. A reconstructed discontinuous Galerkin method for compressible flows in Lagrangian formulation. Computers $\&$ Fluids, 202:104522, 2020.
[53] X. Liu, N. Morgan, and D. Burton. A Lagrangian cell-centered discontinuous Galerkin hydrodynamic method for 2D Cartesian and RZ axisymmetric coordinates. 2018 AIAA Aerospace Sciences Meeting, AIAA-2018-1562, Kissimmee, Florida, 2018.
[54] X. Liu, N. Morgan, and D. Burton. Lagrangian discontinuous Galerkin hydrodynamic methods in axisymmetric coordinates. Journal of Computational Physics, 373:253–283, 2018.
[55] J. Cheng and C-W Shu. A third order conservative Lagrangian type scheme on curvilinear meshes for the compressible euler equations. Communications in Computational Physics, 4:1008–1024, 2008.
[56] D. Kuzmin. A vertex-based hierarchical slope limiter for p-adaptive discontinuous Galerkin methods. Journal of Computational and Applied Mathematics, 233(12):3077–3085, 2010.
[57] D. Kuzmin. Hierarchical slope limiting in explicit and implicit discontinuous galerkin methods. Journal of Computational Physics, 257:1140–1162, 2014.
[58] T. Barth and D.C. Jespersen. The design and application of upwind schemes on unstructured meshes. 27th Aerospace Sciences Meeting, AIAA 1989-366, Reno, NV, 1989.
[59] N. Morgan, X. Liu, and D. Burton. Reducing spurious mesh motion in Lagrangian finite volume and discontinuous Galerkin hydrodynamic methods. Journal of Computational Physics, 372:35–61, 2018.
[60] X. Liu, N. Morgan, and D. Burton. A high-order Lagrangian discontinuous Galerkin hydrodynamic method for quadratic cells using a subcell mesh stabilization scheme. Journal of Computational Physics, 386:110–157, 2019.
[61] X. Liu and N. Morgan. Symmetry-preserving WENO-type reconstruction schemes in Lagrangian hydrodynamics. Computers $\&$ Fluids, 205:104528, 2020.
[62] X. Liu, N. Morgan, and D. Burton. A robust second-order accurate Lagrangian discontinuous Galerkin cell-centered hydrodynamic method on quadratic triangular cells.
[63] N. Morgan, X. Liu, and D. Burton. A Lagrangian discontinuous Galerkin hydrodynamic method for higher-order triangular elements. 2018 AIAA Aerospace Sciences Meeting, AIAA-2018-1092, Kissimmee, Florida, 2018.
[64] E. Lieberman, X. Liu, N. Morgan, D. Luscher, and D. Burton. A higher-order Lagrangian discontinuous Galerkin hydrodynamic method for solid dynamics. Computer Methods in Applied Mechanics and Engineering, 353:467–490, 2019.
[65] E. Lieberman, X. Liu, N. Morgan, and D. Burton. A multiphase Lagrangian discontinuous Galerkin hydrodynamic method for high-explosive detonation physics. Applications in Engineering Science, 4:100022, 2020.
[66] F. Bassi and S. Rebay. High-order accurate discontinuous finite element solution of the 2D Euler equations. Journal of Computational Physics, 138:251–285, 1997.
[67] D. Dubiner. Spectral methods on triangles and other domains. Journal of Scientific Computing, 6:345–390, 1991.
[68] Chi-Wang Shu. Total-Variation-Diminishing time discretizations. SIAM Journal on Scientific and Statistical Computing, 9:1073–1084, 1988.
[69] Chi-Wang Shu. Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics, 77:439–471, 1988.
[70] R.J. Spiteri and S.J. Ruuth. A new class of optimal high-order strong-stability preserving time discretization methods. SIAM Journal on Numerical Analysis, 40:469–491, 2002.
[71] N. Morgan. A dissipation model for staggered grid Lagrangian hydrodynamics. Computers $\&$ Fluids, 83:48–57, 2013.
[72] B. Cockburn, S. Hou, and C-W Shu. The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. IV. The multidimensional case. Mathematics of Computation, 54:545–581, 1990.
[73] C. Michalak and C. Ollivier-Gooch. Accuracy preserving limiter for the high-order accurate solution of the Euler equations. Journal of Computational Physics, 228:8693–8711, 2009.
[74] Per-Olof Persson and J. Peraire. Sub-cell shock capturing for discontinuous Galerkin methods. 44th AIAA Aerospace Sciences Meeting and Exhibit, AIAA 2006-112, Reno, Nevada, 2006.
[75] G. E. Barter and J. Peraire. Shock capturing with higher-order, PDE-based artificial viscosity. 18th AIAA Computational Fluid Dynamics Conference, AIAA 2007-3823, Miami, FL, 2007.
[76] C. Michalak and C. Ollivier-Gooch. A parameter-free generalized moment limiter for high-order methods on unstructured grids. Advances in Applied Mathematics and Mechanics, 1:451–480, 2009.
[77] L.Krivodonova, J. Xin, J.-F. Remacle, N. Chevaugeon, and J.E. Flaherty. Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws. Applied Numerical Mathematics, 48:323–338, 2004.
[78] V. Chiravalle and N. Morgan. A 3D finite element ALE method using an approximate Riemann solution. International Journal for Numerical Methods in Fluids, 83:642–663, 2016.
[79] P-H. Maire. A high-order one-step sub-cell force-based discretization for cell-centered Lagrangian hydrodynamics on polygonal grids. Computers $\&$ Fluids, 46:341–347, 2011.
[80] D. Burton, N. Morgan, T. Carney, and M. Kenamond. Reduction of dissipation in Lagrange cell-centered hydrodynamics CCH through corner gradient reconstruction CGR. Journal of Computational Physics, 299:229–280, 2015.
[81] V. Dobrev, T. Kolev, and R. Rieben. High-order curvilinear finite element methods for Lagrangian hydrodynamics. SIAM Journal on Scientific Computing, 34:B606–B641, 2012.
[82] N. Morgan, J. Waltz, D. Burton, M. Charest, T. Canfield, and J. Wohlbier. A Godunov-like point-centered essentially Lagrangian hydrodynamic approach. Journal of Computational Physics, 281:614–652, 2014.
[83] L. Sedov. Similarity and Dimensional Methods in Mechanics. Academic Press, 1959.
[84] C. Pederson, B. Brown, and N. Morgan. The Sedov blast wave as a radial piston verification test. Journal of Verification, Validation and Uncertainty Quantification, 1:1–9, 2016.
[85] D. Flanagan and T. Belytschko. A uniform strain hexahedron and quadrilateral with orthogonal hourglass control. International Journal for Numerical Methods in Engineering, 17:679–706, 1981.
[86] W.N. Weseloh, S.P. Clancy, and J.W. Painter. Pagosa physics manual. Technical Report LA-14425-M, Los Alamos National Laboratory, 2010.
[87] G. Taylor. The use of flat-ended projectiles for determining dynamic yield stress I. theoretical considerations. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, pages 289–299. 1948.
[88] P.-H. Maire, R. Abgrall, J. Breil, R. Loubère, and B. Rebourcet. A nominally second-order cell-centered Lagrangian scheme for simulating elasticplastic flows on two-dimensional unstructured grid. Journal of Computational Physics, 235:626–665, 2013.

A discontinuous Galerkin method based on a hierarchical orthogonal basis for Lagrangian hydrodynamics on curvilinear grids

Abstract

keywords:

1 Introduction

2 Governing equations

2.1 Constitutive models

3 Discretization

3.1 Basis functions for DG methods

3.2 Discontinuous Galerkin method

3.3 Temporal integration

3.4 Quadrature rules

3.4.1 Riemann solvers

3.5 Subcell Mesh Stabilization (SMS)

3.5.1 Surface integral

3.6 Limiting

4 Test problems

4.1 2D Taylor-Green vortex problem

4.2 2D Gresho vortex

4.3 Sedov blast problem

4.4 Elastic vibration of a Beryllium plate

4.5 Taylor bar impact on a wall

5 Conclusions

6 Acknowledgments

Appendix A Study of the hierarchical basis functions ϕ\boldsymbol{\phi}

Appendix B Gauss-Lobatto quadrature rules

References

Appendix A Study of the hierarchical basis functions $\boldsymbol{\phi}$