Regular Article
LES of Internal Combustion Engine Flows Using Cartesian Overset Grids
^{1}
Institute for Combustion Technology, RWTH Aachen University,
Templergraben 64,
52062
Aachen  Germany
^{2}
Department of Mechanical Engineering, Sogang University,
Mapo,
Seoul
121742  Korea
^{3}
Honda R&D Co., Ltd., Automobile R&D Center,
141 Chuo, Wakoshi,
Saitama
3510193  Japan
^{*} Corresponding author email: t.falkenstein@itv.rwthaachen.de
Received:
18
January
2017
Accepted:
17
July
2017
Accurate computations of turbulent flows using the LargeEddy Simulation (LES) technique with an appropriate SubFilter Scale (SFS) model require low artificial dissipation such that the physical energy cascade process is not perturbed by numerical artifacts. To realize this in practical simulations, energyconserving numerical schemes and highquality computational grids are needed. If unstructured meshes are used, the latter requirement often makes grid generation for complex geometries very difficult. Structured Cartesian grids offer the advantage that uncertainties in mesh quality are reduced to choosing appropriate resolution. However, two intrinsic challenges of the structured approach are local mesh refinement and representation of complex geometries. In this work, the effectiveness of numerical methods which can be expected to reduce both drawbacks is assessed in engine flows, using a multiphysics inhouse code. The overset grid approach is utilized to arbitrarily combine grid patches of different spacing to a flow domain of complex shape during mesh generation. Walls are handled by an Immersed Boundary (IB) method, which is combined with a wall function to treat underresolved boundary layers. A statistically stationary Spark Ignition (SI) engine port flow is simulated at Reynolds numbers typical for engine operation. Good agreement of computed and measured integral flow quantities like overall pressure loss and tumble number is found. A comparison of simulated velocity fields to Particle Image Velocimetry (PIV) measurement data concludes the validation of the enhanced numerical framework for both mean velocity and turbulent fluctuations. The performance of two SFS models, the dynamic Smagorinsky model with Lagrangian averaging along pathlines and the coherent structure model, is tested on different grids. Sensitivity of pressure loss and tumble ratio to the wall treatment and mesh refinement is presented. It is shown that increased wall friction introduced by applying a wall model is overcompensated by some secondary effects, which lead to an overall reduction of pressure loss in the investigated engine geometry. Finally, dynamics of the statistically stationary valve jets are analyzed using Proper Orthogonal Decomposition (POD). Two distinct flow patterns are identified and the relevance for CycletoCycle Variations (CCV) is discussed.
© T. Falkenstein et al., published by IFP Energies nouvelles, 2017
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Introduction
For many practical combustion systems, LES of turbulent reacting flows seems to be an attractive technique compared with the ReynoldsAveraged NavierStokes (RANS) approach, since mixing fields can be computed with high accuracy. This capability is particularly interesting for stochastic effects that are relevant, for instance, for emissions formation or abnormal combustion phenomena. Recently, special interest has been attributed to CCV in SI engines [1]. Although LES is a suitable tool to predict these effects, fully reactive and even motored simulations over multiple cycles which achieve reasonable agreement with experimental data, are still scarce. Possible reasons are the challenges associated with performing accurate LES of high Reynolds number flows in the presence of complex geometries and moving walls.
The discrete representation of complex geometries can be naturally accomplished with unstructured grids. This approach is realized in the AVBP code, which has been used in several engine studies, e.g. [2–4]. Another wellknown example is the OpenFoam code, which has been used to carry out engine LES more recently [5]. However, meshing efforts are usually nontrivial. On the other hand, some engine LES studies have been performed based on structured Cartesian grids [6–8]. Major advantages of this approach are applicability of higher order numerical methods, reduced uncertainty with respect to LES mesh quality, insignificant meshing effort and computational performance. However, intrinsic challenges are the treatment of complexshaped walls and local mesh refinement. For the former, commonly used techniques are IB or Cut Cell methods [9]. Due to the wide range of Reynolds numbers present throughout the engine simulation domain, it is desirable to vary the grid spacing and make LES more affordable. Cartesian grid refinement can be realized by overset or Chimera methods [10], or blockstructured Adaptive Mesh Refinement (AMR) [11].
In this work, the efficacy of such approaches for LES on Cartesian grids is investigated in engine flows. Compared to a previous study [7], an IB and an overset grid method have been introduced into an inhouse code. Numerical accuracy near the walls has been improved and a wall function has been coupled to the IB method. The overset grid method is verified on a simple flow problem with analytic solution. For validation, an engine flowbench configuration is considered as representative test case for intake stroke flow conditions. The absence of moving walls reduces the computational cost to obtain incylinder flow statistics and simplifies the analysis of complex dependencies between flow features in the intake port and further downstream. Hence, it is a suitable case to evaluate the effects of different combinations of numerical methods and grid resolutions on the flow field.
The paper is organized as follows. First, the governing equations and baseline numerical framework are described. Second, recent enhancements in terms of wall treatment and mesh refinement are discussed and verified. Third, validation results and sensitivities for different SFS models, grid resolutions and numerical methods are discussed. Finally, dynamic behavior of the valve jets is analyzed.
1 Computational Framework
1.1 Governing Equations
Typical Mach numbers during intake stroke suggest the occurrence of compressibility effects, i.e. the computational framework should be suitable for at least weakly compressible flows. Here, the fully compressible NavierStokes equations are solved:(1)(2)where the viscous stress tensor is given by(3)Here, constant molecular viscosity μ is considered a reasonable approximation for the flow configurations discussed below.
The energy variable used in this work is total energy(4)where the thermal conductivity λ is computed from the molecular viscosity, assuming Pr = 0.7.
The system of equations is closed by using the ideal gas law as equation of state:(5)
For this study, calorically perfect gas is assumed, i.e. specific heat capacities are constant in all simulations.
1.2 Numerical Methods
In the following, discretization schemes and boundary conditions for the solution of the governing equations are described. Within this work, an existing numerical framework has been extended by three methods. First, the spatial filtering scheme [12], which removes numerical artifacts from the solution, has been modified according to [13]. Consequently, interior and nearboundary points are filtered with the same order of accuracy. Second, an overset grid method has been introduced. Third, an immersed boundary method [14] was added to represent the geometry. All methods have, but are not limited to, secondorder accuracy.
1.2.1 Discretization
The NavierStokes equations are discretized by secondorder finite difference schemes on a staggered mesh as described in [15]. Metrics operators at interior points do satisfy discrete conservation of mass, momentum, and kinetic energy, while at boundaries, only primary conservation is achieved. For the energy equation, a thirdorder accurate Weighted Essentially NonOscillatory (WENO) scheme [16] is utilized. To avoid inconsistencies between discrete mass conservation and scalar convection, which may result from combinations of higherorder finite difference and upwindbiased schemes, a finitevolume mass flux is used in the continuity and scalar transport equations [17].
The nondissipative nature of the spatial discretization used for the NavierStokes equations may require additional procedures to remove spurious oscillations from the solution. Such numerical artifacts might result from finite errors due to explicit time integration of the continuity equation in compressible flows, reduced accuracy near boundaries, discretization errors on nonuniform grids, or dispersive errors. These ‘wiggles’ can be eliminated by locally dissipative numerical schemes, explicitly adding artificial dissipation to the governing equations in parts of the domain, or applying lowpass filters to the solution as a postprocessing step [18,19]. Effective sensors for highfrequency modes are required for the former two approaches, particularly when turbulent flow structures are to be retained. If lowpass filters with sharp spectral cutoff properties are applied, sensors may not be needed. In this work, sixthorderaccurate spectrallike filters according to [13] are used for enginetype flows.
For time integration, an explicit scheme based on a secondorderaccurate, LowDissipation and LowDispersion RungeKutta (LDDRK) method is used [20]. In these types of schemes, the coefficients are optimized to reduce the dissipation and dispersion errors in acoustic wave propagation calculations. This allows for larger time steps than in classical RungeKutta methods to achieve the same accuracy. LDDRK schemes typically require a larger number of stages than the order of accuracy might indicate. Due to this property, such schemes can be implemented with lowstorage requirements [21], which is an attractive feature for LES since computational performance of fluid dynamics applications is mostly bound by memory bandwidth [22]. In all compressible flow simulations discussed below, the secondorderaccurate fivestage scheme (LDD25) [21] is used.
1.2.2 Immersed Boundary Method
In the present work, the Immersed BoundaryApproximated Domain Method (IBADM) as proposed by Kang et al. [14] has been integrated into the inhouse code. The basic concept of enforcing the velocity boundary condition at the actual geometric location of the interface follows the work of Fadlun et al. [23]. In this approach, the discretized equations in the direct vicinity of the immersed boundary are replaced by appropriate numerical boundary conditions, which are interpolated from the nearwall flow solution and the known physical boundary conditions on the immersed interface. This discrete forcing by direct boundary condition imposition from reconstruction eliminates stability issues encountered when continuous or feedback forcing methods are applied to rigid interfaces [9]. An example of the approximated fluid domain where the governing equations are solved is given in Figure 1 for a curved immersed boundary. In this work, simple linear interpolation is used, unless a velocity boundary condition is provided from the wall model (Sect. 1.4.1).
Challenges related to mass conservation arise when boundary conditions are imposed at the closest grid point on the fluid side of the immersed interface. In order to only solve the governing equations on the fluid side of the immersed boundary, and avoid complex treatment to enforce continuity locally in cells crossed by the interface [24], global mass conservation across the whole boundary is used as additional constraint:(6)
This is incorporated into the facevelocity boundary condition as follows. First, a temporary velocity is computed by interpolation from N_{nb} fluid points to the numerical boundary location Γ_{a}, using the velocity of the immersed boundary. Second, a correction is applied to the interpolation coefficients , which are functions of the local immersed boundary geometry and velocities. The final velocity boundary condition at each face m is computed in a local coordinate system, centered on the immersed boundary (no summation over repeated indices):(7)The correction to the interpolation weights is obtained by solving a constrained least square problem to minimize the amount of correction. Note that this redistributes the error in over the numerical boundary Γ_{a}. In flowbench simulations discussed below, the mass flux error to be compensated at each valve was found to be on the order of one percent of the mass flow passing the valve. After applying the correction procedure once, the mass flux error was reduced below 10^{−5}, which could be further decreased by repeated execution. Grid refinement can be used to reduce the amount of correction according to the order of accuracy of the chosen interpolation scheme.
In compressible flow simulations with Neumann boundary conditions for temperature and density, overheating or overcooling effects were observed in some cases near the immersed boundary. Locally, both variables were found to diverge into opposite directions, while maintaining constant pressure. To address this issue, the isobaric fix proposed by Fedkiw et al. [25] is used to calculate temperature and density at the numerical boundary.
Figure 1 Approximated flow domain with immersed boundary. 
1.2.3 Overset Grid Method
The method can be structured into several interdependent procedures. These are, in execution order, overset grid generation and definition of intergrid boundaries, hole cutting, intergrid interpolation and communication setup, as well as coupled time integration. Each step is described in the following.
Cartesian grids are fully parameterized by the bounding box dimension and the grid point distribution in each of the three coordinate directions, which makes grid generation very simple and timeefficient. As preprocessing step, the user first defines a background mesh for the entire flow domain, which can, but does not need to be generated. Overset grids can be defined in arbitrary locations by respective bounding boxes and a refinement/coarsening factor relative to the background mesh. In this way, the flow domain can be composed of a chain of overlapping grids and the base grid can be omitted, if desired. Grid overlap is automatically adjusted. Currently, all grids need to have the same orientation in the global coordinate system. To reduce interpolation errors and for ease of use, integer numbers are used as refinement factors in the present work.
Thus far, Cartesian blocks have been defined without considering the complex geometry. The employed IB method requires a procedure that tags all cells on the fluid side of the actual wall geometry, which is imported into the inhouse code in STereoLithography (STL) format. Cells cut by the geometry are identified by ray tracing, which is conceptually similar to the work in [26]. For this operation, grid and triangulated geometry are mapped to integer space which avoids precisionrelated ambiguities [27] and increases efficiency. Usage of Alternating Digital Trees (ADT) [28] further speeds up the intersection search, which is a desirable feature when applied to handle moving boundaries in the flow solver. The raytracing method is very robust since it is capable of representing the shell of even very thin geometry objects on the grid.
A threedimensional stairstep representation of the actual geometry is determined with a floodfill procedure, starting from userdefined material points on each grid. This geometry recognition approach can be efficiently employed to generate refined boundary layer grids with minimum user input, i.e. specification of the desired wall and refinement layer thickness (Fig. 2).
During the preprocessing step, intergrid coupling boundaries are defined as requested by the user. Finally, the parallel block decomposition as presented in [7], can be applied to all overset grids independently for different numbers of processors to conclude the preprocessing step.
In the flow solver, it is desirable to integrate the discretized equations on the finest available overset grid, with minimum mesh overlap. To accomplish this, holes are cut into coarse mesh regions. The cutting procedure is initiated by comparing resolutions of all available grids at each intergrid coupling boundary. On the finest mesh, the inhouse code internally creates an STL representation of the shape of the coupling boundary, but offset towards the fluid side by a userprescribed number of overlap cells. This internal geometry object needs to be made invisible to the fine mesh, but visible to all coarser meshes. Then, the geometry recognition procedure described in the previous paragraph can be utilized to cut holes into the mesh and create additional intergrid coupling boundaries.
The overlap distance separates locations where interpolation from the coarse to the fine and from the fine to the coarse grid solution is performed. For small overlap, the value interpolated onto one grid appears in the interpolation stencil for the other grid. In such a case, a coupled system needs to be solved which is called implicit interpolation [29]. In this work, grid overlap is chosen sufficiently large so the interpolation is solely explicit.
Once all intergrid coupling boundaries are geometrically defined, adjacent fringe or target points need to be identified in the hole or outer region, where the standard discretization schemes need information from other grids. The actual fringe point number and position depends on the spatial scheme and storage location of each variable. On the staggered mesh with different numerical schemes for velocities and scalars, four different sets of target points need to be considered. For each fringe point, an efficient point searching algorithm based on ADT finds the corresponding location on the donor or source grid, from which the flow solution is to be interpolated. At these donor points, interpolation weights are precomputed during initialization or, for moving boundaries and translated grids, before every time step. Arbitrary order of the polynomialbased interpolation scheme can be requested by the user. According to [29], the order should be one order higher than the order of the spatial discretization scheme.
During time integration with an explicit method, the major differences compared to a singlegrid solver are additional interpolation and communication procedures to be executed after the intragrid updates of variables across processors. One simplification introduced for computational efficiency concerns the dynamic procedure to compute eddy viscosity and diffusivity. To avoid communication of velocity and scalar gradients at intergrid coupling boundaries before the test filter is applied, Neumann boundary condition is assumed for the gradients. Effects on the SFS model results are expected to be very small. Currently, equations on all grids are integrated with the same time step size, which will be changed in the future.
Figure 2 Overset grid for boundary layer refinement in an intake port. 
1.3 Boundary Conditions
At open boundaries, i.e. where mass flux can enter or leave the flow domain, NavierStokes Characteristic Boundary Conditions (NSCBC) in threedimensional formulation [30] are applied. At inlets, no artificial turbulence is introduced, but a mean velocity profile is used as relaxation target. At outlet boundaries, local reverse flow is enabled by switching to the inlet NSCBC formulation in computational cells where mass flux into the domain occurs.
1.4 Physical Models
Details of the recent implementation of a wall model are given below, followed by a brief summary of the SFS models used for the computations of engine port flows.
1.4.1 Wall Model
Due to the high Reynolds numbers expected in engine flows, it is not feasible to fully resolve the boundary layers, even with local mesh refinement. For explicit time marching schemes, the time step would become prohibitively small. Hence, wall models are needed to account for the correct shear stress. While this applies to computations on both bodyfitted and nonbodyconforming grids, the latter approach is even more dependent on a wall function, if a boundary condition away from the wall needs to be known. In the present work, the IB method provides the nearwall velocity by linear interpolation. In regions of developed turbulent flow, the velocity estimate is replaced by loglawbased interpolation, as described in [31]. Underresolved velocity gradients, and accordingly underpredicted wall shear, are compensated by a model viscosity applied to the first cells off the boundary. Two modifications to the model have been introduced within this study. First, instead of using the instantaneous offwall velocity to predict the wall shear stress from the loglaw, the timeaveraged velocity field is used. This seems more meaningful from a physical point of view and is favourable in terms of robustness, as mentioned in [32]. For transient flows, the timeaveraging can be replaced by a timefiltering procedure, as suggested by [33]. The second modification concerns situations, where the boundary location needs to be prescribed at a location inside the viscous sublayer. Since the eddyviscositybased model is not applicable to such cases [31] and the actual wall distance is rather small, a model designed for bodyfitted grids [32] is locally used.
To validate the wall model in a canonical flow configuration, turbulent channel flows at Re_{τ} = 2000 were computed and compared against Direct Numerical Simulation (DNS) data [34]. In this case, the immersed boundary object representing the physical wall Γ_{ib} is parallel to the approximated boundary Γ_{a}. The distance between Γ_{ib} and Γ_{a} was varied between y^{+} = 0 (bodyfitted) and y^{+} = 50. In Figure 3, mean and Root Mean Square (RMS) values of the streamwise velocity component are plotted as function of wall distance. Except for the first cell next to Γ_{a}, the results are in very good agreement with the DNS data, irrespective of the distance to Γ_{ib}.
Figure 3 Mean streamwise velocity (top) and streamwise velocity fluctuation (bottom) as function of wall distance. 
1.4.2 SFS Models
In this study, two eddyviscositybased SFS models are tested in enginetype flows. The baseline model employed is the Dynamic Smagorinsky Model with Lagrangian averaging (DSML) along particle trajectories [35]. A comprehensive overview [36] of several SFS models with respect to the mathematical and physical nature of the NavierStokes equations has shown that the DSM features multiple desirable properties from a theoretical point of view. Here, the DSML is compared against the more recent Coherent structure Smagorinsky Model (CSM) [37]. To the author's knowledge, this comparison has not been done for practical applications.
Both approaches eliminate one of the key shortcomings of the static Smagorinsky model [38], namely the flowproblemdependent selection of a constant model parameter. Instead, the extended models determine the subfilter length scale from local resolved scales. However, the considerations for deriving the two models are quite different. While the dynamic procedure is based on a complex test filter operation and scale similarity assumption, the CSM utilizes a correlation between the second invariant Q of the velocity gradient tensor and turbulent kinetic energy dissipation. Although this correlation is rather weak for LES, promising results for canonical flows were reported [37,39]. A major advantage of the CSM is its simplicity, since the calculation of the local model parameter only requires the vorticity vector as well as velocity strain rate and in fact, the resulting parameter is wellbounded. Consequently, an averaging procedure, required for the dynamic Smagorinsky model is not needed. It can be shown that both models are in some respect conceptually similar. While the CSM employs the second, the DSM actually utilizes the third invariant of the velocity gradient tensor to compute the model parameter which has a similar distribution in turbulent shear flows according to [40]. Both models implicitly recover the correct nearwall behavior and give zero eddy viscosity in the limit of laminar flows.
Regarding the eddy diffusivity of the energy variable, the dynamic procedure [41] is utilized in all LES cases discussed below.
2 Overset Grid Method Verification
Verification of the overset grid method is performed on a simple twodimensional isentropic vortex advection problem. The analytic solution is given by:(8)(9)(10)(11)
Local coordinates with respect to the vortex core are defined by , , and the radius as . The radius where maximum/minimum velocities occur is denoted R, while β is the vortex strength. Vortex convection is determined by U_{∞}, which is obtained from a prescribed freestream Mach number Ma_{∞}. Equations (8)–(11) are equivalent to the test case used in [42], if .
Verification simulations are computed on a streamwiseperiodic domain of length 15R. In spanwise direction, the analytic solution is prescribed in distance ±5R from the centerline. Two equidistant grid spacings are considered, R/Δx = 25 for the coarser, and R/Δx = 50 for the finer meshes. The freestream Mach number is set to 0.2, while the vortex strength is chosen as β = 1.0. Results are obtained after one cycle, i.e. (t·U_{∞})/R = 15.
2.1 Numerical Accuracy
Consistent numerical accuracy of all methods is important to obtain a smooth spatial error distribution and avoid that the solution is polluted locally. In the following, results from the inviscid vortex problem computed on a coarse mesh, an overall refined mesh, and a coarse background mesh with refinement along the centerline using overset grids, are discussed. The latter case is sketched in Figure 4. After the vortex has traveled through the periodic domain back to its initial location, the flow field is compared to the analytic solution. Using the infinity norm L_{∞} allows to identify even small local errors. In Figure 5, the error norms for the velocity and energy variables are plotted on logarithmic scale. Note that for the coarse resolution, all data is obtained from the same simulation. Lines drawn to the single grid and overset grid solutions, respectively, almost perfectly overlap. Since the errors have not been normalized by the respective flow variable, the offset between the lines should be ignored. The slopes prove second order accuracy of the overset grid method.
Figure 4 Overset grid setup and initial pressure field for accuracy test. 
Figure 5 Infinity error norm after one cycle (not normalized). 
2.2 Vortex Preservation
To assess numerical dissipation and dispersion as flow structures cross an intergrid coupling boundary, the overset mesh is placed as shown in Figure 6. This allows the inviscid vortex to be passed from the coarse to the fine mesh and back during one cycle. The overset grid is again refined by a factor of two. In this grid arrangement, coupling the solutions of all staggered flow variables across domains requires interpolation with coefficients different from unity. After one cycle, the vortex shape should be in good agreement with the analytical solution. This is confirmed by the vorticity contour plotted along the centerline in Figure 7.
Figure 6 Overset grid setup and initial pressure field for vorticity test. 
Figure 7 Vorticity along the centerline after one cycle. 
2.3 Conservation
Although interpolation methods exist which enforce at least primary conservation, formal accuracy is often considered sufficient for flows without discontinuities, if densitybased solvers are used [43]. For the vortex advection test considered in the previous section, the domainintegrated mass, normalized by initial mass is given in Figure 8 as function of time. The maximum error is on the order of (2 × 10^{−6}). With respect to kinetic energy, a slight increase is observed as the vortex passes the refined grid (Fig. 9). Interestingly, the initial value is almost perfectly recovered when the vortex arrives back on the base grid. This may indicate that the observed rise in integrated kinetic energy is caused by directly computing the discrete sum on grids with different resolution, instead of applying a spatial filter before. Still, the maximum error is on the order of (4 × 10^{−5}) which appears reasonably accurate. A modified overset grid location such that grid vertices do not coincide with the background mesh changed this result by up to (2 × 10^{−5}).
Figure 8 Mass in flow domain, normalized by initial value. 
Figure 9 Kinetic energy in flow domain, normalized by initial value. 
2.4 Computational Efficiency
One of the major advantages of applying overset grids to practical applications like engine flows, is the ability to afford high spatial resolution only in regions where the solution accuracy can significantly improve and therefore get computational savings as compared to overall refined grids. This requires an efficient implementation, since additional procedures are needed to couple solutions on composite grids. For the cases used to test accuracy in Section 2.1, measured wall times normalized by the fine singlegrid value are summarized in Figure 10. As expected, the coarsegrid solution is approximately four times cheaper than the one computed on the fine grid. The overset grid case reduces the time to solution by more than 30%, while achieving the same accuracy. By normalizing the measured wall time with the number of computational cells, the overhead needed for intergrid coupling can be estimated. For this case, the composite grid solver only takes about 5% longer per cell, than the single domain solution.
The overset grid method integrated into the inhouse code within this work has been verified and validated for a simple test problem. Computational efficiency of the implementation suggests application to more complex flow problems.
Figure 10 Wall time normalized by case ‘1 Grid/Fine’. 
3 Computation of Internal Combustion Engine Port Flows
3.1 Engine Specification
The engine under investigation of incylinder flow is a singlecylinder 4valve SI engine. Engine specifications are given in Table 1.
Engine specification.
3.2 Experimental Procedure
Stereoscopic Particle Image Velocimetry (SPIV) measurements were conducted in a flowbench test rig for a given valve lift to obtain LES validation data (Fig. 11). The intake airflow is guided through turbulence correction filters in order to smooth the turbulent flow produced by typical fan blades. Large plenums mounted upstream of the port and downstream of the cylinder liner damp out pressure fluctuations. For the coordinate system definition, refer to Figure 11. Measurement planes perpendicular to the cylinder axis are located at z = −B/3, z = −B/2, and z = −2/3B below the fire deck (z = 0). Three velocity components are measured by PIV in these planes. A laser beam is formed to a light sheet and the particles are illuminated by a doublepulsed Ng:YAG laser (New Wave Solo PIV III15) for a short time interval (3–6 μs). The scattered light is captured onto consecutive frames of two high speed digital CCD cameras with a resolution of 1600 × 1200 pixels and a frame rate of 10 Hz. From the movement of particles the velocity vectors are calculated by Direct Cross Correlation Method (DCC) using the software FtrPIV (Flowtech Research Inc.). In this study, an interrogation area size of 33 × 33 pixels with an overlap of 16 × 16 pixels was used. The total interrogation area is 100 × 75 pixels (=1600/16 × 1200/16). This resolves the velocity field within 113.978 × 85.479 mm^{2} with a spatial resolution of 1.14 mm. During a SPIV recording, 1000 doubleframed images were acquired for successive statistical analysis. Such measurement was performed at valve lifts from 1 mm to 10 mm in 1 mm intervals. For each condition, the mass flow rate was adjusted so that the pressure difference between inlet (P1) and atmospheric outlet (P2) of the flow domain is 5.88 kPa. To calculate the tumble ratio from the velocity vectors in a single plane, the AVL approach described in [44] was employed.
Figure 11 Experimental setup for PIV measurements. 
3.3 Simulation Setup
Simulations for 8 mm valve lift are carried out on a flow domain spanning from a location just downstream of the inlet plenum, down to 2.6B below the fire deck, which is 1.25B upstream of the outlet plenum. Different combinations of equidistant grids with either 16 or 32 cells per valve lift are considered. Regions near the inlet and outlet boundaries are always discretized with the coarse mesh, while a refined overset grid is used in the port and upper cylinder region for some cases. Singlegrid setups are entirely discretized with coarse cells (case names ending with ‘C’, cf. Tab. 2). Overset grids overlap by five (coarse) grid cells, without a background mesh being present over the full flow domain. A summary of all simulations performed is given in Table 2.
Nonreflecting boundary conditions are used at both inlet and outlet. Total pressure and total temperature are prescribed at the inlet, while the static pressure according to the measurement is set as target value at the outlet. Regarding the wall model, due to the underlying assumption of fully developed turbulent boundary layer, the model may not be applicable in all parts of the domain. Here, it is applied to the port and valves, while other wall boundaries are handled by the standard IB method. Heat transfer is assumed to be negligible at nearatmospheric conditions.
Two different SFS models are employed. The baseline case in focus for most of the analyses, is computed with the DSML. For comparison, other cases have been run with the CSM.
The time step in all simulations is limited by a constant stability limit Courant–Friedrichs–Lewy Number, CFL = 2.0. Once the flow field has reached steady state, statistics are sampled for 40 ms, which is approximately 3.5 times the duration a fluid particle needs to travel from the inlet to the outlet boundary.
Simulation setup summary.
3.3.1 Parallel Performance
Discretization of complexshaped flow domains with Cartesian grids may lead to a large number of solid cells, which can reduce computational efficiency. Hence, the parallel partitioning algorithm needs to balance the amount of active and inactive cells on each processor. In this work, the efficient block decomposition based on a reduced Cartesian communicator which was previously utilized for singlegrid simulations [7], is applied to each overset mesh separately. While the overset grid method introduces additional interpolation and communication efforts, it can be effectively used to cluster more processors in physical regions of high computational load and improve the overall load balance of the simulation. Representing the engine geometry by compact, concatenated grid patches, may also lead to a smaller number of solid cells as compared to flow domains discretized by a single mesh. For the same number of fluid cells per processor, the fraction of solid cells is reduced from 43.3 for Case DSMLC, to 36.6% for Case DSML. In Figure 12, strong scaling results are given for both simulation setups. For each point on the two different xaxes, the number of fluid cells per partition is the same in the single and the multidomain simulation. It can be concluded that the scaling behavior of both cases is very similar. To compute 40 ms simulation time, Case DSMLC needs less than two days on 400 processors.
Figure 12 Strong scaling for cases DSMLC (upper xaxis) and DSML (lower xaxis). 
3.3.2 Convergence of Statistics
Velocity mean and RMS values are computed after every time step. To demonstrate that the sampling time is sufficiently long to obtain converged statistics, mean and RMS of the vertical velocity component, sampled over different intervals are shown in Figure 13. The data is plotted along a line located at z = −B/2, perpendicular to the tumble rotation yaxis (Fig. 14). Mean velocities mostly collapse to the same curve, except for a region 0 < x < 15 mm, where longer simulation times would still change the result. This location is actually very close to the tumble vortex center, as will be shown below. Away from the vortex core, RMS values are well converged while in the center, the agreement of the curves is still reasonable. Since overall trends in all curves are very similar, a significant benefit from even longer run times cannot be expected.
Figure 13 Mean and RMS of the vertical velocity component for different sampling times (DSML). 
Figure 20 Vertical velocity component in plane at z = −B/2, LES case DSML (left) and PIV data (right). 
3.4 Overall Flow Field Characterization
This section is intended to give a broad idea of the physical flow regimes and bulk flow features. Nondimensional numbers are evaluated from instantaneous realizations to highlight the quality and resolution of the LES results. In Figure 14, the local Reynolds number is plotted on a cut plane through the valve center axis. Note that an arbitrary length scale obtained from the sketched cross section in the horizontal port section has been used. Since the port diameter reduces towards the valves, the actual Reynolds number can be smaller than the color scale may indicate. For 8 mm valve lift and the considered pressure gradient, the typical Reynolds numbers are on the order of 100,000.
The second important dimensionless quantity discussed in the following is the Mach number. For inert ideal gases, it is approximately proportional to the velocity magnitude. Accordingly, we expect maximum Mach numbers in the smallest cross section, which is typically the valve gap. Inside that region, values beyond 0.4 can be observed in Figure 15. Consequently, compressibility effects cannot be neglected. In the colour map, small flow structures on the order of the grid spacing can be seen, which is a feature of LES with accurate numerical schemes.
A common quantity of interest in engine intake port design is the overall pressure loss until the flow enters the cylinder. In general, pressure losses can be classified into linear losses, caused by viscous friction, and singular or local losses. The latter may particularly occur in complex wallbounded flows due to flow separation and formation of recirculation bubbles [45]. To estimate the main reason for losses in the present study, the same cut plane as before is colored by total pressure (Fig. 16). Although there is a small separation zone at the bottom side of the port, right after the first bend, it is obvious that the main decrease in pressure occurs very locally in the direct vicinity of the valve. Wall friction does not seem to reduce the total pressure significantly as the flow approaches the valves. However, there might be a secondary effect, since the velocity profile just upstream of the valve stem might be affected by wall friction and can alter the large singular losses, which occur as the flow passes the valves.
To visualize the incylinder flow field, streamlines computed from the timeaveraged velocity data are depicted in Figure 17. Color contours represent the vertical velocity component, red color means upward motion. For the hightumble engine concept subject to the present analysis, the majority of the inlet mass flow passes over the valves and forms a very strong recirculation, even without being redirected by a piston surface. A complex flow pattern is formed by the interaction of the upward tumble flow and the downward flow which directly enters the cylinder from underneath the valves. Figure 17 also shows one of the measurement planes to establish a link between the local validation data discussed below, and the bulk flow. For the coordinate system definition, the reader is referred to the lower left hand side of the figure.
Figure 14 Instantaneous local Reynolds number in cut plane through valve center axis, based on length scale from indicated cross section (DSMLF). 
Figure 15 Instantaneous Mach number in cut plane through valve center axis (DSMLF). 
Figure 16 Instantaneous total pressure in cut plane through valve center axis (DSMLF). 
Figure 17 Streamlines computed from average velocity field and measurement plane, colored by vertical velocity component (DSML). 
3.5 LES Validation
Engine LES should provide both, reliable prediction of overall flow quantities such as pressure loss or tumble ratio, and accurate incylinder flow (and mixing) fields. In the following, the computational framework is validated in these two respects.
3.5.1 Integral Flow Quantities
In order to assess the agreement between computed and measured pressure loss in the present configuration with prescribed pressures at inlet and outlet of the flow domain, mass flows obtained from the LES and the flowbench test rig can be compared. While for the measurement, three timeaveraged data points are available, Figure 18 shows the mass flow time history at both boundaries in the LES, normalized by the experimental mean value. Note that due to the NSCBC formulation applied (Sect. 1.3), boundary pressures are not strictly enforced, but used as relaxation targets to achieve at least partially reflecting properties. Hence, the mass flow can fluctuate to some extent, which is quantified by a normalized RMS value of 1.3% for the inlet and 3.8% for the outlet boundary, respectively. The mean mass flow is overpredicted by 3.2%. It is worth noting that increasing the order of the nearwall filtering operation from the method described in [12] to onesided, higher order operators as proposed in [13], cf. Section 1.2.1, eliminated a systematic overprediction of the pressure loss from the flowbench simulations. In fact, the slight underprediction observed here can be attributed to several simplifications in the simulation setup. The effect of omitting the outlet plenum will be investigated in the future.
A first quantitative comparison of the incylinder flow can be achieved by analyzing the tumble ratio. It is defined as the ratio of the vortex angular velocity and a characteristic engine angular velocity, which is proportional to the vertical convective flow. All required flow variables to quantify tumble motion are evaluated in the measurement plane located at z = −B/2. In Figure 19, the time histories of tumble ratio calculated from the instantaneous and timemean LES velocity fields (Case DSML) are normalized by the value obtained from ensembleaveraged PIV velocity data. LES variables are continuously averaged over the simulation time, i.e. when the ‘Time Avg.’line approaches a constant value, statistics are converged. At the end of the simulation, the tumble ratio is underpredicted by 2.5%. Since parts of the strongly recirculating flow approach the outlet boundary at the bottom of the cylinder, adding the outlet plenum may affect the tumble results.
Considering the simplifications made for computational efficiency in the present model, the agreement in both integral flow quantities can be considered good. How the wall treatment and mesh refinement affect the results, will be discussed in more detail in Section 3.7.
Figure 18 Instantaneous mass flows at inlet and outlet boundary for case DSML, normalized by experimental value. 
Figure 19 Instantaneous and timeaveraged tumble ratios for case DSML, normalized by experimental value. 
3.5.2 Local Velocity Field
For a more detailed validation of the incylinder flow field, data from three measurement planes are available in the present study. A qualitative comparison of the vertical velocity component (zdirection) presented in Figure 20 suggests reasonable agreement. To make the analysis more quantitative, data interpolated along three distinct lines normal to the tumble vortex rotation axis, located in the symmetry plane and underneath each valve's center axis (Figure 20), is used in this work. The cylinder axis is located at x = 0 mm, Intake Valves (IV) are located at x < 0 mm, the closed exhaust valves (EV) are positioned in the region x > 0 mm. For brevity, only the midmeasurement plane located at z = −B/2 and offcenter lines will be discussed for Case DSML in the following. The central line will be focus of analysis in Section 3.6.1.
We begin the quantitative discussion with the vertical velocity component, which is by far largest in magnitude. Since the engine geometry is symmetric, mean data from the two offcenter lines are plotted in Figure 21 for both LES and PIV. Instantaneous velocities from 625 realizations of the computed flow field are depicted in grey color. With respect to the experimental data (red symbols), it is obvious that mean velocities at the two locations differ by up to 30% in the highvelocity region, which can be due to flow instabilities, experimental error or other reasons. A RANS simulation of this configuration actually predicted symmetric velocity fields. Regarding the general velocity distribution, very strong downward flow is present near the liner wall at the exhaust side. The tumble recirculation caused upward flow on the left hand side of the cylinder axis, where the intake valves are located. Overall, the LES mean velocity results are in good agreement with the experimental data, but the discrepancy between the two LES lines occurs at a different location than in the measurement. Interestingly, the deviations are of similar magnitudes in both datasets. Longer averaging time is not expected to remove the asymmetry from the LES data, as discussed in Section 3.3.2. Close to the liner wall on the exhaust side, the computed velocity magnitude rapidly decreases, but there are no measurement points to validate the nearwall behavior of the numerical solution. Since there is very good agreement at the first experimental data point off the wall, the standard IB interpolation applied at the cylinder walls is considered to be a good approximation.
In Figure 22, the mean velocity in tumble rotation (y) axis direction is depicted. From the experimental results it can be concluded that on the intake side (x < 0 mm) liner wall, i.e. where the mean flow recirculates upward, transverse velocities are near zero which indicates a wellformed tumble vortex. Towards the cylinder axis, the yvelocity component increases beyond the magnitude of the vertical velocity and forms a net flow towards the center plane. The largest magnitudes are observed near the exhaustside (x > 0 mm) liner wall, where the strong downward flow is directed towards the symmetry plane by the curvature of the liner walls. As the computed instantaneous data indicate, single flow realizations may differ significantly from the described trends in the mean flow. The timeaveraged LES results reproduce the measured values well, except for the liner wall on the right, where velocity magnitudes are underpredicted.
The horizontal tumble plane velocity (xdirection) is overall smallest in magnitude (Figure 23). In the recirculation region (x < 0 mm), there is almost no flow in this direction. On the exhaust side, the xvelocity is mostly negative, indicating flow towards the intake side. Only in the direct vicinity of the liner, flow still approaches the wall. The computed magnitudes in the mean match the experiment quite well. However, a more intense mean flow towards the intake side is observed near the cylinder axis, which could be due to a too high vertical tumble vortex position or mismatch in bulk vortex shape.
To assess the ability of the numerical framework to predict velocity fluctuations, RMS values will be discussed in the following. In Figure 24, fluctuations in the vertical velocity component are shown along the two offcenter lines. Maximum values occur on the exhaust side in the region of highmomentum flow, as expected. However, in this data set, the observed asymmetry is actually most severe. The measurement shows high fluctuations near the intake side liner wall for one of the two lines as well. In this region, the downward flow from the intake valve along the (left) liner wall interacts with the tumble vortex and may form a very unstable stagnation flow which might cause the large RMS values. In the LES data, fluctuations of very similar magnitudes are found, but the overall distribution seems to be shifted to the left. On the intake side, lower fluctuations compared to the PIV data are predicted for one of the two analysis lines.
In direction of the tumble rotation axis (y), maximum RMS values are of similar order and occur in the same region as the fluctuations in vertical direction, cf. Figure 25. If this agreement is due to a motion of the highvelocity jet tangential to the wall or, due to the presence of highly intermittent turbulent flow structures, or both, is difficult to conclude. Agreement between simulation and experiment is overall very good. Near the cylinder axis, the LES overpredicts the measured fluctuations in one of the two analysis locations.
Fluctuations in the xvelocity component (Fig. 26) are most intense on the exhaust side near the center axis, where mean flow towards the intake side exists. LES results are found to be in good agreement with the measurement, except for some overprediction of the RMS on the intake side in one of the two offcenter lines.
Quantitative analysis of the incylinder velocity field has shown good agreement between LES and PIV results. However, the match is not equally well in all locations which might be due to slightly different tumble vortex locations. Since the flow field is very complex and the free recirculating flow forms several stagnation regions which are unstable by nature, validation of the numerical framework can be considered successful, given the overall agreement in both mean and RMS velocities.
Figure 21 Instantaneous (LES) and mean zvelocity at z = −B/2, y = ±20 mm. 
Figure 22 Instantaneous (LES) and mean yvelocity at z = −B/2, y = ±20 mm. 
Figure 23 Instantaneous (LES) and mean xvelocity at z = −B/2, y = ±20 mm. 
Figure 24 zvelocity RMS at z = −B/2, y = ±20 mm. 
Figure 25 yvelocity RMS at z = −B/2, y = ±20 mm. 
Figure 26 xvelocity RMS at z = −B/2, y = ±20 mm. 
3.6 Robustness
A required property of LES for practical applications is the reproducibility of results for different grids and model combinations. This will be investigated in the following.
3.6.1 SFS Model and Grid Resolution
The same simulation performed for validation purposes and discussed in Section 3.5.2 has been repeated on different grids with all parameters being the same, except for the subfilter eddy viscosity model. In fact, results from the baseline case DSML computed using the dynamic Smagorinsky model with Lagrangian averaging along particle trajectories [35], are compared to simulation data obtained using the coherent structure model [37], Case CSM. For a summary of the two concepts, the reader is referred to Section 1.4.2. Since a priori, none of the two models can be expected to perform better, sensitivity of the LES solution will be discussed in the following. The overall effect of an SFS model on the LES solution is griddependent. Hence, it seems appropriate to perform the comparison on different grids as well. Note that the purpose of the parameter variation is not to identify optimal settings, but to understand sensitivities for the given simulation setup.
As before, the analysis starts with integral quantities. In Figure 27, averaged mass flows from different LES setups, normalized by the experimental value, are plotted for direct comparison. Since the overall pressure loss is governed by singular losses in the valve region, and less dependent on wall friction, both different grid resolution and SFS stresses may have an effect. To quantify the former, Cases DSML and DSMLC, as well as CSM and CSMC can show the effect of refining the port and upper cylinder region by a factor of two (Table 2). The observed difference accounts for 2–3%, which seems reasonable. Comparatively small is the effect of the SFS model on the same grid, while higher resolution seems to increase the sensitivity.
The mass flow distribution over the valve curtain surface as well as the interaction of the unsteady highReynolds number flow with the valves affect the intake jet and consequently the tumble vortex. Accordingly, grid resolution and SFS models can alter the tumble ratio. Since the tumble ratio takes into account the convective flow, it should be independent of the actual mass flow. This is confirmed by comparing Cases DSMLC and CSMC in Figures 27 and 28, respectively. While the mass flow is not affected by the SGS model, the tumble ratio is changed by 3.8%. On the refined grid, the SFS model affects the tumble vortex slightly less. Increased tumble ratios by approximately 5% are observed on the overall coarse as compared to the refined grid.
To show that placing a fine/coarse intergrid coupling boundary in the tumble region does not significantly change the results, the boundary was moved from z = −0.3B for the baseline case (DSML) to z = −B for Case DSMLF, cf. Table 2. From Figures 27 and 28 it can be concluded that the effect is indeed minor.
For the quantitative comparison of velocity data, the discussion is limited to the variation of SFS models on the refined grid, since no systematic difference between solutions computed with and without grid refinement was found. Different to the discussion in Section 3.5.2, velocity data interpolated to a line perpendicular to the tumble axis, but this time located in the symmetry plane (Fig. 14), is used. Considering the mean vertical velocity component (z) plotted in Figure 29, the experimental data are very similar to the offcenter analysis lines discussed before. Also in this location, the LES results from the baseline Case DSML are in very good agreement with the measurement. The CSM simulation overall predicts the mean zvelocity similarly well, except for a region on the exhaust side, close to the cylinder axis where the change in sign is shifted to the left. This is due to a slight offset in tumble vortex location towards the intake side, as will be shown later. Note that in this region, results may slightly change if the averaging time is increased (Sect. 3.3.2).
For symmetry reasons, the mean yvelocity is expected to be close to zero in the stagnation plane, which is confirmed by Figure 30. Both simulations reproduce the measured dataset reasonably well. Locally, the CSM model shows larger discrepancy which again is most pronounced near x = 8 mm.
Regarding the tumble plane horizontal (x) velocity component, peak magnitudes are slightly higher than in the offcenter analysis line locations, which can be explained by flow stagnation in ydirection. While the DSML predicts the measured velocity distribution quite well, the CSM again gives larger mismatch around x = 8 mm and near the exhaustside liner wall (Fig. 31).
Velocity fluctuations for the z, y and xvelocities are plotted in Figures 32–34, respectively. Overall, maximum velocity RMS values along the center line are smaller than in the offcenter locations which is likely due to the absence of the primary jets directly penetrating into the cylinder from the valves. Instead in the symmetry plane, the conical jets form a stagnation flow which penetrates downwards. Similar maximum RMS values are found in all three velocity components and all in the same region close to the exhaust side liner wall. The general agreement of numerical and experimental data is good and mostly similar for both SFS models. Again, the central region close to x = 8 mm has consistently higher velocity fluctuations in the case computed with CSM. However, it should be noted that even longer sampling times for the statistics may change the RMS values in x, y and zdirection in the region of 0 < x < 15 mm by up to 0.5, 1.7 and 1.3 m/s, respectively.
To estimate tumble vortex position, the Γ_{1} vortex identification criterion [46] applicable to twodimensional flows has been computed in three tumble cut planes located at the same ylocations as the line data discussed previously, cf. Figure 35. Note that the analysis is strictly valid if the planenormal velocity component is zero or at least constant, which can be assumed for the symmetry plane. For the comparison, the characteristic vortex length has been set to B/4 and the resulting Γ_{1} fields have been rescaled such that the global maximum is equal to unity in both cases. The tumble vortex core can then be identified by Γ_{1} = 1. In Figure 35, the isocontours for Γ_{1} = 0.8, computed from the mean velocity field, are shown. For orientation purposes, the line at z = −B/2 as well as the cylinder axis and the x = 8 mm location are highlighted in the center plane. Obviously, the vortex cores computed with different SFS models differ by approximately 5 mm, which explains the discrepancies in the velocity field mentioned above. Interestingly, the predicted vortex cores at y = −20 mm (background) are in good agreement while there is an apparent offset at y = 20 mm (foreground).
To give an impression of the differences in the results of the two LES in the region upstream of the tumble vortex, a location with similar discrepancies between the two SFS models as at (x = 8 mm, y = 0, z = −B/2) was identified in the z = −B/3 plane. Upstream of that location, a streamline was computed for the DSML data set. The velocity field computed with the CSM has been interpolated to that streamline. In Figure 36, mean velocity magnitude and RMS of the zvelocity are plotted along the streamline length. On the righthandside axis, the relative difference between the curves is quantified. Regarding the mean velocity, the two solutions are very similar along the intake port. As the flow passes the valve, the curves differ locally by up to 5%. Major discrepancies seem to exist inside the cylinder, but are very likely due to different vortex locations. Accordingly, the maximum relative offset of more than 30% should not be given too much importance, since the flow fields can still be qualitatively very similar.
Similar conclusions can be drawn from the fluctuations in vertical velocity (Figure 36). Although the relative difference upstream of the sharp bend in the port is high, it should be noted that the RMS values are close to zero in that region. More significant variations begin to exist locally near the valve and in the cylinder. However, small offsets in bulk flow structures may cause these differences. Therefore, the performance of both models throughout the flow domain can be considered very similar.
Figure 27 Mass flow computed with different grids and SFS models, normalized by experimental value. 
Figure 28 Tumble ratio computed with different grids and SFS models, normalized by experimental value. 
Figure 29 Mean zvelocity at z = −B/2, y = 0 mm. 
Figure 30 Mean yvelocity at z = −B/2, y = 0 mm. 
Figure 31 Mean xvelocity at z = −B/2, y = 0 mm. 
Figure 32 zvelocity RMS at z = −B/2, y = 0 mm. 
Figure 33 yvelocity RMS at z = −B/2, y = 0 mm. 
Figure 34 xvelocity RMS at z = −B/2, y = 0 mm. 
Figure 35 Isocontour Γ_{1} = 0.8 in symmetry and valve axis center planes. 
Figure 36 Mean velocity magnitude and zRMS along streamline. 
3.7 Sensitivities
Although the quantitative influence of various numerical methods and physical models on the LES result might be casespecific, some trends should be applicable to engine port flows in general. For brevity, only integral flow quantities are considered below.
3.7.1 Effect of Wall Model
The wall model provides the velocity boundary condition to the IB method in the port region and at the valves, as well as the shear stress through a model viscosity. In Figure 37, results from the baseline case (DSML) are compared to a simulation without wall model, but standard IB interpolation (DSMLNW). As mentioned before, the overall loss in total pressure is not governed by wall friction. In fact, the increased wall friction provided by the model is overcompensated by some secondary effects, likely due to the modified velocity profile upstream of the valves. However, the effect accounts for less than 5% increase in mass flow. More profound is the effect of the wall model on tumble ratio. The mass flow distribution across the valves is changed such that the tumble ratio increases by 25%. Therefore, we consider the wall model an indispensable component of a numerical framework for simulation of SI engine port flows.
Figure 37 Mass flow and tumble ratio computed with and without wall model, normalized by experimental values. 
3.7.2 Combined Effects
To quantify the effects of all methods and models implemented into the numerical framework as part of this study, results computed without enhanced nearwall filtering (Sect. 1.2), IB method (Sect. 1.2.2), wall model (Sect. 1.4.1) and overset grids (Sect. 1.2.3), Case DSMLPRVC, are compared to the baseline case in Figure 38. It is obvious that previously, both integral flow measures were significantly underpredicted. With the extended framework, mass flow is increased by 12% and tumble ratio by 31%, respectively. Regarding the overset grid method, for the mesh spacings considered in this study, refinement does not significantly improve the results. However, it can be employed in the future to locally coarsen the mesh for computational efficiency. Systematic errors in the wall treatment have been eliminated, but careful reassessment of remaining discrepancies is necessary to identify shortcomings in the simulation setup.
Figure 38 Mass flow and tumble ratio computed with different methods and models, normalized by experimental values. 
3.8 Valve Jet Dynamics
One of the major advantages of LES is the ability to capture local unsteady effects. Since fluctuations in the conical valve jet are considered as one possible cause of CCV in SI engines [47], it seems meaningful to investigate inflow dynamics in a flowbench configuration. For a qualitative identification of most intermittent regions, the magnitude of the RMS velocity vector can be considered. In Figure 39, streamlines computed from the mean velocity field are colored by the RMS magnitude. The 30 m/s isosurface is sketched to mark the highly fluctuating zone. As the flow penetrates downward into the cylinder, the RMS values decrease. In the regions where a jet interacts with the liner walls or with the other jet, fluctuations seem to exist until further downstream. Note that the RMS does not allow to distinguish between turbulent velocity fluctuations and intermittent motion of bulk flow structures.
For a more quantitative analysis, decomposing the velocity field into spatial and temporal modes can be very useful. Proper Orthogonal Decomposition (POD) [48] is employed on a subsection of the domain (Fig. 40). For Case DSMLF, a total number of 625 realizations have been sampled over 40 ms simulation time and are used for the analysis. In Figure 41, the kinetic energy fraction of the first 20 modes, sorted by energy content, is depicted. Since the first mode corresponds to the mean flow, the kinetic energy is dominant compared to the dynamic modes. Regarding the latter, the energy content decreases very gradually and it is difficult to make a clear distinction into more and less important modes, just by the relative contribution to total kinetic energy. In the following, modes two to six are arbitrarily selected for visualization.
In Figure 40, the basis function of the first mode is represented by streamlines. For the higher modes, only vectors with more than 50% of the maximum magnitude in each mode are plotted as arrows. The color scale indicates the angle between the basis function vector of mode one, i.e. the average flow direction, and the respective dynamic mode direction. Note that the vector length has no meaning. It seems that most of the arrows are clustered near the symmetry plane and are mostly oriented in the transverse direction. This may indicate that the flow dynamics are governed by the interaction of the two valve jets, similar as previously observed in a different engine geometry [49]. The regions of high velocity RMS near the liner walls in Figure 39, are not represented by any of the first six modes. However, some of the selected modes do contribute to flow dynamics downstream of the valve stem. Although directions of these modes deviate less from the bulk flow direction than in the symmetry plane, it is difficult to draw a strong conclusion on the more critical region for incylinder flows in this geometry. Still, kinetic energies contained in the first modes are very similar. Typical frequencies of the visualized modes are on the order of 3000 Hz. Assuming sufficiently large valve lifts for a duration of 125 degrees crank angle in an engine operated at 3000 rpm, approximately 20 periods of the considered POD modes would occur. Under the assumption that similar flow dynamics develop in the running engine, these modes might be relevant for CCV.
To eliminate the effect of the flow instability in the symmetry plane, another LES has been performed. Instead of using a mesh for the full engine geometry, only half of the flow domain is considered. At the central plane, symmetry boundary conditions are applied, imitating two perfectly identical port flows and valve jets. Velocity fields have been sampled with the same frequency as before. In Figure 42, the upper 50% of kinetic energy in each of the modes two to six are again visualized by vector arrows, this time computed from the halfgeometry results. The second mode still has some transverse components, but does not only act in a small region, as in the full geometry. Although a larger number of modes is clustered in the valve wake, the fluctuations still seem to mostly occur with smaller deviation from the mean streamline direction.
To quantify this observation, the distributions of the angles between basis function vectors and the timeaveraged streamlines from Figures 40 and 42 are plotted in Figure 43. While results computed in the full engine geometry have a bimodal distribution with peaks at 30 and 90°, two distinct maxima can be identified at 17 and 38° in the halfgeometry LES data. Hence, dynamic modes which act perpendicular to the mean flow direction seem to mostly exist in the symmetry plane and occur much less frequently, if half of the flow domain is omitted in the simulation. Which of the two effects, the central plane transverse flow intermittency, or pulsation of the bulk tumble flow are more important for CCV, needs to be investigated in a similar LES study of a running engine.
Figure 39 Streamlines colored by velocity RMS vector magnitude. 
Figure 40 Streamlines computed from basis function of mode 1, and vectors with magnitude of at least 50% of the basisfunctionmaximum for POD modes 2–6, colored by angular deviation from mean streamline (LES on full geometry). 
Figure 41 Energy fraction of POD modes. 
Figure 42 Streamlines computed from basis function of mode 1, and vectors with magnitude of at least 50% of the basisfunctionmaximum for POD modes 2–6, colored by angular deviation from mean streamline (LES on half geometry). 
Figure 43 Distribution of angle between timeaveraged streamlines and basis function vectors of POD modes 2–6. 
Conclusion
Challenges intrinsic to LES of wallbounded flows in complex geometries on Cartesian grids have been addressed by the implementation of three numerical methods and a wall model into an inhouse code. An overset grid approach has been demonstrated to be accurate as well as efficient. The enhanced numerical framework was successfully applied to simulate an engine flowbench configuration. For prediction of integral flow quantities like overall pressure loss and tumble ratio, the enhanced wall treatment by a combination of an IB method, higherorder filtering and a wall function, has been shown to be crucial. Local mesh refinement can be used to reduce computational cost. Parameter variations with three different grids and two SFS models have proven robustness of the simulation setup. All five simulations reproduced measured integral flow numbers within a deviation of 5%. The computed velocity field was validated against PIV measurement data. Good agreement in mean and RMS velocities was observed. Direct comparison between two SFS models, the DSML and CSM has shown very similar performance. Due to its simplicity, the CSM seems an attractive alternative for engine flows. Valve jet dynamics have been analyzed by means of POD. Highly unsteady flow motion has been identified in transverse direction in the symmetry plane, and more towards streamwise direction in the valve wake. If similar effects occur during engine operation, needs to be investigated in further studies. The enhanced numerical framework is currently being extended to moving geometries and will be used to compute motored and fired engine operation in the near future.
Acknowledgements
The authors from RWTH Aachen University gratefully acknowledge partial funding by Honda R&D and by the Cluster of Excellence “TailorMade Fuels from Biomass”, which is funded by the Excellence Initiative of the German federal and state governments to promote science and research at German universities.
The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gausscentre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de).
S.K. gratefully acknowledges financial support from a National Research Foundation of Korea (NRF) grant, funded by the Korean government (MSIP) (No. 2012M2A8A4055647).
References
 Fansler T.D., Wagner R.M. (2015) Cyclic dispersion in engine combustionintroduction by the special issue editors, Int. J. Engine Res. 16, 255–259 [CrossRef] [Google Scholar]
 Enaux B., Granet V., Vermorel O., Lacour C., Pera C., Angelberger C., Poinsot T. (2011) LES study of cycletocycle variations in a spark ignition engine, Proc. Combust. Inst. 33, 2, 3115–3122 [CrossRef] [Google Scholar]
 Moureau V., Barton I., Angelberger C., Poinsot T. (2004) Towards large eddy simulation in internalcombustion engines: simulation of a compressed tumble flow, in: SAE Technical Paper, SAE International, 06 [Google Scholar]
 Truffin K., Angelberger C., Richard S., Pera C. (2015) Using largeeddy simulation and multivariate analysis to understand the sources of combustion cyclic variability in a sparkignition engine, Combust. Flame 162, 12, 4371–4390 [CrossRef] [Google Scholar]
 Janas P., Wlokas I., Bohm B., Kempf A. (2017) On the evolution of the flow field in a spark ignition engine, Flow Turbul. Combust. 98, 1, 237–264 [CrossRef] [Google Scholar]
 Nguyen T.M., Proch F., Wlokas I., Kempf A.M. (2016) Large eddy simulation of an internal combustion engine using an efficient immersed boundary technique, Flow Turbul. Combust. 97, 1, 191–230 [CrossRef] [Google Scholar]
 Mittal V., Kang S., Doran E., Cook D., Pitsch H. (2014) LES of gas exchange in IC engines, Oil Gas Sci. Technol.: Rev. d'lFP Energies nouvelles 69, 1, 29–40 [CrossRef] [EDP Sciences] [Google Scholar]
 Verzicco R., MohdYusof J., Orlandi P., Haworth D. (2000) Large eddy simulation in complex geometric configurations using boundary body forces, AIAA J. 38, 3, 427–433 [CrossRef] [Google Scholar]
 Mittal R., Iaccarino G. (2005) Immersed boundary methods, Ann. Rev. Fluid Mech. 37, 1, 239–261 [NASA ADS] [CrossRef] [Google Scholar]
 Meakin R.L. (1998) Composite overset structured grids, in: Handbook of Grid Generation, CRC Press, December 1998 [Google Scholar]
 Berger M.J., Oliger J. (1984) Adaptive mesh refinement for hyperbolic partial differential equations, J. Comput. Phys. 53, 3, 484–512 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Gaitonde D., Shang J., Young J. (1997) Practical aspects of highorder accurate finitevolume schemes for electromagnetics, in: AIAA Paper 970363 [Google Scholar]
 Gaitonde D.V., Visbal M.R. (2000) Padetype higherorder boundary filters for the NavierStokes equations, AIAA J. 38, 2103–2112 [CrossRef] [Google Scholar]
 Kang S., Iaccarino G., Ham F., Moin P. (2009) Prediction of wallpressure fluctuation in turbulent flows with an immersed boundary method, J. Comput. Phys. 228, 18, 6753–6772 [CrossRef] [Google Scholar]
 Desjardins O., Blanquart G., Balarac G., Pitsch H. (2008) High order conservative finite difference scheme for variable density low mach number turbulent flows, J. Comput. Phys. 227, 15, 7125–7159 [CrossRef] [MathSciNet] [Google Scholar]
 Liu X.D., Osher S., Chan T. (1994) Weighted essentially nonoscillatory schemes, J. Comput. Phys. 115, 1, 200–212 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Trisjono P., Kang S., Pitsch H. (2016) On a consistent highorder finite difference scheme with kinetic energy conservation for simulating turbulent reacting flows, J. Comput. Phys. 327, 612–628 [CrossRef] [Google Scholar]
 Poinsot T., Garcia M., Senoner J.M., Gicquel L., Staffelbach G., Vermorel O. (2011) Numerical and physical instabilities in massively parallel les of reacting flows, J. Sci. Comput. 49, 1, 78–93 [CrossRef] [Google Scholar]
 Yee H.C., Sjogreen B. (2002) Designing adaptive lowdissipative high order schemes for longtime integrations, Springer, The Netherlands, Dordrecht, pp. 141–198 [Google Scholar]
 Hu F., Hussaini M., Manthey J. (1996) Lowdissipation and lowdispersion RungeKutta schemes for computational acoustics, J. Comput. Phys. 124, 1, 177–191 [NASA ADS] [CrossRef] [Google Scholar]
 Stanescu D., Habashi W. (1998) 2nstorage low dissipation and dispersion RungeKutta schemes for computational acoustics, J. Comput. Phys. 143, 2, 674–681 [CrossRef] [Google Scholar]
 Waheed A., Yan J. (1998) Workload characterization of CFD applications using partial differential equation solvers, Technical report, Nasa Ames Research Center, Technical Report NAS98011 [Google Scholar]
 Fadlun E., Verzicco R., Orlandi P., MohdYusof J. (2000) Combined immersedboundary finitedifference methods for threedimensional complex flow simulations, J. Comput. Phys. 161, 1, 35–60 [NASA ADS] [CrossRef] [Google Scholar]
 Kim J., Kim D., Choi H. (2001) An immersedboundary finitevolume method for simulations of flow in complex geometries, J. Comput. Phys. 171, 1, 132–150 [CrossRef] [Google Scholar]
 Fedkiw R.P., Marquina A., Merriman B. (1999) An isobaric fix for the overheating problem in multimaterial compressible flows, J. Comput. Phys. 148, 2, 545–578 [CrossRef] [Google Scholar]
 Meakin R. (2001) Object Xrays for cutting holes in composite overset structured grids, in: Fluid Dynamics and Colocated Conferences, American Institute of Aeronautics and Astronautics, June 2001 [Google Scholar]
 Noack, R. (2016) On overset hole cutting and the problems encountered, in: 13th Symposium on Overset Composite Grids And Solution Technology, Mukilteo, WA, USA [Google Scholar]
 Bonet J., Peraire J. (1991) An alternating digital tree (ADT) algorithm for 3D geometric searching and intersection problems, Int. J. Numer. Methods Eng. 31, 1, 1–17 [NASA ADS] [CrossRef] [Google Scholar]
 Chesshire, G., Henshaw W. (1990) Composite overlapping meshes for the solution of partial differential equations, J. Comput. Phys. 90, 1, 1–64 [CrossRef] [MathSciNet] [Google Scholar]
 Lodato G., Domingo P., Vervisch L. (2008) Threedimensional boundary conditions for direct and largeeddy simulation of compressible viscous flows, J. Comput. Phys. 227, 10, 5105–5143 [CrossRef] [Google Scholar]
 Roman F., Armenio V., Frohlich J. (2009) A simple walllayer model for large eddy simulation with immersed boundary method, Phys. of Fluids 21, 10, 101701 [CrossRef] [Google Scholar]
 Lee J., Cho M., Choi H. (2013) Large eddy simulations of turbulent channel and boundary layer flows at high Reynolds number with mean wall shear stress boundary condition, Phys. Fluids 25, 11, 110808 [CrossRef] [Google Scholar]
 Yang X.I.A., Sadique J., Mittal R., Meneveau C. (2015) Integral wall model for large eddy simulations of wallbounded turbulent flows, Phys. Fluids 27, 2, 025112 [CrossRef] [Google Scholar]
 Hoyas S., Jimenez J. (2006) Scaling of the velocity fluctuations in turbulent channels up to Re_{T} = 2003, Phys. Fluids 18, 1, 011702 [CrossRef] [Google Scholar]
 Meneveau C., Lund T.S., Cabot W.H. (1996) A Lagrangian dynamic subgridscale model of turbulence, J. Fluid Mech. 319, 353–385 [CrossRef] [Google Scholar]
 Silvis M.H., Remmerswaal R.A., Verstappen R. (2017) Physical consistency of subgridscale models for largeeddy simulation of incompressible turbulent flows, Phys. Fluids 29, 1, 015105 [CrossRef] [Google Scholar]
 Kobayashi H. (2005) The subgridscale models based on coherent structures for rotating homogeneous turbulence and turbulent channel flow, Phys. Fluids 17, 4, 045104 [CrossRef] [Google Scholar]
 Smagorinsky J. (1963) General circulation experiments with the primitive equations, Monthly Weather Rev. 91, 3, 99–164 [NASA ADS] [CrossRef] [Google Scholar]
 Kobayashi H. (2006) Large eddy simulation of magnetohydrodynamic turbulent channel flows with local subgridscale model based on coherent structures, Phys. Fluids 18, 4, 045107 [CrossRef] [Google Scholar]
 Tanahashi M., Miyauchi T. (1995) Small scale eddies in turbulent mixing layer, in: Proceedings of the Tenth Symposium on Turbulent Shear Flows, Vol. 1, pp. 79–84 [Google Scholar]
 Balarac G., Pitsch H., Raman V. (2008) Development of a dynamic model for the subfilter scalar variance using the concept of optimal estimators, Phys. Fluids 20, 035114 [CrossRef] [Google Scholar]
 Yee H., Sandham N., Djomehri M. (1999) Lowdissipative highorder shockcapturing methods using characteristicbased filters, J. Comput. Phys. 150, 1, 199–238 [CrossRef] [MathSciNet] [Google Scholar]
 Shih T. (2002) Overset grids: Fundamental and practical issues, in Fluid Dynamics and Colocated Conferences, American Institute of Aeronautics and Astronautics, June 2002 [Google Scholar]
 Laimboeck F.J., Glanz R., Modre E., Rothbauer R.J. (1999) AVL approach for small 4stroke cylinderhead, port and combustion chamber layout, in: SAE Technical Paper, SAE International, 09 [Google Scholar]
 Hager W.H. (2010) Wastewater Hydraulics, SpringerVerlag, Berlin, Heidelberg [CrossRef] [Google Scholar]
 Graftieaux L., Michard M., Grosjean N. (2001) Combining PIV, POD and vortex identification algorithms for the study of unsteady turbulent swirling flows, Meas. Sci. Technol. 12, 9, 1422 [Google Scholar]
 J. Heywood, Internal Combustion Engine Fundamentals. Automotive technology series, McGrawHill, 1988 [Google Scholar]
 Lumley J.L. (1967) The structure of inhomogeneous turbulent flows, in: Yaglom A.M., Tatarski V.I. (eds.), Atmospheric turbulence and radio propagation, Nauka, Moscow, pp. 166–178 [Google Scholar]
 Falkenstein T., Bode M., Kang S., Pitsch H., Arima T., Taniguchi, H. (2015) Largeeddy simulation study on unsteady effects in a statistically stationary SI engine port flow, in: SAE Technical Paper, SAE International, 04 [Google Scholar]
All Tables
All Figures
Figure 1 Approximated flow domain with immersed boundary. 

In the text 
Figure 2 Overset grid for boundary layer refinement in an intake port. 

In the text 
Figure 3 Mean streamwise velocity (top) and streamwise velocity fluctuation (bottom) as function of wall distance. 

In the text 
Figure 4 Overset grid setup and initial pressure field for accuracy test. 

In the text 
Figure 5 Infinity error norm after one cycle (not normalized). 

In the text 
Figure 6 Overset grid setup and initial pressure field for vorticity test. 

In the text 
Figure 7 Vorticity along the centerline after one cycle. 

In the text 
Figure 8 Mass in flow domain, normalized by initial value. 

In the text 
Figure 9 Kinetic energy in flow domain, normalized by initial value. 

In the text 
Figure 10 Wall time normalized by case ‘1 Grid/Fine’. 

In the text 
Figure 11 Experimental setup for PIV measurements. 

In the text 
Figure 12 Strong scaling for cases DSMLC (upper xaxis) and DSML (lower xaxis). 

In the text 
Figure 13 Mean and RMS of the vertical velocity component for different sampling times (DSML). 

In the text 
Figure 20 Vertical velocity component in plane at z = −B/2, LES case DSML (left) and PIV data (right). 

In the text 
Figure 14 Instantaneous local Reynolds number in cut plane through valve center axis, based on length scale from indicated cross section (DSMLF). 

In the text 
Figure 15 Instantaneous Mach number in cut plane through valve center axis (DSMLF). 

In the text 
Figure 16 Instantaneous total pressure in cut plane through valve center axis (DSMLF). 

In the text 
Figure 17 Streamlines computed from average velocity field and measurement plane, colored by vertical velocity component (DSML). 

In the text 
Figure 18 Instantaneous mass flows at inlet and outlet boundary for case DSML, normalized by experimental value. 

In the text 
Figure 19 Instantaneous and timeaveraged tumble ratios for case DSML, normalized by experimental value. 

In the text 
Figure 21 Instantaneous (LES) and mean zvelocity at z = −B/2, y = ±20 mm. 

In the text 
Figure 22 Instantaneous (LES) and mean yvelocity at z = −B/2, y = ±20 mm. 

In the text 
Figure 23 Instantaneous (LES) and mean xvelocity at z = −B/2, y = ±20 mm. 

In the text 
Figure 24 zvelocity RMS at z = −B/2, y = ±20 mm. 

In the text 
Figure 25 yvelocity RMS at z = −B/2, y = ±20 mm. 

In the text 
Figure 26 xvelocity RMS at z = −B/2, y = ±20 mm. 

In the text 
Figure 27 Mass flow computed with different grids and SFS models, normalized by experimental value. 

In the text 
Figure 28 Tumble ratio computed with different grids and SFS models, normalized by experimental value. 

In the text 
Figure 29 Mean zvelocity at z = −B/2, y = 0 mm. 

In the text 
Figure 30 Mean yvelocity at z = −B/2, y = 0 mm. 

In the text 
Figure 31 Mean xvelocity at z = −B/2, y = 0 mm. 

In the text 
Figure 32 zvelocity RMS at z = −B/2, y = 0 mm. 

In the text 
Figure 33 yvelocity RMS at z = −B/2, y = 0 mm. 

In the text 
Figure 34 xvelocity RMS at z = −B/2, y = 0 mm. 

In the text 
Figure 35 Isocontour Γ_{1} = 0.8 in symmetry and valve axis center planes. 

In the text 
Figure 36 Mean velocity magnitude and zRMS along streamline. 

In the text 
Figure 37 Mass flow and tumble ratio computed with and without wall model, normalized by experimental values. 

In the text 
Figure 38 Mass flow and tumble ratio computed with different methods and models, normalized by experimental values. 

In the text 
Figure 39 Streamlines colored by velocity RMS vector magnitude. 

In the text 
Figure 40 Streamlines computed from basis function of mode 1, and vectors with magnitude of at least 50% of the basisfunctionmaximum for POD modes 2–6, colored by angular deviation from mean streamline (LES on full geometry). 

In the text 
Figure 41 Energy fraction of POD modes. 

In the text 
Figure 42 Streamlines computed from basis function of mode 1, and vectors with magnitude of at least 50% of the basisfunctionmaximum for POD modes 2–6, colored by angular deviation from mean streamline (LES on half geometry). 

In the text 
Figure 43 Distribution of angle between timeaveraged streamlines and basis function vectors of POD modes 2–6. 

In the text 