High-order variational Lagrangian schemes for compressible fluids

TL;DR

This paper introduces high-order variational Lagrangian finite element methods for compressible fluids that conserve energy and entropy, enabling larger stable time steps especially in challenging flow conditions.

Contribution

The paper develops a novel high-order variational Lagrangian scheme with energy conservation and entropy stability, using fully implicit time stepping for improved stability in compressible fluid simulations.

Findings

Scheme conserves mass, momentum, and energy.

Allows larger time steps due to implicit temporal discretization.

Demonstrates good performance through numerical results.

Abstract

We present high-order variational Lagrangian finite element methods for compressible fluids using a discrete energetic variational approach. Our spatial discretization is mass/momentum/energy conserving and entropy stable. Fully implicit time stepping is used for the temporal discretization, which allows for a much larger time step size for stability compared to explicit methods, especially for low-Mach number flows and/or on highly distorted meshes. Ample numerical results are presented to showcase the good performance of our proposed scheme.

Tables

Table 1. Table 1. Number of quadrature points N k subscript 𝑁 𝑘 N_{k} for the quadrature rule on a simplex that is accurate up to degree 2 ​ k + 1 2 𝑘 1 2k+1 for 0 ≤ k ≤ 6 0 𝑘 6 0\leq k\leq 6 .

	$k=0$	$k=1$	$k=2$	$k=3$	$k=4$	$k=5$	$k=6$
$N_k$ on Triangle	1	6	7	15	19	28	37
$N_k$ on Tetrahedron	1	8	14	36	61	109	171

Table 2. Table 2. Example 5.1 : History of convergence for the BDF [ m ] delimited-[] 𝑚 [m] - P m superscript 𝑃 𝑚 P^{m} scheme for 1 ≤ m ≤ 4 1 𝑚 4 1\leq m\leq 4 .

$m$	mesh	${\Vert {𝒖}_h-𝒖\Vert }_{\Omega }$	order	${\Vert e_h-e\Vert }_{\Omega }$	order
	$8\times 8$	1.052e-01	–	1.337e-01	–
	$16\times 16$	4.131e-02	1.35	7.284e-02	0.88
1	$32\times 32$	1.949e-02	1.08	3.697e-02	0.98
	$64\times 64$	9.710e-03	1.00	1.852e-02	1.00
	$8\times 8$	1.077e-02	–	1.263e-02	–
	$16\times 16$	3.250e-03	1.73	3.578e-03	1.82
2	$32\times 32$	7.933e-04	2.03	8.764e-04	2.03
	$64\times 64$	1.964e-04	2.01	2.192e-04	2.00
	$8\times 8$	5.590e-03	–	2.947e-03	–
	$16\times 16$	5.809e-04	3.27	3.175e-04	3.21
3	$32\times 32$	7.070e-05	3.04	4.181e-05	2.92
	$64\times 64$	8.766e-06	3.01	5.404e-06	2.95
	$8\times 8$	3.986e-03	–	2.510e-03	–
	$16\times 16$	1.361e-04	4.87	7.305e-05	5.10
4	$32\times 32$	5.690e-06	4.58	2.928e-06	4.64
	$64\times 64$	4.755e-07	3.58	3.304e-07	3.15

Equations

$$2\mathcal{D}={\int }_{\Omega }\eta |{\nabla }_s\text{boldmath}u{|}^2+(\xi -\frac{2}{3}\eta )|\nabla \cdot \text{boldmath}u{|}^2\mathrm{dx}$$

$$\psi (\rho ,\theta )=(c_p-c_v)\theta \rho \log (\rho )-c_v\rho \theta \log (\theta ),$$

$$p:=$$

$$e:=$$

$$s:=$$

$$\theta \partial s=\partial e+p\partial (\frac{1}{\rho }),$$

$$\frac{d}{dt}\text{boldmath}x(\text{boldmath}X,t)=\text{boldmath}u(\text{boldmath}x(\text{boldmath}X,t),t),$$

$$\frac{d}{dt}{\int }_{V^t}\rho \mathrm{dx}={\int }_{V^0}\frac{d}{dt}(\rho J)\mathrm{dX}=0,$$

$$\frac{d}{dt}(\rho J)=0,\text{ or }\rho (\text{boldmath}x(\text{boldmath}X,t),t)={\rho }_0(\text{boldmath}X)/J(\text{boldmath}X,t),$$

$$\mathcal{A}={\int }_0^T\left(\mathcal{K}-\mathcal{F}\right)\mathrm{dt}={\int }_0^T{\int }_{{\Omega }^t}\left(\frac{1}{2}\rho |\text{boldmath}u{|}^2-\psi (\rho ,\theta )\right)\mathrm{dx}\mathrm{dt}.$$

$$\text{boldmath}{f}_{\text{inertial}}=\frac{\delta {\int }_0^T\mathcal{K}\mathrm{dt}}{\delta \text{boldmath}x},\text{boldmath}{f}_{\text{cons}}=\frac{\delta {\int }_0^T\mathcal{F}\mathrm{dt}}{\delta \text{boldmath}x}.$$

$$\delta {\int }_0^T\mathcal{K}\mathrm{dt}=$$

$$=$$

$$\delta {\int }_0^T\mathcal{F}\mathrm{dt}=$$

$$=$$

$$\nabla p=\nabla ({\psi }_{\rho }\rho -\psi )=\rho \nabla {\psi }_{\rho }+\underset{=0}{\underbrace{{\psi }_{\rho }\nabla \rho -{\psi }_{\rho }\nabla \rho }}-{\psi }_{\theta }\nabla \theta .$$

$$\text{boldmath}{f}_{\text{diss}}=\frac{\delta \mathcal{D}}{\delta \text{boldmath}u},$$

$$\text{boldmath}{f}_{\text{diss}}=-\nabla \cdot \left(\eta {\nabla }_s\text{boldmath}u\cdot +(\xi -\frac{2}{3}\eta )(\nabla \cdot \text{boldmath}u)\text{boldmath}I\right),$$

$$\frac{\delta \mathcal{A}}{\delta \text{boldmath}x}=\frac{\delta \mathcal{D}}{\delta \text{boldmath}u},$$

$$\rho \dot{\text{boldmath}u}+\nabla p-\nabla \cdot \left(\eta {\nabla }_s\text{boldmath}u+(\xi -\frac{2}{3}\eta )(\nabla \cdot \text{boldmath}u)\text{boldmath}I\right)=0,$$

$${\int }_{{\Omega }^t}\left(\rho \dot{\text{boldmath}u}\cdot \text{boldmath}w-p\nabla \cdot \text{boldmath}w+\eta {\nabla }_s\text{boldmath}u\cdot {\nabla }_s\text{boldmath}w+(\xi -\frac{2}{3}\eta )(\nabla \cdot \text{boldmath}u)(\nabla \cdot \text{boldmath}w)\right)\mathrm{dx}=0,$$

$$\rho \dot{e}=\rho \theta \dot{s}-p\rho \dot{(1/\rho )}=\rho \theta \dot{s}-p\nabla \cdot \text{boldmath}u$$

$$\rho \dot{s}=\nabla \cdot \text{boldmath}j+\Delta .$$

$$\text{boldmath}j\theta =q,$$

$$q=\kappa \nabla \theta$$

$$0=\frac{d}{dt}{\int }_{\Omega }\left(\frac{1}{2}\rho |\text{boldmath}u{|}^2+\rho e\right)\mathrm{dx}=$$

$$0=$$

$$=$$

$$\Delta =\left(\eta |{\nabla }_s\text{boldmath}u{|}^2+(\xi -\frac{2}{3}\eta )(\nabla \cdot \text{boldmath}u{)}^2+\frac{\kappa |\nabla \theta {|}^2}{\theta }\right)/\theta ,$$

$$\rho \theta \dot{s}=\nabla \cdot (\kappa \nabla \theta )+\eta |{\nabla }_s\text{boldmath}u{|}^2+(\xi -\frac{2}{3}\eta )(\nabla \cdot \text{boldmath}u{)}^2.$$

$$\rho \dot{e}=-p\nabla \cdot \text{boldmath}u+\nabla \cdot (\kappa \nabla \theta )+\eta |{\nabla }_s\text{boldmath}u{|}^2+(\xi -\frac{2}{3}\eta )(\nabla \cdot \text{boldmath}u{)}^2.$$

Full text

G. Fu’s research is partially supported by NSF grant DMS-2012031. C. Liu’s research is partially supported by NSF grants DMS-1950868, DMS-2153029 and DMS-2118181.

1. Introduction

We are interested in the numerical simulation of multidimensional compressible fluid flows. There are two typical choices: a Lagrangian framework in which the computational mesh moves with the fluid velocity, and an Eulerian framework in which the fluid flows through a fixed spatial mesh. Lagrangian methods are widely used in many fields for multi-material flow simulations such as astrophysics and inertial confinement fusion (ICF), due to their distinguished advantage of being able to preserve the material interface sharply. In this article, we are concerned exclusively with Lagrangian methods.

To address problems in Lagrangian hydrodynamics, two fundamental approaches are commonly used. The first one is the staggered grid hydrodynamics (SGH) approach which employs a spatial discretization where the kinematic variables such as position and velocity are continuous and defined at the mesh vertices, and the thermodynamic variables including density, pressure, and internal energy are discontinuous and defined at cell centers; see [35, 36]. Artificial viscosity, as originally proposed in [35], is used to generate entropy production and suppress numerical oscillations across shocks. The second approach is the cell-centered hydrodynamics (CCH), which treats all hydrodynamic variables as cell centered and use approximate Riemann solvers to reconstruct a continuous velocity field; see [24, 9, 21, 6]. This process naturally introduces a sufficient level of dissipation at shock boundaries without the need of artificial viscosity.

While the original SGH scheme [35] does not preserve total energy conservation, it can be recovered using the notion of “corner masses” and “corner forces” as in the compatible hydro approach [5]. Energy-conservative high-order finite element schemes were further developed in [10] which can be viewed as a high-order generalization of the SGH approach since the key idea of using a continuous kinematic approximation space and a discontinuous thermodynamic approximation space were adopted. Advantages of the high-order Lagrangian schemes in [10] over the low-order schemes in [5] include the ability to more accurately capture the geometry of the flow and maintain robustness with respect to mesh motion using high-order curved meshes, sharper shock resolution and better symmetry preservation for symmetric flows.

Our spatial discretization follows the SGH approach using a (high-order) continuous finite element space for the kinematic variables, and a (high-order) discontinuous finite element space for the thermodynamic variables, which turns out to be closely related to the high-order finite element scheme [10]. However, we use a completely different derivation. While most of the existing Lagrangian schemes are obtained by directly discretizing the underling PDE system, the starting point of our spatial discretization is an energy dissipation law. In particular, we present a class of variational Lagrangian schemes for compressible flows using a discrete Energetic variational approach, which is an analogue to the energetic variational approach (EnVarA) [19, 14] in a semi-discrete level.

For a given energy-dissipation law and the kinematic (transport) relation, the EnVarA [19, 14] provides a general framework to determine the dynamics of system in a unique and well-defined way, through two distinct variational processes: Least Action Principle (LAP) and Maximum Dissipation Principle (MDP). This approach is originated from pioneering work of Onsager [26, 27] and Rayleigh [33], and has been successfully applied to build up many mathematical models [19, 34, 11, 14]. Most of existing EnVarA literature focuses on the isothermal case where temperature variation is not allowed. Here we adopt the approach used in [20] to model a thermodynamic system with temperature variation using EnVarA. The obtained model is then discretized in space using the discrete EnVarA, which leads to a high-order variational Lagrangian finite element scheme. Main structures of the continuous model including mass/momentum/energy conservation and entropy stability are naturally preserved in the proposed spatial discretization. The resulting ODE system is further discretized in time using high-order implicit time integrators. Due to the special structure of the ODE system, a nonlinear system only for the velocity degrees of freedom (DOFs) needs to be solved in each time step. The allowed time step size for stability is drastically improved over explicit time stepping, especially in the low Mach number regime, at the expense of a nonlinear system solve. Ample numerical examples are used to show the good performance of the proposed scheme.

We summarize the main features of our scheme:

• •

Space-time high-order accuracy.

• •

Mass/momentum/energy conservation and entropy stability for the spatial discretization.

• •

Implicit time stepping allows for large time step size especially for low Mach number flows and/or on highly distorted meshes.

The rest of the paper is organized as follows. In Section 2, we present the ideal gas model using EnVarA. In Section 3, a variational Lagrangian scheme is constructed using discrete EnVarA. The implicit temporal discretization of the resulting ODE system from Section 3 is then presented in Section 4. Numerical examples are presented in Section 5, followed by a summary in Section 6. In the Appendix, we briefly discuss our approach to the compressible isothermal case where temperature is fixed.

2. The energetic variational approach for ideal gas

The EnVarA [14] is a tool to derive the force balance equation starting from a total energy $E^{\text{total}}$ and an energy dissipation rate functional $\mathcal{D}$, where $E^{\text{total}}$ is the sum of the kinetic energy \mathcal{K}=\int_{\Omega}\frac{1}{2}\rho|\hbox{\boldmathu}|^{2}\mathrm{dx} and the Helmholtz free energy $\mathcal{F}=\int_{\Omega}\psi(\rho,\theta)\mathrm{dx}$, and

[TABLE]

is the rate of energy dissipation with dynamic viscosity $\eta$ and bulk viscosity $\xi$, which may depend on both $\rho$ and $\theta$. Here $\nabla_{s}$ denotes the symmetric gradient operator, $\rho$ is the density, $\theta$ is the absolute temperature, and $u$ is the fluid velocity.

Thermodynamics of idea gas is well studied in the literature [1, 2, 3, 7, 8, 12, 22, 29]. In classical thermodynamics, the concept of free energy proves to be useful [1, 29]. The internal energy and pressure of an ideal gas [22, 2] have a linear relationship with temperature and the product of temperature and density, respectively. Using this observation, Liu and Sulzbach [20] proposed a definition for the free energy density of an ideal gas:

[TABLE]

which we utilize in our current work. This definition includes the specific heat at constant volume $c_{v}$ and specific heat at constant pressure $c_{p}$. Associated with the free energy density (2), we define the three thermodynamic variables, namely pressure $p$, internal energy (per unit mass) $e$, and entropy (per unit mass) $s$:

[TABLE]

which are easily verified to satisfy the famous Gibbs equation (by chain rule) [8]:

[TABLE]

where $\partial$ represents any differentiation.

In Lagrangian coordinates, we introduce the flow map: \hbox{\boldmathx}(\hbox{\boldmathX},t):\Omega^{0}\rightarrow\Omega^{t}, which satisfies the trajectory equation

[TABLE]

with initial condition \hbox{\boldmathx}(\hbox{\boldmathX},0)=\hbox{\boldmathX}. Since we are concerning with conventional ideal gas, we postulate the kinematics of the temperature $\theta$ being transported along the trajectory: \frac{d}{dt}\theta=\theta_{t}+\hbox{\boldmathu}\cdot\nabla\theta.

2.1. Kinematics: mass conservation

Within an Lagrangian control volume V^{t}:=\{\hbox{\boldmathx}(\hbox{\boldmathX},t):\;\forall\hbox{\boldmathX}\in V^{0}\}, mass does not change over time:

[TABLE]

where J=\mathrm{Det}(\nabla_{X}\hbox{\boldmathx}) is the Jacobian determinant. Since equality (6) is valid for any control volume $V(t)$, there must hold

[TABLE]

where $\rho_{0}:\Omega^{0}\rightarrow\mathbb{R}^{+}$ is the initial density. The equality (7) represents strong mass conservation. Writing mass conservation (7) back to the Eulerian coordinates, we have \rho_{t}+\nabla\cdot(\rho\hbox{\boldmathu})=0, which is often referred to as the continuity equation. By abuse of notation, from it we can also get \delta\rho=-\nabla\cdot(\rho\delta\hbox{\boldmathx}), where $\delta$ represents the variational derivative. Such relation will be used and made clear in the derivations later in this paper.

2.2. Force balance: LAP and MDP

Here we combine the Least Action Principle (LAP) and Maximum Dissipation Principle (MDP) [26, 27, 33] to derive the force balance equation. The action functional for the system is

[TABLE]

The LAP performs variation on the kinetic and free energies to derive the inertial and conservative forces:

[TABLE]

Taking variation on the kinetic energy and using mass conservation (7), we get

[TABLE]

where we used the short-hand notation \dot{\hbox{\boldmathu}}:=\frac{d}{dt}\hbox{\boldmathu}(\hbox{\boldmathx}(\hbox{\boldmathX},t),t) for the material derivative. This implies that \hbox{\boldmathf}_{\text{inertial}}=-\rho\dot{\hbox{\boldmathu}}. Taking variation on the Helmholtz free energy $\mathcal{F}$ (Hamilton’s principle of virtual work) leads to

[TABLE]

where we used mass conservation and the kinematic transport of the temperature defined after formula (5) in the second row, and the definition of the pressure (3a) in the last row:

[TABLE]

This gives the conservative force \hbox{\boldmathf}_{\text{cons}}=\nabla p.

The MDP performs variation on the energy dissipation rate (1) to get the dissipative force:

[TABLE]

which implies

[TABLE]

where $I$ is the identity matrix.

Combining these, we get the force balance equation [14]

[TABLE]

which takes the following form

[TABLE]

This is the usual momentum equation for the Navier-Stokes equations, with its natural weak formulation

[TABLE]

for any test function $w$ with homogeneous boundary conditions.

2.3. Internal energy and entropy equations

The Gibbs equation (4) naturally gives an update equation for the internal energy:

[TABLE]

The rate of change of the entropy $s$ can be contributed by the entropy flux $j$ and entropy production $\Delta$:

[TABLE]

The second law of thermodynamics states that entropy production is non-negative: $\Delta\geq 0$. If we take the entropy flux $j$ being given by Durhem relation

[TABLE]

and the heat flux

[TABLE]

according to Fourier’s law in which $\kappa$ is the heat conductivity. From here, the explicit expression of $\Delta$ can be derived via conservation of total energy: there holds

[TABLE]

Using equations (12) with test function \hbox{\boldmathw}=\hbox{\boldmathu}, (13), and (14), we get:

[TABLE]

where we applied the chain rule for the heat flux term and used the homogeneous boundary condition \nabla\theta\cdot\hbox{\boldmathn}=0 on $\partial\Omega$. This implies that we can take the entropy dissipation rate as

[TABLE]

which satisfies the second law of thermodynamics as long as $\theta>0$. This implies the following entropy equation:

[TABLE]

Plugging (14) and (18) back to the internal energy equation (13), we obtain:

[TABLE]

By (3b), equation (20) equivalently gives the dynamics of temperature $\theta$, which has the heat equation as the leading term.

2.4. Summary

Combining the above results, we finally obtain the model equations:

[TABLE]

This is nothing but the compressible Navier-Stokes equations of an ideal gas [8]. Moreover, this system further satisfy the entropy equation (19). In the next section we derive a variational Lagrangian scheme for this system, where the spatial derivatives are evaluated by pulling back to the reference configuration (Lagrange coordinates); see, e.g., [36, 6, 10].

3. A discrete energetic variational approach: spatial discretization

In this section, we construct a variational Lagrangian scheme based on a discrete EnVarA. For simplicity, we ignore heat conduction in the model, i.e. we take $\kappa=0$ in this section.

3.1. Notation and the finite element spaces

Our grid-based scheme starts with a conforming triangulation $\mathcal{T}_{h}^{0}=\{T_{\ell}^{0}\}_{\ell=1}^{N_{T}}$ of the initial configuration $\Omega^{0}$ with $N_{T}$ elements, where we assume the element $T_{\ell}^{0}:=\Phi_{T_{\ell}}^{0}(\widehat{T})$ is obtained from a polynomial mapping $\Phi_{T_{\ell}}^{0}$ from the reference element $\widehat{T}$, which, for simplicity, is a simplex or a hypercube.

We denote $\mathcal{P}^{k}(\widehat{T})$ as the polynomial space of degree no greater than $k$ if $\widehat{T}$ is a reference simplex, or the tensor-product polynomial space of degree no greater than $k$ in each direction if $\widehat{T}$ is a reference hypercube, for $k\geq 1$. The mapped polynomial space on a spatial physical element $T^{0}\in\mathcal{T}_{h}^{0}$ is denoted as

[TABLE]

We denote $\{\widehat{\bm{\xi}}_{i}\}_{i=1}^{N_{k}}$ as a set of quadrature points with positive weights $\{\widehat{\omega}_{i}\}_{i=1}^{N_{k}}$ that is accurate for polynomials of degree up to $2k+1$ on the reference element $\widehat{T}$, i.e.,

[TABLE]

Note that when $\widehat{T}$ is a reference square, we simply use the Gauss-Legendre quadrature rule with $N_{k}=(k+1)^{2}$, which is optimal. On the other hand, when $\widehat{T}$ is a reference simplex, the optimal choice of quadrature rule is more complicated; see, e.g., [38, 37] and references cited therein. Table 1 list the number $N_{k}$ for $0\leq k\leq 6$ of the symmetric quadrature rules provided in [38].

The integration points and weights on a physical element $T_{\ell}^{0}$ are simply obtained via mapping: $\{\bm{\xi}_{i}^{\ell}:=\Phi_{T_{\ell}^{0}}(\widehat{\bm{\xi}}_{i})\}_{i=1}^{N_{k}}$, and $\{\omega_{i}^{\ell}:=|\nabla\Phi_{T_{\ell}^{0}}(\widehat{\bm{\xi}}_{i})|\widehat{\omega}_{i}\}_{i=1}^{N_{k}}$. To simplify the notation, we denote the set of physical integration points and weights

[TABLE]

and denote $(\cdot,\cdot)_{h}$ as the discrete inner-product on the mesh $\mathcal{T}_{h}^{0}$ using the quadrature points $\Xi_{h}^{k}$ and weights $\Omega_{h}^{k}$:

[TABLE]

We are now ready to present our continuous and discontinuous finite element spaces:

[TABLE]

where the local space

[TABLE]

is associated with the integration rule in (22) such that $\dim W^{k}(T_{\ell}^{0})=N_{k}$, and the nodal conditions

[TABLE]

in which $\delta_{ij}$ is the Kronecker delta function determines a unique solution $\varphi_{i}^{\ell}\in W^{k}(T_{\ell}^{0})$. This implies that $\{\varphi_{i}^{\ell}\}_{i=1}^{N_{k}}$ is a set of nodal bases for the space ${W}^{k}(T_{\ell}^{0})$, i.e.,

[TABLE]

When $T_{\ell}^{0}$ is a mapped hypercube, we have $N_{k}=(k+1)^{2}$, hence $\delta W_{k}(T_{\ell}^{0})=\emptyset$ and $W^{k}(T_{\ell})$ is simply the (mapped) tensor product polynomial space $\mathcal{P}^{k}(T_{\ell}^{0})$. On the other hand, when $T_{\ell}^{0}$ is a mapped simplex, we have $\dim\delta W_{k}(T_{\ell}^{0})=N_{k}-\dim\mathcal{P}^{k}(T_{\ell}^{0})>0$ for $k\geq 1$. However, we emphasize that the explicit expression of $\delta W_{k}(T_{\ell}^{0})$ or the basis function $\phi_{i}^{\ell}$ does not matter in our discretization, as only their nodal degrees of freedom (DOFs) on the quadrature nodes will enter into the numerical integration. Any function \alpha(\hbox{\boldmathX},t) in $W_{h}^{k}$ (for fixed $t$) can be expressed as

[TABLE]

where $\{{\sf a}_{i}^{\ell}(t)\}$ are the unknown coefficients. We refer to $W_{h}^{k}$ as the (discontinuous) integration rule space, which *only * contains the $N_{T}\times N_{k}$ quadrature points and weights (23), and is easy to implement in practice.

We further denote a set of basis functions for \hbox{\boldmathV}_{h}^{k} as \{\hbox{\boldmath\varphi}_{i}\}_{i=1}^{N_{V}^{k}}, where $N_{V}^{k}$ is the dimension of \hbox{\boldmathV}_{h}^{k}. We use the continuous space \hbox{\boldmathV}_{h}^{k} to approximate the flow map \hbox{\boldmathx}_{h} and velocity \hbox{\boldmathu}_{h}, and the discontinuous space $W_{h}^{k}$ to approximate the density $\rho_{h}$, pressure $p_{h}$, internal energy $e_{h}$, temperature $\theta_{h}$, and entropy $s_{h}$. More specifically, we have

[TABLE]

where ${\sf X}_{h}=[{\sf x}_{1},\cdots,{\sf x}_{N_{V}^{k}}]^{T}$, ${\sf U}_{h}=[{\sf u}_{1},\cdots,{\sf u}_{N_{V}^{k}}]^{T}$, ${\sf R}_{h}=[{\sf\rho}_{1}^{1},\cdots,{\sf\rho}_{N_{k}}^{N_{T}}]^{T}$, ${\sf P}_{h}=[{\sf p}_{1}^{1},\cdots,{\sf p}_{N_{k}}^{N_{T}}]^{T}$, ${\sf E}_{h}=[{\sf e}_{1}^{1},\cdots,{\sf e}_{N_{k}}^{N_{T}}]^{T}$, ${\sf\Theta}_{h}=[{\sf\theta}_{1}^{1},\cdots,{\sf\theta}_{N_{k}}^{N_{T}}]^{T}$, and ${\sf S}_{h}=[{\sf s}_{1}^{1},\cdots,{\sf s}_{N_{k}}^{N_{T}}]^{T}$ are the time dependent coefficient vectors for \hbox{\boldmathx}_{h}, \hbox{\boldmathu}_{h}, $\rho_{h}$, $p_{h}$, $e_{h}$, $\theta_{h}$, and $s_{h}$, respectively. Note that by the thermodynamic relations (3), there are only two independent thermodynamic variables. Here we take density $\rho$ and temperature $\theta$ as the independent variables. The other variables will be updated through the discrete formulas (31) below.

3.2. Trajectory equation and mass conservation

With the notation given in (28), the trajectory equation (5) simply implies that

[TABLE]

We require mass conservation to be satisfied pointwise at the quadrature nodes level, specifically, (7) implies that

[TABLE]

where the discrete Jacobian on the quadrature point \hbox{\boldmath\xi}_{i}^{\ell} is J_{h}(\hbox{\boldmath\xi}_{i}^{\ell},t):=\mathrm{Det}(\nabla_{X}\hbox{\boldmathx}_{h}(\hbox{\boldmath\xi}_{i}^{\ell},t)).

3.3. Thermodynamic relations

We require the relations in (3) for the thermodynamic variables be satisfied on the quadrature points level, which implies

[TABLE]

for all $1\leq i\leq N_{k}$ and $1\leq\ell\leq N_{T}$. It is easy to see that the (pointwise) Gibbs equation (4) is satisfied:

[TABLE]

By the density definition (30) and Jacobi’s formula, we have

[TABLE]

which implies that

[TABLE]

We note that the above pointwise relations also hold for the classical low-order SGH schemes [35, 36, 5] where the thermodynamic variables are approximated via piecewise constants, which, however, does not hold in general for the high order scheme [10] due to the use of a different high-order thermodynamic finite element space.

3.4. Discrete EnVarA and velocity equation

Instead of discretizing the force balance equation (21c), here we discretize the energy law and use the EnVarA to derive the discrete force balance equation directly. We denote the discrete action functional

[TABLE]

where the discrete kinetic energy $\mathcal{K}_{h}$ and the discrete Helmholtz free energy $\mathcal{F}_{h}$ are given as

[TABLE]

and

[TABLE]

Moreover, the discrete dissipation rate is given as

[TABLE]

The discrete force balance equation (10) is then

[TABLE]

Elementary calculation, using (29), (30) and Jacobi’s formula, yields that

[TABLE]

where \dot{\hbox{\boldmathu}}_{h}=\sum_{i=1}^{N_{V}^{k}}{\sf u}_{i}(t)^{\prime}\hbox{\boldmath\varphi}_{i}, and the pressure $p_{h}\in W_{h}^{k}$ satisfies

[TABLE]

according to (31a). We also have

[TABLE]

Plugging (36) and (37) back to (35), and using the definition of the pressure, we get the semi-discrete force balance equation:

[TABLE]

where $\sigma_{h}$ is the viscous stress defined as

[TABLE]

Here the evaluation of spatial derivative terms shall be pulled back to the initial configuration $\Omega^{0}$. In particular,

[TABLE]

where F_{h}=\nabla_{X}\hbox{\boldmathx}_{h} is the deformation tensor. Equation (38) provides an ODE system for the velocity coefficient vector ${\sf U}_{h}(t)$. Taking test function \hbox{\boldmath\varphi}_{j} as a constant vector, we immediately obtain global momentum conservation:

[TABLE]

3.5. Energy conservation and temperature equation

We use energy conservation to get an update equation for the temperature coefficient vector ${\sf\Theta}_{h}(t)$. In the absence of heat conduction ($\kappa=0$), the spatial discretization of the internal energy equation (20) leads to

[TABLE]

Equivalently, the equation (40) has the following pointwise form for the coefficient vector ${\sf E}_{h}$:

[TABLE]

for all $1\leq\ell\leq N_{T},\,1\leq j\leq N_{k}$. Plugging in the relations (31a) and (31b) back to (41), we get the following ODE system for the coefficent vector $\Theta_{h}(t)$:

[TABLE]

One key observation is that (42) is a linear ODE system for ${\sf\Theta}_{h}$. Moreover, taking \hbox{\boldmath\varphi}_{j}=\hbox{\boldmathu}_{h} in (38) and $\phi_{j}^{\ell}=1$ in (40) and adding, we obtain the total energy conservation:

[TABLE]

3.6. The entropy equation

Combining the Gibbs equation (32) with the internal energy equation (41) and simplifying, we obtain the ODE system satisfied by the entropy:

[TABLE]

for all $1\leq\ell\leq N_{T}$, $1\leq j\leq N_{k}$. This implies that

[TABLE]

where positivity of the right hand side, i.e., semi-discrete entropy stability, is guaranteed as long as the temperature $\Theta_{h}>0$. We remark that the entropy stability (44) is a direct consequence of our special choice of (nodal) thermodynamic finite element space (25) and (27).

3.7. Artificial viscosity

To make the scheme robust even in the case of zero physical viscosities with $\eta=\xi=0$, we add artificial viscosity [35] to the system so that shocks can be dissipated. Specifically, we add to the stress term (39) an artificial stress tensor $\sigma_{h}^{av}$ of the following form:

[TABLE]

in which, following [10], the artificial viscosity coefficient $\mu_{av}$ is:

[TABLE]

where $q_{1}$ and $q_{2}$ are linear and quadratic scaling coefficients, $c_{s}=\sqrt{\gamma p_{h}/\rho_{h}}$ is the speed of sound with $\gamma=c_{p}/c_{v}$ being the adiabatic constant, \Delta_{s_{1}}\hbox{\boldmathu}_{h}:={s_{1}\cdot\nabla\hbox{\boldmathu}_{h}\cdot s_{1}} is the directional measure of compression and $\ell_{s_{1}}=\ell_{0}|J_{h}s_{1}|$ is the directional length scale along the direction $s_{1}$, and the two linear switches are \phi_{0}=\frac{|\nabla\cdot\hbox{\boldmathu}_{h}|}{\|\nabla\hbox{\boldmathu}_{h}\|}, and \phi_{1}=\begin{cases}1,&\text{if }\Delta_{s}\hbox{\boldmathu}_{h}<0,\\ 0,&\text{if }\Delta_{s}\hbox{\boldmathu}_{h}\geq 0.\\ \end{cases} Here the direction $s_{1}$ is the unit eigenvector of the symmetric tensor \nabla_{s}\hbox{\boldmathu}_{h} with the smallest eigenvalue $\lambda_{1}$, i.e.,

[TABLE]

With this notation, we have \Delta_{s_{1}}\hbox{\boldmathu}_{h}=\lambda_{1}. Moreover $\ell_{0}=h_{0}/k$ is the mesh size of the initial domain divided by the polynomial degree $k$. We refer interested reader to [10] and references cited therein for more discussion on the choice of the artificial viscosity coefficient.

3.8. Summary

The final form of the semi-discrete scheme is summarized in Algorithm 1 below. This spatial discretization is high-order, mass/momentum/energy conserving, and entropy stable.

4. Temporal discretization

In this section, we focus on the discretization of the ODE system in Algorithm 1. We use fully implicit time discretizations so that the fully discrete scheme is robust for all mach numbers. We refer to [10, 30] for alternative conservative explicit schemes.

4.1. First order energy dissipative scheme

Using implicit Euler for the time derivative terms, we arrive at the following first order scheme: Given data {\hbox{\boldmathx}}_{h}^{n-1},{\hbox{\boldmathu}}_{h}^{n-1}\in\hbox{\boldmathV}_{h}^{k} and ${\theta}_{h}^{n-1}\in W_{h}^{k}$ at time $t^{n-1}$, and time step size $\delta t$, find solution {\hbox{\boldmathx}}_{h}^{n},{\hbox{\boldmathu}}_{h}^{n}\in\hbox{\boldmathV}_{h}^{k} and ${\theta}_{h}^{n}\in W_{h}^{k}$ such that

[TABLE]

where the stress

[TABLE]

in which the artificial viscosity coefficient $\mu_{av}^{n}$ is evaluated at time $t^{n}$, and the Jacobian determinant J_{h}^{n}=|\nabla_{X}\hbox{\boldmathx}_{h}^{n}|. The above system can be solved by first expressing \hbox{\boldmathx}_{h}^{n} and $\theta_{h}^{n}$ in terms of \hbox{\boldmathu}_{h}^{n} using (47a) and (47c):

[TABLE]

and then solve the nonlinear system for \hbox{\boldmathu}_{h}^{n} in (47b) using (48). We use Newton’s method to solve this nonlinear system for the velocity DOFs.

For the scheme (47), positivity of density on the quadrature points is guaranteed as long as the Jacobian $J_{h}^{n}>0$ on these quadrature points. And positivity of the temperature (hence positivity of pressure and internal energy) is satisfied as long as the denominator of right hand side of (48b) stays positive, i.e.,

[TABLE]

Moreover, strong mass conservation is satisfied due to (30), and global momentum conservation is satisfied by taking \hbox{\boldmath\varphi}_{j} to be a global constant in (47b). Finally, taking \hbox{\boldmath\varphi}_{j}=\hbox{\boldmathu}_{h}^{n} in (47b) and combining with (47c), we get

[TABLE]

which implies that

[TABLE]

Hence, the total energy is dissipated over time for the scheme (47).

4.2. Second order energy-conservative scheme

Energy conservation can be recovered from the first order scheme (47) by applying a time filter, which is the same as the midpoint rule; see [4]. This time stepping algorithm has two steps and is recorded in Algorithm 2 for reference.

We notice that the BE step implies

[TABLE]

By the extrapolation relations (49), we have

[TABLE]

Combining these equations, we get total energy conservation for this two-step method:

[TABLE]

4.3. High-order BDF schemes

It is natural to extend the implicit Euler scheme (47) to higher order by replacing the backward Euler time difference terms in (47) using higher-order backward difference formulas (BDFs). For completeness, we record these high-order schemes with uniform time stepping in Algoritm 3 below.

5. Numerical results

We present numerical results in this section using the open-source finite-element software NGSolve [31], https://ngsolve.org/. For all the simulation results, we consider invisid models with zero physical viscosities.

We use a (variable time step size) BDF2 time stepping with high-order spatial discretizations in Algorithm 3 for all the examples, except for Example 5.1 where higher order BDF time steppings are also used to verify the space/time high-order accuracy of the proposed methods. For problems with shocks, we take $q_{1}=0.5$ and $q_{2}=2$ as the default choice of the artificial viscosity parameters in (46) unless otherwise stated.

We take the time step size as

[TABLE]

where the length scale $h_{min}=h_{0}\alpha_{0}/k$ with $\alpha_{0}$ being the minimal singular value of the Jacobian matrix \nabla_{X}\hbox{\boldmathx}_{h}. Here the default choice of the CFL constant is taken to be $\text{CFL}=1$ unless otherwise stated. The automatic time-step control detailed in [10, 7.3] is also used for the examples with shocks. The Newton’s method is used to solve the nonlinear system for velocity DOFs in each time step, where the average iteration counts for all cases are observed to be around 4-8. Most of the linearized systems in each Newton iteration are solved using the sparse Cholesky factorization, with the only exception of the low-Mach number cases in Example 5.6 where a direct pardiso solver is used as the matrix failed to be positive definite therein.

5.1. Accuracy test: 2D Taylor-Green Vortex

We consider the invisid 2D Taylor-Green Vortex problem proposed in [10] to check the high-order convergence of our proposed algorithm on deforming domains. Following [10], we take $\gamma=5/3$ and add a source term

[TABLE]

to the internal energy equation (20). The computational domain is a unit square with wall boundary conditions, and the initial conditions are taken such that the exact solutions are:

[TABLE]

We turn off artificial viscosity in the numerical simulations, and perform mesh convergence studies for the $L^{2}$-errors of velocity and internal energy at final time $t=0.5$ on a sequence of four consecutive uniform rectangular meshes with size $2^{3+l}\times 2^{3+l}$ for $l=0,1,2,3$. We consider the BDF$[m]$ scheme Algorithm 3 with polynomials of degree $m$ used for the spatial discretization for $m=1,2,3,4$. The time step size is taken be $\delta t=0.05/2^{l}$ for $l=0,1,2,3$. The midpoint rule Algorithm 2 (with smaller time steps) is used to generate the starting values for high order BDF schemes. History of convergence for the two $L^{2}$ errors at $t=0.5$ are recorded in Table 2. We observe $m$-th order convergence for both variables for all $1\leq m\leq 4$. Hence, spatial and temporal high order convergence is achieved. This convergence rate is optimal for the BDF time stepping, but suboptimal by one order for the spatial discretization. We find the convergence behavior of the BDF$[m+1]$-$P^{m}$ scheme (not reported here for simplicity) is similar to BDF$[m]$-$P^{m}$ for $1\leq m\leq 3$, which suggests the one-order spatial convergence rate reduction is unavoidable for our scheme.

We remark that our spatial discretization can be made slightly more efficient by taking polynomial degree one order lower for the thermodynamic variables than the flow map approximation, while still maintaining a similar convergence behavior. However, since the major computational cost for our scheme is in the nonlinear system solver in (47b). Such efficiency gain is not significant, and will not be investigated further in this work.

5.2. 1D Shock Tube

We consider a simple 1D Riemann problem, the Sod shock tube on the domain $\Omega=[-5,5]$ with initial condition

[TABLE]

Here $\gamma=1.4$. We apply the BDF2 scheme with polynomial degree $k=4$ on an initial mesh with 20 uniform cells. The results at final time $t=2.0$ on all $20\times(4+1)=100$ quadrature points are shown in Figure 1, along with the deformed cells. The numerical approximation agrees with the exact Riemann solution quite well even on this coarse mesh, where the shock is resolved within 2 cells. As typical of Lagrangian schemes, the contact discontinuity is captured without dissipation. We also observe the “wall heating” phenomenon in the internal energy at the contact.

5.3. 2D Sedov explosion

The Sedov explosion [32] models the expanding wave by an intense explosion in a perfect gas. It is a standard problem to test the ability of codes to preserve the radial symmetry of shocks. The domain is a square $\Omega=[0,1.2]\times[0,1.2]$. The initial condition is set to have unit density and zero velocity, and also to have zero internal energy except at the left bottom corner cell $T_{0}$, where its is a bilinear function whos value is $\frac{0.2448\times 4}{\mathrm{area}(T_{0})}$ on the left/bottom corner vertex and zero on the other three vertices. So the total initial internal energy is $\int_{\Omega_{0}}\rho e\mathrm{dX}=0.2448.$ Symmetry boundary conditions are imposed on the left and bottom boundaries, while free boundary conditions are used on the top and right boundaries. The analytic solution at time $t=1$ gives a shock at radius $r=1$ with a peak density of $6$.

We apply the BDF2 scheme with polynomial degree $k=4$ on initial uniform rectangular meshes of size $N\times N$ with $N=16$ and $N=32$. The density field on deformed meshes, and also the scattered plot of density versus radius $r=\sqrt{x^{2}+y^{2}}$ on all quadrature points are shown in Figure 2. The radial symmetry of the solution is preserved, and the numerical solution agrees quite well with the analytic solution in Figure 2(c).

5.4. 2D Noh explosion

The Noh explosion problem [25] consists of an ideal gas with $\gamma=5/3$, initial density $\rho_{0}=1$, initial internal energy $e_{0}=0$, and initial velocity \hbox{\boldmathu}_{0}=(-\frac{x}{\sqrt{x^{2}+y^{2}}},-\frac{y}{\sqrt{x^{2}+y^{2}}}). The computational domain is $\Omega=[0,1]\times[0,1]$. We use symmetry boundary conditions on left and bottom boundaries, and free boundary conditions on top and right boundaries. Similar to the previous case, this problem has a radial symmetry, and the analytic solution at time $t=0.6$ gives a shock at radius $r=0.2$ with a peak density of $16$.

We apply the BDF2 scheme with polynomial degree $k=4$ on initial uniform rectangular meshes of size $N\times N$ with $N=16$ and $N=32$. For this problem, we observe the default choice of artificial viscosity coefficients with $q_{1}=0.5$ and $q_{2}=2$ leads to quite large post-shock oscillations. So we increase these coefficients to $q_{1}=1$ and $q_{2}=4$ in our numerical experiments reported here. The radial symmetry of the solution is preserved as seen in (a) and (b) of Figure 3, and the numerical solution has a good agreement with the analytic solution for $0.1<r<0.3$ in Figure 2(c), although it is slightly oscillatory.

5.5. Triple point problem

The triple point problem is a multimaterial test case proposed in [17]; see also [13]. The initial data is shown in Fig. 4. The computational domain has a rectangular shape with $7\times 3$ edge ratio. It includes three materials at rest located in $\Omega_{1}=[0,1]\times[0,3]$, $\Omega_{2}=[1,7]\times[0,1.5]$, and $\Omega_{3}=[1,7]\times[1.5,3]$, initially forming a T-junction. The high-pressure material in $\Omega_{1}$ creates a shock wave moving to the right. Due to different material properties in $\Omega_{2}$ and $\Omega_{3}$, the shock wave moves faster in $\Omega_{3}$, which leads to vortex formation around the triple point. For pure Lagrangian methods, there is a limit to how long this problem can be run due to vortex generation. Similar to [10], we run simulation till time $t=3.3$.

We apply the BDF2 scheme with polynomial degree $k=4$ and on initial uniform rectangular meshes of size $28\times 12$ and $56\times 24$. Here we reduce the artificial viscosity coefficients to be $q_{1}=0.25$ and $q_{2}=1$, and take a larger CFL number with $\text{CFL}=3$. A total of 238 time steps is used to drive the solution to final time $t=3.3$ on the fine mesh with $56\times 28$ cells. We note that an explicit scheme would require about two order of magnitude more time steps for a stable simulation. For example, the Laghos code freely available in the github repository https://github.com/CEED/Laghos, which implements the high-order Lagrangian finite element scheme in [10], requires about $90,000$ time steps for $k=4$ on the fine mesh with RK4 time stepping and CFL=1. We further note that the default choice of parameters also work for this problem; we make the modifications to illustrate that the high-order method is still robust with a larger time step size and a smaller artificial viscosity.

The density plot in log-scale on the deformed meshes are shown in Figure 5. We observe the material interfaces are sharply preserved as typical of Lagrangian schemes, and the the shock locations are essentially the same, but the total amount of ”roll-up” at the triple point increases as the mesh refines.

5.6. Gresho vortex

The Gresho vortex problem [15] is an example of a stationary, incompressible rotating flow around the origin in two spatial dimensions, where centrifugal forces are exactly balanced by pressure gradients. It was first applied to the compressible Euler equations in [18]. Here we use the low-Mach setup given in [23]. The domain is a unit square $\Omega=(0,1)\times(0,1)$ with wall boundary conditions, $\gamma=1.4$, and the initial conditions are

[TABLE]

where the radius $r=\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}$, angular velocity $u_{\phi}$ is

[TABLE]

the unit vector \hbox{\boldmathe}_{\phi}=(-(y-0.5)/r,(x-0.5)/r), and $M_{\max}$ is the parameter used to adjust the maximum Mach number of the problem, where M(r)=\frac{|\hbox{\boldmathu}|}{\sqrt{\gamma p/\rho}} is the Mach number. The maximum number of the initial condition is achieved at $r=0.2$ with $M(0.2)=M_{\max}$.

We apply the BDF2 scheme with polynomial degree $k=2$ on a $40\times 40$ mesh, and $k=4$ on a $20\times 20$ mesh. The total number of velocity DOFs are the same for the two cases. Since the problem is smooth, we turn off the artificial viscosity. The period of one rotation for $r=0.2$ is $2\pi r=0.4\pi$. We run simulation till final time $t=\frac{3}{4}\times{0.4\pi}$ so that the internal flow has rotated $\frac{3}{4}\times 180=135$ degrees. The default choice of time step size (51) is linear proportional to the Mach number since $\max c_{s}\approx 1/M_{\max}$. We remove this Mach-number dependency on time step size by multiplying the sound speed in (51) with the maximum Mach number, i.e.,

[TABLE]

We use $M_{\max}=0.1,0.01,0.001$, and take the CFL number to be CFL=0.25 for all cases. The relative Mach number $M/M_{\max}$ at final time $t=\frac{3}{4}\times 0.4\pi$ for all cases are shown in Figure 6. The first row of Figure 6 show results for the $k=2$ simulations, where we clearly observe a locking phenomena as the Mach number decreases. On the other hand, the higher order simulations with $k=4$ leads to almost identical results for all three Mach numbers. This example illustrates the advantage of using a higher order scheme over a low-order scheme in the low-Mach number regime.

Moreover, the total number of time steps for $k=4$ are between 900 to 1000 for all three cases. If an explicit scheme (e.g., RK3) were to be used to solve this problem, the total number of time steps will be about three orders of magnitude larger when $M_{\max}=0.001$ due to sound speed based CFL constraints. Finally, we note that the nonlinear system in each time step is harder to solve as the Mach number $M_{\max}$ decreases. For example, Newton’s method with sparse Cholesky direct solver works for $M_{\max}=0.1$ with CFL=1, but it fails for $M_{\max}=0.001$, where we have to reduce CFL to be 0.25 and replace the Cholesky solver by a pardiso solver, which indicates the linear system for $M_{\max}=0.001$ in the Newton iteration is no longer positive definite. The linear system solver issue in the low Mach number regime will be further investigate in our future work.

5.7. Shock-bubble interaction

This test case corresponds to the interaction of shock wave with a cylindrical Helium bubble surrounded by air at rest [28]. We use the same steup as in [13, Section 8.4]. The initial domain is a rectangular box $\Omega_{0}=(0,L)\times(-H/2,H/2)=(0,0.650)\times(-0.089,0.089)$, which includes a circular bubble with center $(0.320,0)$ and radius $r_{b}=0.025$. Initial data are shown in Figure 7(a). Wall boundary conditions at each boundary is prescribed except at the right boundary, where we impose a piston-like boundary condition with inward velocity $(-124.824,0)$. The left going shock wave hits the bubble at time $t_{i}=668.153\times 10^{-6}$. The final time of simulation is $t_{f}=t_{i}+674\times 10^{-6}=1342.153\times 10^{-6}$, which corresponds to the time where experimental shadow-graph extracted from [16] is displayed in [28].

We use two unstructured triangular meshes that is fitted to the bubble boundary for this problem. The coarse mesh has 8324 triangular cells whose mesh size is $h=H/32$, while the fine mesh has 33526 cells whose mesh size is $h=H/64$; see Figure 7(b)-(c) for the zoom-in view of the two meshes around the bubble. The BDF2 time stepping is used in combination with polynomial degree $k=2$ and $k=4$ on these two meshes. The zoom-in view around the deformed bubble at final time $t_{f}$ are shown in Figure 8. We observe that the location and shape of the deformed bubble is similar for each simulation, with a better resolution being obtained on a finer mesh with a higher order polynomial degree. These shapes are also qualitatively similar to the experimental Schlieren image in Figure 8(e) obtained from [16].

We display in Figure 9 the time evolution of the bubble at times $t=800\times 10^{-6},1100\times 10^{-6},1342.153\times 10^{-6}$ for $k=4$ on the two meshes. We note that the results obtained with both meshes are quite similar.

5.8. Multimaterial implosion in cylindrical geometry

The aim of this example is to access the capability of the implicit high-order Lagrangian scheme to handle a multi-mode implosion in cylindrical geometry. Here we consider a simple 1D multimaterial implosion problem on unstructured 2D triangular meshes. The problem consists of a low-density material with $\rho_{1}=0.05$ in the radial range $r\in[0,1]$ surrounded by a shell of high density material with $\rho_{2}=1.0$ in the radial range $r\in(1.0,1.2]$. Both materials are initially at rest with pressure $p=0.1$ and adiabatic index $\gamma=5/3$. This problem was originally proposed in [13] with a time dependent pressure source on the outer radial surface $r=1.2$. Here we use the modification in [10] that applies a constant radial velocity source of \hbox{\boldmathu}=-5(x,y)/1.2 on the outer boundary, which drives a cylindrical shock wave inwards. Due to the 1D setup (symmetry) , the material interface should be a function of radius only for all time. The movement of the interface radius against time is shown on the left panel of Figure 11. We clearly observe the deceleration of the interface starting around $t=0.12$, and the so-called stagnation phase is reached around $t=0.14$ where the radius obtains its minimum value. The flow becomes Rayleigh-Taylor unstable after this time as a small perturbation of the interface will grow exponentially as function of time due to the fact that the light fluid inside is pushing the heavy fluid outside after the stagnation phase. It is very challenging to preserve the interface symmetry for a Lagrangian scheme on general unstructured meshes, especially when time past the stagnation phase.

Due to symmetry of the problem, we take the computational domain to be a quarter circle $\Omega=\{(x,y):x\geq 0,y\geq 0,x^{2}+y^{2}\leq 1.2^{2}\}$ with wall boundary conditions on the left and bottom boundaries. We apply the BDF2 scheme with polynomial degree $k=2,4$ on two set of unstructured triangular meshes. The coarse mesh has mesh size $h=0.05$ and is used for $k=4$, while the fine mesh is a uniform refinement of the coarse mesh and is used for $k=2$; see Figure 10(d) for the coarse mesh and Figure 10(a) for the fine mesh. We take $q_{1}=1$ and $q_{2}=4$ as the artificial viscosity coefficients and set CFL number equals 0.25. The two simulations have the same number of velocity DOFs, and requires a similar amount of total time steps to reach the final time $t=0.16$. Density plots in log-scale on the deformed meshes are shown in Figure 10 at times $t=0.08$ and $t=0.16$, along with the initial density at $t=0$. We observe while the results for $t=0.08$ are similar for both cases, where radial symmetry of the material interface is preserved quite well. The results for $k=4$ at $t=0.16$ (which past the stagnation phase) is better than that for $k=2$ in terms of the interface symmetry. In Figure 11 we plot the time evolution of the average radius of the material interface, and the normalized standard deviation of this radius at different times as an indication of the symmetry error over time. We observe similar results for the average radius for both simulations, and a smaller symmetry error for the higher order case.

6. Conclusion

We presented a class of high-order variational Lagrangian schemes for compressible flow. The discrete EnVarA approach is used to derive the high-order spatial finite element discretization. Features of our spatial discretization include mass/momentum/energy conservation and entropy stability. Fully implicit time stepping is then applied to the resulting ODE system. Each time step requires a nonlinear system solve for the velocity DOFs only. Ample numerical results are shown to support the good performance of the proposed scheme. We plan to extend this pure Lagrangian scheme to the arbitrary Eulerian Lagrangian framework, which has the potential to address the mesh distortion issue.

Appendix

Here we briefly discuss the isothermal case [14] where the temperature does not change over time. In this case, the free energy (2) is a function of density only: $\psi=\psi(\rho)$. Typical choices include $\psi(\rho)=\alpha\rho^{\gamma}$ or $\psi(\rho)=\alpha\rho\log(\rho)$.

This model only has two thermodynamic variables: density $\rho$ and pressure

[TABLE]

The EnVarA derivation leads to the following model equations; see also [14]:

[TABLE]

where $p=\psi_{\rho}\rho-\psi$.

The spatial discretization for the model (52) is summarized below:

It is clear that the spatial discretization in Algorithm 4 is mass and momentum conservative. Next we prove Algorithm 4 is also entropy stable. Using a similar derivation as in (9) and (36) and the definition of density and pressure, we obtain

[TABLE]

Taking test function \hbox{\boldmath\varphi}_{h}=\hbox{\boldmathu}_{h} in (53b) and use the above relation, we get

[TABLE]

The ODE system in Algorithm 4 can be discretized using the implicit schemes discussed in Section 4. Here we present a variational implicit time discretization similar to the backward Euler scheme in Section 4.1. Given data \hbox{\boldmathx}^{n-1},\hbox{\boldmathu}_{h}^{n-1}\in V_{h}^{k} at time $t^{n-1}$ and time step size $\delta t$, denote the following discrete energy functional for the flow map \hbox{\boldmathx}_{h}:

[TABLE]

where the dissipation

[TABLE]

with artificial viscosity $\mu_{av}^{n-1}$ explicitly evaluated at time level $t^{n-1}$, and J_{h}=|\nabla_{X}\hbox{\boldmathx}_{h}|. The flow map at next time level is obtained by solving the following minimization problem:

[TABLE]

Assuming piecewise constant viscosity coefficients $\eta,\xi$, and $\mu_{av}^{n-1}$, the Euler-Lagrangian equation for this minimization problem is simply the Backward Euler scheme in (47) applied to the ODE system (53). We note that such energy minimization interpretation is not available for the non-isothermal case discussed in Section 4 due to the temperature equation (47c).