Non-Zonal Detached Eddy Simulation coupled with a steady RANS solver in the wall region

Xiao et al. [1] proposed an interesting hybrid LES/RANS method in which they use two solvers and solve the RANS and LES equations in the entire computational domain. In the present work this method is simpliﬁed and used as a hybrid RANS-LES method, a wall-modeled LES. The two solvers are employed in the entire domain. Near the walls, the ﬂow is governed by the steady RANS solver; drift terms are added to the DES equations to ensure that the time-integrated DES ﬁelds agree with the steady RANS ﬁeld. Away from the walls, the ﬂow is governed by the DES solver; in this region, the RANS ﬁeld is set to the time-integrated LES ﬁeld. The disadvantage of traditional DES models is that the RANS models in the near-wall region – which originally were developed and tuned for steady RANS – are used as URANS models where a large part of the turbulence is resolved. In the present method – where steady RANS is used in the near-wall region – the RANS turbulence models are used in a context for which they were developed. In this method, it may be worth while to use an accurate, advanced RANS model. The EARSM model is used in the steady RANS solver in the present work. The new method is called N-Z S-DES .


Introduction
DES (Detached-Eddy Simulation) uses unsteady RANS near walls (URANS region) and LES further away from walls (LES region). The resolved turbulence in the URANS region is often larger than the modeled part. But the RANS models used in the URANS region were originally developed and tuned in steady RANS simulations. Hence the accuracy and the validity of the RANS models in the URANS region can be questioned. In the present work, DES is coupled with steady RANS near the walls. We denote the method N-Z S-DES (Non-Zonal approach using Steady RANS coupled to DES).
Xiao et al. [1] proposed a new method in which they solve both the LES and RANS equations in the entire domain. The flow is in the near-wall region governed by the RANS equations and in the outer region it is governed by the LES equations. This is achieved by adding drift terms in the LES and RANS equations. In the interface region(s), the drift terms are modified by a linear ramp function. Drift terms are used in all equations in the RANS equations (momentum equations, the pressure equation (PISO is used)) and in the modelled turbulent equations (k and ε). Two drift terms are added in the LES momentum equations; one to ensure than the mean velocity fields in the RANS and LES equations are the  same and one to ensure that the total turbulent kinetic energies are the same.
In [2] they extended the method to account for nonconformal meshes. They used a Cartesian mesh for the LES equations and a body-fitted mesh for the RANS equations. They applied the method to fully developed channel flow and the flow over periodic hills. The paper shows how an accurate academic solver -massively parallel -can be combined with industrial, flexible RANS solvers.
Tunstall et al. [3] implement and use the method in [1] and modify it (different ramp function, different constants, reducing the number of case-specific constants etc). They apply it to fully developed channel flow and a rather complex flow consisting of a pipe junction including heat transfer. Hence, they have to introduce drift terms also in the energy equations.
Laage de Meux et al. [4] use forcing to achieve resolved Reynolds stress profiles equal to target modeled RANS stress profiles. The resolved Reynolds stresses are integrated in time in the same way as in [1,3].
Breuer and Schmidt [5] use an advanced RANS turbulence model -the Explicit Algebraic Reynolds Stress Model, EARSM -in an hybrid LES-RANS methodology. However, as mentioned above, the disadvantage is that the RANS equations are solved in transient mode, where a large part of the large-scale turbulence is resolved.
Davidson [6] use the same model as in the present study. In [6], however, the interface is defined along a pre-selected gridline and the grid in the hump flow simulations is much coarser upstream the hump and in the outlet region. Moreover, the present work uses a different timescale in the drift term as well as evaluates different locations of the interface (see Eq. 6).
In the present study, the steady RANS equations are solved. Here it makes sense to use advanced RANS turbulence models, since these models were developed for steady RANS. The EARSM [7] is used in the RANS solver. The present method is in many aspects similar to that proposed in [1,3] but it is simplified: the RANS equations are used in steady mode, a more advanced RANS turbulence model is used and the present method includes fewer drift terms and tuning constants.

Numerical solvers
The momentum equations with an added turbulent viscosity read where the first term on the right side is the driving pressure gradient in the streamwise direction, which is used in the fully-developed channel flow simulations.

DES solver
An incompressible, finite volume code is used [8,9]. The convective terms in the momentum equations are discretized using central differencing.
Hybrid central/upwind is used for the k and ω equations. The Crank-Nicolson scheme is used for time discretization of all equations. The numerical procedure is based on an implicit, fractional step technique with a multigrid pressure Poisson solver [10] and a non-staggered grid arrangement.

RANS solver
An incompressible, finite volume code -CALC-BFC [11] -is used. The transient term in Eq. 1 (the first term on the left side) is not included. The convective terms in the momentum equations are discretized using the second-order bounded scheme, MUSCL [12]. Hybrid central/upwind is used for the k and ω equations. The numerical procedure is based on SIMPLEC and a staggered grid arrangement.

The N-Z S-DES methodology
Two sets of equations are solved (steady RANS solver, see where τ r = max(0.1k/ε, ∆t) following [3]. No drift term is used in the pressure equation. · T indicates integration over time, T , i.e.
where a = 1/(1+∆t/T ) and n denotes the timestep number. Note that although the flow cases in the present work include homogeneous direction(s), no space averaging is made in Eq. 3. It may be noted that although the velocity field in the RANS solver is steady, it is time integrated when used in Eqs. 2 and 4 because it varies slightly in time. This time integration may not be necessary.
In the LES region, the RANS velocities are prescribed as v RAN S i = v LES i T by adding a large source term, i.e.
where ǫ = 10 −10 . The pressure is simply set asp RAN S = p LES T and the pressure correction is set to zero. This means that, in reality, the steady RANS solver needs to be solved only in the wall region. In the LES region the momentum equations in the RANS solver are merely transporting the turbulence quantities, k and ω, to ensure that correct values of k and ω are transported into the RANS region through the LES-RANS interface at y = δ I , see Fig. 1. The pressure, p LES j+1 T , and the streamwise velocity, ū LES j+1 T , at the LES-RANS interface are used as a boundary condition for the RANS equations in the wall region, see Fig. 3. The wall-normal velocity,v RAN S j , is solved for using the pressure at node j + 1.

The k − ω model
The Wilcox k − ω turbulence model reads where d/dt = ∂/∂t +v j ∂/∂x j (∂/∂t = 0 in the RANS solver). The standard coefficients are used, i.e. C ω1 = 5/9, C ω2 = 3/40, σ k = σ ω = 2 and C µ = 0.09. The location of the interface is defined as C I = 1.0 is the standard DES value. Here we use two different values, C I = 1.0 and C I = 1.4. The result of the latter value is that the interface location is moved further away from the wall compared with standard DES.

The k − ω model in the DES solver
The DES equations are solved in the entire region, but they govern the flow only in the LES region, see Fig. 1. In the RANS regions, the lengthscale in Eq. 5 is computed as ℓ t = k 1/2 /(C µ ω) and in the LES region it is taken from the standard DES model, i.e.
ℓ t = C DES ∆ max , ∆ max = max{∆x, ∆y, ∆z} (7) with C DES = 0.67. The location of the interface in the DES solver is defined in the same was as between the RANS solver and the DES solver, i.e. by Eq. 6.

The k − ω EARSM model in the RANS solver
The steady RANS equations are solved in the entire region, but they govern the flow only in the RANS region, see Fig. 1. The Reynolds stresses, v ′ i v ′ j , are computed from the two-dimensional explicit algebraic Reynolds stress model (EARSM) [7].
Vertical black lines show DES interface.

Initialization
The simulations are initialized as follows: first the 2D RANS equations are solved. Anisotropic synthetic fluctuations, (V ′ i ) m , are then superimposed to the 2D RANS field which gives the initial LES velocity field. The initial time integrated fields, v LES i T and v RAN S i T , are also set from the 2D RANS field.
In order to compute (V ′ i ) m , synthetic fluctuations, v ′ i,synt , are computed plane-by-plane (y − z) in the same way as prescribing inlet boundary conditions. The synthetic fluctuations in the y − z planes are coupled with an asymmetric space filter where m denotes the index of the x 1 location and a = exp(−∆x 1 /L int ) and ∆x 1 and L int denote the grid size and the integral length scale, respectively (L int = 0.2).

Results
The first test case is fully developed channel flow with periodic boundary conditions in streamwise (x) and spanwise (z) directions. The Reynolds number, Re τ = u τ h/ν, is 8 000 where h denotes half-channel width. The size of the domain is x max = 3.2, y max = 2 and z max = 1.6. The mesh has 32 × 96 × 32 (x, y, z) cells.
The interface between the wall region and LES is defined by Eq. 6. Four different integration times, T (see Eq. 3), are evaluated. It is important that the sampling time is much larger than the integration time. If it is too small, it often gives an asymmetric time-averaged flow field. : DES solver, modeled turbulence; : DES solver, resolved turbulence; : RANS solver, modeled turbulence.      Figure 4 compares the velocity profiles obtained with N-Z S-DES with that of pure DES. The velocity profiles predicted with N-Z S-DES forT = 1, 10 and 20 agree all very well with Reichardt's law butT = 50 shows a small log-layer mismatch. Pure DES exhibits the usual log-layer mismatch (larger than that forT = 50).
The turbulent viscosity and turbulent kinetic energies are presented in Figure 5. Figure 5a confirms that steady RANS gives a much larger turbulent viscosity than does the DES, i.e. steady RANS give much larger turbulent viscosity than URANS. The reason is that in DES, low k is transported into the URANS region from the LES region (the magnitudes of ω at the interface are similar in steady RANS and DES). It can, however, be noted that close to the wall (for y < 0.01, see zoom), steady RANS and URANS give almost identical modeled kinetic energy, see Fig. 5b. Figure 6a shows the resolved shear stresses. It can be seen that the peak of − ū ′v′ decreases for large (T = 50) and small (T = 1) integration time. In fully developed channel flow, the total (resolved plus modeled) shear stress is given by Normally, τ w = 1 because the driving pressure gradient is equal to one (the first term on the right side of Eq. 1). However, for small and large integration times, it turns out that the drift term in the DES momentum equation increases (see Fig. 6b) because v RAN S i T is slightly larger than v DES i T . Hence, in this case the wall shear stress, τ w , balances not only the driving pressure gradient but also the drift term. As a result, τ w increases. Xiao et al. [2] also report that the drift term affects the resolved shear stresses.
The second test case is the flow over a two-dimensional hump, see PSfrag replacements is refined upstream of the hump and in the outlet region, see Fig. 7a. The spanwise extent of the domain is set to 0.3 so that ∆z = 0.3/64 = 0.0047. The required resolution for an LES away from the wall (in the log-region) is ∆x/δ ≃ 10, ∆z/δ ≃ 20. The streamwise spacing, ∆x, near the wall is shown in Fig. 7b. The inlet boundary layer thickness is δ in = 0.08 which means that the resolution in the inlet region reasonable (∆x/δ in ≃ 10, ∆z/δ in ≃ 17 ). The boundary layer thickness after the recirculation bubble (x > 1.3) is δ ≃ 0.13 (see Fig. 9d) so that ∆x/δ ≃ 3 which is somewhat too small; this region is, however, believed to be less critical than the inlet and hump region.
The inlet profiles are taken from a separate 2D RANS simulation with the same momentum thickness as the experimental velocity profiles. Anisotropic synthetic fluctuations are superimposed to the inlet velocity profile (for more detail, see [14]). Periodic boundary conditions are used in the spanwise direction (z). The interface be- : RANS solver; : DES solver, modeled; : DES solver, resolved.
tween the wall region and the LES domain as well at that between the steady RANS solver and the DES solver is defined by Eq. 6. The pressure coefficient and skin friction are presented in Fig. 8. As can be seen, the agreement with experiments is good for both N-Z S-DES and pure DES except that the pure DES predicts slightly too low a skinfriction upstream of the hump. For the channel flow it was found that for too large and too smallT , the drift term in the DES equation was comparable to the driving pressure gradient. Figure 8a shows no such problems for the hump flow withT = 10; the pressure from the RANS and DES solvers are virtually identical (there are small differences at x ≃ 0 and x ≃ 0.66). Figure 9 compares the predicted velocity profiles with experiments. The S-Z S-DES gives a better agreement with experiments than pure DES; in particular, the latter gives a worse velocity profile of the boundary layer at x = 0.65. The velocity profiles of the RANS solver are also included. The RANS solver velocity profiles in the wall region match those of the DES solver as they should. The locations of the interface is shown with a red plus sign. It may be noted that the location of the interface at x = 1.3 is not seen and it seems that it is located at PSfrag replacements y > 0.2. This will be discussed below.
The turbulent viscosity is shown in Fig. 10. As for the channel flow, the turbulent viscosity predicted by the RANS solver is much larger than that predicted by the DES solver, especially in the recirculation region. Again, as was seen in Fig. 9d, is seems that the interface is located at y > 0.2; this is also seen for the turbulent shear stress in Fig. 11d. Figure 11 presents the predicted and measured shear stresses. The agreement between the resolved shear stresses and the measured is good downstream of the separation point. What is somewhat surprising is that at x = 0.65 the N-Z S-DES models most of the turbulence rather than resolving it (the magnitude of the modeled shear is much larger than the resolved one). This should probably be viewed as a drawback of the N-Z S-DES. : integration times,T (h denotes hump height, see Fig. 2). This is much smaller than for the channel flow. It is believed that the channel flow is particularly sensitive because of the streamwise periodic boundary conditions which lacks a stabilizing inlet boundary condition. Figures 12 -13 show that the flow is only weakly dependent to the integration time,T . The lowest value,T = 5, gives somewhat worse agreement than the other three. The switch between RANS and DES (i.e. the interface) is in Figs. 8 -11 defined by C I = 1.4 (see Eq. 6). In Figs. 12 -13 the switch is defined by C I = 1.0, i.e. the standard switch in pure DES. Figure 14 presents the location of the interface for the hump simulations. When the switch is defined as C I = 1.4 it is seen that there are three regions where the interface is located at y > 0.04: near the inlet (x < −1), at x ≃ 0.12 and at 1.26 < x < 1.8. The interface in these three regions is actually located at the upper boundary, i.e. the entire region is covered by the RANS solver. However, as seen in Figs. 12 -13, the influence of the location is not large (compare the blue solid lines in Figs. 12 -13 with Figs. 8 -9). The largest differences are found in the skinfriction at x ≃ 0.18 and as a result the velocity profile is less full in Fig. 13a than in Fig. 9a.

Conclusions
A new non-zonal model based on a steady RANS solver in the wall region coupled to a DES solver which covers the entire region. The steady RANS solver is called very 10th timestep. The RANS solver dictates the flow in the wall region. A drift term in the wall region in the DES solver forces the time-integrated DES flow to match that of the RANS flow. The steady RANS solver velocity field is set to that of the time-integrated DES flow in the offwall region. The only object of the steady RANS field in the off-wall region is to transport the modeled turbulent quantities (k and ω) into the wall region. They -k and ω -could instead be transported by the time-integrated DES flow.
The new model is evaluated in fully developed channel flow and the hump flow. Very good agreement with experiments is obtained. However, three issues are identified. First, the drift term in the DES solver in the wall region -which is the ratio of the difference between the steady RANS velocity field and that of the time-integrated DES to a turbulent relaxing timescaleis non-zero which acts as an additional driving pressure gradient. No such problem is seen in the hump flow. Second, the time integration is chosen over a certain time, T =T h/U b , where values ofT between one and 50 have been evaluated. For the channel flow simulations, too large (T = 50) or too small (T = 1) values act as an additional driving pressure gradient. For the hump flow, the influence ofT is small. Third, in the attached bound-ary layer over the hump, the modeled turbulence in the DES solver is larger than the resolved one. In the near future, the reason for this will be analyzed by making simulations of developing flat-plate boundary layers.