© T. Nguyen and A.M. Kempf, published by IFP Energies nouvelles, 2017

Introduction

The extensive use of Computational Fluid Dynamics (CFD) in numerous engineering applications leads to many fundamental questions regarding the applied numerical approaches themselves, especially the effects of the numerical schemes on the simulation results, which have been investigated in many studies, as for examples [1-4]. Obviously, reducing the numerical errors is essential for an accurate solution of the transport equations of different physical quantities [5, 6]. It is clear that to obtain more accurate numerical results, fine grids and high order numerical schemes are preferable. However, high simulation costs due to finer computational grids as well as the numerical instability resulting from high order schemes are factors that must be considered. As a matter of fact, a compromise between accuracy, stability and convergence is often necessary. Depending on the nature of the problem, one may try and choose suitable numerical schemes that provide good results.

In order to investigate the influence of numerical schemes applied to engine simulations, several nominally 2nd-order numerical schemes are employed to examine their effect on the simulated in-cylinder flow field, as well as the combustion process. In Large Eddy Simulation (LES), investigators often outline the importance of the filter size and high-order discretization schemes to yield reasonable simulation results. However, it should be stressed that the order of numerical schemes is just one part of the story: It has been observed in many studies [6-9] that two schemes with the same nominal order of accuracy can yield very different results.

In the context of engine simulations, Misdariis et al. [10] investigated the influence of several numerical schemes for the discretization of the convective flux on the flow field and the combustion inside the engine. They observed a rather small difference between the in-cylinder flow fields resulting from the numerical schemes of second and third order of accuracy. In comparison to their findings, we observe a stronger impact of the applied numerical schemes on the small scale features of the simulation results, as illustrated in Figure 1. The Central Differencing Scheme (CDS), Quadratic Upstream Interpolation for Convective Kinematics (QUICK) and the Total Variation Diminishing (TVD) schemes yield significant differences in the velocity field where small structures are found in the whole flow field for the simulation using CDS with 2nd and 4th order, but not in the simulation with the TVD schemes. Interestingly, in spite of distinctive flow structures between different simulations, a rather small discrepancy of the total kinetic energy of the in-cylinder flow between the simulation results of all the employed schemes is achieved. This may not be overly surprising, since most of the turbulent kinetic energy exists on the larger scales. However, significant differences are observed on the smaller scales and hence the modeled sub-grid energy that has much impact on the wrinkling factor and turbulent flame speed.

In the scope of this paper, the commonly used LES quality indices are also discussed. Interestingly, more dissipative schemes (like TVD) tend to cause a smaller turbulent viscosity due to smaller velocity gradients, as the numerical scheme itself is more dissipative. One may argue that with a TVD scheme, much of the dissipation is achieved by the discretization, so that the turbulence model needs to dissipate less energy. This leads to a misleading situation where many commonly used LES quality criteria must fail. Criteria based on the amount of modeled turbulent kinetic energy [11, 12], as well as the criteria based on the ratio between the turbulent and the laminar viscosity (as proposed by Celik et al. [13]), will no longer be adequate.

In this work, simulations of a four-stroke engine [14] with two intakes and two exhaust valves are performed on two different computational grids with cell sizes of 0.8 mm and 0.5 mm. To avoid the complexity of the mesh generation in a moving engine-geometry, the motion of the piston and the valves are described by Lagrangian particles in combination with an Immersed Boundary Method (IBM) [15]. This mesh-free technique has been implemented in our code PsiPhi [16] and is simple and well suited for High-Performance Computing (HPC) applications. Flame propagation is modeled using a Flame Surface Density (FSD) model [17] with iso-octane as surrogate fuel. The flow conditions at the intake and the exhaust ports of the engine are computed according to Navier-Stokes characteristic boundary conditions [18]. Our implementation [19] of Nicoud’s Sigma model [20] is used to calculate the sub-grid scale stresses.

1 Numerical Methodology for Engine Simulation

The Favre-filtered compressible Navier-Stokes equations for mass, momentum and total energy are discretized and solved numerically on an equidistant Cartesian grid using the finite-volume method. To handle the moving parts of the engine such as the piston, the intake and exhaust valves, an approach by Nguyen et al. [16] using Lagrangian particles and an IBM [15] is applied for the description of the geometries as well as the interaction between the moving boundaries and the fluid. In this approach, the moving geometries are presented with a cloud of logic and massless particles instead of the unstructured cells from body-fitted methods. The harmonic motion of the piston and the motion of the valves from the valve lift profiles are governed by the transportation of different groups of Lagrangian particles, which are used to describe the moving objects. Extensive discussion about the advantages, limitations and implementation of the method can be found in [16].

2 Combustion Modeling with Flame Surface Density

The modeling of the combustion process is described by the following equations:$$\begin{array}{l}\frac{\mathrm{ }\mathrm{\partial }\bar{\rho }\tilde{c}}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\bar{\rho }\tilde{c}{\tilde{u}}_i}{\mathrm{\partial }{x}_i}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left[\bar{\rho }\left({\widetilde{\mathrm{cu}}}_i-\tilde{c} {\tilde{u}}_i\right)\right]=\overline{\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(\rho D\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_i}\right)}+{\overline{\stackrel{̇}{w}}}_c=\overline{\rho S_d}{\mathrm{\Sigma }}_{\mathrm{gen}}\\ \end{array}$$(1) $$\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left[\bar{\rho }\left({\widetilde{\mathrm{cu}}}_i-\tilde{c}{\tilde{u}}_i\right)\right]=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(\bar{\rho }\frac{{\nu }_t}{\mathrm{Sc}}\frac{\mathrm{\partial }\tilde{c}}{\mathrm{\partial }{x}_i}\right)+F_{\mathrm{cgt}}$$(2) $$\begin{array}{l}\overline{\rho S_d}{\mathrm{\Sigma }}_{\mathrm{gen}}\approx {\rho }_u{S}_{\mathrm{l}}{\mathrm{\Sigma }}_{\mathrm{gen}}={\rho }_u{S}_{\mathrm{l}}\mathrm{\Xi }\left|\nabla \bar{c}\right|-F_{\mathrm{cgt}}\approx {\rho }_u{S}_{\mathrm{l}}\mathrm{\Xi }\left|\nabla \tilde{c}\right|\\ \end{array}$$(3)

On the LHS of Equation (1), the sub-grid scalar flux is calculated by the classical gradient hypothesis and a counter gradient term F _cgt (Eq. 2). In Equation (3), the generalized flame surface Σ_gen is approximated by replacing the Reynolds-averaged progress variable by the Favre-averaged one , where the counter gradient transport F _cgt is implicitly included (see the Discussion by Ma et al. [21]). The two terms and D in Equation (1) are the mean reaction source term and the molecular diffusivity. In Equations (2)-(3), the turbulent viscosity, the turbulent Schmidt number, the laminar flame speed and the wrinkling factor are represented by ν _t, S _c, S _l and , respectively.

The influence of turbulence on the flame propagation is described by the wrinkling factor , which is modeled using an approach proposed by Muppala et al. [17].$$\mathrm{\Xi }=1+a{{\mathrm{Re}}_{\mathrm{t}}}^{0.25}{\left(\frac{u^{\mathrm{\prime }}}{S_{\mathrm{l}}}\right)}^b{\left(\frac{p}{p_0}\right)}^c$$(4)

The turbulent Reynolds number Re _t is a function of the sub-grid scale velocity fluctuation u′ [22], the cell size Δ and the laminar viscosity ν:$${\mathrm{Re}}_{\mathrm{t}}=\frac{u\mathrm{\text{'}\Delta }}{\nu }$$(5)

In Equation (4), the model parameters are set to b = 0.2, c = 0.2, and a = 0.46 for the iso-octane/air flame [17, 23]. In the scope of this paper, elaborated ignition modeling is not considered. Instead, we set the progress variable c to 1.0 (burned area) in a small region below the spark-plug so that the flame can start developing and propagating. It is clear that a proper ignition modeling can be done in more sophisticated ways. However, this would add additional uncertainty to the present study and be beyond the scope of the present study. We have therefore chosen to just show the effect of the numerical scheme. It is also clear that some shortcomings of numerical schemes can be partially compensated by suitable models or model constants. Further information on combustion modeling can also be found in our previous paper [16].

3 Experiment

In this study, we consider a four-stroke optical engine, which is operated by the Dreizler group at Technical University of Darmstadt [14]. The engine runs at 800 rpm with iso-octane (C₈H₁₈) as a surrogate fuel. Optical access is provided via a quartz glass cylinder liner and a flat piston window. The engine consists of a twin-cam, an overhead-valve pent-roof cylinder head with two intakes and two exhaust valves. The characteristics and operation point are shown in Table 1. For further details of the engine, the readers are referred to the paper of Baum et al. [14].

4 Numerical Setup

The LES in-house code PsiPhi [16, 24-27] is used to solve the Favre-filtered governing equations for a compressible fluid. An equidistant Cartesian grid is used throughout the whole computational domain, which provides good accuracy and avoids the numerical dissipation caused by deformed or stretched cells of an unstructured-grid. In an equidistant Cartesian grid, the numerical error directly results from the applied interpolation scheme without adding the truncation error from the irregular grid. The disadvantage of this approach is that a local refinement is not possible, which leads to an increase of the total number of cells in the computational domain.

Results from different numerical schemes are compared: a 2nd and a 4th order CDS (CDS2, CDS4), the Limited Linear scheme [8] (LL01), the QUICK scheme [28] and two TVD schemes using CHARM [29] (TVD-CH) and Van Leer [30] (TVD-VL) limiters. It is important to note that we use a value β = 0.1 for the LL01 scheme, which makes this scheme less dissipative and supposedly “close to CDS” (Tab. 2).

The TVD scheme with CHARM limiter [29] is used for the discretization of all scalar quantities including density , progress-variable and total energy .

An explicit third-order Runge-Kutta scheme is employed for time integration. The unresolved turbulent viscosity is determined with Nicoud’s σ model [20] with the model constant C _m of 1.5. In our experience, the combination of a sub-grid model, good equidistant Cartesian grids and an explicit low storage Runge-Kutta scheme with a small CFL number is sufficient to achieve stability – some upwinding is only introduced right at the walls and for increased Mach number. This approach has enabled stable – but not always easy to setup.

It is noted that the numerical analyses in the following sections are based on the phase-averaged data of five consecutive engine cycles for each numerical scheme. Nevertheless, we demonstrated some result from a chosen engine cycle to highlight the significant distinctions of the numerical results due to the different numerical schemes.

5 Kinetic Energy Analysis

Figure 1 shows the normal velocity components at −270 and −90 CAD, which result from six different numerical schemes. In comparison to the simulated results obtained by CDS2 and CDS4, where small structures are dominant in the flow field, much smoother flow fields are produced by the simulations using the TVD schemes (TVD-CH, TVD-VL, LL01). On the other hand, only big structures are visible for the simulations with the TVD schemes for the momentum at −90 CAD. It is noticeable that the flow structures obtained by QUICK are somewhat in between CDS and TVD. One should note that the use of non-central schemes for momentum transport in LES is a fairly recent phenomenon, which has emerged in recent years with the availability of “LES models” in robust mainstream codes like ANSYS [31], STAR-CCM+ [32] or OpenFOAM [33].

As observed in Figure 1, the simulations using CDS with 2nd and 4th order schemes resolve more small structures in the in-cylinder flow field than the other schemes. In order to quantitatively analyze the differences, we evaluate the integral of the kinetic energy, , the integral of the velocity gradient and the integral of turbulent viscosity inside the cylinder Ω. Figure 2 shows the comparison of the integrated kinetic energy KE _Ω , integrated velocity gradient and integrated viscosity between six numerical schemes for momentum. As it is shown in Figure 2a, the total resolved kinetic energy during engine operation is quite consistent between different numerical schemes, where the maximum value of KE _Ω is achieved around −270 CAD and gradually decreases towards the compression and expansion strokes. The spikes of KE _Ω observed at −125 CAD and 105 CAD result from closing the intake valves and opening the exhaust valves, which cause a high velocity in the small gap between the valves and port walls. Interestingly, it is hard to notice any significant difference in the integrated kinetic energy KE _Ω of the in-cylinder flow field between the simulations using different numerical schemes. Despite the difference in the flow structures, all the numerical schemes seem to perform very well in resolving the kinetic energy of the flow field. This may not be surprising, since most of the kinetic energy exists on the largest scales.

Since the wide dynamic range of data values obtained during engine simulation makes it difficult to distinguish the differences of KE _Ω , and resulting from different numerical schemes (Fig. 2a-2c), using the normalization of all simulated data to the corresponding data obtained with the CDS2 scheme is a convenient way to highlight the variation between the schemes. As it can be seen in Figures 2a and 2d, the resolved kinetic energy yielded by the six numerical schemes using 0.5 mm grid are not significantly different from each others, whereas on the coarse grid of 0.8 mm, the resolved kinetic energy is clearly lower. As seen in Figure 2d, during the intake stroke, the total resolved kinetic energy of the simulations using more dissipative schemes such as QUICK, TVD-CH, TVD-VL and LL01 is lower than the resolved kinetic energy obtained by simulations using CDS2 and CDS4. However, during the compression stroke, KE _Ω obtained by simulations using QUICK, TVD and LL01 seems to dissipate quite slowly in comparison to KE _Ω obtained by CDS2 and CDS4. This may explain why we see an increasing trend of KE _Ω from QUICK, TVD and LL01 with reference to KE _Ω from simulations using CDS. Despite the significant differences in the flow fields obtained by different numerical schemes (Fig. 1), the resolved kinetic energy KE _Ω of the in-cylinder flow shows no substantial difference for six studied numerical schemes. Obviously, the higher or lower value of resolved kinetic energy KE _Ω give no definite indication whether the applied numerical scheme is more or less accurate and dissipative than other schemes. Hence, it is questionable whether the total resolved kinetic energy can be used to determine the quality of LES – or what LES quality actually is.

6 Velocity Gradient and Turbulent Viscosity

We have shown that the resolved kinetic energy may not be adequate to differentiate the influence of numerical schemes on the in-cylinder flow field, even though clear differences between these schemes are observed in Figure 1. On the other hand, the integrated velocity gradients in Figures 2b and 2f show a clear distinction between numerical schemes. As expected, the s obtained by both CDS2 and CDS4 are comparable and consistently larger than s of QUICK, TVD-CH, TVD-VL and LL01. Intuitively one may conclude that the turbulence level of the in-cylinder flow field obtained by QUICK (Fig. 1) lies somewhere in between TVD and CDS schemes. This is also clearly shown in Figures 2b and 2e where . Since the calculation of sub-grid turbulent viscosity depends on the velocity gradient, obtained by CDS2 and CDS4 are, as expected, greater than from QUICK, TVD-CH, TVD-VL and LL01. As shown in Figure 2f, the total sub-grid turbulent viscosity obtained from the dissipative schemes TVD-CH, TVD-VL and LL01 is just half the corresponding turbulent viscosity from CDS2 and CDS4. The distribution of the turbulent viscosity in the tumble plane in Figure 3 reveals the substantial difference between the less and more dissipative numerical schemes, where large values of the turbulent viscosity are found in the simulations using CDS2, CDS4 and QUICK. This obviously calls into question whether it is reasonable to determine good or bad LES quality based on the sub-grid turbulent viscosity, since the less dissipative numerical schemes would always lead to higher sub-grid turbulent viscosity, hence implying that more turbulence is resolved and less modeled by the more dissipative schemes.

7 Turbulent Viscosity in Dissipative Schemes

To examine further the influence of sub-grid modeling in dissipative schemes such as LL01, TVD-CH and TVD-VL, a simulation without sub-grid modeling using LL01 for the discretization of the momentum is carried out. The main reason for this test is that the scheme is very popular in both LES [34, 35] and Monotone Integrated Large Eddy Simulation (MILES) [36] simulations. The numerical results obtained by two simulations with and without turbulent viscosity can be distinguished in Figure 2.

As shown in Figure 2, the differences in the resolved kinetic energy KE _Ω of the in-cylinder flow field as well as in the absolute velocity gradient with and without turbulent viscosity are insignificant. With LL01, the contribution of turbulent viscosity has become minimal and it is hardly noticeable. The turbulent viscosity ν _t depends on the velocity gradients which are better maintained with the CDS or QUICK schemes than with TVD or LL01 schemes. By destroying the velocity gradient, the influence of sub-grid modeling becomes negligible and the engine simulation with LL01 scheme even without using turbulent viscosity ν_t can still converge and produces similar results as the simulation with the LL01 scheme with a turbulent viscosity ν _t. Figure 4 demonstrates clearly the influence of sub-grid modeling on the simulations using CDS2 and LL01. This finding may imply that LES with a limited linear (or similar dissipative) scheme for momentum are actually MILES simulations (implicit LES) with an ineffectual turbulence model. In other words, such simulations might be “implicitly implicit LES”. In such cases, LES quality criteria based on the turbulent viscosity must fail.

8 How much Kinetic Energy is Contained in Small Flow Structures?

Despite the small structures are dominant in the flow field of the simulations using CDS2 and CDS4, the resolved kinetic energy obtained by these simulations is quite similar to the results obtained by TVD-CH, TVD-VL or LL01 schemes (Fig. 2a). This raises a question of how much kinetic energy is contained in small structures. In order to reasonably estimate the amount of kinetic energy containing in these small vortices, we employ a Gaussian filter with the filter size of 3Δ × 3Δ × 3Δ (Δ denotes the size of the equidistant cells) to filter the small structures out of the flow field. The coefficients of the applied Gaussian filter are presented as a 3 × 3 × 3 dimensional array:$$\frac{1}{64}\left(\left[\begin{array}{ccc}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right], \left[\begin{array}{ccc}2& 4& 2\\ 4& 8& 4\\ 2& 4& 2\end{array}\right], \left[\begin{array}{ccc}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right]\right)$$

The difference in kinetic energy ΔKE _Ω  = KE_original − KE_filtered between the original flow field and the filtered flow field is a good estimation of how much kinetic energy the small structures contain. Figure 5 shows the vertical velocity component and its corresponding high-pass filtered component from different numerical schemes at −180 CAD. The total kinetic energy and the difference ΔKE _Ω between the original and the filtered flow field are plotted in Figure 6. According to our method, the kinetic energy contained in the small flow structures, is less than 10% of the total resolved kinetic energy of the flow field. This finding is interesting because it can explain why a comparable total kinetic energy is obtained by simulations using different numerical schemes for momentum.

From Figures 1, 2 and 6, it is clear that the dominantly small structures in the flow field, obtained by simulations using CDS, comprise just a small amount of kinetic energy in comparison to the big structures which govern the flow field of the simulations using TVD or LL schemes. In spite of holding a limited amount of kinetic energy, the fluctuations in the flow field significantly contribute to the velocity gradient and eventually, to flame propagation during combustion.

9 Dissipation Rate

The resolved dissipation rate is evaluated by the turbulent viscosity ν _t and the phase-averaged strain rate tensor 〈S _ij 〉:$$\epsilon ={\nu }_{\mathrm{t}}\langle {\mathrm{S}}_{\mathrm{ij}}\rangle \langle {\mathrm{S}}_{\mathrm{ij}}\rangle$$(6)

Figure 7 shows the distribution of dissipation rate at −180 CAD in the tumble plane. As expected, the dissipation rate obtained by the simulations with CDS2, CDS4 and QUICK exhibits much higher values in comparison to the ones achieved with TVD-CH, TVD-VL and LL01. The higher dissipation rate resulting in simulations using CDS or QUICK schemes is a quantitative indication of how well the small structures are resolved. It can be clearly seen in Figure 7 that simulations using dissipative schemes (TVD or LL01) can only resolve small structures where strong turbulence exists.

The averaged dissipation rate of the in-cylinder flow for different numerical schemes is shown in Figure 8. The averaged dissipation rate from less dissipative schemes like CDS2 and CDS4 are almost two times larger the dissipation rate resulting from simulations using TVD or LL01 schemes.

10 Estimation of the Integral Length Scale

In order to estimate the integral length scale of the in-cylinder flow, one can compute the spatial correlation function (Eq. 7), in which the two-point correlation is evaluated along two chosen vertical and horizontal sampling lines (Fig. 9). We do not show the integral length-scales themselves as the number of cycles is not sufficient to determine the length-scale with sufficient precision. We show the differences in the correlation functions from which the length scales are determined instead.

The spatial correlation function R of the velocity fluctuation u′(x) reads:$$\mathrm{R}\left(r\right)=\frac{\langle u^{\mathrm{\prime }}\left(x\right)u^{\mathrm{\prime }}\left(x+r\right)\rangle }{\langle u^{\mathrm{\prime }}(x{)}^2\rangle }$$(7)

The velocity fluctuation u′(x) is calculated from the instantaneous velocity u(x) and the phase averaged velocity 〈u(x)〉:$$u^{\mathrm{\prime }}\left(x\right)=u\left(x\right)-\langle u\left(x\right)\rangle$$(8)

The integral length scale L is calculated by integrating R(r) along the sampling line up to the intersection point P between the curve R (Fig. 10) and the horizontal axis (s = OP):$$L={\int }_0^s\mathrm{}\mathrm{R}\left(r\right)\mathrm{d}r$$(9)

As shown in Figure 9, the two-point correlation function is computed along two sampling lines using two velocity components: R _xx is computed along the vertical line using the horizontal velocity u and vice versa. The idea is simple, a change in sign of the correlation R _yy along the horizontal sampling line or R _xx along the vertical sampling line can be considered as a very coarse estimation of the radius of the vortex [37]. More detailed explanation and discussion of the method can be found in the papers [12, 38].

Figure 10 shows the calculated two-point correlation function R(r) along two chosen sampling lines. To avoid confusion by plotting too much data (six schemes and 5 crank angle degrees), only results obtained from two simulations with CDS2 and TVD, are presented. Based on the Equation (9) and the comparison of the two-point coefficient curves R _xx and R _yy between CDS2 and TVD-VL, it would be reasonable to consider that the length scales increase with more dissipative schemes.

11 Influence of Numerical Schemes on Combustion

Figure 11 shows instantaneous volume renderings of the reaction source terms at different crank angle degrees. Since the volume rendering is quite similar between the simulations using CDS2 and CDS4 and between the simulations using TVD-CH, TVD-VL and LL01, we present only the visualization of flame propagation from simulations using CDS2, TVD-CH and QUICK instead of all six simulations. As evidenced by Figure 11, from −8 to 8 CAD, the growth rate of flame area from the CDS2 and QUICK schemes are quite comparable. The flame-wrinkling produced by these less dissipative schemes is much stronger in comparison to the wrinkling from simulations using TVD. Although the total kinetic energy produced by the six numerical schemes is obviously in agreement (Fig. 12), the obtained burning rates behave very differently. This is due to the fact, that much of turbulence is still very well “maintained” by the simulations using CDS in comparison to the ones using TVD schemes, which can be clearly demonstrated by the vertical velocity component at −16 CAD (Fig. 13). Evidently, the integrated turbulent viscosity achieved with the simulations using CDS is more than 40% higher than the turbulent viscosity ν _t from the simulations using TVD, which directly results in a larger discrepancy in the integrated reaction source term and the in-cylinder pressure. Interestingly, the difference in between the simulation using CDS and QUICK (less than 20%) does not lead to strong deviation in the flame propagation and the in-cylinder pressure, especially between the simulated results from CDS2 and QUICK up to 6 CAD.

12 Les Quality Criteria

Determination of good or bad LES quality with a given resolution and numerical scheme has been extensively discussed in many studies [13, 38-42]. Focusing on engine simulations, some authors [12, 37] have applied several LES quality indices [11, 13]. These LES indices are used to estimate, whether a given resolution is sufficient for a “good LES”. However, most of these criteria are likely to fail due to the influence of numerical schemes on the numerical results. We have shown that simulations using highly dissipative schemes could be considered as a “good LES” since they seem to achieve less unresolved kinetic energy or less turbulent viscosity in comparison to the simulations using less dissipative schemes such as CDS2, CDS4 or QUICK. Similar arguments can also be found in the work of others [39, 43]. Obviously, any LES quality criteria, that consider the minimum of turbulent viscosity or unresolved kinetic energy as an indicator for a “good LES”, should be taken with care. To demonstrate the inadequacy of such LES quality criteria, we apply them to the LES engine simulations on the same grid using different numerical schemes.

The first criterion is the ratio of resolved turbulent kinetic energy k _res (x, t) to the amount of total turbulent kinetic energy (k _res (x, t) + k _sgs (x, t)).$$M\left(x,t\right)=\frac{k_{\mathrm{res}}\left(x,t\right)}{k_{\mathrm{res}}\left(x,t\right)+k_{\mathrm{sgs}}\left(x,t\right)}$$(10)

If a LES computation is able to resolve more than 80% of the total kinetic energy (M (x, t) > 0.8), then it is supposed to be a “good LES”. The calculations of resolved kinetic energy k _res (x, t) are outlined below.$$k_{\mathrm{res}}(x,t)=\frac{1}{2}(U_i-\langle U_i\rangle {)}^2$$(11)where 〈U _i 〉 is the phase-averaged velocity component. The sub-grid kinetic energy k _sgs (x, t) can be modeled from the phase averaged strain rate tensor 〈S _ij 〉² according to Equation (12) [12, 44] with the constant C _i set to 0.202:$$k_{\mathrm{sgs}}\left(x,t\right)=2C_i{\Delta }^2{\langle S_{\mathrm{ij}}\rangle }^2$$(12)

Figure 14 shows the distributed value of M (x, t) as computed from the in-cylinder flow at −270, −180, −90 CAD. As discussed, six investigated schemes lead to results that satisfy the “good LES criterion” related to resolved kinetic energy. Obviously, from the PDF of M (x, t) in Figure 14, the number of points exceeding the value of 0.95 by numerical schemes like TVD and QUICK is greater than for CDS2 and CDS4. This implies that the less dissipative schemes like CDS and CDS4 resolve a smaller part of the turbulent kinetic energy in comparison to TVD, which is clearly not true.

A second common criterion for LES quality is based on the ratio of the effective viscosity to the molecular viscosity . A good LES is supposed to obtain a value of IQ _ν above 0.8, while an IQ _ν value of 0.952, is considered a fully resolved Direct Numerical Simulation (DNS) [13, 43].$$\mathrm{I}{\mathrm{Q}}_{\nu }={\left(1.0+0.05 {\left(\frac{\bar{\nu }+{\bar{\nu }}_{\mathrm{t}}}{\bar{\nu }}\right)}^{0.53}\right)}^{-1}$$(13)

According to the PDF of LES IQ _ν at −270, −180 and −90 CAD in Figure 14, six engine simulations with different numerical schemes satisfy the LES IQ _ν index with the value of IQ _ν above 0.8 for all the computational cells. Similar to the energy criterion, a larger number of computational cells from simulations using dissipative schemes such as TVD and LL01 achieves a much higher value of LES IQ _ν in comparison to the LES IQ _ν obtained by simulations with CDS and QUICK. Figure 15 shows the averaged values of the two LES quality criteria, in which simulations using TVD schemes always show a better “LES quality” than the simulations with CDS or QUICK schemes. As stressed in Sections 7 and 8, it is not a good approach to use a quantity, that is strongly influenced by the numerical schemes, to evaluate the quality of a simulation.

Conclusion

This work studies the numerical effects on the flow behavior as well as the prediction of flame propagation in the context of engine simulation using LES. To carry out the investigation, six different numerical schemes (CDS2, CDS4, QUICK, TVD-CH, TVD-VL, LL) were used to discretize the momentum equation over multiple engine cycles using a real engine geometry [14].

As it was expected, small structures were found to be dominant in the flow field from simulations using less dissipative schemes (CDS2, CDS4 and QUICK), while big and smooth structures prevailed in the flow field from simulations using more dissipative schemes (TVD and LL). Interestingly, it is found that the total kinetic energy of the in-cylinder flow of all the simulations (using less or more dissipative schemes) was comparable from intake to exhaust stroke, despite much difference observed in the flow fields.

In this study, we have shown that the contribution of sub-grid modeling may be negligible for simulations using dissipative schemes such as LL or TVD. Since the velocity gradient is not well-resolved by the more dissipative schemes, the calculated turbulent viscosity is much smaller than the one obtained by simulations using less dissipative schemes.

Our study also showed that the numerical schemes have a strong effect on the combustion process since flame propagation is strongly influenced by the level of turbulence that must be well-maintained throughout the cycle. Obtaining a correct wrinkling factor or flame propagation requires a well-resolved flow field with sufficient turbulence.

We also confirm that many common LES quality criteria are inadequate to determine a “good LES”. We have illustrated that simulations using highly dissipative schemes are likely to be considered “good LES” under these criteria, which is clearly not correct.

References

All Tables

All Figures