Development and validation of a hybrid temporal LES model in the perspective of applications to internal combustion engines

CFD simulation tools are increasingly used nowadays to design more fuel-efficient and clean Internal Combustion Engines (ICE). Within this framework, there is a need to benefit from a turbulence model which offers the best compromise between prediction capabilities and computational cost. The Hybrid Temporal LES (HTLES) approach is here retained within the perspective of an application to ICE configurations. HTLES is a hybrid Reynolds-Averaged Navier Stokes/Large Eddy Simulation (RANS/LES) model based on a solid theoretical framework using temporal filtering. The concept is to model the near-wall region in RANS and to solve the turbulent structures in the core region if the temporal and spatial resolutions are fine enough. In this study, a dedicated sub-model called Elliptic Shielding (ES) is added to HTLES in order to ensure RANS in the near-wall region, regardless of the mesh resolution. A modification of the computation of the total kinetic energy and the dissipation rate was introduced as first adaptions of HTLES towards non-stationary ICE configurations. HTLES is a recent approach, which has not been validated in a wide range of applications. The present study intends to further validate HTLES implemented in CONVERGE code by examining three stationary test cases. The first validation consists of the periodic hill case, which is a standard benchmark case to assess hybrid turbulence models. Then, in order to come closer to real ICE simulations, i.e., with larger Reynolds numbers and coarser near-wall resolutions, the method is validated in the case of a channel flow using wall functions and in the steady flow rig case consisting in an open valve at a fixed lift. HTLES results are compared to RANS k-x SST and wall-modeled LES r simulations performed with the same grid and the same temporal resolution. Unlike RANS, satisfactory reproduction of the flow recirculation has been observed with HTLES in the case of periodic hills. The channel flow configuration has underlined the capability of HTLES to predict the wall friction properly. The steady flow rig shows that HTLES combines advantages of RANS and LES in one simulation. On the one hand, HTLES yields mean and rms velocities as accurate as LES since the scale-resolving simulation is triggered in the core region. On the other hand, hybrid RANS/LES at the wall provides accurate pressure drop in contrast with LES performed on the same mesh. Future work will be dedicated to the extension of HTLES to non-stationary flows with moving walls in order to be able to tackle realistic ICE flow configurations.


Introduction
The fight against global warming and the reduction of atmospheric pollution are meaningful environmental concerns in industrial societies and involve pursuing the continuous effort of reducing fuel consumption and control of emissions induced by Internal Combustion Engines (ICE). For these reasons, pollutant standards and CO 2 emission targets are becoming more and more restrictive. In addition to these environmental considerations, the cost and time for the development of new products is a main concern due to the very competitive environment of the automotive industry. In this framework, there is a strong necessity to benefit from simulation tools which offer a good compromise between prediction capabilities and cost. Among all the physical phenomena occurring in internal combustion engines, the motion of the internal flow is one of the key factors as it is known to directly affect the global efficiency of the engine and its emissions. Turbulent flows encountered in ICE are wall-bounded and characterized by high Reynolds numbers of about 10 5 . This induces turbulent scales from an order of magnitude of 10 cm related to the bore diameter down to 10 lm for the Kolmogorov turbulent scale. Direct Numerical Simulation (DNS) resolves the entire scale of motion and the computational cost increases as Re 4 s [1], where Re s denotes the Reynolds number based on the friction velocity u s , and this remains far beyond the capabilities of computers, both at present and in the near future. The reduction of the computational cost by introducing a model for turbulence is thus a necessity, and two main approaches are today commonly employed.
On one hand, Reynolds-Averaged Navier Stokes (RANS), which can also be denoted as URANS for unsteady RANS, in non-stationary flows, is the most used approach and is fully industrialized. Its small CPU requirements, as well as the availability of a large number of models explain its popularity. In the frame of ICE, RANS only provides phase-averaged quantities, which is sufficient for assessing stable and steady operating points. Comparisons between experimental and numerical results exhibit good agreement either for spark ignited engines [2] or for compression ignited engines [3] by using for instance the Extended Coherent Flame Model (ECFM) as combustion model [4]. Nevertheless, due to its averaged nature, RANS can hardly capture fully unsteady phenomena such as Cycle-to-Cycle Variations (CCV), even if recent work suggests that RANS could be relevant for a qualitative assessment [5]. On the other hand, Large Eddy Simulation (LES), based on a spatial filtering of the equations of motion, resolves large scales explicitly and models the small scales. This approach is fully adapted to predict phenomena such as the occurrence of knock or CCV [6]. However, LES requires that the grid size must be sufficient to solve a large amount of the turbulent energy which demands high computational cost, especially near the wall where turbulent structures are rather small [7]. Wall-modeling partially solves this issue by significantly decreasing the computational cost of LES since the grid resolution increases weakly as ln(Re s ) instead of Re 2 s for wall-resolved LES which still remains unfeasible for industrial applications [7]. Setting apart the gain of computational cost, applications of wallmodeled LES in ICE flows have shown limited accuracy for the wall friction assessment [8], potentially due to the inconsistency between resolved LES fields, which relies on the spatial filtering and wall functions that provide statistical averaged quantities rather compatible with the RANS framework.
An alternative to wall-modeled LES consists in achieving the coexistence in the same computation of RANS and LES models. The approach called hybrid RANS/LES modeling has emerged very early since the concept was proposed by Schumann [9], and since then two different families of hybrid methods can be identified according to Sagaut et al.'s classification [10]: -Zonal hybrid methods rely on a discontinuous treatment at the RANS-LES interface. The RANS and LES regions are predefined by the user, and standard RANS and LES models are used in each region. The main difficulty consists in addressing the management of the interfaces: the resolved content has to be reconstructed explicitly at the inlet of a LES region to take into account the lack of resolved fluctuations in the RANS region. -Global or seamless methods are based on one set of equations and continuous treatment at the interface between RANS and LES zones. In particular, this approach usually treats the near-wall region with RANS and operates in LES mode where the grid resolution is suitable by reducing the modeled-stresses or the turbulent viscosity of the original RANS model.
Despite the development of a large number of hybrid methods, few of them were applied to ICE flows: Buhl et al. [11] performed a comparative study of several hybrid approaches, such as Scale Adaptive Simulation (SAS) and Detached Eddy Simulation (DES), with LES sub-grid models. Moreover, Partially-Averaged Navier-Stokes (PANS) has been recently applied to moving IC engine [12].
Despite their attractiveness, many hybrid models suffer from empirical bases, and do not tackle the issue of mutual combination in the same computational domain of averaged RANS and instantaneous filtered LES quantities. In this framework, the recent seamless model called Hybrid Temporal LES (HTLES) method [13] is preferred because of its solid theoretical formalism. Indeed, the multi-scale approach is used to model the unresolved structures, and the introduction of the temporal filtering [14,15] solves the inconsistency for stationary flows that appears between the scale-resolving and RANS regions in conventional hybrid methods. The first part of this article is dedicated to recalling the main hypotheses and set of equations of HTLES, and to introducing the Elliptic Shielding (ES) sub-model which intends to enforce RANS near the wall, regardless of the mesh resolution. Further modifications of the computation of the total turbulent kinetic energy and the dissipation rate are added to HTLES implemented in CONVERGE as a first adaption towards non-stationary flows applications.
Regarding the validation of the HTLES model, Manceau [13] studied the flow around a square-sectioned cylinder and showed that HTLES provides results close to experimental measurements with a substantial computational cost reduction compared to LES. The present paper intends to complete the validation of HTLES and to assess the validity of the proposed modifications of the model. For that purpose, three stationary test cases for which reference data are available are examined. Firstly, the developed HTLES is validated on a flow over periodic hills, which is a common benchmark test case used to assess hybrid turbulence models [16,17]. Nevertheless, this test case significantly differs from the targeted industrial applications since the separation due to adverse pressure gradient on a smooth wall requires the resolution of the boundary layer down to the wall, which is not a common practice for ICE simulations. In order to evaluate the predictive capabilities of the method in the case of a simulation representing the near-wall region by wall functions, two additional cases are computed. The channel flow case makes the assessment of the capability to properly predict the wall friction possible. From an automotive engineering standpoint, the correct prediction of friction is of prior importance as it is directly related to the mass of trapped air, which is strongly coupled to the output power of the engine. The third test case is related to a steady flow rig dedicated to studying the flow around an open valve with a fixed lift whose dimensions are very close to actual engines.
2 Hybrid temporal large eddy simulation

Solution to the inconsistency issue
Whatever the concept applied for making the model sensitive to the grid resolution, many models define the scale-resolving mode empirically, which can lead to wrong dissipation levels in the LES region as observed with the VLES model [10]. Moreover, few of them address the issue of combining in the same computational domain RANS and LES quantities, which are of different nature [18].
The Temporal Partially Integrated Transport Model (TPITM) [19] addresses the above-mentioned issues. On the one hand, TPITM converts any RANS model into a scale-resolving model [20] in a way similar to the Partially Integrated Transport Model (PITM) to define the subfiltered scales derived from the consistent multi-scale approach [21][22][23]. This method uses a spectral cut-off to define the unresolved scales, and makes the simulation of turbulent flows possible on relatively coarse grids when the cut-off wave number can be located in the energycontaining range and as far as the grid-size is suitable to correctly describe the resolved flow. Moreover, TPITM eliminates the inconsistency in the transition zone between LES and RANS in the framework of statistically stationary flows by using a temporal formalism for LES [14], since any temporally filtered quantity tends to the statistically averaged quantity in the limit of an infinite filter width [18].
TPITM transforms any RANS model into a scaleresolving model by reducing the modeled stresses or the turbulent viscosity. In the present work, the sub-filter stress s ij is modeled using the Boussinesq approximation [24] where m sfs is the sub-filter viscosity and S ij the resolved strain rate tensor. The sub-filter turbulent kinetic energy, k sfs , represents the energy of eddies with higher frequencies than the cut-off frequency imposed by the model. The TPITM model reduces the sub-filter viscosity of the RANS model by making the dissipation equation sensitive to the temporal filter width via a modified C e2 coefficient. However, applications of this model have shown difficulties in sustaining the correct level of resolved energy, since the energy partition is indirectly controlled through the e-equation and not directly in the k-equation as it is the case for some hybrid approaches such as DES [25].
As a solution, Tran et al. [24] proposed the HTLES model which transfers the dependence of TPITM on the temporal filter width from the dissipation equation to the sub-filter energy equation to sustain the correct level of resolved energy [26]. To this end, equivalence between TPITM and HTLES is ensured by assuming that both models provide the same level of sub-filter energy for the same filter width and tend to the same RANS model when the filter width approaches infinity. As a consequence, HTLES avoids the main weakness of TPITM of setting the energy partition control on the dissipation equation while keeping all the theoretical advantages of TPITM. The set of equations of HTLES is provided in the following section.

Sub-filter viscosity evaluation
HTLES is based here on the k-x SST model [13,24], in which the turbulence viscosity expresses as where S is the strain rate, F 2 ¼ max ffiffiffiffi ffi k sfs p 0:09xy ; 500m In the same way as the two-equation DES, the resolved scales are controlled by an appropriate modification of the dissipation term in the k-equation, but this modification is based here on a time scale rather than a length scale where q is the density, P sfs is the sub-filter production, D sfs is the sub-filter diffusion and T is the width of the temporal filter defined as [24,27] T ¼ r where the energy ratio r ¼ k sfs ktot is defined as the ratio of the modeled energy k sfs to the total turbulent energy k tot , e ¼ b Ã k sfs x is the dissipation rate, and C e1 = 1.44 and C e2 = 1.92 are constants. To take into account the cutoff frequency, Tran et al. [24] proposed a relation between r and the cut-off frequency, where U s ¼Ũ þ c ffiffiffiffiffiffiffi k tot p is the sweeping velocity [13].Ũ denotes the mean resolved velocity magnitude, and c is a constant.
is the highest cut-off frequency, with dt and D = max(dx, dy, dz) the time and grid steps, respectively, and b is a constant fixed at 0.667, calibrated in order to exhibit the correct amount of dissipation in the case of decaying homogeneous turbulence [24]. HTLES estimates the ratio r based on the analytical integration of the equilibrium Eulerian spectrum, between the cut-off frequency x c and infinity, with C j % 1.5 is the Kolmogorov constant. The ratio r is bounded by unity, which corresponds to the RANS mode, and the scale-resolving mode is triggered for values of r less than one. The total turbulent kinetic energy introduced above is defined as the sum of the resolved kinetic energy k res and the sub-filter kinetic energy k sfs , where k res is defined as, withŨ i the i-th component of the resolved velocity, and h. . .i is the statistical average operator. The parameters used in HTLES are averaged in time. However, k res is a statistical quantity that requires a relatively long time-averaging to converge. It was found (not shown here) that the simplified expression, does not affect the results significantly, although k res is now a fluctuating quantity. Moreover, this simplification is useful for future non-stationary ICE application, in which estimating statistical quantities by time-averaging is not possible. For such cases, using an estimate of U i , based, for instance, on time filtering, will be necessary.

Elliptic shielding
One of the main objectives of hybrid models is to operate in RANS mode in the near-wall region in order to reduce the cost of wall-bounded simulations. However, several hybrid methods do not properly enforce the RANS mode near the wall. For instance, DES is conceptually designed to cover attached boundary layers in RANS mode and to operate in LES mode in detached regions. However, applications to industrial flows show Modeled-Stress Depletion (MSD) and Grid-Induced Separation (GIS) when the model operates in LES mode within the boundary layer with insufficient grid refinement [28]. In the same way as DES, results of the PITM model in a channel flow at Re s = 395 show similar issues with resolved fluctuations occurring too close to the wall and increasing too rapidly when moving away from wall [19].
First applications of HTLES on wall-bounded flow exhibited similar problems as outlined in PITM simulations [19]. The ratio r predicted by the model can be less than unity near the wall (Fig. 1), which means that the model may operate in the scale-resolving mode near the wall which is undesirable.
Several solutions were proposed within this context. Shielding functions are added to DES in an enhanced version called Delayed Detached-Eddy Simulation (DDES) [28] in order to enforce the RANS mode in the boundary layer. Similarly, a shielding was proposed in the PITM framework by adding the wall-distance dependency to the energy ratio [19] using a parameter a based on the Elliptic Blending (EB) framework, such that the RANS mode is prescribed at walls.
In a similar way, the blending factor a is added to the expression for r used in HTLES where a is the solution of the following equation with a wall = 0 and the sub-filter length L sfs expresses as with C L = 0.161, C g = 80.

Near wall treatment
The HTLES model reduces the cost of scale-resolving simulations significantly near the wall by using RANS models in the near-wall region. To solve the boundary layer with RANS, the grid resolution must be such that the walladjacent cell is as small as y + = 2 in the wall-normal direction. This requirement must be verified in the whole domain, which can be challenging in high Reynolds number flows and complex geometries [29]. Further cost reduction for wall-bounded flows mandatory for industrial ICE simulations can be achieved through the use of wall functions, which can be used as an alternative to the explicit resolution of the whole boundary layer. HTLES is here combined with the automatic wall function [29]. This method offers a versatile formulation which gradually switches from the low-Re formulation to the log-wall function formulation depending on the grid resolution using a function of y + : -The specific dissipation rate x is computed from a blending formula of the analytical solution of x in the viscous sub-layer x visc ¼ 6m 0:075y 2 (y is the distance to nearest wall) and in the logarithmic region -The boundary condition for momentum equations is defined from the formulation of the friction velocity in the viscous sub-layer u s;visc ¼ U t y þ and in the logarithmic region u s;log where C = 5.1 and U t is the component of velocity tangential to the wall. The aforementioned equations require the a priori knowledge of y + . An approximation of the friction velocity u s is used to estimate y + where C l = 0.09. This formulation extends the use of standard wall functions, which are initially valid in the logarithmic region only.

CONVERGE CFD
HTLES has been implemented in CONVERGE CFD, which is widely used for ICE simulations. CONVERGE is a cell-centered code using the finite-volume method to solve the equations for compressible flows. The mesh is generated automatically using hexahedra far from the walls, while the cells near the walls are cut by the outer surface of the geometry using the cut-cell technique. The pressure-velocity coupling is achieved using a modified Pressure Implicit with Splitting of Operator (PISO) method [30]. The convection and diffusion terms are approximated by second-order central differencing, and they are advanced in time using the Crank-Nicolson second-order scheme. The Redlich-Kwong equation of state [31] is used for gas to couple the density, pressure and temperature. The Rhie-Chow [32] interpolation scheme is used to maintain collocated variable and to prevent pressure velocity decoupling and associated oscillations. A step-flux limiter is used to maintain physical solutions by restricting fluxes and mimicking a switch to a lower-order spatial discretization to avoid spurious oscillations in the solution.
The next section aims to assess the performance of the HTLES method in three statistically stationary flow cases.

Methodology
The predictions of HTLES are compared in three different cases (periodic hills, channel flow and steady flow rig) for which experimental or numerical reference data are available. The first part discusses the predictive capabilities of HTLES in periodic hills. This test case is known to be challenging for both RANS and LES. On one hand, RANS can predict the separation accurately, but shows difficulties to locate the reattachment correctly. Moreover, a large sensitivity of the results to the turbulence model used was reported in [33]. On the other hand, LES computational cost is large for solving the flow within the vicinity of the lower wall, which is essential to predict the flow correctly in the whole domain. Then, in order to evaluate the predictive capabilities of the method in the case of a simulation representing the near-wall region by wall functions, the channel flow and the steady flow rig cases are computed.
In order to assess the positioning of HTLES, RANS and LES computations are also performed. RANS simulations are performed using the unsteady two-equation k-x SST (Menter and Rumsey [34]) model combined with the automatic wall functions previously described. For LES, the r-model [35] combined with Werner and Wengle wall function is used [36].
Regarding RANS simulations, results of the periodic hills and the channel yielded steady solutions. Nevertheless, for the steady flow rig simulations, small-amplitude fluctuations have been observed and no steady solution could be reached and statistical average is thus performed. These fluctuations may be related to the dispersion error of the central differencing scheme, which can produce spatial oscillations acting as perturbations in the detached shear layer [37].
It is worth noting that the same grid is used for all simulations whatever the turbulence modeling approach (HTLES, RANS or LES) in order to focus on the single impact of the models on the results.

Flow over periodic hills
This configuration is part of the ERCOFTAC database and corresponds to a 2-D separating flow over periodic hills. Due to the hill slope, the flow separates downstream the hill, which results in separation and recirculation followed by reattachment on the flat wall upstream the subsequent hill. The Reynolds number based on the hill height h = 28 mm and the bulk velocity is Re h = 10 595.
The dimensions of the computational domain were examined in [38] and are thus chosen as L x = 9.0 h, L y = 3.035 h and L z = 4.5 h. Well-resolved LES with a grid of 12 million elements was performed in [38], which is taken as the reference data. The flow rate is imposed by a streamwise momentum source term, which is adjusted in order to provide the same mass flow rate as the reference. This configuration has been simulated with the k-x SST model and the HTLES with ES. The same grid, displayed in Figure 2, has been employed for both simulations. The mesh is refined along the bottom wall and in the recirculation region. The grid setup consists of 5.8 million cells. At the bottom wall, the grid resolution is refined such that the thickness of the cells adjacent to the wall is around 2 in wall units. Indeed, this typical test-case for hybrid model is commonly performed by solving the wall boundary layer [16,17] and this best-practice is here also retained.
A total time of 2.5 s was simulated, which approximately corresponds to 67 convective times t x ¼ Lx U bulk , where L x is the length of the channel and U bulk the bulk velocity. The flow was first established during 11 convective time and then the statistical average was performed during 56 convective times. The initial field is defined by a prior steady RANS solution performed using the same setup.
In Figure 3, iso-surfaces of the non-dimensional Q-criterion equal to 0.8 are used to identify the turbulent structures given by the HTLES model. The threshold retained for plotting the Q-isosurfaces highlights the complex turbulent structures occurring in the bottom part of the geometry while less structures can be seen on the top due to their lower turbulent intensity and the coarse grid resolution. Figure 4 represents the statistical average of the ratio r around the hill provided by the HTLES. Setting apart the near-wall regions where the ES enforces the RANS mode and thus r = 1, the domain exhibits very small values of r, which indicates that the scale-resolving method is active almost everywhere.
The streamlines computed from the results of the statistically-averaged velocities for the simulations performed and the reference data are compared in Figure 5. One may observe that the k-x SST model overestimates the recirculation length. This issue has already been encountered in previous studies using RANS models [33] since they fail to asses properly flows featured by massively separated structures. Unlike RANS, the streamline pattern shown by HTLES is very similar to the reference data. The center of the recirculation region is properly located, and the flow reattaches at the same location as the reference. Figure 6 shows the development of the streamwise velocity along the hill channel. As already mentioned, the k-x SST velocity profiles show limited accuracy in this configuration. At the top of the hill, the streamwise velocity is underestimated near the wall because of too late reattachment of the flow upstream the hill. It is observed at x h ¼ 2

Flow conditions
The channel flow consists of a flow at constant flow rate between two infinite parallel walls, as shown in Figure 8. Except for pressure, the flow is periodic in the x-direction. A streamwise source term equal to 176 Pa/m is added to the momentum equation in the x-direction in order to sustain the flow. The Reynolds number based on the friction velocity u s is 1000 and the reference data are extracted from DNS (Lee and Moser [39]).
The flow is initialized using a coarse DNS, performed on the same grid and using a first-order upwind scheme. A physical time of 2 s was simulated, which corresponds to 200 convective times (t x ¼ Lx U bulk ). Mean quantities are collected over a time period of 180 convective times and are also averaged in the spanwise direction over a distance equal to the channel height. The simulations are performed on a grid which contains N x Â N y Â N z = 500 Â 50 Â 100 hexahedra with an uniform distribution. The height of the cells adjacent to walls has been chosen to target y + = 30. This choice is motivated by the today practice for industrial ICE simulations. Further details of the simulations are summarized in Table 1. Figure 9 represents the dimensionless axial velocity as a function of the wall distance in wall units y + along the wall units y + of the channel given by HTLES with and without the Elliptic Shielding (ES). One can notice the impact of ES on the flow rate. Table 2 compares flow rates given by the simulations with DNS data. HTLES without ES shows a relative error of +11% compared to DNS data, while HTLES with ES yields accurate flow rate showing an error of À2% only. Figure 10 represents iso-surfaces of the non-dimensional Q-criterion equal to 0.2, based on U b and d, to identify the turbulent structures given by the HTLES model. Turbulent structures are solved explicitly in the core region whereas unsteady structures can hardly be seen near the walls. This observation shows that HTLES is able to operate in the scale-resolving mode in the near-wall region and sustain the turbulence by means of the direct control of the energy partition directly in the k-equation unlike some hybrid models that degenerate in RANS mode in this situation [13]. Figure 11 represents the statistical mean of the ratio r given by the two simulations. Using ES, the ratio r reaches unity at wall and decreases at a slope grade, which is less steep than in the simulation without ES. Far from wall, ES does not affect the results since the level of r is very close to the one without ES.

Elliptic shielding in HTLES
In-depth analysis is performed by examining the fluctuations given by HTLES with and without the ES as reported in Figure 12, which examines the impact of ES on the streamwise fluctuations. In the near-wall region, the resolved and the sub-filter fluctuations given by HTLES without ES are of the same order of magnitude. The scaleresolving mode is activated too close to the wall due to small values of the ratio r given by the HTLES model near the wall (Fig. 11). In HTLES with ES, the sub-filter fluctuations are increased near the wall and thus the scale-resolving region is shifted away from the wall. Although the targeted r plotted in Figure 11, evaluated by equation (10), goes to 1 at the wall, it is observed in Figure 12 that a significant part of the energy of the streamwise fluctuations is resolved. Indeed, r is only an ingredient of the model   that enters equation (4), and there is no mechanism in the simulation that imposes that the observed partition of energy among resolved and modeled scale satisfy this target ratio. A favorable side effect of the penetration of the resolved fluctuations in the RANS zone is that the level of resolved energy is high in the region where r rapidly falls from 1 to about 0.1, such that no modeled stress depletion is observed. It is observed that the near-wall peak of streamwise energy is underestimated in Figure 12 (right). This is due to the fact that, in this region, the simulation mostly rely on the RANS model, which is not able to reproduce this peak.

Comparison to RANS and LES
RANS and LES are also performed in the channel flow case, and velocity profiles are compared to HTLES and to the reference data in Figure 13. The k-x SST shows very accurate velocity profiles in this configuration, and the flow rate is almost the same as DNS with a relative error of À0.4%, as listed in Table 3. Unlike RANS, the velocity profile given by LES shows the same overestimation of the flow rate as HTLES without ES. Note that similar conclusions were already outlined in [8] when using wall functions in the channel flow. Indeed, LES combined with wall functions is not adapted for these configurations were the wall vicinity plays an essential role for the dynamics of the entire flow domain.

Configuration
The steady flow rig represents a simplified in-cylinder flow around a valve with a fixed lift. The gas is injected in the intake port, which consists of a cylinder and fixed valve placed axisymmetrically inside the jet nozzle. The flow exiting through the valve opening is representative of the intake flow into an IC engine geometry. An overview of the geometry with its dimensions similar to a common passenger car engine are displayed in Figure 14. The Reynolds number based on the diameter of the cylinder and the bulk velocity is 30 000. Laser Doppler Anemometry (LDA) measurements were investigated by Thobois et al. [40], which provide the mean and rms for axial as well as radial velocity profiles located 20 mm and 70 mm from the abrupt expansion (Fig. 14).
The results shown herein compare three simulations of the flow inside the steady flow rig using the same grid composed of 1 402 628 cells. The side view of the hexahedral grid is shown in Figure 15: the grid step in the core of the domain is 2 mm; a refinement of 1 mm is embedded at the inlet pipe and the head of the valve where the jet develops; at the outlet a coarse grid size of 8 mm is used as a sponge layer to avoid the reflection of acoustic waves inside the chamber. At the walls, the grid resolution is refined insofar as the wall distance is less than 100.
At the inlet of the computational domain, a uniform mean axial velocity profile is imposed without fluctuations with the bulk velocity U bulk = 65 m s À1 at 300 K, such that the experimental mass flow rate is recovered (0.06 kg s À1 ). A physical time of 2 s was simulated. The flow is established for a period of t = 0.3 s, equivalent to 8 flow through the intake port; then, the statistics are collected during a period of 44 flow through the intake port which is equivalent to 1.7 s. The static pressure at the outlet is set to 1 bar. The walls are assumed to be hydraulically smooth and adiabatic. The initial field used in the RANS simulation is defined as uniform, equal to the bulk velocity, while HTLES and LES initial fields are defined by the RANS solution.  simulations accurately reproduce the topology of the first recirculation region, with the same center location identified by a dashed line. Nevertheless, the second recirculation in RANS is longer than for LES and HTLES. Q-criterion is examined in Figure 17 in order to help visualizing the structures, only one half of the domain is shown. The grid resolution employed here is purposely insufficient for solving the turbulent structures upstream of the valve. Downstream, a part of the turbulent structures is solved where the mesh is sufficiently fine, far from the boundaries.

Quantitative comparison to RANS and LES
In order to assess the three approaches quantitatively, mean and rms velocity profiles along line 1 and line 2, respectively at 20 mm and 70 mm from the top of the cylinder (Fig. 14), are compared to the experimental measurements.
The mean axial velocity along line 1, given by the simulations, is compared in Figure 18    good agreement with the experiment. LES and HTLES mean velocity profiles are quasi-identical, and the velocity peaks are slightly better predicted compared to the k-x SST.
Discrepancies between the simulations can be observed for the mean radial velocity profiles as shown in Figure 19. HTLES and LES show similar radial profiles and provide a good assessment of radial velocity profiles within the recirculation region (0:5 < r R < 1), whereas to the k-x SST underestimates the radial velocity within the recirculation. The experimental data are distributed along a large deviation at this location and hence no clear conclusion can be drawn regarding velocity magnitudes.
The mean axial velocity profiles along line 2 are compared in Figure 20. Similar axial velocity profiles are obtained in HTLES and LES, which are in a good agreement with the experiments. The k-x SST overestimates the axial velocity in the central region where À0:5 < r R < 0:5, which confirms that the larger recirculation is overestimated by the RANS simulation.
In Figure 21, the radial velocity along line 2 (Fig. 14) is represented. One can notice the small velocity magnitudes at this location. Despite the axisymmetric configuration, both experimental and simulation data are not fully symmetrical. The asymmetry in experimental data can be related to an intrinsic uncertainty of the measurement device due to the small velocity magnitudes, while the     asymmetry in the simulations can be related to insufficiently converged statistics. HTLES and LES still provide results close to experimental data, whereas the k-x SST clearly under-estimates the radial velocity at this location.
Profiles are compared hereafter. Both the sub-filter and the resolved parts are taken into account. The axial rms velocity profiles are compared to LDA measurements in Figure 22. All simulations show acceptable results compared to the experimental data. HTLES provides results similar to RANS in the smaller recirculation region and nearly identical to LES in the larger one. These results are consistent with the previous observation since the RANS mode is enforced close to walls and the grid resolution is sufficient to make a transition toward the scale-resolving method. Finally, k-x SST results exhibit higher levels of fluctuations in the core region than HTLES and LES results.
It is observed in Figure 23 that RANS and LES both overestimate the peak of Ur rms at the location of the     annular jet r R AE 0:5 and mispredicts its location. HTLES accurately reproduces the shape of the profile everywhere except for the peak that is completely missed. Figures 24 and 25 show the profiles along line 2 for the axial and radial velocities, respectively, and the trends are the same as those already noticed for line 1. However, the experimental uncertainties, reflected by the asymmetry of the data, are such that it is difficult to draw any conclusion for rms velocity profiles along this line. Figure 26 compares the axial profile of the static pressure at the wall in order to assess the pressure drop in the different regions. All the simulations reproduce the stiff pressure decrease at x = 0 due to the sudden expansion, and the recirculation zone featured by a slight pressure increase at x/D = 1. Despite this good qualitative agreement, the global pressure drop defined as the difference of pressure between the position downstream the valve and outlet of the domain is different between the three simulations. While RANS and HTLES provide a prediction close to the experiment with a relative error of À3% (Tab. 4), LES is less accurate with a relative error of +11%. As already discussed in the channel flow section, this might be due to the wall modeling employed in the LES.

Analysis of HTLES modeled stress
The ratios r without and with ES along line 1 and line 2 are compared in Figures 27 and 28, respectively. The ES only affects the values of r in the near-wall region where it reaches unity, unlike the simulation without ES where r is much lower than one at the wall along both lines 1 and 2. This is consistent with the results obtained in the channel flow case (Fig. 11).
In order to perform a deeper analysis of the HTLES modeled stress, rms velocity for both axial and radial components are decomposed into the resolved and sub-filter parts in Figures 29 and 30. Regarding the sub-filter rms velocities, the evolution is consistent, for both components and for both lines, with the evolution of r displayed in Figures 27 and 28. The sub-filter fluctuations increase at walls and decrease in the core region for both components but are not negligible compared to the resolved ones.
The plots of the resolved rms velocities show that the axial contribution vanishes at the wall along both lines 1 and 2, as expected. However, the radial component surprisingly does not vanish at the wall along both lines. The same issue is observed in RANS and LES simulations in Figures 23 and 25, which points the difficulty of the code to impose tangential velocities to walls.

Conclusion
A first application of a hybrid RANS-LES method, which is intended to be used for simulations of ICE was presented. The proposed HTLES derived from the TPITM model offers a formalism based on the multi-scale approach. The HTLES model relies on the same set of equations as the k-x SST model and introduces a temporal filter width on the dissipation term of the k-equation, which is made sensitive to the grid and temporal resolutions. ES was adjoined to HTLES to protect the near-wall region from incursions of the scale-resolving method too close to walls. HTLES uses time averaged velocities for assessing the total turbulent kinetic energy which is required for the closing of the dissipation term of the sub-filter energy. A modification from the original HTLES is proposed by applying the averaging process to the resolved kinetic energy instead of the total one, and by removing the averaging process on the dissipation term in the balance equation of the sub-filter kinetic energy. This constitutes a first step towards the application of HTLES towards unsteady cases and ICE flows.
The HTLES method is applied to three cases and the results are compared with reference data as well as with the URANS k-x SST model and LES r model. The main conclusions are listed below: -First validation in the periodic hills has underlined a very good prediction capability of the developed HTLES for this separated flow, whereas the recirculation length obtained with the k-x SST model is severely overestimated. -Regarding the channel flow case, HTLES with ES solves the wall friction and the unsteadiness of the flow properly. The smooth variation with the wall-distance of the energy ratio induces a smooth transition from the RANS region and the scale-resolving region, and thus no modeled stress depletion is observed. RANS mode at wall yields correct prediction of the flow rate  on the one hand and the scale resolving simulation gives satisfactory assessment of the fluctuations on the rest of the domain.
-HTLES has been employed for the steady flow rig, which is a configuration close to ICE. Very good positioning of HTLES versus RANS and LES was noticed. On the one hand, for the same grid, the results of HTLES in terms of mean and rms velocity profiles are as accurate as LES unlike RANS, which exhibits some discrepancies with the reference data. On the other hand, the experimental pressure drop is properly predicted by HTLES in contrast to LES.
The application of HTLES implemented in the CONVERGE code shows promising results in stationary flows. The only remaining restriction for using HTLES in non-stationary flows is related to the statistical averaging of the velocity required for computing the resolved turbulent kinetic energy, and future work will be dedicated to this issue in order to make HTLES applicable to ICE configuration with moving walls.