Boundary Conditions and SGS Models for LES of Wall-Bounded Separated Flows: An Application to Engine-Like Geometries

The implementation and the combination of advanced boundary conditions and subgrid scale models for Large Eddy Simulations are presented. The goal is to perform reliable cold flow LES simulations in complex geometries, such as in the cylinders of internal combustion engines. The implementation of an inlet boundary condition for synthetic turbulence generation and of two subgrid scale models, the local Dynamic Smagorinsky and the Wall-Adapting Local Eddy-viscosity SGS model ( WALE) is described. The WALE model is based on the square of the velocity gradient tensor and it accounts for the effects of both the strain and the rotation rate of the smallest resolved turbulent fluctuations and it recovers the proper y3 near-wall scaling for the eddy viscosity without requiring dynamic pressure; hence, it is supposed to be a very reliable model for ICE simulation. Model validation has been performed separately on two steady state flow benches: a backward facing step geometry and a simple IC engine geometry with one axed central valve. A discussion on the completeness of the LES simulation (i.e. LES simulation quality) is given.


INTRODUCTION
Numerical simulation represents nowadays a very common approach to improve the understanding of the complex physical phenomena occurring in Internal Combustion (IC) Engines. A correct modeling of the turbulent motion is very important to achieve reliable predictions of the gas exchange phase (swirl, tumble), of fuel mixing and combustion, of the cyclic combustion variability.
In CFD engine simulation, the most widely used approach is the Reynolds Averaged Navier-Stokes (RANS) method. In RANS, the turbulence effects are taken into account without the direct solution of the turbulent structures but only by adopting specific turbulence models acting on the flow viscosity [1,2]. If the turbulent flow is modeled without directly solving the turbulence cascade, the computational cost results to be drastically reduced and, at the same time, an acceptable accuracy to reproduce the mean engine flow characteristics may be achieved. The main drawback of the RANS approach lies in the impossibility to correctly predict the intrinsically unsteady behavior of the turbulent flow into the engine. In this context, the numerical representation of the unsteady flow structures is a fundamental aspect for the overall quality of an IC-Engine CFD simulation, since they strongly affect complex phenomena like fuel-air mixing and cyclic combustion variability. From this point of view, Large Eddy Simulation (LES) applied to ICE is potentially a reliable tool to simulate turbulent motion in the engine cycle, without the loss of detail of the standard RANS approach. In LES, the Navier-Stokes equations are filtered in space, allowing a separation of the turbulent scales between Grid-Scales (GS) and Sub-Grid Scales (SGS). The GS part corresponds to the resolved large-scales of motion responsible for the momentum, heat, and mass transfer in separated or free shear flow regions. The SGS part corresponds to the more statistically isotropic and universal small not resolved residual motions.
Over the last decade the development of massively parallel computers and the improvement of CFD codes allowed to apply LES to industrial configurations of IC engines. However, before applying LES to real engine cases, this technique was firstly applied on simplified engine configurations. Thobois et al. [3,4] performed methodological studies about the LES applications to steady flow rings characterized by valves at fixed lift. The same steady state engine flow bench proposed in [3] was also carefully studied by the authors [5,6] adopting specific methodologies to evaluate the LES quality. Moureau et al. [7] presented the LES analysis of the generation and disruption of tumble vortexes inside a square piston engine. The results collected on six consecutive engine cycles showed encouraging results about the possibility to reproduce the mean velocity experimental profiles. The research work described in [8][9][10] demonstrated the potential of LES to capture the cycle-to-cycle engine variability. These studies represented only first steps towards LES of complex IC engine geometries; at the same time. the results showed were fundamental to carry out and highlight important insight that helped researchers to apply LES to more realistic engine configurations [11][12][13][14].
Despite the application of LES to IC-Engine simulation looks promising, still today there are some known issues that cannot be neglected. First, boundary conditions in LES must be able to properly handle both acoustic waves and turbulence properties. Then, the discrete solution of the Navier Stokes equations on unstructured meshes, typically used for real-world IC-engine applications, may lead to a numerical dissipation that is comparable to the eddy viscosity m sgs . Moreover, when simulations of complex geometries such as in internal combustion engines are performed, also the near-wall behavior of the eddy-viscosity represents a further difficulty [15]. This paper is divided in two main parts. In the first part, a condition for synthetic turbulence generation at inlet boundary is presented together with the implementation in an open-source CFD code [16] of two subgrid scale models, namely the classical Dynamic Smagorinsky SGS model [17] and the WALE SGS model [18]. In the second part, the implemented LES models are used to predict the flow conditions over two test cases: a backward facing step geometry [19] and a like IC-Engine configuration [3]. For both cases, LES results are compared to experimental LDA measurements. The quality of the performed LES simulations are evaluated by the Length Scale Resolution (LSR) parameter defined by the authors [20]. All the models described have been implemented in Lib-ICE Ò , a C+ + library based on the OpenFOAM Ò technology [21][22][23][24].

SYNTHETIC TURBULENCE INLET BC
In a LES, the velocity field at the inflow must be as close as possible to a ''real'' turbulent flow. In particular, it must reflect the salient characteristics of turbulence (like randomness, time-and space-correlation, and solenoidality), and it must correspond to the real flow with respect to some statistics (mean velocity, energy spectrum, lengthscales, etc.) [15,25]. One way to obtain ''turbulent-like'' velocity at the inlet is to artificially generate synthetic fluctuations with the desired characteristics, and superpose them onto a mean velocity profile [26]. As first step, the fluctuating component of the velocity is generated as a sum of Fourier modes: whereû n , w n are the amplitude and phase of Fourier mode n; j n j is the corresponding wavenumber vector, that is oriented in space according to angles h n and u n , and it ranges from j 1 to j N . Finally, r n j is the unit velocity vector of j n j (Fig. 1). In order to avoid that the generated fluctuations are destroyed by the Navier-Stokes solver, solenoidality of the velocity field must be enforced; this can be achieved by ensuring that the velocity vector r n j lies in a plane orthogonal to j n j [26]. The angle a n is the orientation of r n j on the ðn 1 ; n 2 Þ plane, as shown in Figure 1. Geometric angles / n , a n and h n and phase angle / n are random variables; they are generated for each mode n according to a given statistical distribution (Tab. 1). Spectral mode amplitudeû n is calculated according to a prescribed spectrum shape. In this work, a modified version of the Von Ka´rma´n energy spectrum has been used: The highest wave number of the artificial spectrum is taken as the spatial filter cutoff frequency, that is a function of the filter width: j max ¼ 2p=ð2DÞ. Conversely, the smallest wave number is defined from j 1 ¼ j e =p, where j e corresponds to the energy-carrying eddies lengthscale. Factor p should be larger than one to make the largest scales larger than those corresponding to j e . In the present work, p ¼ 2 [26]. The wavenumber space, j max À j 1 is divided into N modes (typically 150-600), equally spaced, of size Dj. If turbulence has to be generated for channel flow case, j g can be calculated as

Þ.
A fluctuating velocity field is then generated for each time step as described above. Fields generated at different timesteps are however independent of each other, and their time correlation is thus zero: this is unphysical. To enforce correlation in time, the new fluctuating velocity field, u 0 i is weighted with the previous one on the basis of an asymmetric time filter [27].
where (T is the integral timescale): In the simulations, the mean inlet profiles are either set from experiments or, for example, from the law of the wall. The synthetic fluctuations created yield isotropic turbulence in the inlet plane: u rms , v rms and w rms are constant in time, though they vary with the distance from solid walls. The algorithm for synthetic turbulence generation is applied at each temporal integration step, leading to a small increase of the total simulation time.

SGS STRESS TENSOR MODEL
In this work, the dynamic Smagorinsky [17] and the WALE (Wall-Adapting Local Eddy-viscosity) model Illustration of the symbols used in Equation (1). The wavenumber velocity r n i must lie in a plane orthogonal to the wavenumber vector j n j to ensure solenoidality of the velocity field [26].

TABLE 1
Statistical distribution of angles u n , w n , h n and a n Pðu n Þ ¼ 1=ð2pÞ 0 u n 2p Pðw n Þ ¼ 1=ð2pÞ 0 w n 2p Pðh n Þ ¼ 1= sinðhÞ 0 h n p Pða n Þ ¼ 1=ð2pÞ 0 a n 2p F. Piscaglia et al. / Boundary Conditions and SGS Models for LES of Wall-Bounded Separated Flows: An Application to Engine-Like Geometries [18] have been implemented in an open-source code [16] and were used in the simulations. Both the models are based on the eddy-viscosity assumption: where s ij is the subgrid scale tensor, defined as: and is the resolved rate-of-strain tensor. The overbar denotes the filtering operation (implicitly performed by the computational grid), and incompressibility is assumed. For most of the eddy-viscosity models, the eddy-viscosity can be written in the general formulation: where C m is the constant of the model, D is the subgrid characteristic length scale and OP is an operator of space and time, homogeneous to a frequency, defined from the resolved fields.

The Dynamic Smagorinsky Model
In the Dynamic Smagorinsky model [17], a test filterD bigger than the initial filter D (usuallyD ¼ 2D) is considered. The tensor operator OP in Equation (7) is based on the resolved rate-of-strain tensor like the classic Smagorinsky model: but the model constant C m is computed dynamically from the resolved scales: where L ij is the Leonard tensor (computable): and The superscript ''+'' in Equation (9) denotes a positive clipping of all the negative values to zero and the sign hi is a stabilization method that consists in a local volume averaging.

The WALE SGS Model
If compared to the dynamic formulation of the Smagorinsky model [17], the WALE model includes in the formulation of the operator OP the traceless symmetric part of the square of the velocity gradient tensor s d ij : where g 2 ij ¼ g ik Á g kj and d ij is the Kronecker symbol, X is the anti-symmetric part of g (or, the vorticity tensor): By construction, the trace of s d is zero and its second invariant remains finite and proportional to s d ij s d ij , so that for the incompressible case: where Since with the WALE model the invariant S ij S ij is zero in the case of pure shear, the model is able to reproduce the laminar to turbulent transition.
The formulation of the WALE model is based on the operator s d ij s d ij and this represents its main advantage with respect to the dynamic Smagorinsky model, since it is sensitive to both the strain and the rotation rate of the small turbulent structures. Because of this, it is well suited for LES in complex geometries with structured or unstructured methods because no explicit filtering is needed and only local information is required to build the eddy-viscosity. This is achieved by the definition of an operator OP with the following properties: -invariant to any coordinate translation or rotation; -easily assessed on any kind of computational grid; -function of both the strain and the rotation rates; -it goes naturally to zero at the wall so that neither damping function nor dynamic procedure are needed to reproduce the effect of the no-slip condition.
The resulting formulation for the eddy viscosity then is: In Equation (15), the ratio OP 1 =ðS ij S ij Þ 5=2 ensures to keep the y 3 behavior of the whole operator near the wall and to have a formulation of operator consistent with Equation (7). At the same time, OP 1 =ðS ij S ij Þ 5=2 is not well conditioned numerically since the denominator can locally tend to zero while OP 1 remains finite. The way used in [18] to avoid this situation is to scale OP 1 by ðS ij S ij Þ 5=2 þ ðs d ij s d ij Þ 5=4 ; the second term in the denominator is negligible near the wall but it avoids numerical instabilities because OP 2 does not go to zero for pure shear or irrotational strain. Although the dynamic procedure could also be applied to the WALE model, in this work C w has been considered as a true constant, assessed from the case of isotropic homogeneous turbulence [18].

ANALISYS AND DISCUSSION OF THE RESULTS
A validation study for the implementation of the proposed models and of the simulation methodology has been applied on the Backward-Facing Step (BFS) testcase [19], for which an extensive set of experimental measurements was available. Moreover, it has features that make it a good validation test case for engines, since it includes phenomena that are similar to the ones occurring during the intake stroke of the engine, when air flows to the combustion chamber through the valve ports. Finally, to explore the suitability of the proposed approach for reciprocating internal combustion engines, a simplified IC geometry with one axis-centered valve [3] has been simulated.
The main goal of a LES simulation is to directly solve all the non-isotropic turbulent scales over the inertial sub-range. Therefore, the quality of a LES simulation can be directly referred to the evaluation of the flow turbulent energy budged directly solved during the simulation itself.
A method to evaluate the turbulence resolution for a LES simulation is to use the Length Scale Resolution (LSR) parameter defined by authors as: where D is the local filter size and l di is the lower limit of the inertial sub-range, that is usually estimated as [28]: where g is the Kolmogorov scale. The LSR parameter is proportional to the deviation between the actual resolved energy level and the corresponding lower limit of the inertial sub-range. Where the LSR value is equal to 1 all the turbulent scales up to the viscosity range are resolved. By the LSR definition, the evaluation of the actual resolved energy level is directly linked to the local filter size. It helps to clearly match the adopted mesh size to the local energy resolution all over the computational domain. In [20,29], the authors found that a LSR value of 3-5 is the upper limit to guarantee a reasonable LES resolution at an affordable computational cost.

Backward Facing Step Geometry
The single-side expansion of Figure 2 was considered in this study. The results obtained by LES for mean velocity profiles and for velocity fluctuations have been quantitatively compared to the experimental measurements by hot-wire anemometry performed by Eaton et al. [19] in terms of both mean and RMS velocities. In the experiments, the turbulent boundary layer was generated in the inlet channel of the backward facing step by positioning two trips, that were placed at the same distance from the step: the first one was located on the bottom wall and the second one on the top wall. Trips heights were 2.5 mm and 1.25 mm respectively; the boundary layer thickness on the walls was therefore different. The inlet stream, at Re ¼ 40 000 (based on the bulk velocity U 0 at the inflow) was turbulent and fully developed at the step. In the simulations, only a small portion of the domain of the BFS was modeled, as shown in Figure 3; the inlet channel length L 1 was chosen to have fully developed flow near the step, while the outlet channel length L 2 was chosen long enough to not have the numerical solution affected by the vicinity of the outlet boundary. Experimental setup of the backward facing step geometry by Eaton et al. [19].

Computational Procedure for BFS
In the present case a two-block orthogonal, fully-structured mesh system was used, with one block covering the upstream channel and the other one covering the downstream channel of the BFS. Simulations were carried out by an open-source CFD code [16], based on the finite volume approach. Second-order central differencing schemes in space for advection and diffusion were blended with linear-upwind schemes to stabilise solutions while maintaining second-order behavior [30]. The schemes applied result to be fully conservative and since the coefficients are always positive they are unconditionally bounded. Also, they satisfy the transportiveness requirement for large values of local (cell) Peclet numbers. Unfortunately, no higher order schemes can be applied in the FV context when unstructured grids are used [31]. The numerical setup described in this paragraph will be kept for all the simulations carried out in this work. Time-marching was implicit, with the time derivative being discretised by a second-order backward approximation. The flux terms were advanced explicitly; the provisional velocity field was then corrected via the pressure gradient by a projection onto a divergence-free velocity field; pressure was computed as a solution to the pressure-Poisson problem. In order to improve the global convergence of the simulation, the operator splitting technique for pressure-velocity coupling was implemented as a transient-SIMPLE algorithm [2], where convergence of pressure-velocity solution is enforced by iterating the coupling procedure for each time step; the resulting fluid-dynamic CFL number, defined as CFL ¼ jujDt=Dx, has been set to 10, with significant advantages in terms of required wall time.
The code and the developed libraries were fully parallelised. The fluid-dynamic information recorded by a RANS simulation on a precursor domain has been used to reconstruct the turbulent fluctuations at the boundary inlet by the synthetic turbulence inlet bc described in Section 1.
An unsteady convective boundary condition has been used for the outlet: where hU i is the normal velocity derived by the fluxes on the cell face located at the boundary end. Walls were set as adiabatic. The subgrid stress tensor has been modeled by the dynamic Smagorinsky model [17] and by the WALE model [18]. While the first model does not make use of any fixed-value parameter, the WALE model makes use of a constant C w that is flow-dependent. In the simulations presented in this work, an average value of C w ¼ 0:58 has been taken from the literature, since it has been demonstrated that it is not significantly varying for very different types of flows [18]. Since both models are self-adapting in the vicinity of solid walls, no wall functions were applied and flow equations were solved up to y þ ! 0. The BFS domain was discretised by two different meshes, having 1.2 M (case ''coarse'') and 2.2 M cells (case ''fine'') respectively. This is in agreement with the studies of other authors [32].
Grids are the outcome of precursor testing. Local refinement was used near the walls to ensure the maximum value of y þ within the limits reported in Table 2. In both cases, the width of the computational domain in the spanwise direction was set to 2H and cell-aspect ratio along the spanwise direction was maintained constant for the two cases. The spanwise direction was treated as statistically homogeneous, with periodic conditions prescribed at the boundaries. Figure 4 shows distributions of the auto-correlation coefficient R uu along several lines at different locations along the x-axis within the separated shear layer. As it can be seen, the spanwise extent is sufficient to ensure uncorrelation along that direction. It is important to note that mesh resolution near the walls was not high enough to capture the full wall turbulence cycle, which requires also Dx þ < 10 Schematic of the backward facing step geometry simulated [19]; and Dz þ < 5 [15]. This was a deliberate choice of the authors, since the final goal of the work was to find the best setup to perform reliable LES simulations of internal combustion engines, where the correct mesh resolution near the walls along non-normal directions implies a number of computational cells that is not affordable with the common computational resource limits.
In the authors' experience, when generating a grid for axis-symmetric geometries (like intake/exhaust ducts or combustion chambers), constraining the first off-wall point to be at y þ < 1 has a lower impact on the overall mesh size, with respect to axial (Dx) or circumferential (Dz) refinement. In fact, the latter affect also the cell size near the cylinder axis, lowering the maximum timestep allowed by the CFL criterion, while the former leads to a greater increase in total cell number since the axial dimension is usually dominant over the other two.
The main characteristics of the grids used for the simulation of the BFS are reported in Table 2.

Results for Backward-Facing Step
The backward-facing step case is an example of a wallbounded flow, it involves reattachment of separated turbulent shear layers and it is a very good test-case to validate turbulence models near the walls. The main quantities to be analized are the first and the second moment of the velocity field, namely, the mean velocity hU i and the RMS streamwise fluctuations hu 0 u 0 i. For each of them, experimental measurements by hot wire anemometry were available for comparisons. In Figure 5, the predicted mean velocity profiles are shown at different locations over the domain. Velocity has been averaged both in time and along the spanwise direction. Time averaging has been performed over five flow-through times after the steady state condition was reached. The agreement between simulations and experiments looks satisfactory for all the cases considered, expecially for the separated region after the step (x=H ! 4); moreover, velocity profiles downstream of the step seem quite independent from both the mesh and the SGS model. As expected, there are no significant differences between the statistical properties of the solutions obtained by the dynamic Smagorinsky and the WALE model, at least in the bulk region of flow. As the post-separation and reattachment processes are dominated by large-scale vortices, the influence of subgrid modelling is expected to be generally weak, provided that the grid is reasonably fine. When RMS streamwise velocity fluctuations hu 0 u 0 i are considered (Fig. 6), the influence of the subgrid model and of the mesh is then apparent. Simulations using the dynamic Smagorinsky model reproduce the experimental curves quite well, even when a coarse mesh is used (Fig. 6, first  row), while the WALE model seems to be more sensitive to the level of refinement of the mesh used. When the fine mesh is adopted, both models exhibit predictions that are in good agreement with experimental data (Fig. 6, second  row). The coordinate of the reattachment point of the recirculation bubble cannot be inferred from the plots of Figure 5, since the experimental dataset used for comparison did not include any point located at y=H < À0:85.
The reattachment point, however, can be precisely determined by looking at the coordinate where the wall friction coefficient:  changes its sign. The wall friction coefficient C f calculated along the step wall is plotted in Figure 7.
Although the trend of C f is qualitatively captured by both models, significant differences exist. The dynamic Smagorinsky model (Fig. 7a) seems to better predict the reattachment length (Fig. 7b), even though the WALE model better predicts the absolute values in the vicinity of x=H ¼ 5. All configurations show a significant difference in the predicted values of C f for x=H > 10, and this discrepancy increases towards the outlet end. A possible explanation can be found in the shorter length of the simulated channel with respect to the experimental configuration; as a consequence, a stronger pressure gradient is generated at the boundary and this influences the upstream flow. The region upstream of the step (x=H < 0) can be regarded as a plane channel and studied accordingly. Here, the main turbulent dynamics are represented by near-wall streaks and ejections, that have to be completely resolved to fully account for turbulence production [15]. Grid resolution at solid walls is of foremost importance and a rough guideline must be given in terms of wall units. In particular, the following values are thought as the maximum allowable cell size for a complete representation of near-wall turbulence [33]: y þ < 1, Dx þ % 100, Dz þ % 30. Whereas the requirement on the wall-normal resolution (y þ ) can be easily achieved without increasing significantly the computational load, the other two requirements are seldom fulfilled. The near-wall cell size used in the present work is reported in Table 2.
In Figure 8, the near-wall profiles for three stations are plotted in scaled coordinates and compared against the theoretical Spalding's law-of-the-wall [34]. The coarser mesh has only one grid point in the viscous sublayer (y þ < 1) along the wall-normal direction. However, the laminar sublayer is correctly resolved for all cases, though it extends to a wall-distance greater than expected. On the other hand, both models overestimate the kinetic energy in the outer region (Fig. 8a-c), and the log-law behavior develops for y þ > 100. This can be probably ascribed to the insufficient streawise grid resolution that leads to an incomplete representation of the near-wall turbulence dynamics. In addition, it might be possible that the length of the inlet channel is too short to allow for a full development of turbulence, even though a synthetic turbulent inlet is applied at the inflow. In this region, the refinement of the mesh does not lead to a significant improvement in the level of accuracy: while the WALE model seems to get closer to the theoretical curve, predictions obtained by the dynamic Smagorinsky model become worse (Fig. 8d-f). The different near-wall behavior of the WALE model with respect to the dynamic Smagorinsky can be investigated by examining Figure 9, where the subgrid viscosity is plotted against wall-normal scaled distance. Both models exhibit a decreasing trend while approaching to the solid wall, but the slope of the WALE curve is steeper, since it asymptotically tends to m sgs / ðy 3 Þ. This is consistent with the theory [18]. Conversely, the dynamic Smagorinsky curve tends to flatten as the wall distance decreases, and this effect is greater when a coarse mesh is used. Figures 10-13 show some results about the evaluation of the LES simulation quality; the vorticity field and LSR parameter calculated from simulations running with the dynamic Smagorinsky and the WALE SGS models are showed. The vorticity field does not seem to be much sensitive to mesh resolution when the dynamic Smagorinsky model is applied: the contour plots of Figure 10a and b look very similar, suggesting that turbulence is advected and dissipated in the same way by the coarse and the fine grid; this is also confirmed by the mean velocity fluctuations Comparison of computed and experimental friction coefficient C f for different SGS models. a) Dynamic Smagorinsky; b) WALE. À À Experimental, -coarse mesh, ÀÀ fine mesh.  Comparison of near-wall velocity profiles by SGS model upstream of the step (x=H < 0). First row (a-c): coarse mesh; second row (d-f) fine mesh. À À Dynamic Smagorinsky, ÀDÀ WALE; À Á À law of the wall. Comparison of near-wall SGS viscosity for different meshes and models. First row (a-c): coarse mesh; second row (d-f) fine mesh. À À Dynamic Smagorinsky, ÀDÀ WALE; ---y 3 slope. profiles (Fig. 6), where the experimental data are matched quite well by both grids. On the other hand, simulations performed with the WALE model show higher differences between the coarse grid and the fine grid solution, as shown in Figure 11a and b. In particular, a larger cell size causes too strong a dissipation of turbulent structures, as it can be seen near the outlet of Figure 11a; this reflects in a wrong estimate of fluctuating velocity some lengths after the step, as shown in Figure 6.
Independently by the SGS model applied, the maximum value of LSR in the inlet channel decreases in the vicinity of the inlet boundary, from % 5 (coarse mesh) to % 3 (fine mesh). This means that the filter size of the fine mesh in the inlet channel is closer to the lower limit of the inertial subrange: hence, the flow energy budget directly solved during the simulation is higher. Figures  12 and 13 clearly show also that a better mesh resolution (and an improvement of the quality parameter LSR) before the step has a direct influence on the flow resolution after the step. In detail, by the streamwise velocity fluctuation profiles showed in Figure 6 it is possible to recognize the improvement in the agreement between predictions and experiments over the channel height of the BFS, when the fine mesh is used together with the WALE SGS model.

Cold Flow IC Simulation
The second part of the work concerns with LES simulation of a simplified IC engine configuration, that is shown in Figure 14. The geometry consists of a circular pipe with a sudden expansion from d ¼ 34 mm to D ¼ 120 mm (diameters ratio is about 3.53). A single axis-centered poppet valve is positioned across the expansion with a fixed lift of 10 mm, to originate a circular jet that expands inside the larger cylinder. The side of the larger pipe opposite to the valve is an open end, and the flow is driven by the difference of total pressure between the two pipe ends. Air enters the smaller pipe with a mean Reynolds number of about 30 000, which corresponds to a bulk velocity of 65 m/s. The flow inside the smaller annular duct can be assimilated to a circular pipe flow; however, this region is of less interest since most turbulence production takes place at the shear layer between the circular jet and the air into the cylinder that is, initially, at rest. Two large toroidal recirculation zones originate inside the cylinder: one is located between the jet and the top wall, while the other one is close to the cylinder axis. Although simplified, this configuration can be consid-ered as representative of the main flow types that occur in a real engine during the induction stroke, where a complex jet coming from the intake valve enters the cylinder at high speed and it originates large-scale motion of the charge. LDA measurements of mean flow velocity and of RMS fluctuations (along the radial and tangential directions) are available on two planes located at a distance of 20 mm and 70 mm from the cylinder top, respectively. Similar studies on this configuration have already been done in [33].

Computational Procedure for the Cold Cylinder Flow
Two computational grids were used for the same geometry: the former had 1.4 million of cells, the latter had about 13 million of cells. In both cases the mesh was a block-structured fully hexhaedral, with near-wall refinement in the inlet duct and on the valve sides, as shown in Figure 15.
The fluid is air at ambient temperature and Mach number at the inlet is 0.19, which is below the usual limit for compressible flows (Ma ¼ 0:3). Since the curtain area   of the valve is greater than the inlet pipe cross section, it can be assumed that the incompressible hypothesis does preserve its validity.
Simulations were performed by the WALE SGS model with C w ¼ 0:58 as in the BFS case. Timeresolved velocity field has been averaged in time and   along the circumferential direction to calculate first and second statistical moments. Due to the high magnitude of the advective fluxes in the momentum equation, a limited scheme had to be applied to avoid spurious oscillation in the velocity field. However, this introduces artificial viscosity, which is not a desirable effect in LES since it is a source of errors. Moreover, when using a limited scheme, the mesh influences the numerical dissipation, as shown in [35]. Both effects have to be taken into account when analyzing the results.
In this work, a particular flavour of the Van Leer limiter, particularly suited for vector fields [36] has been applied. Time derivatives were discretized with a second-order backward finite difference. The same algorithm as the BFS case has been used to solve the discretized equations.

Results for the Cold Cylinder Flow
Comparisons of the current LES with LDA measurements are shown in Figure 14. Quantities of graphs a-c are referred to the upper plane (20 mm below the cylinder top), while plots d-f are referred to the lower plane (70 mm below the cylinder top). Both on the upper (Fig. 14a) and on the lower plane (Fig. 14a), mean  Comparison of the computational grids used for cold cylinder flow simulations. On the left is the "coarse" grid (1.4 M cells), on the right is the "fine" grid (13 M cells). The cylinder length has not been represented in its entirety to highlight the valve region.
velocities are captured rather well by both the coarse and the fine mesh, with a slight superiority of the latter with respect to the peak value around r=R ¼ 0:5.
On the other hand, RMS of the fluctuating component of the velocity exhibits a significant deviation from the measured values. The LES with the finer mesh overestimates the axial fluctuations on both planes, whereas the coarser mesh seems to adhere better to the measured points. With respect to the circumferential fluctuations ( Fig. 14c and f), there are no significant differences between the meshes, but a noticeable overestimation of the values on the upper plane still exists (Fig. 14c). Conversely, on the bottom plane the circumferential fluctuations (Fig. 14f) are captured quite well, although their value is almost uniform across that surface.
In Figure 14c-f, RMS fluctuations look better predicted when the coarse mesh is used. This can be explained by the reciprocal cancellation of modeling and discretization errors; in fact, a too low value of the model constant C w , which has been tuned for the channel-flow case, may cause the overestimation of fluctuating velocity observed for the fine mesh, as a consequence of the underestimation of the subgrid viscosity m sgs . On the other hand, when solving the filtered equations on a non-Cartesian mesh, a first-order (dissipative) numerical error, that acts like an artificial viscosity added to the eddy viscosity, arises [37]: m eff ¼ m lam þ m sgs þ m num . This extra dissipation may compensate the aforementioned modeling error, so that the correct value of m eff might eventually be restored. An extra argument supporting this consideration might come from the resolved vorticity field, Figure 16. The increasing dissipation of turbulent structures arising when a large cell size is used, that has been already noted with respect to the BFS case, is clearly shown in Figure 16a.
The LSR parameter, shown in Figure 17 for both meshes, can be used to estimate the completeness of the scale resolution. It is clearly visible that the fine mesh has lower values of LSR almost everywhere: this means that the smaller resolved scales are closer to the lower limit of the inertial subrange. However the global quality of the simulation, evaluated in terms of match between simulations and experiments, cannot be estimated by the LSR parameter only: in the case shown, the coarse mesh results lead to better predictions, despite the LSR for that case is higher. There are other factors, like cell skewness, mesh non-orthogonality and cell shape that influence the final result. The conclusion that can be drawn at this point is that LSR can only be used to compare two meshes generated with the same strategy (like in the case of the BFS) or to detect those region that need to be refined in order to improve the results (like in the case presented in this section).

CONCLUSIONS
In this paper, evaluation of SGS models to apply to the simulation of IC engines has been carried out. The WALE SGS model [18] and the classic dynamic Smagorinsky model [17] have been implemented in OpenFOAM Ò and applied to the simulation of standard test cases, that have some features in common with IC engine flow. Moreover, a synthetic-turbulence inlet based on the procedure of Davidson [26] has been implemented by the author and used in all simulations as inflow boundary condition.
The first case considered is the backward-facing step geometry [19]. Two computational meshes were generated for this case, with different cell sizes, and simulations were run with each of the SGS models. Results on this geometry show that both subgrid models are able to capture the mean velocity with good accuracy. There is a sligth superiority of the WALE model for channel flows, where it better reproduces the near-wall velocity profile. On the other hand, the dynamic Smagorinsky model proved to be less sensitive to mesh resolution, since it is able to correctly capture the RMS fluctuations in coarse meshes, where the WALE model fails. With a small cell size, both model are able to correctly predict the RMS value of velocity fluctuations, with little differences between the models. The LSR parameter has been computed in order to estimate the grid quality in terms of cell size, and it has proven useful to point out those mesh regions that need local refinement.
The second case is the simplified engine geometry by [3]. For this case, only the WALE model has been adopted, again with two different computational grids. There is a good correspondance between experimental and computed mean velocities, on both the coarse and the fine mesh. RMS fluctuations are solved with less accuracy with respect to the previous case, probably because of uncertainties in the value of the model constant; unexpectedly, turbulence is better predicted by the coarser mesh, probably due to a fortuitous reciprocal cancellation of modeling and discretization errors. Tests on engine-like geometry evidenced that LSR parameter does not provide reliable comparisons between computational grids having different block structures.
As stated in the introduction, this work represents a preparatory stage towards the application of LES to IC engines. Some of the basic assumptions of this work will no longer hold when simulating real engine cases, that require both moving piston and models for compressible flows. However, a significant part of the methodology shown here (choice of the mesh resolution, applicability of the SGS models based on the WALE operator) will be still valid when increasing the complexity of the simulation. Work on compressible flows will be focused on the implementation of the dynamic WALE model [38] in its compressible formulation and of a non-reflecting inflow with synthetic turbulence generation.