# An active hybrid Reynolds-Averaged Navier-Stokes/Large Eddy Simulation approach for grey area mitigation

Mahitosh Mehta, Remi Manceau, Vladimir Duffal, Benoit de Laage de Meux

► To cite this version:

Mahitosh Mehta, Remi Manceau, Vladimir Duffal, Benoit de Laage de Meux. An active hybrid
Reynolds-Averaged Navier-Stokes/Large Eddy Simulation approach for grey area mitigation. Physics
of Fluids, 2023, 35, 10.1063/5.0174381. hal-04287062

HAL Id: hal-04287062

https://inria.hal.science/hal-04287062v1

Submitted on 15 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An active hybrid RANS/LES approach for grey area mitigation

Mahitosh Mehta and Rémi Manceau*

University of Pau & Pays Adour, E2S UPPA, CNRS, Inria, CAGIRE project-team,
Department of mathematics and applied mathematics, Pau, France

Vladimir Duffal and Benoît de Laage de Meux

EDF R&D, MFEE Dept., 6 quai Watier, 78400 Chatou, France

(*Corresponding author: remi.manceau@univ-pau.fr)

(Dated: November 15, 2023)

During the transition from the RANS (Reynolds-Averaged Navier-Stokes) mode to the LES (Large Eddy Simulation) mode, i.e., in the so-called grey area, continuous hybrid RANS/LES approaches suffer from the well-known problem of excessively slow generation of resolved structures. Indeed, when the mesh is refined in the direction of the flow, the model is designed to reduce the modeled energy, but there is no mechanism to transfer the equivalent amount of energy into the resolved motion. Hence, the total turbulent energy and turbulent stresses are underestimated, which strongly affects the prediction of the mean flow. This also constitutes a violation of the conservation of mechanical energy, which can only be corrected by an *active* approach, i.e., an approach that allows the injection of resolved energy. The aim of this work is to develop such an active approach based on the introduction of a fluctuating volume force into the resolved momentum equation, similar to the anisotropic linear forcing (ALF) method proposed previously. The major difference with ALF is that the new method does not require target statistics obtained by a RANS computation, but is based on a simple analysis of the rate of energy transfer related to variations in resolution, enabling the forcing to be extended to continuous hybrid RANS/LES. Application of the new method to the cases of a channel with or without periodic constriction shows a drastic improvement over the case without forcing. Although the method is applied herein to a particular hybrid RANS/LES approach (HTLES, hybrid temporal LES), it can easily be extended to any other approach, as long as a parameter identifies variations in resolution, and thus offers vast application prospects.

## I. INTRODUCTION

Computational fluid dynamics (CFD) is a powerful tool for predicting complex turbulent flows. Although RANS (Reynolds-Averaged Navier-Stokes) computations are the industrial standard practice, they often do not provide accurate enough results in separated flow regions and also do not provide time-dependent information.

For instance, it is of great importance to describe the unsteady behavior of large-scale energetic structures of turbulence to predict aerodynamic loads and thermal fatigue in industrial applications. LES (Large-Eddy Simulation) can provide these quantities but at a high computational cost, which motivates the use of hybrid RANS/LES models which offer the best potential to achieve the balance between LES and RANS advantages. Indeed, LES can capture large-scale structures in regions of interest, while RANS models are applied to regions where LES is not needed or too CPU-demanding. Over the last couple of decades many continuous hybrid RANS/LES approaches have been proposed, such as Detached Eddy Simulations (DES) and variants,¹⁻³ Partially Averaged Navier-Stokes (PANS⁴), Scale-Adaptive Simulations (SAS⁵), Partially Integrated Transport model (PITM⁶), Continuous Eddy Simulation (CES⁷), to name a few (see, for instance Chaouat⁸ or Heinz,⁹ for recent reviews).

The main purpose of any continuous hybrid RANS/LES method is to reduce the modeled tur-

bulent stresses in the LES region either by modifying the eddy-viscosity formulation, or by manipulating terms in the equations. While transitioning from RANS to LES, we go from fully modeled to partially resolved quantities and by doing so, the reduction in modeled energy should be compensated by the increase in resolved energy. However, the growth of the resolved fluctuations relies entirely on natural instabilities, and is generally not fast enough to compensate for the stress decay imposed in the model, a phenomenon known as Modeled Stress Depletion.¹⁰ This is illustrated by Fig. 1, which shows Q-isocontours for the periodic hill case, where upon transitioning from RANS to LES, it is seen that it takes a long time to generate the resolved structures. This lack of resolved energy in the transition zone from RANS to LES is also known as the *grey area* problem, and significantly degrades the mean flow predictions. This problem can be alleviated by decreasing to zero the subgrid viscosity in the initial region of a detached shear layer, in order to avoid this viscosity delaying the development of instabilities.¹¹ However, this solution relies essentially on the presence of inflectional velocity profiles and does not inject energy into the resolved motion, such that, as shown below, mechanical energy is still lost.

The problem of the absence of energy transfer between modeled and resolved parts can be clearly illustrated by a thought experiment, the simplified case of a RANS-LES transition in grid turbulence, as shown in Fig. 2:
<figure><img src="image_2.png" /><figcaption>FIG. 1: Isocontours $Q = 0.2U_b^2/h^2$ colored with the velocity magnitude and, in the background, modeled-to-total turbulent energy ratio $r$, indicating the regions solved in RANS ($r = 1$) and in LES ($r < 1$). Here the HTLES model of Duffal et al.¹² is used.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2: Illustration in grid turbulence of the problem of the absence of turbulent energy transfer between the modeled part and the resolved part. (a) Schematic view of a grid turbulence experiment. (b) Schematic evolution of the turbulent energy partition during a RANS-to-LES transition in hybrid RANS/LES.</figcaption></figure>

this corresponds to famous grid turbulence experiments, as measured for instance by Comte-Bellot and Corrsin.¹³ In such a flow, turbulence is generated by the grid and convected downstream by a uniform flow. Due to the absence of a velocity gradient, there is no turbulence production, and the turbulent energy decreases as a power law of $x$. If an eddy-viscosity model is used for the sub-filter stress, the molecular and turbulent diffusion terms in the resolved momentum equation write

$$ \frac{\partial}{\partial x_j} \left[ (\nu + \nu_{sfr}) \frac{\partial \tilde{u}_i}{\partial x_j} \right] \quad (1) $$

(the tilde stands for the filtering operator as defined in section II). In the RANS region, the resolved velocity is the mean velocity, i.e., the uniform velocity field. Consequently, in the momentum equation, the velocity gradient $\partial \tilde{u}_i / \partial x_j$, the convection term $\tilde{u}_j \partial \tilde{u}_i / \partial x_j$ and the diffusion term (1) are zero. Therefore, the resolved momentum equation reduces to

$$ \frac{1}{\rho} \frac{\partial \tilde{p}}{\partial x_i} = 0 \quad (2) $$

whose solution is $0 = 0$. When at a given position the mesh is refined to switch to LES mode, the hybrid

RANS/LES model strongly decreases the turbulent viscosity. However, since the turbulent viscosity does not appear in Eq. (2) as explained above, the RANS solution $0 = 0$ remains a solution of the equation. Resolved structures will appear only slowly, due to the natural instability of the flow. However, as the meshes used are much coarser than those required for DNS, this natural transition can be strongly delayed, or completely suppressed. Fig. 2 shows the desired result: the modeled energy $k_m$ is strongly reduced due to the refinement, which is compensated by an increase of the resolved energy $k_r$, so that the total energy $k = k_r + k_m$ remains insensitive to the refinement. We will in practice obtain $k_r \simeq 0$ and the energy $k$ will decrease drastically, following $k_m$.

Two fundamental conclusions can be drawn from this thought experiment:

i) It is obvious that no modification of the model for the unresolved part can correct this problem, since the model is not involved in Eq. (2).

ii) There is actually here a fundamental problem of violation of the conservation of mechanical energy.

It will be seen below that the decrease of the modeled en-
ergy corresponds to a modification of the partition mod-
eled/resolved energy, which generates terms in the mo-
mentum equation that are not taken into account. This
loss of turbulent energy, and thus the decrease of the to-
tal turbulent stresses which drive the mean momentum,
affects the mean velocity profiles in turbulent shear flows.
Since the missing terms are actually forces, the goal of
this paper is develop a volume forcing method to inject
energy into the resolved motion.

It is legitimate to wonder whether it is necessary to cor-
rect this problem in all cases. Indeed, many applications
of continuous hybrid RANS/LES models are available in
the literature that do not introduce this type of forcing:
if fluctuations appear naturally at the transition from
RANS to LES, which is often the case in the presence of
a detached shear layer, the solution is acceptable, and it
is possible to promote the appearance of resolved struc-
tures, for example by radically decreasing the subgrid
viscosity in this zone.¹¹ However, it is important to note
that in this case there is still a violation of the conser-
vation of mechanical energy: there is no mechanism that
transfers modeled energy to resolved energy. The growth
of the resolved fluctuating energy is therefore necessar-
ily obtained by taking from the mean flow this missing
part of the mechanical energy, which consequently alters
the mean flow. It is therefore preferable, in all cases, to
correct this mechanical energy conservation issue.

To address this issue, a few methods have been de-
veloped that introduce, in one way or another, resolved
fluctuations. A comparative study of several methods
was performed by Keating et al..¹⁴ Synthetic turbulence
methods are generally based on Fourier reconstruction
techniques, where spatial and temporal correlations can
be introduced.¹⁵,¹⁶ Other approaches introduce virtual
vortices of simple shape, either by directly superimpos-
ing them to the velocity field or by introducing a restoring
force in the momentum equation.¹⁷,¹⁸ Spille-Kohoff and
Kaltenbach¹⁹ proposed an approach based on a volume
forcing of the wall-normal velocity in order to control
the shear stress in a set of planes close to the LES inlet.
Laraufie et al.²⁰ modified this approach by controlling
the diagonal Reynolds stress in the wall-normal direction
rather than the shear stress. A common feature of all
these approaches is that the interface between the RANS
and LES zones must be geometrically simple (typically a
plane normal to a shear layer) and the flow direction well
identified. Indeed, when synthetic vortices are generated,
a great deal of information is required to adapt their size,
impose a possible periodicity in the transverse direction,
and advect the vortices. Typically, the size and energy
of the vortices must be obtained from existing informa-
tion, such as upstream statistics from the RANS region
or distance to the wall. Approaches based on volume
forcing are also linked to a particular component of the
Reynolds stress tensor and therefore to the relative ori-
entations of the flow and the wall. For future application
in complex industrial geometries, and for possible adap-
tation of the RANS and LES zones during the course of

the computation, it is highly desirable to have a more
flexible approach, requiring only quantities at the point
under consideration and no limitation as to the orienta-
tion or geometry of the interface. One such approach is
the ALF method (anisotropic linear forcing), which is a
generalization of the above-mentioned volume forcing ap-
proaches, based on the imposition of a complete target
Reynolds-stress tensor in a volume at the beginning of
the LES region.²¹,²² This approach yields an appreciable
reduction of the adaptation distance downstream of the
inlet. However, it is still restricted to zonal RANS/LES
approaches, and also requires target statistics.

In the present paper, a volume forcing method is
proposed for continuous hybrid RANS-LES approaches.
This method is a modification of the ALF to avoid the
need for target statistics provided by a RANS calcula-
tion, which is at the origin of the restriction of the ALF
to zonal hybrid RANS/LES with an overlap between the
RANS and LES domains. Its essence is to evaluate the
energy transfer between the modeled part of the turbu-
lent spectrum and the resolved part, in order to inject
the necessary energy into the resolved motion via a vol-
ume force. It is first evaluated in comparison with the
ALF in the case of a series of hills with idealized inlet
conditions (i.e., given by a preliminary periodic compu-
tation), for which the ALF has already been successfully
tested.²³ Then it is further extended to a more realistic
case consisting in a succession of inter-hill regions treated
in RANS and LES, in order to investigate the influence
of the new methodology and analyse the generation of
the turbulent fluctuations. Finally, it is also validated
for the case of a simple plane channel flow, which does
not present a shear layer instability and thus is a more
challenging case for the new method.

In the following section, the continuous hybrid
RANS/LES method used herein (HTLES, Hybrid Tem-
poral LES) is briefly introduced and the results obtained
in the periodic hill case are compared with the reference
LES. These results in the periodic configuration will be
used in the remainder of the article as reference results
for the simulations in spatially developing cases. Section
III presents the theoretical arguments and analytical de-
velopments leading to the proposal of the new volume
forcing method. Section IV A is then devoted to the val-
idation of this new approach in a case in LES mode with
stationary inlet conditions, enabling comparison with the
ALF method, and section IV B to the case of a transi-
tion from RANS to LES followed by a transition from
LES to RANS, which prefigures the use of this new ac-
tive approach in practical applications. Section IV C is
dedicated to the extension of the validation to a plane
channel flow and the investigation of the influence of the
length of the RANS-to-LES transition region.
## II. HYBRID RANS/LES FRAMEWORK

The active approach developed in this paper can potentially be applied to any hybrid RANS/LES method, as long as one can identify a parameter that drives the transition from RANS to LES, be it simply the mesh size or any other parameter that can be calculated during the computation. Here, the method used is the Hybrid Temporal LES (HTLES)¹²,²⁴⁻²⁶, which is briefly presented below.

In filtered approaches, any flow variable $f$ can be decomposed into two parts: a filtered part, $\tilde{f}$, and a sub-filtered part, $f''$,

$$f = \tilde{f} + f'', \quad (3)$$

using a general spatio-temporal filtering operator denoted as $\tilde{\cdot}$,

$$\tilde{f}(\mathbf{x}, t) = \int_{\mathcal{D}} \int_{-\infty}^{t} \mathcal{G}(\mathbf{x}, \mathbf{x}', t, t') f(\mathbf{x}', t') d\mathbf{x}' dt' \quad (4)$$

where $\mathcal{D}$ is the domain and $\mathcal{G}$ the filter kernel. This definition encompasses standard LES based on spatial filtering, and temporal LES on which the HTLES approach is built, hence the name Hybrid Temporal LES.²⁴ Whatever the type of filter, the equation of the filtered momentum is:

$$\frac{\partial \tilde{u}_i}{\partial t} + \tilde{u}_j \frac{\partial \tilde{u}_i}{\partial x_j} = -\frac{1}{\rho} \frac{\partial \tilde{p}}{\partial x_i} + \nu \frac{\partial^2 \tilde{u}_i}{\partial x_j \partial x_j} - \frac{\partial \tau_{ij\text{sfs}}}{\partial x_j} \quad (5)$$

where the subfilter scale (sfs) tensor is defined as the generalized second-order moment $\tau_{ij\text{sfs}} = \tilde{u}_i \tilde{u}_j - \tilde{u}_i \tilde{u}_j$.²⁷ This equation is of the same form as the one obtained by using the Reynolds decomposition based on the statistical average $\cdot$, so that the cornerstone of all continuous hybrid RANS/LES methods is to use a model making it possible to transition the closure for this stress tensor from a RANS model to a LES model in some regions of the flow. It can be shown that this is equivalent to making the characteristic size of a spatial filter (resp. a temporal filter) in homogeneous turbulence (resp. statistically stationary turbulence) tend to infinity in Eq. (4).²⁸

### A. Hybrid model

In HTLES,¹² the underlying RANS model is sensitized to the filter width via the time scale $T_m(r)$, dependent on the resolution parameter $r$, which is the modeled-to-total energy ratio

$$r = \frac{k_m}{k}, \quad (6)$$

where $k = k_m + k_r$ is the total turbulent energy, decomposed into a modeled part $k_m = \overline{k_{\text{sfs}}} = \overline{\tau_{ii_{\text{sfs}}}}/2$ and a resolved part $k_r = (\tilde{u}_i \tilde{u}_i - \tilde{u}_i \tilde{u}_i)/2$, were the

overbar denotes statistical averaging, approximated by an exponentially-weighted average²⁹ performed during the computation. Here, following Duffal et al.,¹² this approach is applied to the underlying $k-\omega$-SST RANS model, leading to the definition of a resolution-dependent turbulent time-scale

$$T_m(r) = \frac{r}{\psi(r)} \frac{k_m + c_r k_r}{C_\mu k_m \omega_m^*}, \quad (7)$$

where the hybridization function writes

$$\psi(r) = \frac{\beta_{\omega}}{C_{\mu}\gamma_{\omega} + r(\beta_{\omega} - C_{\mu}\gamma_{\omega})} \quad (8)$$

and the subfilter-scale eddy-viscosity

$$\nu_{\text{sfs}} = \frac{a_1 k_{\text{sfs}}}{\max[a_1 \psi(r) \omega_{\text{sfs}}^*, F_2 \tilde{S}]} \quad (9)$$

The transport equations of the $k-\omega$-SST model then become

$$\begin{align}
\frac{\mathrm{d}k_{\mathrm{sfs}}}{\mathrm{d}t} &= P_k + \frac{\partial}{\partial x_j} \left[ \left( \nu + \frac{\nu_{\mathrm{sfs}}}{\sigma_k} \right) \frac{\partial k_{\mathrm{sfs}}}{\partial x_j} \right] - \frac{k_{\mathrm{sfs}}}{T_m(r)}, \nonumber \\
\frac{\mathrm{d}\omega_{\mathrm{sfs}}^*}{\mathrm{d}t} &= \gamma_{\omega} \frac{1}{\psi(r)} \tilde{S}^2 + \frac{\partial}{\partial x_j} \left[ \left( \nu + \frac{\nu_{\mathrm{sfs}}}{\sigma_{\omega}} \right) \frac{\partial \omega_{\mathrm{sfs}}^*}{\partial x_j} \right] - \beta_{\omega} \omega_{\mathrm{sfs}}^* \nonumber \\
&\quad + (1 - F_1) 2 \frac{1}{\sigma_{\omega 2}} \frac{1}{\psi(r)} \omega_{\mathrm{sfs}}^* \frac{\partial \omega_{\mathrm{sfs}}^*}{\partial x_j} \frac{\partial k_{\mathrm{sfs}}}{\partial x_j}, \tag{10}
\end{align}$$

where $\gamma_\omega$, $\beta_\omega$, $a_1$, $a_2$, $\sigma_k$, $\sigma_\omega$, $F_1$ and $F_2$ are the usual coefficients and functions of the RANS model.³⁰ The variable $\omega_{\text{sfs}}^*$, solution of the second equation, is marked by a superscript star to emphasize that it is not the specific dissipation rate except at the RANS limit: in general, the specific dissipation rate is $\psi(r) \omega_{\text{sfs}}^*$, which appears in Eq. (9). When $r=1$ (RANS mode), the equations of the RANS model are recovered. Note that the coefficient $c_r$ in Eq. (7) is necessary to reach this limit properly (called the *internal consistency constraint*): if resolved structures penetrate from the LES region to the RANS region, their energy is not counted in order to avoid modifying the RANS model, since $c_r$ is defined as

$$\left\{ \begin{array}{ll} c_r = 0 & \text{if } r = 1, \\ c_r = f_s & \text{if } r < 1, \end{array} \right. \quad (11)$$

where $f_s$ is the shielding function defined below. When $r < 1$ (LES mode), $T_m$ is reduced, which leads to an increase of the dissipation $k_{\text{sfs}}/T_m$ in the $k_{\text{sfs}}$ equation, which in turn decreases $k_{\text{sfs}}$ and thus the turbulent viscosity $\nu_{\text{sfs}}$, a mechanism similar to that operating in two-equation DES.

HTLES being based on temporal filtering, assuming an equilibrium Eulerian spectrum, the energy ratio $r$ can be related to the highest frequency resolvable by the numerical scheme

$$\omega_c = \min \left( \frac{\pi}{\tau}, \frac{U_s \pi}{\Delta} \right) \quad (12)$$
through the relation

$$r_K = \frac{1}{\beta_0} \left( \frac{U_s}{\sqrt{k}} \right)^{2/3} \left( \omega_c \frac{k}{\varepsilon} \right)^{-2/3}, \quad (13)$$

where the index K indicates that this evaluation is based on the Kolmogorov assumptions.²⁴ In these relations, τ is the time step, Δ the grid step (estimated herein by the cubic root of the cell volume), β₀ = 0.48 and U_s = U + √(2k/3) is the sweeping velocity.³¹ Note that if the time step τ is smaller than Δ/U_s (or the sweeping CFL number τU_s/Δ is less than unity) everywhere in the domain, which will be the case in practice in the present paper, the cut-off frequency ω_c is related to the grid step Δ and Eq. (13) reduces to the standard relations used, for instance, in PITM,⁶

$$k_m = \int_{\kappa_c}^{\infty} E(\kappa) d\kappa = \frac{\Delta^{2/3} \varepsilon^{2/3}}{\beta_0 \pi^{2/3}}, \quad (14)$$

$$r_K = \frac{1}{\beta_0} \left( \frac{\pi k^{3/2}}{\Delta \varepsilon} \right)^{-2/3}. \quad (15)$$

Following Duffal et al.,¹² a shielding function f_s is introduced in the model as

$$r = (1 - f_s) + f_s \min(1, r_K) \quad (16)$$

$$f_s = 1 - \tanh[\max(\xi_K^{p_1}, \xi_D^{p_2})] \quad (17)$$

This function serves two objectives. The first is to avoid grid-induced separation when the mesh is ambiguous near the wall, by forcing the RANS mode independently of the mesh up to a certain distance from the wall, based on the dimensionless parameter

$$\xi_K = C_1 \frac{(\nu^3 / \varepsilon)^{1/4}}{d_w}. \quad (18)$$

With the values used herein, C₁ = 45 and p₁ = 8, this parameter activates the shielding up to y⁺ ≈ 100, as shown by Duffal et al.¹² The second objective is to mitigate the log-layer mismatch,³² by accounting for the fact that near-wall cells are generally too elongated for LES in the near-wall region and postponing accordingly the RANS-to-LES transition based on the parameter

$$\xi_D = C_2 \frac{\Delta_{\max}}{d_w}, \quad (19)$$

with C₂ = 1.2 and p₂ = 6. Δmax is the longest edge of the local cell. For details about HTLES, its derivation and calibration, the reader is referred to Duffal et al.¹²

## B. Comparison with highly-resolved LES

In section IV, the active hybrid RANS/LES approach will be tested in cases of the spatially developing flows.

As illustrated in Fig. 1, the main issue here is to ensure a fast transition from a RANS behavior to a developed LES behavior (with RANS regions along the walls). In the present section, we will therefore show a short validation of the fully developed behavior of the HTLES model described above, by performing periodic simulations that are compared with the reference refined LES.³³ A more detailed validation is available in Duffal et al.¹² Once validated, this periodic HTLES solution will in turn serve as a reference for the simulations in spatial development presented in the following sections.

For hybrid RANS/LES approaches, a periodic arrangement of 2D hills in a plane channel constitutes a standard benchmark case. Because low-frequency phenomena govern the characteristics of the mean flow, including moving separation and reattachment points, and large-scale structures generated by the detached shear layer have a strong influence on the recirculation region, reproducing the flow statistics is a real challenge for RANS models and hybrid RANS/LES generally bring a significant improvement.

The flow is characterized by the bulk Reynolds number, Re_b = U_bh/ν = 10595, where U_b is the bulk velocity and h is the height of the hill. The computational domain is L_x × L_y × L_z = 9h × 3.035h × 4.5h, and periodicity is applied in the spanwise and streamwise directions. Using a source term in the streamwise momentum equation, a constant mass flow is imposed. The periodic hill flow is studied using a coarse mesh M1 (N_x × N_y × N_z = 80 × 80 × 40) and a fine mesh M2 (N_x × N_y × N_z = 120 × 120 × 60) containing more than three times more cells, which both satisfy Δy⁺ ≈ 1 at the wall. As a reference, the highly resolved LES produced by Breuer et al.³³ with 13.1 million cells is considered. Computations are performed using code-saturne,³⁴ the open-source general purpose CFD software developed by EDF, based on a co-located finite-volume approach that accepts meshes with any type of cells and any type of grid structure. The velocity-pressure system is solved using a SIMPLEC algorithm. A hybrid second-order upwind (SOLU)/centered difference (CDS) scheme is utilized, similar to the one proposed by Travin et al.,³⁵ using c_r from Eq. (11) as a blending coefficient.¹²

The distribution of skin friction on the lower wall highlights the behavior of HTLES and LES. In Fig. 3, it is seen that HTLES reproduces results similar to the reference data, using both grids, although M1 appears slightly too coarse at the top of the hill. Focusing on profiles extracted at different streamwise locations (x/h = 0, 2, 4, 6 and 8) for the mean velocity and total turbulent energy (resolved + modeled), the reference data are well reproduced for both grids as shown in Fig. 4, although the turbulent energy is overestimated. The discrepancies between the two grids are minor considering the large difference in the grid resolution.
<figure><img src="image_4.png" /><figcaption>FIG. 3: Comparison of HTLES and LES for the periodic hill. Skin friction along the lower wall.</figcaption></figure>

### III. ACTIVE HYBRID RANS/LES APPROACH

As mentioned in the introduction, the objective of this paper is to improve the transition between an upstream RANS zone and a downstream LES zone, i.e., to mitigate the grey area problem. Firstly, Fig. 1 illustrates the baseline situation: we perform a channel simulation with periodic constriction, this time not in a periodic domain, but in a spatially developing domain, with inlet and outlet boundaries. The computational domain contains $n$ inter-hill subdomains bounded by $n+1$ hills (the entrance and exit of the domain being located at the tops of the hills, the first and the last hills are in fact half-hills). There will be $n = 3$ subdomains in section IVA and $n = 5$ subdomains in section IVB. Note that the term subdomain is used here for convenience to describe a region between two successive hills, but there is only one computational domain. In order to decouple the numerical issue from the modeling issue and also to avoid problems of connections between areas refined differently by unstructured or non-conforming meshes, all the simulations will be performed with the meshes used in section II B, duplicated in each subdomain. The transition between upstream RANS mode and downstream LES mode will be imposed by prescribing a spatial evolution of the resolution parameter $r$ in the model. On the other hand, in the LES area, $r$ will still be given by Eq. (16). As illustrated by Fig. 1, the resolved content in the LES region takes a long time to develop. In order to remedy this issue, we will now consider a new method to inject fluctuating energy in the resolved motion.

#### A. Theoretical framework

Although the active approach developed herein is applicable to any continuous hybrid RANS/LES model, it is convenient to illustrate its rationale from a filtering formalism. Indeed, the effect of a change of resolution on the energy partition between the resolved motion and the modeled motion is illustrated in Fig. 5, here for a refinement, but the opposite is true for a coarsening. When going from a coarse mesh to a fine mesh, the cutoff wavenumber of the filter (or frequency for a temporal filter) is modified. As shown in Fig. 5a, when the mesh is refined, the resolved energy increases and the modeled energy decreases. Consequently, the energy that appears shaded in Fig. 5a corresponds to an energy transfer from the modeled part to the resolved part. It is this energy transfer that we will try to represent in the active model. The interpretation in physical space is illustrated by Fig. 5b, in which the typical evolution of the turbulent energy along a streamline is plotted. The solution on a coarse mesh, shown in dashed lines, exhibits some partitioning of the total turbulent energy between the modeled $k_m$ and resolved $k_r$ parts. On a finer mesh, the solution, represented in dash-dotted lines, exhibits a different energy partition, with more resolved energy and less modeled energy. In theory, however, this different distribution should not affect the total turbulent energy. If a fluid particle passes along a streamline through a zone of gradual mesh refinement, then the resolved and modeled energies should follow the evolution plotted in solid lines, and the total turbulent energy should remain constant.

This observation is related to the issue, shown by several authors,²⁸,³⁶⁻⁴⁰ that when the equations for the filtered quantities are derived, terms arise due to the variation of the filter size in the equation due to the fact that the filter does not commute with the derivatives. These terms, also called commutation error, correspond to the transfer of energy which is missing in the model. For instance, following Chaouat and Schiestel,³⁹ one can note that if the filter width of the kernel $\mathcal{G}$ defined by Eq. (4) is variable in space, the spatial derivative of a filtered quantity $\tilde{f}$ generates the commutation error

$$ \frac{\partial \tilde{f}}{\partial x_i} - \frac{\partial \tilde{f}}{\partial x_i} = \iiint_{\mathcal{D}} \int_{-\infty}^{t} \frac{\partial \mathcal{G}}{\partial x_i} f \, dx' dy' dz' dt'. \quad (20) $$

Here, a different and simpler way of viewing this problem is used, which has the additional advantage of not referring to a particular formalism and thus emphasizes the generality of the approach. It simply consists in considering, as illustrated by Fig. 5b, that when the mesh is refined, the variables change (typically, the modeled turbulent energy is reduced): it is not a question here of a dependence on the numerical scheme, but of a dependence of the model on the grid step. Thus, it is proposed to express any filtered field $\phi$ as depending not only on
<figure><img src="image_5.png" /><figcaption>FIG. 4: Comparison of HTLES and LES for the periodic hill. (a) Streamwise velocity profiles and (b) total turbulent kinetic energy profiles.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5: Illustration of the transfer of energy from the modeled to the resolved part in hybrid RANS/LES due to a change in resolution. (a) In spectral space. (b) In physical space, with an arbitrary spatial evolution of $k$ ($s$ is the curvilinear coordinate along the streamline).</figcaption></figure>

the coordinates $t$ and $x_i$, but also on the grid step $\Delta$,

$$ \phi(t, x_i, \Delta). \qquad (21) $$

Thus, the variable is no longer a function of four parameters, three of space and one of time, but also of a fifth parameter representing the dependence of the variable on the grid step, illustrated in Fig. 5b. In other terms, the variable is now part of a hyperspace of dimension five. Note that one could, in a completely equivalent way, express $\phi$ not as a function of $\Delta$ but as a function of the cut-off wavenumber $\kappa_c$ or frequency $\omega_c$, or of the energy ratio $r$ as will be shown later on in the paper, but it is

easier to visually represent a dependence on the grid step.

Fig. 6d provides a representation of this hyperspace, in which here, for obvious reasons, the two dimensions $z$ and $t$ are not represented. If we follow a fluid particle (blue trajectory) crossing a domain meshed with a constant grid step $\Delta$, Fig. 6a, the fluid particle evolves at a constant $\Delta$, i.e., remains in an horizontal plane in Fig. 6d. On the other hand, if we follow the fluid particle (red trajectory) in Fig. 6b, where the mesh is suddenly refined, the variables that are dependent on $\Delta$ are now modified by this refinement. Instead of remaining in the horizontal plane of constant $\Delta$ in Fig. 6d, the solution goes from a coarse-mesh solution (horizontal plane at
<figure><img src="image_7.png" /><figcaption>FIG. 6: Fluid particle trajectories in: (a,b,c) physical space; (d) the hyperspace $x_i-\Delta$ (the z-direction is omitted). (a) Uniform mesh; (b) Mesh with a sudden refinement; (c) Mesh with a gradual refinement; (d) View in the (x-y-Δ) space of the three trajectories shown in figures (a), (b) and (c).</figcaption></figure>

large Δ) to a fine-mesh solution (horizontal plane at a smaller Δ). In the region where the mesh is refined, the solution thus rapidly moves downward in the direction of small Δ's in the hyperspace. This is just a different way of representing the change in solution linked to a mesh refinement shown in Fig. 5b. The purple trajectory of Figs.6c and d corresponds to the solution obtained on a gradually refined mesh.

This description in a five-dimensional hyperspace allows us to distinguish, in the evolution of the different physical quantities, what is due to the dynamics of the flow from what is due to the changes in grid step. Indeed, the total variation $d\phi$ of the variable $\phi$ is expressed as

$$d\phi = \frac{\partial\phi}{\partial t}dt + \frac{\partial\phi}{\partial x_i}dx_i + \frac{\partial\phi}{\partial\Delta}d\Delta, \quad (22)$$

such that

$$\begin{aligned} \frac{d\phi}{dt} &= \frac{\partial\phi}{\partial t} + \frac{\partial\phi}{\partial x_i} \frac{dx_i}{dt} + \frac{\partial\phi}{\partial\Delta} \frac{d\Delta}{dt} \\ &= \underbrace{\frac{\partial\phi}{\partial t} + \nabla\phi \cdot U}_{\text{Material derivative}} + \underbrace{\frac{\partial\phi}{\partial\Delta} \frac{d\Delta}{dt}}_{\text{Extra term}} \end{aligned} \quad (23)$$

which is nothing but the summation of the classical material derivative and the extra term due to mesh variations. For the particular case of a constant grid step, this extra term disappears. Hence the material derivative of the filtered quantity can be decomposed in two parts:

$$\frac{d\phi}{dt} = \frac{d\phi}{dt}\Big|_{\Delta} + \frac{d\phi}{dt}\Big|_{t,x_i}. \quad (24)$$

Following the standard notation used in thermodynamics,

$$\left. \frac{d\phi}{dt} \right|_{\Delta} = \frac{\partial\phi}{\partial t} + \nabla\phi \cdot U \quad (25)$$

denotes the derivative of the variable $\phi$ at constant $\Delta$, i.e., the classical material derivative. It is the derivative taken in an horizontal plane in Fig. 6d. Now,

$$\left. \frac{d\phi}{dt} \right|_{t,x_i} = \frac{\partial\phi}{\partial\Delta} \frac{d\Delta}{dt} \quad (26)$$

is the derivative at constant $t$ and $x_i$, i.e., the part due only to variations of $\Delta$, in others words, the derivative taken in the vertical direction in Fig. 6d.

It is easy to interpret this term by looking at the quantities $k$, $k_r$ and $k_m$. Indeed, as illustrated by Fig. 5b, the total turbulent energy is independent of the grid step, and thus

$$\frac{\partial k}{\partial \Delta} = 0. \quad (27)$$

Since, $k = k_r + k_m$,

$$\frac{\partial k_r}{\partial \Delta} = -\frac{\partial k_m}{\partial \Delta}, \quad (28)$$

which reflects the fact that a variation of the grid step leaves the total energy unchanged but changes the energy partition. For example, if we follow a fluid particle as it passes through a mesh refinement zone, the case illustrated by Fig. 5 and Fig. 6b, applying Eq. (23) to $k_m$ and $k_r$ shows that their evolution along the streamline can be separated into two contributions: their mesh-independent dynamics

$$\left. \frac{dk_m}{dt} \right|_{\Delta} = \frac{\partial k_m}{\partial t} + U_k \frac{\partial k_m}{\partial x_k} ; \quad \left. \frac{dk_r}{dt} \right|_{\Delta} = \frac{\partial k_r}{\partial t} + U_k \frac{\partial k_r}{\partial x_k} \quad (29)$$

and their variation due to the grid step refinement

$$\left. \frac{dk_m}{dt} \right|_{t,x_i} = \frac{\partial k_m}{\partial \Delta} \frac{d\Delta}{dt} ; \quad \left. \frac{dk_r}{dt} \right|_{t,x_i} = \frac{\partial k_r}{\partial \Delta} \frac{d\Delta}{dt}. \quad (30)$$
As for the total turbulent energy $k$, it is independent of the grid step, so that

$$ \left. \frac{dk}{dt} \right|_{t,x_i} = 0. \qquad (31) $$

Since $k = k_r + k_m$, from Eqs. (30) and (31), we can conclude that the term

$$ \frac{\partial k_r}{\partial \Delta} \frac{d\Delta}{dt} = - \frac{\partial k_m}{\partial \Delta} \frac{d\Delta}{dt} \qquad (32) $$

represents the rate of energy transfer from the modeled part to the resolved part due the refinement of the grid, i.e., the shaded area in Fig. 5a. It is therefore clear that this term must be accounted for in the evolution equations of these quantities. If numerous models take into account the variations of $k_m$, in order to reduce $\nu_t$, there is no mechanism to add energy in the resolved part. Interestingly, by deriving the subfilter energy equations in the case of a variable filter, Chaouat$^{41}$ showed that different terms appear related to the non-commutativity of the filter with the differential operators, and that the term (32), originating from the material derivative, is the dominant term. Our analysis shows that this additional term arises from the variations of the grid step, independent of any underlying formalism and is therefore valid for any scale-resolving approach in which the partition of energy is related to the grid step (or any other resolution parameter), such as the widely spread DDES or PANS approaches, and even pure LES.

In the framework of the HTLES model, the size of the filter can, possibly, be related to the time step $\tau$, according to Eq. (12). Thus, there is also a sensitivity of the solution to the time step, and the extra term must be replaced by

$$ \left. \frac{d\phi}{dt} \right|_{t,x_i} = \frac{\partial \phi}{\partial \Delta} \frac{d\Delta}{dt} + \frac{\partial \phi}{\partial \tau} \frac{d\tau}{dt} \qquad (33) $$

Note that, in most of the usual simulations, the time step is constant and thus $d\tau/dt = 0$. However, since the material derivative inherently involves the time derivative, Eq. (33) encompasses cases where the time step or the grid step vary in time, such as with automatic mesh refinement methods.

The example of grid turbulence illustrated in Fig. 2 shows that the hybrid model takes indeed into account the variations of the modeled energy as a function of the grid step. In other words, the transport equation of the modeled energy contains, as for all $\Delta$ dependent variables, the two contributions of Eq. (24), and the model is build in order to represent these two terms, i.e., both the natural evolution and the grid sensitivity of the modeled energy. For instance, in two-equation DES approaches, the variation of the modeled energy with the grid step is taken into account by replacing the turbulent length scale $\ell$ by $\min(\ell, C_{DES}\Delta)$ in the dissipation rate. On the other hand, the momentum equation, and consequently

the resolved turbulent energy, remains insensitive to variations of $\Delta$: this is what we will try to model in the next section.

## B. Modeling of the energy transfer due to variations of the resolution

In the situation described in Fig. 5, where the mesh is refined, the HTLES turbulence model is sensitized to the variation of $\Delta$ because the transport equation of the subfilter energy

$$ \frac{dk_{sfs}}{dt} = P_{sfs} + D_{sfs} - \frac{k_{sfs}}{T_m} \qquad (34) $$

involves $T_m$ which depends on the energy ratio $r$, which is linked to $\Delta$ by Eqs. (12) and (13). In other words, the model is intended to represent the total variation of $k_{sfs}$, i.e., both terms of Eq. (24),

$$ \frac{dk_{sfs}}{dt} = \left. \frac{dk_{sfs}}{dt} \right|_{\Delta} + \left. \frac{dk_{sfs}}{dt} \right|_{t,x_i} . \qquad (35) $$

As mentioned in section III A, the second term of Eq. (35) corresponds to the transfer of energy between the modeled part and the resolved part.

However, the resolved momentum equation does not contain a corresponding term. The only way to introduce a source of resolved energy is to add a force to the resolved momentum equation. Following De Laage de Meux et al.,$^{21}$ it is proposed to introduce a fluctuating force of the form

$$ f_i = A_{ij} \tilde{u}_i + B_i, \qquad (36) $$

which linearly depends on the resolved velocity $\tilde{u}_i$. Upon constructing the transport equation for the resolved stress $u'_i u'_j$, where $u'_i = \tilde{u}_i - \tilde{\overline{u}}_i$ is the resolved fluctuating velocity, it can be easily shown that the introduction of this force introduces the production term

$$ T_{ij}^{f} = \overline{f'_{i} u'_{j}} + \overline{f'_{j} u'_{i}} = A_{ik} \overline{u'_{j} u'_{k}} + A_{jk} \overline{u'_{i} u'_{k}}, \qquad (37) $$

where $f'_i$ is the fluctuating part of the force. A mean force $\overline{f}_i$ also appears in the mean momentum equation, which must be equal to zero, because one does not want to influence the mean velocity. Hence,

$$ \overline{f}_i = A_{ij} \overline{u_j} + B_i = 0 \qquad (38) $$

The production term of the resolved turbulent kinetic energy $k_r$ being half the trace of that of the resolved stresses, it can be written as

$$ T^f = A_{ik} \overline{u'_i u'_k}, \qquad (39) $$

which is nothing else than the average work done by the force per unit time. If we want this power to be equal
to the energy decrease in the modeled part due to the variations of $\Delta$, the relation to satisfy is

$$A_{ik} \overline{u'_i u'_k} = - \frac{\partial k_m}{\partial \Delta} \frac{d\Delta}{dt}. \quad (40)$$

Individually, each component must satisfy

$$T_{ij}^f = A_{ik} \overline{u'_j u'_k} + A_{jk} \overline{u'_i u'_k} = - \frac{\partial \tau_{ij}^m}{\partial \Delta} \frac{d\Delta}{dt}, \quad (41)$$

where $\tau_{ij}^m$ is the modeled part of the Reynolds stress, i.e., the Reynolds average of the subfilter-stress tensor, $\tau_{ij}^m = \overline{\tau_{ij} \text{sfs}}$. Anisotropy is simply accounted for by assuming

$$\frac{\partial \tau_{ij}^m}{\partial \Delta} \approx \frac{\tau_{ij}^m}{k_m} \frac{\partial k_m}{\partial \Delta}. \quad (42)$$

Now, the modeled turbulent kinetic energy $k_m$ is given by Eq. (14), such that

$$\frac{1}{k_m} \frac{\partial k_m}{\partial \Delta} = \frac{2}{3\Delta}. \quad (43)$$

The conclusion of this analysis is that, at each point of the domain, the following system of equations must be resolved to determine the values of $A_{ij}$ and $B_i$, which define the fluctuating force $f_i$,

$$A_{ik} \overline{u'_j u'_k} + A_{jk} \overline{u'_i u'_k} = -\frac{2}{3} C \tau_{ij}^m \frac{1}{\Delta} \frac{d\Delta}{dt}, \quad (44a)$$

$$B_i = -A_{ij} \overline{u_j}. \quad (44b)$$

The statistical operator $\overline{\cdot}$ is approximated by exponentially-weighted averaging²⁹ to evaluate the necessary quantities during the computation. Note that a coefficient $C$ has been introduced in the right hand side of Eq. (44a). The analytical derivation above shows that $C = 1$, but in practice, this coefficient introduces a degree of freedom to adjust the intensity of the forcing, as will be discussed later.

Finally, it should be pointed out that the implementation of the method is relatively straightforward. No additional transport equation compared to the standard HTLES method need to be considered. The resolved energy equation has only been introduced into the theoretical derivation to analyze energy transfers, and does not need to be solved. The only modification to the computational code is the introduction of the force $f_i$ given by Eq. (36) as an explicit source term in the resolved momentum equation. This implies a very small increase of the computational cost, due to the solving of the system of Eqs. (44) by a direct method (LU decomposition) in each point of the domain. For the applications considered in the present article, the increase in CPU cost is less than 1%.

### C. Accounting for the shielding functions

The forcing given by Eq. (44) is general and can be applied to any mesh-dependent hybrid RANS/LES approach, since Eq. (43) simply assumes a Kolmogorov spectrum.

However, since many approaches introduce shielding functions, such as HTLES used herein, the dependence of the solution on the grid step $\Delta$ is disrupted near the walls. Therefore, it makes more sense in this case to consider the dependence of the solution on the energy ratio $r$ directly,

$$\phi(t, x_i, r). \quad (45)$$

Upon following the same analytical procedure as in section III B,

$$T^f = - \frac{\partial k_m}{\partial r} \frac{dr}{dt}. \quad (46)$$

Now, $k_m = rk$ and $k$ is independent of the energy ratio $r$, such that

$$T^f = -k \frac{dr}{dt} \quad (47)$$

The system of equations to be solved becomes

$$A_{ik} \overline{u'_j u'_k} + A_{jk} \overline{u'_i u'_k} = -C \tau_{ij}^m \frac{1}{r} \frac{dr}{dt} \quad (48a)$$

$$B_i = -A_{ij} \overline{u_j}, \quad (48b)$$

which is equivalent to Eq. (44) in the absence of shielding functions, since in this case Eq. (15) leads to

$$\frac{1}{r} \frac{dr}{dt} = \frac{2}{3} \frac{1}{\Delta} \frac{d\Delta}{dt}.$$

This expression thus makes it possible to take into account a more complex formulation than a simple dependence on the grid step $\Delta$. It can be extended to any other resolution parameter.

## IV. RESULTS AND DISCUSSION

The volume forcing given by Eq. (48) is a modification of the original ALF (Anisotropic Linear Forcing) proposed by De Laage de Meux et al.²¹ (see below). The ALF method was originally developed to generate fluctuations at the entrance of a LES domain in a zonal hybrid RANS/LES approach, with an overlap of the RANS and LES zones: in the overlap zone, the intensity of the forcing in LES is driven by the statistics given by the RANS. In order to validate that the modified ALF developed in the previous section, which does not use RANS statistics, works similarly to the original ALF, a configuration that makes it possible to compare the two approaches is first considered. In a second step, the validation of the new method, the continuous active hybrid RANS/LES approach based on this modified ALF, will be performed in a case that represents its real objective, namely a continuous RANS/LES hybrid simulation with transition from RANS to LES and from LES to RANS. Finally, the validation will be extended to the case of a plane channel, which is more challenging since it does not exhibit an inflexional velocity profile likely to easily generate structures.
A. Comparison with the original ALF

The configuration to compare the two approaches is the case $n = 3$ (3 inter-hill subdomains) described at the beginning of section III, in which the inlet conditions are independent of time. The aim of the forcing is to achieve rapid spatial development of the LES solution, i.e., to tend as quickly as possible towards the periodic HTLES solution presented in section II B. If we were to use the periodic RANS solution as an inlet condition, the simulation would have to perform a double adaptation in order to tend towards the periodic HTLES solution: correcting the RANS velocity field and generating resolved structures. In order to decouple these two issues, we introduce here as inlet condition at $x/h = 0$ the mean profile resulting from the periodic HTLES solution, so as to focus solely on the issue of generating resolved structures, leaving the coupled problem for later (section IV B). The same mesh (M1) and numerical parameters as for the periodic HTLES are used, and a simulation without forcing is also performed.

In the original ALF, the forcing is designed in order to drive both the mean and fluctuating fields towards a target mean velocity $U_i^\dagger$ and Reynolds-stress tensor $\tau_{ij}^\dagger$ obtained from a RANS computation in an overlapping zone. By identifying the work done by the force as in section IIIB, De Laage de Meux et al.²¹ showed that this can be achieved by using the system of equations

$$A_{ik} \overline{u_j' u_k'} + A_{jk} \overline{u_i' u_k'} = \frac{1}{\tau_r} \left( \tau_{ij}^\dagger - \overline{u_i' u_j'} \right), \quad (49a)$$

$$A_{ij} \overline{u_j} + B_i = \frac{1}{\tau_v} \left( U_i^\dagger - \overline{u_i} \right), \quad (49b)$$

where $\tau_v$ and $\tau_r$ are the relaxation times

$$\tau_r = \max \left( 2\tau, 0.01 \frac{k^\dagger}{\varepsilon^\dagger} \right) \quad \text{and} \quad \tau_v = \frac{5h}{U_b}, \quad (50)$$

with $\tau$ the time step and $U_b$ the bulk velocity. The computation using this original ALF consists of two overlapping zones: the upstream zone, which only contains the flow statistics and provides the target mean velocity $U_i^\dagger$ and Reynolds-stress tensor $\tau_{ij}^\dagger$; and the downstream zone, treated with HTLES in LES mode (with shielding functions). In the overlap, which covers the region $x/h \in [0, 2]$, the ALF is activated. Again, to get rid of the detrimental influence of a RANS model, here the target solution is obtained from averaging in time the periodic HTLES, which can be viewed as an ideal RANS computation, since it indeed provides the targeted solution. Therefore, the computation in the upstream zone corresponds to the periodic HTLES computation described in section III and is in practice performed in advance. The inlet conditions for the downstream computation is also the time-averaged periodic HTLES. Moreover, in order not to bias the comparison, it is ensured that, as for the modified ALF presented in previous section, the force

does not influence the mean flow: the rhs of Eq. (49b) is replaced by 0.

Moreover, in order to make a fair comparison between the different approaches, the following procedure is used for the other spatially-developing HTLES computations, without forcing (standard HTLES) and with forcing (active HTLES). The same inlet conditions are used, with the average velocity profile obtained in periodic HTLES, and also $k_{sfs}$ equal to the total turbulent energy. Just at the inlet, the RANS mode is activated, by imposing $r = 1$ in the model. A transition from RANS to LES is imposed gradually from $x/h = 0$ to $x/h = 2$, by enforcing a modified $r$

$$r_{\text{mod}} = (1 - f) + f r \quad \text{with} \quad f\left(\frac{x}{h}\right) = \frac{1}{2}\frac{x}{h} \\ \text{for} \quad \frac{x}{h} \in [0; 2], \quad (51)$$

such that $r = 1$ at the inlet and gradually reduces to the usual value in HTLES given by Eq. (16). Since the mesh is the same mesh as used for the periodic HTLES and for the HTLES with the ALF, after $x/h = 2$ the LES mode is activated.

From Fig. 7, it can be seen that in the case of the standard HTLES, without forcing, the recovery of the total turbulent kinetic energy is very slow. Just after the inlet ($x/h \simeq 0$), the turbulent energy is virtually equal to the one of the periodic case, since it is entirely modeled ($r = 1$, RANS mode) and the periodic profile of the total turbulent energy is imposed as the boundary condition for the modeled energy $k_{sfs}$. In the RANS-to-LES transition region ($x/h \in [0, 2]$), $k_{sfs}$ is gradually decreased due to Eq. (51), and the growth of resolved turbulent energy is supposed to compensate for this decrease. However, Fig. 7 shows that this compensation mechanism does not work. Indeed, the total turbulent energy remains very small in the first subdomain ($x/h \in [0, 9]$); starts increasing, essentially in the detached shear layer, only in the second subdomain ($x/h \in [9, 18]$); and eventually recovers the expected level at the end of the third subdomain ($x/h > 24$) (it is recalled that there is only one computational domain and that the term subdomain is used to identify the regions between two successive hills). This underestimation is due to the very slow development of the turbulent structures and thus of the resolved energy. This has a great influence on the velocity profiles, visible in Fig. 7: the lack of mixing of momentum by the turbulent structures is at the origin of a too weak entrainment and thus of an underestimation of the backward velocity and an overestimation of the size of the recirculation zone. Consequently, the friction coefficient $C_f$ is very poorly reproduced, until the end of subdomain 2 (Fig. 12a). When the turbulent energy returns to a correct level in the third subdomain, the velocity profiles are much better predicted and $C_f$ tends towards the value obtained in the periodic simulation when approaching the last hill.

In contrast, it can be seen that the ALF is very efficient in generating resolved energy quickly and thus maintain-
<figure><img src="image_8.png" /><figcaption>(a) Subdomain 1</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d) Subdomain 1</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Subdomain 2</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(e) Subdomain 2</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) Subdomain 3</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(f) Subdomain 3</figcaption></figure>

FIG. 7: Profiles of total turbulent kinetic energy (a,b,c) and mean velocity (d,e,f) in the different subdomains. Comparison of periodic HTLES, HTLES without forcing, HTLES with the ALF and Active HTLES, both with a forcing in the region $x/h \in [0, 2]$.
ing a total energy level very close to that obtained in the periodic simulation. This is achieved by forcing at each point in the region $x/h \in [0, 2]$ towards the target turbulent energy according to Eq. (49a). However, Fig. 7a shows that at $x/h = 4$, after the end of the forcing region, a part of the generated fluctuations have been dissipated, and the turbulent energy is slightly underestimated. However, at the end of the first subdomain, the energy, the velocity profile (Fig. 7d) and the friction coefficient (Fig 12a) are very well reproduced. The results tend almost perfectly to the periodic results in the middle of the second subdomain (Figs. 7b, 7e and 12a).

The challenge for the new forcing method, which has no target field but is simply based on the estimate Eq. (48a) of the energy rate to be injected into the resolved part, is to approach the very good results obtained with the ALF. It can be seen in Fig. 7a that in the middle of the RANS-to-LES transition zone, at $x/h = 1$, the total turbulent energy is underestimated: the modeled energy has been reduced by the HTLES model, and the fluctuating energy generated by the forcing is not sufficiently developed yet. The activation of the forcing is related to the transition from RANS to LES driven by Eq. (51), which imposes a decrease of $r$ along the streamlines. On the other hand, at the end of the transition zone ($x/h = 2$), we observe a slight overshoot of the turbulent energy in the detached shear layer and in the recirculation zone, and a stronger overshoot in the core of the flow. The velocity profiles (Fig. 7d) are then slightly less accurate than those given by the ALF method. However, it can be seen that from the end of subdomain 1, the results given by the active HTLES are very close to those given by the ALF and the periodic HTLES, both for the turbulent energy, the velocity profiles and the friction coefficient.

It is to be noted here that the coefficient $C = 1$ is used in Eq. (48a). A study of the influence of this coefficient has of course been performed (not shown here). An increase of $C$ accelerates the production of resolved energy, at the cost of a stronger overshoot at the end of the transition region; a decrease of $C$ suppresses the overshoot but finally degrades the results at the exit of the transition region. Similar to what was observed for the ALF, a part of the resolved energy generated by the forcing is dissipated, and it can be considered that a small overshoot at the exit of the transition region for the active HTLES is favorable, such that the theoretical coefficient $C = 1$ is kept for the following.

In conclusion, both the ALF associated with HTLES in zonal version and the new continuous approach called active HTLES achieve a fast transition from RANS to LES in this configuration. As the ALF was previously validated by De Laage de Meux et al.²¹ and Duffal,²³ this comparison was used to support the fact that it is possible to avoid the need for a target field in the forcing, by replacing the restoring term in Eq. (49a) with the estimate of the rate of energy to be transferred from the modeled part to the resolved part used in Eq. (48a).

This modification is decisive, since it makes it possible to adapt to continuous hybrid approaches the forcing method initially developed for zonal hybrid approaches with an overlapping zone. The following section is therefore devoted to the validation of the new approach in a configuration for which it was developed, with transition from RANS to LES and from LES to RANS.

## B. Validation in a fully realistic configuration

The active approach is now tested in a more realistic scenario, where the simulation transitions from a fully developed RANS solution to LES, and back from LES to RANS at the end of the domain. The geometrical representation of this case $n = 5$ (5 inter-hill subdomains) is shown in Fig. 8. Two subdomains are treated in RANS mode, one at the entrance and one at the exit; and the other three are treated in LES mode (once again, it is recalled that there is only one computational domain and that the term subdomain is used for convenience). The RANS-to-LES transition is imposed in the region $x/h \in [0, 2]$ using Eq. (51), and for the LES-to-RANS transition, the modified $r$ is imposed in a similar way in the region $x/h \in [25, 27]$ using $f(x/h) = (27 - x/h)/2$.

Three cases are compared: (i) a standard HTLES without forcing; (ii) Active HTLES with forcing only activated in case of RANS-to-LES transition; (iii) Active HTLES with forcing both for RANS-to-LES and for LES-to-RANS transition. The results are compared with the periodic RANS solution in RANS regions and with the periodic HTLES transition in LES regions. The computations use the same mesh M1 as in sections IIB and IVA and the same numerical parameters.

Figs. 12b, 9, 10 and 11 show, for the three cases, the comparisons of the friction coefficient, Q-isocontours, turbulent energy profiles and mean velocity profiles, respectively. Note that in all these figures, the reference solution is the periodic RANS solution in subdomains 1 and 5, and the periodic HTLES solution in subdomains 2, 3 and 4. In the first subdomain, treated in RANS mode ($r = 1$), since the inlet conditions correspond to the periodic RANS solution, the three approaches preserve this solution almost perfectly. The small discrepancies are only due to the fact that at the end of the RANS region at $x/h = 0$, instead of having a periodicity condition, the transition zone from RANS to LES starts, which has a small influence on the upstream region.

From the isocontours $Q = 0.2U_b^2/h^2$ plotted in Fig. 9, it is seen that when no force is applied, the resolved structures take a long time to develop after the RANS-to-LES transition located between $x/h = 0$ and $x/h = 2$. In contrast, when the forcing is applied, one can observe a rapid appearance of resolved structures.

This makes a major difference for the turbulent energy in the LES zone shown in Figs. 10b, 10c and 10d. Without forcing, the turbulent energy is very strongly underestimated until the end of subdomain 4, in a man-
<figure><img src="image_14.png" /><figcaption>FIG. 8: Geometrical representation of the case $n = 5$, showing the isocontours of the modeled-to-total energy ratio $r$, indicating the regions solved in RANS ($r = 1$) and in LES ($r < 1$).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 9: Isocontours $Q = 0.2U_b^2/h^2$ colored with the velocity magnitude and, in the background, modeled-to-total turbulent energy ratio $r$, indicating the regions solved in RANS ($r = 1$) and in LES ($r < 1$). Top: Without forcing, Middle: Forcing at the RANS-to-LES transition, Bottom: Forcing at both the RANS-to-LES and the LES-to-RANS transitions.</figcaption></figure>

ner very similar to the case studied in section IV A, in which the essential difference is that at $x/h = 0$, the periodic HTLES solution was used as a boundary condition, whereas here, the values at $x/h = 0$ are advected from the RANS zone. The velocity profiles in Fig. 11 and especially the comparison of $C_f$ in Figs. 12a and 12b show that this modification does not change much the results obtained in the LES zone, which are equally inaccurate.

When the forcing is introduced (regardless of whether it is activated or not in the downstream LES-to-RANS transition zone), the generation of resolved turbulent energy occurs rapidly (Fig. 9). However, a larger overshoot is observed compared to the situation discussed in section IV A (Fig. 10b). This difference arises from the fact that the velocity profiles and turbulent variables at the location $x/h = 0$ are no longer the ideal periodic HTLES profiles, but rather those advected from the upstream RANS domain, which are quite different. Fig. 12 illustrates that in this case, the reproduction of the friction coefficient is less accurate in the region $x/h \in [0, 9]$ compared to the case in section IV A, because the LES mode of HTLES, in addition to developing resolved structures, needs to correct the mean field originating from the RANS region. Fig. 10f displays, at the end of the transition zone ($x/h = 2$), the breakdown of turbulent energy into modeled ($k_m$) and resolved ($k_r$) components. In the absence of forcing, the resolved energy is nearly zero throughout. As mentioned earlier, the reverse flow velocity is very low in the recirculation region (Fig. 11b), causing complete flow relaminarization, with the modeled energy also

going to zero. In contrast, with active HTLES, the resolved energy develops significantly, leading to the aforementioned overshoot. In this case, similar to the situation discussed in section IV A, the investigation of the influence of the coefficient $C$ (not presented here) indicates that it is preferable to generate an overshoot at the end of the transition zone, thereby maintaining $C = 1$. The recovery of the reference (periodic) HTLES solution is slower than in the idealized scenario discussed in the previous section, where periodic profiles were imposed at the inlet. The turbulent energy and velocity profiles still exhibit slight deviations from the periodic solution in the middle of subdomain 4 (Figs. 10d and 11d). However, Fig. 12b shows that the friction coefficient is very satisfactory starting from the end of subdomain 2, i.e., the first subdomain in LES mode.

An interesting observation in this context emerges when the forcing is applied also at the LES-to-RANS transition, i.e., at the end of subdomain 4 ($x/h \in [25, 27]$). Fig. 13 illustrates the behavior of the half-trace of the rhs of Eq. (48a), which represents the average production of resolved energy due to the volume forcing. It reveals that the situations in the RANS-to-LES and LES-to-RANS transition zones are reversed. Specifically, in Fig. 13a, the forcing contributes energy to the resolved motion, except near the lower wall at the location $x/h = 2$, where the flow reenters from the LES region into the near-wall RANS region imposed by the shielding function. Conversely, in Fig. 13b, except for near the wall at the beginning of the hill, one predom-
<figure><img src="image_16.png" /><figcaption>(a) Subdomain 1</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Subdomain 2</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(c) Subdomain 3</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(d) Subdomain 4</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(e) Subdomain 5</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(f) Breakdown of total energy.</figcaption></figure>

FIG. 10: (a) to (e): Profiles of total turbulent kinetic energy in the different subdomains, (f): Breakdown of total turbulent kinetic energy into modeled turbulent kinetic energy and resolved turbulent kinetic energy at the end of the region of transition from RANS to LES in the case with forcing (active HTLES) in the first subdomain.
<figure><img src="image_22.png" /><figcaption>(a) Subdomain 1</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(b) Subdomain 2</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(c) Subdomain 3</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(d) Subdomain 4</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(e) Subdomain 5</figcaption></figure>

<caption>FIG. 11: Profiles of averaged velocity in different subdomains.</caption>
<figure><img src="image_27.png" /><figcaption>FIG. 12: Skin friction coefficient along the lower wall</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>FIG. 13: Profiles of the half-trace of the rhs of the forcing (48) in the transition region (mean production of resolved energy). Case of active HTLES with forcing activated at both RANS-to-LES and LES-to-RANS transitions.</figcaption></figure>

inantly observes a flow transitioning from the LES region ($r < 1$) to the RANS region ($r = 1$), resulting in a positive material derivative of $r$ and consequently the destruction of resolved energy. However, Fig. 9 reveals that resolved structures persist in all cases at the beginning of the RANS domain. Nonetheless, as shown in Fig. 10e, turbulent energy is indeed reduced in the core of the flow due to the negative forcing, but not significantly in the separated shear layer and recirculation zone. Despite this modest effect, it is still visible that the velocity profiles within the recirculation zone are closer to the RANS profiles (Fig. 11e). In Figs. 10e and 11e, the reference periodic HTLES solution is also shown, in order to emphasize that it is not a serious problem if resolved structures survive in the RANS zone after the RANS-to-LES transition,

since the solution then remains close to the solution obtained in LES mode, which is better than the RANS solution. However, it may seem more consistent for the solution to return as rapidly as possible to the RANS solution. To this end, it is possible to increase the effect of the forcing by identifying the RANS-to-LES transition and the LES-to-RANS transition from the sign of $dr/dt$ in Eq. (48a), and using a larger coefficient $C = 1.75$ in the second case instead of $C = 1$. Fig. 14 shows that this increase in coefficient indeed speeds up the transition to the RANS solution: in particular, even though there is still excessive resolved energy, the velocity profile of the RANS solution is recovered as the last hill is approached.

Finally, the robustness of the method to mesh refinement is verified using the two meshes M1 and M2 of sec-
tion II B. It can be seen from Fig. 15, which compares streamwise velocities in subdomain 2 and friction coefficients at the lower wall, that mesh sensitivity is slightly more pronounced than for periodic HTLES, but does not significantly influence the results. In the transition zone at the beginning of subdomain 2, the velocity profiles are virtually identical, and only a moderate difference is observed in the second part of subdomain 2. However, the friction coefficients at the lower wall show no more difference between the M1 and M2 meshes than in the case of the periodic HTLES studied in section II B. These results therefore show that the active method is relatively robust to mesh refinement.

### C. Channel flow

The previous section successfully showed that the proposed active method can efficiently generate resolved structures at the RANS-to-LES transition in the case of the hill flow. This flow, representative of many detached flows, is a relatively favorable case. Indeed, resolved structures are generated, even if insufficiently, in the detached sheared layer, and the forcing has the role of amplifying them. Many other flows, in particular boundary layers, do not exhibit linear instability promoting the generation of resolved structures. It is also desirable in general to be able to generate these structures in a boundary layer upstream of an obstacle, for example, to be able to resolve the flow around the obstacle in LES mode. This is why the present section focuses on the question of the RANS-to-LES transition in a channel flow.

The flow is characterized by the friction Reynolds number, $Re_τ = u_τh/ν = 590$, where $u_τ$ is the friction velocity and $h$ is the half-height of the channel. The periodic solution, which will serve as a reference here, obtained in a domain of size $L_x \times L_y \times L_z = 6.4h \times 2h \times 3.2h$, was studied in detail by Duffal et al.¹² and shown to be satisfactorily close to the DNS data.⁴² The mesh used here, $N_x \times N_y \times N_z = 64 \times 96 \times 64$, satisfying $Δy^+ \simeq 1$ at the wall, was shown to be sufficiently fine for HTLES.¹²

Here, simulations in spatial development are performed, with inlet and outlet conditions. The computational domain and the mesh used for the periodic calculation are simply duplicated ten times, so that the overall size is $L_x \times L_y \times L_z = 64h \times 2h \times 3.2h$, discretized using the mesh given by $N_x \times N_y \times N_z = 640 \times 96 \times 64$. The initial part of the channel, $(x/h \in [0, 6.4])$, is treated in RANS mode by imposing the energy ratio $r = 1$ in the model. The gradual transition from RANS-to-LES occurs in $x/h \in [6.4, 12.8]$, by enforcing a modified energy ratio $r_{\text{mod}}$, as done in previous section:

$$r_{\text{mod}} = f + (1 - f)r \quad \text{with} \quad f\left(\frac{x}{h}\right) = 2 - \frac{1}{6.4} \frac{x}{h} \quad \text{for } x/h \in [6.4, 12.8] \quad (52)$$

and the rest of the channel is treated in LES mode

($f = 0$). Periodic RANS computations are used as inlet boundary conditions.

Fig. 16c shows that without forcing, the total turbulent energy is severely underestimated, all along the channel, since, as can be seen in Fig. 16b, virtually no resolved energy is produced. Note that in the computation without forcing, denoted simply as “HTLES” in the figures, the RANS-to-LES transition also occurs in the region $x/h \in [6.4, 12.8]$ and is imposed by Eq. (52). This leads to completely wrong velocity profiles (Fig. 17) and consequently, a much underestimated friction coefficient in Fig. 18. On the other hand, when the forcing is applied, the resolved turbulent energy shown in Fig. 16b rapidly tends towards the correct level, at the end of the forcing region, located at $x/h = 12.8$. It is seen in Fig. 19 that despite the absence of inflectional velocity profile, the active approach is able to rapidly generate resolved structures (note that only the first half of the computational domain ($x/h \in [0, 32]$) is shown in this figure). However, at the end of the forcing region ($x/h = 12.8$), it is apparent that a part of the structures generated by the forcing rapidly disappear in the central region of the channel. In Fig. 17, an improved velocity profile is observed compared to the case without forcing. Although far from perfect, the prediction of the friction coefficient is significantly improved (Fig. 18), and tends to the periodic value at the end of the domain.

It is interesting to look in more detail at Fig. 16, which shows the modeled and resolved contributions to the total turbulent energy, in the transition zone ($x/h \in [6.4, 12.8]$). In the case without forcing, the gradual decrease of $k_m$ is observed (Fig. 16a), which tends towards the profile of the periodic solution at the exit of the transition zone. This decrease is not compensated by an increase in $k_r$, which remains virtually zero (Fig. 16b), which causes the strong underestimate of $k = k_m + k_r$ (Fig. 16). When the forcing is activated, $k_r$ tends towards the periodic profile at the end of the transition zone. On the other hand, $k_m$ no longer decreases as quickly as without forcing, which can be understood by looking at Eq. (7): the reduction of $r$ tends to decrease the value of the time scale $T_m$, but this effect is partially compensated by the appearance of modeled energy $k_r$. Quite paradoxically, when the forcing is introduced, the overshoot observed in Fig. 16c in the total energy is not due to the rapid increase in resolved energy but due to too slow a decrease in modeled energy.

So far in this section, the transition zone has been arbitrarily set to the interval $x/h \in [6.4, 12.8]$, which appears relatively long for practical applications. It is then legitimate to wonder if this length can be reduced: the intensity of the forcing (i.e., the work of the force) being controlled by the material derivative of $r$ in the rhs of Eq. (48a), a shortening of the transition zone induces an intensification of the force, so that at the end of the transition zone, the right amount of energy has been transferred between the modeled part and the resolved part. Two additional simulations were therefore carried out, in
<figure><img src="image_29.png" /><figcaption>FIG. 14: Improvement of the LES-to-RANS transition (subdomain 5) by using different coefficients in the two transition regions.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>FIG. 15: Comparison of the solutions obtained with two different grids.</figcaption></figure>

which the transition zone starts at the same place but is shortened by a factor of two ($x/h \in [6.4, 9.6]$) and by a factor of four ($x/h \in [6.4, 8]$). Figs. 17 and 18 show that this actually has a very detrimental effect on the mean flow prediction. To understand this issue, the evolution of the modeled and resolved contributions in Fig. 16 must be examined. Fig. 16a shows that, as expected, when the RANS-to-LES transition is more abrupt, the modeled energy $k_m$ decreases more quickly. Fig. 16b indicates that indeed, the forcing being more intense, the resolved energy is produced more quickly, but unfortunately it does not survive at the end of the transition zone, when the forcing stops. For the case of a transition zone of length four times shorter, the first profile, located at $x/h = 6.45$, i.e., in the transition zone, shows that the resolved energy

grows very quickly. However, the second profile, located at $x/h = 8.35$, just after the end of the transition zone, shows that this resolved energy is very quickly dissipated. The same phenomenon is observed for the case of a transition zone reduced only by a factor of two, simply shifted further downstream. Even in the case of the longest transition zone, some of the resolved structures are dissipated in the central region, as can be seen in Fig. 19, but enough remain for them to be regenerated quickly.

It can be concluded from this study that the forcing makes it possible to generate resolved structures efficiently, even in the absence of an inflectional velocity profile allowing structures to be generated naturally. Even in a channel, when the turbulent viscosity is reduced, very
<figure><img src="image_31.png" /><figcaption>FIG. 16: Channel flow: Profiles of turbulent energy (modeled, resolved and total) in the region $x/h \in [6.4; 12.8]$ extracted at $x/h = 6.45; 8.35; 10.35$ and $12.65$.</figcaption></figure>

small fluctuations appear which the forcing manages to amplify during the iterations, until reaching the expected level. Indeed, the amplification of disturbances must not be considered only in the direction of the flow, but also over the iterations: at each point, the small initial fluctuations are amplified by the forcing at each iteration, until reaching the equilibrium solution. However, the study also shows that the RANS-to-LES transition should not be too rapid, otherwise the generated fluctuations are dissipated when leaving the forcing zone. We can compare this remark to the conclusions of the article by Druault et al.:⁴³ when the forcing is too intense, the fluctuations generated no longer resemble turbulent structures, and are then dissipated when the forcing stops. In other words, the introduction of the forcing modifies the

equations of motion, and the structures generated are no longer solutions to the equations when the forcing stops. It is therefore appropriate to have a transition zone that is not too short so that the forcing term is not too strong compared to the other terms of the resolved momentum equation at the exit from the transition zone. In a general case, the question of the length of the minimum transition zone remains an open question.

V. CONCLUSIONS

In hybrid RANS/LES methods, when the flow transitions from a RANS zone to a LES zone, there is a shift from fully modeled to partially resolved quantities. This
<figure><img src="image_32.png" /><figcaption>FIG. 17: Channel flow: Profiles of the streamwise mean velocity extracted every 6.4h from x = 3.5h. In the RANS zone (x/h ∈ [0; 6.4]), the reference is the periodic RANS solution; in the LES zone (x/h > 6.4), it is the periodic HTLES solution.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>FIG. 18: Channel flow: skin friction coefficient along the lower wall.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>FIG. 19: Channel flow: Q-isocontours obtained for a RANS-to-LES transition in the region x/h ∈ [6.4; 12.8]. Note that only the first half of the computational domain is shown (x/h ∈ [0; 32]).</figcaption></figure>

transition necessitates compensating for the reduction
in modeled energy with an increase in resolved energy.
Unfortunately, the growth of resolved fluctuations relies
solely on natural instabilities and is not sufficiently rapid
to counterbalance the stress decay imposed by the model.
A simple analysis in spectral space reveals that the miss-
ing element is the transfer of energy from the modeled
component to the resolved component, since the momen-
tum equation remains unaffected by variations in grid
resolution.

To address this issue, a method is proposed to model
the terms resulting from grid step variations, which are
responsible for this energy transfer but are commonly dis-

regarded, through the use of a fluctuating volume forcing.
Although the existence of these terms, due to commuta-
tion errors between the filter and the differential opera-
tors, has been noted by various authors using filter-based
formalisms, the analytical approach used here is simpler
and applies to any hybrid RANS/LES model, whether it
explicitly refers to a particular formalism or not, as long
as the turbulent energy partition between modeled and
resolved scales depends on an identifiable parameter. It
is demonstrated that the intensity of this force can be
determined by equating the work done by the force per
unit time to the decrease in modeled energy caused by
grid step variations.
This novel approach, referred to as *active* hybrid RANS/LES, is validated in a channel flow configuration with periodic constriction. In this setup, the flow transitions from RANS to LES, and the purpose is to rapidly reach a fully developed periodic solution. The results indicate that in the absence of forcing, the development of resolved structures is slow, and it takes three times the inter-hill distance to approach the periodic solution. In contrast, when the forcing is introduced, the solution becomes immediately physically acceptable, with resolved structures developing rapidly. After just one inter-hill distance, the solution closely resembles the fully developed state. Furthermore, the forcing is capable of handling the LES-to-RANS transition by changing the energy transfer direction. This reduces the resolved energy when the flow reenters a RANS region.

The case of a plane channel finally shows that the forcing makes it possible to generate resolved structures even in a situation without a detached sheared layer which naturally creates fluctuations. The amplification during the iterations of fluctuations initially of very small amplitude makes it possible to reach a final equilibrium state with the desired level of resolved energy. However, the study of the influence of the spatial length of the transition zone shows that it is necessary not to introduce too brutal a forcing into the equations, otherwise it would generate structures that are not sufficiently realistic and which are dissipated at the exit from the forcing zone.

The active approach presented in this study has the potential to be applied to any hybrid RANS/LES method, provided a parameter driving the RANS-to-LES transition can be identified. It solves the problem of the appearance of resolved structures at the RANS-to-LES transition in the flow direction, but could also improve the normalwise transition and avoid the log-layer mismatch, in the spirit of the work of Hamba,⁴⁴ by considering, in addition to the material derivative (Eq. 23),

the diffusion terms. Importantly, this method is highly versatile as it does not require information about the interface between RANS and LES regions; it relies solely on local quantities, particularly the material derivative of the resolution parameter. This is a significant advantage over approaches that introduce synthetic vortices, which require scales evaluated from a RANS calculation or from the upstream RANS zone. Since it is based on a first principle, conservation of energy, it is to be assumed that the method is not fundamentally dependent on flow characteristics, and can be extended to other geometries and higher Reynolds numbers. The promising results obtained pave the way for a flexible utilization of continuous hybrid approaches in industrial configurations.

The authors have no conflicts to disclose.

## ACKNOWLEDGMENTS

This study was funded by E2S UPPA (ANR-16-IDEX-0002) in the framework of the ASTURIES project (grant E2S-20-ScientificChallenges-01).

The computing resources were provided by MCIA (Mésocentre de Calcul Intensif Aquitain) and GENCI-IDRIS (Grants 2022-A0122A10980 and 2023-A0142A10980).

**Mahitosh Mehta:** Methodology; Software; Investigation; Visualization; Writing/Original Draft Preparation. **Rémi Manceau:** Conceptualization; Methodology; Supervision; Writing/Review & Editing; Funding Acquisition; Project Administration. **Vladimir Dufal:** Methodology; Software; Investigation; Visualization; Writing/Review & Editing. **Benoît de Laage de Meux:** Methodology; Supervision; Writing/Review & Editing; Funding Acquisition; Project Administration.

<footnote n="1">¹ P. R. Spalart. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. In *Proceedings of first AFOSR international conference on DNS/LES*. Greyden Press, 1997.</footnote>

<footnote n="2">² P. R. Spalart, S. Deck, M. L. Shur, K. D. Squires, M. K. Strelets, and A. Travin. A new version of detached-eddy simulation, resistant to ambiguous grid densities. *Theor. Comput. Fluid Dyn.*, 20(3):181–195, 2006.</footnote>

<footnote n="3">³ M. L. Shur, P. R. Spalart, M. K. Strelets, and A. K. Travin. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities. *Int. J. Heat Fluid Fl.*, 29(6):1638–1649, 2008.</footnote>

<footnote n="4">⁴ S. S. Girimaji. Partially-Averaged Navier-Stokes Model for Turbulence: A Reynolds-Averaged Navier-Stokes to Direct Numerical Simulation Bridging Method. *J. Appl. Mech.*, 73(3):413–421, 2005.</footnote>

<footnote n="5">⁵ F. Menter and Y. Egorov. A scale adaptive simulation model using two-equation models. In *43rd AIAA aerospace sciences meeting and exhibit*, page 1095, 2005.</footnote>

<footnote n="6">⁶ B. Chaouat and R. Schiestel. A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows. *Phys. Fluids*, 17(6): 065106, 2005.</footnote>

<footnote n="7">⁷ S. Heinz, R. Mokhtarpoor, and M. Stoellinger. Theory-based Reynolds-averaged Navier-Stokes equations with large eddy simulation capability for separated turbulent flow simulations. *Phys. Fluids*, 32(6), 2020.</footnote>

<footnote n="8">⁸ B. Chaouat. The State of the Art of Hybrid RANS/LES Modeling for the Simulation of Turbulent Flows. *Flow Turbul. Combust.*, 99(2):279–327, 2017.</footnote>

<footnote n="9">⁹ S. Heinz. A review of hybrid RANS-LES methods for turbulent flows: Concepts and applications. *Prog. Aerosp. Sci.*, 114:100597, 2020.</footnote>

<footnote n="10">¹⁰ P. Spalart. Detached-eddy simulation. *Annu. Rev. Fluid Mech.*, 41:181–202, 2009.</footnote>

<footnote n="11">¹¹ C. Mockett, W. Haase, and D. Schwamborn. Go4hybrid: Grey area mitigation for hybrid RANS-LES methods. *Notes on Numerical Fluid Mechanics and Multidisciplinary*</footnote>
Design, 134, 2018.

12 V. Duffal, B. de Laage de Meux, and R. Manceau. Development and validation of a new formulation of hybrid temporal large eddy simulation. *Flow Turbul. Combust.*, 108(1):1–42, 2022.

13 G. Comte-Bellot and S. Corrsin. The use of a contraction to improve the isotropy of grid-generated turbulence. *J. Fluid Mech.*, 25:657–682, 1966.

14 A. Keating, G. De Prisco, and U. Piomelli. Interface conditions for hybrid RANS/LES calculations. *Int. J. Heat Fluid Fl.*, 27(5):777–788, 2006. ISSN 0142-727X. Special issue of the 6th International Symposium on Engineering Turbulence Modelling and Measurements – ETMM6.

15 M. L. Shur, P. R. Spalart, M. K. Strelets, and A. K. Travin. Synthetic turbulence generators for RANS-LES interfaces in zonal simulations of aerodynamic and aeroacoustic problems. *Flow Turbul. Combust.*, 93(1):63–92, 2014.

16 J. Janin, F. Duval, C. Friess, and P. Sagaut. A new linear forcing method for isotropic turbulence with controlled integral length scale. *Phys. Fluids*, 33:045127, 2021.

17 N. Jarrin, S. Benhamadouche, D. Laurence, and R. Prosser. A synthetic-eddy-method for generating inflow conditions for large-eddy simulations. *Int. J. Heat Fluid Fl.*, 27(4):585–593, 2006.

18 S. W. Haering, T. A. Oliver, and R. D. Moser. Active model split hybrid RANS/LES. *Phys. Rev. Fluids*, 7: 014603, 2022.

19 A. Spille-Kohoff and H.-J. Kaltenbach. Generation of turbulent inflow data with a prescribed shear-stress profile. In C. Liu, L. Sakell, and T. Beutner, editors, *DNS/LES Progress and Challenges, Proc. third AFOSR International Conference on DNS/LES Arlington, TX, USA*. Greyden Press, Columbus, OH, USA, 2001.

20 R. Laraufie, S. Deck, and P. Sagaut. A dynamic forcing method for unsteady turbulent inflow conditions. *J. Comput. Phys.*, 230(23):8647–8663, 2011.

21 B. De Laage de Meux, B. Audebert, R. Manceau, and R. Perrin. Anisotropic linear forcing for synthetic turbulence generation in large eddy simulation and hybrid RANS/LES modeling. *Phys. Fluids*, 27(3):035115, 2015.

22 D. Kempf and C.-D. Munz. Zonal direct-hybrid aeroacoustic simulation of trailing edge noise using a high-order discontinuous Galerkin spectral element method. *Acta Acustica*, 6, 2022.

23 V. Duffal. *Développement d'un modèle hybride RANS-LES pour l'étude des efforts instationnaires en paroi*. PhD thesis, Université de Pau et des Pays de l'Adour, 2020.

24 R. Manceau. Progress in Hybrid Temporal LES (invited keynote paper). In Y. Hoarau, S.-H. Peng, D. Schwamborn, and A. Revell, editors, *Papers contributed to the 6th Symp. Hybrid RANS-LES Methods, 26–28 September 2016, Strasbourg, France*, volume 137 of *Notes on Numerical Fluid Mechanics and Multidisciplinary Design*, pages 9–25. Springer, 2018.

25 M. Mays, S. Lardeau, and S. Laizet. Capturing the drag crisis in the flow around a smooth cylinder using a hybrid RANS-LES model on coarse meshes. *Int. J. Heat Fluid Fl.*, 103:109203, 2023.

26 G. Hyde-Linaker, P. Hall Barrientos, S. Stoumpos, D. Kingsmore, and A. Kazakidi. Patient-specific computational haemodynamics associated with the surgical creation of an arteriovenous fistula. *Med. Eng. Phys.*, page 103814, 2022.

27 M. Germano. Turbulence: the filtering approach. *J. Fluid Mech.*, 238:325–336, 1992.

28 A. Fadai-Ghotbi, C. Friess, R. Manceau, T. B. Gatski, and J. Borée. Temporal filtering: A consistent formalism for seamless hybrid RANS-LES modeling in inhomogeneous turbulence. *Int. J. Heat Fluid Fl.*, 31(3):378–389, 2010.

29 C. D. Pruett, T. B. Gatski, C. E. Grosch, and W. D. Thacker. The temporally filtered Navier-Stokes equations: Properties of the residual stress. *Phys. Fluids*, 15(8):2127–2140, 2003.

30 F. R. Menter. Two-equation eddy-viscosity turbulence models for engineering applications. *AIAA J.*, 32(8):1598–1605, 1994.

31 H. Tennekes. Eulerian and lagrangian time microscales in isotropic turbulence. *J. Fluid Mech.*, 67(3):561–567, 1975.

32 N. Nikitin, F. Nicoud, B. Wasistho, K. Squires, and P. R. Spalart. An approach to wall modeling in large-eddy simulations. *Phys. Fluids*, 12(7):1629–1632, 2000.

33 M. Breuer, N. Peller, C. Rapp, and M. Manhart. Flow over periodic hills—numerical and experimental study in a wide range of reynolds numbers. *Comput. Fluids*, 38(2): 433–457, 2009.

34 F. Archambeau, N. Méchitoua, and M. Sakiz. Code Saturne: A Finite Volume Code for the computation of turbulent incompressible flows - Industrial Applications. *Int. J. on Finite Volume*, 1(1), 2004.

35 A. Travin, M. Shur, M. Strelets, and P. Spalart. *Advances in LES of Complex Flows*, chapter Physical and Numerical Upgrades in the Detached-Eddy Simulation of Complex Turbulent Flows, pages 239–254. New York: Kluwer Acad., 2002.

36 S. Ghosal and P. Moin. The Basic Equations for the Large Eddy Simulation of Turbulent Flows in Complex Geometry. *J. Comput. Phys.*, 118:24–37, 1995.

37 M. Germano. Properties of the hybrid RANS/LES filter. *Theor. Comput. Fluid Dyn.*, 17(4):225–231, 2004.

38 S. Wallin and S. Girimaji. Commutation error mitigation in variable-resolution pans closure: Proof of concept in decaying isotropic turbulence. In *6th AIAA Theoretical Fluid Mechanics Conference*, page 3105, 2011.

39 B. Chaouat and R. Schiestel. Partially integrated transport modeling method for turbulence simulation with variable filters. *Phys. Fluids*, 25, 2013.

40 L. Davidson. Zonal pans: evaluation of different treatments of the RANS-LES interface. *J. Turbul.*, 17(3):274–307, 2016.

41 B. Chaouat. Commutation errors in pitm simulation. *Int. J. Heat Fluid Fl.*, 67:138–154, 2017.

42 R. D. Moser, J. Kim, and N. N. Mansour. Direct numerical simulation of turbulent channel flow up to $Re_τ = 590$. *Phys. Fluids*, 11(4):943–945, 1999.

43 P. Druault, S. Lardeau, J.-P. Bonnet, F. Coiffet, J. Delville, E. Lamballais, J. Largeau, and L. Perret. Generation of Three-Dimensional Turbulent Inlet Conditions for Large-Eddy Simulation. *AIAA J.*, 42(3):447–456, 2004.

44 F. Hamba. Log-layer mismatch and commutation error in hybrid RANS/LES simulation of channel flow. *Int. J. Heat Fluid Fl.*, 30(1):20–31, 2009.