On the extension of a Riemann solver for RANS simulations

Abstract

Purpose

Reynolds-averaged Navier–Stokes (RANS) models often perform poorly in shock/turbulence interaction regions, resulting in excessive wall heat load and incorrect representation of the separation length in shockwave/turbulent boundary layer interactions. The authors suggest that this can be traced back to inadequate numerical treatment of the inviscid fluxes. The purpose of this study is an extension to the well-known Harten, Lax, van Leer, Einfeldt (HLLE) Riemann solver to overcome this issue.

Design/methodology/approach

It explicitly takes into account the broadening of waves due to the averaging procedure, which adds numerical dissipation and reduces excessive turbulence production across shocks. The scheme is derived based on the HLLE equations, and it is tested against three numerical experiments.

Findings

Sod’s shock tube case shows that the scheme succeeds in reducing turbulence amplification across shocks. A shock-free turbulent flat plate boundary layer indicates that smooth flow at moderate turbulence intensity is largely unaffected by the scheme. A shock/turbulent boundary layer interaction case with higher turbulence intensity shows that the added numerical dissipation can, however, impair the wall heat flux distribution.

Originality/value

The proposed scheme is motivated by implicit large eddy simulations that use numerical dissipation as subgrid-scale model. Introducing physical aspects of turbulence into the numerical treatment for RANS simulations is a novel approach.

Citation

Buck, A. and Mundt, C. (2024), "On the extension of a Riemann solver for RANS simulations", International Journal of Numerical Methods for Heat & Fluid Flow, Vol. ahead-of-print No. ahead-of-print. https://doi.org/10.1108/HFF-11-2023-0708

Nomenclature

c_s

= wave speed factor [−];

F

= mean flux vector;

k

= turbulence kinetic energy [m²/s²];

P

= turbulence production [m²/s³];

q

= heat flux [w/m²];

S

= wave speed [m/s];

t

= time [s];

Tu

= turbulence intensity [−];

u_i

= velocity components [m/s];

W

= mean state vector;

x_i

= spacial coordinates [m];

ρ

= density [kg/mm/³];

^∼

= Favre average; and

″

= Favre fluctuation.

Subscripts

L, R

= Left and right value;

T

= turbulent; and

W

= wall.

1. Introduction

Most real-world flows are in the turbulent regime, so treating the effects of turbulence correctly is critical for obtaining accurate results. Reynolds-averaged Navier–Stokes (RANS) turbulence models have a persisting relevance for engineering applications due to their efficiency compared with higher order approaches such as large eddy simulations (LES) and direct numerical simulations (DNS). Furthermore, engineers are often only interested in time-averaged solutions, rather than resolving the turbulent fluctuations in time, so computing the mean solution directly is practical.

Turbulence and shock/turbulence interaction (STI) is a critical phenomenon for many aerospace vehicles such as launchers, airplanes and hypersonic vehicles. While RANS models perform robustly in certain flow conditions, some physical aspects, such as STI, are usually not captured properly. Turbulence production in STI is generally over-predicted compared to DNS data. Sinha et al. (2003) showed that post-shock turbulence kinetic energy is excessive with standard two equation models. The result of this depends on the chosen turbulence model formulation with some models producing excessive heat load and decreased separation length (Knight et al., 2002), whereas separation length increases with other models (Brown, 2013). For compressible flows, a combination of Reynolds decomposition ϕ=ϕ¯+ϕ′ with ϕ¯=lim⁡τ→∞1/τ∫tt+τϕ(x→,t)dt for the pressure and density, and Favre decomposition ϕ=ϕ˜+ϕ″ with ϕ˜=ρϕ¯/ρ¯ for the velocities yields the most compact set of equations. Note that ρ¯ is the Reynolds-averaged density. The exact production term P for the turbulence kinetic energy (TKE) in the RANS equations can then be written as:

(1) P=−ρui″uj″¯∂u˜i∂xj.

This is modeled in RANS eddy-viscosity models as (Einstein notation is applied):

(2) PEVM=[μT(∂u˜i∂xj+∂u˜j∂xi−23∂u˜k∂xkδij)−23ρ¯k]∂u˜i∂xj,

where u_i is the velocity in x_i direction, δ_ij is the Dirac delta, k=1/2ui″ui″¯ is the turbulence kinetic energy. µ_T is the eddy viscosity, which must be determined by the turbulence model. For a 1D case, we find:

(3) PEVM,1D=μT43(∂u˜∂x)2∝(∂u˜∂x)2.

In shocks, the magnitude of the mean velocity gradient ∂u˜/∂x is large and so is the turbulence production. Turbulent flows are inherently unsteady, i.e. they are characterized by small-scale fluctuations, which also affect the location of shock waves. RANS solutions are statistical averages, and while shocks are instantaneously discontinuous, in a time-averaged solution they appear spread out in space (Sinha et al., 2003). Because of that, the mean velocity gradients are reduced compared with laminar or inviscid flows, which also decreases the turbulence production.

To improve existing turbulence models, Sinha et al. adapted the production term to account for shock unsteadiness, arguing that the Boussinesq hypothesis breaks down near shock waves (Sinha et al., 2003). While this proved to be effective in their test cases, it is a fundamental change in the model. The production term in equation (1), is exact and, given that the Boussinesq hypothesis is correct, equation (2) is an exact extension. We suggest that the shortcomings of existing models can also be interpreted as the result of excessively steep mean flow gradients in the numerical solution of RANS models. We propose that the current production term can give acceptable results when paired with adequate numerical methods that take into account the statistical nature of RANS solutions, i.e. spread shocks in space sufficiently.

For LES, the idea of combining the effect of shock-capturing numerics and subgrid-scale models has been around for a while. It was realized that methods for capturing steep gradients and subgrid-scale models both add dissipative fluxes, so it seemed reasonable to combine the two. Boris et al. (1992) were early proponents of the implicit LES (ILES) approach that aims to eliminate the need for an explicit subgrid-scale model by applying a numerical scheme that dissipates energy on the small scales. Adams and Stolz (2002) further added to the method of combining the effects of subgrid-scale modeling and shock-capturing by means of approximate deconvolution. More recently, Sousa and Scalo (2022) proposed to unify the artificial viscosity of shock-capturing schemes with the subgrid-scale model using their quasi-spectral viscosity method.

For RANS simulations, this consolidation of turbulence modeling and numerical treatment is not common, usually the two are treated separately. Compared to LES, where the subgrid-scale model is primarily required to dissipate kinetic energy, RANS models are supposed to emulate the effect of all turbulent scales. Therefore, the same level of inclusion of turbulence model and numerics seems challenging for RANS simulations. However, there has been some work adapting existing numerical schemes for the application in RANS simulations to improve accuracy and stability. Chuvakhov (2014) extended the Roe scheme for two-dimensional 2-eq q-ω models. The author derived Roe matrices for the convective fluxes of q-ω and k-ω models, resulting in better convergence and stability.

The RANS equations can be derived from the Navier–Stokes equations by decomposing the field variables into a mean and a fluctuating part, and then averaging the equations. Inviscid and viscous fluxes can be distinguished in the resulting set of equations. Generally, these are discretized using different methods; flux-difference splitting methods (Riemann solvers) are popular for the inviscid fluxes. These methods were developed for the Euler equations and not for the statistical representation of RANS solutions. Because the unsteadiness of turbulent flows is not a viscous effect, it should also be included in the treatment of the inviscid fluxes. Furthermore, with existing shock-capturing schemes, shock waves are captured over a fixed number of cells. Because of that, mesh independence cannot be achieved: as the mesh is refined, shocks become sharper and turbulence production is increased (Lacombe et al., 2021).

In this publication, we extend the popular Harten, Lax, van Leer, Einfeldt (HLLE) Riemann solver (Harten et al., 1983) to better represent the averaged solution in turbulent flows, which should result in a more physical representation of the flow features in RANS solutions. Due to the incoming turbulent fluctuations, the positions of all waves are unsteady. The time-averaged result of the waves are therefore spread out in space. Their size depends on the local turbulence intensity. By treating the waves as fans with finite width in the Riemann solver, additional numerical dissipation is added to the solution. This decreases mean flow gradients and hence reduces turbulence production across steep gradients, such as shocks.

This publication is structured as follows:

In Section 2, we give an introduction to Riemann solvers and derive the proposed alterations for RANS models, before showing the governing equations and numerical methods of the implementation in Section 3. The cases that are used to investigate the scheme and the results are presented in Section 4.

2. Riemann solvers

When using the finite-volume method to discretize any conservation equation in space, approximations for the flux between cells must be found. See Figure 1 for an illustration of the location of known values of the state vector W and the required flux vector F in cell-centered finite-volume schemes. Following the method originally proposed by Godunov and Bohachevsky (1959), this inter-cell flux can be found by solving the local Riemann problem, where the adjacent cell values represent the left and right initial conditions. In addition to finite-volume schemes, the proposed Riemann solver may also be used with other types of spacial discretization, such as finite-difference and discontinuous Galerkin schemes.

Exact solutions of the Riemann problem exist, but due to the nonlinearity of the Euler equations, they must be found iteratively. A more efficient alternative are approximate Riemann solvers, such as the HLLE (Harten et al., 1983), or Harten, Lax, van Leer, Contact (HLLC) (Toro et al., 1994) schemes. Instead of solving the exact Riemann problem iteratively, an approximate solution is computed directly.

2.1 Standard Riemann solvers

Riemann solvers are a numerical method to solve the Riemann problem. They are generally used for systems of hyperbolic conservation laws to correctly represent shock waves in the solution. A generic 1D conservation law may be written in conservative form:

(4) ∂W∂t+∂F(W)∂x=0.

As the solution W may be discontinuous, the integral form of equation (4) over a closed area Ω bounded by ∂Ω in x − t space is more useful. Using the divergence theorem it can be written as:

(5) ∬Ω(∂W∂t+∂F(W)∂x)dtdx=∮∂Ω(Wdx−Fdt)=0.

The local Riemann problem is constructed by placing the cell interface at x = 0. We now want to find the solution W* = W(x = 0, t > 0) and F* = F(x = 0, t > 0) of the local Riemann problem for the initial distribution:

(6) W(x,t=0)={WL, for x<0WR, for x>0.

Without any further assumptions, the integral in equation (5) can only be solved for the inter-cell flux F^* iteratively. Approximate Riemann solvers, such as the HLLE scheme, are derived by first assuming that the solution consists of one left-running and one right-running wave with corresponding wave speeds S_L and S_R. Because it forms the basis of the proposed scheme, its main equations shall be presented briefly. Two cases must be treated separately: left and right supersonic flow and subsonic flow. The solution for the supersonic cases is simple:

(7) W∗=WL and F∗=FL if SL>0,

(8) W∗=WR and F∗=FR if SR<0.

Figure 2(a) shows the wave pattern for the subsonic case, where the dashed line shows the integration domain in the x − t space and the solid lines represent the left and right characteristic waves. Because the state vector is piecewise constant along quadrilateral ABCD, integration simply yields:

(9) W∗=SRWR−SLWL−(FR−FL)SR−SL.

The inter-cell flux F^* can either be determined as F(W^*), or directly by integrating equation (5) over one the quadrilaterals AEFD or EBCF, which gives:

(10) F∗=SRFL−SLFR+SLSR(WR−WL)SR−SL.

2.2 Proposed Riemann solver for turbulent flows

To enhance the capability of the HLLE scheme to capture the mean flow gradients in statistically averaged turbulent flows, both waves are treated as compression or expansion fans. Figure 2(b) shows the wave pattern in an averaged turbulent flow, where the inherent unsteadiness spreads all waves in space. Instead of discontinuous jumps, the fans represent gradual changes in W and F. The left and right fans have the respective limiting mean wave speeds S_j,min and S_j,max with the mean wave speed variance ΔS_j = S_j,max − S_j,min (for j = L or R). The states outside the fans are still denoted as:

(11) W=WL and F=FL for x<tSL,min

and:

(12) W=WR and F=FR for x>tSR,max.

Similar to inviscid Riemann solvers, multiple flow cases must be distinguished. Left and right supersonic cases are treated analogously:

(13) W∗=WL and F∗=FL if SL,min>0,

(14) W∗=WR and F∗=FR if SR,max<0.

The subsonic case is shown in Figure 2(b), where S_L,max < 0 and S_R,min > 0. To find an explicit solution to the Riemann problem in this case, we assume that W varies linearly in the fans along t = const:

(15) W(x,t)=W(tSj,min,t)+W(tSj,max,t)−W(tSj,min,t)tSj,max−tSj,min(x−tSj,min)  for  tSj,min<x<tSj,max,

where j = L, R. The result of the integration with this assumption is also correct for all distributions that are symmetric about the center of the fans. More complex, nonsymmetric distributions could be investigated in future studies. Integrating W on a line with t = const across a fan gives:

(16) ∫tSj,mintSj,maxW(x,t)dx=12[W(tSj,min,t)+W(tSj,max,t)](tSj,max−tSj,min).

We can now integrate equation (5) over quadrilateral ABCD to find the inter-cell state analogously to the HLLE scheme:

(17) W∗=−2FL+2FR+WL(SL,max+SL,min)−WR(SR,max+SR,min)SL,max+SL,min−SR,max−SR,min.

We find the inter-cell flux by integrating over AEFD:

(18) F∗=12(2FL+(W∗−WL)(SL,max+SL,min))

and inserting equation (17):

(19) F∗=2FR(SL,max+SL,min)−(SR,max+SR,min)(2FL−(SL,max+SL,min)(WL−WR))2(SL,max+SL,min−SR,max−SR,min).

Note that equation (19) reduces to equation (10) as intended, when S_L,min = S_L,max = S_L and S_R,min = S_R,max = S_R. The integration in the x − t space is equivalent for both cases except for the sections DD’ and C’C.

In addition to subsonic and supersonic cases, a right transonic case with S_L,min < 0 and S_L,max > 0 and a left transonic case with S_R,min < 0 and S_R,max > 0 exist for turbulent flows. Figure 3 shows the wave pattern for the two cases. For the sake of brevity, we use W^* for the state between the fans and F^* = F(0, t > 0) for the flux along the positive t-axis. W^* in the transonic case is equivalent to the one in the subsonic case. The inter-cell flux in the right transonic case is found by integrating over quadrilateral ABCD. As we assume that the state vector varies linearly along DC, see equation (15), this gives:

(20) ∫ΔtSL,min0W(x,Δt)dx=W∗−WLSL,max−SL,minSL,min22Δt−WLSL,minΔt.

After performing the integration over ABCD, the inter-cell flux is:

(21) F∗=2FLSL,max−2FLSL,min−SL,min2W∗+SL,min2WL2(SL,max−SL,min).

Inserting equation (17) results in:

(22) F∗=2FL(SL,max2+(SL,min−SL,max)(SR,max+SR,min))+SL,min2(−2FR−(SR,max+SR,min)(WL−WR))2(SL,max−SL,min)(SL,max+SL,min−SR,max−SR,min).

We can treat the left transonic case in the same way by integrating over A’B’C’D’, where the integral over section D’C’ is:

(23) ∫0ΔtSR,maxW(x,Δt)dx=W∗+WRSR,max−SR,minSR,max22Δt−SR,maxSR,minWRSR,max−SR,minΔt.

The integral over the quadrilateral gives the inter-cell flux:

(24) F∗=2FRSR,max−2FRSR,min+SR,max2W∗−SR,max2WR2(SR,max−SR,min)=−2FR(SR,min2+(SL,min+SL,max)(SR,max−SR,min))+SR,max2(2FL−(SL,max+SL,min)(WL−WR))2(SR,max−SR,min)(SR,max+SR,min−SL,min−SL,max).

After determining the limiting wave speeds, the inter-cell flux can directly be computed from equations (13), (14), (19), (22) and (24). In the remainder of this publication, we will refer to the proposed scheme as Harten, Lax, van Leer, Einfeldt Shock Unsteadiness (HLLESU).

2.3 Wave speed estimates

Toro et al. (1994) gave multiple possible estimates for the left and right wave speeds. The estimates that are implemented in the applied code follow Einfeldt (1988). Based on the velocity u and speed of sound a, we use:

(25) SL=min(uL−aL,u∗−a∗)

(26) SR=max(uR+aR,u∗+a∗),

With a^* = max (a_L, a_R) and:

(27) u∗=uLρL+uRρRρL+ρR.

The apparent thickness of the shock is generally a function of the mean flow upstream and the intensity of the turbulent fluctuations. Lacombe et al. derived a relation for the turbulent shock thickness from the Rankine–Hugoniot conditions (Lacombe et al., 2021), a simpler estimate relates the thickness to the magnitude of the turbulent velocity fluctuations ΔS∝u″∝k. We may therefore give estimates for the left and right wave speed variance as:

(28) ΔSL=c′u″u″¯L=cS23kL

(29) ΔSR=c′u″u″¯R=cS23kR.

Note that the turbulence is assumed to be isotropic, where u″u″¯:=u1″u1″¯=u2″u2″¯=u3″u3″¯ and k=32u″u″¯. The parameter c_s has been added to tune the scheme.

Two main options exist for applying the wave speed variance to the left and right mean wave speeds. We can either choose the outer limits of the fans as the left and right wave speeds:

(30) SL,min=SL, SL,max=SL+ΔSLSR,max=SR, SR,min=SR−ΔSR,

or the inner limits:

(31) SL,min=SL−ΔSL, SL,max=SLSR,max=SR+ΔSR, SR,min=SR,

See Figure 4 for a schematic representation. When using equation (30), there may be cases with overlapping fans. If we want to treat these correctly, the integration becomes more complicated. For this publication, the overlapping cases were neglected. When equation (31) is used, this problem is avoided because overlapping fans are impossible. There is also a third, centered option:

(32) Sj,min=Sj−ΔSj2

(33) Sj,max=Sj+ΔSj2.

This form reduces exactly to the HLLE scheme and is therefore not investigated further.

2.4 Numerical dissipation

To reduce mean flow gradients, the proposed scheme must add a component to the flux vector, which has the type:

(34) Fadd∝−∂W∂x.

Three terms can be distinguished in the enumerator of the subsonic flux equation (19). The first term 2F_R(S_L,max + S_L,min) only depends on the left wave speeds and the second term −2F_L(S_R,max + S_R,min) only depends on the right wave speeds. The third term (−S_R,max − S_R,min) (−S_L,max − S_L,min)(W_L − W_R) is proportional to the product of the left and right wave speeds. Hence, the magnitude of the third term increases when the limiting left and right wave speeds S_L,min and S_R,max are chosen as described in equation (31). In the subsonic case, S_L,min < S_L,max < 0 and the additional flux is:

(35) Fadd>0 if WL>WR and

(36) Fadd<0 if WL<WR,

which satisfies equation (34).

The enumerator in the transonic cases equations (24) and (22) has a similar structure. There are two terms proportional to the left and right flux and the square of the left or right wave speeds. The third term is proportional to (W_L − W_R) and the third power of left/right wave speeds. Again, this term increases in magnitude compared with the other two when equation (31) is used. In the right transonic case, the last term is:

(37) −SL,min2(SR,max+SR,min)2(SL,max−SL,min)(SL,max+SL,min−SR,max−SR,min)︸A(WL−WR).

Since S_L,min < 0, S_L,max > 0 and S_L,max < S_R,min, A > 0. Equation (34) is therefore fulfilled. In the left transonic case, the last term is:

(38) −SR,max2(SL,max+SL,min)2(SR,max−SR,min)(SR,max+SR,min−SL,min−SL,max)︸A(WL−WR).

Again, since S_R,min < 0, S_R,max > 0 and S_R,min > S_L,max, A > 0 and equation (34) is fulfilled.

This shows that by using equation (31) for the wave speed limits, an additional diffusive flux is introduced. Using equation (30), the additional flux component may have either sign, depending on the magnitude of the wave speeds. This results in unpredictable behavior and can even introduce an anti-diffusive flux component.

3. Numerical implementation

The scheme was implemented in Navier–Stokes multi-block (Leyland et al., 1995; Hoarau et al., 2016); Goebel et al., 2012), which is a cell-centered finite volume solver that uses block-structured meshes with hexagonal cells. It solves the compressible RANS equations, which, using Einstein summation, can be written as (Wilcox, 2006):

(39) ∂ρ¯∂t+∂∂xi(ρ¯u˜i)=0,

(40) ∂∂t(ρ¯u˜i)+∂∂xj(ρ¯u˜ju˜i+δijp)=∂∂xj[t¯ji+ρ¯τji],

(41) ∂∂t(ρ¯E˜)+∂∂xj(ρ¯u˜jE˜+u˜jp¯)=∂∂xj[−qj−qT,j+tjiui″¯−ρuj″1/2ui″ui″¯]+∂∂xj[u˜i(t¯ij+ρ¯τji)],

p is the mean pressure, E˜=e˜+1/2u˜iu˜i+k is the mean total energy, t¯ij=μ[∂u˜i∂xj+∂u˜j∂xi−23∂u˜k∂xkδij] is the mean laminar shear stress tensor and ρ¯τij=−ρui″uj″¯ denotes the Reynolds stress tensor. The components of the laminar and turbulent heat flux q_i and q_i,T are proportional to the gradient of the mean temperature T˜: q(T,)j=−λ(T)∂T˜/∂xj. An eddy diffusivity model with constant turbulent Prandtl number Pr_T = (c_pµ_T)/λ_T = 0.9 was used. Note that the second terms on the left-hand side of equations (39)–(41) represent the inviscid fluxes that are treated with the proposed scheme.

The system of equations is closed with the thermodynamic equation of state:

(42) p¯ρ¯=RT˜

and caloric equation of state:

(43) e˜=cvT˜,

where R is the specific gas constant. With the ratio of specific heats γ, the specific heat capacity at constant volume c_v = R/(γ − 1) can be defined.

To keep the influence of the turbulence model consistent between schemes, the Wilcox k-ω 2-equation model (Wilcox, 2006) without compressibility correction was used for all tests. This model has previously shown good performance in attached flows, even for hypersonic flows (Buck and Mundt, 2021); its equations for the turbulence kinetic energy and specific dissipation can be written as (Wilcox, 2006):

(44) ∂ρ¯k∂t+u˜j∂ρ¯k∂xj=ρ¯τij∂u˜j∂xj−β∗ρ¯kω+∂∂xj[(μ+σ∗μT)∂k∂xj],

(45) ∂ρ¯ω∂t+u˜j∂ρ¯ω∂xj=γωkτij∂u˜j∂xj−βρ¯ω2+∂∂xj[(μ+σμT)∂ω∂xj].

The eddy viscosity can then be computed as μT=α∗(ρ¯k)/ω. The model parameters β^*, σ^*, γ, β, σ and α^* can be found in Wilcox (2006).

An implicit scheme was used to generate steady-state solutions, where the system of equations is solved with an lower-upper symmetric Gauss–Seidel method each timestep. Unsteady simulations were integrated in time with an explicit five-stage Runge–Kutta scheme. Before the Riemann solver is applied to compute the inter-cell fluxes, the primitive variables are extrapolated to the left and right side of the cell interface using MUSCL extrapolation (Monotonic Upstream-centered Scheme for Conservation Laws) with the van Leer limiter (van Leer, 1974), or the van Albada limiter (van Albada et al., 1982), which results in second order accuracy in space for smooth regions.

4. Results

To investigate the performance of the proposed scheme, three test cases were chosen that should examine the scheme in shock/turbulence interaction cases. The simplest STI case is the steady interaction of a statistically planar shock wave with homogenous isotropic turbulence, which can be simulated on an orthogonal mesh. However, this case is numerically difficult because the shock is not geometrically fixed in space, so the outlet conditions must be chosen iteratively to achieve a stationary shock. To remedy this, Sod’s unsteady shock tube problem (Sod, 1978) was chosen. This case is equivalent to the steady planar STI in a nonmoving reference frame. To suppress the dissipation of pre-shock turbulence, sustaining terms (Spalart and Rumsey, 2007) are added to the turbulent governing equations (44) and (45). The standard Sod problem with a shock Mach number of M_shock ≈ 1.66 and one with M_shock = 3 were investigated. An exact solution exists for the inviscid Sod problem.

While the proposed scheme should add numerical dissipation near shocks, the performance in viscous cases should not be impaired. A second test case was therefore chosen to examine the performance in the simplest viscosity-dominated case: a flat plate boundary layer. Supersonic conditions were chosen because of the intended use cases for the scheme. They are also more demanding because of the large velocity gradient near the wall. The M_∞ = 2 and M_∞ = 5 cases by Jespersen et al. (2016) were used.

Finally, the goal of the new scheme is to improve the prediction in realistic STI cases, so the shock/turbulent boundary layer interaction case by Schülein et al. (1996) was chosen. It consists of an oblique shock wave that impinges on a flat plate boundary layer at M_∞ = 5.

4.1 Viscous shock tube problem

Sod’s shock tube test case consists of a 1-m tube that initially has a membrane at x = l/2, see Figure 5. The initial conditions for the left (L) and right (R) half of the tube are given in Table 1; the standard conditions produce a shock Mach number of around 1.66. Upon initialization, the membrane is removed instantaneously resulting in a right-running shock wave and contact discontinuity and a left-running expansion fan. Because the low pressure results in unrealistic flow properties for the standard case, it was run with synthetic constant wave speed variances ΔS, whereas the M_shock = 3 case produces relevant conditions, so ΔS∝k was used.

The tube was discretized with 1,000 cells in axial direction and all walls were modeled with the slip condition. Both cases were run until the inviscid shock reached a similar final position, resulting in t_max = 0.25 s for the standard case and t_max = 0.00044 s for the M_shock = 3 conditions. For the second case, the initial pressure in the right section of the tube was chosen to be the atmospheric pressure, which necessitates the high initial pressure in the left section.

Because the initial velocity is zero in the shock tube, initial turbulence must be added artificially. The initial turbulence kinetic energy was set to kinitial=Tu umax2 with a turbulence intensity of Tu = 1%, where u_max is the maximum velocity at t_max in the inviscid case. The gas was air with R = 287 J/kgK and γ = 1.4. Using turbulence models also requires the inclusion of the laminar viscous terms, so the laminar dynamic viscosity was set to the sea level value of air at T = 0°C: µ = 1.736 × 10⁻⁵ Pas = const for both cases.

The velocity distribution and shock dilatation for Sod’s problem with standard initial conditions are shown in Figure 6. The first option, equation (30), of choosing the wave speeds is referred to as S_outer, the second option, equation (31), as S_inner. As the velocity distribution on the left shows, the proposed scheme does not affect the overall accuracy of the solution. The position of all waves is nearly identical to the exact and HLLE solutions. The shock dilatation on the right shows that choosing the outer definition of equation (30) results in a slight increase in shock speed and shock strength (dashed lines). Furthermore, solutions with ΔS ≥ 0.5 m/s show spurious oscillations and ΔS = 2 m/s results in entirely nonphysical behavior (not shown). The reason might be the overlap of waves in subsonic and transonic cases, which is not accounted for. The inner wave speed definition equation (31), however, does not impair the shock speed and succeeds in reducing the shock dilatation (solid lines).

This shows the impact of the observations in Section 2.4. Using the inner wave speed definition of equation (31) increases the wave speeds and adds numerical dissipation. Mean flow gradients are decreased, resulting in lower shock dilatation. The outer choice of wave speeds in equation (30), on the contrary, introduces a counter-gradient flux that can produce spurious oscillations.

The diffusive component of the flux is:

(46) Fadd=(SR,max+SR,min)(SL,max+SL,min)2(SL,max+SL,min−SR,max−SR,min)(WL−WR)

in the subsonic case and:

(47) Fadd=−SL,min2(SR,max+SR,min)2(SL,max−SL,min)(SL,max+SL,min−SR,max−SR,min)(WL−WR)

in the right transonic case. The three components of the flux vector as a function of ΔS for the standard Sod case at the right-running shock wave are shown in Figure 7 for the inner and outer wave speed definitions. For the inner wave speed definition, the additional flux is positive for all ΔS, resulting in a diffusive property and stable behavior. Using the outer wave speed definition, the additional flux decreases as ΔS is increased, making the simulation less stable. At ΔS = S_R,max − S_L,min ≈ 2.528 m/s, the flux function is undefined. When ΔS > S_R,max − S_L,min, the additional flux component is negative, resulting in counter-gradient flux, which is unstable. To achieve the desired level of numerical dissipation, the inner wave speed definition equation (31) must therefore be used. It was applied for the rest of this publication. Using equation (30) highlights numerical issues near strong gradients and does not appear useful.

Figure 8 shows the velocity distribution and shock dilatation for the M_shock = 3 case. Compared with the standard conditions, the turbulent distributions differ more from the exact inviscid results. The higher turbulence kinetic energy appears to increase the impact of the turbulent terms in the momentum equation and energy equation, which results in a faster and stronger shock wave in the turbulent case. The post-shock velocity in the turbulent case is therefore higher than the exact solution. Both the standard HLLE and the proposed scheme give similar results.

The shock dilatation in Figure 8 also shows that the shock position moves further to the right with the proposed scheme, an effect that is amplified when c_s is increased. The graph also shows a reduced dilatation for c_s = 100; the stronger dilatation with c_s = 10 might be a result of the position of the shock wave in the c_s = 1 case.

To investigate the impact of the lower shock dilatation on the evolution of turbulence, the TKE amplification across the shock wave is shown in Figure 9. Figure 9(a) shows the data for the standard conditions over the wave speed variance ΔS. As the wave speed increases, the turbulence amplification decreases. The distribution appears close to linear, with beginning saturation at ΔS > 5 m/s. Figure 9(b) shows the same data for the M = 3 case over different values of the scaling factor c_s, where ΔS was treated as a function of the turbulence kinetic energy. Similar to the standard case, a larger wave speed variance decreases the turbulence amplification across the shock. Interestingly, the shock-induced turbulence amplification is much smaller in the M = 3 case compared to the standard conditions. These results indicate that the proposed scheme can reduce the steep velocity gradients across shock waves and limit the resulting turbulence production.

4.2 Supersonic flat plate boundary layer

The previous results show that the proposed scheme can reduce the impact of shock waves on turbulence. The supersonic flat plate boundary layer case is used to determine whether the added numerical dissipation affects the performance in cases that are not dominated by shock/turbulence interaction. A supersonic flat plate with a length of 2 m was chosen, see Figure 10. The boundary conditions from Jespersen et al. (2016) are given in Table 2. The chosen wall temperatures are approximately equal to the respective recovery temperatures.

Based on recommendations from the NASA Langley turbulence modeling resources, the plate was discretized with 544 × 384 cells, with 448 cells along the solid wall and 96 cells in the freestream. The first grid point has a wall distance of 5 × 10⁻⁴mm, which results in a maximum nondimensional wall distance of y⁺ ≈ 0.1 near the leading edge of the plate. The fluid was air with γ = 1.4 and R = 287 J/kg K. Sutherland’s model was used for the viscosity, and the conductivity was determined with a constant Prandtl number Pr = 0.72.

The top row of Figure 11 shows the skin friction cf=τw/(1/2ρ∞U∞2) over the momentum thickness Reynolds number Re_Θ = (ρ_∞U_∞Θ)/µ_∞, where Θ is the momentum thickness:

(48) Θ=∫0δ99.5ρρ∞uU∞(1−uU∞)dy.

Note that the 99.5% velocity boundary layer edge δ_99.5, rather than the maximum extent of the domain, was chosen as the upper limit of the integral. Because the flow state downstream of the shock differs slightly from the freestream conditions, the numerical integration would otherwise produce large errors (Jespersen et al. (2016). The theoretical compressible skin friction distribution is obtained from the incompressible Karman-Schoenherr relation (Huang et al., 1993), to which the van Driest II transformation (Gatski, 2013) is applied.

For M_∞ = 2 all skin friction profiles are similar, only the one computed with HLLESU and c_s = 100 is slightly increased. The skin friction predicted by theory is lower than the numerical data. It seems that the added numerical dissipation has a small impact on the boundary layer. The M = 5 case also shows nearly identical skin friction for the numerical results. Again, the theory predicts lower skin friction.

Finally, the bottom row of Figure 11 shows the nondimensional boundary layer profiles at Re_Θ = 10,000 compared to theoretical data from Jespersen et al. (2016). The nondimensional velocity u⁺ = u/u_τ is shown over the nondimensional wall distance y⁺ = (yu_τ)/v, where uτ=τw/ρ. The simulation results for both cases show a favorable performance of all applied schemes. The agreement to theory is satisfactory. These results indicate that the added numerical dissipation in the proposed scheme does not impair the performance in attached, shock-free boundary layers.

4.3 Shock/turbulent boundary layer interaction

To investigate the performance of the proposed scheme in a shock/turbulence boundary layer interaction (STBLI) case, the experiment by Schülein et al. (1996) was simulated. This case has proven challenging for RANS models in previous studies (Roy et al., 2018), which failed to reproduce the interaction of the shock with the incoming turbulence and consequently over-predicted the wall heating. Because the wall heat flux is significantly lower in laminar simulations of this case, it is reasonable to assume that excessive turbulence levels are the root cause of the problem. The proposed scheme should reduce turbulence amplification across the shock, resulting in lower post-shock turbulence levels and lower wall heat flux.

Figure 12 shows the numerical domain, the operating conditions are given in Table 3. The fluid is air with γ = 1.4 and R = 287 J/kg K, the resulting unit Reynolds number is Re_∞,1 = 40 × 10⁶ 1/m. A 10^° shock generator creates an oblique shock wave that interacts with a flat plate boundary layer. The shock generator is positioned, so that the shock impinges on the bottom wall 0.35 m downstream of the leading edge for an inviscid flow. The rarefaction fan generated by the end of the shock generator impinges on the lower wall downstream of x = 0.44 after interacting with the reflected shock. Experimental heat flux and shear stress data were measured by Schülein et al. (1996) and Schülein (2006).

Other studies of this case, e.g. Wallin and Johansson (2000), often limit the computational domain to a rectangular region where the effect of the incident shock and the rarefaction fan are represented by defining boundary profiles. In this study, the entire plate and the shock generator were simulated, which removes the need for precursor simulations and accurate interpolation of boundary profiles.

Five grids were compared in a grid convergence study, their properties are shown in Table 4. Grids 1, 2 and 3, all have the same cell distribution along the channel with varying cell counts in wall normal direction. Grids 4 and 5 are refined in the interaction region (IR) (0.34 m < x < 0.37m). This also required a minor change in the number of axial nodes outside the IR. Grid 5 was created by refining grid 4 by a factor of 2 in x and y direction. The first node in the IR has a wall distance of 0.77 mm for grids 1 to 4, which results in a nondimensional wall distance of ymax+<0.3 in the IR. The lower wall distance of grid 5 results in ymax+<0.15. Because the exact distance between the shock generator and the lower wall were not specified in the original paper by Schülein et al. (1996), the shock generator was modeled as a slip wall with some minor clustering to better resolve the attached shock wave.

The grid convergence study was conducted with the HLLE scheme. Unlike for the other cases, the modified van Albada limiter (van Albada et al., 1982) was used, which does not clip the solution near smooth extrema and alleviates known convergence problems in regions with near-uniform flow.

Figure 13 shows the mesh convergence of the wall heat flux and skin friction in the IR. The Stanton number St = q/(ρ_∞U_∞c_p(T₀ − T_w)) reported in Schülein (2006) was made dimensional with the data given in Table 3.

As the grid resolution increases, the peak heat flux and skin friction decrease. The solutions with grid 4 and grid 5 are very similar even though grid 5 doubles the number of cells. Therefore, the grid convergence with grid 4 seems satisfactory and grid 4 will be used for the rest of the simulations.

Figure 14(a) shows the heat flux distribution along the wall for the different schemes compared to experimental data by Schülein (2006). Upstream of the shock interaction, the agreement to the experiment is good. The experimental data show a more upstream position of the separation shock in both the heat flux and skin friction distributions. Downstream of the reattachment, the numerical results strongly overpredict the wall heat flux, whereas the skin friction data lie within the experimental margin of error, see Figure 14(b).

Comparing the HLLESU scheme with different values of c_s shows an increase in wall heat flux for increasing c_s. While c_s = 1 and c_s = 10 are overall very similar to the HLLE result, c_s = 100 and especially c_s = 1,000 overpredict the heat flux even more strongly than the HLLE result. Similarly, increasing c_s also increases the skin friction, but the difference between schemes is smaller and the c_s = 1,000 case even improves the prediction at x > 0.4 m.

Investigating the numerical Schlieren images in Figure 15 shows that the HLLESU scheme indeed adds numerical dissipation resulting in wider shock waves. This effect is especially noticeably at c_s = 1,000, but can also be observed at lower c_s, where the oscillations around the reflected shock are reduced. The maximum value of the turbulence kinetic energy and its location are also given in Figure 15. Interestingly, the maximum TKE does not change drastically for c_s = 1 and c_s = 100 and even decreases appreciably for c_s = 1,000. This indicates that the scheme successfully reduces the mean flow gradients and turbulence amplification across shocks also for STBLI flows given appropriate scaling of c_s.

The increased wall heat flux can therefore not be explained by excessive turbulence amplification alone and appears to be a result of the added numerical diffusion of the scheme. A first-order solution with the HLLE scheme, which also adds numerical dissipation, similarly over-predicts the heat flux significantly, albeit at even higher levels. These diffusive terms effectively increase the viscous momentum and energy transfer, which results in higher wall heat flux and wall shear stress. The fact that this effect is not visible for the attached flat plate boundary layer is likely a result of the much higher turbulence intensity Tu_max ≈ 7,700 in the IR of the STBLI case compared to Tu_max ≈ 0.075 in the flat plate case. These results seem to indicate that a fixed value of parameter c_s may be inadequate. Future studies should also consider a variable c_s that could, e.g. be a function of the local pressure gradient, reducing the added dissipation in smooth regions.

It should also be noted that the axial Mach number M_x = u_x/a is greater than one across the incident shock. Because the supersonic fluxes are not changed from the HLLE scheme, no numerical dissipation is added in x-direction. In y-direction, the flow is subsonic in the entire domain and the HLLESU flux function are active. Numerical dissipation is therefore only added in wall-normal direction, which appears to affect the boundary layer profiles significantly.

5. Conclusion and outlook

The aim of this publication was to show the derivation and performance of an extended Riemann solver for RANS simulations. The concept is based on the observation that time-averaged solutions of instantaneously discontinuous waves in turbulent flows are continuous in space. With that, we can derive equations for the inter-cell flux based on the left and right mean flow and the velocity variance computed from the turbulence kinetic energy. A mathematical analysis of the resulting flux functions shows that the new scheme adds numerical dissipation in the subsonic and transonic cases.

The shock tube test case indicates that the proposed extension of the HLLE scheme decreases turbulence amplification in shock turbulence interaction given an appropriate choice of the limiting wave speeds. The results show that the wave speed variance must be applied to the outside of the wave pattern to introduce an additional dissipative flux. Based on the observation that RANS simulations overpredict turbulence amplification across shocks, the results generated with the proposed HLLESU scheme, which show lower postshock turbulence levels and are an improvement over the ones generated with the HLLE scheme. It is also revealed that the scheme must be calibrated thoroughly to find ideal relations for the scaling parameter c_s.

The flat plate test case gives very similar results for the standard HLLE scheme and the proposed extension. This demonstrates that the added numerical dissipation does not impair the results in cases without strong shock waves. Recalibration of existing models should not be necessary.

The shock/turbulent boundary layer interaction case shows that the scheme may add numerical dissipation excessively in cases with locally high turbulence intensity. While the more dissipative numerical scheme reduces the velocity gradients and therefore the postshock turbulence, the required numerical dissipation impairs the results further downstream. Because of that the wall heat flux and wall shear stress are increased. Switching to a less dissipative scheme in regions without strong gradients could be a possible solution to the problem. This might reduce the shock-induced turbulence amplification without deteriorating the boundary layer. Additional numerical experiments have to be performed to investigate the performance of the proposed scheme in different test cases.

While this study succeeded in proving the effectiveness of the proposed scheme in reducing mean flow gradients, it also shows that the performance of the scheme depends strongly on the free parameter c_s. Choosing an appropriate value for this parameter introduces additional complexity and additional work is required to determine an appropriate value. It may be difficult to find a single value that is satisfactory for all cases, so c_s ≠ const may also be investigated. A practical implementation of this could scale the parameter c_s by the magnitude of some mean flow gradient, e.g. the pressure gradient |∇p|, which results in c_s = 0 in smooth regions, recovering the HLLE scheme. As the different test cases seem to show different sensitivities to the value of c_s, a different definition of the wave speed variance with ΔS(k)∝k could also be investigated.

To further improve the performance of the scheme, the contact wave could be included in the wave pattern, extending the method to a complete Riemann solver, such as HLLC (Toro et al., 1994). If the idea of combining turbulence modeling and numerical treatment should be applied to Riemann-solver-free schemes, additional artificial viscosity terms, based on some turbulent property, e.g. the TKE, could be included in other numerical methods as well.

Figures

Figure 1.

Inter-cell flux in cell-centered finite-volume schemes

Figure 2.

Wave diagram for the derivation of Riemann solvers for inviscid and mean turbulent flows

Figure 3.

Transonic cases for turbulent flows: the width of the transonic fans is increased to improve readability

Figure 4.

Selection of the possible estimates for the minimum/maximum wave speed based on the wave speed variance ΔS

Figure 5.

Sod’s shock tube problem

Figure 6.

Velocity distribution and dilatation across the shock for Sod’s problem at standard conditions at t_max

Figure 7.

Diffusive component of the flux vector F_add at the right-running shock wave in the standard Sod case, units on the vertical axis are omitted for clarity

Figure 8.

Velocity distribution and shock dilatation for Sod’s problem at M_shock = 3 conditions at t_max, wave speeds with the inner limits from equation (31)

Figure 9.

Turbulence amplification across the shock for Sod’s problem, equation (31) was used for the wave speeds

Figure 10.

Domain of supersonic flat plate cases

Figure 11.

Supersonic flat plate boundary layer flow

Figure 12.

Domain of the STBLI case by Schülein (2006)

Figure 13.

Mesh convergence of the STBLI case

Figure 14.

Wall heat flux and shear stress for the shock/boundary layer interaction case with different schemes

Figure 15.

Numerical Schlieren images of the interaction region in the STBLI case: ∂ρ/∂x is shown

Table 1.

Initial conditions for the shock tube test cases

Case	pL [Pa]	pR [Pa]	pL [kg/m3]	pR [kg/m3]	uL [m/s]	uR [m/s]
Standard	1	0.1	1	0.125	0	0
Mshock = 3	64 100 968	101 325	817.676	1.293	0	0

Source: Created by the authors

Table 2.

Boundary conditions for the supersonic flat plate boundary layer

Case	M∞ [−]	T∞ [K]	Tw [K]	Re∞,1 [1/m]	Tu∞ [%]	µT/µ [−]
1	2	300	513.6	15 × 106	0.004	0.009
2	5		1,635		0.002

Source: Created by the authors

Table 3.

Environmental conditions for the STBLI case by Schülein et al. (1996)

M∞ [−]	p0 [Pa]	T0 [K]	Tw [K]
5	2.12 × 106	410	300

Source: Created by the authors

Table 4.

Mesh properties for the grid convergence study

Grid	Nodes	IR cells	Description
1	100 × 435	24	Baseline mesh
2	150 × 435	24	Fine wall normal
3	200 × 435	24	Finer wall normal
4	200 × 696	300	Refined in IR
5	399 × 1,389	600	Refinement level 2

Source: Created by the authors