A Hybrid Large Eddy Simulation Algorithm Based on the Implicit Domain Decomposition

Abstract

The resolution of small near-wall eddies encountered in high-Reynolds number flows using large eddy simulation (LES) requires very fine meshes that may be computationally prohibitive. As a result, the use of wall-modeled LES as an alternative is becoming more popular. In this paper, the near-wall domain decomposition (NDD) approach that was originally developed for Reynolds-averaged Navier–Stokes simulations (RANSs) is extended to the hybrid RANS/LES zonal decomposition. The algorithm is implemented in two stages. First, the solution is computed everywhere with LES on a coarse grid using a new non-local slip boundary condition for the instantaneous velocity at the wall. The solution is then recomputed in the near-wall region with RANS. The slip boundary conditions used in the first stage guarantee that the composite solution is smooth at the inner/outer region interface. Another advantage of the model is that the turbulent viscosity in the inner region is computed based on the corresponding RANS velocity. This shows improvement over those hybrid models that have only one velocity field in the whole domain obtained from LES. The model is realized in the open source code OpenFOAM with different approximations of turbulent viscosity and is applied to the planar channel flow at frictional Reynolds numbers of ${\mathrm{Re}}_{\tau }=950$, 2000, and 4200. Mean streamwise velocity and Reynolds stress intensities are predicted reasonably well in comparison to the solutions obtained with unresolved LES and available DNS benchmarks. No additional forcing at the interface is required, while the log–layer mismatch is essentially reduced in all cases.

1. Introduction

Numerical simulation of turbulent flows has been widely used over recent decades thanks to the advancement of computational power and resources. The direct numerical simulation (DNS) for large-Reynolds-number flows is still limited to simple geometries due to the very large number of grid points required to capture all eddies up to the Kolmogorov scale [1]. In DNS, the grid size required to resolve all structures varies as $N\approx O\left({\mathrm{Re}}^{2.64}\right)$ [2]. In the large eddy simulation (LES) of free-streaming flows, the required grid density increases only modestly with the growth of $\mathrm{Re}$ as $N\approx O\left({\mathrm{Re}}^{0.4}\right)$ [3]. Therefore, LES has evolved as an advanced engineering tool in the analysis of turbulent flows.

The so-called resolved LES, i.e., direct resolution of all energetic scales, is a powerful tool for free shear flows such as jets, mixing layers, and wakes. However, its application to wall-bounded flows is still limited due to the presence of fine eddies encountered in the near-wall area, which necessitates the use of fine meshes for LES [4]. In close proximity to walls, the length scale of dynamically important motions decreases, and energetic eddies are not large any more [5]. It is shown in [2] that the grid resolution requirement for a resolved LES in the boundary layer flows is proportional to ${\mathrm{Re}}^{1.85}$, which is almost as costly as DNS. For a moderate frictional Reynolds number of ${\mathrm{Re}}_{\tau }=2000$, Temmerman et al. [6] showed that the unresolved LES produces a large discrepancy of around 40% from the DNS benchmark and that only resolved LES with a grid more than an order of magnitude finer in each direction can produce acceptable results. Such fine meshes are still a significant computational burden for real engineering applications.

In recent years, the use of wall-modelled LES (WMLES) [7] instead of resolved LES has received more attention. WMLES employs near-wall approximations that permit a much coarser resolution to be used close to the walls without a significant loss of accuracy. The idea is to model the near-wall turbulent structures rather than to resolve them. In WMLES, two approaches are widely used: wall-stress and RANS/LES models. It should be noted that the classification and names of the methods in the literature may vary.

The wall-stress models (WSMs) can be classified as equilibrium or non-equilibrium models [8]. The equilibrium WSM approximates the wall-shear stress (WSS) based on the logarithmic law or power law and applies this approximation in the form of a Neumann velocity boundary condition to the wall. The use of equilibrium WSMs dates back to the papers by Deardorff [9] and Schumann [10]. This approximation disregards convection and pressure-gradient terms in the boundary layer and is therefore primarily valid for simple configurations but is inaccurate with separated flows or recirculating zones. The deficiency of wall laws and the importance of streamwise grid resolution in predicting separation zones are already addressed by Temmerman et al. [11].

In the two-layer WSM (non-equilibrium), an auxiliary fine one-dimensional (1D) grid is embedded between an LES grid point and the wall. Then, a simplified set of the Reynolds-averaged Navier–Stokes equations (RANS) based on the thin boundary layer equations (TBLEs) is solved on the embedded mesh [12,13,14,15]. The outer-layer LES provides the velocity boundary condition for the inner-layer RANS model. In return, the inner-layer wall-model predicts the wall shear stress required by the LES at the wall. Non-equilibrium effects, such as unsteady advection and pressure gradient, can be accounted for using the TBLE model. Nonetheless, the coupling between the LES and RANS models is rather weak and only occurs through the boundary conditions at the wall ($y=0$) and the interface boundary ($y=y_{int}$). Cabot and Moin [14] demonstrated the importance of the pressure gradient and advection terms in the simulation of a backward-facing step. Several modifications have been made to the model over the course of recent years (see, e.g., [8,16,17]).

In the other approach, often referred to as hybrid RANS/LES models, RANS and LES solutions are usually (although not always) blended [7,18,19]. For this purpose, a blending function is introduced [20] for the modeled turbulent viscosity [7] or the length scale as in detached eddy simulations (DESs) [21,22,23]. DES has shown very good efficacy for separated flows but not for attached flows [22]. In addition, blending RANS and LES models always has an intermediate region where the blended model is not justified.

In the user-defined-interface hybrid methods, the inner region is covered by RANS and separated from the LES area by an interface. The influence of near-wall turbulent structures is modeled by RANS. In this way, the near-wall RANS region can be superimposed onto the LES domain [6,24]. Then, the boundary conditions (BCs) for RANS are obtained from LES. Alternatively, switching from LES to RANS can be undertaken by selecting the turbulent viscosity ${\nu }_T$ and kinetic energy k from either LES or RANS models [25,26]. A detailed description of BCs at the RANS/LES interface is provided in [27]. In all cases, attention is needed to maintain turbulence during the transition between RANS and LES.

Despite recent advancements in hybrid RANS/LES methods, there are still issues that must be resolved. This is a result of the difficulties associated with managing interface boundary conditions (IBCs). LES equations predict fluctuations, whereas RANS equations do not. The addition of artificial fluctuations introduces its own uncertainties [18]. The log–layer mismatch (LLM) is another significant issue in hybrid RANS/LES simulations. This is evident even in planar channel flow simulations [20,28,29] and leads to overestimation or underestimation of the skin friction. The formation of an artificial buffer layer and consequently unrealistic physical structures can contaminate the LES area.

To avoid the LLM, Keating and Piomelli [30] recommend stochastic forcing at the interface area. This can expedite the formation of resolved eddies and lead to improved outcomes. Hamba [20] introduces additional filtering. Zhong and Tucker [31] use the $k-l$ method to approximate the eddy viscosity. Davidson and Peng [25] utilize the $k-\omega$ for their eddy viscosity model in the RANS region. Davidson and Dahlström [32] increase the resolved shear stress by introducing fluctuations at the interface. Temmerman et al. [6] determine the eddy viscosity in the RANS region using the $k-l$ or $k-\omega$ method; the eddy viscosity of the RANS interface is matched with the LES subgrid viscosity. Davidson and Billson [26] add fluctuations at the interface to increase the resolved shear stress. Larsson et al. [5] use additional forcing to induce resolved motions in the LES region. They discuss that matching the turbulent viscosity at the interface has a minimal effect on the final solution.

However, these solutions require additional data from a preliminary benchmark solution or contain user-defined parameters (such as the forcing amplitude) to increase resolved turbulent fluctuations at the interface. The model accuracy varies considerably based on the values of these free parameters, and robust selection criteria for the forcing-method parameters are unavailable, to the best of our knowledge. Consequently, LLM still occurs for the majority of channel flow simulations when hybrid RANS/LES models are used. This illustrates the importance of accurate enough modeling the inner region and the connection between the inner and outer regions. According to Piomelli [7], the RANS/LES methods are least accurate when applied to attached and thin shear layers. They are seemingly better suited for separated flows over curved surfaces and with adverse pressure gradients.

Even with a near-wall model, the grid resolution in the outer LES region should be good enough. Kawai and Larsson [33] highlight the importance of grid independence studies in WSMs or hybrid models with an auxiliary RANS model. LES must be well-resolved at the interface point to provide accurate enough information to the RANS model. They also show that the auxiliary RANS model should not necessarily be located between the first off-wall grid point and the wall because the LES model is under-resolved in the near-wall region and is therefore contaminated by truncation and subgrid-scale (SGS) model errors [34]. The effect of decreasing the interface boundary height on reducing accuracy also occurs in [28]. Keating and Piomelli [30] obtain better results by setting the interface far enough after the buffer region. This helps them generate the RANS logarithmic layer before the start of the LES region.

The near-wall domain decomposition (NDD) was first developed for high-Reynolds RANS [35,36,37]. The original approach to NDD presumes the computation domain to be split into the inner and outer regions. In the inner region, a simplified locally one-dimensional TBLE model with a prescribed turbulent viscosity can be used [38]. For the outer region, the effect of the inner region is replaced by nonlocal IBCs, which are obtained by the transfer of the original BCs from the wall to the interface boundary [37]. As demonstrated in [39], the NDD can be realized with the full set of RANS equations, although this affects the rate of convergence. A principally new approach called the implicit NDD is proposed in [40]. In this case, the problem can be solved on a relatively coarse grid in the entire region with specially derived computational slip BCs at the wall. Then, the solution is recomputed in the inner region with the use of the TBLE model. The derived slip BCs at the wall guarantee that the composite solution is smooth up to the first derivative. A significant advantage of this approach is that it does not require any explicit split of the computational domain into two sub-domains.

In the current paper, the implicit NDD (INDD) [40] is extended to the RANS/LES zonal decomposition. As a result, in the first (predictor) stage, only LES is used with the slip BCs at the wall. Then, in the second (corrector) stage, the solution in the inner region is recomputed. In this way, the RANS region is embedded onto the LES domain. The computing time needed for the RANS solution is almost negligible since a one-dimensional (TBLE) model is considered in the inner region. In this work, in contrast to the original NDD, the turbulent viscosity is determined from the $k-l$ model rather than prescribed.

The remainder of this paper is structured as follows. Section 2 describes the model used in the inner region, including the slip BCs. Then, in Section 3, the INDD algorithm is provided. In turn, the configuration and characteristics of the channel flow test case are described in Section 4. The obtained results are presented there and compared with other approaches, including DNS, unresolved LES, and conventional hybrid methods. Finally, conclusions are drawn, and recommendations and future prospects are discussed in Section 5.

2. Near-Wall Model

2.1. Slip Boundary Conditions

In the implicit NDD, the solution is initially obtained on a relatively coarse mesh using LES and slip BCs at the wall. Then, a TBLE is solved in the inner region on a one-dimensional mesh to obtain the velocity with the no-slip BC at the wall and the velocity taken from the LES at the opposite end. In this way, the grid used for RANS is independent from the grid used for LES. At the end, the complete solution consists of the RANS solution in the inner region and the LES solution beyond the interface (see Figure 1). The red dashed line represents the preliminary velocity profile in the inner region resulting from the slip boundary condition. The slip BCs for LES guarantee that the composite time-averaged solution is smooth at the interface. The same procedure is applied to both tangential components of the velocity.

Let us consider a one-dimensional RANS region embedded into the LES mesh between $y=0$ and $y=y^{*}$ and follow the TBLE format of the boundary-layer equation:

$${\displaystyle \frac{\partial }{\partial y}}\left(\mu {\displaystyle \frac{\partial u}{\partial y}}\right)=R_u\left(y\right),$$

(1)

where u is the streamwise (wall-parallel) velocity, $\mu \left(y\right)$ is the total viscosity as the sum of the laminar (${\mu }_l$) and turbulent (${\mu }_T$) viscosities, and $R_u\left(y\right)$ includes non-equilibrium terms such as the pressure gradient and convective terms. As shown in [40], for a small enough $y^{*}$ the convective terms can be neglected with low effect on the accuracy. Thus, in the current study, $R_u$ corresponds to the inverse pressure gradient along the wall taken at $y^{*}$.

For Equation (1), the Dirichlet BCs are set:

$$\left\{\begin{array}{cc}u=0\hfill & \mathrm{at}y=0,\hfill \\ u=u^{*}\hfill & \mathrm{at}y=y^{*}.\hfill \end{array}\right.$$

According to [35,36,37], the BC can be exactly transferred from the wall $y=0$ to the interface boundary $y=y^{*}$:

$$u\left(y^{*}\right)=f_1{\left.{\displaystyle \frac{du}{dy}}\right|}_{y^{*}}-\widetilde{f_2},$$

(2)

where

$$f_1={\int }_0^{y^{*}}{\displaystyle \frac{\mu \left(y^{*}\right)}{\mu \left(y\right)}}dy,{\tilde{f}}_2={\int }_0^{y^{*}}{\displaystyle \frac{1}{\mu \left(y\right)}}\left({\int }_y^{y^{*}}{R}_{u}d\xi \right)dy.$$

(3)

It is worth noting that (2) cannot be obtained by double-integrating the governing Equation (1) from the wall since the derivative of the solution at the wall is unknown.

For the outer region, Equation (2) is fully equivalent to the governing Equation (1) and the no-slip BC at the wall [36]. To simplify the decomposition algorithm, at the implicit NDD [40] it is proposed to transfer the IBC (2) back to the wall with $\mu$ and $R_u$ taken at $y^{*}$. This leads to the computational slip BC at the wall [40]:

$$u\left(0\right)=f_{w1}{\left.{\displaystyle \frac{du}{dy}}\right|}_{y=0}+f_{w2},$$

(4)

with

$$f_{w1}=f_1-y^{*},\mathrm{and}{f}_{w2}=\widetilde{f_2}+{\displaystyle \frac{y^{*}{R}_u^{*}}{{\mu }^{*}}}\left(f_1-{\displaystyle \frac{y^{*}}{2}}\right),$$

(5)

where symbol * corresponds to the values taken at $y^{*}$.

The aim of the slip BC (4) is two-fold. It allows the LES solution to be formally obtained up to the wall on a relatively coarse grid. In addition, it guarantees that the time-averaged composite RANS/LES solution is smooth up to the first derivative [40]. The solution obtained in the preliminary (LES) stage is schematically shown in the inner region by the red curve in Figure 1.

As soon as the LES solution with the slip BC has been obtained, the RANS solution in the inner region can be updated. To update the shear stress (${\tau }_w$), Equation (1) is integrated twice using the velocity at $y^{*}$ (already known from LES) which leads to

$${\tau }_w={\displaystyle \frac{u\left(y^{*}\right)-{\int }_0^{y^{*}}{\displaystyle \frac{{\int }_0^y{R}_{u}d\xi }{\mu }}dy}{{\int }_0^{y^{*}}{\displaystyle \frac{dy}{\mu }}}}.$$

(6)

It is to be noted here that $u\left(y^{*}\right)$ corresponds to an instantaneous velocity from LES. Next, the RANS velocity in the entire inner region is updated from the exact solution of the governing Equation (1):

$$u\left(y_i\right)={\tau }_w{\int }_0^{y_i}{\displaystyle \frac{dy}{\mu }}+{\int }_0^{y_i}{\displaystyle \frac{{\int }_0^y{R}_{u}d\xi }{\mu }}dy.$$

(7)

2.2. Turbulent Viscosity Estimation

In the aforementioned equations, the turbulent viscosity must be evaluated. It can be estimated from algebraic correlations or from k-based methods. In this work, we examine three approaches, namely the mixing length model, the $k-\omega$ method, and the $k-l$ method by Wolfshtein [41,42]. In further consideration, by default, we imply the Cartesian coordinate system and Einstein notation for repeated indices, which are supposed to vary from 1 to 3.

Model 1: The mixing length model by Cabot [43] including a damping function is used to estimate the turbulent viscosity:

$${\nu }_T=\kappa y{u}_{\tau }{D}^2,D=1-\exp (-y{u}_{\tau }/\nu A^{+}),$$

(8)

where $\kappa =0.41$, $A^{+}=19$, $u_{\tau }=\sqrt{|{\tau }_w|/\rho }$, and ${\nu }_T={\mu }_T/\rho$.

Model 2: We use the two-equation $k-\omega$ model by Wilcox [44]. The k equation in this model is given by

$${\displaystyle \frac{\partial k}{\partial t}}+{\displaystyle \frac{\partial u_{j}k}{\partial x_j}}={\displaystyle \frac{\partial }{\partial x_j}}\left[\left({\nu }_l+{\sigma }_k{\nu }_T\right)\left({\displaystyle \frac{\partial k}{\partial x_j}}\right)\right]+P_k-{\beta }^{*}k\omega ,$$

(9)

where ${\sigma }_k=1/2$ is the diffusion Prandtl number, ${\beta }^{*}=0.09$, and $u_j$ is the RANS velocity, while $P_k$ is the turbulence production term that reads

$$P_k={\tau }_{ij}^{turb}{\displaystyle \frac{\partial u_i}{\partial x_j}}.$$

(10)

Here, ${\tau }_{ij}^{turb}$ corresponds to the Reynolds stress tensor:

$${\tau }_{ij}^{turb}=\frac{{\nu }_T}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_i}{\partial x_j}\right),(i, j=1, 2, 3).$$

The transport equation for $\omega$ reads

$${\displaystyle \frac{\partial \omega }{\partial t}}+{\displaystyle \frac{\partial u_j\omega }{\partial x_j}}={\displaystyle \frac{\partial }{\partial x_j}}\left[\left({\nu }_l+{\sigma }_{\omega }{\nu }_T\right){\displaystyle \frac{\partial \omega }{\partial x_j}}\right]+P_{\omega }-\beta {\omega }^2.$$

(11)

Here, $\beta =3/40$, ${\sigma }_{\omega }=1/2$, $ϵ={\beta }^{*}\omega k$ and the production term $P_{\omega }$ reads

$$P_{\omega }=\alpha {\displaystyle \frac{\omega }{k}}{\tau }_{ij}^{turb}{\displaystyle \frac{\partial u_i}{\partial x_j}},$$

(12)

where $\alpha =5/9$. Turbulent viscosity is then evaluated as

$${\nu }_T={\displaystyle \frac{k}{\omega }}.$$

(13)

As discussed below, k and $\omega$ equations are solved in the approximation of TBLE.

Model 3: We also adapt the $k-l$ method [45] based on the transport of turbulent kinetic energy. In this case, we have a slightly different form of the k equation:

$${\displaystyle \frac{\partial k}{\partial t}}+{\displaystyle \frac{\partial u_{j}k}{\partial x_j}}=\frac{\partial }{\partial x_j}\left[\left({\nu }_l+{\displaystyle \frac{{\nu }_T}{{\sigma }_k}}\right){\displaystyle \frac{\partial k}{\partial x_j}}\right]+P_k-ϵ.$$

(14)

Here, ${\sigma }_k=1$. The boundary value for k at the interface ($y^{*}$) is equal to $k^{*}=k^{res}+k^{sgs}$, where $k^{res}$ is the time-averaged resolved and dominant part of the turbulent kinetic energy in LES, and $k^{sgs}$ is its sub-grid scale part; $k^{res}=0.5\overline{u_i^{\prime }{u}_i^{\prime }}$, where $u_i^{\prime }(i=1, 2, 3)$ are turbulent pulsations.

Under the TBLE assumption, Equation (14) can be rewritten in the 1D form along the wall normal direction (y) as

$${\displaystyle \frac{\partial }{\partial y}}\left[\left({\nu }_l+{\displaystyle \frac{{\nu }_T}{{\sigma }_k}}\right){\displaystyle \frac{\partial k}{\partial y}}\right]+{\nu }_T{\left({\displaystyle \frac{\partial u}{\partial y}}\right)}^2-ϵ=0.$$

(15)

This is a reasonable approximation for a low enough y. As noted in [24], low-$\mathrm{Re}$ RANS models often yield acceptable solutions when the normal to the wall direction is predominant. Therefore, high aspect ratio cells (order of 100–1000) can be used in the wall vicinity. This gives a superiority over resolved LESs, where fine streamwise meshes are required. It is worth noting that the embedded mesh for any fixed coordinate x along the wall is entirely local. In the aforementioned $k-l$ equation, we have

$$ϵ=C_{ϵ}{k}^{3/2}/l_{ϵ},$$

(16)

and turbulent viscosity is evaluated through

$${\nu }_T=C_{\mu }{l}_{\mu }{k}^{1/2},$$

(17)

where $C_{ϵ}=1$ and $C_{\mu }=0.09$. The length scales in this RANS model are defined as

$$l_{ϵ}=2.4y\left(1-\exp (-0.263y^{\times })\right),$$

(18)

$$l_{\mu }=2.4y\left(1-\exp (-0.016y^{\times })\right),$$

(19)

where $y^{\times }=y\rho k^{1/2}/\mu$.

It is clear that the $k-l$ method is computationally less expensive than the $k-\omega$ model. We use the $k-l$ model in all cases, unless otherwise stated.

3. Implicit Near-Wall Domain Decomposition

3.1. Algorithm

To summarize, the implicit NDD algorithm for the zonal RANS/LES model is implemented as follows:

• Initialize the coarse grid for LES.
• Initialize the interface boundary $y^{*}$ and a sub-grid for RANS.
• Initialize the flow fields for both the unresolved LES grid and RANS sub-grid.
• Compute the turbulent viscosity.
• Compute $f_1$ and ${\tilde{f}}_2$ (Equation (3)) based on the chosen turbulent viscosity model of Section 2.2.
• Compute the coefficients $f_{w1}$ and $f_{w2}$ (Equation (5)) to impose the slip boundary condition at the wall.
• Solve the LES governing equations on the coarse grid with the slip boundary condition (Equation (4)).
• Transfer the LES streamwise velocity and turbulent kinetic energy at $y^{*}$ (i.e., $u^{*}$ and $k_{tot}^{*}$) to the embedded RANS model.
• Compute the wall shear stress in the inner region using $u^{*}$ value (Equation (6)).
• Compute the RANS velocity solution in the inner region (Equation (7)) with the updated ${\tau }_w$.
• Compute the turbulent kinetic energy (Equation (14)), if needed.
• Repeat the procedure from step 4.

3.2. Discussion

The INDD effectively combines the zonal and wall-stress approaches and has some advantages over conventional methods, as discussed below.

(I) The solution at the interface is smooth because the slip boundary condition is transferred to the wall. In addition, the composite time-averaged solution is also smooth. This follows from the design of the computational slip BCs [40].

(II) Since we use RANS equations in the inner region, we avoid the issues regarding the selection of a proper $k^{*}$ value as a BC for the RANS region (see, e.g., [6]). If $u^{LES}$ is used in the inner region for the evaluation of k from Equation (14), the contribution of resolved turbulent kinetic energy ($k_{res}$) into the inner region is already taken into account; therefore, $k^{*}$ should only include the modeled part, i.e., $k_{sgs}$. However, Temmerman et al. [6] argued that the issue is not that simple and other assumptions are needed regarding the correct $k^{*}$ value. In the current work, we take $k^{*}=k_{sgs}+k_{res}$ at the boundary of the RANS region because $u^{RANS}$ is used in the inner region to predict k.

(III) In the current approach, since both the RANS velocity and turbulent kinetic energy are obtained on a fine 1D mesh in the inner region, the grid independence of the RANS solution can be easily ensured. Moreover, using u from the RANS solution in the inner region can help prevent the formation of an “artificial buffer layer”, that is, a transition zone between modeled and resolved turbulence. Larsson et al. [5] discussed that the artificial buffer layer is an effect of the merging of two different modeling regions. We additionally show that the current algorithm does not need any enforcement of external force or synthetic turbulence at the interface.

(IV) In comparison to WSMs, according to Larsson et al. [8], while LESs can impart flow structures onto the wall-model, it is problematic for the flow structures in the wall-model to penetrate into the LES region. Therefore, in WSMs, the RANS/LES coupling is rather weak. In the INDD, this coupling is enhanced via the RANS solution in the entire inner region and non-local IBCs.

(V) One shortcoming of the current approach against WSMs is that it is computationally more expensive.

4. Test Cases

4.1. Setup

Simulations are carried out with the open source and finite-volume-based CFD software OpenFOAM distributed by the OpenFOAM Foundation [46]. Planar channel flow at three frictional Reynolds numbers, namely ${\mathrm{Re}}_{\tau }=950$, 2000, and 4200, is considered, where the friction Reynolds number is defined as

$${\mathrm{Re}}_{\tau }={\displaystyle \frac{u_{\tau }h}{\nu }}.$$

(20)

Here, h is the half channel height, $\nu$ is the fluid kinematic viscosity, and $u_{\tau }=\sqrt{|{\tau }_w|/\rho }$ is the friction velocity.

Domain size, grid resolution, and the interface location $y^{*}$ are listed in

Table 1. Characteristics of studied test cases.

Case	${\mathbf{Re}}_{\mathit{\tau }}$	${\mathit{L}}_{\mathit{x}}\times {\mathit{L}}_{\mathit{y}}\times {\mathit{L}}_{\mathit{z}}$	Resolution	$\Delta {\mathit{x}}^{+}$	$\Delta {\mathit{y}}^{+}$	$\Delta {\mathit{z}}^{+}$	$\Delta {\mathit{y}}_{\mathit{j}+1}/\Delta {\mathit{y}}_{\mathit{j}}$	${\mathit{y}}^{*}/\mathit{h}$	${{\mathit{y}}^{*}}^{+}$
C950	950	$2\pi h\times 2h\times \pi h$	$40\times 64\times 36$	148	$1.57\to 126$	83	1.15	0.063	60
C2000	2000	$2\pi h\times 2h\times \pi h$	$40\times 72\times 36$	312	$1.78\to 267$	174	1.15	0.065	129
C4200	4200	$3\pi h\times 2h\times \pi h$	$60\times 84\times 42$	659	$1.82\to 547$	314	1.15	0.0476	200

, where x, y, and z represent the streamwise, wall-normal, and spanwise directions, respectively. The streamwise and spanwise directions have periodic boundary conditions. At the wall, $y=0$, we originally set the no-flux and no-slip BCs. The channel half-height is $h=1$ m. The channel flow is driven by a driving force determined from the DNS-computed bulk velocity.

The mesh is uniform in the x and z directions, and refined in the y direction with a specified stretch factor. In the inner region, the embedded mesh usually consists of 20−30 nodes.

The results of different approaches, including unresolved LES, implicit NDD, and benchmark DNS are compared. The first two models are implemented in the current study and compared to the available benchmark DNS data. In each test case, all studied models (except for DNS) have identical time and space resolutions. In addition, the effects of three turbulent viscosity models of Section 2.2 are compared with each other.

The WALE model [47] is used for the subgrid-scale stress modeling in the LES. The model is originally implemented in OpenFOAM. Each simulation starts with a fully developed turbulent channel flow and allows for enough running time before capturing the output values.

The interface boundary in the simulations is situated at $0.05h\lesssim y^{*}\lesssim 0.1h$ as recommended in [30], where h is the channel half-height.

4.2. Simulation Results

In this section, the simulation results for the mean streamwise velocity and Reynolds normal stress intensities obtained by INDD are compared to DNS and unresolved LES (i.e., LES computed on a relatively coarse grid without any special near-wall modeling). In the inner region, the $k-l$ eddy viscosity model is used.

The normalised mean streamwise velocity $u^{+}=\overline{u}/u_{\tau }$ for ${\mathrm{Re}}_{\tau }=950$ is plotted in Figure 2 against the normalised distance from the wall $y^{+}=y{u}_{\tau }/\nu$ and compared to the DNS data from Hoyas and Jiménez [48]. It is evident that the unresolved LES simulation significantly overestimates the mean streamwise velocity. The unresolved LES can only predict well enough the velocity profile near the viscous sub-layer ($y^{+}<10$), whereas in the logarithmic region the velocity exhibits a strong upward lift. This situation is substantially improved by using the INDD, as the agreement with DNS in both the inner and outer regions is satisfactory. Here and further, the location of the interface is indicated by a vertical dashed line. Intriguingly, no LLM occurs with the current approach, despite avoiding additional force terms in the momentum equation. That can be attributed to the appropriate exchange of flow physics between the RANS and LES regions. With the TBLE approximation and focus on the wall-normal gradients, we are able to use wall-parallel cells, which are much coarser than the local eddies. They cannot be appropriately resolved with the conventional SGS models.

The normalised Reynolds stress intensities (${\overline{u^{\prime }{u}^{\prime }}}^{+}$, ${\overline{v^{\prime }{v}^{\prime }}}^{+}$, and ${\overline{w^{\prime }{w}^{\prime }}}^{+}$) for ${\mathrm{Re}}_{\tau }=950$ are depicted in Figure 3. In turn, Figure 4 illustrates the normalized Reynolds shear stress for C950. The overall agreement with the DNS is rather good. It is to be noted that the prediction of the Reynolds stress intensities with INDD in the inner region is rather formal. It corresponds to LES with the slip boundary conditions in the preliminary stage. The LES is not intended to be resolved in this region. The overall prediction of INDD is rather accurate, with the exception of the interface region. Near the RANS region, we cannot expect a high accurate prediction of the Reynolds stress intensities. A similar issue also occurs in other works related to RANS/LES models (see, e.g., [30]). This happens because RANS cannot predict turbulent velocity fluctuations. As a result, the Reynolds stresses in the interface area cannot be well predicted.

The results for the mean streamwise velocity at ${\mathrm{Re}}_{\tau }=2000$ are displayed in Figure 5. They are compared to DNS, unresolved LES, and those from Davidson and Dahlström [32] (shown as DD2005) both without and with a force term. The mean velocity of the unresolved LES is significantly overpredicted, which is due to the substantial underprediction of the friction velocity (around 40%). Without incorporating additional force terms into the momentum equation, the results obtained with INDD are in a rather good agreement with the DNS data. Our prediction almost coincides with that obtained by Davidson and Dahlström [32] when they included a force term. However, the nature and magnitude of the added force term can change their profile considerably. As depicted in the figure, when no forcing is applied in [32], a large deviation from the DNS data occurs.

The Reynolds normal stress intensities for ${\mathrm{Re}}_{\tau }=2000$ are shown in Figure 6. The INDD results correspond to ${y^{*}}^{+}=129$, while ${y^{*}}^{+}=62$ in [32]. Again, away from the interface, the agreement of INDD with the DNS data is rather good. Despite taking a larger $y^{*}$ value than DD2005, INDD demonstrates superior agreement with the benchmark data for all components. The velocity fluctuations at ${y^{*}}^{+}=129$ are not available in [32] but we anticipate for DD2005 a larger shift to the right and, thus, a greater difference from DNS in this case. In turn, Figure 7 shows the relative values of $INDD-DNS$. The main difference is near the interface area and it fades away when approaching the channel center. The hybrid model predicts the maximum normal stress to be at the interface, while in DNS it is essentially closer to the wall. This happens because LES is under-resolved in the inner region. As a result, the pulsations $u^{\prime }$ sharply drop there.

The shear stress profile is plotted in Figure 8. Considering the coarse mesh used here and the deficiency of the RANS model in predicting fluctuations, the overall agreement with DNS is reasonably good. The hybrid model predicts the maximum point close to the interface.

Figure 9 demonstrates the results obtained with the INDD and our conventional hybrid RANS/LES simulation without forcing (like [5]) on the same mesh. With the hybrid approach, the LES solution is computed in the entire region in such a way that, in the inner region, the subgrid eddy viscosity is replaced by the eddy viscosity predicted from Model 3. It can be seen that a mismatch at the interface is reproduced with this approach. As noted in [49], the implementation of a blending function for the eddy viscosity is able to reduce the mismatch at the interface, although it does not remove it beyond the interface. In contrast, the INDD profile is more in line with the DNS benchmark, and the velocity profile is smoother even without the use of any blending function. INDD significantly mitigates the LLM and does not have a velocity jump.

The mean streamwise velocity for ${\mathrm{Re}}_{\tau }=4200$ is illustrated in Figure 10. Overall, the agreement between the current study and the DNS [50] for both inner and outer regions is rather good. The results of Keating and Piomelli [30] (shown as KP2006), when no forcing is applied, are overpredicted. The addition of the force term improves the KP2006 results, but it requires the selection of free forcing parameters that are not known a priori. This issue is addressed in a recent review by Larsson et al. [8]. They note that the log–layer mismatch is a linear function of the forcing amplitude. Finally, it is concluded that, since there is no robust general method to determine the forcing amplitude, the concept of adding force terms to initiate the transition from RANS to LES is not a comprehensive strategy.

As shown in Figure 11, the Reynolds stress intensities predicted by INDD are closer to the benchmark data with the increase in the distance from the interface. Specifically with regard to the streamwise component, the prediction with the INDD is comparable with that of KP2006, where a force term is used. The result of INDD for the shear stress in Figure 12 shows a reasonable agreement with the benchmark data considering the coarse mesh utilized in the current study.

As known, the LLM results in streak-like areas in the regions in the vicinity of the interface [30,51]. This is different from the low-speed streaks that trigger the turbulence generation cycle. Figure 13 depicts a snapshot of the streamwise velocity at $y/H=0.04$. Clearly, the “super-streaks” that are distinct to the LLM in the conventional RANS/LES model are not observed here.

Overall, for all Reynolds numbers considered, either minor or negligible LLM has been observed in cases where INDD is used.

4.3. Effect of Eddy Viscosity Model

Next, the influence of the eddy viscosity model on the mean velocity profile for ${\mathrm{Re}}_{\tau }=4200$ is studied. Comparisons are made in Figure 14. The results obtained with the $k-l$ method, $k-\omega$, and mixing length models are compared against the DNS data. As can be seen, the $k-l$ model has a clear superiority. The mixing length and $k-\omega$ models significantly underestimate the mean velocity in the core region. This result is consistent with the statement by Temmerman et al. [6], who argue that unresolved LES cannot provide a reliable prediction for $ϵ$ or $\omega$ because the appropriate evaluation of the dissipation rate requires a fine-enough mesh as dissipation occurs at small scales.

5. Conclusions

In this paper, the implicit near-wall domain decomposition approach has been extended to hybrid RANS/LES simulations. In the approach, the solution is obtained in two stages. First, an LES solution is computed for the entire region with specially derived non-local slip boundary conditions at the wall. Then, in the second stage, the solution is recomputed in an inner region near the wall based on a TBLE approximation. The slip BCs at the wall guarantee that the time-averaged composite solution is smooth. To evaluate the eddy viscosity in the inner region, it is recommended to use the $k-l$ turbulence model. The obtained results demonstrate that the approach is able to provide rather accurate predictions in a wide range of Reynolds numbers without the implementation of any additional forcing term at the interface. The velocity profile corresponds to the DNS data reasonably well, such that the LLM effect is practically negligible. The conventional hybrid models often suffer from a large LLM. They may include extra force terms or synthetic turbulence fluctuations at the interface. The amplitude of the force term is case-dependent and there is not a universal approach to determine it. In contrast, INDD does not need any additional force term. In the current study, we applied the model to frictional Reynolds numbers of 950, 2000, and 4200. In all cases, the mean velocity profile is well predicted. However, better results are obtained for higher Reynolds numbers. The current work demonstrates that the developed approach can be successfully used to reduce or eliminate LLM. In future work, the approach will be tested on flows around complex geometries.