Entropy stable spectral collocation schemes for the 3-D Navier-Stokes equations on dynamic unstructured grids

TL;DR

This paper introduces high-order entropy stable spectral collocation schemes for 3-D Navier-Stokes equations on dynamic unstructured grids, ensuring conservation, accuracy, and stability on moving and deforming meshes.

Contribution

It develops a novel entropy conservative flux and a scheme that preserves geometric conservation laws on dynamic unstructured grids.

Findings

Achieves the desired order of accuracy in numerical tests.

Ensures entropy stability on moving and deforming grids.

Maintains freestream preservation properties.

Abstract

New entropy stable spectral collocations schemes of arbitrary order of accuracy are developed for the unsteady 3-D Euler and Navier-Stokes equations on dynamic unstructured grids. To take into account the grid motion and deformation, we use an arbitrary Lagrangian-Eulerian (ALE) formulation. As a result, moving and deforming hexahedral grid elements are individually mapped onto a cube in the fixed reference system of coordinates. The proposed scheme is constructed by using the skew-symmetric form of the Navier-Stokes equations, which are discretized by using summation-by-parts spectral collocation operators that preserve the conservation properties of the original governing equations. Furthermore, the metric coefficients are approximated such that the geometric conservation laws (GCL) are satisfied exactly on both static and dynamic grids. To make the scheme entropy stable, a new entropy conservative flux is derived for the 3-D Euler and Navier-Stokes equations on dynamic unstructured grids. The new flux preserves the design order of accuracy of the original spectral collocation scheme and guarantees the entropy conservation on moving and deforming grids. We present numerical results demonstrating design order of accuracy and freestream preservation properties of the new schemes for both the Euler and Navier-Stokes equations on moving and deforming unstructured grids.

Tables

Table 1. Table 1: Face numbering convention

index	face number
$i=1$	1
$i=N$	2
$j=1$	3
$j=N$	4
$k=1$	5
$k=N$	6

Table 2. Table 2: Error convergence on dynamic grids for the isentropic vortex problem.

$p=3$	$L_2$ error	$L_2$ rate	$L_{\infty }$ error	$L_{\infty }$ rate
$6\times 6$	$3.92\times {10}^{-5}$	–	$4.26\times {10}^{-4}$	–
$12\times 12$	$6.50\times {10}^{-6}$	2.59	$9.99\times {10}^{-5}$	2.09
$24\times 24$	$7.85\times {10}^{-7}$	3.05	$1.55\times {10}^{-5}$	2.69
$48\times 48$	$7.85\times {10}^{-8}$	3.32	$1.58\times {10}^{-6}$	3.29
$p=4$
$6\times 6$	$1.71\times {10}^{-6}$	–	$2.63\times {10}^{-5}$	–
$12\times 12$	$1.88\times {10}^{-7}$	3.18	$3.27\times {10}^{-6}$	3.00
$24\times 24$	$1.26\times {10}^{-8}$	3.89	$3.39\times {10}^{-7}$	3.27
$48\times 48$	$7.00\times {10}^{-10}$	4.17	$2.23\times {10}^{-8}$	3.92
$p=5$
$6\times 6$	$6.69\times {10}^{-8}$	–	$1.00\times {10}^{-6}$	–
$12\times 12$	$3.86\times {10}^{-9}$	4.12	$9.74\times {10}^{-8}$	3.36
$24\times 24$	$1.26\times {10}^{-10}$	4.93	$3.34\times {10}^{-9}$	4.82
$48\times 48$	$3.19\times {10}^{-12}$	5.31	$7.73\times {10}^{-11}$	5.48

Table 3. Table 3: Error convergence on dynamic randomly perturbed grids for the viscous shock problem.

$p=3$	$L_2$ error	$L_2$ rate	$L_{\infty }$ error	$L_{\infty }$ rate
$6\times 6\times 6$	$1.99\times {10}^{-3}$	–	$8.33\times {10}^{-2}$	–
$12\times 12\times 12$	$1.95\times {10}^{-4}$	3.36	$1.37\times {10}^{-2}$	2.60
$24\times 24\times 24$	$1.85\times {10}^{-5}$	3.39	$1.22\times {10}^{-3}$	3.49
$48\times 48\times 48$	$2.56\times {10}^{-6}$	2.86	$2.12\times {10}^{-4}$	2.52
$p=4$
$6\times 6\times 6$	$4.05\times {10}^{-4}$	–	$1.52\times {10}^{-2}$	–
$12\times 12\times 12$	$2.55\times {10}^{-5}$	3.99	$2.09\times {10}^{-3}$	2.86
$24\times 24\times 24$	$1.19\times {10}^{-6}$	4.43	$1.35\times {10}^{-4}$	3.95
$48\times 48\times 48$	$7.04\times {10}^{-8}$	4.07	$9.71\times {10}^{-6}$	3.80
$p=5$
$6\times 6\times 6$	$9.62\times {10}^{-5}$	–	$5.36\times {10}^{-3}$	–
$12\times 12\times 12$	$3.53\times {10}^{-6}$	4.77	$6.53\times {10}^{-4}$	3.04
$24\times 24\times 24$	$7.42\times {10}^{-8}$	5.57	$1.24\times {10}^{-5}$	5.72
$48\times 48\times 48$	$2.28\times {10}^{-9}$	5.02	$6.19\times {10}^{-7}$	4.32

Equations

$${\mathbf{\xi }}_1={\left[{\xi }_{1,1},\dots ,{\xi }_{1,N_1}\right]}^{\mathrm{T}},$$

$$U\left({\xi }_1,{\xi }_2,{\xi }_3,t\right)\in L^2\left(\left[{\alpha }_1,{\beta }_1\right]\times \left[{\alpha }_2,{\beta }_2\right]\times \left[{\alpha }_3,{\beta }_3\right]\times \left[0,T\right]\right)$$

$$\mathbf{u}={\left[U\left({\xi }_{1,1},{\xi }_{2,1},{\xi }_{3,1},t\right),\dots ,U\left({\xi }_{1,N_1},{\xi }_{2,N_2},{\xi }_{3,N_3},t\right)\right]}^{\mathrm{T}}.$$

$$\mathsf{A}\otimes\mathsf{B}\equiv\left[\begin{array}[]{ccc}a_{1,1}\mathsf{B}&\dots&a_{1,m}\mathsf{B}\\ \vdots&&\vdots\\ a_{n,1}\mathsf{B}&\dots&a_{n,m}\mathsf{B}\end{array}\right].$$

$$\mathsf{A}\circ\mathsf{B}\equiv\left[\begin{array}[]{ccc}a_{11}b_{11}&\dots&a_{1n}b_{1n}\\ \vdots&&\vdots\\ a_{n1}b_{n1}&\dots&a_{nn}b_{nn}\end{array}\right].$$

$${\mathbf{\xi }}_l^k\equiv {\left[{\xi }_{l,1}^k,\dots ,{\xi }_{l,N_l}^k\right]}^{\mathrm{T}},$$

$${\left({\mathbf{\xi }}_l^i\right)}^{\mathrm{T}}{\mathsf{E}}_{{\xi }_l}^{(1D)}{\mathbf{\xi }}_l^j={\beta }_l^{i+j}-{\alpha }_l^{i+j},i,j=0,1,\dots ,r,r\ge p.$$

$${\mathsf{E}}_{{\xi }_l}^{(1D)}={\mathbf{t}}_{{\beta }_l}{\mathbf{t}}_{{\beta }_l}^{\mathrm{T}}-{\mathbf{t}}_{{\alpha }_l}{\mathbf{t}}_{{\alpha }_l}^{\mathrm{T}}.$$

$${\mathbf{t}}_{{\alpha }_l}^{\mathrm{T}}{\mathbf{\xi }}_l^k={\alpha }_l^k,{\mathbf{t}}_{{\beta }_l}^{\mathrm{T}}{\mathbf{\xi }}_l^k={\beta }_l^k,k=0,1,\dots ,r.$$

$$\begin{array}[]{ll}\overline{\mathsf{P}}\equiv\mathsf{P}_{\xi_{1}}^{(1D)}\otimes\mathsf{P}_{\xi_{2}}^{(1D)}\otimes\mathsf{P}_{\xi_{3}}^{(1D)}\otimes\mathsf{I}_{5},&\\ \\ \overline{\mathsf{D}}_{\xi_{1}}\equiv\overline{\mathsf{P}}^{-1}\overline{\mathsf{Q}}_{\xi_{1}},&\overline{\mathsf{Q}}_{\xi_{1}}\equiv\mathsf{Q}_{\xi_{1}}^{(1D)}\otimes\mathsf{P}_{\xi_{2}}^{(1D)}\otimes\mathsf{P}_{\xi_{3}}^{(1D)}\otimes\mathsf{I}_{5},\\ \end{array}$$

$$\begin{array}[]{lll}\overline{\mathsf{E}}_{\xi_{1}}\equiv\overline{\mathsf{E}}_{\beta_{1}}+\overline{\mathsf{E}}_{\alpha_{1}},&\overline{\mathsf{E}}_{\alpha_{1}}\equiv-\overline{\mathsf{R}}_{\alpha_{1}}^{\mathrm{T}}\overline{\mathsf{P}}_{\perp\xi_{1}}\overline{\mathsf{R}}_{\alpha_{1}},&\overline{\mathsf{E}}_{\beta_{1}}\equiv\overline{\mathsf{R}}_{\beta_{1}}^{\mathrm{T}}\overline{\mathsf{P}}_{\perp\xi_{1}}\overline{\mathsf{R}}_{\beta_{1}},\\ \\ \overline{\mathsf{R}}_{\alpha_{1}}\equiv\bm{t}_{\alpha_{1}}^{\mathrm{T}}\otimes\mathsf{I}_{\xi_{2}}\otimes\mathsf{I}_{\xi_{3}}\otimes\mathsf{I}_{5},&\overline{\mathsf{R}}_{\beta_{1}}\equiv\bm{t}_{\beta_{1}}^{\mathrm{T}}\otimes\mathsf{I}_{\xi_{2}}\otimes\mathsf{I}_{\xi_{3}}\otimes\mathsf{I}_{5},&\overline{\mathsf{P}}_{\perp\xi_{1}}\equiv\mathsf{P}_{\xi_{2}}^{(1D)}\otimes\mathsf{P}_{\xi_{3}}^{(1D)}\otimes\mathsf{I}_{5},\end{array}$$

$${\mathsf{D}}_{{\xi }_1}\equiv {\mathsf{D}}_{{\xi }_1}^{(1D)}\otimes {\mathsf{I}}_{{\xi }_2}\otimes {\mathsf{I}}_{{\xi }_3}.$$

$$\begin{array}{rl}& \frac{\partial U}{\partial t}+\sum _{m=1}^3\frac{\partial F_{x_m}}{\partial x_m}=\sum _{m=1}^3\frac{\partial F_{x_m}^{(v)}}{\partial x_m},\forall \left(x_1,x_2,x_3\right)\in \Omega ,t\ge 0,\\ & U\left(x_1,x_2,x_3,t\right)=G^{(B)}\left(x_1,x_2,x_3,t\right),\forall \left(x_1,x_2,x_3\right)\in \Gamma ,t\ge 0,\\ & U\left(x_1,x_2,x_3,0\right)=G^{(0)}\left(x_1,x_2,x_3,0\right),\forall \left(x_1,x_2,x_3\right)\in \Omega ,\end{array}$$

$$U={\left[\rho ,\rho V_1,\rho V_2,\rho V_3,\rho E\right]}^{\mathrm{T}},$$

$$F_{x_m}={\left[\rho V_m,\rho V_m{V}_1+{\delta }_{m,1}P,\rho V_m{V}_2+{\delta }_{m,2}P,\rho V_m{V}_3+{\delta }_{m,3}P,\rho V_{m}H\right]}^{\mathrm{T}},$$

$$F_{x_m}^{(v)}={\left[0,{\tau }_{1,m},{\tau }_{2,m},{\tau }_{3,m},\sum _{i=1}^3{\tau }_{i,m}{U}_i-\kappa \frac{\partial T}{\partial x_m}\right]}^{\mathrm{T}}.$$

$${\tau }_{i,j}=\mu \left(\frac{\partial V_i}{\partial x_j}+\frac{\partial V_j}{\partial x_i}-{\delta }_{i,j}\frac{2}{3}\sum _{n=1}^3\frac{\partial V_n}{\partial x_n}\right),$$

$$H=c_{P}T+\frac{1}{2}{v}^{\mathrm{T}}v,P=\rho RT,R=\frac{R_u}{M_w},$$

$$s=\frac{R}{\gamma -1}\log \left(\frac{T}{T_{\infty }}\right)-R\log \left(\frac{\rho }{{\rho }_{\infty }}\right),\gamma =\frac{c_p}{c_p-R},$$

$$\begin{array}[]{l l}t&=\tau,\\ \bm{x}&=\bm{x}(\tau,\bm{\xi}),\end{array}$$

$$\frac{\partial(t,\bm{x})}{\partial(\tau,\bm{\xi})}=\left[\begin{array}[]{c c c c}1&0&0&0\\ \frac{\partial x_{1}}{\partial\tau}&\frac{\partial x_{1}}{\partial\xi_{1}}&\frac{\partial x_{1}}{\partial\xi_{2}}&\frac{\partial x_{1}}{\partial\xi_{3}}\\ \frac{\partial x_{2}}{\partial\tau}&\frac{\partial x_{2}}{\partial\xi_{1}}&\frac{\partial x_{2}}{\partial\xi_{2}}&\frac{\partial x_{2}}{\partial\xi_{3}}\\ \frac{\partial x_{3}}{\partial\tau}&\frac{\partial x_{3}}{\partial\xi_{1}}&\frac{\partial x_{3}}{\partial\xi_{2}}&\frac{\partial x_{3}}{\partial\xi_{3}}\end{array}\right],\quad J=\left|\frac{\partial(t,\bm{x})}{\partial(\tau,\bm{\xi})}\right|.$$

$$\begin{array}{rl}& \frac{\partial }{\partial t}=\frac{\partial }{\partial \tau }+\sum _{l=1}^3\frac{\partial {\xi }_l}{\partial t}\frac{\partial }{\partial {\xi }_l},\\ & \frac{\partial }{\partial x_m}=\sum _{l=1}^3\frac{\partial {\xi }_l}{\partial x_m}\frac{\partial }{\partial {\xi }_l},m=1,2,3,\end{array}$$

$$J\frac{\partial U}{\partial \tau }+\sum _{l=1}^{3}J\frac{\partial {\xi }_l}{\partial t}\frac{\partial U}{\partial {\xi }_l}+\sum _{m,l=1}^{3}J\frac{\partial {\xi }_l}{\partial x_m}\frac{\partial F_{x_m}}{\partial {\xi }_l}=\sum _{m,l=1}^{3}J\frac{\partial {\xi }_l}{\partial x_m}\frac{\partial F_{x_m}^{(v)}}{\partial {\xi }_l}.$$

$$\begin{array}[]{l}\frac{\partial J\bm{\mathit{U}}}{\partial\tau}+\sum\limits_{l=1}^{3}\frac{\partial}{\partial\xi_{l}}\left(J\frac{\partial\xi_{l}}{\partial t}\bm{\mathit{U}}\right)+\sum\limits_{m,l=1}^{3}\frac{\partial}{\partial\xi_{l}}\left(J\frac{\partial\xi_{l}}{\partial x_{m}}\left(\bm{\mathit{F}}_{x_{m}}-\bm{\mathit{F}}_{x_{m}}^{(v)}\right)\right)\\ -\bm{\mathit{U}}\left[\frac{\partial J}{\partial\tau}+\sum\limits_{l=1}^{3}\frac{\partial}{\partial\xi_{l}}\left(J\frac{\partial\xi_{l}}{\partial t}\right)\right]-\left(\bm{\mathit{F}}_{x_{m}}-\bm{\mathit{F}}_{x_{m}}^{(v)}\right)\sum\limits_{m,l=1}^{3}\frac{\partial}{\partial\xi_{l}}\left(J\frac{\partial\xi_{l}}{\partial x_{m}}\right)=0.\end{array}$$

$$\begin{array}[]{l}\frac{\partial J}{\partial\tau}+\sum\limits_{l=1}^{3}\frac{\partial}{\partial\xi_{l}}\left(J\frac{\partial\xi_{l}}{\partial t}\right)=0,\\ \sum\limits_{l=1}^{3}\frac{\partial}{\partial\xi_{l}}\left(J\frac{\partial\xi_{l}}{\partial x_{m}}\right)=0,\quad m=1,2,3,\end{array}$$

$$\frac{\partial JU}{\partial \tau }+\sum _{l=1}^3\frac{\partial }{\partial {\xi }_l}\left(J\frac{\partial {\xi }_l}{\partial t}U\right)+\sum _{m,l=1}^3\frac{\partial }{\partial {\xi }_l}\left(J\frac{\partial {\xi }_l}{\partial x_m}\left(F_{x_m}-F_{x_m}^{(v)}\right)\right)=0.$$

$$\begin{array}[]{l}\mathcal{S}=-\rho s,\\ \mathcal{F}=-\rho s\bm{V},\end{array}$$

$${\mathbf{W}}^{\mathrm{T}}\equiv \frac{\partial \mathcal{S}}{\partial U}={\left[\frac{h}{T}-s-\frac{{\mathbf{V}}^{\mathrm{T}}\mathbf{V}}{2T},\frac{V_1}{T},\frac{V_2}{T},\frac{V_3}{T},-\frac{1}{T}\right]}^{\mathrm{T}}.$$

$${\mathbf{W}}^{\mathrm{T}}\frac{\partial F_{x_m}}{\partial x_m}={\mathbf{W}}^{\mathrm{T}}\frac{\partial F_{x_m}}{\partial U}\frac{\partial U}{\partial x_m}=\frac{\partial {\mathcal{F}}_{x_m}}{\partial U}\frac{\partial U}{\partial x_m}=\frac{\partial {\mathcal{F}}_{x_m}}{\partial x_m},m=1,2,3,$$

$$\begin{array}{rl}& \frac{\partial U}{\partial W}\frac{\partial W}{\partial t}+\sum _{m=1}^3\frac{\partial F_{x_m}}{\partial W}\frac{\partial W}{\partial x_m}=\sum _{l,m=1}^3\frac{\partial }{\partial x_l}\left({\mathsf{C}}_{l,m}\frac{\partial W}{\partial x_m}\right),\end{array}$$

Full text

{AMS}

65M12, 65M70

1 Introduction

In many fluid dynamics applications including fluid-structure interaction, maneuvering, free-surface flows, and others, the domain boundaries undergo large displacements and deformations. One of the most widely used approaches for solving this class of problems is based on the arbitrary Lagrangian-Eulerian (ALE) framework that maps the Navier-Stokes equations in a reference frame wherein the computational grid can move and deform independently of the fluid flow [8, 16]. It should be noted that the grid motion and deformation may destroy both the conservation and stability properties of numerical methods based on the ALE formulation [20], thus indicating that a rigorous analysis is required for constructing stable ALE-based numerical schemes for the Euler and Navier-Stokes equations in moving and deforming domains.

In the ALE methodology, it is desirable to represent the Euler and Navier-Stokes equations in strong conservation form to satisfy the conservation laws  in time–dependent curvilinear coordinates. The continuous Euler and Navier-Stokes equations can be recast in strong conservation form if the geometric conservation law (GCL) equations are satisfied, which is always the case if a mapping between the physical and time-dependent curvilinear coordinates is diffeomorphism (e.g., see [11]). Though the GCL equations are satisfied exactly at the continuous level, this is not necessarily the case at the discrete level. A concept of the discrete geometric conservation laws was first introduced by Thomas and Lombard in [29]. To satisfy the discrete GCL equations exactly and improve the accuracy of the numerical solution on moving and deforming grids, the metric coefficients in [29] are first recast in conservation law form and then discretized by the same finite difference scheme used for approximating the Navier-Stokes equations.

It should be noted that satisfying the discrete GCL equations improves not only the accuracy but also stability of a numerical scheme on moving grids. Formaggia and Nobile show that for a linear advection-diffusion equation, the first-order discrete geometric conservation law is a sufficient condition for the backward Euler implicit scheme to be unconditionally stable. In [9], it is proven that for a nonlinear scalar hyperbolic conservation law equation, the discrete GCL equations are necessary and sufficient condition to guarantee that the 1st- and 2nd-order total variation diminishing (TVD) finite volume schemes are nonlinearly stable in the sense that they satisfy the discrete maximum principle.

Recently, provably stable ALE formulations have been developed for high-order finite difference, discontinuous Galerkin, and spectral collocation methods. Nikkar and Nordström use SBP operators both in time and space to construct high-order, fully discrete, conservative, energy stable finite-difference schemes for the skew-symmetric form of linear constant-coefficient hyperbolic systems in deforming domains [22]. This approach is extended to semi-discrete discontinuous Galerkin spectral collocation discretizations in [19]. For nonlinear conservation laws, the current state of the art is entropy stable schemes on static meshes. In [11, 2, 3], high-order entropy stable finite-difference and spectral collocation schemes based on summation-by-parts (SBP) operators that satisfy the discrete GCL equations are constructed for the 3-D unsteady Navier-Stokes equations on static curvilinear grids. These entropy stable methods have recently been extended to entropy stable spectral collocation weighted essentially non-oscillatory (WENO) [31] and fully discrete schemes [13] on static grids.

Despite the progress in ALE formulations, so far, no high-order provably stable numerical schemes have been developed for the unsteady compressible Euler and Navier-Stokes equations on dynamic grids. To address this problem, we develop new high-order spectral collocation schemes that are entropy stable, conservative, and satisfy the discrete GCL equations exactly on moving and deforming unstructured grids. The proposed scheme is constructed by discretizing the skew-symmetric form of the Navier-Stokes equations using SBP operators. A new entropy conservative flux is derived for the 3-D Euler and Navier-Stokes equations. The new flux preserves the design order of accuracy of the baseline scheme and provides entropy stability on moving and deforming grids.

The paper is organized as follows. The notation used in the paper is laid out in Section 2. The continuous geometric conservation law (GCL) equations and entropy stability analysis of the Navier-Stokes equations in moving and deforming domains are reviewed in Section 3. The numerical algorithm and its analysis are presented in Section 4. Implementational details are discussed in Section 5, while results are presented in Section 6. Finally, conclusions are drawn in Section 7.

2 Notation and definitions

Partial differential equations (PDEs) under consideration are discretized on cubes having Cartesian computational coordinates denoted by the 3-tuple $(\xi_{1},\xi_{2},\xi_{3})$, where the physical coordinates are denoted by the 3-tuple $(x_{1},x_{2},x_{3})$. Vectors are represented by lowercase bold font, for example

[TABLE]

while matrices are represented using sans-serif font, for example, $\mathsf{A}$. Continuous scalar functions on a space-time domain are denoted by capital letters in script font. For example,

[TABLE]

represents a square integrable function, where $t$ is the temporal coordinate. The restriction of such functions onto a set of mesh nodes is denoted by lower case bold font. For example, the restriction of $\mathit{U}$ onto a grid of $N_{1}\times N_{2}\times N_{3}$ nodes is given by the vector

[TABLE]

One-dimensional SBP operators are used to discretize derivatives and are extended to multiple dimensions using tensor products.

Definition 2.1**.**

The tensor product between the matrix $\mathsf{A}\in\mathbb{R}^{n,m}$ and $\mathsf{B}\in\mathbb{R}^{p,q}$ is denoted by

[TABLE]

For basic properties of tensor products, the reader is referred to [15]. The Hadamard product is used in the construction of entropy conservative/stable discretizations and is defined by

Definition 2.2**.**

The Hadamard product between the matrix $\mathsf{A}\in\mathbb{R}^{n\times n}$ and $\mathsf{B}\in\mathbb{R}^{n\times n}$ is denoted by

[TABLE]

For basic properties of Hadamard products, we refer the reader to [15].

SBP operators are matrix difference operators that approximate derivatives at mesh nodes and are characterized by their exactness in differentiating monomials, thus, an order $p$ SBP operator is one that differentiates exactly monomials up to degree $p$. Oftentimes monomials are discussed and the following notation is used:

[TABLE]

with the convention that $\bm{\xi}_{l}^{k}=\bm{0}$ for $k<0$.

We now give the definition of a one-dimensional SBP operator in the $\xi_{l}$ direction, $l=1,2,3$ (also, see [6, 7, 27]).

Definition 2.3**.**

**Summation-by-parts (SBP) operator for the first

derivative: A matrix operator, $\mathsf{D}_{\xi_{l}}^{(1D)}\in\mathbb{R}^{N_{l}\times N_{l}}$, is an SBP operator of order $p$ approximating the derivative $\frac{\partial}{\partial\xi_{l}}$ on the domain $\xi_{l}\in\left[\alpha_{l},\beta_{l}\right]$ with nodal distribution $\bm{\xi}_{l}$ having $N_{l}$ nodes, if**

$\mathsf{D}_{\xi_{l}}^{(1D)}\bm{\xi}_{l}^{k}=\left(\mathsf{P}_{\xi_{l}}^{(1D)}\right)^{-1}\mathsf{Q}_{\xi_{l}}^{(1D)}\bm{\xi}_{l}^{k}=k\bm{\xi}_{l}^{k-1}$, $k=0,1,\dots,p$; 2. 2.

$\mathsf{D}_{\xi_{l}}^{(1D)}=\left(\mathsf{P}_{\xi_{l}}^{(1D)}\right)^{-1}\mathsf{Q}_{\xi_{l}}^{(1D)}$, where the norm matrix, $\mathsf{P}_{\xi_{l}}^{(1D)}$, is symmetric and positive definite; 3. 3.

$\mathsf{Q}_{\xi_{l}}^{(1D)}=\left(\mathsf{S}_{\xi_{l}}^{(1D)}+\frac{1}{2}\mathsf{E}_{\xi_{l}}^{(1D)}\right)$, where $\mathsf{S}_{\xi_{l}}^{(1D)}=-\left(\mathsf{S}_{\xi_{l}}^{(1D)}\right)^{\mathrm{T}}$, $\mathsf{E}_{\xi_{l}}^{(1D)}=\left(\mathsf{E}_{\xi_{l}}^{(1D)}\right)^{\mathrm{T}}$, and $\mathsf{E}_{\xi_{l}}^{(1D)}$ satisfies

[TABLE]

In order to impose boundary conditions or inter-element coupling using simultaneous-approximation-term (SAT) penalty type conditions [4], it is useful to decompose $\mathsf{E}_{\xi_{l}}^{(1D)}$ as the outer product of interpolation/extrapolation operators [6]:

[TABLE]

The interpolation/extrapolation operators $\bm{t}_{\alpha_{l}}$ and $\bm{t}_{\beta_{l}}$ are $(r+1)$th–order accurate and therefore, satisfy the following accuracy conditions:

[TABLE]

Here, the three-dimensional Euler and Navier-Stokes equations are considered and the required derivative operators are given as

[TABLE]

where $\mathsf{I}_{5}$ is a $5\times 5$ identity matrix and $\overline{\mathsf{D}}_{\xi_{2}}$ and $\overline{\mathsf{D}}_{\xi_{3}}$ are defined similarly. The surface matrices are given as

[TABLE]

and the operators in the other two computational directions are defined similarly.

Remark 2.4**.**

It is at times necessary to consider scalar SBP operators and in such a case, the overbar notation is dropped, for example,

[TABLE]

Remark 2.5**.**

*While tensor-product SBP operators are used, all the theoretical results are presented in terms of multidimensional SBP operators and are therefore general. Along these lines, throughout the text, general nodal distributions are referenced as $C_{\kappa}\equiv\left\{\bm{\xi}_{i}\right\}_{i}^{N}$. *

The physical domain $\Omega\subset\mathbb{R}^{3}$, with boundary $\Gamma\equiv\partial\Omega$, with Cartesian coordinates $\left(x_{1},x_{2},x_{3}\right)\subset\mathbb{R}^{3}$, is partitioned into $K$ nonoverlapping hexahedral elements. The domain of the $\kappa^{\text{th}}$ element is denoted by $\Omega_{\kappa}$ and has boundary $\Gamma_{\kappa}$. Numerically, PDEs are solved in computational coordinates, $\left(\xi_{1},\xi_{2},\xi_{3}\right)\subset\mathbb{R}^{3}$, where each $\Omega_{\kappa}$ is locally transformed to $\hat{\Omega}_{\kappa}$, with boundary $\hat{\Gamma}\equiv\partial\hat{\Omega}_{\kappa}$, under the following assumption: {assume} Each element in physical space is transformed using a local, smooth, one-to-one, and invertible curvilinear coordinate transformation that is compatible at shared interfaces; thus, a grid generated by these local curvilinear coordinate transformations is watertight.

Associating the index 3-tuple $(i,j,k)$ with the computational coordinate 3-tuple

$(\xi_{1},\xi_{2},\xi_{3})$, the face numbering convention is given in Table 1.

3 Continuous analysis

3.1 Governing equations

The 3-D compressible Navier-Stokes equations in conservation law form in the Cartesian coordinates $(x_{1},x_{2},x_{3})$ are given by

[TABLE]

where $\bm{\mathit{U}}$ is a vector of conservative variables, and $\bm{\mathit{F}}_{x_{m}}$, and $\bm{\mathit{F}}_{x_{m}}^{(v)}$ are the inviscid and viscous fluxes associated with the $x_{m}$ coordinate, respectively. The boundary data, $\bm{\mathit{G}}^{(B)}$ and the initial condition, $\bm{\mathit{G}}^{(0)}$, are assumed to be bounded in $L^{2}\cap L^{\infty}$. In addition, $\bm{\mathit{G}}^{(B)}$ is assumed to contain boundary data that are stable in the entropy sense.

The vector of conservative variables is given as

[TABLE]

where $\rho$ denotes the density, $\bm{\mathit{V}}=\left[\mathit{V}_{1},\mathit{V}_{2},\mathit{V}_{3}\right]^{\mathrm{T}}$ is the velocity vector, and $\mathit{E}$ is the specific total energy. The inviscid fluxes, $\bm{\mathit{F}}_{x_{m}},m=1,2,3$, are given by

[TABLE]

where $\mathit{P}$ is the pressure, $\mathit{H}$ is the specific total enthalpy and $\delta_{i,j}$ is the Kronecker delta (we have broken our convention for continuous scalar functions in the definition of $\rho$).

The viscous fluxes, $\bm{\mathit{F}}_{x_{m}}^{(v)},m=1,2,3$, are defined as

[TABLE]

The viscous stresses are given by

[TABLE]

where $\mu(T)$ is the dynamic viscosity and $\kappa(T)$ is the thermal conductivity (not to be confused with the choice of parameter for element numbering).

To close the Navier-Stokes equations, Eq. (3), the following constituent relations are used:

[TABLE]

where $\mathit{T}$ is the temperature, $R_{u}$ is the universal gas constant, $M_{w}$ is the molecular weight of the gas, and $c_{\mathit{P}}$ is the specific heat capacity at constant pressure. Finally, the specific thermodynamic entropy is given as

[TABLE]

where $\mathit{T}_{\infty}$ and $\rho_{\infty}$ are reference temperature and density, respectively (the stipulated convention has been broken here and $s$ has been used rather than ${S}$ for reasons that will be discussed in Section 3.3).

3.2 Geometric Conservation Law (GCL) Equations

To solve the Navier-Stokes equations in complex geometries with moving and deforming boundaries, we use the arbitrary Lagrangian-Eulerian (ALE) formulation. A dynamic unstructured grid in the physical domain is generated by individually mapping a reference domain $(\xi_{1},\xi_{2},\xi_{3})\in\hat{\Omega}=[-1,1]^{3}$ with time $\tau$ onto each grid element in the physical domain $(x_{1},x_{2},x_{3})\in\Omega$ with time $t$. Assuming that each individual transformation

[TABLE]

is a diffeomorphism, it can be described by the following Jacobian matrix:

[TABLE]

Taking into account that

[TABLE]

and multiplying Eq. (3) by the metric Jacobian, $J$, the Navier-Stokes equations can be recast in the time-dependent curvilinear coordinates $(\xi_{1},\xi_{2},\xi_{3},\tau)$ as follows:

[TABLE]

In contrast to the original Navier-Stokes equations in the Cartesian coordinates, Eq. (7) is not written in divergence from. Using the product rule, the Navier-Stokes equations in the curvilinear coordinates can be represented in conservation law from as follows:

[TABLE]

Based on the definition of the metric coefficients, the following identities hold:

[TABLE]

which are called the geometric conservation law (GCL) equations. Taking into account the above identities, the last two terms in Eq. (8) vanish and the Navier-Stokes equations can be rewritten in the following fully conservative form:

[TABLE]

Note that the GCL equations (9) guarantee that any physically meaningful constant vector of conservative variables $\bm{\mathit{U}}=\bm{const}$ is a solution of the Navier-Stokes equations (10). Though, the GCL equations (9) are satisfied exactly at the continuous level, this is not the case at the discrete level [29]. A detailed analysis on how the Navier-Stokes equations and the corresponding metric coefficients should be discretized to satisfy the GCL equations is presented in Section 4.

3.3 Continuous Entropy Analysis

A necessary condition for selecting a unique, physically relevant solution among possibly many weak solutions of Eq. (3) is the entropy inequality. It is well known that the entropy inequality holds for the Navier-Stokes equations in the Cartesian and curvilinear coordinates that are independent of time (e.g., see [3]). In the present analysis, we show that the entropy inequality also holds for the Navier-Stokes equations in time-dependent curvilinear coordinates.

The compressible Navier-Stokes equations are equipped with a convex scalar entropy function $\mathcal{S}$ and the corresponding entropy flux $\mathcal{F}$, which are given by

[TABLE]

where $s$ is the thermodynamic entropy defined by Eq. (6) and $\bm{V}$ is the velocity vector. Note that the mathematical entropy $\mathcal{S}$ has the opposite sign from the thermodynamic entropy. Thus, the mathematical entropy across a shock decreases rather than increases. This nomenclature is used throughout the paper.

The entropy function $\mathcal{S}$ satisfies the following properties:

$\mathcal{S}(\bm{\mathit{U}})$ is convex and its Hessian matrix, $\frac{\partial^{2}\mathcal{S}}{\partial U^{2}}$, is positive definite provided that $\rho>0$ and $T>0$ $\forall\bm{x}\in\Omega$, thus yielding a one-to-one mapping from the conservative to entropy variables that are defined as follows:

[TABLE] 2. 2.

The entropy variables satisfy the following compatibility relations for all inviscid fluxes of the compressible Navier-Stokes equations:

[TABLE]

where $\mathcal{F}_{x_{m}}$ is the entropy flux in the $m$-th spatial direction. 3. 3.

The entropy variables symmetrize the compressible Navier-Stokes equations, which can be recast in terms of $\bm{W}$ as follows:

[TABLE]

with the symmetry conditions $\frac{\partial U}{\partial W}=\left(\frac{\partial U}{\partial W}\right)^{\mathrm{T}}$, $\frac{\partial\bm{\mathit{F}}_{x_{m}}}{\partial U}=\left(\frac{\partial\bm{\mathit{F}}_{x_{m}}}{\partial U}\right)^{\mathrm{T}}$, and $\mathsf{C}_{l,m}=(\mathsf{C}_{l,m})^{\mathrm{T}}$. Furthermore, $\frac{\partial U}{\partial W}$ is positive definite, and the matrices $\mathsf{C}_{l,m}$ are positive semi-definite, provided that $\rho>0$ and $T>0$ $\forall\bm{x}\in\Omega$. Note that the term on the right–hand side of Eq. (14) is a recast of the viscous fluxes in terms of entropy variables, that is,

[TABLE]

It has been proven by Godunov in [14] that if Eq. (3) is symmetrized by introducing new variables $\bm{W}$ and $\bm{W}$ is a convex function of $\varphi$, then the entropy function and the corresponding entropy flux satisfy the following equations:

[TABLE]

[TABLE]

where the functions $\varphi$ and $\bm{\psi}$ are called the entropy potential and entropy potential flux, respectively.

We now show that the entropy inequality holds for the compressible Navier-Stokes equations in the time-dependent curvilinear coordinates. Contracting Eq. (10) with the entropy variables given by Eq. (12) yields

[TABLE]

The matrices $\hat{\mathsf{C}}_{l,n}$ on the right-hand side of Eq. (18) are symmetric semi-definite matrices which are given by

[TABLE]

For further details on how $\mathsf{C}_{m,j}$ and $\hat{\mathsf{C}}_{l,n}$ are constructed, we refer the reader to [11].

Using $\bm{\mathit{W}}^{\mathrm{T}}=\frac{\partial\mathcal{S}}{\partial\bm{\mathit{U}}}$, the term $I$ in Eq. (18) can be manipulated as follows:

[TABLE]

Similarly, the term $II$ is recast in the following form:

[TABLE]

Using the compatibility relations (Eq. (13)), the term $III$ is reduced to

[TABLE]

The last term in Eq. (18) can be manipulated as follows:

[TABLE]

Integrating Eq. (18) over the physical domain and taking into account Eqs. (20–23), we have

[TABLE]

The last terms in the first two lines of Eq.(24) are identically equal to zero because of the GCL equations given by Eq. (9). Using the integration-by-parts (IBP) formula, the above equation can be recast in the following form:

[TABLE]

Taking into account that the matrices $\hat{\mathsf{C}}_{l,m}$ are positive semi-definite and assuming that the boundary conditions are entropy stable, Eq. (25) becomes

[TABLE]

Equation (26) represents the entropy inequality in the domain. Note that for the Euler equations with smooth solutions, Eq. (25) becomes an equality. The entropy inequality (25) is only a necessary condition, which is not by itself sufficient to guarantee the convergence to a physically relevant weak solution of the Navier-Stokes equations.

4 Discrete analysis

In this section, the various components of the entropy stable semi-discrete form for the compressible Navier-Stokes equations are presented. The approach taken is to first construct an entropy conservative semi-discrete form for the Euler equations, Section 4.1, then add appropriate interface dissipation to make the scheme entropy stable, Section 4.2, and finally, discretize the viscous components such that entropy stability is not lost, Section 4.3.

4.1 Entropy conservative semi-discrete form for the Euler equations

The skew-symmetric splitting of the Euler equations that is discretized on the $\kappa^{\rm{th}}$ element is given as

[TABLE]

The skew–symmetric form is derived by averaging the conservative and nonconservative forms. Note that the conservative form is derived by taking advantage of the GCL conditions outlined in Section 3.2:

[TABLE]

[TABLE]

In order to achieve entropy conservation, the approximations to the spatial derivatives are constructed such that the continuous entropy stability analysis is mimicked in a one-to-one fashion. To do so, the approximation to the spatial derivatives needs to be mimetic of a special case of integration by parts (IBP), namely

[TABLE]

(this is referred to as a nonlinear IBP relation) where $\bm{\mathit{V}}$ and $\bm{\mathit{G}}\left(\bm{\mathit{U}}\right)$ are vector valued functions, $\mathit{h}\left(\bm{\mathit{U}}\right)$ is a scalar function, $\zeta$ is some independent variable over the domain $D$ with boundary $\Gamma$, and $n_{\zeta}$ is the component of the outward facing unit normal in the $\zeta$ direction.

Remark 4.1**.**

In the current context, the two nonlinear IBP formulas that are of interest are

[TABLE]

*where the GCL conditions (28) and (29) have been used. *

To introduce nonlinear SBP operators and their associated notation, the generic nonlinear IBP property (30) and the derivative $\frac{\partial\bm{\mathit{G}}\left(\bm{\mathit{U}}\right)}{\partial\zeta}$ are used; the derivative is approximated using the following nonlinear SBP operator:

[TABLE]

The symbol $\circ$ is the Hadamard product, i.e., element-wise multiplication while $C_{\kappa}$ is the set of nodes on the $\kappa^{\rm{th}}$ element; therefore, the notation

[TABLE]

denotes the vector constructed by evaluating $\frac{\partial\bm{\mathit{G}}\left(\bm{\mathit{U}}\right)}{\partial\zeta}$ at the nodes $C_{\kappa}$. The vector $\bm{u}_{\kappa}$ is a vector of vectors corresponding to the discrete solution of the $5$ conservation equations at each node ordered with the $\xi_{3}-$direction varying most rapidly, then $\xi_{2}$, and then $\xi_{1}$. The vector $\overline{\bm{1}}$ is defined as $\overline{\bm{1}}\equiv\bm{1}\otimes\bm{1}_{5}$, where $\bm{1}$ is a vector of ones of size $N$, the number of nodes in each element, and $\bm{1}_{5}$ is a vector of ones of size $5$. It is necessary to discuss the components of $\bm{u}_{\kappa}$ at each node and the values of conservative variables at all nodes in an element. For this purpose a matrix like notation, $\bm{u}_{\kappa}\left(i,\,j\right)$, is introduced to index the vector such that $\bm{u}_{\kappa}\left(i,\,:\right)=\bm{u}_{\kappa}((i-1)5:i5)$ and $\bm{u}_{\kappa}\left(:,\,j\right)$ is the vector constructed from the $j^{\rm{th}}$ conserved variable at all nodes. Furthermore, the two-argument matrix function $\mathsf{G}^{sc}\left(\bm{u}_{\kappa},\bm{u}_{\kappa}\right)$ is a block matrix of diagonal blocks, that is

[TABLE]

and it is of size $\left(5\,N\right)\times\left(5\,N\right)$. The five element column vectors,

[TABLE]

are constructed from two-point flux functions that are symmetric in their arguments and consistent, thus,

[TABLE]

where the superscript denotes that this is an entropy conservative two-point flux function. This implies that

[TABLE]

Finally, in order to construct entropy conservative schemes, the two-point flux functions must satisfy the Tadmor shuffle condition [28]:

[TABLE]

The next theorem demonstrates that the nonlinear SBP operator is mimetic of the nonlinear IBP relation (30) to high order [10]:

Theorem 4.2**.**

Consider a two-argument matrix flux function $\mathsf{G}^{sc}\left(\bm{u}_{\kappa},\bm{u}_{r}\right)$ constructed from the two point flux function $\bm{\mathit{G}}^{sc}\left(\bm{u}_{\kappa}\left(i,\,:\right),\bm{u}_{r}\left(j,\,:\right)\right)$ that is symmetric, consistent, and satisfies the Tadmor shuffle condition (33). Then the nonlinear SBP operator is an approximation to the derivative as below,

[TABLE]

and satisfies the following nonlinear SBP property,

[TABLE]

*where the LHS and RHS are high-order approximations to the LHS and RHS of the nonlinear IBP property (30). *

Proof 4.3**.**

*The proof is given in [10]. *

Theorem 4.2 is proven using accuracy estimates on the nonlinear SBP operators and the bilinear forms that result from its constituent components. In addition to Theorem 4.2, the following theorem (which will be used to prove entropy conservation) is employed:

Theorem 4.4**.**

Consider the matrix $\mathsf{A}$ of size $N\times N$ and a two-argument matrix flux function $\mathsf{G}^{sc}\left(\bm{u}_{\kappa},\bm{u}_{r}\right)$ constructed from the two point flux function $\bm{\mathit{G}}^{sc}\left(\bm{u},\bm{v}\right)$ that satisfies the Tadmor shuffle condition (33) and is symmetric, then

[TABLE]

*where $\bm{v}^{\mathrm{T}}\bm{g}$ is to be understood as $\left(\bm{v}^{\mathrm{T}}\bm{g}\right)(i)\equiv\bm{v}(:,\,i)^{\mathrm{T}}\bm{g}(:,\,i)$. *

Proof 4.5**.**

*The proof is given in [10]. *

The discrete operators that are used to discretize the terms in (27) and their error properties are given in the next theorem.

Theorem 4.6**.**

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\overline{\left[\mathit{J}\frac{\partial\xi_{l}}{\partial x_{m}}\right]}_{\kappa}$ and $\overline{\left[\mathit{J}\frac{\partial\xi_{l}}{\partial t}\right]}_{\kappa}$ are diagonal matrices with the metric terms along their diagonals and $\mathsf{F}_{x_{m}}^{sc}$ and $\mathsf{U}^{sc}$ are constructed from two-point flux functions that are symmetric, consistent, and satisfy the following shuffle conditions:

[TABLE]

[TABLE]

Proof 4.7**.**

*The proofs of (35) and (36) are given in [10], while the proofs of (37) and (38), follow identically. *

Using the operators in Theorem 4.6, the semi-discrete form of (27) is given as

[TABLE]

where the SATs on the right-hand side are defined as

[TABLE]

[TABLE]

where $\bm{u}_{(2l-1)}$ and $\bm{u}_{2l}$ are the solution vectors for the elements touching the $(2l-1)$ and $2l$ faces, respectively.

For analysis, it is convenient to separate out the discretization of each scalar conservation law in (27) and for this purpose the following is introduced:

[TABLE]

where the SATs on the right-hand side are defined as

[TABLE]

[TABLE]

where

[TABLE]

Similarly,

[TABLE]

In order to be entropy conservative, the approximation to the metric terms used in (41) need to satisfy the conditions contained in the next theorem.

Theorem 4.8**.**

The discretization (41) has a telescoping entropy on interior elements if the metric terms satisfy

[TABLE]

and

[TABLE]

This implies that for the periodic problem, the scheme is entropy conservative, i.e.,

[TABLE]

*and entropy conservative/stable for non-periodic problems under the assumption of appropriate SATs are used for the weak imposition of boundary conditions. *

Proof 4.9**.**

*The proof is contained in Appendix A *

Next, let us show that the proposed scheme is freestream preserving.

Theorem 4.10**.**

*Under the same conditions on the metric terms as in Theorem 4.8 the discretization (41) is freestream preserving. *

Proof 4.11**.**

The proof follows similarly to that given in Ref. [10]. Here, we sketch the required steps. It is convenient to use the scalar version of the semi-discrete form (44). Consider a physically meaningful constant state, $\bm{\mathit{U}}_{c}$, then we denote the evaluation of this constant state on the mesh nodes of the $\kappa^{\rm{th}}$ element as $\bm{u}_{\kappa}^{(c)}$; then, via the consistency of the two-point flux functions, the derivative operators in the scalar version of semi-discrete form reduce to

[TABLE]

and

[TABLE]

*In a similar manner, it is easy to show that the SATs vanish and the resulting conditions on the metric terms are the same as in Theorem 4.8. *

Remark 4.12**.**

*There is a misconception that the freestream preservation is sufficient to guarantee stability on moving and deforming grids. In fact, the freestream preservation property is only a necessary condition for stability. *

In the present analysis, the interest is in schemes that are discretely conservative so that when they converge, they converge to a weak solution of the conservation law equations. A proof that the scheme has the desired property is based on a semi-discrete analog of the Lax-Wendroff theorem. The conservation property can be demonstrated by showing that the scheme has a telescoping flux form. Because the nonlinear SBP operators are used in the construction of the semi-discrete form (41), a proof of a “stronger” discrete conservation statement, namely subcell conservation, is given herein.

Theorem 4.13**.**

*The semi-discrete scheme (41) is subcell conservative. *

Proof 4.14**.**

The proof is sketched here. The nonlinear SBP operators can be recast in telescoping flux form. For example, for the one-dimensional nonlinear operator on $N$ nodes

[TABLE]

where the rectangular matrix $\mathsf{\Delta}$ of size $(N+1)\times N$ is given as

[TABLE]

*The vector $\bm{\mathit{f}}_{x_{m}}$ is an interpolation of $\mathit{F}_{x_{m}}$ to a set of nodes that interdigitate the nodal distribution of the element. The nonlinear SBP operators constructed globally over the mesh are then recast in a telescoping flux form over each mesh line. The resulting finite-volume like scheme is in the form required by the Lax-Wendroff theorem and therefore, the discretization is subcell conservative (for more details, see Refs. [12, 11]). *

4.2 Entropy stable semi-discrete form for the Euler equations: interface dissipation

The weak imposition of boundary conditions using SATs can be reinterpreted in the sense of discontinuous Galerkin methods by discretely integrating against a test function $\mathit{V}$; thus,

[TABLE]

In the context of the entropy conserving SATs thus far introduced, for example on the $2l$ interface, the weak form of the SAT can be reinterpreted as

[TABLE]

and

[TABLE]

In this section, the entropy conservative SATs are augmented with dissipative terms motivated by the upwinding used in the Roe approximate Riemann solver, which has generic form

[TABLE]

where the superscript $+$ denotes quantities evaluated on a side of an interface in the positive direction of the unit normal, and the superscript $-$ denotes quantities in the negative direction. Since the Roe flux is not entropy consistent, the central flux,

[TABLE]

is replaced by the numerical fluxes given in the above entropy conservative SAT terms. To provide entropy dissipation at element interfaces, the dissipative term, i.e.,

[TABLE]

should also be modified.

The approach is to construct a dissipation term that enables entropy stability, while remaining design order accurate and conservative. This is accomplished by using the flux Jacobian with respect to the entropy variables rather than the conservative variables. The eigenvectors of the conservative variable flux Jacobian can be scaled such that

[TABLE]

Thus, for the mesh velocity terms such as $\mathit{J}_{\kappa}\frac{\partial\xi_{l}}{\partial t}\frac{\partial\bm{\mathit{U}}}{\partial\xi_{l}}$, the added dissipative term is of the form

[TABLE]

For the remaining spatial terms, note that

[TABLE]

this property was first derived by Merriam [21] (for more details, see [11]). Thus, the added dissipative terms are

[TABLE]

where

[TABLE]

and the matrices $\mathsf{Y}$ and $\mathsf{\Lambda}_{\xi_{l}}$ are constructed pointwise from the Roe-average of the two states on either side of the interface. The next theorem characterizes the addition of the dissipation to the scheme.

Theorem 4.15**.**

*The added dissipative terms (49) are design order, lead to an entropy stable scheme, and retain subcell conservation. *

Proof 4.16**.**

*The proofs of design order and entropy stability follow similarly to those given in Ref. [10]; subcell conservation follows by recasting the dissipation terms in telescoping flux form. *

4.3 Extension to the Navier-Stokes equations

A local discontinuous Ga-lerkin (LDG)-type approach is used to discretize the viscous terms (see [2, 5, 24, 23]). Thus, Eq.(10) is discretized as follows:

[TABLE]

where

[TABLE]

The terms $\bm{f}_{I}$ and $\bm{SAT}_{I}$ are the discretization of the inviscid portion of the governing equations. The interior penalty term, $\bm{I}_{\mathrm{P}}^{\kappa}$, adds interface dissipation and is given as

[TABLE]

where the notation, for example $\left[\hat{\mathsf{C}}_{l,l}\right]_{\kappa,2l}$, implies that $\hat{\mathsf{C}}$ is constructed from the Roe average of the solutions $\bm{u}_{\kappa}$ and $\bm{u}_{2l}$.

The next theorem demonstrates that the discretization of the viscous terms leads to an entropy stable formulation.

Theorem 4.17**.**

*Assuming that the inviscid portion is entropy conservative, then the semi-discretization (50) is entropy stable and freestream preserving. *

Proof 4.18**.**

The proof of entropy stability follows similarly to that given in Ref. [10]. Here, we sketch the steps in the proof:

• •

Contract the viscous portions of Eq. (50) with $\bm{w}_{\kappa}^{\mathrm{T}}\overline{\mathsf{P}}$, contract Eq. (51) with $\left(\bm{\theta}_{n}^{\kappa}\right)^{\mathrm{T}}\left[\hat{\mathsf{C}}_{l,n}\right]_{k}\overline{\mathsf{P}}$ for $n=1,2,3$, and add the result.

• •

After simplification, the result is of the form

[TABLE]

where the first term on the right-hand side can be shown to be negative semi-definite and $CT$ represents the coupling terms from the SATs.

• •

Summing over the elements, it can be shown that the interface coupling terms cancel and assuming appropriate boundary SATs, the remaining terms at the boundaries are consistent with the continuous analysis.

*The proof of freestream preservation is identical to that of Theorem 4.10, taking into account the fact that all viscous fluxes and SAT terms vanish if the state vector is constant. *

5 Implementation

5.1 Approximation of the metrics

As follows from Theorem 4.8, the discrete GCL equations (47–48) must be satisfied exactly. Note that this property is not immediately satisfied. Several approaches have been proposed in the literature to address this problem [29, 30, 1, 26]. These studies suggest that the following three properties should hold to guarantee that the GCL equations are satisfied exactly at the discrete level:

The discrete operators used to approximate the derivatives of the inviscid fluxes should also be used to approximate the metric coefficients; 2. 2.

The metric tensor $\bm{\xi}_{\bm{x}}$ cannot be calculated directly by inverting $\bm{x}_{\bm{\xi}}$, but instead it should be discretized so that it satisfies the GCL equation (48); 3. 3.

The grid Jacobian should be computed by integrating the first discrete GCL equation (47) using the same time integrator employed for approximating the governing equations.

To satisfy the first two properties, the metrics in the present analysis are discretized by using the approach developed in [29], which is selected because of its computational efficiency. Using the SBP spectral collocation operators defined in Section 2, the metric tensor $\bm{x}_{\bm{\xi}}$ is discretized as follows:

[TABLE]

The inverse metric coefficients are then calculated using the following formulas:

[TABLE]

where the Jacobian $J$ is the determinant of the metric tensor $\bm{x}_{\bm{\xi}}$ given by Eq. (53). The metric coefficients given by Eq. (54) satisfy the second discrete GCL equation (48) exactly. Though there are other formulas for calculating the metric coefficients that identically satisfy Eq. (48) (e.g., see [1, 26]), these alternative approaches are more computationally expensive and have therefore not been considered in the present analysis.

As has been noted above, the grid Jacobian $J$ should satisfy the GCL equations (47–48) exactly. There are at least three possible options for calculating $J$, which may provide this property.

The first option is to discretely integrate Eq. (47) by using the same temporal and spatial derivative operators used for discretization of the governing equations. This numerically integrated Jacobian is then used to calculate the inverse metric coefficients given by Eq. (54). Though this approach satisfies Eq. (48) exactly, Eq. (47) may not be satisfied identically, because the integrated and determinant-based Jacobians differ from each other by the magnitude of the order of the truncation error. Our numerical results show that this integrated Jacobian approach provides the freestream preservation only with the accuracy that is 2–3 orders of magnitude larger than the machine zero. This is the main reason why this approach is not used in the present analysis.

The second approach uses the integrated and determinant-based Jacobians simultaneously. The key difference of this approach from the previous one is that both the spatial terms in the discretized Navier-Stokes equations (41) and the mesh velocity terms in Eq. (47) are calculated by using the determinant-based Jacobian, while the Jacobian in the time derivative term is computed by integrating Eq. (47). This method uses the conservative discretization in time, satisfies both GCL equations exactly, and provides the freestream preservation to machine accuracy. This dual-Jacobian approach is used in all numerical experiments presented in Section 6.

The last option considered is based on a nonconservative approximation of the time derivative term in Eq. (54). Using the product rule, the time derivative term can be represented as follows:

[TABLE]

In this approach, the grid Jacobian is evaluated by calculating the determinant of the metric tensor $\bm{x}_{\bm{\xi}}$. The time derivative of the Jacobian, $\frac{\mathrm{d}\overline{\mathsf{J}}_{\kappa}}{\mathrm{d}\tau}$, is computed using the GCL equation (47). Though this approach satisfies both GCL equations exactly, it is nonconservative in time and therefore not used in the present analysis.

5.2 Entropy conservative fluxes

As follows from Theorem 4.6, the two-point entropy conservative flux functions $\bm{\mathit{U}}^{sc}$ and $\bm{\mathit{F}}^{sc}_{x_{m}}$ should be symmetric, consistent, and satisfy the shuffle conditions (39) and (40), respectively. For the Euler and Navier-Stokes equations, a computationally affordable entropy consistent flux function $\bm{\mathit{F}}^{sc}_{x_{m}}$ satisfying all the required constraints has been derived by Ismail and Roe in [17]. It should be noted that this two-point entropy consistent flux function $\bm{\mathit{F}}^{sc}_{x_{m}}$ is not unique. We refer the reader to [25] for a review of alternative entropy consistent flux functions satisfying (40).

Using the same technique developed in [17] for constructing $\bm{\mathit{F}}^{sc}_{x_{m}}$, the following two-point entropy consistent fluxes satisfying the shuffle condition (39) can be derived:

[TABLE]

[TABLE]

where $\mathit{Z}=\left[\frac{1}{\sqrt{T}},\frac{V_{1}}{\sqrt{T}},\frac{V_{2}}{\sqrt{T}},\frac{V_{3}}{\sqrt{T}},\rho\sqrt{T}\right]^{\mathrm{T}}$ for $\bm{\mathit{U}}_{1}^{sc}$ and $\mathit{Z}=\left[\rho,V_{1},V_{2},V_{3},\frac{1}{T}\right]^{\mathrm{T}}$ for $\bm{\mathit{U}}_{2}^{sc}$. In the above equations, $\overline{Z}$ and $\overline{Z}_{\rm{log}}$ are the arithmetic and logarithmic averages:

[TABLE]

where the subscripts $L$ and $R$ correspond to two arbitrary, yet physically meaningful, states of the flow variables.

Lemma 5.1**.**

*The two-point fluxes given by Eqs. (55) and (56) are symmetric, consistent, and satisfy the shuffle condition (39). *

Proof 5.2**.**

*The symmetry of the entropy consistent fluxes $\bm{\mathit{U_{1}}}^{sc}$ and $\bm{\mathit{U_{2}}}^{sc}$ follows immediately form the symmetry of the arithmetic and logarithmic averages given by Eq. (57). The consistency can be directly verified by substituting $Z_{L}=Z_{R}=Z$ in Eqs. (55–55) and using the L’Hopital’s rule for evaluating the logarithmic averages. The shuffle condition (39) can be verified directly by substituting Eqs. (55–55) in Eq. (39). *

Remark 5.3**.**

*As evident from Eqs. (55–55), the two-point entropy consistent flux $\bm{\mathit{U}}^{sc}$ is not unique and other two-point flux functions satisfying Eq. (39) can be readily constructed by using different auxiliary flow variables $Z$. *

6 Results

The accuracy and freestream preservation properties of the proposed high-order entropy stable spectral collocation schemes for the 3-D Euler and Navier-Stokes equations on moving and deforming grids are tested using standard benchmark problems including isentropic vortex and viscous shock flows. In all simulations presented herein, the low-storage, five-stage, fourth-order Runge-Kutta scheme developed in [18] is used to advance the semi-discretization in time. Note that this scheme violates the entropy stability property of the semi-discrete operator by a factor proportional to the local temporal truncation error. In the present analysis, the temporal error is evaluated and monitored by using the Runge-Kutta built-in error estimator. For all test problems considered, the temporal error is at least 2-3 orders of magnitude less than the corresponding spatial error. The spatial derivatives in the Euler and Navier-Stokes equations are discretized using the entropy stable spectral collocation schemes with two-point entropy conservative flux given by Eq. (56) and the polynomial bases of various degrees ($p=3,4,5$). It should be emphasized that no additional artificial dissipation other than that introduced by the entropy consistent Roe-type flux (see Eq. (49)) at element interfaces is used for all numerical experiments presented in this paper.

6.1 Inviscid flows

To evaluate the accuracy of the new entropy stable spectral collocations schemes for the Euler equations, we consider the propagation of an isentropic vortex in a moving deforming domain. The isentropic vortex is an exact solution of the Euler equations, which is given by

[TABLE]

[TABLE]

For this test problem, we set $M_{\infty}=0.5,V_{\infty}=0.25,\epsilon=5.0,x_{0}=-0.25,y_{0}=0,\alpha=0$.

The dynamic grid is constructed by mapping vertices of each grid cell of the uniform Cartesian grid in the computational domain $\hat{\Omega}=[-1,1]\times[-1,1]\times[-1,1]$ onto the corresponding vertices of a moving deforming hexahedral grid element by using the following coordinate transformation:

[TABLE]

where $K$ is the number of elements in the $\xi-$ and $\eta-$directions and the parameters $A_{1}$, $A_{2}$, $\nu$, and $\omega$ are set equal to $0.09$, $0.06$, $3\pi/4$, and $\pi$, respectively. Note that only the vertices of each hexahedral grid element are moved according to Eq. (59), while the element edges are constructed by connecting the corresponding vertices with straight lines and the interior grid points are generated by using the 3-D transfinite interpolation. As a result, the grid is nonsmooth across all element interfaces.

Though this test problem is two-dimensional in nature, a 3-D code is used to solve this isentropic vortex flow for testing the proposed entropy stable spectral collocation schemes on dynamic unstructured grids. In this study, the Dirichlet boundary conditions are imposed at the moving boundaries along the $\xi-$ and $\eta-$directions by using the exact solution given by Eq. (58) and the periodic boundary conditions are used in the $\zeta$-direction. The mesh velocity vector appearing in both the proposed numerical scheme (41) and the GCL equation (28) is computed analytically by differentiating the time-dependent coordinate transformation given by Eq. (59).

At the initial moment of time $t=0$, the vortex is centered at $x_{0}=-0.25,y_{0}=0$ and the grid is not deformed. The final time is set equal to $2.5$, so that the grid undergoes $1\frac{1}{4}$ deformation cycles and reaches its maximum deformation, while the vortex propagates throughout the entire domain. Density contours obtained using the new entropy stable spectral collocation scheme for $p=4$ on the $12\times 12\times 1$ dynamic grid at $t=0$ and $t=2.5$ are shown in Fig. 1. As one can see in the figure, the grid elements are highly deformed and skewed at $t=2.5$. Despite this large unsteady grid deformations, the numerical solution is free of spurious oscillations and preserves the symmetry of the vortex.

The error convergence rate of the new high-order entropy stable spectral collocation schemes is evaluated on a sequence of moving deforming grids with $6\times 6\times 1$, $12\times 12\times 1$, $24\times 24\times 1$, $48\times 48\times 1$, elements. The refined dynamic grids are obtained by constructing a sequence of nested grids in the computational domain and then mapping them onto the physical domain by using Eq. (59). As mentioned above, all 3-D grids used in this grid refinement study are one–element thick in the $\zeta$ direction and their thickness is proportionally reduced as the grid is globally refined.

The $L_{2}$ and $L_{\infty}$ errors of the entire solution vector obtained with the new entropy stable schemes for polynomial bases of degree $p=3,4,5$ on globally refined dynamic grids are presented in Table 2. The convergence rate in the range of $p$ to $p+0.5$ is achieved on highly deformed, skewed moving grids for all polynomial degrees considered. It should be noted that the $(p+1)$th-order convergence or so-called superconvergence is not attainable on highly deformed skewed moving meshes.

We now corroborate our theoretical results presented in Theorem 4.10 and show that the proposed spectral collocation schemes for the Euler equations preserve the freestream solution on the same dynamic grids used for the isentropic vortex problem. Figure 2 shows the $L_{\infty}$ error of the entire vector of conservative variables computed using the new entropy stable schemes (41) for $p=3,4,5$ on the $12\times 12\times 1$ grid given by Eq. (59) for the constant inviscid flow at $M_{\infty}=0.5$. As evident from this comparison, the constant solution is preserved on the moving deforming grid with the accuracy of $O(10^{-14})$ for all polynomial degrees considered over the entire interval of integration. As one can see in Fig. 2, there is a slight increase in error as the polynomial degree increases. This increase can be explained by accumulation of the roundoff error, since the higher order schemes require larger number of operations per degree of freedom.

6.2 Viscous flows

We now evaluate the accuracy and freestream preservation properties of the proposed entropy stable spectral collocation schemes for the Navier-Stokes equations on moving deforming unstructured hexahedral grids. To assess the effect of grid motion and deformation on the error convergence of the high-order spectral collocation schemes, the propagation of a viscous shock is considered. For Prandtl number equal to $3/4$, the 1-D compressible Navier-Stokes equations reduce to the single momentum equation, which can be solved exactly and used to measure the accuracy of the numerical scheme under consideration (e.g., see [11]). The 3-D exact solution is then obtained by rotating the 1-D viscous shock profile, which is initially located in the middle of the domain, by $45^{\circ}$ with respect to both $x-$ and $z-$axes. For this test problem, the Dirichlet boundary conditions are used at all physical boundaries and the Reynolds and Mach numbers are set equal to $10$ and $2.5$, respectively.

Moving grids for the viscous shock problem are generated by perturbing the vertices of a uniform Cartesian grid by $25\%$ of the grid spacing in a randomly chosen direction and mapping them onto vertices of a dynamic grid in the physical domain by using the following coordinate transformation:

[TABLE]

where the parameters $A_{1}$, $A_{2}$, $A_{3}$, and $\omega$ are set equal to $0.4$, $-0.2$, $0.3$, and $\pi$, respectively. Similar to the previous test problem, all moving grids for the viscous shock problem are straight-sided and the mesh velocities are calculated by differentiating Eq. (60).

The 3-D unsteady Navier-Stokes equations are integrated until $t=0.5$, which corresponds to an instant in time when the grid reaches the maximum deformation. Figure 3 shows density contours computed using the proposed entropy stable spectral collocation scheme for $p=4$ on the randomly perturbed $12\times 12\times 12$ grid at $t=0$ and $t=0.5$. The numerical solution obtained on this moving and deforming grid is free of spurious oscillations and accurately approximates the shock propagation in the domain.

To evaluate the accuracy of the new entropy stable schemes for the Navier-Stokes equations, we measure the $L_{2}$ and $L_{\infty}$ errors for the viscous shock problem on a sequence of the globally refined randomly perturbed moving grids. Since each grid is individually perturbed, the grids are not nested. As follows from our numerical results presented in Table 3, the proposed entropy stable schemes demonstrate in average $p$th-order convergence for all polynomial degrees considered. It should be noted that the error convergence is nonmonotonic because of the random perturbations introduced in dynamic grids used in this grid refinement study.

Similar to the Euler equations, a constant flow is an exact solution of the unsteady Navier-Stokes equations. To demonstrate that the new entropy stable spectral collocation schemes preserve a freestream viscous flow, we integrate the unsteady Navier-Stokes equations on the same randomly perturbed dynamic grids generated for the viscous shock problem. The $L_{\infty}$ errors of the entire solution vector computed with the new spectral collocation scheme given by Eq. (50) for $p=3,4,5$ on the $12\times 12\times 12$ randomly perturbed moving grid are presented in Fig. 4. Our numerical results show that the proposed entropy stable scheme preserves the constant solution with the machine accuracy for all polynomial degrees considered. Similar to the inviscid flow test case, there is a slight growth in the roundoff error accumulation as the polynomial degree increases.

7 Conclusions

New provably stable spectral collocation schemes of arbitrary order of accuracy are developed for the unsteady 3-D Euler and Navier-Stokes equations on dynamic unstructured grids. Using the arbitrary Lagrangian-Eulerian formulation, the governing equations are mapped onto a fixed reference system of coordinates. To provide the design order of accuracy, entropy stability, and conservation on moving and deforming grids, the skew-symmetric form of the Euler and Navier-Stokes equations are discretized by using high-order summation-by-parts spectral collocation operators. In addition to that, the entropy stability proof requires that the geometric conservation laws (GCL) must be satisfied exactly at the discrete level and the flux function has to satisfy an additional shuffle condition. To satisfy the GCL equations exactly, the metric coefficients are computed by using formulas that guarantee conservation, which are discretized by the same SBP operators employed to approximate the flux derivatives. Furthermore, the spatial derivative terms in the Euler and Navier-Stokes equations are calculated using the determinant-based Jacobian, while the time derivative term is calculated using the Jacobian obtained by integrating the corresponding GCL equation. This dual-Jacobian formulation satisfies all GCL equations simultaneously. Two new two-point entropy conservative flux functions satisfying the required shuffle condition are derived for the 3-D Euler and Navier-Stokes equations. Our numerical results show that the developed entropy stable spectral collocation schemes are fully conservative, design-order accurate, and provide freestream preservation for both the Euler and Navier-Stokes equations on highly deforming, moving unstructured grids.

Appendix A Entropy conservation of the semi-discrete form of the Euler equations

In this appendix, it is proven that the semi-discrete from (41) is entropy conservative for periodic domains. If the scheme is entropy conservative on a periodic domain, then it telescopes the entropy to the boundaries on non-periodic domains and the analysis of entropy conservation/dissipation reduces to the analysis of the weak imposition of boundary conditions.

The analysis follows in a one-to-one fashion the continuous analysis. The first step is to multiply by the entropy variables and discretely integrate in space; this is accomplished by multiplying (41) by $\bm{w}_{\kappa}^{\mathrm{T}}\overline{\mathsf{P}}$ which gives

[TABLE]

First, the term $I$ in (61) reduces to

[TABLE]

where $\bm{J}_{\kappa}$ is a vector containing the metric Jacobian at the nodes.

Using the symmetry of $\mathsf{F}_{x_{m}}^{sc}\left(\bm{u}_{\kappa},\bm{u}_{\kappa}\right)$, the SBP property $\mathsf{Q}+\mathsf{Q}=\mathsf{E}$, Theorem 4.4, the consistency of the derivative operator, $\mathsf{D}\bm{1}=\bm{0}\rightarrow\bm{1}^{\mathrm{T}}\mathsf{Q}=\bm{1}^{\mathrm{T}}\mathsf{E}$, and the relation, $\varphi=\bm{\mathit{W}}^{\mathrm{T}}\bm{\mathit{U}}-\mathcal{S}$, the term $II$ reduces to

[TABLE]

Using similar relations as for term $II$, term $III$ reduces to

[TABLE]

Inserting (62), (63), (64) into (61) and rearranging gives

[TABLE]

Ignoring for now the terms on the right-hand side of (65), since the term

[TABLE]

and is therefore mimetic of the appropriate continuous term, this implies that

[TABLE]

which are nothing more than discretizations of the continuous GCL (28) and (29).

Assuming the above hold,

[TABLE]

Next, the term $IV$ reduces to

[TABLE]

where the fact that $\overline{\mathsf{E}}_{\xi_{l}}$ is diagonal has been used throughout. Next, the term $V$ reduces to

[TABLE]

Substituting (67) and (68) into (66) and rearranging and summing over all elements gives

[TABLE]

The coupling terms each have counter-parts from abutting elements, thus, for example

[TABLE]

where Theorem 4.4 has been used. Moreover, advantage has been taken of 1) the fact that terms like $\mathsf{R}_{\alpha_{l}}^{\mathrm{T}}\mathsf{P}_{\perp\xi_{l}}\mathsf{R}_{\beta_{l}}$ are diagonal and therefore commute with the diagonal matrices containing the metric terms, 2) such terms only pick of metrics at the interface, and 3) the metric terms match at the interface.

The result of the above analysis is that the surface terms in (69) cancel out at interior elements; thus, on a periodic domain

[TABLE]

For non-periodic domains, the remaining terms are surface terms at the boundary of the domain and appropriate entropy conservative/stable SATs need to be employed.

Acknowledgments

N. Yamaleev and J. Lou gratefully acknowledge the support provided by Army Research Office through grant W911NF-17-1-0443. We would also like to thank Old Dominion University for the use of the Turing High Performance Computing cluster.