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

Nonlinear entropy stability analysis is used to derive entropy stable no-slip wall boundary conditions for the Eulerian model proposed by Sv\"{a}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.


Introduction
The classical compressible Navier-Stokes (CNS) equations can be derived based on the material (Lagrangian) derivative formulation [30].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 [33], 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 [12,32].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 [33] and their semidiscrete counterpart.As done in [26] 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-approximationterm (SAT) penalty technique.At the semidiscrete level, the work in [26] was sharpened in [13] 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 [26,13] 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 [33], Dolejší and Svärd [15], 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, fullydiscrete entropy stable SSDC solver described, validated, and verified in [27].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.

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.

General form of the Eulerian model
Svärd [33] arrives at the following form of the Eulerian model for viscous and heat conducting compressible flows: where i, j = 1 : N DIM (in MATLAB notation), and ρ, U i , p, and 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 where R is the gas constant.The specific total energy, E, is given by the following equation 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, ν, is given by where µ is the dynamic viscosity, α ∈ [1, 4  3 ], and β is an additional diffusion coefficient.Finally, the speed of sound, c, is defined as where γ = cp cv .Remark 1 Equation ( 4) is the most general expression of the kinematic viscosity of the Eulerian model.However, in practice, Svärd [33] chooses β = 0 and thus, equation ( 4) reduces to a scaled kinematic viscosity coefficient

Entropy analysis
We cover the entropy analysis of the Eulerian model ( 1) by rewriting it using the following compact form, where f i , and f are the inviscid and viscous fluxes in the x i direction, respectively.
The next theorem ensures that the Eulerian model is entropy stable.

Theorem 1
The following boundary integral bounds the time derivative of the entropy function, S, of the Eulerian model (1).
Proof Given an entropy pair (S, F i ) that symmetrizes the Eulerian model (see [33]), we define the entropy variable w = dS dq .Then, using the matrix dq dw , we change the variables of system (1) to q := q(w).These steps yield the symmetric form dq dw ∂w ∂t + ∂g i , and g 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 dq dw Further, as shown by Tadmor [36], the contribution of the symmetrized inviscid flux term yields w ∂g i.e., the divergence of the entropy flux.Lastly, the viscous term contribution can be manipulated to obtain the following expression: Therefore, using Equations ( 10), ( 11) and ( 12), the scalar partial differential equation for the entropy function, S, integrated over a generic domain Ω with boundary Γ reads Now, using the fundamental theorem of calculus on ∂Fi ∂xi , and the integrationby-parts rule on w where Because the kinematic viscosity, ν, is positive and the matrix dq dw is symmetric positive-definite, the product ν dq dw is also a positive-definite matrix.Thus, DT is positive and appropriately bounded (See Svärd [32]).We obtain the following expression 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.

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 [26,13].
The entropy-entropy flux pair (S, F i ) = (−ρs, −ρU i s), used in [13,26], is the only pair that admits a diffusive entropy flux for the CNS model [33].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 [17,9] where h = c pr T is the enthalpy.

Remark 2
The entropy-entropy flux pair for the Eulerian model is not necessarily restricted to (S, F i ) = (−ρs, −ρU i s).Svärd [33] 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 From the proof of Theorem 1, we arrive at Integrating the previous expression over a cubical domain, Ω = (0, 1)×(0, 1)× (0, 1), and, without loss of generality, considering a solid wall located at x i = 0, we get the following contributions [26] Ω where j, k = i.From Theorem 1, we get the bound 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 [33].Four comes from the linear analysis of the CNS model [3,34], in addition to a boundary condition that rises from density diffusion [33].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 [3,34,33].These five boundary conditions lead to linear stability of both the CNS and Eulerian model [3,34,33].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.

Inviscid contribution
The inviscid contribution to the time rate of change of the entropy function, appearing on the RHS of (21), is treated as in [13, Theorem 2.2].However, for clarity of presentation and completeness, we report it here.
Theorem 2 The no-slip boundary conditions U i = 0 and U j = U wall j , where j = i, bound the inviscid contribution of the time derivative of the entropy function, i.e.
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 [13, Theorem 2.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., The main result is provided in the following theorem.
Theorem 3 The boundary conditions 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 that contributes positively to the time derivative of the entropy function.Then, setting 1 ρ ∂ρ ∂xi = 0 yields the boundary condition 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 κ 1 T ∂T ∂xi [26].However, for the Eulerian model, we require the impostion of two boundary conditions, i.e. g(t) = µ 1 T ∂T ∂xi and 1 ρ ∂ρ ∂xi = 0.In practice with the above boundary conditions, the time derivative of the entropy function satisfies the following relation Bound ( 28) can be translated into a bound for the conserved quantities [12,32], and hence, for the primitive variables.

Semidiscrete entropy stable framework
Herein, using summation-by-parts (SBP) operators [35] and the simultaneousapproximation-technique (SAT) [8,23], we provide an entropy stable framework of any order for the semidiscretization of the Eulerian model (1) using unstructured grids.

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, D ∈ R N ×N , is a linear SBP operator of degree p approximating the derivative ∂ ∂xi on the domain 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 onedimensional SBP operators used in this work are explicitly constructed in [10].Their extension to two-and three-dimensions is achieved using tensor product operations [10,26]: where D N , Q N , B N , ∆ N and P N are the one-dimensional SBP operators, and I N is the identity operator. 1In this context, we chose diagonal P N , and by the definition of SBP operators The matrices B (•) pick off the interface terms in the respective directions.For the spectral element discretization considered in this paper, the B (•) matrices take on a particularly simple form; as an example, consider B x1 , which is given as 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.P = (P N ⊗ P N ⊗ P N ) .
We finally note that in the present work, the quadrature nodes and solution nodes are collocated.

Semidiscretization of the Eulerian model
Using the operators shown in (29), we can write the semidiscretization of (1) as where vectors g xi enforce the boundary conditions, while g patches interfaces together using a SAT approach [26].The bolded letters represent quantities, q, and functions, f and g for all nodes in an element.
Following [17,9], we use the telescoping property of an SBP operator, to re-write the semidiscrete counterpart to the equation (30) as The vector f (I) and the one-dimensional telescoping operator, ∆ N , is given by The operator I StoF interpolates the value at the solution nodes to the interfaces between nodes as shown in Figure 1.Additionally, we define dq dw as a block diagonal matrix applied to all LGL points in an element.The kinematic viscosity ν is a scalar and we can bring it inside the matrix dq dw , i.e., ν dq dw .The block diagonal matrix ν dq dw commutes with the operator D xi .Thus, we can re-write the viscous flux as where the symbol Θ xi is the gradient of the entropy variables in the x i direction.
Using the procedure based on local discontinuous Galerkin (LDG) and interior penalty approach (IP) described in [26,13], the semidiscretization (32) can be recast as The terms g , and g , g (In),Θ are the SAT penalty boundary (b) and interface (In) terms on the conservative variables, q, and the gradient of the entropy variables, Θ, 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 I StoF is replaced with a nonlinear interpolation operator [17,9].Thus, the operator i,ss for the entropy stable case) is a nonlinear operator that is a discrete counterpart to the term that appears in Theorem 2. Thus, we arrive at the following relation (35) The flux vector f (I) i,sc is constructed using an entropy conservative two-point flux [36], and therefore, it satisfies the relation (11).In the case of an entropy stable flux, f (I) i,ss , the resulting term is entropy dissipative [9].

Time derivative of entropy function
Following the entropy analysis detailed in [9,26,25], the semidiscretization (34) yields the following expression for the time derivative of the entropy function, S: where and 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., g (In),q xi = g (In),Θ xi = 0. Therefore, equation (38) reduces to 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: 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 [26].
The first term enforces Euler no-penetration wall where f (b,I) i,sc is the entropy conservative flux and v represent the primitive variables.The manufactured inviscid states are defined as in [26] for i = 1, without loss of generality, The viscous boundary term is given as where f (b,V ) i will be defined later in this section.Together with (43), the following term enforce the no-slip boundary condition weakly.The analysis of the previous viscous terms match the analysis described by [26,13] since the entropy variables are the same, and the matrix ν dq dw 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 where the manufactured viscous boundary primitive variables, v (b,V ) , is defined, as in [13], by The change of variables matrix, dq dw (v), as a function of the primitive variables, v, is defined in Appendix B. The term Θ xi is the manufactured gradient of the entropy variables.The construction of Θ xi 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 ) , where β 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 [13].
We summarize the results for the RHS of (36) in the following three theorems.The first theorem is a result from [26] and is reprinted here for completeness Theorem 4 The penalty inviscid flux contribution, 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 [26,Theorem 5.1].

Theorem 5
The penalty terms for the viscous flux in the conserved varaibles, g (b,V ),q , and the gradient of entropy variables, g (b,V ),Θ , are entropy conservative if the wall is adiabatic, i.e. g(t) = 0, entropy stable in the presence of a heatflux, g(t) = 0, where g(t) is a given L 2 function.
Proof Substituting the expressions for g (b,V ),q , and, g (b,V ),Θ into (39), yields 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.
Proof By expanding the penalty terms with respect to the conservative variables, g (b,V ),q , and focusing only on the dissipative IP term, M (b,V ) , in (36), we arrive at Thus, 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.

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 [25].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 [25]) ∂q r ∂t 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, g xi,r in equation (34), which are replaced by the interior penalty interface coupling terms, g (In),• xi,r in equations (51).

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 [27].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) [2], its mesh topology abstraction (DMPLEX) [21], and its scalable ODE/DAE solver library [1].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 α = 1.We use the two-point entropy conservative flux presented by Chandrashekar [11].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 [4] with adaptive time-step and both relative and absolute tolerances set to 10 − 8.The meshes are generated using Gmsh [19], 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.Table 1: Convergence study for the flow in a pipe with annular cross-section at Re = 1 and Ma = 1e − 3. axial velocity error; solution polynomial degree: p = 2.
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 [13], this case is considered for two reasons.First, it has an analytic solution for incompressible flow that cannot be represented in polynomial space [28].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: where R o = 0.5, R o /R i = 4, and G/µ = 1.We consider uniform density and temperature, and zero nonaxial velocities.The parameters used are Re = 1, Ma = 1e−3 and α = 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 [20].
The error calculation uses the following discrete norms where Ω represents the volume of the computational domain, 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).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 âr = ar ar , where a r = (1, 1, 1  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 Re = 1, and 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 4a shows the time derivative of the entropy function of the semidiscretization of the Eulerian model ( 30), d dt 1 PS, and the dissipation term, DT, as a function of time, t (see Equation ( 36)).As shown in Figure 4b, 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.In this section, we present the two-dimensional flow around a cylinder of diameter D = 1 at Re = 40, and 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 steadystate.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

SSDC Eulerian SSDC CNS
Park et al. [24]   Fig. 8: Wall pressure coefficient for the flow past a cylinder at Re = 40 and Ma = 0.07.The reference data are obtained from [24].Vorticity SSDC Eulerian SSDC CNS Park et al. [24]   Fig. 9: Vorticity at the wal for the flow past a cylinder at Re = 40 and Ma = 0.07.The reference data was obtained from [24].
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 [7] and Russell and Wang [29].
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 7a.The difference observed is accumulated at the leading edge and remains consistent along with the tailwind, as shown in Figure 7b.
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. [24].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 • where they both undershoot the reference results.

Blast-wave
Svärd [33] 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 Re = 10 and Ma = 0.07, and the same initial conditions as in [33].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 [33] 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 [33] stable.As we will see shortly, we are unable to reproduce the "lagging" effect observed in [33].
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 [33], 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.

Supersonic flow around a cylinder
In this last test case, we present some results for the flow past a two-dimensional cylinder enclosed between two solid walls [6].The similarity parameters are  Fig.12: Supersonic flow past a cylinder enclosed between two plates at Re = 10, 000 and Ma = 3.5.Re = 10, 000 and Ma = 3.5.The initial condition is a uniform flow with unit density, ρ, unit temperature, T , and only a non-zero unit velocity component in the x 1 direction, 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 [16].The solution polynomial degree distribution is shown in Figure 11, and it is the same setup used in [27].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, ρ, and the velocity component in the x 1 direction, U 1 , of the developed solution.
A quantitative analysis of the numerical results can be performed using the oblique shock wave theory [31].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, θ, as a function of the shock wave angle, β, and local Mach number, Ma [22]: The pair θ − β 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 ∆ = 0.299 and ∆ = 0.308 for the Eulerian and CNS models, respectively.These values are in good agreement with the analytical value, which reads ∆ = 0.293 [5].

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 [33].In that context, using summationby-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.
Since P is diagonal, we have Thus, Equation (55a) is recast in the following form: B Change of variables matrix, 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 where

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, M = 0, the RHS of (36) becomes Define W = diag(w) and write where WPx j ,x k = Px j ,x k W, since W is diagonal, and χ is a vector on a node.We use (71) for both terms in the RHS of (70).The first term reads

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.

Fig. 3 : 1
Fig. 3: Spinning sphere enclosed in a cubical box; Eulerian model with Re = 1, and 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.

Fig. 4 :
Fig. 4: Verification of the entropy conservation for the spinning sphere enclosed in a cubical box; Eulerian model with Re = 1, and Ma = 0.05.

Fig. 6 :
Fig. 6: Wake region of the flow past a cylinder at Re = 40 and Ma = 0.07.The background is colored by the magnitude of the pointwise velocity vector.

Fig. 7 :
Fig. 7: The density profile for the flow past a cylinder at Re = 40 and Ma = 0.07.

Fig. 11 :
Fig. 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 [27].

Table 2 :
Convergence study for the flow in a pipe with annular cross-section at Re = 1 and Ma = 1e − 3. axial velocity error; solution polynomial degree: p = 3.

Table 3 :
Convergence study for the flow in a pipe with annular cross-section at Re = 1 and Ma = 1e − 3. axial velocity error; solution polynomial degree: p = 4.

Table 4 :
Bubble length, L bubble , and drag coefficient, C D .
. Our discretization of the Eulerian and CNS models are entropy 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.
dt 1 PS + 1 Px j ,x k Bx i Fx i + (Dx i w) P ν dq dw Θx j = w Px j ,x k Bx i ν dq dw Θx j + w Px j ,x k g Px j ,x k Bx i Fx i + w Px j ,x k Bx i ν dq dw Θx j + w Px j ,x k g Px j ,x k Bx i Fx i + w Px j ,x k Bx i ν dq dw Θx j + w Px j ,x k g