Conservative and entropy stable solid wall boundary conditions for the compressible Navier-Stokes equations: Adiabatic wall and heat entropy transfer

TL;DR

This paper introduces a new method for imposing entropy conservative and stable solid wall boundary conditions in compressible Navier-Stokes simulations, ensuring stability and accuracy for complex flow scenarios.

Contribution

It generalizes previous boundary condition techniques using summation-by-parts and penalty flux methods, applicable across various discretization schemes.

Findings

The method ensures a semi-discrete entropy estimate for the entire domain.

Numerical tests demonstrate robustness and non-linear stability in 3D flows.

Compatible with multiple high-order discretization methods.

Abstract

We present a novel technique for the imposition of non-linear entropy conservative and entropy stable solid wall boundary conditions for the compressible Navier-Stokes equations in the presence of an adiabatic wall, or a wall with a prescribed heat entropy flow. The procedure relies on the formalism and mimetic properties of diagonal-norm, summation-by-parts, and simultaneous-approximation-term operators, and is a generalization of previous works on discontinuous interface coupling [1] and solid wall boundary conditions [2]. Using the method of lines, a semi-discrete entropy estimate for the entire domain is obtained when the proposed numerical imposition of boundary conditions are coupled with an entropy-conservative or entropy-stable discrete interior operator. The resulting estimate mimics the global entropy estimate obtained at the continuous level. The boundary data at the wall are weakly imposed using a penalty flux approach and a simultaneous-approximation-term technique for both the conservative variables and the gradient of the entropy variables. Discontinuous spectral collocation operators (mass lumped nodal discontinuous Galerkin operators), on high-order unstructured grids, are used for the purpose of demonstrating the robustness and efficacy of the new procedure for weakly enforcing boundary conditions. Numerical simulations confirm the non-linear stability of the proposed technique, with applications to three-dimensional subsonic and supersonic flows. The procedure described is compatible with any diagonal-norm summation-by-parts spatial operator, including finite element, finite difference, finite volume, discontinuous Galerkin, and flux reconstruction schemes.

Tables

Table 1. Table 1 : Convergence study for the flow in a pipe with annular section; p = 2 𝑝 2 p=2 ; error in the axial velocity.

Grid	$L^1$	Rate	$L^2$	Rate	$L^{\infty }$	Rate
4	8.77e-02	-	1.79e-01	-	5.69e-01	-
8	1.28e-02	-2.77	3.45e-02	-2.38	1.50e-01	-1.92
16	2.03e-03	-2.66	5.47e-03	-2.66	3.12e-02	-2.27
32	2.68e-04	-2.92	7.51e-04	-2.86	5.07e-03	-2.62
64	3.33e-05	-3.01	9.66e-05	-2.96	7.23e-04	-2.81

Table 2. Table 2 : Convergence study for the flow in a pipe with annular section; p = 3 𝑝 3 p=3 ; error in the axial velocity.

Grid	$L^1$	Rate	$L^2$	Rate	$L^{\infty }$	Rate
4	2.04e-02	-	3.87e-02	-	1.54e-01	-
8	2.71e-03	-2.91	5.68e-03	-2.77	3.17e-02	-2.28
16	2.15e-04	-3.65	5.69e-04	-3.32	4.43e-03	-2.84
32	1.25e-05	-4.10	4.40e-05	-3.69	4.63e-04	-3.26
64	6.65e-07	-4.23	2.87e-06	-3.94	4.05e-05	-3.51

Table 3. Table 3 : Convergence study for the flow in a pipe with annular section; p = 4 𝑝 4 p=4 ; error in the axial velocity.

Grid	$L^1$	Rate	$L^2$	Rate	$L^{\infty }$	Rate
4	4.54e-03	-	9.55e-03	-	4.70e-02	-
8	4.48e-04	-3.34	9.58e-04	-3.32	6.31e-03	-2.90
16	1.89e-05	-4.57	5.68e-05	-4.08	5.14e-04	-3.66
32	5.98e-07	-4.98	2.34e-06	-4.60	2.76e-05	-4.22
64	2.22e-08	-4.75	7.65e-08	-4.95	1.13e-06	-4.61

Table 4. Table 4 : Subsonic flow past a cylinder: mean drag coefficient and Strouhal number for different polynomial orders p 𝑝 p and uniform refinements r ​ l 𝑟 𝑙 rl , R ​ e ∞ = 300 𝑅 subscript 𝑒 300 Re_{\infty}=300 , M ∞ = 0.1 subscript 𝑀 0.1 M_{\infty}=0.1 . Reference values are also provided.

Poly. order	$rl=0$	$rl=1$	$rl=2$
$p=2$ $p=3$ $p=4$ $p=5$ [50]	${\overline{c}}_d$ $St$ 1.579 0.1995 1.439 0.1997 1.458 0.1998 1.414 0.1999 1.381 0.2000	${\overline{c}}_d$ $St$ 1.541 0.1998 1.439 0.1999 1.386 0.1999 1.382 0.1999 1.381 0.2000	${\overline{c}}_d$ $St$ 1.443 0.1999 1.383 0.1999 1.384 0.1999 1.386 0.1999 1.381 0.2000

Table 5. Table 5 : Subsonic flow past a sphere: mean drag coefficient c ¯ d subscript ¯ 𝑐 𝑑 \bar{c}_{d} , mean lift coefficient c ¯ l subscript ¯ 𝑐 𝑙 \bar{c}_{l} and its orientation in the x 1 = 0 subscript 𝑥 1 0 x_{1}=0 plane θ 𝜃 \theta , and Strouhal number S ​ t 𝑆 𝑡 St for different polynomial orders p 𝑝 p and uniform refinements l 𝑙 l , R ​ e ∞ = 300 𝑅 subscript 𝑒 300 Re_{\infty}=300 , M ∞ = 0.1 subscript 𝑀 0.1 M_{\infty}=0.1 . Reference values are also provided.

Poly. order	$rl=0$	$rl=1$	$rl=2$
$p=2$ $p=3$ $p=4$ $p=5$ [52]	${\overline{c}}_d$ ${\overline{c}}_l$ $\theta$ $St$ 1.125 0.095 -0.34 0.060 0.815 0.074 -0.18 0.120 0.688 0.069 0.14 0.140 0.650 0.066 0.39 0.140 0.656 0.069 - 0.137	${\overline{c}}_d$ ${\overline{c}}_l$ $\theta$ $St$ 0.843 0.076 -0.10 0.140 0.671 0.068 0.24 0.140 0.656 0.067 1.80 0.140 0.659 0.067 1.68 0.140 0.656 0.069 - 0.137	${\overline{c}}_d$ ${\overline{c}}_l$ $\theta$ $St$ 0.691 0.069 0.45 0.140 0.659 0.068 1.71 0.140 0.659 0.068 1.70 0.140 0.659 0.068 1.70 0.140 0.656 0.069 - 0.137

Equations

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

$$\mathcal{Q}={\left[\rho ,\rho {\mathcal{U}}_1,\rho {\mathcal{U}}_2,\rho {\mathcal{U}}_3,\rho \mathcal{E}\right]}^{\top },$$

$${\mathcal{F}}_{x_m}^{(I)}={\left[\rho {\mathcal{U}}_m,\rho {\mathcal{U}}_m{\mathcal{U}}_1+{\delta }_{m,1}\mathcal{P},\rho {\mathcal{U}}_m{\mathcal{U}}_2+{\delta }_{m,2}\mathcal{P},\rho {\mathcal{U}}_m{\mathcal{U}}_3+{\delta }_{m,3}\mathcal{P},\rho {\mathcal{U}}_m\mathcal{H}\right]}^{\top },$$

$${\mathcal{F}}_{x_m}^{(V)}={\left[0,{\tau }_{1,m},{\tau }_{2,m},{\tau }_{3,m},\sum _{i=1}^3{\tau }_{i,m}{\mathcal{U}}_i-\kappa \frac{\partial \mathcal{T}}{\partial x_m}\right]}^{\top },$$

$${\tau }_{i,j}=\mu \left(\frac{\partial {\mathcal{U}}_i}{\partial x_j}+\frac{\partial {\mathcal{U}}_j}{\partial x_i}-{\delta }_{i,j}\frac{2}{3}\sum _{n=1}^3\frac{\partial {\mathcal{U}}_n}{\partial x_n}\right),$$

$$\mathcal{H}=c_{\mathcal{P}}\mathcal{T}+\frac{1}{2}{\mathcal{U}}^{\top }\mathcal{U},\mathcal{P}=\rho R\mathcal{T},R=\frac{R_u}{M_w},$$

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

$$\mathcal{S}=-\rho s,$$

$$\bm{\mathcal{W}}=\left[\frac{\mathcal{H}-1/2\,{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\bm{\mathcal{U}}^{\mathrm{\top}}\bm{\mathcal{U}}}}}{\mathcal{T}}-s-\frac{{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\bm{\mathcal{U}}^{\mathrm{\top}}\bm{\mathcal{U}}}}}{\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]^{\mathrm{\top}}.$$

$$\frac{\partial \mathcal{S}}{\partial \mathcal{Q}}\frac{\partial {\mathcal{F}}_{x_m}^{(I)}}{\partial x_m}=\frac{\partial \mathcal{S}}{\partial \mathcal{Q}}\frac{\partial {\mathcal{F}}_{x_m}^{(I)}}{\partial \mathcal{Q}}\frac{\partial \mathcal{Q}}{\partial x_m}=\frac{\partial {\mathcal{F}}_{x_m}}{\partial \mathcal{Q}}\frac{\partial \mathcal{Q}}{\partial x_m}=\frac{\partial {\mathcal{F}}_{x_m}}{\partial x_m},m=1,2,3,$$

$$\frac{\partial \mathcal{Q}}{\partial \mathcal{W}}\frac{\partial \mathcal{W}}{\partial t}+\sum _{m=1}^3\frac{\partial {\mathcal{F}}_{x_m}^{(I)}}{\partial \mathcal{W}}\frac{\partial \mathcal{W}}{\partial x_m}=\sum _{m,j=1}^3\frac{\partial }{\partial x_m}\left({\mathsf{C}}_{m,j}\frac{\partial \mathcal{W}}{\partial x_j}\right),$$

$${\mathcal{F}}_{x_m}^{(V)}=\sum _{j=1}^3{\mathsf{C}}_{m,j}\frac{\partial \mathcal{W}}{\partial x_j}.$$

$${\mathcal{Q}}^{\top }=\frac{\partial \Phi }{\partial \mathcal{W}},{\left({\mathcal{F}}_{x_m}^{(I)}\right)}^{\top }=\frac{\partial {\Psi }_{x_m}}{\partial \mathcal{W}}.$$

$$\Phi ={\mathcal{W}}^{\top }\mathcal{Q}-\mathcal{S},$$

$${\Psi }_{x_m}={\mathcal{W}}^{\top }{\mathcal{F}}_{x_m}^{(I)}-{\mathcal{F}}_{x_m}.$$

$$\frac{\partial \mathcal{S}}{\partial \mathcal{Q}}\frac{\partial \mathcal{Q}}{\partial t}+\sum _{m=1}^3\frac{\partial \mathcal{S}}{\partial \mathcal{Q}}\frac{\partial {\mathcal{F}}_{x_m}^{(I)}}{\partial x_m}$$

$$\sum _{m=1}^3\left[\frac{\partial \mathcal{S}}{\partial t}+\frac{\partial {\mathcal{F}}_{x_m}}{\partial x_m}\right]$$

$$=\sum _{m=1}^3\left(\frac{\partial }{\partial x_m}\left({\mathcal{W}}^{\top }{\mathcal{F}}_{x_m}^{(V)}\right)-\sum _{j=1}^3{\left(\frac{\partial \mathcal{W}}{\partial x_m}\right)}^{\top }{\mathsf{C}}_{m,j}\frac{\partial \mathcal{W}}{\partial x_j}\right).$$

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\Omega }\mathcal{S}\mathrm{d}\Omega & \le \sum _{m=1}^3{\int }_{\Gamma }\left({\mathcal{W}}^{\top }{\mathcal{F}}_{x_m}^{(V)}-{\mathcal{F}}_{x_m}\right)n_{x_m}\mathrm{d}\Gamma -\sum _{m,j=1}^3{\int }_{\Omega }{\left(\frac{\partial \mathcal{W}}{\partial x_m}\right)}^{\top }{\mathsf{C}}_{m,j}\frac{\partial \mathcal{W}}{\partial x_j}\mathrm{d}\Omega \\ & =\sum _{m=1}^3{\int }_{\Gamma }\left({\mathcal{W}}^{\top }{\mathcal{F}}_{x_m}^{(V)}-{\mathcal{F}}_{x_m}\right)n_{x_m}\mathrm{d}\Gamma -DT,\end{array}$$

$$DT=\sum _{m,j=1}^3{\int }_{\Omega }{\left(\frac{\partial \mathcal{W}}{\partial x_m}\right)}^{\top }{\mathsf{C}}_{m,j}\frac{\partial \mathcal{W}}{\partial x_j}\mathrm{d}\Omega .$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\Omega }\mathcal{S}\mathrm{d}{x}_1\mathrm{d}{x}_2\mathrm{d}{x}_3=-DT$$

$$+{\int }_{x_1=0}\left[+{\mathcal{F}}_{x_1}-{\mathcal{W}}^{\top }\left({\mathsf{C}}_{1,1}\frac{\partial \mathcal{W}}{\partial x_1}+{\mathsf{C}}_{1,2}\frac{\partial \mathcal{W}}{\partial x_2}+{\mathsf{C}}_{1,3}\frac{\partial \mathcal{W}}{\partial x_3}\right)\right]\mathrm{d}{x}_2\mathrm{d}{x}_3$$

$$+{\int }_{x_1=1}\left[-{\mathcal{F}}_{x_1}+{\mathcal{W}}^{\top }\left({\mathsf{C}}_{1,1}\frac{\partial \mathcal{W}}{\partial x_1}+{\mathsf{C}}_{1,2}\frac{\partial \mathcal{W}}{\partial x_2}+{\mathsf{C}}_{1,3}\frac{\partial \mathcal{W}}{\partial x_3}\right)\right]\mathrm{d}{x}_2\mathrm{d}{x}_3$$

$$+{\int }_{x_2=0}\left[+{\mathcal{F}}_{x_2}-{\mathcal{W}}^{\top }\left({\mathsf{C}}_{2,1}\frac{\partial \mathcal{W}}{\partial x_1}+{\mathsf{C}}_{2,2}\frac{\partial \mathcal{W}}{\partial x_2}+{\mathsf{C}}_{2,3}\frac{\partial \mathcal{W}}{\partial x_3}\right)\right]\mathrm{d}{x}_1\mathrm{d}{x}_3$$

$$+{\int }_{x_2=1}\left[-{\mathcal{F}}_{x_2}+{\mathcal{W}}^{\top }\left({\mathsf{C}}_{2,1}\frac{\partial \mathcal{W}}{\partial x_1}+{\mathsf{C}}_{2,2}\frac{\partial \mathcal{W}}{\partial x_2}+{\mathsf{C}}_{2,3}\frac{\partial \mathcal{W}}{\partial x_3}\right)\right]\mathrm{d}{x}_1\mathrm{d}{x}_3$$

$$+{\int }_{x_3=0}\left[+{\mathcal{F}}_{x_3}-{\mathcal{W}}^{\top }\left({\mathsf{C}}_{3,1}\frac{\partial \mathcal{W}}{\partial x_1}+{\mathsf{C}}_{3,2}\frac{\partial \mathcal{W}}{\partial x_2}+{\mathsf{C}}_{3,3}\frac{\partial \mathcal{W}}{\partial x_3}\right)\right]\mathrm{d}{x}_1\mathrm{d}{x}_2$$

$$+{\int }_{x_3=1}\left[-{\mathcal{F}}_{x_3}+{\mathcal{W}}^{\top }\left({\mathsf{C}}_{3,1}\frac{\partial \mathcal{W}}{\partial x_1}+{\mathsf{C}}_{3,2}\frac{\partial \mathcal{W}}{\partial x_2}+{\mathsf{C}}_{3,3}\frac{\partial \mathcal{W}}{\partial x_3}\right)\right]\mathrm{d}{x}_1\mathrm{d}{x}_2.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\Omega }S\mathrm{d}{x}_1\mathrm{d}{x}_2\mathrm{d}{x}_3=-DT$$

$$+{\int }_{x_1=0}\left[{\mathcal{F}}_{x_1}-{\mathcal{W}}^{\top }\left({\mathsf{C}}_{1,1}\frac{\partial \mathcal{W}}{\partial x_1}+{\mathsf{C}}_{1,2}\frac{\partial \mathcal{W}}{\partial x_2}+{\mathsf{C}}_{1,3}\frac{\partial \mathcal{W}}{\partial x_3}\right)\right]\mathrm{d}{x}_2\mathrm{d}{x}_3.$$

$${\mathcal{F}}_{x_1}={\mathcal{W}}^{\top }{\mathcal{F}}_{x_1}^{(I)}-{\Psi }_{x_1}=-\rho s{\mathcal{U}}_{1}R,{\Psi }_{x_1}=-\rho {\mathcal{U}}_{1}R.$$

$$\mathtt{g}(t)=\kappa \frac{\partial \mathcal{T}}{\partial n}\frac{1}{\mathcal{T}},$$

$$\kappa \left(\frac{\partial \mathcal{T}}{\partial x_1}\frac{1}{\mathcal{T}}\right)={\mathcal{W}}^{\top }\left({\mathsf{C}}_{1,1}\frac{\partial \mathcal{W}}{\partial x_1}+{\mathsf{C}}_{1,2}\frac{\partial \mathcal{W}}{\partial x_2}+{\mathsf{C}}_{1,3}\frac{\partial \mathcal{W}}{\partial x_3}\right),$$

$$\begin{gathered}{\mathsf{D}}_{x_{1}}=\left(\mathsf{D}_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right),\quad\cdots\quad{\mathsf{D}}_{x_{3}}=\left({\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes\mathsf{D}_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right),\\ \\ {\mathsf{Q}}_{x_{1}}=\left(\mathsf{Q}_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right),\quad\cdots\quad{\mathsf{Q}}_{x_{3}}=\left({\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes\mathsf{Q}_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right),\\ \\ {\mathsf{B}}_{x_{1}}=\left(\mathsf{B}_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right),\quad\cdots\quad{\mathsf{B}}_{x_{3}}=\left({\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes\mathsf{B}_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right),\\ \\ {\Delta}_{x_{1}}=\left(\Delta_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right),\quad\cdots\quad{\Delta}_{x_{3}}=\left({\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes\Delta_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right){\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{,}}\\ \\ {\mathsf{P}}_{x_{1}}=\left(\mathsf{P}_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right),\quad\cdots\quad{\mathsf{P}}_{x_{3}}=\left({\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes\mathsf{P}_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right),\quad\\ \\ {\mathsf{P}}_{x_{1},{x_{2}}}=\left(\mathsf{P}_{N}\otimes\mathsf{P}_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right),\quad\cdots\quad{\mathsf{P}}_{x_{2},{x_{3}}}=\left({\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}}\otimes\mathsf{P}_{N}\otimes\mathsf{P}_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right),\\ \\ \mathsf{P}={\mathsf{P}}_{x_{1},{x_{2}},{x_{3}}}=\left(\mathsf{P}_{N}\otimes\mathsf{P}_{N}\otimes\mathsf{P}_{N}\otimes{\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{5}}}\right),\\ \\ \end{gathered}$$

$${\mathsf{B}}_{x_1}={\mathsf{B}}_{x_1}^{+}-{\mathsf{B}}_{x_1}^{-},{\mathsf{B}}_{x_1}^{-}=\mathrm{diag}\left(1,0,\dots ,0\right)\otimes {\mathsf{I}}_N\otimes {\mathsf{I}}_N\otimes {\mathsf{I}}_5,{\mathsf{B}}_{x_1}^{+}=\mathrm{diag}\left(0,\dots ,0,1\right)\otimes {\mathsf{I}}_N\otimes {\mathsf{I}}_N\otimes {\mathsf{I}}_5.$$

Full text

Using the method of lines, a semi-discrete entropy estimate for the entire domain is obtained when the proposed numerical imposition of boundary conditions are coupled with an entropy-conservative or entropy-stable discrete interior operator. The resulting estimate mimics the global entropy estimate obtained at the continuous level. The boundary data at the wall are weakly imposed using a penalty flux approach and a simultaneous-approximation-term technique for both the conservative variables and the gradient of the entropy variables.

Discontinuous spectral collocation operators (mass lumped nodal discontinuous Galerkin operators), on high-order unstructured grids, are used for the purpose of demonstrating the robustness and efficacy of the new procedure for weakly enforcing boundary conditions. Numerical simulations confirm the non-linear stability of the proposed technique, with applications to three-dimensional subsonic and supersonic flows. The procedure described is compatible with any diagonal-norm summation-by-parts spatial operator, including finite element, finite difference, finite volume, discontinuous Galerkin, and flux reconstruction schemes.

1 Introduction

Next-generation numerical algorithms for use in large eddy simulations and direct numerical simulations of computational fluid dynamics will rely on efficient, high-order formulations, that are able to deliver better accuracy per degree of freedom than low-order methods, and that feature much smaller numerical errors both in terms of dispersion and dissipation [3, 4]. While these properties make high-order methods well suited for time-dependent simulations, these techniques are more prone to instability when compared to their lower-order counterparts. This is because numerical instabilities may occur if the flow contains discontinuities or under-resolved physical features. Various stabilization strategies (e.g. filtering [5], artificial viscosity, over-integration, and slope limiting [4] to cite a few) are commonly used to address these issues. However, such stabilization techniques possess several drawbacks since i) they reduce accuracy [4], ii) they usually require tuning parameters for each problem configuration, and iii) they do not guarantee that solvers designed to be high-order accurate in space will not crash.

A very promising and mathematically rigorous alternative consists in focusing on discrete operators that are non-linearly stable777We use the term rigorous because, as we will see in the next sections, these operators can mimic at the discrete level the results of the non-linear stability analysis at the continuous level. or, as in the case of the compressible Navier–Stokes equations, entropy stable. These operators simultaneously conserve mass, momentum, and total energy. In addition, they satisfy a discrete analogue to the conservation or dissipation of entropy which, with positivity assumptions on temperature and density, guarantees an $L^{2}$ bound on the conservative variables [6, 7]. We remark that the idea of enforcing entropy stability in numerical methods is old and commonly used for low-order operators, see e.g. [8, 9]. For extensions to high-order accurate operators see [10, 11, 12, 13, 14].

Until recently, fully discrete entropy stability was mostly established for implicit time stepping schemes. However, Ranocha and colleagues [15] developped and applied new explicit Runge–Kutta schemes (i.e., relaxation Runge–Kutta schemes) to entropy conservative or entropy dissipative semi-discretizations of any order for the compressible Euler and Navier–Stokes equations. The new time integration schemes can conserve or dissipate any solution properties with respect to any convex functional by the additional of a relaxation parameter that multiplies the Runge–Kutta update at each step. The general technique is not limited to the compressible Euler and Navier–Stokes equations setting but can be applied to many ordinary differential equations, and to both explicit and implicit Runge–Kutta methods.

However, issues remain on the path towards complete entropy stability for the compressible Navier–Stokes equations, e.g. shock capturing and bound-preserving limiter for high-order accurate discretizations. One major obstacle is the need for boundary conditions that preserve the entropy conservation or stability property of the interior operator. Practical experience indicates that numerical instabilities frequently originate at domain boundaries; the interaction of shocks with these physical boundaries is particularly challenging for high-order formulations. An important step towards entropy stable wall boundary conditions for the compressible Euler and Navier–Stokes equations appears in [2, 16, 17]. More specifically, non-linearly stable wall boundary conditions for the compressible Navier–Stokes equations are presented in [2]. Therein, it is shown that entropy stability requires two conditions to be satisfied: i) Euler no-penetration, and ii) a prescribed value for the product of temperature and the gradient of the temperature in the normal direction to the wall. An additional term providing a controllable numerical dissipation has to be introduced to impose a zero relative velocity at the wall, i.e. the no-slip condition. Therefore, the solid wall boundary conditions proposed in [2] are entropy stable, but not entropy conservative. Note that in [17] it is shown that demanding a bound on velocity gradients necessitates the use of the full no-slip conditions, i.e. the thermal and the relative velocity boundary conditions.

In this work we present a general procedure for the development of point-wise entropy conservative boundary conditions representing either an adiabatic solid wall or a wall with a prescribed heat entropy flow for the compressible Navier–Stokes equations, discretized by using diagonal-norm, summation-by-parts (SBP) and simultaneous-approximation-term (SAT) operators (i.e. SBP-SAT operators). Entropy conservation is obtained by penalizing, using a SAT penalty, both the entropy variables and their gradients in the normal direction to the wall, as in the local discontinuous Galerkin approach [18]. The overall algorithm closely follows the treatment of the discontinuous interior interfaces coupling presented in [1]; a single implementation, with different inputs, can be used for interface penalization and imposition of boundary conditions. A controllable amount of dissipation can be added to make the boundary conditions entropy stable. The new procedure can be immediately applied to a moving wall, as will be shown in the theoretical and numerical results sections.

The manuscript is organized as follows. A brief review concerning the derivation of continuous entropy inequalities and the entropy analysis of the viscous wall boundary conditions for the compressible Navier–Stokes equations is provided in Section 2. The weak, point-wise, imposition of entropy conservative and entropy stable boundary conditions is carried out in Section 3 for an adiabatic solid wall and for a wall with a prescribed heat entropy transfer. Section 5 presents numerical results which confirm the accuracy and stability of the proposed boundary conditions. Conclusions are drawn in Section 6. Finally, in A a Python script is provided that symbolically verifies all proofs for curvilinear grids, while in B a simple and dimension-agnostic implementation of the entropy stable solid wall boundary condition coded in FORTRAN is presented.

2 A brief review of entropy stability theory

In this Section, we review the continuous entropy theory for the compressible Navier–Stokes equations and the solid wall boundary conditions by closely following [19, 2].

2.1 The compressible Navier–Stokes equations

To keep the presentation simple but without loss of generality, we consider the three-dimensional compressible Navier–Stokes equations in Cartesian coordinates $\left(x_{1},x_{2},x_{3}\right)$ for an ideal gas in a bounded domain $\Omega$ with boundary $\Gamma$

[TABLE]

The vectors $\bm{\mathcal{Q}}$, $\bm{\mathcal{F}}_{x_{m}}^{(I)}$, and $\bm{\mathcal{F}}_{x_{m}}^{(V)}$ respectively denote the conserved variables, the inviscid ($I$) fluxes, and the viscous ($V$) fluxes. The boundary data, $\bm{\mathcal{G}}^{(B)}$, and the initial condition, $\bm{\mathcal{G}}^{(0)}$, are assumed to be in $L^{2}(\Omega)$, with the further assumption that $\bm{\mathcal{G}}^{(B)}$ will be set to coincide with linear well posed boundary conditions and such that entropy conservation/stability is achieved.

The vector of conserved variables is given as

[TABLE]

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

[TABLE]

where $\mathcal{P}$ is the pressure, $\mathcal{H}$ is the specific total enthalpy, and $\delta_{i,j}$ is the Kronecker delta. The viscous flux $\bm{\mathcal{F}}_{x_{m}}^{(V)}$ is given as

[TABLE]

where $\kappa=\kappa(T)$ is thermal conductivity, and the viscous stresses is given by

[TABLE]

where $\mu=\mu(T)$ is the dynamic viscosity.

The required constitutive relations are

[TABLE]

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

[TABLE]

with $\mathcal{T}_{\infty}$ and $\rho_{\infty}$ the reference temperature and density, respectively.

It is well known that the compressible Navier–Stokes equations given in (1) possess a convex extension that, when integrated over the physical domain $\Omega$, only depends on the boundary data. Such an extension yields the entropy function

[TABLE]

which is a useful tool for proving stability in the $L^{2}$ norm [6, 7]. We can then define the entropy variables $\bm{W}=\partial\mathcal{S}/\partial\bm{\mathcal{Q}}$, which for the compressible Navier–Stokes equations given in (1) are

[TABLE]

We remark that the convexity of $\mathcal{S}$ guarantees the invertibility of the mapping between conservative and entropy variables, provided that the temperature $\mathcal{T}$ and the density $\rho$ are positive. In what follows, we always assume that such positivity is preserved.

The vector of entropy variables simultaneously contracts all of the inviscid spatial fluxes $\bm{\mathcal{F}}_{x_{m}}^{(I)}$ as

[TABLE]

where the scalar $\mathcal{F}_{x_{m}}(\bm{\mathcal{Q}})$ denotes the entropy flux in the $m$-th direction. By letting $\bm{\mathcal{W}}$ take the role as of a new set of independent variables, i.e. $\bm{\mathcal{Q}}=\bm{\mathcal{Q}}(\bm{\mathcal{W}})$, the entropy variables (7) symmetrize the system (1) as [20]

[TABLE]

where the viscous fluxes $\bm{\mathcal{F}}_{x_{m}}^{(V)}$ have been recast in term of the entropy variables as

[TABLE]

For the definition of the symmetric and semi-definite $\mathsf{C}_{m,j}$ matrices see [21, 2].

Due to the symmetric nature of $\partial\bm{\mathcal{Q}}/\partial\bm{\mathcal{W}}$ and $\partial\bm{\mathcal{F}}_{x_{m}}^{(I)}/\partial\bm{\mathcal{W}}$, there exist scalar functions in entropy variables whose Jacobians represent the conservative variables $\bm{\mathcal{Q}}$ and the inviscid fluxes $\bm{\mathcal{F}}_{x_{m}}^{(I)}$ as

[TABLE]

$\mathcal{\varPhi}$ is the potential, whereas the $\mathcal{\varPsi}_{x_{m}}$ functions are the potential fluxes in the $x_{m}$ direction, with $\left(\mathcal{\varPhi},\mathcal{\varPsi}_{x_{m}}\right)$ the potential-potential flux pair [22]. A close relation between entropy and the potential-potential flux pair is summarized in the following Theorem [23], which is due to Godunov (see also [24]).

Theorem 2.1**.**

If a system of conservation laws can be symmetrized by introducing new variables $\bm{\mathcal{W}}$, and $\bm{\mathcal{Q}}$ is a convex function of $\mathcal{\varPhi}$, then an entropy function $\mathcal{S}=\mathcal{S}(\bm{\mathcal{Q}})$ is given by

[TABLE]

and the entropy fluxes $\mathcal{F}_{x_{m}}(\bm{\mathcal{Q}})$ satisfy

[TABLE]

By contracting the system of equations (1) with the entropy variables,

[TABLE]

and applying the relations given in (8), (9), and (10), we arrive at the differential form of the (scalar) entropy equation

[TABLE]

To obtain a global conservation statement for the entropy function $\mathcal{S}$, we then integrate equation (15) over the domain $\Omega$

[TABLE]

where $n_{x_{m}}$ is the $m$-th component of the outward facing unit normal and

[TABLE]

We remark that viscous dissipation always introduces a negative rate of change in entropy, since the $-DT$ term in (16) is negative semi-definite. An increase in entropy within the domain can only result from data that convects or diffuses through the boundaries $\Gamma$. For smooth flows, we finally note that the inequality sign in (16) becomes an equality.

2.2 No-slip wall boundary conditions

For simplicity, we let the domain of interest be $\Omega=[0,1]^{3}$ and we only consider entropy conservation (i.e., the equality relation in (16)). Thus, expanding the notation in equation (16) yields

[TABLE]

Note that the plus and minus signs within the integrand terms of (17) account for the direction of the outward facing normals on the six faces of the unit cube $\Omega$.

Furthermore, without loss of generality, we consider the case of a wall placed at $x_{1}=0$ such that the normal vector is $\bm{n}=(-1,0,0)^{\top}$, and we assume that all the other boundaries terms are entropy conservative, which allows us to neglect their contributions. Then, estimate (17) reduces to

[TABLE]

Within the context of linear analysis for a solid viscous wall, the wall behaves like a subsonic outflow [25], and four independent boundary conditions must be imposed to prove energy (linear) stability [26, 27, 28] (see also [2] and the references therein). The first three correspond to the no-slip boundary conditions $\mathcal{U}_{1}=\mathcal{U}_{2}=\mathcal{U}_{3}=0$ that impose a zero relative velocity with respect to the wall. The fourth condition can be either imposed on the gradient of the temperature normal to the wall $(\partial\mathcal{T}/\partial n)_{wall}$ (Neumann boundary condition, e.g. the adiabatic wall), or to the temperature at the wall $\mathcal{T}_{wall}$, (the Dirichlet or isothermal wall boundary condition), or a mixture of these two (the Robin boundary condition) [27, 28].

In the non-linear case, entropy conservation and entropy stability in the adiabatic solid wall case or a wall with a prescribed heat entropy flow are attained by means of the next two theorems. These theorems provide the conditions that result in a bound on the time rate of change of the entropy function in (18), and are point-wise valid [2]. The first theorem is a generalization of Theorem 3.1 presented in [2] to a moving wall.

Theorem 2.2**.**

The no-slip boundary conditions $\mathcal{U}_{1}=0$ and $\mathcal{U}_{m}=\mathcal{U}_{m}^{wall},\,m=2,3$ bound the inviscid contribution to the time derivative of the entropy in equation (18).

Proof.

Equation (13) provides the following relation for the entropy flux

[TABLE]

Substituting the no-slip conditions, into the definition of the inviscid flux, $\bm{\mathcal{F}}_{x_{1}}^{(I)}$, (equation (3)) and the condition $\mathcal{U}_{1}=0$ into the definition of {\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathcal{\varPsi}}_{x_{1}}}}, yields the desired result $\mathcal{F}_{x_{1}}=0$. ∎

Remark 2.1**.**

In a general setting, the no-slip boundary conditions read as $\bm{\mathcal{U}}=\bm{\mathcal{U}}^{wall}$ with $\bm{\mathcal{U}}^{wall}\leavevmode\nobreak\ \cdot\leavevmode\nobreak\ {\bm{n}}=0$, where $\bm{\mathcal{U}}$, $\bm{\mathcal{U}}^{wall}$ and ${\bm{n}}$ denote the velocity vector of the fluid, the velocity vector of the wall and the unit normal vector, respectively.

Theorem 2.3**.**

The boundary condition

[TABLE]

where $\partial\mathcal{T}/\partial n$ denotes the normal derivative of $\mathcal{T}$, bounds the viscous contribution to the time derivative of the entropy (18).

Proof.

See Theorem 3.2 in [2]. ∎

Remark 2.2**.**

The scalar value,

[TABLE]

accounts for the change in entropy due to the wall heat flux at $x_{1}=0$ [2] and is often denoted as heat entropy transfer or heat entropy flow [29].

3 Entropy conservative and entropy stable solid wall boundary conditions

for the semi-discrete system

To discretize in space, we partition the physical domain $\Omega$ into non-overlapping hexahedral elements and we semi-discretize the system (1) using a multi-dimensional SBP operator, constructed from a one-dimensional SBP operator by way of tensor products. The nodal distribution within each element is based on $N^{3}$ Legendre-Gauss-Lobatto (LGL) points [30, 2, 12, 31], where $N$ is the number of LGL point in one direction.

Here, we summarize the relevant SBP operators used to discretize (1), and to derive the new procedure to impose the solid wall boundary conditions.

[TABLE]

$\mathsf{D}_{N}$, $\mathsf{Q}_{N}$, $\mathsf{B}_{N}$, $\Delta_{N}$ and $\mathsf{P}_{N}$ are the one-dimensional SBP operators, and {\leavevmode\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathsf{I}_{N}}} is the identity matrix of dimension $N$. The matrices $\mathsf{P}_{(\cdot)}$ may be thought of as mass matrices in the context of the discontinuous Galerkin finite element method. Herein, the focus is exclusively on diagonal-norm SBP operators, based on fixed element-based polynomials of order $p$ ($p=N-1$). The matrices $\mathsf{D}_{(\cdot)}$ are used to approximate the first derivatives and are defined as $\mathsf{P}^{-1}_{(\cdot)}\mathsf{Q}_{(\cdot)}$. The nearly skew-symmetric matrices $\mathsf{Q}_{(\cdot)}$ are undivided differencing operators where all rows sum to zero, and the first and last column sum to $-1$ and $1$ respectively. The matrices $\mathsf{B}_{(\cdot)}$ pick off the interface terms in the respective directions. For the spectral element discretization considered in this paper, the $\mathsf{B}_{(\cdot)}$ matrices take on a particularly simple form; as an example, consider $\mathsf{B}_{x_{1}}$, which is given as

[TABLE]

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.

SBP operators can be recast in telescoping flux form [21]. For example,

[TABLE]

where $\bar{\bm{f}}_{x_{1}}^{(I)}=\mathsf{I}_{S\text{to}F}^{x_{1}}\bm{f}_{x_{1}}^{(I)}$ and the one-dimensional telescoping operator, $\Delta_{N}$, is defined as

[TABLE]

The operator $\mathsf{I}_{S\text{to}F}^{x_{1}}$ interpolates the flux at the solution nodes on to a set of flux nodes (see Figure 1).

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.

[TABLE]

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

Using an SBP operator and its equivalent telescoping form, the semi-discrete form of the three-dimensional compressible Navier–Stokes equations (1) for Cartesian grids in each hexahedral element reads (see, e.g. [30, 2])

[TABLE]

We remark that we have omitted any term involved with the geometrical mapping of the elements in order to simplify the discussion. The extension of the following analysis to curvilinear grids is straightforward and all the theorems presented herein remain valid. The interested reader is referred to the supplementary python script (see A) which symbolically verifies all the proofs for the general case of curvilinear grids. Furthermore, a simple, dimension-agnostic implementation of the entropy stable solid wall boundary conditions for the general case is provided in B.

The vectors $\bm{{g}}_{x_{m}}^{(B)}$ in (22) enforce the boundary conditions, while $\bm{{g}}_{x_{m}}^{(In)}$ patches interfaces together using a SAT approach [1]. The derivatives appearing in the viscous fluxes $\bm{f}_{x_{m}}^{(V)}$ are also computed using the operators $\mathsf{D}_{x_{m}}$ defined in (21).

The discrete no-slip wall boundary conditions constructed herein follows a local discontinuous Galerkin-type approach [18], where also the gradient of the entropy variables is penalized. The penalization of the gradient represents one of the main novelties in this work within the context of SBP-SAT discretizations for the imposition of boundary conditions. In fact, this is a key difference with respect to previous work [2] which introduced for the first time an entropy stable approach to impose solid wall boundary conditions.

Following the same procedure based on local discontinuous Galerkin (LDG) and interior penalty approach (IP) described in [32, 1], equation (22) can be recast as

[TABLE]

where $\bm{\Theta}_{x_{m}}$ is the gradient of the entropy variables in the $m$-th direction, whereas $\bm{{g}}_{x_{m}}^{(B),q}$, $\bm{{g}}_{x_{m}}^{(B),\Theta}$ and $\bm{{g}}_{x_{m}}^{(In),q}$, $\bm{{g}}_{x_{m}}^{(In),\Theta}$ are the SAT penalty boundary ($B$) and interface ($I$) terms on the conservative variable $\bm{q}$, and the gradient of the entropy variable $\bm{\Theta}$, respectively [32]. The contributions of the interface penalty terms are non-zero only in the normal direction to the interface. The matrices $\left[\mathsf{C}_{m,j}\right]$ are block diagonal matrices with blocks of size 5, corresponding to the viscous coefficients at each solution point.

Remark 3.1**.**

In order to build a high-order accurate entropy-conservative/stable spatial discretization, the linear interpolation operation, $\mathsf{I}_{S\text{to}F}^{x_{m}}$ is replaced with a non-linear interpolation operator (or equivalently, the linear SBP operator is replaced with a non-linear SBP operator) [22, 21, 30, 2]. Critically, the resulting non-linear operator, ${\mathsf{P}}_{x_{m}}^{-1}\Delta_{x_{m}}\bar{\bm{f}}_{x_{m},sc}^{(I)}$ or ${\mathsf{P}}_{x_{m}}^{-1}\Delta_{x_{m}}\bar{\bm{f}}_{x_{m},ssr}^{(I)}$ (entropy conservative or entropy stable, respectively), has the property that when contracted with the entropy variables and discretely integrated over the domain, the result is a discrete surface integral with respect to the entropy flux; that is, for the entropy-conservative formulation, the operator telescopes in the entropy flux. For example

[TABLE]

*where $\bar{\bm{f}}_{x_{1},sc}^{(I)}$ is the vector that results from the entropy-conservative non-linear interpolation of the flux (this flux is replaced with $\bar{\bm{f}}_{x_{1},ssr}^{(I)}$ for the entropy-stable version), $\bm{1}$ is a vector of ones of appropriate size, and $\bm{F}_{x_{1}}$ is a vector of the entropy flux in the $x_{1}$ direction evaluated on the mesh nodes. *

To obtain an equation for the entropy of the system, we follow the entropy stability analysis presented in [30, 32]. Therefore, multiplying the two discrete equations by $\bm{w}^{\top}\mathsf{P}$ and $\left(\left[\mathsf{C}_{m,j}\right]\bm{\Theta}_{x_{j}}\right)^{\top}\mathsf{P}$, respectively, the expression for the time derivative of the entropy function in each element is

[TABLE]

where

[TABLE]

is a positive quadratic term in the approximation of the first derivative of the solution [2], $\mathbf{1}$ is the unit vector of appropriate size, and $\bm{F}_{x_{m}}$ the vector of the entropy flux in the $m$-th direction.

Equation (24) can be conveniently rewritten as

[TABLE]

This compact notation will be useful in Section 5.

As in the continuous analysis, we assume that we have a hexahedral element with edge length equal to one, and consider only the plane $(0,x_{2},x_{3})$ as wall boundary face. Thus, we have $\bm{{g}}_{x_{m}}^{(In),q}=\bm{{g}}_{x_{m}}^{(In),\Theta}=0,\leavevmode\nobreak\ \forall m$. We also assume that all the points that lie on the other faces of the cube are treated in an entropy stable fashion such that their contribution can be neglected. With these assumptions, equation (24) reduces to

[TABLE]

The penalty source term $\bm{{g}}_{x_{1}}^{(B),q}$ is composed of three design-order terms plus a source boundary term

[TABLE]

where

[TABLE]

$\bm{f}^{sc}\left(\bm{v},\bm{v}^{(B,I)}\right)$ is the entropy conservative flux vector, which is a function of the two arguments $\bm{v}$ and $\bm{v}^{(B,I)}$, and

[TABLE]

For clarity of presentation, the expressions of $\bm{\mathcal{M}}^{(B,V)}$, which is used to add provably dissipation, and of $\bm{\mathcal{L}}^{(B,V)}$ will be given and analyzed at the end of the Section.

The penalty $\bm{{g}}_{x_{1}}^{(B),\Theta}$ contains a single design-order term

[TABLE]

In each of the contributions, the first component (the numerical state) is constructed from the numerical solution, while the second component (the boundary state) is constructed from a combination of the numerical solution and four independent components of physical boundary data.

In what follows, we analyze each of these contributions from (26) separately. We further restrict our analysis to a single solution point lying on the wall boundary face. Therefore, we will no longer make use of the bold notation, which will be replaced by italics to denote the vector of five components at the collocated point.

The term given in (28) enforces the Euler no-penetration wall condition through the inviscid flux of the compressible Euler equations. The boundary state is formed by constructing an entropy conservative flux based on the numerical state in primitive variables at the face point, $v$, and a manufactured boundary state given by the vector of the primitive variables

[TABLE]

The term defined in (29), together with the penalty (30), allow a weak imposition of the no-slip condition

[TABLE]

in an entropy conservative way, where $\mathcal{U}^{wall}_{m}$ is the $m$-th component of the wall velocity. Furthermore, entropy flow is enforced if $\kappa\left(\frac{\partial\mathcal{T}}{\partial x_{1}}\frac{1}{\mathcal{T}}\right)\neq 0$.

Using the following primitive variables

[TABLE]

we define the point-wise boundary viscous flux ${\mathcal{F}}_{x_{1}}^{(B,V)}$, the term ${\mathcal{M}}^{(B,V)}$, and source term ${\mathcal{L}}^{(B,V)}$ as

[TABLE]

[TABLE]

[TABLE]

where $\mathtt{g}(t)=\kappa\frac{\partial\mathcal{T}}{\partial x_{1}}\frac{1}{\mathcal{T}}$ is a given $L^{2}$ function. Thus, in case of a boundary condition with imposed non-zero heat flux, $\mathtt{g}(t)$ is non-zero. We remark that the vector $v^{(B,V)}$ is used to evaluate the matrix of the viscous coefficients $\mathsf{C}_{1,j}^{(B,V)}$, as well as to compute the penalty term in (30).

The manufactured gradient of the entropy variables at the boundary $\Theta_{x_{j}}^{(B,V)}$ is constructed using the following procedure:

Rotate the gradient of the entropy variables $\Theta_{x_{j}}$ to the gradient of the primitive variables, $v$,

[TABLE]

where $\frac{\partial{V}}{\partial{W}}$ is the Jacobian of the transformation in primitive variables with respect to the entropy variables evaluated at the face point.

• 2.

Construct the gradient of the primitive variables at the wall boundary point

[TABLE]

This choice follows the imposition of Neumann boundary conditions in the context of the nodal discontinuous Galerkin method [5].

• 3.

Rotate the gradient of the primitive variables $\Pi_{x_{j}}^{(B,V)}$ to the gradient of the entropy variables

[TABLE]

where \frac{\partial{W}}{\partial{V}}\big{\rvert}_{(B)} is the Jacobian of the entropy variables with respect to the primitive variables evaluated using the state $v^{(B,V)}$ defined in (32).

Finally, the entropy variables $w^{(B,V)}$ needed in (30) are computed from the primitive variables defined in expression (32) by using the relations in (7).

The matrix $L$ in (34) is a negative semi-definite 5x5 matrix which is defined as

[TABLE]

where $\mathsf{C}_{1,1}$ and $\mathsf{C}_{1,1}^{(B,V)}$ are the positive semi-definite viscous coefficient matrices in the normal direction ($m=j=1$), respectively evaluated using the states $v$ and $v^{(B,V)}$, and $\beta$ is a positive coefficient that modulates the strength of the penalty term. This coefficient has to scale as the inverse of the typical element length888This is done such that the scaling remains dimensionally consistent..

Summarizing, the penalty at the face point for the conservative variables $\bm{q}$ is the sum of two terms:

the difference between inviscid and entropy conservative fluxes in the normal direction,

• 2.

the difference between internal viscous and boundary viscous fluxes in the normal direction.

The penalty on the gradient of the entropy variables $\bm{\Theta}$ is instead given by the difference between the solution at the node and the data imposed at the boundary expressed in terms of entropy variables.

The entropy conservation and stability of the penalty source terms (27) and (30) is demonstrated in the following three theorems. The first theorem, which ensures entropy conservation for the inviscid SAT penalty in (27), and which enforces the no penetration condition, is Theorem 5.1 in [2]. The second, is a new theorem and it ensures entropy conservation or stability for the viscous SAT penalty. The third theorem is also new and ensures that the term $\bm{\mathcal{M}}^{(B,V)}$ is a dissipative entropy contribution.

Theorem 3.1**.**

The penalty inviscid flux contribution in equation (27) is entropy conservative if the vector $v^{(B,I)}$ is defined as in (31).

Proof.

See Theorem 5.1 in [2]. ∎

Theorem 3.2**.**

The penalty terms for the viscous flux on the conserved variables (27), together with the viscous penalty on the gradient of the entropy variables (30), are

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

• 2.

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

Proof.

By substituting into (26) the expressions for $\bm{{g}}_{x_{1}}^{(B),q}$ with $\bm{\mathcal{M}}^{(B,V)}=0$ (i.e. no dissipation) and $\bm{{g}}_{x_{1}}^{(B),\Theta}$ given in (27) and (30), respectively, yields

[TABLE]

For an adiabatic wall $\mathtt{g}(t)=0$ and therefore the proposed boundary conditions are entropy conservative. For $\mathtt{g}(t)\neq 0$ the boundary conditions are entropy stable because the contribution to the time rate of change of the entropy function is only a function of the data $\mathtt{g}(t)$. ∎

Theorem 3.3**.**

The interior penalty term

[TABLE]

added to the SAT (27) is entropy dissipative.

Proof.

By expanding the contraction $\bm{w}^{\top}{\mathsf{P}}_{x_{2},{x_{3}}}\bm{{g}}_{x_{1}}^{(B),q}$ in (26), and by focusing on the dissipation term only, we arrive at the following point-wise contribution to the time-rate of change of the entropy function

[TABLE]

Note that we have omitted an extra positive scaling factor corresponding to the entry of the matrix ${\mathsf{P}}_{x_{2},{x_{3}}}$ associated with a wall boundary point. Plugging in the definitions of the matrix $L$ (39), $w$ and $w^{(B,V)}$ we obtain

[TABLE]

which completes the proof. ∎

Remark 3.2**.**

In order to construct an entropy stable solid wall boundary condition for the Euler equations, it is necessary to add provable entropy dissipation. One solution consists of replacing the entropy conservative flux $\bm{f}^{sc}$ in (27) with an entropy stable flux $\bm{f}^{ssr}$ in (27) as described in [2].

4 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 allows for a SAT implementation which is identical to the interface treatment shown in [1]. We can use a single subroutine with different inputs corresponding to the imposition of the interior interface couplings, or of the adiabatic solid wall or of the wall with a prescribed heat entropy flow. In fact, the interior interface coupling can be written as (see equations (16a-16d) in [1])

[TABLE]

which have exactly the same structure as LDG-IP approach used for the imposition of the solid wall boundary conditions except for the boundary penalty interface terms, $\bm{{g}}_{x_{m},r}^{(B),\cdot}$ in equation (23), which are replaced by the interior penalty interface coupling terms, $\bm{{g}}_{x_{m},r}^{(In),\cdot}$ in equations (42).

5 Numerical results

In this section we present four three-dimensional test cases which demonstrate the robustness of the new wall boundary conditions coupled with the family of high-order accurate entropy-stable interior SBP-SAT algorithms developed in [30, 2, 1, 12]. The systems of ordinary differential equations arising from the spatial discretizations are integrated using the fourth-order accurate Dormand–Prince method [33] endowed with an adaptive time stepping technique based on digital signal processing [34, 35]. We note that small enough tolerances are always used to make the temporal error negligible.

The unstructured grid solver used herein has been developed at the Extreme Computing Research Center (ECRC) at KAUST on top of the Portable and Extensible Toolkit for Scientific computing (PETSc) [36], its mesh topology abstraction (DMPLEX) [37] and scalable ordinary differential equation (ODE)/differential algebraic equations (DAE) solver library [38], and the Message Passing Interface (MPI). Additionally, the numerical solver is based on the algorithms proposed in [30, 2, 1, 12]. It uses a transformation from computational to physical space that satisfies both the entropy conservation and the geometric conservation law at the semi-discrete level [21]. Unless otherwise stated, the meshes used in this work have been generated using the GMSH package [39].

We present the numerical results for five test cases:

Laminar flow in a pipe with annular section to verify the accuracy of procedure for the imposition of the solid wall boundary conditions;

• 2.

Laminar flow in a lid-driven cavity to validate the entropy conservation and entropy stability properties of the interior discretization operator coupled with the boundary procedure;

• 3.

Laminar flow past a three-dimensional cylinder and a sphere to demonstrate the engineering capabilities of the interior discretization operator coupled with the boundary procedure;

• 4.

Supersonic turbulent flow past a three-dimensional rod with square section to demonstrate the robustness of the solver and solid wall boundary conditions (“standard” SBP-SAT operators and procedure for imposing solid wall boundary conditions based on linear analysis crash).

5.1 Flow in pipe with annular cross-section

In this section we investigate the accuracy of the solid wall boundary conditions. The proposed entropy stable no-slip wall boundary conditions do not force the numerical solution to exactly fulfill the boundary conditions. Instead the effect can be described as a rubber-band pulling the solution towards the boundary conditions. The computed boundary value (or numerical state) typically deviates slightly from the prescribed value but the deviation is reduced as the grid is refined. To verify the accuracy of the new procedure, we perform a grid convergence study for the flow in a full three-dimensional pipe with annular cross-section. For an incompressible flow, this test case has an analytical solution [40] for both the velocity distribution and the volume flux through the annular pipe. The expression for the axial velocity distribution, $\mathcal{U}_{1}$, as a function of the radial coordinate, $r$, is

[TABLE]

where $G$ is the pressure gradient forcing term and $R_{i}$ and $R_{o}$ are the inner and outer radii of the pipe, respectively. Herein, we set $R_{o}=0.5$, $R_{o}/R_{i}=4$ and $G/\mu=1$. We highlight that we have chosen this test problem because i) it has an analytical solution which cannot be represented exactly by the polynomial space of the numerical solution, and ii) it requires the use of curved boundary element faces to capture accurately the geometry of the annular section of the pipe.

The code that is used is a compressible code and in order to obtain results that are very close to those found for the incompressible equations, a Mach number of $M=1.0e-05$ is considered. Periodic boundary condition are used in the axial direction.

We run a grid convergence study for $p=2,3,4$ with a sequence of nested grids generated using rational Bezier basis functions such that the geometrical description of the pipe is preserved exactly for each refinement level.

The error in the axial velocity profile are computed using discrete norms as follows:

[TABLE]

where $\mathsf{J}_{\kappa}$ 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 grid convergence study are shown in Tables 1, 2 and 3 where the numbering the first column indicates the number of elements in the radial, angular and axial coordinates. It can be seen that the computed order of accuracy is very close to the formal value of $\sim(p+1)$.

5.2 Lid-driven cavity

Next, we validate the algorithm on the simple problem of the three-dimensional lid-driven cavity with adiabatic solid walls. The domain is a cube with sides of length $l$ discretized using a Cartesian grid composed of eight elements in each direction. A velocity field is imposed on one of the walls, corresponding to a rigid body rotation about the center of the wall at a speed $\omega$ (see figure 2). Based on the rotation velocity and the length of the cavity, this example is characterized by a Reynolds number $Re=l^{2}\omega/\nu=100$ and a Mach number $M=l\omega/c=0.05$. All the dissipation terms used for the interface coupling [1] and the imposition of the boundary conditions are turned off, including upwind and interior-penalty SAT terms. The two terms on the left-hand side of equation (25), $\frac{d}{dt}\mathbf{1}^{\top}\widehat{\mathsf{P}}\,\bm{\mathcal{S}}+\mathbf{DT}$, are monitored at every time step. Note that because we do not include any dissipation terms for the interface couplings and adiabatic wall boundary conditions are used, $\mathbf{\Xi}$ in equation (25) is zero. The monitored values are reported in Figure 3, left panel, together with the error committed in entropy conservation, right panel, which is below machine (double) precision. Therefore, we have confirmed numerically that the newly developed wall boundary conditions together with the interior discretization operator are entropy conservative when all the entropy dissipative terms are turned off.

We also present the results of a modified version of the lid-driven cavity example by considering non-adiabatic boundaries, and by imposing a nonzero entropy flux in a face adjacent to the rotating face, given as $\mathtt{g}(t)=10^{-4}\sin(4\pi t)$ (see figure 4). Figure 5 shows the terms on the entropy balance for this example, including the boundary contribution (left panel), and the error in the entropy conservation (right), which is below machine (double) precision.

5.3 Subsonic flow past a cylinder

We further explore the engineering capabilities of the entropy-stable SBP-SAT operators and the numerical procedure for the weak imposition of the solid wall boundary conditions by simulating the flow around a cylinder, a canonical example of external flows with important applications such as particle transport, fluid-structure interaction and bluff body aerodynamics, that has been extensively studied both numerically [41, 42, 43, 44] and experimentally [45, 46, 47, 48, 49].

We described the flow in a Cartesian coordinate system ($x_{1}$,$x_{2}$,$x_{3}$), with the free-stream velocity aligned in the $x_{1}$ direction. A circle of diameter $d$ is centered at the origin, with the domain of interest delimited by a rectangular box that respectively extends $20d$ and $60d$ upstream and downstream of the flow direction, and $30d$ in the $x_{2}$ direction. Such a 2D domain is then extruded a distance $d$ in the $x_{3}$ direction. We prescribe an adiabatic no-slip wall boundary condition on the surface of the cylinder; periodic boundaries are applied in the $x_{3}$ direction. The remaining faces of the box are treated as a far field.

The domain is discretized as follows: we first mesh the $(x_{1},x_{2})$ plane with second-order quadrilateral elements, and include a boundary layer around the cylinder wall as pictured in Figure 6. We then extrude this mesh using three layers of elements in the $x_{3}$ direction over the span of the cylinder.

Considering the free-stream velocity $u_{\infty}$ and the diameter of the cylinder $d$, the free-stream flow is characterized by a Mach number $M_{\infty}=0.1$ and a Reynolds number $Re_{\infty}=300$. Under these conditions, the flow developed behind the cylinder is three-dimensional and results in vortex shedding at a constant frequency.

In Table 4 we report the time average drag coefficient, $\bar{c}_{d}$, and the Strouhal number, $St=fd/u_{\infty}$ based on the frequency of the vortex shedding, resulting from uniform $p$-refinements (with polynomial orders for the solution space ranging from $p=2$ to $p=5$) and three levels of uniform $h$-refinement of the initial mesh described below, resulting in 714 (denoted by $rl=0$ in Table 4), 5,712 ($rl=1$), and 45,969 ($rl=2$) hexahedral elements respectively. Comparisons with results reported in the literature for these aerodynamic coefficients [50] are also provided. From these tables, it can be seen that in all cases the accuracy of the results improve by increasing the order of accuracy of the scheme and the grid resolution. The fifth- ($p=4$) and sixth-order ($p=5$) accurate entropy-stable schemes perform very well on the last level of $h-$refinement considered ($rl=2$), which is coarser compared to the typical grids used with second-order finite volume and finite differences schemes.

5.4 Subsonic flow past a sphere

We then test our implementation within a more complicated setting represented by the flow around a sphere. In this case, a sphere of diameter $d$ is centered at the origin, and a box is respectively extended $20d$ and $60d$ upstream and downstream of the flow direction; the box size is $30d$ in both the $x_{2}$ and $x_{3}$ directions. As boundary conditions, we consider adiabatic solid walls at the surface of the sphere and far field on all faces of the box.

The surface of the sphere is first triangulated using second-order simplices and a boundary layer composed of triangular prisms is extruded from the sphere surface for a total length of $3d$. The rest of the domain is meshed with an unstructured tetrahedral mesh. We then obtain an unstructured conforming hexahedral mesh by uniformly splitting each tetrahedron in to four hexahedral elements, and each prism in to three hexahedral elements, resulting in a total of 4,328 hexahedral elements. A cut of the final mesh is illustrated in Figure 7, together with a representative splitting of a tetrahedral and a prismatic cell.

The free-stream flow is characterized by a Mach number $M_{\infty}=0.1$ and a Reynolds number $Re_{\infty}=300$. Under these conditions, the flow developed behind the sphere is non-axisymmetric, with hairpin vortices shedding from the wake at a constant rate [51], inducing a total non-zero lift force on the sphere. Figure 8 depicts these hairpin vortices with isocontours of the Q-criterion colored by the vorticity magnitude at a given time instant.

As in the previous test case, we analyze the convergence of some quantities of interest under $p$- and $h$-refinement. The number of elements in the sequence of nested grids considered are 4,328, 34,624 and 276,992 respectively; the polynomial order of the solution ranges from $p=2$ to $p=5$. For this test, we monitor the time average drag coefficient of the sphere, $\bar{c}_{d}$, the average lift coefficient, $\bar{c}_{l}$ (with its orientation $\theta$ in the $x_{3}-x_{2}$ plane), and the Strouhal number $St$ based on the frequency of the vortex shedding. We remark that $\bar{c}_{l}$ is nonzero because of the asymmetric flow features [51]. The results provided in Table 5 are in good agreement with those reported in the literature [52] for sufficiently refined meshes (i.e. $rl=1$ and $rl=2$), and polynomial orders greater than or equal to 3.

Figure 8 shows the 0.02 isocontour of the Q-criterion at $t=150$ for the solution of the flow around a sphere with $rl=2$, and $p=5$. At this time, the flow has reached a constant periodic state. From this figure, it can be seen the asymmetric pattern of the flow, consistent with a nonzero $\bar{c}_{l}$. Notice the heads of the hairpin vortices are almost aligned with the $x_{2}$ axis, in agreement with our reported value of $\theta$ for this solution.

5.5 Supersonic flow past a square cylinder

We finally provide further evidence of the robustness of the algorithm in the context of supersonic flow around a square cylinder with $Re_{\infty}=10^{4}$ and $M_{\infty}=1.5$, which features shocks, expansion regions and three-dimensional vortical structures [2]. We start with a square of side $s$ placed in the $x_{1}-x_{2}$ plane, and construct an unstructured mesh around it, manually refined in order to capture the main features of the flow (see Figure 9). We then extrude the mesh for a total size $s$ in the $x_{3}$ direction using four elements; the final three-dimensional mesh used in the study consists of 87,872 hexahedral elements. The boundary conditions imposed are adiabatic solid wall on the square cylinder surfaces, periodic boundary conditions in the $x_{3}$ direction, and far field at the remaining boundaries. The problem is solved using a fourth-order accurate ($p=3$) discretization.

Figure 10 show the results for the supersonic square cylinder at $t=100$. At this point in time, the flow is fully unsteady and the shock in front of the cylinder has reached its final position. The flow is characterized by the shock in front of the square cylinder and those in the near wake region. There is also an unsteady wake populated by three-dimensional vortices shedding from the body.

We finally remark that the small oscillations near the shock region are caused by discontinuities in the solution and are expected for this scheme. In fact, we are not using any shock capturing method or reducing the order of scheme at the discontinuity. Nevertheless, the simulation remains stable at all time, and the oscillations are always confined to small regions near the discontinuities. This is a feat unattainable with several alternative approaches to wall boundary conditions based on linear analysis which for this test problem lead to numerical instabilities and an almost immediate crash of the solver.

6 Conclusions

We have used entropy stability and the summation-by-parts framework to derive entropy stable wall boundary conditions for the three-dimensional compressible Navier–Stokes equations in the presence of an adiabatic wall, or a wall with a prescribed heat entropy flow.

A point-wise entropy-stable numerical procedure has been presented for weakly enforcing these solid wall boundary conditions at the semi-discrete level combining a penalty flux approach with a simultaneous-approximation-term technique for the conservative variables and the variables representing the gradient of entropy. The resulting semi-discrete operator mimics exactly the behavior at the continuous level, and the proposed non-linear boundary treatment provides a mechanism for ensuring non-linear stability in the $L^{2}$ norm of the continuous and semi-discretized compressible Navier–Stokes equations.

The design order properties of the scheme are validated in the context of laminar flow in a pipe with an annular section. Detailed viscous numerical computations in a three-dimensional subsonic lid-driven cavity flow have been presented to assess the accuracy of the proposed numerical techniques. The error in the entropy function balance showed an excellent agreement with the theory with or without heat entropy flux.

Unsteady laminar flow past a cylinder and a sphere have been presented to highlight the efficacy in computing aerodynamic forces; numerical simulations considering both $p$- and $h$-refinements showed very good agreement with results available from the literature.

The robustness of the complete semi-discrete operator (i.e. the entropy-stable interior operator coupled with the new boundary treatment) was demonstrated for the supersonic flow past a three-dimensional square cylinder at $Re_{\infty}=10^{4}$ and $M_{\infty}=1.5$, as proposed in [2]. This test has been successfully computed with a fourth-order accurate method without the need of introducing artificial dissipation, limiting techniques, or filtering, for the purpose of stabilizing the computations, a feat unattainable with several alternative approaches based on linear analysis only.

Although the robustness and efficacy of the techniques presented in this work have been validated using discontinuous spectral collocation operators on unstructured grids, the new boundary conditions can be applied to a very broad class of spatial discretizations and they are compatible with any diagonal-norm summation-by-parts spatial operator, including finite element, finite difference, finite volume, discontinuous Galerkin, and flux reconstruction schemes.

Acknowledgments

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.