Impact of wall modeling on kinetic energy stability for the compressible Navier-Stokes equations
IImpact of wall modeling on kinetic energy stability forthe compressible Navier-Stokes equations
Vikram Singh a, ∗ , Steven Frankel a , Jan Nordstr¨om b,c a Faculty of Mechanical Engineering, Technion - Israel Institute of Technology, HaifaIsrael b Department of Mathematics, Computational Mathematics, Link¨oping University,SE-58183 Link¨oping, Sweden c Department of Mathematics and Applied Mathematics, University of Johannesburg,P.O. Box 524, Auckland Park 2006, South Africa
Abstract
Affordable, high order simulations of turbulent flows on unstructuredgrids for very high Reynolds’ number flows require wall models for efficiency.However, different wall models have different accuracy and stability prop-erties. Here, we develop a kinetic energy stability estimate to investigatestability of wall model boundary conditions. Using this norm, two wall mod-els are studied, a popular equilibrium stress wall model, which is found to beunstable and the dynamic slip wall model which is found to be stable. Theseresults are extended to the discrete case using the Summation-by-parts (SBP)property of the discontinuous Galerkin method. Numerical tests show thatwhile the equilibrium stress wall model is accurate but unstable, the dynamicslip wall model is inaccurate but stable.
Keywords:
Discontinuous Galerkin, Summation-by-parts, Wall modelling, Stability,Skew-symmetric form ∗ Corresponding Author, Present Address: Max Planck Institute for Meteorology, Ham-burg, Germany
Email addresses: [email protected] (Vikram Singh), [email protected] (Steven Frankel), [email protected] (Jan Nordstr¨om) a r X i v : . [ phy s i c s . f l u - dyn ] F e b . Introduction High-order methods have been shown to provide more accurate solutionsthan low-order one’s for the same number of degrees of freedom [14]. Here,we focus on the discontinuous Galerkin method (DG) [27, 2, 7], which hasbecome a popular tool for high order simulations for fluids over the pastfew decades. It can be made conservative, stable and arbitrarily high-orderaccurate on unstructured grids for linear problems.There are studies indicating that DG is a suitable method for implicitlarge eddy simulations (ILES) [3]. However, for wall resolved LES the totalnumber of grid points required are of the order of ∼ Re / L , where Re L is theReynolds’ number based on streamwise length [6]. The resources requiredto do simulations on such large grids are prohibitive. One way of avoidingthis requirement is to use wall modeled LES (WMLES), where the no-slipboundary condition is replaced by a suitable model. In the case of WMLES,the grid points requirement scales as ∼ Re L which is much more affordable.Many different approaches exist for wall modeling [19]. They can, in gen-eral, be categorized into algebraic and ordinary differential equation (ODE)based models. In the ODE based models, two grids are required. One gridis used to solve the ODE governing the wall model which provides the wallstress boundary conditions for the full set of equations. This procedure in-troduces an additional level of complexity. The algebraic models are simplerto implement requiring only solutions of algebraic equations. Recent studieshave shown the effectiveness of using such wall models for DG [10]. An-other more recent approach to wall modelling is the dynamic slip wall model[5, 1]. It relies on deriving boundary conditions for the filtered Navier-Stokesequations. This method requires a length scale as a regularization param-eter which ensures that the no-slip boundary condition is recovered as theresolution approaches direct numerical simulation (DNS) levels.Stability is highly desirable for any numerical method. Various meth-ods exist to stabilize the discrete nonlinear equations ; see [33, 23] and thereferences therein. However stability issues arising from the use of the wallmodel boundary conditions are not well studied. If the boundary conditionsof the initial boundary value problem (IBVP) are not well-posed, there isno guarantee that the simulation will be stable or even meaningful. Variousstudies have been done to determine boundary conditions that satisfy stabil-ity estimates [32, 4, 22]. However, these studies usually deal with standardboundary conditions such as no-slip, no-penetration or inflow/outflow type2f boundary conditions.In this paper, we consider two wall models. The first one is an equilibriumstress model and the second one a dynamic slip model. Next, we derivekinetic energy stability estimates for the wall boundary conditions. Usingthis bound, the equilibrium stress model is shown to be unstable while thedynamic slip wall model is stable. We then introduce the DG method andthe associated SBP property. The stability estimates for the wall models arethen extended to the discrete case using SBP. Finally numerical results arepresented that compare the accuracy and stability of the two wall models.
2. Wall models
There are various techniques for wall modelling in LES (WMLES) [19].Here we will explore an equilibrium wall stress model [10], and the morerecent dynamic slip wall model [5].
For the compressible Navier Stokes equations, four boundary conditionsare required at the wall [32]. For example, the no-slip boundary conditionis imposed by requiring the three components of the velocity to be zero andimposing either an isothermal or adiabatic condition for temperature. Inthis paper, the isothermal boundary condition will be considered. Note thatisothermal implies temperature T = C where C is the given boundary data.In the following sections on the wall model, the notation () ∗ , for example, u ∗ will be used for boundary data. For example, the no slip boundary conditionis u ∗ i = 0 , i = 1 , , T ∗ = C . The two wall models shown here use eithera Neumann or Robin boundary condition. The other boundary conditionsare the no-penetration slip and the isothermal wall. In a rectangular domain { ( x, y ) : 0 ≤ x ≤ , ≤ y ≤ , ≤ z ≤ } or [0 , , with y = 0 as the wallthe boundary conditions will be implemented as u ∗ i = 0 , i = 1 , T ∗ = C and a boundary value for τ ∗ where τ is the shear stress. The equilibrium wall stress model assumes that the flow is fully turbulentand in equilibrium. Under this assumption the wall stress is determined usinga log-law like profile [29] for the near wall velocity. To present the technique,we first define the so called friction velocity as u τ = (cid:112) τ w /ρ . Then the scaledwall normal co-ordinate is defined as y + = yu τ /ν and the scaled velocity as3 + = u/u τ . Using these scaled variables, the algebraic Reichardt function [28] is u + = 1 κ ln(1 + κy + ) + ( C − κ ln κ )(1 − e − y +11 − y + e − y +3 ) , (1)where we choose C = 4 . κ = 0 .
38. Given some location y wm away fromthe wall and a velocity u wm at this location, (1) defines an implicit scalarequation for τ w . Given the velocity u wm = [ u, v, w ] wm at y wm , we simply let u wm = | u wm | with || implying magnitude. The choice of the location y wm isimportant for the performance of this method. The convention is to choosea location sufficiently far away from the wall and yet at a fraction of thelocal boundary layer thickness y wm ∼ . δ . Here, following [10], we takethe solution point farthest away from the wall boundary point in the walladjoining element, see figure 1. Using this velocity, wall stress is determinedsuch that (1) holds. In the case of Gauss Lobatto solution points, the wallboundary points are also solution points. Figure 1: Illustration of wall modelled element. The solution points are in blue and black,while the red points are the nodes where wall boundary condition is imposed. The velocity u wm used for the wall model is the velocity at the black solution points. Given y wm , u wm , we solve for the wall stress using (1). It is an implicitequation in τ w which we solve numerically. The result has to be imposed asthe viscous flux, hence we need to convert the scalar value τ w to a vector.We partition the stress using the velocity vector. For example, in a domain40 , the expression for the shear stress at y = 0 wall becomes, τ ∗ xy = u wm,x τ w | u wm | , τ ∗ yz = u wm,z τ w | u wm | (2)where u wm,x , u wm,z are the tangential components of the velocity in the x, z directions respectively at y = y wm and τ ∗ xy , τ ∗ yz is the shear stress boundarydata. The process to implement the wall model is summarized below for eachwall adjoining element.1. Loop through boundary points2. For each boundary point, find the surface normal pointing into theflow. Using this normal find the solution point inside the element thatis furthest away from the boundary point along this normal. Take theheight of this point from the boundary point as y wm .3. Let the velocity magnitude at that point with height y wm be u wm .4. Use ρ, τ w at the boundary as initial guess for u τ = (cid:112) τ w /ρ and getinitial values for u + , y + .5. Get updated value for u + , y + by substituting previous values in (1).Iterate to convergence for u + .6. The value of u + thus obtained provides u τ and thus τ w . Use (2) to getthe shear stress boundary data and impose this as a Neumann boundarycondition using simultaneous approximation terms (SAT) [33].7. The other boundary conditions used are the no-penetration conditionfor normal velocity and a fixed temperature isothermal wall. Impositionof these boundary conditions are also done using SAT. The dynamic slip wall model was introduced in [5]. The boundary condi-tion was derived for the filtered Navier Stokes equations and reduces to theno-slip condition as the grid is refined. It is given by u i − l p ∂u i ∂y = 0 , i = 1 , , , (3)at the wall. Here u i indicates the velocity components. This is a Robin typeboundary condition that reduces to the no-slip boundary condition as l p goesto zero. We introduce two changes to the wall model (3). First, note that (3)allows penetration. In a rectangular domain { ( x, y ) : 0 ≤ x ≤ , ≤ y ≤ }
5r [0 , , (3) gives us u − l p ∂u ∂y = 0. Here u is the velocity in y direction. Wemodify the boundary condition to directly impose no-penetration, u = 0.The modified boundary condition is now mass conserving. Second, we donot use the dynamic procedure from [5] to obtain l p . We will treat it as aparameter with a constant positive value following [9]. We will study theeffects of this parameter on the accuracy of the wall model, in the numericalsection.Note that the Robin boundary condition (3) needs to be converted toboundary data such as (2). We do this by using penalization. To illustratethis, consider again the domain [0 , . At y = 0 wall, let the local values forthe velocity, shear stress etc. be u , τ xy . Using these values the boundarydata will be given by τ ∗ xy = τ xy + σ ( u − l p ∂u ∂y ) .τ ∗ yz = τ yz + σ ( u − l p ∂u ∂y ) . (4)Here τ ∗ xy , τ ∗ yz is the shear stress boundary data and τ xy , τ yz , u , u , ∂u ∂y , ∂u ∂y arelocal shear stress, wall parallel velocity in the x, z directions and gradient ofwall parallel velocity evaluated at the wall. σ is a parameter, which we willdetermine based on the kinetic energy stability, analysed later. The dynamicslip wall model implementation is summarized below.1. Loop through boundary points2. Choose a constant positive value for the parameter l p . Here we willtest three cases with l p = 0 . , . , . u , u , τ xy , τ yz , ∂u ∂y , ∂u ∂y .4. Using these values, determine the boundary data from (4) and imposeusing a SAT term.5. The other boundary conditions used are the no-penetration conditionfor normal velocity and a fixed temperature isothermal wall. These arealso imposed using SAT terms.
3. Stability for wall model: continuous estimate
The wall models (2) and (3) are Neumann and Robin type boundaryconditions, respectively. We will now look at the stability properties of the6all models. We use the compressible Navier-Stokes (CNS) equations in 2Dto derive stability estimates. They are given by ∂U∂t + ∇ · F Ix + ∇ · G Iy − ∇ · F Vx − ∇ · G Vy = 0with F I = (cid:2) ρu, ρu + p, ρuv, u ( p + e ) (cid:3) T G I = (cid:2) ρv, ρuv, ρv + p, v ( p + e ) (cid:3) T F V = [0 , τ xx , τ xy , uτ xx + vτ xy + q x ] T G V = [0 , τ xy , τ yy , uτ xy + vτ yy + q y ] T . with τ ij = µ (cid:18) ∂u i ∂x j + ∂u j ∂x i − δ ij ∂u k ∂x k (cid:19) , q i = κ ∂T∂x i . (5)Using (5), we have the following identities1 µ ( τ xx + 12 τ yy ) = ∂u∂x µ ( τ yy + 12 τ xx ) = ∂v∂y µ τ xy = ∂u∂x + ∂u∂x . (6)We will now look at the kinetic energy contribution only from the viscousterms, and assume that all other terms are bounded. The kinetic energy canbe derived, by multiplying the continuity equation with − u , momentumequation with u and adding them together. Since the continuity equationhas 0 viscous terms, we will only consider the momentum terms. The aboveprocedure gives us ∂k∂t = u ∂τ xx ∂x + u ∂τ xy ∂y + v ∂τ xy ∂x + v ∂τ yy ∂y + conv, (7)where conv refers to the convective terms. For the convective terms, wewill use a kinetic energy conserving scheme. In the following, we will onlyconsider the viscous terms.Assume now that we have a rectangular domain. Integrating (7) over thedomain yields: (cid:90) ∂k∂t dxdy = (cid:90) (cid:18) u ∂τ xx ∂x + u ∂τ xy ∂y + v ∂τ xy ∂x + v ∂τ yy ∂y (cid:19) dxdy. (8)7sing the chain rule, for example, u ∂τ xx ∂x = ∂uτ xx ∂x − τ xx ∂u∂x in (8), results in (cid:90) ∂k∂t dxdy = (cid:90) (cid:18) ∂uτ xx ∂x + ∂uτ xy ∂y + ∂vτ xy ∂x + ∂vτ yy ∂y (cid:19) dxdy − (cid:90) (cid:18) τ xx ∂u∂x + τ xy ∂u∂y + τ xy ∂v∂x + τ yy ∂v∂y (cid:19) dxdy. (9)Consider the second term on the right hand side. The identity (6) yields τ xx ∂u∂x + τ xy ∂u∂y + τ xy ∂v∂x + τ yy ∂v∂y = 1 µ τ xx ( τ xx + 12 τ xx ) + 1 µ τ yy ( τ xx + 12 τ yy )+ 1 µ τ xy τ xy = 1 µ ( τ xx + τ xx τ yy + τ yy + τ xy ) . (10)Now we also have the inequality a + b + ab = ( a + b ) + a + b ≥ . (11)Using (11) in (10) yields τ xx ∂u∂x + τ xy ∂u∂y + τ xy ∂v∂x + τ yy ∂v∂y = 1 µ ( τ xx + τ xx τ yy + τ yy + τ xy ) ≥ . (12)Consequently (12) and (9) leads to (cid:90) ∂k∂t dxdy ≤ (cid:90) (cid:18) ∂uτ xx ∂x + ∂uτ xy ∂y + ∂vτ xy ∂x + ∂vτ yy ∂y (cid:19) dxdy. (13)Note that (13) only contains divergence terms. Therefore, integratingover the domain leaves only boundary terms. We consider a rectangulardomain [0 , and investigate the result at the boundary y = 0. This gives (cid:90) ∂k∂t dxdy ≤ − (cid:90) (cid:0) u ∗ τ ∗ xy + v ∗ τ ∗ yy (cid:1) dx (14)Using the slip boundary condition u ∗ = u, v ∗ = 0 at the wall leads to (cid:90) ∂k∂t dxdy ≤ − (cid:90) u ∗ τ ∗ xy dx. (15)8xpanding (15) with u ∗ = u and τ ∗ xy = µ ( ∂u∂y + ∂v∂x ) (cid:90) ∂k∂t dxdy ≤ − (cid:90) uτ ∗ xy dx = − (cid:90) uµ ( ∂u∂y + ∂v∂x ) dx = − (cid:90) µu ∂u∂y dx. (16)Therefore, to bound kinetic energy, ∂u∂y and u must have opposite signs. If we use the equilibrium wall model, the boundary data is given by τ ∗ xy = u wm,x τ w | u wm | This would lead to (cid:90) ∂k∂t dxdy ≤ − (cid:90) uτ ∗ xy dx = − (cid:90) uu wm,x τ w | u wm | dx. (17)Since, u, u wm,x in general do not have the same sign, the wall model is notkinetic energy stable. This is particularly serious for highly sheared flows. For the dynamic slip wall model, the stability estimate using (16) becomes (cid:90) ∂k∂t dxdy ≤ − (cid:90) µl p u dx. (18)The dynamic slip wall model is therefore kinetic energy stable, since thelength l p >
0. The implementation (4) in the continuous setting using (15)becomes (cid:90) ∂k∂t dxdy ≤ − (cid:90) uτ ∗ xy dx = − (cid:90) u (cid:18) µ ∂u∂y + σ ( u − l p ∂u∂y ) (cid:19) dx. (19)Specifying the penalty parameter σ as σ = µl p in (19) results in the stabilityestimate (cid:90) ∂k∂t dxdy ≤ − (cid:90) µl p u , dx. (20)
4. Numerical Method
In this section we review our specific discontinuous Galerkin method, theassociated SBP property and how it leads to energy stability.9 .1. Notation
An underlined lower case letter denotes a vector quantity, e.g. u . Adoubly underlined lower case letter denotes a vector imposed on the diagonalof a matrix, i.e. u = diag ( u ). A doubly underlined upper case letter denotesany other matrix, e.g. U . A variable with no underline denotes a scalarquantity. Quantities marked as () ∗ , for example, u ∗ denotes boundary data. Consider the linear advection equation. ∂u∂t + ∂u∂x = 0 . (21)If we multiply (21) by u and integrate over the domain, we get (cid:90) Ω uu t dx = ddt (cid:90) Ω u dx = − (cid:90) Ω u ∂u∂x dx = − u | d Ω (22)Therefore, the energy integral over the domain only depends on the valuesof u at the boundaries.Following [12, 25], we can discretize this using the strong form of DG atGauss Lobatto nodes also known as the discontinuous Galerkin collocationspectral element method (DGSEM) as ∂u∂t + D u + M − R T B ( u ∗ − R u ) = 0 . (23)The operators D, R represent the derivative and restriction matrices respec-tively. The restriction matrix can be defined as
R u = (cid:20) u L u R (cid:21) where u L , u R are the left and right values of u at the edge of an element and u ∗ is the numerical flux (eg. Rusanov flux [34]).Summation-By-Parts (SBP) is a discrete analogue of integration by partsthat allows extension of continuous stability estimates to the discrete one’ssee [33]. First, we define a derivative and quadrature matrix D, M such that Du ≈ ∂ x u, u T M v ≈ (cid:90) uvdx. M = diag ( ω , ω , . . . , ω P ), where ω j , j = 0 , . . . , P denotes quadra-ture weights. Next define the boundary matrix as B = diag ( − , R then satisfies R T B R = diag ( − , , .., , . With these definitions, the SBP operator satisfies the following property:
M D + D T M = R T B R. (24)The symmetric positive definite matrix M defines a norm || u || M = u T M u .Integration by parts can be mimicked using these operators since: u T M D v + v T M Du = (
R u ) T B ( R v ) ≈ (cid:90) − u∂ x v + (cid:90) − v∂ x u = uv | − Next, we consider a single element and for a moment ignore the simultaneousapproximation term (SAT). Multiply (23) on the left by 2 u T M which by using(24) leads to 2 u T M ddt u = ddt u = − ( u T M D u + u T D T M u )= − u T R T B R u = u L − u R . (25)Clearly, (25) is the exact discrete equivalent of (22). Suitable boundaryconditions can be imposed using the SAT term. Therefore, we discretelymimicked the continuous stability estimate using the SBP property.The solution of the compressible Navier Stokes(CNS) equations requiresdiscretization of viscous terms. This will be illustrated using the advectiondiffusion equation which can be written as ∂u∂t + ∂u∂x = (cid:15) ∂ u∂x , or ∂u∂t + ∂u∂x = (cid:15) ∂q∂x ,q = ∂u∂x . (26)11quation (26) can be discretized as ∂u∂t + D u + M − R T B ( u ∗ − R u ) = (cid:15) (cid:0)
D q + M − R T B ( q ∗ − R q ) (cid:1) q = D u + M − R T B ( u ∗ − R u ) (27)The advection diffusion equation can have 3 types of boundary conditions.They are Dirichlet u = g , Neumann u x = g or Robin αu − βu x = g where g , g , g is boundary data. Dirichlet boundary condition is imposedusing u ∗ = g and q ∗ = R q . Neumann boundary condition is imposed using q ∗ = g and u ∗ = R u . Robin boundary condition is imposed using u ∗ = R u and q ∗ = R q + σ ( αR u − βR q − g ) . Here σ is a suitably chosen penaltyparameter such as the one chosen in (20).In the interest of keeping the details to a minimum, we end this section bymentioning that we will use the Kennedy and Gruber split form [16] whichis kinetic energy conserving [13], in the Navier Stokes calculations.
5. Stability for wall model: discrete estimate
To get discrete estimates in d dimensions, we will use the tensor productformulation. We use the notation m, d, r, b to denote the 1 dimensional mass,derivative, restriction and boundary matrices. Next we need the Kroneckerproduct. Consider two matrices a ∈ R k × l , b ∈ R m × n Then we define theKronecker product as a ⊗ b := a b a b . . . a l b ... . . . ... a k b a k b . . . a kl b (28)Define the required two dimensional matrices as M = m ⊗ m, D = I ⊗ d, D = d ⊗ IR = I ⊗ r, R = r ⊗ I, B = I ⊗ b B = b ⊗ I (29)In this case the SBP formulation (24) becomes (see [26]) M D i + D iT M = R iT B i R i (30)12e now have the derivative operators D i acting in the i th , i = 1 , .., d di-rection. These operators satisfy the multi-dimensional SBP conditions (see[24]). Next, we look at the discretization of the viscous terms. The stress isdiscretized as τ xx = 43 µD u − µD v + 43 µM − R T B ( u ∗ − R u ) − µM − R T B ( v ∗ − R v ) τ yy = 43 µD v − µD u + 43 µM − R T B ( v ∗ − R v ) − µM − R T B ( u ∗ − R u ) τ xy = µ (cid:16) D u + D v (cid:17) + µM − R T B ( v ∗ − R v ) + µM − R T B ( u ∗ − R u ) (31)where u ∗ , v ∗ are boundary data. To impose the no-penetration boundarycondition, we use u ∗ = R u, v ∗ = 0.The temperature wall boundary condition is imposed in the energy equa-tion which does not contribute to the kinetic energy bound. The isothermalwall boundary condition is imposed as T ∗ = C where C is the desired con-stant value. From (31) it is obvious that τ xx + 12 τ yy = µD u + µM − R T B ( u ∗ − R u ) τ yy + 12 τ xx = µD v + µM − R T B ( v ∗ − R v ) (32)Stress gradients are given by τ xx,x = D τ xx + M − R T B ( τ ∗ xx − R τ xx ) τ xy,x = D τ xy + M − R T B ( τ ∗ xy − R τ xy ) τ xy,y = D τ xy + M − R T B ( τ ∗ xy − R τ xy ) τ yy,y = D τ yy + M − R T B ( τ ∗ yy − R τ yy ) . (33)where ( · ) ,x denotes the derivative with respect to x . We can impose theboundary data (2), (4) using τ ∗ xy in (33). For τ ∗ xx , τ ∗ yy we use τ ∗ xx = R τ xx , τ ∗ yy =13 τ yy . Recall that we aim for the discrete version of kinetic energy (8). ddt T M k = u T M τ xx,x + u T M τ xy,y + v T M τ xy,x + v T M τ yy,y . (34)Using the SBP property and velocity boundary data, we get ddt T M k ≤ τ xyT R T B R u + u T R T B ( τ ∗ xy − R τ xy )= u T R T B τ ∗ xy . (35)which is the discrete equivalent of (15). The full derivation is presented inAppendix A and the conclusion is that the discrete kinetic energy stabilityestimates mimic the continuous kinetic energy estimates.
6. Results
In this section, we will present numerical results. The properties of theunderlying DGSEM and split form is shown in [30]. First, the isentropicvortex test case is used to show the order of convergence of the numericalmethod. Next, we use the turbulent channel flow to compare the accuracyof the wall models. We then consider a simple 2D airfoil to show that kineticenergy instability of the wall model can lead to solution blow-up.
The 2D isentropic vortex problem is an exact solution for the compressibleEuler equations. It can be used to test the order of accuracy of a numericalmethod. There are various initializations reported for this problem, suchas shown in [31]. In this study, the domain [ − , is considered with thefollowing initial condition: ρ = (cid:20) − β ( γ − γπ exp(1 − r ) (cid:21) u = M − β ( y − y c )2 π exp( 1 − r )2 v = β ( y − y c )2 π exp( 1 − r r = (cid:112) ( x − x c ) + ( y − y c ) . Also, p = ρ γ , β = 5 , γ = 1 .
4, where β determines vortex strength. Further, ( x c , y c ) = (0 , , M = 0 .
5, and periodic14 able 1: L error and order of convergence for the isentropic vortex. P = 3 P = 4N L error OC L error OC6 1.03E-2 - 1.92E-3 -12 8.19E-4 3.65 1.04E-4 4.2124 7.49E-5 3.45 5.32E-6 4.2948 7.34E-6 3.35 2.48E-7 4.4296 5.00E-7 3.88 1.26E-8 4.30boundary conditions are imposed. All simulations are performed at CF L =0 . N elements, after one flow through. Two polynomial order, P = 3 , P but slightly less than the expected orderof P + 1. This is due to the fact that Gauss Lobatto points generally leadsto lower accuracy than Gauss Legendre points [18]. One of the objectives in developing high order methods for unstructuredgrids is to enable simulation of wall bounded turbulent flows. The simplestpossible test case for this is turbulent channel flow [17]. This flow has aperiodic domain in the streamwise and spanwise directions with walls at thetop and bottom. The flow is driven by a forcing term to get a constant massflow rate. The test case is characterized by the parameter Re τ given by Re τ = u τ Lν , u τ = (cid:114) τ w ρ . DNS studies have been done for this case at Re τ = 5200 [17, 20]. Studies havealso been done using spectral type methods such as flux reconstruction (FR)and DG [35, 21]. We use a domain length 2 π in the streamwise direction and π in the spanwise direction. The channel height is H = 2. It is discretized into22 × ×
12 uniform elements. Figure 2 shows the mesh. The computationswere compared to DNS data from [20].The initialization and forcing is similar to [21], with the streamwise ve-locity initialized using the profile, u = 95 u (1 − ( y − ) . igure 2: Uniform mesh for channel flow.Figure 3: Instantaneous snapshot showing horizontal velocity contours for flow with equi-librium wall model The vertical velocity in initialized using v = 0 . u exp (cid:34) − (cid:18) x − π π (cid:19) (cid:35) exp (cid:20) − (cid:16) y (cid:17) (cid:21) cos(4 z ) . The initial density is taken as ρ = 1 and the pressure chosen such thatthe Mach number for velocity u is M = 0 .
2. To this uniform profile, per-turbations are added to trigger transition. They are of the form u (cid:15) = 0 . λ πy λ πz igure 4: Time averaged mean profiles of streamwise velocity at Re τ = 5200 for wallmodelled ILES (WMILES) at P=3 using equilibrium stress wall model compared withDNS from [20]. where λ is a suitable parameter, eg. λ = 20 , ,
40. To enable constant massflow rate, a source term is added to the right hand side of (3) of the form S = (0 , s , , , s ) T , s = F w V + 0 . dt ( m Re − m + m ) , s = u b s (36)Here F w , m , m, m Re are the shear at the wall, mass flow at the previoustime step, mass flow at the current time step and mass flow required for theprescribed Re τ respectively. The first term in s balances the total shearstress generated at the wall, while the second term accelerates convergence.17 igure 5: Time averaged mean profiles of second moments at Re τ = 5200 for wall modelledILES (WMILES) at P=3 using equilibrium stress wall model compared with DNS from[20]. The term s is there simply to be energy consistent with the forcing where u b is the average streamwise velocity over the domain. Figure 3 shows horizontalvelocity contours for this flow with the equilibrium stress wall model. Theflow is clearly highly turbulent. In addition, the boundary layer is very thinas expected for a very high Reynolds’ number flow. Figure 4 shows the mean velocity profile averaged over time for Re τ =5200 with the equilibrium stress wall model. Notice first that the first two18 igure 6: Time averaged mean profiles of streamwise velocity at Re τ = 5200 for wallmodelled ILES (WMILES) using dynamic slip wall model compared with DNS from [20]. points are slightly off from the DNS curve. Recall that in the wall model,the tangential velocity at the top solution point in the wall adjacent elementis taken as the input for the shear stress imposed at the wall. Therefore,the top point is assumed to be the first point that is within the resolvedregion of the flow. This indicates that the first 3 points probably do notmatch the DNS data. We notice that from the third point onwards, ourwall modeled data compares well with the DNS data. This is despite thefact that our first solution point is at y + ∼
200 which is in the log-layer.Therefore, with the wall model we see that the statistics of the first momentare well captured. Figure 5 shows the second moments. Similar to the meanvelocity, the simulations depart near the boundary but away from the wall,they compare well with DNS data. 19 igure 7: Time averaged mean profiles of second moments at Re τ = 5200 for wall modelledILES (WMILES) using dynamic slip wall model compared with DNS from [20]. Next we use the dynamic slip wall model and compare it to DNS data.We use l p = C ∆ y where C, ∆ y denotes a constant and the height of the firstelement respectively. We test three cases with C = 0 . , . , .
50. Figure 6compares the mean velocity profile with DNS data. We note that there is alarge mismatch between wall modeled results and DNS data. However, as thevalues of C is being reduced, the results gradually improve but the mismatchpersists. Figure 7 shows the second moments, which also deviates from DNSdata. This indicates that the choice of a constant l p is too simplistic. The NACA4412 airfoil is a popular test case for various experimental andnumerical studies [8, 15, 11]. Flow over an airfoil is characterized by changing20 igure 8: Horizontal velocity contours over NACA4412 with equilibrium stress wall modelboundary condition. Note the unstable flow over both surfaces.Figure 9: Horizontal velocity contours over NACA4412 with dynamic slip boundary con-dition. velocity and pressure gradients unlike in the channel flow. We use this testcase to demonstrate that there are stability problems with the wall model asour analysis indicates. We have previously shown the effectiveness of usingthe equilibrium wall model for this airfoil [30]. The wall model was used forflow at Re ≈ . × , M = 0 . , α = 12 ◦ . In this paper, we want to show21 igure 10: Close-up of horizontal velocity contours over NACA4412 with equilibrium stresswall model boundary condition. Note the high negative velocities near the wall.Figure 11: Close-up horizontal velocity contours over NACA4412 with dynamic slip wallmodel boundary condition. here that with increasing Reynolds’ and Mach number, stability issues canmake the equilibrium wall model unsuitable. To show this, we take a 2Dunstructured grid and simulate flow at Re ≈ . × , M = 0 . , α = 5 ◦ ,where α is the angle of attack. The mach number and angle of attack areonly moderately high and so should not be challenging for a stable solver.22e ran the simulation to test robustness. With the equilibrium stresswall model, the simulation only ran for 4 flow-through times and then blew-up. Whereas for the dynamic slip wall model, the simulation did not blowup even after 30 flow through time periods. Figures 8 and 9 show horizontalvelocity over the airfoil for the two wall models. For the equilibrium stresswall model, we see irregular flow features all over the surface in contrast tothe dynamic slip wall model. A more detailed illustration of these issuescan be seen in figures 10 and 11. For the equilibrium stress wall model, weobserve very high negative velocity near the wall. On the other hand, for thedynamic slip wall model, the near wall velocities are close to 0 .The high negative velocity in the equilibrium stress wall model occurssince the shear stress is not guaranteed to oppose the velocity at the wall,but rather a velocity that is away from the wall. This leads to kinetic energyincrease and eventually to blow up. We have also investigated a hybrid wall model which choses one of thetwo wall models based on a particular condition. This condition guaranteeskinetic energy stability. τ ∗ xy = (cid:40) u wm,x τ w | u wm | , if u · u wm,x > τ xy + σ ( u − l p ∂u∂y ) , otherwise (37)where u, u wm,x are the components of the velocity in the x direction at y =0 , y = y wm respectively. Note that equation (37) selects (3) or (2) dependingon the condition for the velocities. Shear stress boundary values for τ ∗ yz canbe similarly chosen.Figure 12 compares the mean velocity profiles for the three wall modelsagainst DNS data. We select the parameter l p = 0 . y since that was theclosest to the DNS profile for the dynamic slip wall model. The hybrid wallmodel curve is much closer to the DNS curve than the dynamic slip wallmodel. However, the equilibrium wall model is still more accurate. Figure13 compares the second moments. Again, we observe similar features as themean velocity. This suggests that using a hybrid wall model is a possiblesolution to a stable and accurate wall model. For the sheared flow test case,this wall model behaves very similar to the dynamic slip wall model.23 igure 12: Time averaged mean profiles of streamwise velocity at Re τ = 5200 for wallmodelled ILES (WMILES) using all three wall models compared with DNS from [20].
7. Conclusion
In the present study we investigate stability of wall modelled DG for thecompressible Navier Stokes equations. First, the formulations and the algo-rithm for implementing the wall model have been introduced. Next, kineticenergy stability estimates in the continuous case were developed, which showthat the wall shear stress must oppose the velocity at the wall. Further, weshow that this estimate can be extended to the discrete DGSEM formula-tion because of its SBP property. This framework is used to investigate twowall models. The popular equilibrium stress wall model is not kinetic energystable, whereas the dynamic slip wall model is stable.Next, two test cases are studied numerically to contrast the accuracy andstability properties of the wall models. First, we use the turbulent channelflow to show that at high Reynolds’ number, the first moment is accurately24 igure 13: Time averaged mean profiles of second moments at Re τ = 5200 for wall mod-elled ILES (WMILES) using all three wall models compared with DNS from [20]. captured by the equilibrium stress wall model away from the wall. However,the dynamic slip wall model is not very accurate. Next, the NACA4412 testcase is run at high Reynolds; number but at moderate mach number andangle of attack. With the equilibrium stress wall model, the flow becomesvery irregular and blows-up very quickly. Near the wall, there are very largenegative velocities. This occurs because the wall model does not satisfy thekinetic energy stability estimate. However, with the dynamic slip wall model,the flow is smooth and stable. Finally, a hybrid wall model is shown thatchoses the equilibrium wall model when it is stable and the dynamic slip wallmodel otherwise. This is shown to behave much better in terms of accuracyto the dynamic slip wall mdoel. This shows a possible way forward towardsaccurate and stable wall modeling. To enable wider use of the kinetic energystable dynamic slip wall model, better parameterizations are required for25ore accuracy. Acknowledgments
Vikram Singh and Steven Frankel would like to acknowledge that thiswork was partially supported by Israel Science Foundation ISF-NSFC jointresearch program (ISF Grant No. 2232/15). Jan Nordstr¨om was supportedby Vetenskapsr˚adet (award number 2018-05084 VR), Sweden.
Appendix A. Discrete kinetic energy estimate
Here, we derive a discrete estimate for the kinetic energy for the viscousterms. Using (33) in (34) yields ddt T M k = u T M D τ xx + u T R T B ( τ ∗ xx − R τ xx )+ u T M D τ xy + u T R T B ( τ ∗ xy − R τ xy )+ v T M D τ xy + v T R T B ( τ ∗ xy − R τ xy )+ v T M D τ yy + v T R T B ( τ ∗ yy − R τ yy )= τ xxT D T M u + u T R T B ( τ ∗ xx − R τ xx )+ τ xyT D T M u + u T R T B ( τ ∗ xy − R τ xy )+ τ xyT D T M v + v T R T B ( τ ∗ xy − R τ xy )+ τ yyT D T M v + v T R T B ( τ ∗ yy − R τ yy ) . (A.1)since terms like u T M D τ xx are scalars implying u T M D τ xx = τ xxT D T M u.
If we multiply the first equation in (32) with τ xxT M and the second with τ yyT M from the left, we get τ xxT M ( τ xx + 12 τ yy ) = µτ xxT M D u + µτ xxT R T B ( u ∗ − R u ) τ yyT M ( τ yy + 12 τ xx ) = µτ xxT M D v + µτ xxT R T B ( v ∗ − R v ) . (A.2)Using the SBP property (30) in (A.2) yields1 µ τ xxT M ( τ xx + 12 τ yy ) = − τ xxT D T M u + τ xxT R T B u ∗ µ τ yyT M ( τ yy + 12 τ xx ) = − τ xxT D T M v + τ xxT R T B v ∗ . (A.3)26imilarly, multiplying the expression for τ xy on the left by τ xyT M in (31)yields1 µ τ xyT M τ xy = τ xyT M D u + τ xyT M D v + τ xyT R T B ( v ∗ − R v )+ τ xyT R T B ( u ∗ − R u ) (A.4)Again using the SBP property (30) in (A.4) yields1 µ τ xyT M τ xy = − τ xyT D T M u − τ xyT D T M v + τ xyT R T B v ∗ + τ xyT R T B u ∗ (A.5)Using (A.3) and (A.5) in the expression for kinetic energy (A.1) yields ddt T M k = − µ τ xyT M τ xy + τ xyT R T B v ∗ + τ xyT R T B u ∗ − µ τ xxT M ( τ xx + 12 τ yy ) − µ τ yyT M ( τ yy + 12 τ xx )+ τ xxT R T B u ∗ + τ yyT R T B v ∗ + u T R T B ( τ ∗ xx − R τ xx ) + u T R T B ( τ ∗ xy − R τ xy )+ v T R T B ( τ ∗ xy − R τ xy ) + v T R T B ( τ ∗ yy − R τ yy ) . (A.6)Using the property that M is a positive definite matrix and the identity (11)yields1 µ τ xxT M ( τ xx + 12 τ yy ) + 1 µ τ yyT M ( τ yy + 12 τ xx ) + 1 µ τ xyT M τ xy =1 µ ( τ xyT M τ xy + τ xxT M τ xx + τ xxT M τ yy + τ yyT M τ yy ) ≥ ddt T M k ≤ τ xyT R T B v ∗ + τ xyT R T B u ∗ + τ xxT R T B u ∗ + τ yyT R T B v ∗ + u T R T B ( τ ∗ xx − R τ xx ) + u T R T B ( τ ∗ xy − R τ xy )+ v T R T B ( τ ∗ xy − R τ xy ) + v T R T B ( τ ∗ yy − R τ yy ) . (A.8)27s before, if we look at only the y boundary (A.8) becomes ddt T M k ≤ τ xyT R T B u ∗ + τ yyT R T B v ∗ + u T R T B ( τ ∗ xy − R τ xy )+ v T R T B ( τ ∗ yy − R τ yy ) . (A.9)Substituting the boundary condition for slip velocity v ∗ = 0 , u ∗ = R u in(A.10) yields ddt T M k ≤ τ xyT R T B R u + u T R T B ( τ ∗ xy − R τ xy )= u T R T B τ ∗ xy . (A.10) References [1]
H. J. Bae, A. Lozano-Dur´an, S. T. Bose, and P. Moin , Dynamicslip wall model for large-eddy simulation , Journal of Fluid Mechanics,859 (2018), pp. 400–432.[2]
F. Bassi and S. Rebay , A high-order accurate discontinuous fi-nite element method for the numerical solution of the compressibleNavier–Stokes equations , Journal of Computational Physics, 131 (1997),pp. 267 – 279.[3]
A. D. Beck, T. Bolemann, D. Flad, H. Frank, G. J.Gassner, F. Hindenlang, and C.-D. Munz , High-order discon-tinuous Galerkin spectral element methods for transitional and turbulentflow simulations , International Journal for Numerical Methods in Fluids,76 (2014), pp. 522–548.[4]
J. Berg and J. Nordstr¨om , Stable Robin solid wall boundaryconditions for the Navier–Stokes equations , Journal of ComputationalPhysics, 230 (2011), pp. 7519 – 7532.[5]
S. T. Bose and P. Moin , A dynamic slip boundary condition for wall-modeled large-eddy simulation , Physics of Fluids, 26 (2014), p. 015104.[6]
H. Choi and P. Moin , Grid-point requirements for large eddy sim-ulation: Chapman’s estimates revisited , Physics of Fluids, 24 (2012),p. 011702. 287]
B. Cockburn, G. E. Karniadakis, and C.-W. Shu , The devel-opment of discontinuous galerkin methods , in Lecture Notes in Com-putational Science and Engineering, Springer Berlin Heidelberg, 2000,pp. 3–50.[8]
D. Coles and A. J. Wadcock , Flying-hot-wire study of flow past anNACA 4412 airfoil at maximum lift , AIAA Journal, 17 (1979), pp. 321–329.[9]
C. C. de Wiart and S. M. Murman , Assessment of wall-modeledLES strategies within a discontinuous-galerkin spectral-element frame-work , in 55th AIAA Aerospace Sciences Meeting, American Institute ofAeronautics and Astronautics, Jan. 2017.[10]
A. Fr`ere, C. C. de Wiart, K. Hillewaert, P. Chatelain,and G. Winckelmans , Application of wall-models to discontinuousGalerkin LES , Physics of Fluids, 29 (2017), p. 085111.[11]
A. Fr`ere, K. Hillewaert, P. Chatelain, and G. Winckelmans , High Reynolds number airfoil: From wall-resolved to wall-modeled LES ,Flow, Turbulence and Combustion, 101 (2018), pp. 457–476.[12]
G. J. Gassner , A skew-symmetric discontinuous Galerkin spectral el-ement discretization and its relation to SBP-SAT finite difference meth-ods , SIAM Journal on Scientific Computing, 35 (2013), pp. A1233–A1253.[13]
G. J. Gassner, A. R. Winters, and D. A. Kopriva , Split formnodal discontinuous Galerkin schemes with summation-by-parts propertyfor the compressible Euler equations , Journal of Computational Physics,327 (2016), pp. 39 – 66.[14]
J. S. Hesthaven and T. Warburton , Nodal Discontinuous GalerkinMethods , Springer New York, 2008.[15]
S. Hosseini, R. Vinuesa, P. Schlatter, A. Hanifi, and D. Hen-ningson , Direct numerical simulation of the flow around a wing sectionat moderate reynolds number , International Journal of Heat and FluidFlow, 61 (2016), pp. 117–128. 2916]
C. A. Kennedy and A. Gruber , Reduced aliasing formulations of theconvective terms within the Navier–Stokes equations for a compressiblefluid , Journal of Computational Physics, 227 (2008), pp. 1676 – 1700.[17]
J. Kim, P. Moin, and R. Moser , Turbulence statistics in fully devel-oped channel flow at low reynolds number , Journal of Fluid Mechanics,177 (1987), p. 133.[18]
D. A. Kopriva and G. Gassner , On the quadrature and weak formchoices in collocation type discontinuous Galerkin spectral element meth-ods , Journal of Scientific Computing, 44 (2010), pp. 136–155.[19]
J. Larsson, S. Kawai, J. Bodart, and I. Bermejo-Moreno , Large eddy simulation with modeled wall-stress: recent progress andfuture directions , Mechanical Engineering Reviews, 3 (2016), pp. 15–00418–15–00418.[20]
M. Lee and R. D. Moser , Direct numerical simulation of turbulentchannel flow up to Re τ ≈ G. Lodato, P. Castonguay, and A. Jameson , Discrete filter oper-ators for large-eddy simulation using high-order spectral difference meth-ods , International Journal for Numerical Methods in Fluids, 72 (2012),pp. 231–258.[22]
M. Parsani, M. H. Carpenter, and E. J. Nielsen , Entropystable wall boundary conditions for the three-dimensional compressibleNavier–Stokes equations , Journal of Computational Physics, 292 (2015),pp. 88 – 113.[23]
S. Pirozzoli , Numerical methods for high-speed flows , Annual Reviewof Fluid Mechanics, 43 (2011), pp. 163–194.[24]
H. Ranocha , SBP operators for CPR methods: Master’s thesis , (2016).[25]
H. Ranocha, P. ¨Offner, and T. Sonar , Summation-by-parts op-erators for correction procedure via reconstruction , Journal of Compu-tational Physics, 311 (2016), pp. 299 – 328.3026] ,
Extended skew-symmetric form for summation-by-parts operatorsand varying Jacobians , Journal of Computational Physics, 342 (2017),pp. 13 – 28.[27]
W. Reed and T. Hill , Triangular mesh methods for the neutron trans-port equation , (1973).[28]
H. Reichardt , Vollst¨andige darstellung der turbulentengeschwindigkeitsverteilung in glatten leitungen , ZAMM - Zeitschrift f¨urAngewandte Mathematik und Mechanik, 31 (1951), pp. 208–219.[29]
H. Schlichting and K. Gersten , Boundary-Layer Theory , SpringerBerlin Heidelberg, 2017.[30]
V. Singh and S. Frankel , On the use of split forms and wall modelingto enable accurate high-reynolds number discontinuous galerkin simula-tions on body-fitted unstructured grids , Computers & Fluids, 208 (2020),p. 104616.[31]
S. C. Spiegel, J. R. DeBonis, and H. Huynh , Overview of theNASA Glenn flux reconstruction based high-order unstructured grid code ,in 54th AIAA Aerospace Sciences Meeting, American Institute of Aero-nautics and Astronautics, jan 2016.[32]
M. Sv¨ard and J. Nordstr¨om , A stable high-order finite differencescheme for the compressible Navier–Stokes equations , Journal of Com-putational Physics, 227 (2008), pp. 4805–4824.[33]
M. Sv¨ard and J. Nordstr¨om , Review of summation-by-partsschemes for initial–boundary-value problems , Journal of ComputationalPhysics, 268 (2014), pp. 17 – 38.[34]
E. F. Toro , Riemann Solvers and Numerical Methods for Fluid Dy-namics , Springer Berlin Heidelberg, 2009.[35]