∎

11institutetext: M. Sayyari 22institutetext: 22email:
<EMAIL_ADDRESS>33institutetext: L. Dalcin 44institutetext:
44email<EMAIL_ADDRESS>55institutetext: M. Parsani
66institutetext: 66email<EMAIL_ADDRESS>
King Abdullah University of Science and Technology (KAUST), Computer
Electrical and Mathematical Science and Engineering Division (CEMSE), Extreme
Computing Research Center (ECRC), 23955-6900, Thuwal, Saudi Arabia

# Development and analysis of entropy stable no-slip wall boundary conditions
for the Eulerian model for viscous and heat conducting compressible flows

Mohammed Sayyari  Lisandro Dalcin  Matteo Parsani

(Received: date / Accepted: date)

###### Abstract

Nonlinear entropy stability analysis is used to derive entropy stable no-slip
wall boundary conditions for the Eulerian model proposed by Svärd (_Physica A:
Statistical Mechanics and its Applications, 2018_). and its spatial
discretization based on entropy stable collocated discontinuous Galerkin
operators with the summation-by-parts property for unstructured grids. A set
of viscous test cases of increasing complexity are simulated using both the
Eulerian and the classic compressible Navier–Stokes models. The numerical
results obtained with the two models are compared, and differences and
similarities are then highlighted.

††journal: SN PDEs with Applications

##  1 Introduction

The classical compressible Navier–Stokes (CNS) equations can be derived based
on the material (Lagrangian) derivative formulation Schneiderbauer and Krieger
(2013). In the Lagrangian sense, diffusion between gas pockets is non-
existent, and thus, the continuity equation is hyperbolic. On the other hand,
in the Eulerian model of Svärd (2018), air molecules diffuse into other parts
of the domain, and thus, the continuity equation is modeled as a parabolic
equation.

Generally speaking, entropy conservation and entropy stability are used to
preserve the second law of thermodynamics in the mathematical sense, i.e., the
mathematical entropy function, $S$, decreases monotonically outside the
equilibrium. This yields entropy estimates and bounds on $S$, which can be
translated into bounds on the conservative variables, $q$, of the underlying
model Dafermos et al. (2005); Svärd (2015). In this work, nonlinear entropy
stability and the summation-by-parts (SBP) framework are used to derive
entropy stable wall boundary conditions for the Eulerian model for the viscous
and heat-conducting compressible flows proposed by Svärd (2018) and their
semidiscrete counterpart. As done in Parsani et al. (2015b) for the classical
CNS equations, a semidiscrete entropy estimate for the entire domain is
achieved when the new boundary conditions are coupled with an entropy stable
discrete interior operator. The data at the boundary are weakly imposed using
a penalty flux approach and a simultaneous-approximation-term (SAT) penalty
technique. At the semidiscrete level, the work in Parsani et al. (2015b) was
sharpened in Dalcin et al. (2019) by first constructing entropy conservative
wall boundary conditions and then adding a precise interior penalty term.
Here, we follow the semidiscrete analysis and derivation presented in Parsani
et al. (2015b); Dalcin et al. (2019) by using a collocated discontinuous
Galerkin framework based on SBP operators constructed at the
Legendre–Gauss–Lobatto (LGL) points. We verify the entropy conservation
property of the baseline boundary conditions by simulating the flow around a
rotating sphere placed in a cubic domain. In addition, to verify the accuracy
of the proposed boundary condition implementation, we present the convergence
study for a three-dimensional test case constructed using the method of
manufactured solutions.

In Svärd (2018); Dolejší and Svärd (2021), numerical simulations were
presented to highlight the difference between the classical CNS equations and
the Eulerian model. Here, we also use a set of test cases of increasing
complexity simulated using both the CNS and Eulerian models. We use the $hp$,
fully-discrete entropy stable SSDC solver described, validated, and verified
in Parsani et al. (2020). The numerical results of the two models in regions
near solid walls are compared, and differences and similarities are then
highlighted.

The paper is organized as follows. In Section 2, we present the Eulerian model
in a general form, and we show its entropy stability analysis. Then, in
Section 3, we derive the wall boundary conditions at the continuous. Later, in
Section 4, we discretize the system using SBP-SAT operators and present a
discrete entropy analysis of the Eulerian model. The latter step sets the
context for the construction of discrete entropy conservative and entropy
stable wall boundary conditions. Furthermore, in Section 5, a common SAT
procedure is presented to allow the use of a single subroutine for imposing
interface coupling and wall boundary conditions. In Section 6, we present a
set of numerical results that demonstrate the efficacy and accuracy of the new
boundary conditions. In addition, we present a few more test cases that
highlight similarities and differences between the two models slightly away
from solid walls. Finally, conclusions are drawn in Section 7.

##  2 Entropy analysis of the Eulerian model

To analyze the Eulerian model, we begin by presenting its general form and
then its entropy analysis. The latter step will set the context for deriving
solid wall boundary conditions that preserve the nonlinear entropy stability.

###  2.1 General form of the Eulerian model

Svärd Svärd (2018) arrives at the following form of the Eulerian model for
viscous and heat conducting compressible flows:

$\displaystyle\frac{\partial\rho}{\partial
t}+\frac{\partial\rho\mathcal{U}_{j}}{\partial x_{j}}$
$\displaystyle=\frac{\partial}{\partial
x_{j}}\left(\nu\frac{\partial\rho}{\partial x_{j}}\right),$ (1a)
$\displaystyle\frac{\partial\rho\mathcal{U}_{i}}{\partial
t}+\frac{\partial\rho\mathcal{U}_{i}\mathcal{U}_{j}}{\partial
x_{j}}+\frac{\partial p}{\partial x_{i}}$
$\displaystyle=\frac{\partial}{\partial
x_{j}}\left(\nu\frac{\partial\rho\mathcal{U}_{i}}{\partial x_{j}}\right),$
(1b) $\displaystyle\frac{\partial\rho\mathcal{E}}{\partial
t}+\frac{\partial\rho\mathcal{E}\mathcal{U}_{j}}{\partial
x_{j}}+\frac{\partial p\mathcal{U}_{j}}{\partial x_{j}}$
$\displaystyle=\frac{\partial}{\partial
x_{j}}\left(\nu\frac{\partial\rho\mathcal{E}}{\partial x_{j}}\right),$ (1c)

where $i,j=1:N_{\text{DIM}}$ (in MATLAB notation), and $\rho$,
$\mathcal{U}_{i}$, $p$, and $\mathcal{E}$ are the density, the velocity
component in the $x_{i}$ direction, the thermodynamic pressure and the
specific total energy, respectively.

In what follows, we assume thermodynamically perfect (or ideal) gas. Thus, the
thermodynamic pressure, $p$, is given by

$p=\rho R\mathcal{T},$ (2)

where $R$ is the gas constant. The specific total energy, $\mathcal{E}$, is
given by the following equation

$\mathcal{E}=c_{v}\mathcal{T}+\frac{1}{2}\mathcal{U}_{i}\mathcal{U}_{i},$ (3)

where $c_{v}$ is the heat capacity at constant volume and $T$ is the
temperature. The heat capacity at constant pressure, $c_{p}$, is related to
$c_{v}$ through the gas constant, i.e., $R=c_{p}-c_{v}$. Furthermore, the
generalized form of the kinematic viscosity, $\nu$, is given by

$\displaystyle\nu=\frac{\alpha\mu}{\rho(\textbf{x},t)}+\beta(\rho,\mathcal{T}),$
(4)

where $\mu$ is the dynamic viscosity, $\alpha\in[1,\frac{4}{3}]$, and $\beta$
is an additional diffusion coefficient. Finally, the speed of sound, $c$, is
defined as

$c=\sqrt{\gamma R\mathcal{T}},$ (5)

where $\gamma=\frac{c_{p}}{c_{v}}$.

###### Remark 1

Equation (4) is the most general expression of the kinematic viscosity of the
Eulerian model. However, in practice, Svärd (2018) chooses $\beta=0$ and thus,
equation (4) reduces to a scaled kinematic viscosity coefficient

$\displaystyle\nu=\frac{\alpha\mu}{\rho(\textbf{x},t)}.$ (6)

###  2.2 Entropy analysis

We cover the entropy analysis of the Eulerian model (1) by rewriting it using
the following compact form,

$\displaystyle\frac{\partial{q}}{\partial
t}+\frac{\partial{f}^{(I)}_{i}}{\partial
x_{i}}=\frac{\partial{f}^{(V)}_{i}}{\partial x_{i}},$ (7)

where ${f}^{(I)}_{i}$, and ${f}^{(V)}_{i}$ are the inviscid and viscous fluxes
in the $x_{i}$ direction, respectively.

The next theorem ensures that the Eulerian model is entropy stable.

###### Theorem 2.1

The following boundary integral

$\displaystyle\int_{\Gamma}\left({w}^{\top}\left(\nu\frac{dq}{dw}\right)\frac{\partial{w}}{\partial
x_{i}}-\mathcal{F}_{i}\right)\ d\Gamma$ (8)

bounds the time derivative of the entropy function, $S$, of the Eulerian model
(1).

###### Proof

Given an entropy pair $(\mathcal{S},\mathcal{F}_{i})$ that symmetrizes the
Eulerian model (see Svärd (2018)), we define the entropy variable
$w^{\top}=\frac{dS}{dq}$. Then, using the matrix $\frac{d{q}}{d{w}}$, we
change the variables of system (1) to $q:=q(w)$. These steps yield the
symmetric form

$\frac{d{q}}{d{w}}\frac{\partial{w}}{\partial
t}+\frac{\partial{{g}}^{(I)}_{i}}{\partial
x_{i}}=\frac{\partial{{g}}^{(V)}_{i}}{\partial x_{i}},$ (9)

where ${g}^{(I)}_{i}$, and ${g}^{(V)}_{i}$ are the symmetrized inviscid and
viscous fluxes, respectively. Next, we multiply the symmetric system (9) from
the left by the entropy variables. Thus, using the chain rule, the term
associated with the time derivative reads

${w}^{\top}\frac{d{q}}{d{w}}\frac{\partial{w}}{\partial
t}=\frac{d\mathcal{S}}{d{q}}\frac{d{q}}{d{w}}\frac{\partial{w}}{\partial
t}=\frac{\partial\mathcal{S}}{\partial t}.$ (10)

Further, as shown by Tadmor (2003), the contribution of the symmetrized
inviscid flux term yields

${w}^{\top}\frac{\partial{{g}}^{(I)}_{i}}{\partial
x_{i}}=\frac{\partial\mathcal{F}_{i}}{\partial x_{i}},$ (11)

i.e., the divergence of the entropy flux. Lastly, the viscous term
contribution can be manipulated to obtain the following expression:

$w^{\top}\frac{\partial{{g}}^{(V)}_{i}}{\partial
x_{i}}=w^{\top}\frac{\partial}{\partial
x_{i}}\left(\left(\nu\frac{dq}{dw}\right)\frac{\partial{w}}{\partial
x_{i}}\right).$ (12)

Therefore, using Equations (10), (11) and (12), the scalar partial
differential equation for the entropy function, $\mathcal{S}$, integrated over
a generic domain $\Omega$ with boundary $\Gamma$ reads

$\int_{\Omega}\frac{\partial\mathcal{S}}{\partial
t}+\frac{\partial\mathcal{F}_{i}}{\partial x_{i}}\
dx=\int_{\Omega}w^{\top}\frac{\partial}{\partial
x_{i}}\left(\left(\nu\frac{dq}{dw}\right)\frac{\partial{w}}{\partial
x_{i}}\right)\ dx.$ (13)

Now, using the fundamental theorem of calculus on
$\frac{\partial\mathcal{F}_{i}}{\partial x_{i}}$, and the integration-by-parts
rule on $w^{\top}\frac{\partial{{g}}^{(V)}_{i}}{\partial x_{i}}$ gives

$\int_{\Omega}\frac{\partial\mathcal{S}}{\partial
t}=\int_{\Gamma}\left(w^{\top}\left(\nu\frac{dq}{dw}\right)\frac{\partial{w}}{\partial
x_{i}}-\mathcal{F}_{i}\right)\ d\Gamma-DT,$ (14)

where

$DT=\int_{\Omega}\frac{\partial w^{\top}}{\partial
x_{i}}\left(\nu\frac{dq}{dw}\right)\frac{\partial{w}}{\partial x_{i}}dx.$ (15)

Because the kinematic viscosity, $\nu$, is positive and the matrix
$\frac{dq}{dw}$ is symmetric positive-definite, the product
$\left(\nu\frac{dq}{dw}\right)$ is also a positive-definite matrix. Thus, $DT$
is positive and appropriately bounded (See Svärd (2015)). We obtain the
following expression

$\int_{\Omega}\frac{\partial\mathcal{S}}{\partial t}\
dx\leq\int_{\Gamma}\left(w^{\top}\left(\nu\frac{dq}{dw}\right)\frac{\partial{w}}{\partial
x_{i}}-\mathcal{F}_{i}\right)\ d\Gamma$ (16)

which completes the proof.

In the next section, we derive solid wall boundary conditions such that the
boundary integral on the left-hand side of (16) is bounded by data.

##  3 Entropy stable wall boundary condition

In this section, we derive the boundary conditions for preserving the entropy
stability of the system (1) for a solid wall. The analysis presented next
closely follows the works of Parsani et al. (2015b); Dalcin et al. (2019).

The entropy-entropy flux pair $(\mathcal{S},\mathcal{F}_{i})=(-\rho
s,-\rho\mathcal{U}_{i}s)$, used in Dalcin et al. (2019); Parsani et al.
(2015b), is the only pair that admits a diffusive entropy flux for the CNS
model Svärd (2018). Therefore, it is the only suitable pair for a comparison
between the CNS and Eulerian systems of equations. This pair is also used in
this work. Thus, the entropy variables are computed as Fisher (2012);
Carpenter et al. (2014)

$w=\left(\frac{h}{\mathcal{T}}-s-\frac{\mathcal{U}_{i}\mathcal{U}_{i}}{2\mathcal{T}},-\frac{\mathcal{U}_{1}}{\mathcal{T}},-\frac{\mathcal{U}_{2}}{\mathcal{T}},-\frac{\mathcal{U}_{3}}{\mathcal{T}},-\frac{1}{\mathcal{T}}\right)^{\top},$
(17)

where $h=c_{pr}\mathcal{T}$ is the enthalpy.

###### Remark 2

The entropy-entropy flux pair for the Eulerian model is not necessarily
restricted to $(\mathcal{S},\mathcal{F}_{i})=(-\rho s,-\rho\mathcal{U}_{i}s)$.
Svärd (2018) mentions that the Eulerian model admits a diffusive entropy flux
for all Harten’s generalized entropies.

To simplify, we rewrite the Eulerian form (7) as

$\frac{\partial{q}}{\partial t}+\frac{\partial{f}^{(I)}_{i}}{\partial
x_{i}}=\frac{\partial}{\partial x_{i}}(\nu\frac{\partial{q}}{\partial
x_{i}}).$ (18)

From the proof of Theorem 2.1, we arrive at

$\frac{\partial\mathcal{S}}{\partial
t}+\frac{\partial\mathcal{F}_{i}}{\partial
x_{i}}={w}^{\top}\frac{\partial}{\partial
x_{i}}\left(\left(\nu\frac{dq}{dw}\right)\frac{\partial{w}}{\partial
x_{i}}\right).$ (19)

Integrating the previous expression over a cubical domain,
$\Omega=(0,1)\times(0,1)\times(0,1)$, and, without loss of generality,
considering a solid wall located at $x_{i}=0$, we get the following
contributions Parsani et al. (2015b)

$\int_{\Omega}\frac{\partial\mathcal{S}}{\partial t}\
dx=-\int_{x_{i}=0}\left[{w}^{\top}\left(\nu\frac{dq}{dw}\right)\frac{\partial{w}}{\partial
x_{i}}-\mathcal{F}_{i}\right]dx_{j}dx_{k}-DT,$ (20)

where $j,k\neq i$. From Theorem 2.1, we get the bound

$\int_{\Omega}\frac{\partial\mathcal{S}}{\partial t}\
dx\leq-\int_{x_{i}=0}\left[{w}^{\top}\left(\nu\frac{dq}{dw}\right)\frac{\partial{w}}{\partial
x_{i}}-\mathcal{F}_{i}\right]dx_{j}dx_{k}.$ (21)

To bound the time derivative of the entropy, the right-hand-side of (21)
requires boundary data. For a solid viscous wall, assuming linear analysis,
five independent boundary conditions are required Svärd (2018). Four comes
from the linear analysis of the CNS model Berg and Nordström (2011); Svärd and
Nordström (2008), in addition to a boundary condition that rises from density
diffusion Svärd (2018). Four boundary conditions are the three no-slip
boundary conditions, as the CNS case, and one on the density flux presented
below. The fifth condition is the gradient of the temperature normal to the
wall (Neumann boundary condition; e.g., the adiabatic wall), or the
temperature of the wall (the Dirichlet or isothermal wall boundary condition),
or a mixture of Dirichlet and Neumann conditions (the Robin boundary
condition). The last condition is the same as for the CNS model Berg and
Nordström (2011); Svärd and Nordström (2008); Svärd (2018). These five
boundary conditions lead to linear stability of both the CNS and Eulerian
model Berg and Nordström (2011); Svärd and Nordström (2008); Svärd (2018). In
the remainder of this section, we will show the type and the form of the wall
boundary conditions that have to be imposed to bound the estimate (21) and,
hence, to attain entropy stability.

###  3.1 Inviscid contribution

The inviscid contribution to the time rate of change of the entropy function,

$\int_{x_{i}=0}\mathcal{F}_{i}\ dx_{j}dx_{k}.$ (22)

appearing on the RHS of (21), is treated as in (Dalcin et al., 2019, Theorem
2.2). However, for clarity of presentation and completeness, we report it
here.

###### Theorem 3.1

The no-slip boundary conditions $\mathcal{U}_{i}=0$ and
$\mathcal{U}_{j}=\mathcal{U}_{j}^{\text{wall}}$, where $j\neq i$, bound the
inviscid contribution of the time derivative of the entropy function, i.e.

$\mathcal{F}_{i}=0,$ (23)

for an inviscid solid wall with an outfacing normal vector pointing in the
$x_{i}$ direction.

###### Proof

The proof of this theorem can be found in (Dalcin et al., 2019, Theorem 2.2).

###  3.2 Viscous contribution

In this subsection, we derive the viscous boundary conditions for a no-slip
wall associated to the first term of the RHS of (21), i.e.,

$-\int_{x_{i}=0}{w}^{\top}\left(\nu\frac{dq}{dw}\right)\frac{\partial{w}}{\partial
x_{i}}\ dx_{j}dx_{k}.$ (24)

The main result is provided in the following theorem.

###### Theorem 3.2

The boundary conditions

$\displaystyle
g(t)=\mu\frac{1}{\mathcal{T}}\frac{\partial\mathcal{T}}{\partial
x_{i}},\qquad\textrm{and}\qquad\frac{1}{\rho}\frac{\partial\rho}{\partial
x_{i}}=0,$ (25)

bound the viscous contribution to the time derivative of the entropy function
for a no-slip solid wall with an outfacing normal vector pointing in the
$x_{i}$ direction.

###### Proof

The explicit evaluation of the viscous contribution gives

$-{w}^{\top}\left(\nu\frac{dq}{dw}\right)\frac{\partial{w}}{\partial
x_{i}}=\mu\frac{1}{\mathcal{T}}\frac{\partial\mathcal{T}}{\partial
x_{i}}-\mu(s+(\gamma-1))\frac{1}{\rho}\frac{\partial\rho}{\partial x_{i}}$
(26)

that contributes positively to the time derivative of the entropy function.
Then, setting $\frac{1}{\rho}\frac{\partial\rho}{\partial x_{i}}=0$ yields the
boundary condition

$-{w}^{\top}\left(\nu\frac{dq}{dw}\right)\frac{\partial{w}}{\partial
x_{i}}=\mu\frac{1}{\mathcal{T}}\frac{\partial T}{\partial x_{i}}=g(t)$ (27)

which completes the proof.

For the CNS model, the viscous contribution to the time derivative of the
entropy function can be bounded by only setting
$\kappa\frac{1}{\mathcal{T}}\frac{\partial T}{\partial x_{i}}$ Parsani et al.
(2015b). However, for the Eulerian model, we require the impostion of two
boundary conditions, i.e.
$g(t)=\mu\frac{1}{\mathcal{T}}\frac{\partial\mathcal{T}}{\partial x_{i}}$ and
$\frac{1}{\rho}\frac{\partial\rho}{\partial x_{i}}=0$.

In practice with the above boundary conditions, the time derivative of the
entropy function satisfies the following relation

$\int_{\Omega}\frac{\partial\mathcal{S}}{\partial t}\
dx\leq\mu\frac{1}{\mathcal{T}}\frac{\partial\mathcal{T}}{\partial
x_{i}}=g(t)=\text{DATA}.$ (28)

Bound (28) can be translated into a bound for the conserved quantities
Dafermos et al. (2005); Svärd (2015), and hence, for the primitive variables.

##  4 Semidiscrete entropy stable framework

Herein, using summation-by-parts (SBP) operators Svärd and Nordström (2014)
and the simultaneous-approximation-technique (SAT) Carpenter et al. (1994);
Nordström and Carpenter (1999), we provide an entropy stable framework of any
order for the semidiscretization of the Eulerian model (1) using unstructured
grids.

###  4.1 SBP operators

The one-dimensional SBP operator for the first derivative in the direction
$x_{i}$ is defined as the following.

###### Definition 1

Summation–by–parts (SBP) operator for the first derivative: A matrix operator
with constant coefficients, $\mathcal{D}\in\mathbb{R}^{N\times N}$, is a
linear SBP operator of degree $p$ approximating the derivative
$\frac{\partial}{\partial x_{i}}$ on the domain
$x_{i}\in\left[a,b\right]$ with nodal distribution $x$ having $N$ nodes, if

  1. 1.

$\mathcal{D}x_{i}^{j}=jx_{i}^{j-1}$, $j=0,1,\cdots,p;$

  2. 2.

$\mathcal{D}=\mathcal{P}^{-1}\mathcal{Q}$, where the norm matrix $\mathcal{P}$
is symmetric positive-definite;

  3. 3.

$\mathcal{Q}+\mathcal{Q}^{\top}=\mathcal{B}$, where
$\mathcal{B}=\mathrm{diag}\left[-1,0,\cdots,0,1\right]$.

In other words, an SBP operator of degree $p$ is one that exactly
differentiates monomials up to degree $p$ ($p=N-1$).

In this work, a collocated discontinuous Galerkin approach is used.
Specifically, diagonal norm SBP operators are constructed on the LGL nodes.
The one-dimensional SBP operators used in this work are explicitly constructed
in Carpenter et al. (2015). Their extension to two- and three-dimensions is
achieved using tensor product operations Carpenter et al. (2015); Parsani et
al. (2015b):

$\begin{gathered}\mathcal{D}_{x_{1}}=\mathcal{D}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{5},\quad\cdots\quad\mathcal{D}_{x_{3}}=\mathcal{I}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{D}_{N}\otimes\mathcal{I}_{5},\\\
\mathcal{Q}_{x_{1}}=\mathcal{Q}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{5},\quad\cdots\quad\mathcal{Q}_{x_{3}}=\mathcal{I}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{Q}_{N}\otimes\mathcal{I}_{5},\\\
\mathcal{B}_{x_{1}}=\mathcal{B}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{5},\quad\cdots\quad\mathcal{B}_{x_{3}}=\mathcal{I}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{B}_{N}\otimes\mathcal{I}_{5},\\\
\Delta_{x_{1}}=\Delta_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{5},\quad\cdots\quad\Delta_{x_{3}}=\mathcal{I}_{N}\otimes\mathcal{I}_{N}\otimes\Delta_{N}\otimes\mathcal{I}_{5},\\\
\mathcal{P}_{x_{1}}=\mathcal{P}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{5},\quad\cdots\quad\mathcal{P}_{x_{3}}=\mathcal{I}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{P}_{N}\otimes\mathcal{I}_{5},\\\
\mathcal{P}_{x_{1},x_{2}}=\mathcal{P}_{N}\otimes\mathcal{P}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{,}5\quad\cdots\quad\mathcal{P}_{x_{2},x_{3}}=\mathcal{I}_{N}\otimes\mathcal{P}_{N}\otimes\mathcal{P}_{N}\otimes\mathcal{I}_{5},\\\
\mathcal{P}_{x_{1},x_{2},x_{3}}=\mathcal{P}_{N}\otimes\mathcal{P}_{N}\otimes\mathcal{P}_{N}\otimes\mathcal{I}_{5},\end{gathered}$
(29)

where $\mathcal{D}_{N}$, $\mathcal{Q}_{N}$, $\mathcal{B}_{N}$, $\Delta_{N}$
and $\mathcal{P}_{N}$ are the one-dimensional SBP operators, and
$\mathcal{I}_{N}$ is the identity operator.111In this work, we use $N=5$
because in three dimensions the number of partial differential equations is
five for both the CNS and Eulerian model. In this context, we chose diagonal
$\mathcal{P}_{N}$, and by the definition of SBP operators
$\mathcal{D}_{N}=\mathcal{P}_{N}^{-1}\mathcal{Q}_{N}$ and
$\mathcal{B}_{N}=\mathcal{Q}_{N}^{\top}+\mathcal{Q}_{N}$. The matrices
$\mathcal{B}_{(\cdot)}$ pick off the interface terms in the respective
directions. For the spectral element discretization considered in this paper,
the $\mathcal{B}_{(\cdot)}$ matrices take on a particularly simple form; as an
example, consider $\mathcal{B}_{x_{1}}$, which is given as

$\begin{gathered}\mathcal{B}_{x_{1}}=\mathcal{B}_{x_{1}}^{+}-\mathcal{B}_{x_{1}}^{-},\\\
\mathcal{B}_{x_{1}}^{-}=\operatorname{diag}\left(1,0,\dots,0\right)\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{5},\\\
\mathcal{B}_{x_{1}}^{+}=\operatorname{diag}\left(0,\dots,0,1\right)\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{N}\otimes\mathcal{I}_{5}.\end{gathered}$

For a high-order accurate scheme on a tensor product cell, they pick off the
values of whatever vector they act on (typically the solution or the flux) at
the nodes of the two opposite faces multiplied by the orthogonal component of
the unit normal.

When applying any of these operators to the scalar entropy equation in space,
a hat will be used to differentiate the scalar operator from the full vector
operator, e.g.

$\widehat{\mathcal{P}}=\left(\mathcal{P}_{N}\otimes\mathcal{P}_{N}\otimes\mathcal{P}_{N}\right).$

We finally note that in the present work, the quadrature nodes and solution
nodes are collocated.

###  4.2 Semidiscretization of the Eulerian model

Using the operators shown in (29), we can write the semidiscretization of (1)
as

$\frac{\partial\mathbf{q}_{i}}{\partial
t}=\mathcal{D}_{x_{i}}\left(\mathbf{f}^{(V)}_{i}-\mathbf{f}^{(I)}_{i}\right)+\mathcal{P}^{-1}_{x_{i}}\left(\mathbf{g}^{(b)}_{x_{i}}+\mathbf{g}^{(In)}_{x_{i}}\right),$
(30)

where vectors $\mathbf{g}^{(b)}_{x_{i}}$ enforce the boundary conditions,
while $\mathbf{g}^{(In)}_{x_{i}}$ patches interfaces together using a SAT
approach Parsani et al. (2015b). The bolded letters represent quantities,
$\mathbf{q}$, and functions, $\mathbf{f}$ and $\mathbf{g}$ for all nodes in an
element.

Following Fisher (2012); Carpenter et al. (2014), we use the telescoping
property of an SBP operator,

$\mathcal{D}_{x_{i}}\mathbf{f}^{(I)}_{i}=\mathcal{P}^{-1}_{x_{i}}\mathcal{Q}_{x_{i}}\mathbf{f}^{(I)}_{i}=\mathcal{P}^{-1}_{x_{i}}\Delta_{x_{i}}\mathcal{I}_{StoF}\mathbf{f}^{(I)}_{i}=\mathcal{P}^{-1}_{x_{i}}\Delta_{x_{i}}\bar{\mathbf{f}}^{(I)}_{i},$
(31)

to re-write the semidiscrete counterpart to the equation (30) as

$\frac{\partial\mathbf{q}}{\partial
t}=\left(\mathcal{D}_{x_{i}}\mathbf{f}^{(V)}_{i}-\mathcal{P}^{-1}_{x_{i}}\Delta_{x_{i}}\bar{\mathbf{f}}^{(I)}_{i}\right)+\mathcal{P}^{-1}_{x_{i}}\left(\mathbf{g}^{(b)}_{x_{i}}+\mathbf{g}^{(In)}_{x_{i}}\right).$
(32)

The vector $\bar{\mathbf{f}}^{(I)}_{i}$ is defined as
$\bar{\mathbf{f}}^{(I)}_{i}=\mathcal{I}_{StoF}\mathbf{f}^{(I)}_{i}$ and the
one-dimensional telescoping operator, $\Delta_{N}$, is given by

$\Delta_{N}=\left(\begin{array}[]{cccccc}-1&1&0&0&0&0\\\ 0&-1&1&0&0&0\\\
0&0&\ddots&\ddots&0&0\\\ 0&0&0&-1&1&0\\\ 0&0&0&0&-1&1\end{array}\right).$

The operator $\mathcal{I}_{StoF}$ interpolates the value at the solution nodes
to the interfaces between nodes as shown in Figure 1.

Additionally, we define $\left[\frac{dq}{dw}\right]$ as a block diagonal
matrix applied to all LGL points in an element. The kinematic viscosity $\nu$
is a scalar and we can bring it inside the matrix
$\left[\frac{dq}{dw}\right]$, i.e., $\left[\nu\frac{dq}{dw}\right]$. The block
diagonal matrix $\left[\nu\frac{dq}{dw}\right]$ commutes with the operator
$\mathcal{D}_{x_{i}}$. Thus, we can re-write the viscous flux as

$\mathbf{f}^{(V)}_{i}=\nu\mathcal{D}_{x_{i}}\mathbf{q}_{i}=\nu\mathcal{D}_{x_{i}}\left[\frac{dq}{dw}\right]\mathbf{w}_{i}=\left[\nu\frac{dq}{dw}\right]\mathcal{D}_{x_{i}}\mathbf{w}_{i}=\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta},$
(33)

where the symbol $\mathbf{\Theta}_{x_{i}}$ is the gradient of the entropy
variables in the $x_{i}$ direction.

Figure 1: Onedimensional discretization using LGL points of order $p=4$.
$\cdot$ and $\times$ denote solution and flux points, respectively. Reprinted
from Carpenter et al. (2014) with permission.

Using the procedure based on local discontinuous Galerkin (LDG) and interior
penalty approach (IP) described in Parsani et al. (2015b); Dalcin et al.
(2019), the semidiscretization (32) can be recast as

$\displaystyle\frac{\partial\mathbf{q}_{i}}{\partial
t}+\mathcal{P}^{-1}_{x_{i}}\Delta_{x_{i}}\bar{\mathbf{f}}^{(I)}_{i,sc}-\mathcal{D}_{x_{i}}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{i}}=\mathcal{P}^{-1}_{x_{i}}\left(\mathbf{g}^{(b),q}_{x_{i}}+\mathbf{g}^{(In),q}_{x_{i}}\right),$
(34a)
$\displaystyle\mathbf{\Theta}_{x_{i}}-\mathcal{D}_{x_{i}}\mathbf{w}=\mathcal{P}^{-1}_{x_{i}}\left(\mathbf{g}^{(b),\Theta}_{x_{i}}+\mathbf{g}^{(In),\Theta}_{x_{i}}\right).$
(34b)

The terms $\mathbf{g}^{(b),q}_{x_{i}}$, $\mathbf{g}^{(In),q}_{x_{i}}$, and
$\mathbf{g}^{(b),\Theta}_{x_{i}}$, $\mathbf{g}^{(In),\Theta}$ are the SAT
penalty boundary ($b$) and interface ($In$) terms on the conservative
variables, $q$, and the gradient of the entropy variables, $\Theta$,
respectively. The contributions of the interface penalty terms are non-zero
only in the normal direction to the interface.

###### Remark 3

To build a high-order entropy conservative (or, equivalently, entropy stable)
semi discretization, the linear interpolation operator $\mathcal{I}_{StoF}$ is
replaced with a nonlinear interpolation operator Fisher (2012); Carpenter et
al. (2014). Thus, the operator
$\mathcal{P}^{-1}_{x_{i}}\Delta_{x_{i}}\bar{\mathbf{f}}^{(I)}_{i,sc}$ (or,
$\mathcal{P}^{-1}_{x_{i}}\Delta_{x_{i}}\bar{\mathbf{f}}^{(I)}_{i,ss}$ for the
entropy stable case) is a nonlinear operator that is a discrete counterpart to
the term that appears in Theorem 3.1. Thus, we arrive at the following
relation

$\mathbf{1}^{\top}\mathcal{P}\mathcal{P}^{-1}_{x_{i}}\Delta_{x_{i}}\bar{\mathbf{f}}^{(I)}_{i,sc}=\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\widehat{\mathcal{B}}_{x_{i}}\mathbf{F}_{x_{i}}\simeq\int_{x_{i}=1}\mathcal{F}_{i}\
dx_{j}dx_{k}-\int_{x_{i}=0}\mathcal{F}_{i}\ dx_{j}dx_{k}.$ (35)

The flux vector $\bar{\mathbf{f}}^{(I)}_{i,sc}$ is constructed using an
entropy conservative two-point flux Tadmor (2003), and therefore, it satisfies
the relation (11). In the case of an entropy stable flux,
$\bar{\mathbf{f}}^{(I)}_{i,ss}$, the resulting term is entropy dissipative
Carpenter et al. (2014).

###  4.3 Time derivative of entropy function

Following the entropy analysis detailed in Carpenter et al. (2014); Parsani et
al. (2015b, a), the semidiscretization (34) yields the following expression
for the time derivative of the entropy function, $\mathcal{S}$:

$\displaystyle\frac{d}{dt}\mathbf{1}^{\top}\widehat{\mathcal{P}}\mathbfcal{S}+\mathbf{DT}=\mathbf{\Xi}$
(36)

where

$\mathbf{DT}=\left\|\left[\nu\frac{dq}{dw}\right]^{\frac{1}{2}}\mathbf{\Theta}_{x_{i}}\right\|^{2}_{\mathcal{P}},$
(37)

and

$\displaystyle\begin{split}\mathbf{\Xi}&=-\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\widehat{\mathcal{B}}_{x_{i}}\mathbf{F}_{x_{i}}\\\
&+\mathbf{w}^{\top}\mathcal{P}_{x_{j},x_{k}}\mathcal{B}_{x_{i}}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{i}}\\\
&+\mathbf{w}^{\top}\mathcal{P}_{x_{j},x_{k}}\left(\mathbf{g}^{(b),q}_{x_{i}}+\mathbf{g}^{(In),q}_{x_{i}}\right)\\\
&+\left(\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{i}}\right)^{\top}\mathcal{P}_{x_{j},x_{k}}\left(\mathbf{g}^{(b),\Theta}_{x_{i}}+\mathbf{g}^{(In),\Theta}_{x_{i}}\right).\end{split}$
(38)

The full derivation of (36) is detailed in Appendix A.

Now, consider a cubic element of length 1 with a solid wall with a normal
vector in the $x_{i}$ direction, we have that the interface penalty terms are
zero, i.e., $\mathbf{g}^{(In),q}_{x_{i}}=\mathbf{g}^{(In),\Theta}_{x_{i}}=0$.
Therefore, equation (38) reduces to

$\displaystyle\begin{split}\mathbf{\Xi}&=-\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\widehat{\mathcal{B}}_{x_{i}}\mathbf{F}_{x_{i}}\\\
&+\mathbf{w}^{\top}\mathcal{P}_{x_{j},x_{k}}\mathcal{B}_{x_{i}}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{i}}\\\
&+\mathbf{w}^{\top}\mathcal{P}_{x_{j},x_{k}}\mathbf{g}^{(b),q}_{x_{i}}\\\
&+\left(\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{i}}\right)^{\top}\mathcal{P}_{x_{j},x_{k}}\mathbf{g}^{(b),\Theta}_{x_{i}}.\end{split}$
(39)

###  4.4 Entropy stable wall boundary conditions for the semidiscrete system

The boundary condition penalty term with respect to the conservative variables
is split into three design-order terms plus one source boundary term:

$\mathbf{g}^{(b),q}_{x_{i}}=\mathbf{g}^{(b,I),q}_{x_{i}}+\mathbf{g}^{(b,V),q}_{x_{i}}+\mathbfcal{M}^{(b,V)}+\mathbfcal{L}^{(b,V)}.$
(40)

The first component of each term is computed from the numerical solution, and
the second component is constructed from a combination of the numerical
solution and five independent boundary data, as done in Parsani et al.
(2015b).

The first term enforces Euler no-penetration wall

$\mathbf{g}^{(b,I),q}_{x_{i}}=-\frac{1}{2}\mathcal{B}^{-}_{x_{i}}\left(\mathbf{f}^{(I)}_{i}-\mathbf{f}^{(b,I)}_{i,sc}(v^{(E)},v^{(b,I)})\right),$
(41)

where $\mathbf{f}^{(b,I)}_{i,sc}$ is the entropy conservative flux and $v$
represent the primitive variables. The manufactured inviscid states are
defined as in Parsani et al. (2015b) for $i=1$, without loss of generality,

$\displaystyle
v^{(E)}=\operatorname{diag}([1,-1,1,1,1]),\quad\textrm{and}\quad$
$\displaystyle
v^{(b,I)}=\left(\rho,-\mathcal{U}_{1},\mathcal{U}_{2},\mathcal{U}_{3},\mathcal{T}\right)^{\top}.$
(42)

The viscous boundary term is given as

$\mathbf{g}^{(b,V),q}_{x_{i}}=\frac{1}{2}\mathcal{B}^{-}_{x_{i}}\left(\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{i}}-\mathbf{f}^{(b,V)}_{i}\right),$
(43)

where $\mathbf{f}^{(b,V)}_{i}$ will be defined later in this section. Together
with (43), the following term

$\mathbf{g}^{(b,V),\Theta}_{x_{i}}=\frac{1}{2}\mathcal{B}^{-}_{x_{i}}\left(\mathbf{w}-\mathbf{w}^{(b,V)}\right),$
(44)

enforce the no-slip boundary condition weakly. The analysis of the previous
viscous terms match the analysis described by Parsani et al. (2015b); Dalcin
et al. (2019) since the entropy variables are the same, and the matrix
$\left[\nu\frac{dq}{dw}\right]$ is symmetric positive-definite.

The remaining terms are defined in point-wise notation, and thus the bold
notation is dropped. The analysis of their contributions is summarized next.

  * •

The manufactured viscous flux is defined as

$f^{(b,V)}_{i}=\left(\nu\frac{dq}{dw}\left(v^{(b,V)}\right)\right)\widetilde{\Theta}_{x_{i}},$
(45)

where the manufactured viscous boundary primitive variables, $v^{(b,V)}$, is
defined, as in Dalcin et al. (2019), by

$v^{(b,V)}=\left(\rho,-\mathcal{U}_{1}+2\mathcal{U}_{1}^{\text{wall}},-\mathcal{U}_{2}+2\mathcal{U}_{2}^{\text{wall}},-\mathcal{U}_{3}+2\mathcal{U}_{3}^{\text{wall}},\mathcal{T}\right)^{\top}.$
(46)

The change of variables matrix, $\frac{dq}{dw}(v)$, as a function of the
primitive variables, $v$, is defined in Appendix B. The term
$\widetilde{\Theta}_{x_{i}}$ is the manufactured gradient of the entropy
variables. The construction of $\widetilde{\Theta}_{x_{i}}$ is summarized in
Appendix C.

  * •

The dissipative IP term is an averaged state of the viscous flux evaluated at
the primitive state, $v$, and the manufactured primitive state, $v^{(b,V)}$,

$\mathcal{M}^{(b,V)}=-\beta\frac{f^{(V)}(v)+f^{(V)}(v^{(b,V)})}{2},$ (47)

where $\beta$ is a positive constant that is scaled with the inverse of the
element length in the normal direction, controlling the strength of the
penalty term Dalcin et al. (2019).

  * •

The source term

$\mathcal{L}^{(b,V)}=-\left(0,0,0,0,1\right)^{\top}g(t),$ (48)

imposes the heat flux boundary condition, as done in Parsani et al. (2015b).

We summarize the results for the RHS of (36) in the following three theorems.
The first theorem is a result from Parsani et al. (2015b) and is reprinted
here for completeness

###### Theorem 4.1

The penalty inviscid flux contribution, $\mathbf{g}^{(b,I),q}$, is entropy
conservative if the vector $v^{(b,I)}$ is defined as in (42).

###### Proof

The proof of this theorem can be found in (Parsani et al., 2015b, Theorem
5.1).

###### Theorem 4.2

The penalty terms for the viscous flux in the conserved varaibles,
$\mathbf{g}^{(b,V),q}$, and the gradient of entropy variables,
$\mathbf{g}^{(b,V),\Theta}$, are

  * •

entropy conservative if the wall is adiabatic, i.e. $g(t)=0$,

  * •

entropy stable in the presence of a heatflux, $g(t)\neq 0$, where $g(t)$ is a
given $L^{2}$ function.

###### Proof

Substituting the expressions for $\mathbf{g}^{(b,V),q}$, and,
$\mathbf{g}^{(b,V),\Theta}$ into (39), yields

$\frac{d}{dt}\mathbf{1}^{\top}\widehat{\mathcal{P}}\mathbfcal{S}+\mathbf{DT}=\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}g(t).$
(49)

The function $g(t)$ for an adiabatic wall is $g(t)=0$. Thus the RHS is zero,
and the scheme is entropy conservative. Otherwise, the term $g(t)$ is bounded
by data. Hence, the RHS is bounded, completing the proof.

The full details on the computation of (49) are reported in Appendix C.

###### Theorem 4.3

The IP term, $\mathbfcal{M}^{(b,V)}$ , added to (40) is entropy dissipative.

###### Proof

By expanding the penalty terms with respect to the conservative variables,
$g^{(b,V),q}$, and focusing only on the dissipative IP term,
$\mathcal{M}^{(b,V)}$, in (36), we arrive at

$w^{\top}\mathcal{M}^{(b,V)}=-\frac{2\beta\alpha\mu}{R\mathcal{T}^{2}}\|\Delta\mathcal{U}\|^{2}\left(\|\Delta\mathcal{U}\|^{2}+R\mathcal{T}\right).$
(50)

Thus, $\mathcal{M}^{(b,V)}$ is entropy dissipative because all the parameters
and variables appearing on the RHS of (50) are positive. This completes the
proof.

The full details on the compution of (50) are reported in Appendix D.

##  5 A common SAT procedure for the imposition of wall boundary conditions
and interior interface coupling

The proposed approach for imposing the solid wall boundary conditions allow
for a SAT implementation that is identical to the interface treatment shown in
Parsani et al. (2015a). We can use a single subroutine with different inputs
corresponding to the imposition of the interior interface couplings, or the
adiabatic solid wall, or the wall with a prescribed heat entropy flow. In
fact, the interior interface coupling can be written as (see equations
(16a-16d) in Parsani et al. (2015a))

$\displaystyle\frac{\partial\mathbf{q}_{l}}{\partial
t}+\left(\mathcal{P}^{-1}_{x_{i,l}}\Delta_{x_{i,l}}\bar{\mathbf{f}}^{(I)}_{i,l,sc}-\mathcal{D}_{x_{i,l}}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{i,l}}\right)=\mathcal{P}^{-1}_{x_{i,l}}\mathbf{g}_{x_{i,l}}^{(In),q},$
(51a)
$\displaystyle\mathbf{\Theta}_{x_{i,l}}-\mathcal{D}_{x_{i,l}}\mathbf{w}=\mathcal{P}^{-1}_{x_{i,l}}\mathbf{{g}}_{x_{i,l}}^{(In),\Theta},$
(51b) $\displaystyle\frac{\partial\mathbf{q}_{r}}{\partial
t}+\left(\mathcal{P}^{-1}_{x_{i,r}}\Delta_{x_{i,r}}\bar{\mathbf{f}}^{(I)}_{i,r,sc}-\mathcal{D}_{x_{i,r}}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{i,r}}\right)=\mathcal{P}^{-1}_{x_{i,r}}\mathbf{{g}}_{x_{i,r}}^{(In),q},$
(51c)
$\displaystyle\mathbf{\Theta}_{x_{i,r}}-\mathcal{D}_{x_{i,r}}\mathbf{w}=\mathcal{P}^{-1}_{x_{i,r}}\mathbf{g}_{x_{i},r}^{(In),\Theta}.$
(51d)

The above equations have exactly the same structure as the LDG-IP approach
used for the imposition of the solid wall boundary conditions except for the
boundary penalty interface terms, $\mathbf{{g}}_{x_{i,r}}^{(b),\cdot}$ in
equation (34), which are replaced by the interior penalty interface coupling
terms, $\mathbf{{g}}_{x_{i,r}}^{(In),\cdot}$ in equations (51).

##  6 Numerical results

In this section, we numerically investigate the proposed entropy stable wall
boundary conditions. The numerical experiments reported in this paper are
performed with the entropy stable collocated Discontinuous Galerkin algorithm
and relaxation Runge–Kutta schemes implemented in the $hp$-adaptive,
unstructured, curvilinear grid framework SSDC Parsani et al. (2020). SSDC is
developed in the AANSLab, which is part of the Extreme Computing Research
Center at King Abdullah University of Science and Technology (KAUST). The core
entropy stable adaptive algorithms of SSDC are built on top of the Portable
and Extensible Toolkit for Scientific computing (PETSc) Balay et al. (2021),
its mesh topology abstraction (DMPLEX) Knepley and Karpeev (2009), and its
scalable ODE/DAE solver library Abhyankar et al. (2018). The SSDC framework
uses a non-dimensional formalism; thus, all quantities are scaled to units.
For the Eulerian model, the kinematic viscosity is scaled using $\alpha=1$. We
use the two-point entropy conservative flux presented by Chandrashekar (2013).
Furthermore, all the simulations have been performed in double (machine)
precision. For all the cases, we use the Runge–Kutta scheme of Bogacki and
Shampine (1989) with adaptive time-step and both relative and absolute
tolerances set to $10^{-}8$. The meshes are generated using Gmsh Geuzaine and
Remacle (2009), and Pointwise V18.3 released in September, 2019\. The SSDC
Eulerian and SSDC CNS data computed for this section is available in
http://doi.org/10.5281/zenodo.5041436.

###  6.1 Convergence study

Figure 2: Axial velocity, $\mathcal{U}_{1}$, in a pipe with annular section; Eulerian model with $\text{Re}=1$, $\text{Ma}=1e-3$. Grid | $L^{1}$ | Rate | $L^{2}$ | Rate | $L^{\infty}$ | Rate
---|---|---|---|---|---|---
$4$ | $8.36276e-04$ | - | $1.77874e-03$ | - | $5.70764e-03$ | -
$8$ | $1.25039e-04$ | -2.742 | $3.44481e-04$ | -2.368 | $1.49599e-03$ | -1.932
$16$ | $1.99325e-05$ | -2.649 | $5.46416e-05$ | -2.656 | $3.12259e-04$ | -2.260
$32$ | $2.64237e-06$ | -2.915 | $7.50755e-06$ | -2.864 | $5.07517e-05$ | -2.621
$64$ | $3.31631e-07$ | -2.994 | $9.65089e-07$ | -2.960 | $7.22735e-06$ | -2.812
$128$ | $4.12357e-08$ | -3.008 | $1.21565e-07$ | -2.989 | $9.60653e-07$ | -2.911
Table 1: Convergence study for the flow in a pipe with annular cross-section at $\text{Re}=1$ and $\text{Ma}=1e-3$. axial velocity error; solution polynomial degree: $p=2$. Grid | $L^{1}$ | Rate | $L^{2}$ | Rate | $L^{\infty}$ | Rate
---|---|---|---|---|---|---
$4$ | $2.00203e-04$ | - | $3.85016e-04$ | - | $1.54031e-03$ | -
$8$ | $2.71540e-05$ | -2.882 | $5.67563e-05$ | -2.762 | $3.16887e-04$ | -2.281
$16$ | $2.15610e-06$ | -3.655 | $5.69093e-06$ | -3.318 | $4.42491e-05$ | -2.840
$32$ | $1.24768e-07$ | -4.111 | $4.40714e-07$ | -3.691 | $4.63001e-06$ | -3.257
$64$ | $6.57969e-09$ | -4.245 | $2.86670e-08$ | -3.942 | $4.05012e-07$ | -3.515
$128$ | $2.95557e-10$ | -4.477 | $1.63672e-09$ | -4.131 | $3.12241e-08$ | -3.697
Table 2: Convergence study for the flow in a pipe with annular cross-section at $\text{Re}=1$ and $\text{Ma}=1e-3$. axial velocity error; solution polynomial degree: $p=3$. Grid | $L^{1}$ | Rate | $L^{2}$ | Rate | $L^{\infty}$ | Rate
---|---|---|---|---|---|---
$4$ | $4.54267e-05$ | - | $9.55211e-05$ | - | $4.70479e-04$ | -
$8$ | $4.47553e-06$ | -3.343 | $9.58936e-06$ | -3.316 | $6.31346e-05$ | -2.898
$16$ | $1.88611e-07$ | -4.569 | $5.68142e-07$ | -4.077 | $5.14343e-06$ | -3.618
$32$ | $5.97100e-09$ | -4.981 | $2.34227e-08$ | -4.600 | $2.75719e-07$ | -4.221
$64$ | $1.85738e-10$ | -5.007 | $7.90368e-10$ | -4.889 | $1.13276e-08$ | -4.605
$128$ | $5.76894e-12$ | -5.009 | $2.49773e-11$ | -4.984 | $3.99322e-10$ | -4.826
Table 3: Convergence study for the flow in a pipe with annular cross-section
at $\text{Re}=1$ and $\text{Ma}=1e-3$. axial velocity error; solution
polynomial degree: $p=4$.

In this section, we use the method of manufactured solutions (MMS) to verify
the accuracy and correct implementation of the proposed boundary conditions.
We use a pipe with an annular section and length $1$. Similarly to Dalcin et
al. (2019), this case is considered for two reasons. First, it has an analytic
solution for incompressible flow that cannot be represented in polynomial
space Rosenhead (1963). Second, it allows exercising the high-order mesh
capabilities and better approximate the circular geometry of the pipe. The
following solution is used for the axial velocity:

$\mathcal{U}_{1}=\frac{G}{4\mu}\left((R_{1}^{2}-r^{2})+(R_{2}^{2}-R_{1}^{2})\frac{\log(r/R_{1})}{\log(R_{2}/R_{1})}\right),$
(52)

where $R_{o}=0.5$, $R_{o}/R_{i}=4$, and $G/\mu=1$. We consider uniform density
and temperature, and zero nonaxial velocities. The parameters used are
$\text{Re}=1$, $\text{Ma}=1e-3$ and $\alpha=1$. No-slip adiabatic wall
boundary conditions are used on the outer and inner cylinder whereas, periodic
boundary conditions are used on the remaining boundaries. We use Mathematica
to compute the source term of the Eulerian model Inc. (2020).

The error calculation uses the following discrete norms

$\begin{split}&\text{Discrete
}L^{1}:\|\mathbfcal{U}\|_{L^{1}}=\|\Omega\|^{-1}\sum\limits_{k=1}^{K}\mathbf{1}_{N_{k}}^{\top}\mathcal{P}_{k}\mathbfcal{J}_{k}\text{abs}\left(\mathbfcal{U}_{k}\right),\\\
&\text{Discrete
}L^{2}:\|\mathbfcal{U}\|_{L^{2}}^{2}=\|\Omega\|^{-1}\sum\limits_{k=1}^{K}\mathbfcal{U}_{k}^{\top}\mathcal{P}_{k}\mathbfcal{J}_{k}\mathbfcal{U}_{k},\\\
&\text{Discrete
}L^{\infty}:\|\mathbfcal{U}\|_{L^{\infty}}=\max\limits_{k=1\dots
K}\text{abs}\left(\mathbfcal{U}_{k}\right),\end{split}$ (53)

where $\|\Omega\|$ represents the volume of the computational domain,
$\mathbfcal{J}_{k}$ is the metric Jacobian of the curvilinear transformation
from physical space to computational space of the $k$-th hexahedral element,
and $K$ is the total number of non-overlapping hexahedral elements in the
mesh.

The results of the convergence study are shown in Tables 1, 2 and 3, where the
first column reprisents the number of elements in the radial and angular
coordinates. We observe that the computed order of accuracy is approximately
$(p+1)$.

###  6.2 Spinning sphere: verification of entropy conservation

(a) Overview of the spinning sphere.

(b) A closer look at the velocity field near the spinning sphere.

Figure 3: Spinning sphere enclosed in a cubical box; Eulerian model with
$\text{Re}=1$, and $\text{Ma}=0.05$. The magnitude of the pointwise velocity
is used to scale the arrow glyphs vectors and to color the surface of the
sphere.

(a) Evolution of the time derivative of the entropy function and the
dissipation term.

(b) Evolution of the sum of the time derivative of the entropy function and
the dissipation term.

Figure 4: Verification of the entropy conservation for the spinning sphere
enclosed in a cubical box; Eulerian model with $\text{Re}=1$, and
$\text{Ma}=0.05$.

To verify the entropy conservation property of the new boundary conditions
without IP term, we simulate a spinning sphere enclosed in a cubic domain. The
sphere rotates at a constant angular velocity around a unit vector given by
$\hat{a_{r}}=\frac{a_{r}}{\|a_{r}\|}$, where $a_{r}=(1,1,1)^{\top}$. The
domain size is $2\times 2\times 2$ with a sphere of diameter $D=0.6$ located
at the center of it. Solid wall boundary conditions are imposed on the sphere
surface and all six faces of the cubic box. The mesh is composed of $4,374$
hexahedral elements. The sphere surface is discretized with quadratic boundary
elements. A solution polynomial degree $p=5$ is used. We run the Eulerian
model with $\text{Re}=1$, and $\text{Ma}=0.05$. In Figure 3, the velocity
vector field near the surface of the sphere is shown. The sphere is also
colored based on the module of its pointwise velocity vector.

Figure 4(a) shows the time derivative of the entropy function of the
semidiscretization of the Eulerian model (30),
$\frac{d}{dt}\mathbf{1}^{\top}\widehat{\mathcal{P}}\mathbfcal{S}$, and the
dissipation term, $\mathbf{DT}$, as a function of time, $t$ (see Equation
(36)). As shown in Figure 4(b), the two terms cancel out, up to machine
precision. Therefore, the procedure proposed to impose the boundary conditions
is entropy conservative if the IP term in (40) is set to zero. This simulation
is a strong numerical verification of what is proven in Section 4.4.

###  6.3 Laminar flow around a cylinder $\text{Re}=40$

\begin{overpic}[width=238.49231pt,tics=5]{figures/vis_cyl_domain.png} \put(50.0,1.0){$x_{1}$} \put(0.0,27.5){$x_{2}$} \put(40.0,54.0){\rotatebox{0.0}{Inviscid slip wall}} \par \put(-9.0,25.0){\rotatebox{90.0}{Inlet}} \put(114.0,22.0){\rotatebox{90.0}{Outlet}} \put(27.5,29.0){$D$} \put(29.0,28.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\vector(1,0){4.5}} \put(29.0,28.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\vector(-1,0){4.5}} \put(65.0,-3.0){\rotatebox{0.0}{$30D$}} \put(74.0,-1.75){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\vector(1,0){24.0}} \put(64.0,-1.75){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\vector(-1,0){34.5}} \put(15.0,-3.0){\rotatebox{0.0}{$10D$}} \put(24.0,-1.75){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\vector(1,0){4.5}} \put(14.0,-1.75){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\vector(-1,0){6.0}} \put(103.0,22.0){\rotatebox{90.0}{$20D$}} \put(105.0,30.5){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\vector(0,1){20.0}} \put(105.0,21.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\vector(0,-1){15.0}} \par\end{overpic} Figure 5: Cylinder domain with $D=1$. \begin{overpic}[width=411.93767pt,tics=5]{figures/vis_cyl_40_p7_wake_bw.png} \put(48.5,33.5){\rotatebox{0.0}{\small\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}Eulerian}} \put(51.0,24.75){\rotatebox{0.0}{\small\color[rgb]{0.4765625,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.4765625,0,0}CNS}} \end{overpic} Figure 6: Wake region of the flow past a cylinder at $\text{Re}=40$ and $\text{Ma}=0.07$. The background is colored by the magnitude of the pointwise velocity vector. | $L_{\text{bubble}}$ | $C_{D}$
---|---|---
Park et al. (1998) | $2.25$ | $1.51$
Ye et al. (1999) | $2.27$ | $1.52$
Fornberg (1980) | $2.24$ | $2.5$
Dennis and Chang (1970) | $2.35$ | $1.52$
Calhoun (2002) | $2.18$ | $1.62$
Russell and Wang (2003) | $2.29$ | $1.6$
SSDC CNS | $2.254$ | $1.608$
SSDC Eulerian | $2.256$ | $1.609$
Table 4: Bubble length, $L_{\text{bubble}}$, and drag coefficient, $C_{D}$.

(a) Leading edge.

(b) Trailing edge.

Figure 7: The density profile for the flow past a cylinder at $\text{Re}=40$
and $\text{Ma}=0.07$. \begin{overpic}[width=433.62pt,trim=0.0pt 0.0pt
-5.01874pt 0.0pt,clip]{figures/cp_40.pdf}
\put(89.75,52.25){\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{park1998numerical}{\@@citephrase{(}}{\@@citephrase{)}}}}
\end{overpic} Figure 8: Wall pressure coefficient for the flow past a cylinder
at $\text{Re}=40$ and $\text{Ma}=0.07$. The reference data are obtained from
Park et al. (1998). \begin{overpic}[width=433.62pt,trim=0.0pt 0.0pt -5.01874pt
0.0pt,clip]{figures/vort_40.pdf}
\put(89.75,10.75){\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{park1998numerical}{\@@citephrase{(}}{\@@citephrase{)}}}}
\end{overpic} Figure 9: Vorticity at the wal for the flow past a cylinder at
$\text{Re}=40$ and $\text{Ma}=0.07$. The reference data was obtained from Park
et al. (1998).

In this section, we present the two-dimensional flow around a cylinder of
diameter $D=1$ at $\text{Re}=40$, and $\text{Ma}=0.07$. Figure 5 shows the
computational domain and summarizes the boundary conditions. The initial
condition is a uniform flow in the $x_{1}$ direction. The case is run until
convergence to a steady-state. Specifically, we stop the simulation when the
residual reaches $10^{-14}$ and plateau around that value. Entropy stable
adiabatic no-slip wall boundary conditions are enforced on the cylinder
surface. Far-field boundary conditions are used for the remaining boundaries.
The mesh is composed of $1,140$ elements with quadratic edges for representing
the surface of the cylinder. We use a solution polynomial degree $p=7$.
Therefore, the total number of degrees of freedom is $72,960$. The solution
computed with this setup is denoted here as a “numerically converged solution”
in the sense that if we increase the order of accuracy of the method, $p$, or
refine the mesh, the difference between two consecutive numerical solutions
for all the primitive variables is machine precision.

First, we compute the length of the recirculation bubble in the wake region.
Figure 6 shows the streamlines of the circulation bubble downstream of the
sphere. The background is the magnitude of the pointwise velocity vector. The
top part of the plot corresponds to the solution computed with the Eulerian
model, whereas the bottom part is the solution obtained using the CNS model.
Qualitatively, the two solutions look identical. Furthermore, both models show
a wake length of $2.25$. The first column in Table 4 compares the length of
the bubble with several results published in the literature. The agreement is
good. The second column of Table 4 compares the drag coefficient, $C_{D}$. The
values obtained with both models are very close to each other and agree well
with that reported by Calhoun (2002) and Russell and Wang (2003).

Next, we present some measurements in the boundary layer of the cylinder wall.
Figure 7 shows the density profile for both models, where for the Eulerian
model, we enforce a zero normal density gradient at the wall. Furthermore, for
both models, we use adiabatic solid wall boundary conditions. Because in our
simulations, we use a non-dimensional formalism, the difference between the
two solutions is not negligible, especially near the leading edge; see Figure
7(a). The difference observed is accumulated at the leading edge and remains
consistent along with the tailwind, as shown in Figure 7(b).

Figure 8 shows the wall pressure coefficient, $C_{p}$. The curves obtained
with the Eulerian model and the CNS model are very close to each other and
follow relatively well the DNS results published by Park et al. (1998). A
similar comparison is made in Figure 9 for the vorticity at the wall. Again,
the curves obtained with the Eulerian model and the CNS model are very close
to each other and follow well the DNS results except around $50^{\circ}$ where
they both undershoot the reference results.

###  6.4 Blast-wave

Figure 10: Density distribution for the blastwave at time $t=0.01$ computed
with the Eulerian and the CNS models. The base figure uses a grid of $32,768$
elements.

Svärd (2018) presents the one-dimensional blastwave as an example where the
implementation of the CNS model converges to an unphysical solution. Here, we
investigate the same test case. We use $\text{Re}=10$ and $\text{Ma}=0.07$,
and the same initial conditions as in Svärd (2018). All the simulations are
performed with $p=1$. We present the results for a sequence of five uniform
nested grids, with respectively $2,048$, $4,096$, $8,192$, $16,384$, and
$32,768$ elements. The solutions obtained with the finest grid are denoted
here as “numerically converged solutions”. Solid wall adiabatic no-slip
boundary conditions with zero density gradient for the Eulerian model are
imposed at the left and right boundaries of the domain.

Svärd (2018) observes that since the density flux is not bounded, the solution
can converge to a state where the wall dissipates density. This is a non-
physical solution and is computed with a non-entropy stable implementation of
the CNS model Svärd (2018). Our discretization of the Eulerian and CNS models
are entropy stable. As we will see shortly, we are unable to reproduce the
“lagging” effect observed in Svärd (2018).

The main plot in Figure 10 compares the density profiles computed with the
Eulerian and the CNS models on the finest grid with $8,192$ elements. We can
easily observe that the two curves do not lie on top of each other.
Specifically, we notice a clear difference in the region near the maximum
peak, whose position along the $x$ axis also appears to be influenced by the
model. Furthermore, differently from Svärd (2018), we do not notice any
difference at the left and right boundaries of the domain between the solution
computed with the Eulerian and the CNS models. In Figure 10, we also report
two zoom views of the maximum peaks. There, we plot the density profiles
obtained with the five nested grids. We can observe how the solutions approach
the “numerically converged solution”. In particular, it appears that the
solution obtained with the Eulerian model reaches the “converged state” faster
than the solution obtained with the CNS model. Furthermore, we notice that the
value of the maximum peak obtained with the Eulerian model is slightly smaller
than the same value computed by the CNS model.

###  6.5 Supersonic flow around a cylinder

\begin{overpic}[width=238.49231pt,tics=5]{figures/cylinder_mesh.pdf}
\put(19.0,32.0){$x_{1}$} \put(9.0,40.0){$x_{2}$}
\put(40.0,18.5){\rotatebox{90.0}{\small{Cylinder}}}
\put(34.0,29.0){$\mathrm{D}$}
\put(-4.75,11.25){\rotatebox{90.0}{${\mathcal{U}}=\left(\mathcal{U}_{\infty,~{}0}\right)$}}
\put(100.0,22.0){\rotatebox{90.0}{$6\mathrm{D}$}}
\put(68.0,47.0){\rotatebox{0.0}{Inviscid slip wall}}
\par\put(65.0,-3.0){\rotatebox{0.0}{$30\mathrm{D}$}}
\put(18.0,-3.0){\rotatebox{0.0}{$10\mathrm{D}$}} \par\put(85.0,24.0){$p=4$}
\put(85.0,40.0){$p=3$} \put(85.0,6.0){$p=3$} \put(10.0,7.0){$p=2$}
\put(10.0,7.0){$p=2$} \put(12.0,16.0){$p=2$} \end{overpic} Figure 11:
Computational domain, $hp-$nonconforming grid, and solution polynomial degree
distribution, $p$, for the simulation of the supersonic flow past a circular
cylinder. The underlying figure was adapted from Parsani et al. (2020).

(a) Density, $\rho$.

(b) Velocity component in the $x_{1}$ direction, $\mathcal{U}_{1}$.

Figure 12: Supersonic flow past a cylinder enclosed between two plates at
$\text{Re}=10,000$ and $\text{Ma}=3.5$. Figure 13: $\theta$-$\beta-\text{Ma}$
plot for the flow past a cylinder at $\text{Re}=10,000$, and $\text{Ma}=3.5$.
\begin{overpic}[width=411.93767pt,tics=5]{figures/normalized_distance.png}
\definecolor{cnscolor}{HTML}{7a0000}
\put(39.0,22.0){\rotatebox{0.0}{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}SSDC
Eulerian $\Delta=0.299$}}
\put(39.0,16.5){\rotatebox{0.0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}Analytic
$\Delta=0.293$}}
\put(39.0,11.0){\rotatebox{0.0}{\color[rgb]{0.4765625,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.4765625,0,0}SSDC
CNS $\Delta=0.308$}} \end{overpic} Figure 14: Normalized distance between the
shock and the cylinder leading edge for the flow past a cylinder at
$\text{Re}=10,000$, and $\text{Ma}=3.5$. The background is colored by the
pointwise velocity magnitude.

In this last test case, we present some results for the flow past a two-
dimensional cylinder enclosed between two solid walls Boukharfane et al.
(2018). The similarity parameters are $\text{Re}=10,000$ and $\text{Ma}=3.5$.
The initial condition is a uniform flow with unit density, $\rho$, unit
temperature, $T$, and only a non-zero unit velocity component in the $x_{1}$
direction, $\mathcal{U}_{1}$. The flow is computed with a non-uniform
distribution of the solution polynomial degree. Therefore, the solution is
computed using the entropy stable $p$-nonconforming interface technology
Fernández et al. (2020). The solution polynomial degree distribution is shown
in Figure 11, and it is the same setup used in Parsani et al. (2020). Entropy
stable adiabatic no-slip wall boundary conditions are enforced on the cylinder
surface. Inviscid (slip) wall boundary conditions are imposed on the top and
bottom horizontal boundaries, whereas far-field boundary conditions are used
for the (left) inlet and the (right) outlet boundaries. The mesh consists of
$5,067$ quadrangles. The solution computed with this setup is denoted here as
a “numerically converged solution” in the sense that if we increase the order
of accuracy of the method, $p$, or we refine the mesh, the difference in the
solution for all the primitive variables is machine precision. Figure 12 shows
the contour plot of the density, $\rho$, and the velocity component in the
$x_{1}$ direction, $\mathcal{U}_{1}$, of the developed solution.

A quantitative analysis of the numerical results can be performed using the
oblique shock wave theory Shapiro (1953). At an oblique shock, the flow
changes direction, and there are three directions of interest: the upstream
and downstream flow directions and the shock wave direction itself. The most
useful relation of oblique shock wave theory is the one providing the
deflection angle, $\theta$, as a function of the shock wave angle, $\beta$,
and local Mach number, Ma Liepmann and Roshko (2001):

$\tan(\theta)=2\cot(\beta)\left[\frac{\text{Ma}^{2}\sin^{2}(\beta)-1}{\text{Ma}^{2}(\gamma+\cos(2\beta))+2}\right].$
(54)

The pair $\theta-\beta$ obtained by postprocessing the solutions computed with
the Eulerian and CNS models are shown in Figure 13. We can see that the
results for both models are in good agreement with the theoretical curve.
Finally, we report the normalized distance between the shock and the leading
edge of the cylinder, as shown in Figure 14. The value of the latter quantity
is $\Delta=0.299$ and $\Delta=0.308$ for the Eulerian and CNS models,
respectively. These values are in good agreement with the analytical value,
which reads $\Delta=0.293$ Boiron et al. (2009).

##  7 Conclusion

Guided by the entropy stability analysis, we have derived solid wall boundary
conditions that preserve the entropy stability of the Eulerian model for
viscous compressible fluids proposed by Svärd (2018). In that context, using
summation-by-parts operators and the simultaneous-approximation term
technique, we have constructed entropy conservative and entropy stable solid
wall boundary conditions for the semidiscretized system, which mimics the
continuous entropy analysis results. The proposed boundary conditions have
been validated in terms of accuracy, entropy conservation, and entropy
stability for a set of test cases of increasing complexity. The numerical
results obtained with the Eulerian model have been compared with the solutions
computed using the classic compressible Navier–Stokes equations also
discretized with an entropy stable methodology constructed using summation-by-
parts operators. Differences and similarities have then been highlighted. This
study is a first attempt to understand better the effects of the viscous flux
term introduced in the conservation of mass equation.

  * Data Availability

The datasets generated during and/or analysed during the current study are
available in the Zenodo repository, http://doi.org/10.5281/zenodo.5041436

###### Acknowledgements.

The research reported in this paper was funded by King Abdullah University of
Science and Technology. We are thankful for the computing resources of the
Supercomputing Laboratory and the Extreme Computing Research Center at King
Abdullah University of Science and Technology.

## References

  * Abhyankar et al. (2018) Abhyankar S, Brown J, Constantinescu EM, Ghosh D, Smith BF, Zhang H (2018) Petsc/ts: A modern scalable ode/dae solver library. arXiv preprint arXiv:180601437
  * Balay et al. (2021) Balay S, Abhyankar S, Adams MF, Brown J, Brune P, Buschelman K, Dalcin L, Dener A, Eijkhout V, Gropp WD, Karpeyev D, Kaushik D, Knepley MG, May DA, McInnes LC, Mills RT, Munson T, Rupp K, Sanan P, Smith BF, Zampini S, Zhang H, Zhang H (2021) PETSc users manual. Tech. Rep. ANL-95/11 - Revision 3.15, Argonne National Laboratory
  * Berg and Nordström (2011) Berg J, Nordström J (2011) Stable robin solid wall boundary conditions for the Navier–Stokes equations. Journal of Computational Physics 230(19):7519–7532
  * Bogacki and Shampine (1989) Bogacki P, Shampine L (1989) A 3(2) pair of runge - kutta formulas. Applied Mathematics Letters 2(4):321–325
  * Boiron et al. (2009) Boiron O, Chiavassa G, Donat R (2009) A high-resolution penalization method for large mach number flows in the presence of obstacles. Computers & fluids 38(3):703–714
  * Boukharfane et al. (2018) Boukharfane R, Ribeiro FHE, Bouali Z, Mura A (2018) A combined ghost-point-forcing/direct-forcing immersed boundary method (ibm) for compressible flow simulations. Computers & Fluids 162:91–112
  * Calhoun (2002) Calhoun D (2002) A cartesian grid method for solving the two-dimensional streamfunction-vorticity equations in irregular regions. Journal of computational physics 176(2):231–275
  * Carpenter et al. (1994) Carpenter MH, Gottlieb D, Abarbanel S (1994) Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: Methodology and application to high-order compact schemes. Journal of Computational Physics 111(2):220–236
  * Carpenter et al. (2014) Carpenter MH, Fisher TC, Nielsen EJ, Frankel SH (2014) Entropy stable spectral collocation schemes for the Navier–Stokes equations: Discontinuous interfaces. SIAM Journal on Scientific Computing 36(5):B835–B867
  * Carpenter et al. (2015) Carpenter MH, Parsani M, Fisher TC, Nielsen EJ (2015) Entropy stable staggered grid spectral collocation for the Burgers’ and compressible Navier-Stokes equations. NASA/TM–2015–218990
  * Chandrashekar (2013) Chandrashekar P (2013) Kinetic energy preserving and entropy stable finite volume schemes for compressible Euler and Navier-Stokes equations. Communications in Computational Physics 14(5):1252–1286
  * Dafermos et al. (2005) Dafermos CM, Dafermos CM, Dafermos CM, Mathématicien G, Dafermos CM, Mathematician G (2005) Hyperbolic conservation laws in continuum physics, vol 3. Springer
  * Dalcin et al. (2019) Dalcin L, Rojas D, Zampini S, Fernández DCDR, Carpenter MH, Parsani M (2019) Conservative and entropy stable solid wall boundary conditions for the compressible Navier–Stokes equations: Adiabatic wall and heat entropy transfer. Journal of Computational Physics 397:108775
  * Dennis and Chang (1970) Dennis S, Chang GZ (1970) Numerical solutions for steady flow past a circular cylinder at reynolds numbers up to 100. Journal of Fluid Mechanics 42(3):471–489
  * Dolejší and Svärd (2021) Dolejší V, Svärd M (2021) Numerical study of two models for viscous compressible fluid flows. Journal of Computational Physics 427:110068
  * Fernández et al. (2020) Fernández DCDR, Carpenter MH, Dalcin L, Fredrich L, Winters AR, Gassner GJ, Parsani M (2020) Entropy-stable p-nonconforming discretizations with the summation-by-parts property for the compressible Navier–Stokes equations. Computers & Fluids 210:104631
  * Fisher (2012) Fisher TC (2012) High-order $l^{2}$ stable multi-domain finite difference method for compressible flows. PhD thesis, Purdue University
  * Fornberg (1980) Fornberg B (1980) A numerical study of steady viscous flow past a circular cylinder. Journal of Fluid Mechanics 98(4):819–855
  * Geuzaine and Remacle (2009) Geuzaine C, Remacle JF (2009) Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities. International journal for numerical methods in engineering 79(11):1309–1331
  * Inc. (2020) Inc WR (2020) Mathematica, Version 12.1. URL https://www.wolfram.com/mathematica, champaign, IL, 2020
  * Knepley and Karpeev (2009) Knepley MG, Karpeev DA (2009) Mesh algorithms for PDE with Sieve I: Mesh distribution. Scientific Programming 17(3):215–230, DOI 10.3233/SPR-2009-0249
  * Liepmann and Roshko (2001) Liepmann HW, Roshko A (2001) Elements of gasdynamics. Courier Corporation
  * Nordström and Carpenter (1999) Nordström J, Carpenter MH (1999) Boundary and interface conditions for high-order finite-difference methods applied to the Euler and Navier–Stokes equations. Journal of Computational Physics 148(2):621–645
  * Park et al. (1998) Park J, Kwon K, Choi H (1998) Numerical solutions of flow past a circular cylinder at reynolds numbers up to 160. KSME international Journal 12(6):1200–1205
  * Parsani et al. (2015a) Parsani M, Carpenter MH, Nielsen EJ (2015a) Entropy stable discontinuous interfaces coupling for the three-dimensional compressible Navier–Stokes equations. Journal of Computational Physics 290:132–138
  * Parsani et al. (2015b) Parsani M, Carpenter MH, Nielsen EJ (2015b) Entropy stable wall boundary conditions for the three-dimensional compressible Navier–Stokes equations. Journal of Computational Physics 292:88–113
  * Parsani et al. (2020) Parsani M, Boukharfane R, Nolasco IR, Fernández DCDR, Zampini S, Hadri B, Dalcin L (2020) High-order accurate entropy-stable discontinuous collocated galerkin methods with the summation-by-parts property for compressible cfd frameworks: Scalable ssdc algorithms and flow solver. Journal of Computational Physics 424:109844
  * Rosenhead (1963) Rosenhead L (1963) Laminar boundary layers; an account of the development, structure, and stability of laminar boundary layers in incompressible fluids, together with a description of the associated experimental techniques., 1st edn. Oxford [England] Clarendon Press
  * Russell and Wang (2003) Russell D, Wang ZJ (2003) A cartesian grid method for modeling multiple moving objects in 2d incompressible viscous flow. Journal of Computational Physics 191(1):177–205
  * Schneiderbauer and Krieger (2013) Schneiderbauer S, Krieger M (2013) What do the Navier–Stokes equations mean? European Journal of Physics 35(1):015020
  * Shapiro (1953) Shapiro AH (1953) The Dynamics and Thermodynamics of Compressible Fluid Flow. No. 1 in Engineering special collection, Ronald Press Company
  * Svärd (2015) Svärd M (2015) Weak solutions and convergent numerical schemes of modified compressible Navier–Stokes equations. Journal of Computational Physics 288:19–51
  * Svärd (2018) Svärd M (2018) A new Eulerian model for viscous and heat conducting compressible flows. Physica A: Statistical Mechanics and its Applications 506:350–375
  * Svärd and Nordström (2008) Svärd M, Nordström J (2008) A stable high-order finite difference scheme for the compressible Navier–Stokes equations: no-slip wall boundary conditions. Journal of Computational Physics 227(10):4805–4824
  * Svärd and Nordström (2014) Svärd M, Nordström J (2014) Review of summation-by-parts schemes for initial–boundary-value problems. Journal of Computational Physics 268:17–38
  * Tadmor (2003) Tadmor E (2003) Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems. Acta Numerica 12:451–512
  * Ye et al. (1999) Ye T, Mittal R, Udaykumar H, Shyy W (1999) An accurate cartesian grid method for viscous incompressible flows with complex immersed boundaries. Journal of computational physics 156(2):209–240

##  Appendix A Derivation of the time derivative of the entropy

In this section, we provide the derivation of equation (36). To obtain the
expression of the time derivative of the entropy function, $\mathbfcal{S}$, we
multiply (34a) and (34b) by $\mathbf{w}^{\top}\mathcal{P}$ and
$\left(\left[\nu\frac{dq}{dw}\right]\Theta_{x_{j}}\right)^{\top}\mathcal{P}$,
respectively. This step yields

$\displaystyle\begin{split}&\frac{d}{dt}\mathbf{1}^{\top}\widehat{\mathcal{P}}\mathbfcal{S}+\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\widehat{\mathcal{B}}_{x_{i}}\mathbf{F}_{x_{i}}-\mathbf{w}^{\top}\mathcal{P}\mathcal{D}_{x_{i}}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}\\\
&=\mathbf{w}^{\top}\mathcal{P}_{x_{j},x_{k}}\left(\mathbf{g}^{(b),q}_{x_{i}}+\mathbf{g}^{(In),q}_{x_{i}}\right),\end{split}$
(55a)
$\displaystyle\begin{split}&\left(\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}\right)^{\top}\mathcal{P}\mathbf{\Theta}_{x_{i}}-\left(\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}\right)^{\top}\mathcal{P}\mathcal{D}_{x_{i}}\mathbf{w}\\\
&=\left(\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}\right)^{\top}\mathcal{P}_{x_{j},x_{k}}\left(\mathbf{g}^{(b),\Theta}_{x_{i}}+\mathbf{g}^{(In),\Theta}_{x_{i}}\right).\end{split}$
(55b)

Using the properties of SBP operators, we can write

$\mathbf{w}^{\top}\mathcal{P}\mathcal{D}_{x_{i}}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}=\mathbf{w}^{\top}\mathcal{P}\mathcal{P}^{-1}_{x_{i}}\mathcal{B}_{x_{i}}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}-\mathbf{w}^{\top}\mathcal{P}\mathcal{P}^{-1}_{x_{i}}\mathcal{Q}^{\top}_{x_{i}}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}.$
(56)

Since $\mathcal{P}$ is diagonal, we have

$\mathcal{P}\mathcal{P}^{-1}_{x_{i}}\mathcal{Q}^{\top}_{x_{i}}=\mathcal{P}\mathcal{Q}^{\top}_{x_{i}}\mathcal{P}^{-1}_{x_{i}}=\mathcal{P}\mathcal{D}^{\top}_{x_{i}}=\mathcal{D}^{\top}_{x_{i}}\mathcal{P}.$
(57)

Thus, Equation (55a) is recast in the following form:

$\displaystyle\begin{split}&\frac{d}{dt}\mathbf{1}^{\top}\widehat{\mathcal{P}}\mathbfcal{S}+\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\widehat{\mathcal{B}}_{x_{i}}\mathbf{F}_{x_{i}}+\left(\mathcal{D}_{x_{i}}\mathbf{w}\right)^{\top}\mathcal{P}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}\\\
&=\mathbf{w}^{\top}\mathcal{P}_{x_{j},x_{k}}\mathcal{B}_{x_{i}}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}+\mathbf{w}^{\top}\mathcal{P}_{x_{j},x_{k}}\left(\mathbf{g}^{(b),q}_{x_{i}}+\mathbf{g}^{(In),q}_{x_{i}}\right).\end{split}$
(58)

Additionally, because $\left[\nu\frac{dq}{dw}\right]$ is SPD, we have

$\left(\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}\right)^{\top}\mathcal{P}\mathbf{\Theta}_{x_{i}}=\left(\left[\nu\frac{dq}{dw}\right]^{\frac{1}{2}}\mathbf{\Theta}_{x_{j}}\right)^{\top}\mathcal{P}\left(\left[\nu\frac{dq}{dw}\right]^{\frac{1}{2}}\mathbf{\Theta}_{x_{i}}\right)=\left\|\left[\nu\frac{dq}{dw}\right]^{\frac{1}{2}}\mathbf{\Theta}_{x_{i}}\right\|^{2}_{\mathcal{P}}.$
(59)

Now, substituting (55b) into (58) gives the following expression

$\displaystyle\begin{split}&\frac{d}{dt}\mathbf{1}^{\top}\widehat{\mathcal{P}}\mathbfcal{S}+\left\|\left[\nu\frac{dq}{dw}\right]^{\frac{1}{2}}\mathbf{\Theta}_{x_{i}}\right\|^{2}_{\mathcal{P}}\\\
&=-\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\widehat{\mathcal{B}}_{x_{i}}\mathbf{F}_{x_{i}}\\\
&+\mathbf{w}^{\top}\mathcal{P}_{x_{j},x_{k}}\mathcal{B}_{x_{i}}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}\\\
&+\mathbf{w}^{\top}\mathcal{P}_{x_{j},x_{k}}\left(\mathbf{g}^{(b),q}_{x_{i}}+\mathbf{g}^{(In),q}_{x_{i}}\right)\\\
&+\left(\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}\right)^{\top}\mathcal{P}_{x_{j},x_{k}}\left(\mathbf{g}^{(b),\Theta}_{x_{i}}+\mathbf{g}^{(In),\Theta}_{x_{i}}\right).\end{split}$
(60)

By defining

$\mathbf{DT}=\left\|\left[\nu\frac{dq}{dw}\right]^{\frac{1}{2}}\mathbf{\Theta}_{x_{i}}\right\|^{2}_{\mathcal{P}},$
(61)

and

$\displaystyle\begin{split}\Xi&=-\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\widehat{\mathcal{B}}_{x_{i}}\mathbf{F}_{x_{i}}\\\
&+\mathbf{w}^{\top}\mathcal{P}_{x_{j},x_{k}}\mathcal{B}_{x_{i}}\left[\nu\frac{dq}{dw}\right]\mathbf{\Theta}_{x_{j}}\\\
&+\mathbf{w}^{\top}\mathcal{P}_{x_{j},x_{k}}\left(\mathbf{g}^{(b),q}_{x_{i}}+\mathbf{g}^{(In),q}_{x_{i}}\right)\\\
&+\left(\left[\nu\frac{dq}{dw}\right]\Theta_{x_{j}}\right)^{\top}\mathcal{P}_{x_{j},x_{k}}\left(\mathbf{g}^{(b),\Theta}_{x_{i}}+\mathbf{g}^{(In),\Theta}_{x_{i}}\right),\end{split}$
(62)

Equation (60) can be further simplified to obtain (36), i.e.,

$\displaystyle\frac{d}{dt}\mathbf{1}^{\top}\widehat{\mathcal{P}}\mathbfcal{S}+\mathbf{DT}=\Xi.$
(63)

##  Appendix B Change of variables matrix, $\frac{dq}{dw}(v)$

The Jacobian of the conservative variables, $q$, in terms of the entropy
variables, $w$, as a function of the primitive variables, $v$, is given by

$\frac{dq}{dw}(v)=\begin{bmatrix}\frac{\rho}{R}&\frac{\mathcal{U}_{1}\rho}{R}&\frac{\mathcal{U}_{2}\rho}{R}&\frac{\mathcal{U}_{3}\rho}{R}&\frac{dq}{dw}(v)_{1,5}\\\
\frac{\mathcal{U}_{1}\rho}{R}&\left(\mathcal{T}+\frac{\mathcal{U}_{1}^{2}}{R}\right)\rho&\frac{\mathcal{U}_{1}\mathcal{U}_{2}\rho}{R}&\frac{\mathcal{U}_{1}\mathcal{U}_{3}\rho}{R}&\frac{dq}{dw}(v)_{2,5}\\\
\frac{\mathcal{U}_{2}\rho}{R}&\frac{\mathcal{U}_{1}\mathcal{U}_{2}\rho}{R}&\left(\mathcal{T}+\frac{\mathcal{U}_{2}^{2}}{R}\right)\rho&\frac{\mathcal{U}_{2}\mathcal{U}_{3}\rho}{R}&\frac{dq}{dw}(v)_{3,5}\\\
\frac{\mathcal{U}_{3}\rho}{R}&\frac{\mathcal{U}_{1}\mathcal{U}_{3}\rho}{R}&\frac{\mathcal{U}_{2}\mathcal{U}_{3}\rho}{R}&\left(\mathcal{T}+\frac{\mathcal{U}_{3}^{2}}{R}\right)\rho&\frac{dq}{dw}(v)_{4,5}\\\
\frac{dq}{dw}(v)_{1,5}&\frac{dq}{dw}(v)_{2,5}&\frac{dq}{dw}(v)_{3,5}&\frac{dq}{dw}(v)_{4,5}&\frac{dq}{dw}(v)_{5,5}\end{bmatrix},$
(64)

where

$\displaystyle\frac{dq}{dw}(v)_{1,5}=\frac{\rho\left(-2R\mathcal{T}+\|\mathcal{U}\|_{L^{2}}^{2}+2\mathcal{T}c_{pr}\right)}{2R},$
(65)

$\displaystyle\frac{dq}{dw}(v)_{2,5}=\frac{\mathcal{U}_{1}\rho\left(\|\mathcal{U}\|_{L^{2}}^{2}+2\mathcal{T}c_{pr}\right)}{2R},$
(66)

$\displaystyle\frac{dq}{dw}(v)_{3,5}=\frac{\mathcal{U}_{2}\rho\left(\|\mathcal{U}\|_{L^{2}}^{2}+2\mathcal{T}c_{pr}\right)}{2R},$
(67)

$\displaystyle\frac{dq}{dw}(v)_{4,5}=\frac{\mathcal{U}_{3}\rho\left(\|\mathcal{U}\|_{L^{2}}^{2}+2\mathcal{T}c_{pr}\right)}{2R},$
(68)

$\displaystyle\frac{dq}{dw}(v)_{5,5}=\frac{\rho\left(\left(\|\mathcal{U}\|_{L^{2}}^{2}\right)^{2}+4\mathcal{T}c_{pr}\left(-R\mathcal{T}+\|\mathcal{U}\|_{L^{2}}^{2}+\mathcal{T}c_{pr}\right)\right)}{4R}.$
(69)

##  Appendix C Viscous boundary penalty

This section constructs the viscous boundary penalty term for an adiabatic
wall and a wall with heat entropy flux. We then show that the resulting
expressions lead to an entropy conservative and entropy stable penalty terms,
respectively.

Assuming a dissipative internal penalty term, $\mathcal{M}=0$, the RHS of (36)
becomes

$\frac{d}{dt}\mathbf{1}^{\top}\widehat{\mathcal{P}}\mathbfcal{S}+\mathbf{DT}=w^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\left(g^{(b,V),q}_{x_{i}}\right)+\Theta_{x_{i}}^{\top}\left(\nu\frac{dq}{dw}\right)\widehat{\mathcal{P}}_{x_{j},x_{k}}g^{(b),\Theta}_{x_{i}}.$
(70)

Define $\mathbf{W}=\text{diag}(\mathbf{w})$ and write

$\mathbf{w}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\mathbf{\chi}=\mathbf{1}^{\top}\mathbf{W}\widehat{\mathcal{P}}_{x_{j},x_{k}}\mathbf{\chi}=\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\mathbf{W}\mathbf{\chi},$
(71)

where
$\mathbf{W}\mathcal{P}_{x_{j},x_{k}}=\mathcal{P}_{x_{j},x_{k}}\mathbf{W}$,
since $\mathbf{W}$ is diagonal, and $\mathbf{\chi}$ is a vector on a node. We
use (71) for both terms in the RHS of (70). The first term reads

$\frac{1}{2}\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\widehat{\mathcal{B}}^{-}_{x_{i}}\mathbf{W}\left(\nu\frac{dq}{dw}\right)\Theta_{x_{i}}-\frac{1}{2}\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\widehat{\mathcal{B}}^{-}_{x_{i}}\mathbf{W}\left(\nu\frac{dq}{dw}\left(v^{(b,V)}\right)\right)\widetilde{\Theta}_{x_{i}},$
(72)

where both $\frac{dq}{dw}\left(v^{(b,V)}\right)$ and $\widetilde{\Theta}$ are
evaluated using the manufactured primitive state

$v^{(b,V)}=\left(\rho,-\mathcal{U}_{1}+2\mathcal{U}_{1}^{\text{wall}},-\mathcal{U}_{2}+2\mathcal{U}_{2}^{\text{wall}},-\mathcal{U}_{3}+2\mathcal{U}_{3}^{\text{wall}},\mathcal{T}\right)^{\top}.$
(73)

The kinematic viscosity, $\nu$, depends only on the density $\rho$. Thus, its
evaluation is the same for $v$ and $v^{(b,V)}$. The procedure for computing
$\widetilde{\Theta}_{x_{i}}$ is detailed in Dalcin et al. (2019) and is
summarized in the following equation. The manufactured entropy variables
gradient is defined as

$\widetilde{\Theta}_{x_{i}}=\frac{dw}{dv}\left(v^{(b,V)}\right)\tilde{\Pi}=\frac{dw}{dv}\left(v^{(b,V)}\right)\text{diag}(-1,1,1,1,-1)\Pi,$
(74)

where $\Pi$ is the gradient of primitive variables. The second term on the RHS
of (70) is written as

$\frac{1}{2}\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\widehat{\mathcal{B}}^{-}_{x_{i}}\mathbf{W}\left(\nu\frac{dq}{dw}\right)\Theta_{x_{i}}-\frac{1}{2}\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\widehat{\mathcal{B}}^{-}_{x_{i}}\mathbf{W}\left(v^{(b,V)}\right)\left(\nu\frac{dq}{dw}\right)\Theta_{x_{i}}.$
(75)

The second terms in (72) and (75) cancel out. On the other hand, the first
terms add up to

$\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}\widehat{\mathcal{B}}^{-}_{x_{i}}\mathbf{W}\left(\nu\frac{dq}{dw}\right)\Theta_{x_{i}}=\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}g(t).$
(76)

Therefore, we arrive at (49), i.e.,

$\frac{d}{dt}\mathbf{1}^{\top}\widehat{\mathcal{P}}\mathcal{S}+DT=\mathbf{1}^{\top}\widehat{\mathcal{P}}_{x_{j},x_{k}}g(t).$
(77)

For an adiabatic wall $g(t)=0$, and thus, the RHS is zero yielding an entropy
conservative term. For a wall with non-zero heat entropy flux, $g(t)$ is data,
and the RHS is bounded.

##  Appendix D Interior penalty term, $\mathcal{M}^{(b,V)}$

In this section, we compute the contribution to time derivative of the entropy
function, $\mathcal{S}$, from the interior penalty term $\mathcal{M}^{(b,V)}$.
Assuming a uniform hexahedral element (i.e., a cube) of length one and a wall
in the $x_{i}$ direction, we use the manufactured state of the vector of
primitive variables at the wall (46) to construct the jump in the entropy
variables in the normal direction to the wall, $\hat{n}$,

$\Delta w^{(b,V)}=(w-w^{(b,V)})\cdot\hat{n}^{\top},$ (78)

and the viscous boundary flux state

$\displaystyle
f^{(V)}(v^{(b,V)})=\left(\nu\frac{dq}{dw}(v^{(b,V)})\right)\Delta
w^{(b,V)}\cdot\hat{n}.$ (79)

Then, we use an average of the viscous flux of the numerical solution state at
the wall and the manufactured state

$\mathcal{M}^{(b,V)}=-\beta\frac{f^{(V)}(v)+f^{(V)}(v^{(b,V)})}{2},$ (80)

to compute the dissipative term with Mathematica Inc. (2020)

$w^{\top}\mathcal{M}^{(b,V)}=-\frac{2\beta\alpha\mu}{R\mathcal{T}^{2}}\|\Delta\mathcal{U}\|^{2}\left(\|\Delta\mathcal{U}\|^{2}+R\mathcal{T}\right),$
(81)

where $\Delta\mathcal{U}$ is the jump in the velocity of the numerical
solution state at the wall and the manufactured state. The term is negative
for any non-zero jump in velocity, and thus, is entropy dissipative.