Numerical Methods and Turbulence Modeling for LES of Piston Engines: Impact on Flow Motion and Combustion

Re´sume´ — Me´thodes nume´riques et mode`les de turbulence pour la LES de moteurs a` pistons : impact sur l’ae´rodynamique et la combustion — Cet article pre´sente une e´valuation de l’impact du set-up nume´rique sur l’ae´rodynamique interne et la combustion pre´dite par Simulation aux Grandes E´chelles (LES, Large Eddy Simulations) dans les moteurs a` combustion interne. Du fait de la complexite´ et du couˆt de calcul important associe´s a` ce type de simulation, le set-up le plus classique consiste a` utiliser des sche´mas d’ordre faible (typiquement premier ou second ordre en temps et en espace) et des mode`les de turbulence de sous-maille simples (comme le mode`le de Smagorinsky (Smagorinsky J. (1963) Mon. Weather Rev. 91 , 99-164)). L’objectif de ce travail Abstract — Numerical Methods and Turbulence Modeling for LES of Piston Engines: Impact on Flow Motion and Combustion — In this article, Large Eddy Simulations (LES) of Spark Ignition (SI) engines are performed to evaluate the impact of the numerical set-up on the predicted ﬂow motion and combustion process. Due to the high complexity and computational cost of such simulations, the classical set-up commonly includes “low” order numerical schemes (typically ﬁrst or second-order accurate in time and space) as well as simple turbulence models (such as the well known constant coefﬁcient Smagorinsky model (Smagorinsky J. (1963) Mon. Weather Rev. 91 , 99-164). The scope of this paper is to evaluate the feasibility and the potential beneﬁts of using high precision methods for engine simulations, relying on higher order numerical methods and state-of-the-art Sub-Grid-Scale (SGS) models. For this purpose, two high order convection schemes from the Two-step Taylor Galerkin (TTG) family (Colin and Rudgyard (2000) J. Comput. Phys. single cycle computations don’t permit to distinguish clear improvements on macroscopic parameters such as resolved kinetic energy, heat release or mean in-cylinder pressure. Concerning numerical schemes, TTG schemes also allow a slighlty better resolution of small scale vortices but global quantities such as resolved kinetic energy and SGS viscosity are comparable. Nevertheless, clear differences appear between the different schemes in the combustion stroke. This is attributed to a better resolution of the ﬂame-turbulence interaction process during the free ﬂame propagation period, leading to an increase of the resolved part of heat release. It is also shown in this paper that an adjustment of the efﬁciency constant in the Thickened Flame (TF) model is compulsory to account for the over dissipation of the smallest resolved structures if LW is used. In the light of these conclusions an hybrid set-up, called ES O 2 (Engine Stroke Optimal Order), which consists in using TTGC during combustion and LW elsewhere is proposed and applied to the two engines conﬁgurations. Results are in good agreement with the ones obtained in the case of a full TTGC simulation, while the CPU (Central Processing Unit) cost increase is only about 10% compared to LW. The accuracy of LW seems therefore to be sufﬁcient for pure aerodynamic phases, while the use of TTGC only during combustion permits an improvement in the LES quality. The hybrid ES O 2 method thus appears as an attractive approach to improve further calculations accuracy without being greatly penalized by additional CPU costs in multi-cycle simulations.

single cycle computations don't permit to distinguish clear improvements on macroscopic parameters such as resolved kinetic energy, heat release or mean in-cylinder pressure. Concerning numerical schemes, TTG schemes also allow a slighlty better resolution of small scale vortices but global quantities such as resolved kinetic energy and SGS viscosity are comparable. Nevertheless, clear differences appear between the different schemes in the combustion stroke. This is attributed to a better resolution of the flame-turbulence interaction process during the free flame propagation period, leading to an increase of the resolved part of heat release. It is also shown in this paper that an adjustment of the efficiency constant in the Thickened Flame (TF) model is compulsory to account for the over dissipation of the smallest resolved structures if LW is used. In the light of these conclusions an hybrid setup, called ES O 2 (Engine Stroke Optimal Order), which consists in using TTGC during combustion and LW elsewhere is proposed and applied to the two engines configurations. Results are in good agreement with the ones obtained in the case of a full TTGC simulation, while the CPU (Central Processing Unit) cost increase is only about 10% compared to LW. The accuracy of LW seems therefore to be sufficient for pure aerodynamic phases, while the use of TTGC only during combustion permits an improvement in the LES quality. The hybrid ES O 2 method thus appears as an attractive approach to improve further calculations accuracy without being greatly penalized by additional CPU costs in multi-cycle simulations.

INTRODUCTION
In the past few years, Large Eddy Simulation (LES) has been a subject of growing interest from the automotive community because of its unique potential to reproduce unsteady and sporadic phenomena like Cycle-to-Cycle Variations (CCV) or abnormal combustions [1][2][3][4][5][6]. However, Internal Combustion Engines (ICE) still remain a recent and complex field of application for LES: the flow is highly unsteady and governed by a strong interaction between numerous physical phenomena (turbulence, mixing, combustion, multi-phase flows, acoustics, etc.), the geometry is mobile and initial and boundary conditions are generally badly characterized or defined (wall temperature, inlet/outlet pressure or velocity signal). When dealing with sporadic or erratic phenomena such as CCV or abnormal combustions, an additional difficulty arises: since reliable trends and statistics may only be obtained by computing numerous cycles (typically 50), the solver must necessarily be robust and fast. A direct consequence of these issues (complexity, novelty, cost and robustness) is that classical ICE simulations do not generally use "high order" set-up compared to more academic LES configurations, such as turbulent pipe flows for instance [7], for which highly accurate set-up are the standard. In particular, ICE simulations usually use low order and/or dissipative numerical schemes as well as simple Sub-Grid-Scale (SGS) turbulence models. Concerning numerical schemes, two classes of convective schemes are generally found in the literature: upwind-biased schemes are used by Jhavar and Rutland [8], Dugue et al. [9] or Goryntsev et al. [10] for instance while the Lax-Wendroff total discretization approach [11] (second order accurate in space and time) is used by Richard et al. [12], Vermorel et al. [6] or Granet et al. [13]. Both of them are known to be very dissipative and a priori not well suited for LES [10]. Concerning turbulence modeling, the most popular closure is the constant coefficient Smagorinsky model [14]. In spite of its well-known drawbacks [15], this model is used in many ICE LES, such as Goryntsev et al. [10], Celik et al. [16], Vermorel et al. [6] or Granet et al. [13]. The combined use of a low order numerical scheme and a simple turbulence model does not mean that the results of these computations are systematically wrong or dubious. Many other parameters have also to be taken into account (the resolution and the quality of the grid especially) to state if a methodology is adequate or not. Besides, many promising results have been obtained with such a set-up. For example, previous works conducted at IFP Energies nouvelles and CERF-ACS using the Lax-Wendroff scheme and the classical constant coefficient Smagorinsky model have demonstrated the great potential of LES for solving complex problems in piston engines, notably its ability to help understanding CCV sources [4,13,17] and to build phenomenological models for engine control development [18].
However, it is well known that both numerics and turbulence models can have a huge impact on the simulation results. In the perspective of using LES for predicting engine operation far upstream in its design process, highly precise methods should then be used to provide more reliable results. This objective is laudable but not necessarily realistic: high order nondissipative schemes generally exhibit high-frequency artificial oscillations, which may cause instability in the computations and jeopardize the robustness of the solver. High order set-up may also lead to prohibitive computational costs for engine simulations, which require tens of cycles for generating reliable converged statistics. The objective of this work is therefore to evaluate the feasibility and the potential benefits (or drawbacks) of using high precision numerical schemes and state-of-the-art SGS turbulence models in terms of precision, robustness and cost. For this purpose, two convective schemes from the Two-step Taylor-Galerkin (TTG) family are tested and evaluated in combination with several SGS models (Smagorinsky, dynamic Smagorisnky, sigma). This evaluation is performed considering two different engine configurations of IFPEN: a naturally aspirated engine and a highly supercharged engine. Thanks to these comparisons, an optimal setup is highlighted and a new approach called ES O 2 (Engine Stroke Optimal Order) method, based on the choice of the best compromise in terms of CPU cost and precision, is proposed.

TEST CASE DESCRIPTION AND APPROACH
The aim of the present study is to propose an appropriate numerical set-up which would be able to handle piston engine simulations with the highest fidelity. For this purpose, three different numerical schemes and three different SGS models are evaluated and compared on two different engine configurations. The three numerical schemes are the classical finite volume Lax-Wendroff (LW) scheme [19] (2 nd order accurate in space and time) and two Finite Element (FE) schemes from the 2-step Taylor-Galerkin family: TTGC (3 rd order accurate in space and time) and TTG4A (4 th order accurate in time and 3 rd order accurate in space) [20]. For academic configurations (convection of a vortex, HIT), these two FE schemes give much better results than the LW scheme thanks to their better dispersive and dissipative properties [20]. They have also been widely and successfully used in complex configurations such as gas turbine combustion chambers in recent LES studies [8,21]. Concerning the SGS models, in addition to the classical constant coefficient Smagorinsky model [14], two more recent closures are also evaluated: the dynamic Smagorinsky [15] and the sigma [22] models (see Appendix). Many studies have already highlighted the conceptual and actual advantages of the dynamic procedure compared to the constant coefficient method [23,24]. The sigma model is also supposed to give better results than the classical Smagorinsky model, especially in the vicinity of solid boundaries or in case of pure shear and solid rotation for instance [22]. The two configurations are the F7P and the Ecosural single cylinder engines from IFPEN. The F7P configuration has been widely studied experimentally and numerically in previous works.
In [4,13,25], multicycle LES of different motored and reactive operating points are presented. The numerical set-up includes a LW convection scheme and a constant coefficient Smagorinsky model. Despite its supposed "low" accuracy, this set-up exhibited promising results: for various operating points, experimental results in terms of aerodynamics and combustion cycle-to-cycle variability were correctly reproduced not only qualitatively but also quantitatively. The Ecosural engine is a modern single cylinder engine, specifically designed to support research activities. Experimental measurements are currently performed at IFPEN on the Ecosural engine to provide LES dedicated data, but were not available for this study. For this reason, only qualitative comparisons between the different numerical set-up are achieved, with the aim to support or infer conclusions drawn on the F7P engine. Such results are nevertheless particularly interesting since the engine design and operating conditions (engine load, speed) are very different from those of the F7P engine.

Methodology
Since SI engine configurations can exhibit important levels of cycle-to-cycle variations, the best way to evaluate the above cited numerical schemes and turbulence models would be to perform multicycle computations and to compare statistical results over several tens of cycles. Unfortunately, the computational cost associated to the simulation of these numerous cycles prevents the use of this strategy for all the numerical tests. Thereby it has been chosen to perform a unique cycle calculation for each numerical scheme and turbulence model. Each computation starts at Intake Valve Opening (IVO) from the same initial conditions and ends after the combustion process. The comparison between the different set-up is performed at three levels of interest: trapped mass, flow field and combustion process. It is worth noting that this single-cycle strategy introduces a severe difficulty when comparing the results: it may be difficult to separate the differences due to a change in the numerical set-up from the differences due to "natural" cycle-to-cycle variations. In that sense, the objective here is not to establish a definitive hierarchy between the different numerical schemes and SGS models, but only to extract firsts trends and new elements for future computations and to point out some possible unphysical behaviors as well.

Engine Configurations
The first single cylinder configuration is the F7P engine [26,27]. It is fully equipped with sensors and optical accesses and benefits from a full experimental and numerical characterization on several operating points and a extensive database is then available for comparison. This naturally aspirated configuration consists in a single-cylinder four-valve spark-ignition engine fueled with gaseous propane (Fig. 1). Its main specifications are given in Table 1. The operating point chosen for this study is the one named unst_dil in Granet et al. [13]. This condition is called "unstable" because of its high degree of CCV (COV IMEP ¼ 7:2%, Tab. 2). Compared to a stable operating point, it is expected to exhibit larger differences between the different set-up.
As mentioned in Section 1.1, the comparisons of the different set-up are only based on single-cycle computations due to their high computational cost, which may complicate the analysis and the conclusions. In order to enforce (or not) these conclusions, a second test configuration is also studied. This second configuration is a highly downsized spark ignition engine (Fig. 2) recently developed at IFPEN for the ICAMDAC project [5]. It is equipped with direct injection and is characterized by a high tumble ratio aiming at generating important levels of turbulence in the combustion chamber. It can therefore support elevated boost levels and IMEP (Indicated Mean Effective Pressure) of the order of 30 bar. Table 3 summarizes the main specifications of the Ecosural engine. In order to bring complementary elements compared to the F7P case, very different operating conditions are voluntary chosen, namely both engine speed and load are increased (Tab. 4) to mimic near knocking conditions called "knock" in the following. In practice, gasoline and iso-octane will be experimentally tested in a) Sketch of the experimental F7P engine test bench and b) view of a typical tetrahedral mesh during the intake stroke [25]. this engine, but gaseous propane is used for this qualitative study in order to avoid fuel stratification effects and to facilitate comparisons with the F7P simulations.

Numerical Set-Up
All the computations are performed with the AVBP LES code [12,28], which solves the compressible multispecies Navier-Stokes equations. The energy deposition [29] and Thicken and Flame for LES (TFLES) models [6] are respectively used to simulate spark ignition and flame propagation, as in [13]. Simulation grids are made of tetrahedra, allowing to refine the mesh in specific areas such as the valve seats or the spark plug, while using coarsened meshes in the intake and exhaust pipes or plenum  ( Fig. 1, 2). Due to the piston movement, the number of cells highly evolves during the computation and ranges from about 2 millions at Top Dead Center (TDC) to around 10 millions at Bottom Dead Center (BDC) depending on the engine (Tab. 1, 3). Computational ressources needed to simulate one cycle are about 30 hours on 400 processors of a SGI Altix ICE 8200 for each configuration when using the Lax-Wendroff scheme and the Smagorinsky model.

HIGH ORDER NUMERICAL SCHEMES IMPACT ON LES OF ICE
In this section, the different numerical schemes (LW, TTGC and TTG4A) are evaluated on the unst_dil operating point of the F7P engine. In all computations, the SGS model is the constant coefficient Smagorinsky (Cs ¼ 0:18) and all other parameters (grid, combustion model, CFL, etc.) are unchanged. Results are compared to experimental data and to the reference numerical results (LW-Smagorinsky) reported in [13]. A qualitative study of the Ecosural simulations is also used to bring complementary elements to the F7P LES analysis. In order to identify possible improvements brought by these new numerical schemes, comparisons are first made on trapped mass, then on the flow properties during the intake and compression strokes and finally on the combustion process.

Trapped Mass
The first macroscopic quantity a piston engine computation should be able to predict is the mass trapped in the cylinder after IVC because it has a first order impact on the engine thermodynamic cycle. Most of the time the first cycle of a multicycle LES does not allow to get this quantity with precision due to the influence of initial conditions. The present F7P simulations get round the problem by starting from the end of cycle 21 of the unst_dil database [13], i.e. from fully realistic initial conditions. For the Ecosural engine, two cycles without combustion are computed using the LW scheme and the Smagorinsky model to generate initial conditions for the combustion simulations. Table 5 summarizes the computed trapped mass for the three numerical schemes. The differences are very limited (% 1%) and this quantity is thus almost insensitive to the numerical scheme.

Intake and Compression Aerodynamics
During the intake and compression strokes, high aerodynamic cycle-to-cycle variations were observed by Enaux et al. [17] for the motored engine case and Granet et al. [13] for the reactive cases of the F7P engine. As an illustration, Figure 3 shows in grey the statistical envelope of the x-velocity along the cylinder axisz for the 50 cycles reported in [13] at four different crank angles (see Fig. 1 for a definition of the axis). Here, the statistical envelope of a quantity Q delineates the zone where 95% of the cycles are included and is defined as Q mean ðtÞ AE 2r Q with Q mean the mean value of Q and r Q its standard deviation. The instantaneous profiles obtained with LW, TTG4A and TTGC are also plotted for the four distinct crank angles.
Whatever the crank angle, all TTG profiles differ from the LW ones but remain very similar and lie within the statistical envelope. This observation holds true for the two other velocity components (not shown). As well, kinetic energy evolutions are similar for all schemes as shown in Figure 4. In particular, the same rapid resolved energy drop at the end of the intake stroke and during compression is retrieved. It is particularly interesting to notice that the same observation can be made for both F7P and Ecosural configurations even if the flow energy levels differ a lot during the intake strokes.
At this point, it is of course impossible to claim that a numerical scheme is better than another: no pathologic behavior is noticed whatever the numerical scheme and  all the results remain within the range of cycle-to-cycle variations. However, little differences can be identified if a closer look is made on velocity fields (Fig. 5, 6).
Here again, the overall flow motion and the biggest structures resolution are similar for all schemes. Nevertheless, one should notice that the higher order feature F7P engine: x-velocity profiles along the cylinder axis for the LW, TTGC and TTG4A schemes at a) À280 CAD, b) À240 CAD, c) À180 CAD and d) À100 CAD. The experimental envelope is extracted from [6]. The abscissa z = 0 corresponds to the bottom of the cylinder head.  Ecosural engine: velocity fields at the cylinder center (y = 0) for the LW, TTG4A and TTGC schemes at a) À240 CAD, b) À180 CAD and c) À50 CAD.
of TTG schemes has a direct impact on the convection of the smallest structures. Where LW tends to smooth and dissipates these small structures, TTG schemes and especially TTGC are able of a better accuracy because of less dissipation and diffusion. This behavior can be retrieved on turbulent viscosity levels shown in Figure 7. During the intake and compression phases, a higher level of turbulent viscosity is identifiable for TTG schemes and mostly for TTGC even if the Smagorinsky SGS model is used in all cases with the same constant C s ¼ 0:18. These differences can then only come from variations in the Reynolds tensor used by the model i.e. in the flow resolution. It results that a higher turbulent viscosity can be linked to more velocity gradients and in this case to less dissipation/diffusion of the smallest scales.

Combustion Process
Regarding the combustion process, three phases may be distinguished, as illustrated in Figure 8. The first one is the free flame propagation period during which the flame kernel generated after spark ignition evolves without being directly constrained by wall effects. During this period, the flame expands rapidly due to the very low density of burned gases compared to the fresh mixture and its wrinkling progressively grows under the action of small scale vortices. The kernel is also convected by large scale motions leading to cyclic variations of its localization. This first phase is crucial since it plays a main role in the combustion event phasing in the cycle.
In the second phase of the combustion process, the flame starts interacting with the piston and cylinder head due to the very low height of the combustion chamber around TDC. The reaction zone then propagates towards the periphery of the cylinder while being strongly affected by confinement effects, which vary a lot from one cycle to another. Indeed, the beginning of this phase highly depends on the kernel convection during the free propagation period since the flame can be more or less moved near the walls. The last phase of the combustion process is characterized by flame extinction at the cylinder liner due to the decrease of both flame wrinkling and laminar flame velocity. Flameturbulence interactions are thus not of first order in this part of the combustion process, which is greatly piloted by heat losses and flow kinetic energy dissipation. The distinction between these three phases of the combustion stroke is very important because comparisons should be limited to the free flame propagation period only. One main reason for this statement is that LES analyses are only based on single cycle simulations in this study. Since a modification of the numerical scheme (or of the SGS model) leads to a slightly different engine cycle realization in terms of flow motion, flame kernel convection towards the piston or cylinder head is also affected and the walls influence may highly perturb the comparison during the second phase. In order to separate effects linked to flame-wall interactions from those directly related to the flow resolution quality itself, Mean turbulent viscosity evolution for the F7P a) and the Ecosural b) engines. For this reason, graphs related to combustion will only be analyzed from Ignition to TDC in the following. During this free expansion phase, while TTG schemes predict similar resolved flame surface, LW exhibit a surprisingly low resolved surface as shown in Figure 9a and the same behavior is also retrieved at the SGS level (Fig. 9b). As a result, the combustion process predicted by LW is very slow and the in-cylinder pressure is close to the motored one for this scheme as shown in Figure 10. For the Ecosural engine, the same trend is obtained as shown in Figure 11, i.e. the flame surface is higher with TTGC compared with LW even if the results are not as discriminating as for the F7P engine.
In order to understand the behavior of LW compared to TTG schemes and especially TTGC, the focus is now put on the TFLES model and the efficiency function used to account for the SGS combustion in this study.

TFLES and Efficiency Function
In the previous section, it was found that the two TTG schemes were able to predict a similar combustion, in good agreement with the experimental envelope in the  In-cylinder pressure for F7P engine. The experimental envelope is extracted from [13]. Ecosural engine: resolved a) and SGS b) flame surfaces.
F7P case. On the contrary, combustion predicted by LW was extremely slow which was noticed on resolved and SGS flame surfaces and confirmed on the in-cylinder pressure in both F7P and Ecosural cases. For all these tests, the TFLES approach [6] was used to model the combustion process. In this approach, the flame is artificially thickened in order to ensure an appropriate resolution of the flame front on the LES grid. However, when the flame is thickened, the combustion-turbulence interaction is affected, reducing flame wrinkling at the resolved level and obliging to model the lost part at the SGS level. This role is fulfilled by the efficiency function E, which is given by Equation (36) in the paper of Colin et al. [30]: where Nðd 1 l Þ and Nðd 0 l Þ are respectively the wrinkling factors corresponding to thickened and non thickened flames estimated by: with C a function of the SGS strain rate, u 0 the fluctuating velocity, D e the filter size, s 0 l the laminar flame speed and a a constant given by: In this expression, b is the TFLES model constant and c ms another constant fixed to 0:28 using the Direct Numerical Simulation (DNS) by Yeung et al. [31]. In this study where LW, TTG4A and TTGC are compared, all parameters are the same (mesh, chemistry, operating conditions, etc.). The only parameter which changes depending on the numerical scheme used is the fluctuating velocity u 0 since this velocity is computed using the resolved velocity field [6]: When using the TFLES model, users have to fix the b constant (Eq. 3) which is recommended to be of the order of unity by Colin et al. [30]. A classical value of 0.3 is used here for b according to CERFACS and IFPEN experience of high order schemes [31]. This value was kept for the three tested numerical schemes in order to ensure fair comparisons. If one wants the combustion process to be independent of the numerical scheme used, it obviously means that SGS contributions should exactly complement resolved contributions. In other words, if a numerical scheme allows to account for a larger part of the flame-turbulence interaction at the resolved scale, one expects the SGS part to be reduced accordingly. In Section 2.1.2, it was shown that, even if similar large scale flow motions were retrieved for the three schemes, a better resolution and finer structures were identified with TTG schemes and especially with TTGC. According to Equation (4), a better resolution induces a more intense fluctuating velocity and thus a higher efficiency according to Equation (1) and (2). As a consequence, if one compares TTGC and LW simulations, LW exhibits both lower resolved and lower SGS contributions compared to TTGC.
This unexpected effect is attributed to the construction of the efficiency function and the b constant determination based on very high order schemes and DNS results [6]. In such highly precise simulations, all the turbulent spectrum is resolved with high fidelity until the highest wave numbers associated to the filter size. Taking into account those resolved scales, the efficiency function is built to provide an estimation of the impact of the non-resolved scales (i.e. SGS) on the combustion process. With low order and dissipative schemes, all flow structures larger than the filter size are not precisely resolved, especially the smallest ones. It results that some information is missing for the efficiency function to estimate the SGS contribution. It means that for a given resolved flow, the efficiency function underestimates the SGS contribution if low order schemes are used. It suggest that the b constant can not be scheme independent and has to be adjusted depending on the numerical scheme used.
Actually, this efficiency function adaptation was performed in a very empirical way in several former computations on different engine configurations, in particular the F7P engine considered in this work and the XU10 engine described in [9,17] where authors found that b ¼ 2 gives very satisfactory results with LW. Indeed, Figure 12 shows that with this value, the expected behavior is retrieved: the resolved part is still lower with LW compared to TTGC or TTG4A but this loss is balanced by an increased contribution of SGS surface.
To get an idea of the combustion velocity, Table 6 presents the classical indicators CA2, CA50 and CA90, where CAX is defined as the timing for which the burnt mass fraction reaches X %. With b ¼ 2, the LW computation provides good estimations of the different combustion times compared to experimental observations. It can also be noticed from this table that the standard value of 0.3 with TTGC gives correct results.
Finally, it can be seen in Figure 13 that with b ¼ 2 the LW in-cylinder pressure gets back in the experimental envelope. However it is worth noting that, if the pressure traces are now similar with LW, TTGC and TTG4A, this result is clearly due to an increased modeling contribution for LW as shown in Figure 14. While the maximum SGS contribution never exceeds 12% for TTG schemes, the SGS flame surface represents up to 35% of the total flame surface for the LW simulation. In that sense, the use of higher order schemes such as TTGC seems highly preferable since it clearly allows to lower the impact of the modeling on the simulation. However and even if this conclusion is similar for the two engines tested here, this tendency should be confirmed on multicycle simulations.

Restitution Times
Restitution times for the different convection schemes are presented in Table 7. As expected, the high order of TTG schemes comes with additional simulation costs by a factor close to 2. As mentioned in the introduction, such an increase in the computational resources consumed by the simulation is a non negligible drawback because of the numerous cycles which have to be simulated.

EVALUATION OF SUB-GRID-SCALE MODELS
In this section, the constant coefficient Smagorinsky model [14] (Cs ¼ 0:18) is compared to its dynamic version [15] and to the sigma closure [22] for both F7P and Ecosural configurations. For CPU time reduction purposes, the LW scheme is used in this section with the efficiency constant b = 2. All other parameters (grid, combustion model, CFL, etc.) are kept unchanged.
Results are here again compared to the experimental data and to the reference numerical results (LW-Smagorinsky) reported in [13] for the F7P engine. As for the scheme comparison study, Ecosural simulations are also used to bring complementary elements. The comparison is made on the trapped mass, the aerodynamic fields and combustion process.

In-cylinder Trapped Mass
As shown in Table 8, for both engine configurations, the differences in trapped mass between the three SGS models are very slight, which means that the influence of the turbulence models on this very macroscopic quantity is almost negligible. Figure 15 shows instantaneous x-velocity profiles along the cylinder axis obtained with the constant Smagorinsky, dynamic Smagorinsky and sigma SGS models at four distinct crank angles as well as the experimental envelope extracted from [13] for the F7P engine. From the aerodynamic point of view, no clear trend appears when comparing the three simulations. All profiles are different, with variations reaching several meters per second, but they remain in the statistical envelope. For the two other velocity components (not shown) the same observation can be made. The analysis of the velocity profiles thus do not allow us to draw conclusions since variations remain in CCV envelopes. Cut-planes shown in Figures 16 and 17 permit a slightly different analysis. Large scale motions are very close for the three models and the overall resolved kinetic energy decays presented in Figure 18 are comparable all along the intake and compression strokes, but slight differences appear with the sigma closure. This    This behavior can be attributed to a lower turbulent viscosity level of the sigma model, as illustred in Figure 19. This low viscosity was expected since sigma has been built to avoid an over-estimation of the SGS turbulence in shear layers and in solid rotation zones such as tumble in this case [22]. All these statements are confirmed by the Ecosural simulations, which exhibit the same qualitative trends for all physical quantities although operating conditions and engine design highly differ from the F7P.

Intake and Compression In-Cylinder Flow
Regarding the dynamic Smagorinsky model, it should be noticed that the mean in-cylinder constant, presented in Figure 20, is very close to the reference value 0.18 of the classical formulation during the intake stroke and progressively differs along the compression, where it increases up to 0.22 for the F7P and 0.25 for the Ecosural due to the effect of the walls. The fact that this increase is lower for the F7P configuration explains why the two Smagorinsky models present more comparable velocity fields than for the Ecosural engine (Fig. 16, 17). In this last case, all the fields have indeed very different shapes. Therefore, it should not be considered as a general conclusion that Smagorinsky and dynamic Smagorinsky simulations result in a comparable resolved flow before ignition.    Mean turbulent viscosity in the cylinder for the three SGS models: a) F7P engine, b) Ecosural engine.
Another important finding of this study is that even if the SGS viscosity is lower for sigma than for the Smagorinsky models (Fig. 19), the corresponding resolved kinetic energy is not necessarily higher (Fig. 18). An explanation is that smaller structures induced by the use of sigma are more easily dissipated, while models characterized by a higher SGS viscosity may lead to more organized large scale motions with a greater life-time. In practice both phenomena (small scale and SGS dissipation) will compete and it is difficult to know a priori which one will dominate.
A last result concerns the dynamic Smagorinsky model, which exhibits very high levels of SGS viscosity at the walls compared to other closures, as shown in Figure 21. This observation remains true for both engines and is coherent with the mean in-cylinder constant level at the end of compression. Such a behavior is not physical and may be related to a lack of grid resolution at the walls. Its influence on combustion will be discussed in the following section.

Combustion Stroke
The flame surface with sigma is higher than those from the Smagorinsky models for the two engines as shown in Figures 22 and 23. This statement remains true at both resolved and SGS levels, even if the resolved surface is more important. This behavior may first be related to the velocity field resolution, which presents more and smaller turbulent structures for sigma. This closure also allows to reduce diffusive fluxes within the reaction zone, which potentially allows a better resolution of flame-turbulence interactions even if the velocity field is the same. This last statement was however not verified in this study.
A further analysis of the flame surface evolution of the F7P engine shows that the dynamic Smagorinsky model rapidly leads to resolved and SGS quantities close to    those of sigma as soon as the flame starts interacting with the walls. This may be linked to the high (and unphysical) values of SGS viscosity in these regions, implying a rapid increase of the SGS wrinkling through the efficiency function of the TFLES model. This tendency is not recovered with the Ecosural engine, for which the heat release remains weaker for the dynamic model than for sigma all along the cycle. Therefore, it should not be considered that these two models have in general comparable behaviors. Despite these discrepancies between the SGS models in terms of flame surface, the associated in-cylinder pressure curves exhibit limited differences in the free propagation phase and globally stay within the experimental envelop of the F7P all along the cycle, as shown in Figure 24. Only the sigma model leads to a slight overestimation of the pressure during flame kernel expansion, but its evolution is then close to the dynamic model one, which is coherent with previous findings on flame surfaces. Cylinder pressure curves from the Ecosural engine are not presented since no experimental data, and especially pressure envelopes, are available. This quantity thus does not bring additionnal information compared to flame surface evolutions.
To conclude this section, it can be stated that SGS models do not have a huge impact on global quantities such as the flow kinetic energy during intake and compression or the global heat release in the free flame propagation period. Sigma allows a slightly better resolution of flame-turbulence interactions and could be preferred to the classical Smagorinsky model. On the contrary, the dynamic model shows non-physical behaviors at the walls at the end of compression and during combustion. It may thus be avoided to guaranty that flame propagation is not perturbed by this phenomenon in the two last phases of the combustion process.

Restitution Times
In terms of numerical costs, the SGS turbulence models used have no major impact as reported in Table 9. Only a small increase of the order of 5% is found for the dynamic Smagorinsky model which was expected because of the extra calculation needed to compute the constant. The choice of the SGS model may thus be more guided by stability and physical behavior criteria than computational ressources considerations.

EVALUATION OF A NEW HYBRID APPROACH
Numerical schemes were evaluated in Section 2 and presented a limited influence on the flow evolution during the intake and compression strokes, even if TTG schemes allowed a slightly better resolution of the velocity field. This may partly be explained by the well-refined LES grids used in this study allowing to directly capture the main part of the flow energy on the mesh. This is also because large structures are found during intake and compression, and the tumble motion breakdown process, which produces small scale turbulence, mainly occurs at the end of compression and in the first instant of combustion. On the contrary, huge discrepancies were observed during combustion, where higher order schemes permitted a higher resolution of the flame wrinkling on the grid, due to an improved description of flame-turbulence interactions. Nevertheless, the computational cost of TTG schemes (almost twice the cost of LW) makes the use of these high order schemes unrealistic in practice. Starting from these conclusions, a natural idea is to split the engine cycle in two separate parts, namely intake-compression and combustion, and to use dedicated numerical set-up for each one to optimize both precision and CPU cost. Therefore, TTGC may be used during combustion, while LW may be retained for the rest of the cycle (intake, compression, end of expansion and exhaust).
However it is important to keep in mind that this hybrid approach is limited to homogeneous cases. In very heterogeneous configurations like direct injection of fuel or controlled auto-ignition, TTGC could be activated during engine compression to keep a good resolution of temperature and mixture stratifications, phenomena which are of first order in such cases.
The proposed hybrid approach is called ES O 2 (Engine Stroke Optimal Order), and is analyzed in the following using once again the two engine configurations presented in Section 1.2. For each engine, all simulations are performed keeping the same constant b ¼ 0:3 in the efficiency function of the TFLES model in order to perform fair comparisons. The SGS closure used in this part is the constant Smagonrinsky model.

Flow Motion During Combustion
For an ES O 2 calculation, TTGC is imposed only 2 CAD before ignition (i.e. from À52 CAD for the F7P engine and À22 CAD for the Ecosural engine). This period prior to ignition corresponds to a too short time interval for affecting the velocity field at ignition as shown in Figure 25 (respectively Fig. 26) where ES O 2 and LW velocity fields are very similar at À50 CAD (respectively À15 CAD) for the F7P engine (respectively Ecosural engine). Engine flow motion differences are thus only analyzed during the combustion phase in this section.    During the whole free propagation phase, a memory effect is clearly visible on the velocity field of the ES O 2 computation. Indeed, the ES O 2 field looks like the LW one during the first instant of the flame propagation and the TTGC simulation still presents more turbulent structures than the two other cases. During the flame growth, combustion has an important effect on the flow because of the fast thermal expansion of burnt gases and higher differences are obtained between LW and ES O 2 . At the end of the free propagation phase, close to TDC, the ES O 2 velocity field even presents small eddies and is more comparable to the TTGC one. This statement remains true for both engines and is confirmed by Figures 25 and 26 which show a slightly lower dissipation during combustion with the hybrid method than for LW and ES O 2 levels gradually join those of the full TTGC run.

Flame Propagation Process
Evolutions of the resolved and SGS flame surfaces are plotted in Figures 27 and 28 for the two engine configurations. During the free propagation phase, the resolved surface associated to ES O 2 is higher than the LW one and is more comparable to a full TTGC simulation. However, as described previously, the LW and ES O 2 cases have flow fields which differ from the TTGC one, especially close to the spark plug. Flame-turbulence interactions are thus deeply affected and comparisons between ES O 2 and TTGC can not be directly performed. The ES O 2 cycle may indeed be considered as a numerical perturbation of a full TTGC engine cycle during the first part of the cycle (intake and compression).
At the SGS level, TTGC and ES O 2 also give higher surfaces than LW. This is notably due to their higher resolved flame surface which is used to compute the efficiency level. Nevertheless, whatever the scheme used, the SGS surface stays lower than the resolved surface, confirming the good quality of the computational grid used in the two configurations.
The analysis of progress variable iso-surfaces based on fuel concentrations, plotted at different crank angles in Figures 29 and 30   field history effect at spark timing on flame-turbulence interactions. However, the propagation is faster for ES O 2 than for LW, and the flame expansion velocity is close to the TTGC case leading to similar cylinder pressure curve evolutions as shown in Figure 31. This observation corroborates the idea that the ES O 2 calculation can be view as a "perturbed" TTGC engine cycle.

CPU Time
All the achieved simulations include intake, compression and combustion phases. In other words, no more than about half a cycle is simulated (from À355 CAD to 50 CAD for the F7P configuration, and from À360 CAD to 88 CAD for the Ecosural configuration). In order to evaluate the potential CPU costs improvements associated to ES O 2 , an extrapolation is performed over a whole cycle assuming a TTGC-LW switch when the combustion ends, i.e. at 50 CAD (resp. 90) for the F7P (resp. Ecosural) configuration. The estimated restitution times are given in Table 10.
The CPU cost of ES O 2 is reduced by about 80 to 90% compared to TTGC due to the low combustion duration. Thus, the obtained CPU time remains in line with previous engine LES, while an increased precision can be expected. The proposed hybrid approach finally appears promising for LES studies of industrial configurations.

CONCLUSION
The objective of this study was to evaluate the feasability and the possible benefits of using high order schemes and state-of-the-art SGS models in LES of ICE. For this purpose a first part was dedicated to the comparison of two convective schemes from the Taylor-Galerkin family,     namely TTGC and TTG4A with the classical Lax-Wendroff one. Numerical findings indicate that all set-up lead to similar evolutions of the flow properties (kinetic energy, SGS viscosity) during compression and intake even if slightly more turbulent structures were found for TTG schemes. On the other hand, the free flame propagation phase exhibited completely different behaviors between the schemes. A good flame-turbulence interaction prediction with TTG schemes enables the use of the classical low efficiency constant (b ¼ 0:3) while the higher dissipation of LW had to be balanced by an increased of this constant (b ¼ 2). Because of this lower importance of the modeling contribution, TTG schemes are preferred to increase the precision in further studies. However, the counterpart of this increased precision is a higher restitution time (almost twice the cost of LW). In a second part, several SGS turbulence models were compared, namely the constant coefficient and dynamic Smagorinsky and the Sigma closure, keeping the Lax-Wendroff scheme for all simulations. It was shown that the sigma model allows a slightly better resolution of the velocity field because of a lower SGS vicosity without additionnal CPU time. On the contrary the dynamic Smagorinsky model generates high viscosity levels, especially at the walls, which is not a physical behavior, with a small increase (about 5%) of computational times. Nevertheless, discrepancies between the models were low and all quantities remained in the CCV envelopes, suggesting that the change in SGS closure mainly acts like a numerical perturbation of the flow.
Finally a new hybrid approach is evaluated to take advantage of the better precision of TTGC without being penalized by too large CPU times. This new approach, called ES O 2 (Engine Stroke Optimal Order) consists in using LW during intake and compression, where differences with TTGC are low, and TTGC during combustion to better resolve the flame. Results were very close to TTGC, with a reasonable increase of the CPU cost (about 10%) compared to LW. The fact that similar findings were obtained for both engine configurations, operating in very different conditions, may confer them a general character, even if statistics should be necessary to draw definitive conclusions. Finally, ES O 2 appears as a promising method for future LES studies of ICE. This approach may be ideally used with the sigma model. However, first tests indicate that combining sigma with TTGC may lead to numerical instabilities because of the very low levels of viscosity introduced in this case. As the constant coefficient Smagorinsky model gave results quite close to those of Sigma, it could also be retained for future practical applications in complex geometries.