# 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 spectral collocation schemes developed for the one-dimensional compressible Navier-Stokes equations in [1,2] to three spatial dimensions. The new high-order schemes are provably L2 stable, design-order accurate for smooth solutions, and guarantee the pointwise positivity of thermodynamic variables for 3-D compressible viscous flows. Similar to the 1-D counterpart, the proposed schemes for the 3-D Navier-Stokes equations are constructed by using a flux-limiting technique that combines a positivity-violating entropy stable method of arbitrary order of accuracy and a novel first-order positivity-preserving entropy stable finite volume-type scheme discretized on the same Legendre-Gauss-Lobatto grid points used for constructing the high-order discrete operators. The positivity preservation and excellent discontinuity-capturing properties are achieved by adding an artificial dissipation in the form of the low- and high-order Brenner-Navier-Stokes diffusion operators. To our knowledge, this is the first family of positivity-preserving, entropy stable schemes of arbitrary order of accuracy for the 3-D compressible Navier-Stokes equations.

## Authors

• 2 publications
• 2 publications
11/05/2021

### 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...
01/27/2022

### A positivity preserving strategy for entropy stable discontinuous Galerkin discretizations of the compressible Euler and Navier-Stokes equations

High-order entropy-stable discontinuous Galerkin methods for the compres...
03/23/2020

### High-order accurate entropy stable finite difference schemes for the shallow water magnetohydrodynamics

This paper develops the high-order accurate entropy stable (ES) finite d...
03/19/2020

### Fully-Discrete Explicit Locally Entropy-Stable Schemes for the Compressible Euler and Navier-Stokes Equations

Recently, relaxation methods have been developed to guarantee the preser...
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...
06/20/2022

### Monolithic parabolic regularization of the MHD equations and entropy principles

We show at the PDE level that the monolithic parabolic regularization of...
##### 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

A new family of high-order entropy stable spectral collocation schemes that provide pointwise positivity of thermodynamic variables for the 1-D compressible Navier-Stokes equations has recently been introduced by the authors of the present paper in UY_1Dlow ; UY_1Dhigh . The entropy stability is achieved by using summation-by-parts operators and entropy consistent fluxes for discretizing both the inviscid and viscous terms of the symmetrized form of the compressible Navier-Stokes equations. We have proven in UY_1Dhigh that the new high-order flux-limiting schemes are both pointwise positivity-preserving and stable for the Navier-Stokes equations in one spatial dimension.

Herein, we generalize and extend the 1-D positivity-preserving entropy stable methodology of UY_1Dlow ; UY_1Dhigh to the three-dimensional compressible Navier-Stokes equations on static unstructured hexahedral grids. Similar to the 1-D high-order schemes in UY_1Dhigh , the new schemes for the 3-D compressible Navier-Stokes equations are constructed by combining a positivity-violating entropy stable method of arbitrary order of accuracy and a novel first-order positivity-preserving entropy stable method developed in the companion paper UY_3Dlow . In contrast to positivity–preserving methods developed in Svard ; GHKL ; GMPT , which are at most 2nd-order accurate, the proposed schemes provide an arbitrary order of accuracy for sufficiently smooth solutions of the 3-D compressible Navier-Stokes equations. Unlike the positivity-preserving high-order discontinuous Galerkin (DG) method developed for the Navier-Stokes equations in Zhang , the new high-order schemes guarantee not only the so-called weak positivity of thermodynamic variables, but also the pointwise positivity at individual collocation points that are directly used for approximation of the governing equations. Furthermore, the positivity-preserving DG scheme developed in Zhang imposes very severe constraints on the time step which is about an order of magnitude less than that of the baseline method for high-order polynomial bases. Note that the actual time step constraint in Zhang is much stiffer, because the lower bound on the artificial viscosity coefficient required for positivity may grow dramatically, as the velocity gradients increase.

Another distinctive feature of the proposed methodology is that the new high-order positivity-preserving schemes satisfy the discrete entropy inequality, thus facilitating a rigorous -stability proof for the symmetric form of the discretized Navier-Stokes equations. To our knowledge, this is the first family of high-order schemes that provide both the pointwise positivity of thermodynamic variables and entropy stability for the 3-D compressible Navier-Stokes equations.

## 2 The 3-D Navier-Stokes and Brenner-Navier-Stokes equations

The 3-D compressible Navier-Stokes equations in curvilinear coordinates are written in conservation law form as follows:

 (1)

where , 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,mUi−κ∂T∂xm]T. (2)

In the above equations,

is a vector of the conservative variables,

is the pressure, and is the specific total enthalpy. The viscous stresses in Eq. (2) are given by

 τ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. The governing equations (1) are subject to boundary conditions that are assumed to satisfy the entropy inequality.

It is well known that there are no theoretical results on positivity of thermodynamic variables for the compressible Navier-Stokes equations. To overcome this problem, 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 Brenner-Navier-Stokes equations are given by

 ∂U∂t+3∑m=1∂Fxm∂xm=3∑m=1∂F(B)xm∂xm,∀(x1,x2,x3)∈Ω,t≥0, (4)

where is the volume diffusivity and the viscous fluxes, , are defined as

 F(B)xm=F(v)xm+σ∂ρ∂xm[1VE]⊤. (5)

Despite that Eqs. (1) and (4) are very similar to each other, the Brenner-Navier-Stokes equations possess some remarkable properties that are not available for the Navier-Stokes equations. In contrast to Eq. (1), the Brenner-Navier-Stokes equations guarantee existence of a weak solution and uniqueness of a strong solution, ensure global-in-time positivity of the density and temperature, satisfy a large class of entropy inequalities, and is compatible with a minimum entropy principle FV ; GUERPOP2014 . We rely on these properties of the Brenner-Navier-Stokes equations and regularize the Navier-Stokes equations by adding the following dissipation term to Eq. (1):

where the the artificial dissipation flux can be obtained from the viscous flux of the Brenner-Navier-Stokes equations, , by setting , , and . The coefficient is an artificial viscosity and and are positive tunable coefficients, which are set equal to and for all test problems considered in this paper.

A necessary condition for selecting a unique, physically relevant solution among possibly many weak solutions of Eqs. (1) and (6) is the entropy inequality. Both the Navier-Stokes and regularized Navier-Stokes equations are equipped with the same convex scalar entropy function and entropy flux , where is the thermodynamic entropy. It can be shown that the following inequality holds for Eqs. (1) and (6) assuming the corresponding boundary conditions are entropy stable (e.g., see GCLSS2019 ):

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

Along with the entropy inequality given by Eq. (7), the regularized Navier-Stokes equations (Eq. (6)) preserve some other key properties of the Brenner-Navier-Stokes equations including the positivity of thermodynamic variables. Herein, we propose to develop new numerical schemes that mimic these properties of the 3-D regularized Navier-Stokes equations at the discrete level.

## 3 3-D high-order SBP operators

### 3.1 High–order Diagonal-Norm Summation-by-parts Operators

The derivatives in Eq. (6) are discretized by spectral collocation operators that satisfy the summation-by-parts (SBP) property CFNF . In the one-dimensional setting, this mimetic property is achieved by approximating the first derivative with a discrete operator, , and using local mass and stiffness matrices satisfying the following properties:

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

Only diagonal-norm SBP operators are considered herein, which is critical for proving the entropy inequality at the discrete level CFNF .

In one spatial dimension, the physical domain is divided into non-overlapping elements , so that . The discrete solution inside each element is defined on Legendre-Gauss-Lobatto (LGL) points, . These local points are referred to as solution points. 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. (8), 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 , these discrete SBP derivative operators can be recast in the following telescopic flux form:

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

where is a th-order flux vector defined at the flux points.

Using tensor product arithmetic, the above one-dimensional SBP operators naturally extend to two and three spatial dimensions. The multidimensional tensor product operators are defined as

 Dξ1=(DN⊗IN⊗IN⊗I5), Pξ1=(PN⊗IN⊗IN⊗I5), (10) 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 . We also use the following notation and where is the scalar th diagonal entry of .

## 4 Baseline 3-D high-order spectral collocation schemes

With the 3-D SBP operators discussed in the previous section, a baseline 3-D semi-discrete spectral collocation scheme of arbitrary order of accuracy for the 3-D Navier-Stokes equations (Eq. (1)) can be written as follows:

 ^Ut+3∑l=1P−1ξlΔξl^¯fl−Dξl^f(v)l=3∑l=1P−1ξl^g(BC)l+P−1ξl^g(Int)l (11)

where , and are boundary and interface penalty terms.

The contravariant inviscid fluxes, , defined at flux points are given by

 ^¯fm(→ξ¯i)=N∑j=i+1i∑l=12ql,j¯f(S)(U(→ξl),U(→ξj))^→am(→ξl)+^→am(→ξj)2 for1≤i≤N−1,^¯fm(→ξ¯i)=¯f(S)(U(→ξ¯i),U(→ξ¯i))^→am(→ξ¯i) for¯i∈{0,N}, (12)

where and is a two-point, consistent, entropy conservative flux that satisfies

 (w1−w2)⊤¯f(S)(U1,U2)=→ψ1−→ψ2 (13)

for any two admissible states and TAD2003 . For all test problems considered, we use the entropy conservative flux developed in Chand . In UY_3Dlow , we show that the proposed method for ensuring positivity is independent of a particular choice of .

In the above equation and hereafter, is a th-order discrete approximation of at the solution point , which is constructed such that it satisfied the geometric conservation law (GCL) equations. These approximations are not unique and the specific formulas used for in the present work can be found elsewhere (e.g. (TLGCL, )).

The contravariant viscous fluxes, , are constructed as follows:

 ^f(v)l =3∑m=1[^alm]f(v)xm,f(v)xm=3∑j=1[c(v)m,j]Θxj. (14)

For each , is a block-diagonal matrix with blocks, such that , and , i.e., the full viscous tensor is symmetric positive semi-definite (SPSD).

The gradient of the entropy variables, , is discretized by using an approach that closely resembles the local discontinuous Galerkin (LDG) method developed in CS , which can be written as

 Θxj=3∑l=1[^alj][J−1](Dξlw+P−1ξl^gΘl)^gΘ1(→ξijk)=12(δ1iΔ1w(→ξi−1jk)+δNiΔ1w(→ξijk))Δ1w(→ξijk)=w(→ξi+1jk)−w(→ξijk), (15)

where is the Kronecker delta. Note that similar discretizations are used in each computational direction.

In CFNF , it has been proven that the baseline high-order spectral collocation scheme given by Eqs. (1115) satisfies the discrete entropy inequality if the corresponding boundary conditions are entropy stable. However, entropy stability alone does not guarantee the positivity of thermodynamic variables, if strong discontinuities are present in the domain. One of the key objectives of this paper is to modify the baseline scheme to address this shortcoming.

## 5 Baseline 3-D spectral collocation scheme with high-order artificial dissipation

To control the amount of entropy production in regions where the discrete solution is under-resolved, we generalize the method developed in UY to three spatial dimensions and add artificial dissipation in the form of the Brenner diffusion operator to the baseline high-order scheme (Eq. (11)) presented in the foregoing section. The baseline 3-D spectral collocation scheme with the high-order artificial dissipation is given by

where . The high-order artificial dissipation terms, and , are discretized similarly to the viscous terms of Eq. (11), as discussed in Section 4. The Brenner fluxes, , are constructed as follows:

where , are block–diagonal matrices with blocks, , and , i.e., the full artificial dissipation tensor is symmetric positive semi-definite (SPSD).

To ensure consistency, maintain design-order accuracy for smooth resolved solutions, and control the amount of dissipation added in regions where the solution is under-resolved or discontinuous, we use the artificial viscosity, , that is described in Section 6. The mass and heat viscosity at each solution point are set as , and (see Section 2).

The high-order spectral collocation scheme given by Eq. (16) is conservative and stable in the entropy sense. The conservation follows immediately from the telescopic flux form of the inviscid terms and the SBP form of the viscous and artificial dissipation terms. The entropy stability of the discretized Navier-Stokes terms in Eq. (16) is proven in CFNF . The entropy dissipation properties of the artificial dissipation terms follow immediately form Eq. (17) and the symmetric positive semi-definiteness of the artificial viscous tensor , .

## 6 Artificial Viscosity

The scalar artificial viscosity, , is used for both the high- and low-order artificial dissipation operators. Details on how the artificial viscosity coefficient is constructed are presented in UY_3Dlow . Herein, we only briefly outline its key elements. The artificial viscosity coefficient is constructed based on the finite element residual of the entropy equation and the physical properties of the fluid. In the -th grid element, is defined as follows:

where where is a sensor function () and is the magnitude of the artificial viscosity in the -th grid element.

To detect grid elements where the solution loses its regularity or becomes under-resolved, the sensor is constructed as a function of the finite element residual of the entropy equation, which is given by

 Snk={Snk0,if Snk0≥max(0.2,δ),0,otherwiseSnk0=max(rk)max(1,p−1p−1.5), (18)

where is the polynomial order and is a pointwise normalized entropy residual. To take into account the physics of a problem, we also augment the entropy residual-based sensor with compression and pressure gradient sensors. These sensors are introduced to identify those regions where the amount of artificial viscosity can be reduced without sacrificing the solution accuracy. We refer the reader to UY_3Dlow for further details.

In each element, the upper bound of the artificial viscosity, , is set to be proportional to the maximum value of local velocity and pressure jumps between neighboring solution points UY_3Dlow . The result is that we minimize the amount of artificial dissipation at contact discontinuities and make proportional to the discontinuity strength, such that the velocity and pressure jumps act as a limiter, if spurious oscillations are present in the solution.

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.

## 7 High-order positivity–preserving flux-limiting scheme

Following an approach developed in UY_1Dhigh , we construct a new high-order positivity–preserving flux-limiting scheme for the 3-D Navier-Stokes equations by combining the corresponding positivity-violating high-order spectral collocation scheme (Eq. (11)) and the first-order positivity–preserving finite volume scheme presented in the companion paper UY_3Dlow . This methodology is presented next.

### 7.1 Positivity

We first consider the 1st-order explicit Euler approximation of the time derivative term in Eq. (1), so that on a given element

 ^Un+1p=^Un+τ(d^Udt)p,^Un+11=^Un+τ(d^Udt)1,

where and are th-order and first-order numerical solutions defined on the same Legendre-Gauss-Lobatto (LGL) grid elements with the same high-order metric terms. In the above equation, is obtained by the first-order positivity-preserving entropy stable scheme presented in UY_3Dlow . Therefore, at every th solution point of each element and , where is the internal energy associated with the 1st-order solution .

To combine the 1st- and th-order schemes, we use the flux-limiting technique developed in UY_1Dhigh , which is in fact equivalent to limiting the low- and high-order solution vectors of the conservative variables:

 ^Un+1(θf)=^Un+τ[(1−θf)(d^Udt)1+θf(d^Udt)p]=(1−θf)^Un+11+θf^Un+1p=^Un+11+θf[^Un+1p−^Un+11], (19)

where the flux limiter , , is a constant on a given high-order element.

At each solution point, local lower bounds of density and internal energy are defined as follows:

 ϵρi=(ρ1)n+1iℵ,ϵIEi=IE((^U1)n+1i)ℵ, (20)

where , , is a function that is bounded from below by a small positive number (e.g., ), which approaches to its lower bound if the solution is smooth and goes to 1 if the solution loses its regularity. In the present analysis, is defined as follows:

 ℵk=max(10−8,Lk),Lk=Snkmaxi(|ΔP|2PA), (21)

where is the residual-based sensor given by Eq. (18) and is one half of the maximum relative two–point pressure jump (including jumps at the interfaces) on the th element. Note that and because .

We now prove the following two lemmas.

###### Lemma 1.

For every -th solution point, define a set

 Hρi={θf∈[0,1]|ρn+1i(θf)≥ϵρi}.

Then, the set can be written as where . Furthermore, the following statements hold: 1) if , then and 2) if , then .

###### Proof.

This follows directly from the fact that given by Eq. (19) is a linear equation in the variable with . ∎

###### Lemma 2.

For every -th solution point, define a set

 HIEi={θf∈Hρi|IE(^Un+1i(θf))≥ϵIEi},

where is defined in Lemma 1. Then, the set can be written as where . Furthermore, the following statements hold: 1) if , then and 2) if , then .

###### Proof.

For each th solution point, if , then we set . Assume that there is a solution point such that . Since , is a continuous function with respect to for . Taking into account that and , it follows by the intermediate value theorem that there exists such that . Let (note that there is only one such that ). Now we show that for all , we have . By definition of , . For , we have

 ^Un+1i(θf) =(1−θf)(^U1)n+1i+θf(^Up)n+1i (22) =θfθIEi[θ% IEi((^Up)n+1i−(^U1)n+1i)+(^U1)n+1i]+⎛⎝1−θfθIE% i⎞⎠(^U1)n+1i =θfθIEi^Un+1i(θIEi)+⎛⎝1−θfθIEi⎞⎠(^U1)n+1i.

Hence, due to the concavity of internal energy

 IE(^Un+1i(θf))≥θfθIEiIE(^Un+1i(θIEi))+(1−θfθIEi)%IE((^U1)n+1i)>θfθIEiϵIEi+(1−θfθIEi)ϵIEi=ϵIEi. (23)

###### Remark 1.

Note that in Lemma 2 can readily be found by solving the quadratic equation for internal energy, which is analogous to the one presented in the companion paper UY_3Dlow .

For a given element, we define . By construction, and for every solution point on the element. The solution at the th time level is set equal to , which preserves the pointwise positivity of both density and internal energy.

###### Remark 2.

The above limiting is not immediately conservative for general and . We refer the reader to Section 7.3 which presents an implementation of this limiting procedure in a way that preserves conservation.

### 7.2 Design order of accuracy

In this section, we prove that the proposed limiting scheme is design-order accurate for smooth solutions and sufficient grid resolutions. Without loss of generality, we assume that the grid resolution depends on a single parameter , such that all element edges are directly proportional to . In this section, denotes the Euclidean norm. Let be the smooth exact solution at the th solution point at . For each solution point, we define a local admissible set

 Aϵi={ui=[ρρ→VρE]⊤ | IE% (ui)≥ϵIEi ,ρi≥ϵρi}

and assume that . Note that and are positive user-defined parameters that can be made arbitrarily small by selecting a sufficiently small value of the parameter for a given element. In the present analysis, , which is given by Eq. (21), is set such that it becomes smaller when the regularity of the numerical solution increases. We also assume that the solution is sufficiently smooth, so that , as .

Let us show that for all solution points. If on a given element, then , and the result follows.

We now assume that . In this case, to prove the consistency of the limiting procedure, it is sufficient to show that . Indeed, if , then for every solution point we have

 ∥^Un+1i(θIE)−^Uexi(tn+1)∥≤(1−θIE)∥(^U1)n+1i−^Uexi(tn+1)∥  +θIE∥(^Up)n+1i−^Uexi(tn+1)∥=(1−θIE)O(hx)+θIEO((hx)p)=O((hx)p). (24)

To prove that , it is sufficient to show that if (which is only possible if ), then . Assume that at the th solution point . Since , we only have to consider the following two cases: 1) and , 2) .

Case 1. For , the following inequalities hold , which lead to . From Lemma 1 it follows that satisfies

 ρn+1i(θρi)=(ρ1)n+1i+θρi((ρp)n+1i−(ρ1)n+1i)=ϵρi. (25)

Thus,

 1−θρi=ϵρi−(ρp)n+1i(ρ1)n+1i−(ρp)n+1i=O((hx)p)O(hx)=O((hx)p−1). (26)

Taking into account that , we also have . Using and Eq. (26) yield

 ∥^Un+1i(θρi)−^Uexi(tn+1)∥≤∥^Un+1i(θρi)−(^Up)n+1i∥+∥(^Up)n+1i−^Uexi(tn+1)∥=(1−θρi)∥((^U1)n+1i−(^Up)