# First-order positivity-preserving entropy stable spectral collocation scheme for the 3-D compressible Navier-Stokes equations

In this paper, we extend the positivity-preserving, entropy stable first-order finite volume-type scheme developed for the one-dimensional compressible Navier-Stokes equations in [1] to three spatial dimensions. The new first-order scheme is provably entropy stable, design-order accurate for smooth solutions, and guarantees the pointwise positivity of thermodynamic variables for 3-D compressible viscous flows. Similar to the 1-D counterpart, the proposed scheme for the 3-D Navier-Stokes equations is discretized on Legendre-Gauss-Lobatto grids used for high-order spectral collocation methods. The positivity of density is achieved by adding an artificial dissipation in the form of the first-order Brenner-Navier-Stokes diffusion operator. Another distinctive feature of the proposed scheme is that the Navier-Stokes viscous terms are discretized by high-order spectral collocation summation-by-parts operators. To eliminate time step stiffness caused by the high-order approximation of the viscous terms, the velocity and temperature limiters developed for the 1-D compressible Navier-Stokes equations in [1] are generalized to three spatial dimensions. These limiters bound the magnitude of velocity and temperature gradients and preserve the entropy stability and positivity properties of the baseline scheme. Numerical results are presented to demonstrate design-order accuracy and positivity-preserving properties of the new first-order scheme for 2-D and 3-D inviscid and viscous flows with strong shocks and contact discontinuities.

## Authors

• 2 publications
• 2 publications
11/16/2021

### High-order Positivity-preserving L2-stable Spectral Collocation Schemes for the 3-D compressible Navier-Stokes equations

This paper extends a new class of positivity-preserving, entropy stable ...
09/30/2020

### Stability analysis of a singlerate and multirate predictor-corrector scheme for overlapping grids

We use matrix stability analysis for a singlerate and multirate predicto...
11/18/2019

### Entropy stable discontinuous Galerkin approximation for the Relativistic Hydrodynamic Equations

This paper presents the higher-order discontinuous Galerkin entropy stab...
09/13/2020

### Second-order invariant domain preserving approximation of the compressible Navier–Stokes equations

We present a fully discrete approximation technique for the compressible...
05/17/2020

### An energy-dissipative level-set method for the incompressible two-phase Navier-Stokes equations with surface tension using functional entropy variables

This paper presents the first energy-dissipative level-set method for th...
11/21/2019

### On the robustness and performance of entropy stable discontinuous collocation methods for the compressible Navie-Stokes equations

In computational fluid dynamics, the demand for increasingly multidiscip...
11/21/2019

### On the robustness and performance of entropy stable discontinuous collocation methods for the compressible Navier-Stokes equations

In computational fluid dynamics, the demand for increasingly multidiscip...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

We have recently developed a novel first–order entropy stable finite volume scheme that provides pointwise positivity of thermodynamic variables for the 1-D compressible Navier-Stokes equations on Legendre-Gauss-Lobatto (LGL) grids used for high-order spectral collocation methods UY_1Dlow . The positivity preservation and entropy stability properties are achieved by introducing the first-order artificial dissipation operator that mimics the corresponding diffusion operator of the Brenner-Navier-Stokes equations Brenner1 . It has been proven that the proposed first-order scheme is conservative, entropy stable, and positivity preserving for the 1-D compressible Euler and Navier-Stokes equations. This positivity-preserving methodology has recently been extended to entropy stable spectral collocation schemes of arbitrary order of accuracy for the 1-D Navier-Stokes equations in the companion paper UY_1Dhigh . Herein, we generalize and extend this 1-D positivity-preserving entropy stable finite volume scheme developed in UY_1Dlow to the three-dimensional compressible Navier-Stokes equations on fully unstructured static hexahedral grids.

There are very few papers available in the literature on positivity-preserving methods for the compressible Navier-Stokes equations especially in three spatial dimensions. A first-order positivity-preserving finite difference scheme for the 3-D compressible Navier-Stokes equations on Cartesian uniform grids has been developed in Svard . The positivity proof in this paper is based on some memetic properties of the 1st-order finite difference operators on uniform grids, which are not available for curvilinear or unstructured grids. In GHKL , an implicit first-order positivity-preserving scheme is constructed for the compressible Navier-Stokes equations on staggered grids. This scheme is unconditionally stable and solves the internal energy equation instead of the equation for conservation of total energy. Another positivity-preserving scheme for the compressible Navier-Stokes equations has been proposed in GMPT . This approach relies on the invariant domain preserving approximation of the Euler equations and the Strang’s operator splitting technique that is at most 2nd-order accurate. Note that no entropy stability proof is currently available for the numerical schemes developed in GHKL ; GMPT . Recently, a high-order positivity-preserving discontinuous Galerkin (DG) scheme for the compressible Navier-Stokes equations is presented in Zhang . This explicit in time method provides only so-called weak positivity of density and pressure and imposes very severe constraints on the time step. Note that the actual time step constraint is much stiffer, because the lower bound on the artificial viscosity coefficient, which is required for density and pressure positivity, may grow dramatically, as the velocity gradients increase.

The main objective of the present paper is to construct a first-order positivity-preserving entropy stable scheme defined on LGL grids which are used for high-order spectral collocation methods (e.g., CFNF ; CFNPSY ). To provide the positivity of thermodynamic variables, we construct new first-order artificial dissipation operators that are based on the Brenner regularization of the compressible Navier-Stokes equations UY_1Dlow ; UY . In contrast to the existing positivity-preserving schemes that rely on monotonicity properties of the Rusanov-type dissipation, the proposed method minimizes the amount of artificial dissipation required for pointwise positivity of density and temperature and uses novel velocity and temperature limiters to eliminate the time step stiffness for viscous flows with strong shock waves and contact discontinuities.

## 2 3-D regularized Navier-Stokes equations

We consider the 3-D compressible Navier-Stokes equations in curvilinear coordinates , which are written in conservation law form as follows:

 (1)

where

is a vector of the conservative variables and

, and are the inviscid and viscous fluxes associated with the Cartesian coordinates , which are given by

 F(v)xm=[0,τ1,m,τ2,m,τ3,m,3∑i=1τi,mVi−κ∂T∂xm]T. (2)

The viscous stresses in Eq. (2) are defined as follows:

 τi,j=μ(∂Vi∂xj+∂Vj∂xi−δi,j233∑n=1∂Vn∂xn), (3)

where is the dynamic viscosity, is the thermal conductivity, and is the Kronecker delta. To close the Navier-Stokes equations, the following constituent relations are used:

 H=cPT+12VTV,P=ρRgT,Rg=RuMw,

where is the temperature, is the universal gas constant, is the molecular weight of the gas, and is the specific heat capacity at constant pressure. Equation (1) is subject to boundary conditions that are assumed to satisfy the entropy inequality. Note that only static curvilinear unstructured grids are considered in the present analysis.

Since the Navier-Stokes equations have no inherent mechanism for producing admissible solutions with positive thermodynamic variables, we regularize Eq. (1) by adding artificial dissipation in the form of the diffusion operator of the Brenner-Navier-Stokes equations introduced in Brenner1 . The regularized Navier-Stokes equations in curvilinear coordinates are given by

where is an artificial viscosity and and are positive tunable coefficients. Equations (1) and (4) are derived by using the following identies:

 3∑l=1∂∂ξl(J∂ξl∂xm)=0,m=1,2,3, (5)

which are called the geometric conservation laws (GCL) TLGCL . Though, the GCL equations (5) are satisfied exactly at the continuous level, this is not the case at the discrete level TLGCL . A discussion on how the corresponding metric coefficients should be discretized to satisfy the GCL equation is presented elsewhere (e.g., see (TLGCL, ; GCLSS2019, )).

Both the Navier-Stokes and regularized Navier-Stokes equations are equipped with the same scalar entropy function and the corresponding entropy flux , where is the thermodynamic entropy and is the velocity vector. The mathematical entropy, , is convex and its Hessian matrix, , is positive definite provided that and , thus yielding a one-to-one mapping from the conservative to entropy variables, which are defined by . For entropy stable boundary conditions, the following entropy inequality holds for both Eqs. (1) and (4):

 ∫^Ω∂(JS)∂τd^Ω=ddτ∫^ΩJSd^Ω≤0. (6)

Note that the entropy inequality (6) 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. In contrast to the conventional Navier-Stokes equations (1), the regularized Navier-Stokes equations (4) provide global-in-time positivity of thermodynamic variables FV . Herein, we develop a new numerical scheme that replicates this positivity property of the regularized Navier-Stokes equations (4) at the discrete level.

## 3 Discrete operators

Similar to the first-order scheme developed in one spatial dimension in UY_1Dlow , the proposed first-order scheme for the 3-D Navier-Stokes equations is discretized on the same Legendre-Gauss-Lobatto (LGL) points used for high-order spectral collocation operators. In the one-dimensional setting, the physical domain is divided into non-overlapping elements with nonuniformly distributed points, so that . The discrete solution inside each element is defined on the LGL points, , associated with the Lagrange polynomial basis of degree . These local points are referred to as solution points.

The derivatives of the viscous fluxes in (4) are discretized by high-order spectral collocation operators that satisfy the summation-by-parts (SBP) property CFNF ; SN1 . This mimetic property is achieved by approximating the first derivative by using the following discrete operator, :

 D=P−1Q,P=P⊤,v⊤Pv>0,∀v≠0,Q=B−Q⊤,B=diag(−1,0,…,0,1), (7)

where is a diagonal mass matrix and is a stiffness matrix. Only diagonal-norm SBP operators are considered herein, which is critical for proving the entropy inequality at the discrete level CFNF .

Along with the solution points, we also define a set of intermediate points prescribing bounding control volumes around each solution point. These points referred to as flux points form a complementary grid whose spacing is precisely equal to the diagonal elements of the positive definite matrix in Eq. (7), i.e., , where is a vector of flux points, , and is an matrix corresponding to the two-point backward difference operator CFNF ; YC3 . As has been proven in FCNYS , all discrete SBP derivative operators can be recast into the following telescopic flux form:

 P−1Qf=P−1Δ¯f, (8)

where is a high-order flux vector defined at the flux points. Note that the above telescopic flux form also satisfies the generalized SBP property FCNYS .

For unstructured hexahedral grids, the above one-dimensional SBP operators defined on each grid element naturally extend to three spatial dimensions by using tensor product arithmetic. The multidimensional tensor product operators are defined as follows:

 Dξ1 =(DN⊗IN⊗IN⊗I5), Pξ1 =(PN⊗IN⊗IN⊗I5), (9) Pξ1,ξ2 =(PN⊗PN⊗IN⊗I5), P =(PN⊗PN⊗PN⊗I5), ˆP =(PN⊗PN⊗PN), P⊥,ξ1 =(IN⊗PN⊗PN⊗I5),

with similar definitions for other directions and operators , and . Also, the following notation is used hereafter: and where is the scalar -th diagonal entry of .

The metric coefficients are also discretized by using the high-order SBP operators (Eq. (7))such that the GCL equations given by Eq. (5) are satisfied exactly at the discrete level TLGCL . The discrete approximation of the scalar metric coefficient at the solution point is denoted . The block diagonal matrix contains blocks with entries where

is the identity matrix of size 5. The specific formulas for

are recorded elsewhere (e.g., see TLGCL ; GCLSS2019 ). Note that is continuous at element interfaces and satisfies the following GCL equations exactly:

 3∑l=1Dξl[^alm]15=05,m=1,2,3. (10)

## 4 Artificial Viscosity

The artificial viscosity coefficient, , in Eq. (4) is constructed as a function of the entropy equation residual and consists of three major components: 1) entropy residual, 2) sensor functions, and 3) upper bound of the artificial viscosity coefficient. In this section, we make use of the following globally defined parameters: and , where is the spatial dimensionality, is the total number of elements used, and is the volume on the -th element.

### 4.1 Entropy residual

We directly generalize the 1-D entropy residual presented in UY to three spatial dimensions. The finite element residual of the entropy equation on the -th element is approximated as follows:

 R=3∑l=1(P−1ξlΔξl^¯fl−Dξl^f(v)l−^gl)−3∑l=1(−ˆDξl^Fl+ˆDξl(w⊤^f(v)l)−(Θξl)⊤^f(v)l). (11)

where , , , , and are evaluated by using the 1st-order discrete solution defined on the LGL solution points. Note that no high-order discrete solution is required to compute the entropy residual given by Eq. (11). We refer the reader to UY which discusses the key advantages of approximating the entropy residual by Eq. (11).

### 4.2 Sensor functions

We detect regions where the discrete solution is under-resolved and minimize the amount of artificial dissipation in regions where the solution is smooth by using the following sensor functions: 1) an entropy residual-based sensor, 2) a compression sensor, and 3) a pressure gradient sensor.

To construct the residual-based sensor, we first define an auxiliary pointwise sensor as follows:

 r(→ξijl)=R(→ξijl)J(→ξijl)max(R(→ξijl)J(→ξijl),η(→ξijl)), (12)

where

 (13)

where is the magnitude of a vector. In the above equation, where is the -th component of at and is the -order approximation of the gradient of the entropy variables CFNPSY , and is the speed of sound. The entropy residual sensor for the -th element, , is then defined as follows:

 Snk0=max(rk)max(1,p−1p−1.5),Snk={Snk0,if Snk0≥max(0.2,δ)0,otherwise. (14)

To identify regions where the amount of artificial dissipation can be reduced without generating spurious oscillations, we augment the entropy residual sensor with the compression and pressure gradient sensors, which are used only if the residual sensor on the -th element.

To minimize the amount of artificial dissipation near expansion waves, the compression sensor, , on the -th element is defined by using the integral of the divergence of the velocity filed over this element

 Cnk=(Cnk0)barctan[a(Cnk0−Cn∗)]+π2arctan[a(1−Cn∗)]+π2,Cnk0=max(−JˆP(∇h⋅→V)JˆP∣∣(∇h⋅→V)∣∣+ϵ,0), (15)

where is the approximation of the divergence of the velocity vector, and , , and , are tunable parameters that are set to be , , and for all problems considered.

The pressure, , sensor aims to reduce the amount of artificial dissipation near weak shock waves and is defined as follows:

 Pnk=max(0,Pnk0Pnkd),Pnk0=−∑Ni,j,k=1PijkJ(→ξijk)→V(→ξijk)⋅(∇hP)(→ξijk),Pnkd=ϵ+∑Ni,j,k=1PijkJ(→ξijk)∥→V(→ξijk)∥∥(∇hP)(→ξijk)∥. (16)

### 4.3 Local reference grid spacing

On each -th element where , we define a reference grid spacing, , that is used to compute an upper bound of the artificial viscosity coefficient, . The technique presented in this section for calculating is used for all test problems considered in Section 7.

First, we compute the following array of reference lengths, , defined at each solution point of the -th element:

 Lk(→ξijl) =2[3∏m=1∥Dmx∥¯Em]→ξijl,¯Em(→ξijl)=Em(→ξijl)∑3n=1En(→ξijl), (17) Em(→ξijl) =∥∥ ∥∥[d1→Vd1x]Dmx∥Dmx∥∥∥ ∥∥→ξijl+ϵ,m=1,2,3,

where is a -order approximation of the tangential derivative in the direction and is a block diagonal matrix with entries

 [d1→Vd1x](→ξijl)=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣d1V1d1x1d1V1d1x2d1V1d1x3d1V2d1x1d1V2d1x2d2V2d1x3d1V3d1x1d1V3d1x2d3V3d1x3⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦→ξijl. (18)

Then, an average grid spacing on the -th element is evaluated as follows:

 ^hk=1⊤1ˆPLk1⊤1ˆP11. (19)

Finally, the local reference grid spacing, , is given by

 hk=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩[∏i∈Nkhvi]1|Nk|,if |Nk|>0,0,otherwisehvi=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩ll[∏j∈Ii^hj]1|Ii|,if |Ii|>00,otherwise, (20)

where is a set of indices of all elements that touch the -th global vertex and have a nonzero , and is a set of indices of all vertices that touch the -th element and have a nonzero .

### 4.4 Artificial viscosity coefficient

As in UY , the artificial viscosity for each element is constructed by first finding its physics-based upper bound, . The elementwise scalar function,

, is constructed such that it is proportional to two-point jumps in the maximum eigenvalue of the inviscid flux Jacobian. If

on the -th element, we set . If , then is given by

 μkmax=(hk)2p3(γ+1)32γzSnkmax1≤i,j,l≤N⎡⎣zPnk,Cnk√¯ρ3∑m=1(d1√γPd1xm)2+¯ρ⎛⎜ ⎜⎝P1,12  ⎷3∑n,m=1,n≠m(d1Vnd1xm)2+min(Ma,zCnk)∣∣∣3∑a=1d1Vad1xa∣∣∣⎞⎟ ⎟⎠⎤⎥ ⎥⎦→ξijl,¯ρ(→ξijl)=(ρ(→ξijl)ρ(→ξi+1jl)ρ(→ξi−1jl)ρ(→ξij+1l)ρ(→ξij−1l)ρ(→ξijl+1)ρ(→ξijl−1))17,Ma(→ξijl)=∥→V(→ξijl)∥c(→ξijl),zSnk=min(0.5,1.25(Snk−0.2))≥0,zCnk=P1,12(1−Cnk)+Cnk,zPnk,Cnk=min(P1,12(1−Pnk)+Pnk,zCnk), (21)

where is a reference grid spacing defined in Section 4.3, is the smallest distance between flux points of the 1-D computational element, and is the speed of sound at . Note that is constructed such that it achieves its maximum value at shocks.

The velocity gradient components in Eq. (21) are approximated by using the following two-point discretizations:

 d1Vd1ξm(→ξi)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩V(→ξi+1)−V(→ξi)ξmi+1−ξmi,if i≤p+12V(→ξi+1)−Va(→ξi−1)ξmi+1−ξmi−1,if p+12

for . At element interfaces, we also include velocity jumps formed by the collocated states of elements that share the entire face. The gradient of the term in Eq. (21) is calculated in a similar fashion.

The globally continuous artificial viscosity is then constructed by using the following smoothing procedure. At each element vertex, we from a unique vertex viscosity coefficient, , where contains indices of all elements that share the

-th grid vertex. After that, the globally continuous artificial viscosity is obtained by using the tri-linear interpolation of 8 vertex viscosities,

, of the given hexahedral element.

## 5 First-order positivity-preserving scheme

We now present a positivity-preserving, entropy stable first-order scheme for the regularized -D compressible Navier-Stokes equations (4). Hereafter, is used to denote the vector of primitive variables at the -th solution point. We also make substantial use of the logarithmic, harmonic, arithmetic, and geometric averages, which are denoted for quantities and by using the following subscript notation: , , , and , respectively. Note that the following inequalities hold for any , : .

### 5.1 First-order scheme

The first-order scheme on a given element is approximated on the same LGL points used for the high-order scheme. The first-order element treats solution points in a finite volume manner with the flux points acting as control volume edges and can be written as

where and , are first-order inviscid and artificial dissipation fluxes, , are the high-order physical fluxes, and are high-order penalties that are identical to those used in CFNPSY . The way how first-order inviscid and artificial dissipation fluxes are constructed is discussed next.

### 5.2 First-order inviscid term

The inviscid fluxes in Eq. (23) are represented as the sum of entropy conservative and entropy dissipative terms: . The exact form of is presented in Section 5.4.

The entropy conservative flux, , is defined as follows:

 (24)

where and is any two-point, entropy conservative inviscid flux. For any two admissible states and , this two–point flux satisfies the following condition TAD2003 :

 (w1−w2)⊤¯f(S)(u1,u2)=→ψ1−→ψ2. (25)

In the present analysis, we use the entropy conservative flux developed in Chand . Comparing Eq. (24) with the high-order entropy stable flux in CFNPSY , we note that they are equivalent at the element faces () and only differ at the interior points. The high-order interpolation of the metric terms to the flux points is sufficient for proving Lemma 1.

###### Lemma 1.

The inviscid flux given by Eq. (24) is freestream preserving and entropy conservative, so that the following equation holds:

 3∑l=1w⊤PP−1ξlΔξl^¯f(EC)l=3∑l=11⊤1ˆP⊥,ξlˆBξl^Fl. (26)

Hence, given by Eq. (24) has the same total entropy contribution on each element as the high-order entropy consistent flux in CFNPSY .

###### Proof.

To prove the freestream preservation, we show that the following equation holds for any admissible constant state, : Note that , for any two solution points and on this element. Let us evaluate at point on the element:

 [3∑l=1P−1ξlΔξl^¯f(EC)l](→ξijk)=f(u0)⎡⎣^¯→a1(→ξ¯ijk)−^¯→a1(→ξ¯¯¯¯¯¯i−1jk)Pi,i+^¯→a2(→ξi¯jk)−^¯→a2(→ξi¯¯¯¯¯¯¯j−1k)Pj,j+^¯→a3(→ξij¯k)−^¯→a3(→ξij¯¯¯¯¯¯¯k−1)Pk,k⎤⎦=f(u0)⎡⎢ ⎢⎣N∑n=1qi,n^→a1(→ξnjk)Pi,i+N∑n=1qj,n^→a2(→ξink)Pj,j+N∑n=1qk,n^→a3(→ξijn)Pk,k⎤⎥ ⎥⎦=[0,…,0]⊤, (27)

where the last equality follows from the fact that the metric coefficients satisfy the discrete GCL equation given by Eq. (10).

Let us now show that Eq. (26) holds.

 3∑l=1w⊤PP−1ξlΔξl^¯f(EC)l=3∑l=1w⊤P⊥,ξlΔξl^¯f(EC)l=3∑l=1ˆ1⊤1ˆP⊥,ξl[w⊤Δξl^¯f(EC)l]=3∑l=1ˆ1⊤1ˆP⊥,ξl[w⊤~Bξl^¯f(EC)l−w⊤~Δξl^¯f(EC)l]=3∑l=1ˆ1⊤1ˆP⊥,ξl[ˆBξl(^ψl+^Fl)−w⊤~Δξl^¯f(EC)l]. (28)

Hence, it only remains to show that .

 (29)

where . The last equality in the above equation again follows from the discrete GCL, Eq. (10). ∎