Identification of Large-Scale Structure Fluctuations in IC Engines using POD-Based Conditional Averaging

— Cycle-to-Cycle Variations (CCV) in IC engines is a well-known phenomenon and the de ﬁ nition and quanti ﬁ cation is well-established for global quantities such as the mean pressure. On the other hand, the de ﬁ nition of CCV for local quantities, e.g. the velocity or the mixture distribution, is less straightforward. This paper proposes a new method to identify and calculate cyclic variations of the ﬂ ow ﬁ eld in IC engines emphasizing the different contributions from large-scale energetic (coherent) structures, identi ﬁ ed by a combination of Proper Orthogonal Decomposition (POD) and conditional averaging, and small-scale ﬂ uctuations. Suitable subsets required for the conditional averaging are derived from combinations of the the POD coef ﬁ cients of the second and third mode. Within each subset, the velocity is averaged and these averages are compared to the ensemble-averaged velocity ﬁ eld, which is based on all cycles. The resulting difference of the subset-average and the global-average is identi ﬁ ed as a cyclic ﬂ uctuation of the coherent structures. Then, within each subset, remaining ﬂ uctuations are obtained from the difference between the instantaneous ﬁ elds and the corresponding subset average. The proposed methodology is tested for two data sets obtained from scale resolving engine simulations. For the ﬁ rst test case, the numerical database consists of 208 independent samples of a simpli ﬁ ed engine geometry. For the second case, 120 cycles for the well-established


INTRODUCTION
The development of Internal Combustion (IC) engines is characterized by a large number of conflicting targets such as reducing the fuel consumption and the pollutant emissions while at the same time, the engine performance and the combustion stability has to be increased.Cycle-to-Cycle Variations (CCV) have a strong impact on these targets and their reduction has become a major issue in engine development.However, up to now, the complex non-linear cause-and-effect chain starting from the intake flow and finally leading to cyclic variations of the combustion process is not yet fully understood.
Numerical simulations using URANS, often in combination with experiments, is an established method in the engine development process.However, URANS is not capable to describe CCV and only Scale Resolving Simulations (SRS), e.g.Large-Eddy Simulations (LES), have the potential to capture cyclic fluctuations.One major drawback for the application of SRS are the significant computational demands of such investigations.The high spatial and temporal resolution for each individual cycle and the necessity to calculate a statistically relevant number of cycles pushes the computational requirements of SRS orders of magnitudes beyond what is needed for URANS.Thus, SRS for IC engine flows is still a topic of research rather than a standard method in the engine design and development process.
One important issue, especially in the context of numerical simulations, is the characterization and sometimes even the definition of CCV.While in experimental investigations CCV can be identified clearly for example by global quantities such as the cyclic variations of the peak and/or the mean pressure, such definitions are less straightforward when looking at local numerical results based on a limited number of cycles (being usually smaller than the number of measured cycles).In the following, a heuristic approach for the identification of CCV based on the velocity field in IC engines is proposed, which allows to identify fluctuations of large-scale (coherent) and smaller structures.It is important to note that there is an active discussion concerning the separation between large-scale variability and turbulence already mentioned in [1].This study is not meant to enter into this debate or to provide a clear definition of the separation or the associated length or time scales.As already mentioned by Lumley [2], there is no scale separation in engine flows.Rather it is our aim to use the well-known capability of POD to identify large scale energetic structures to quantify the cyclic fluctuations of the coherent motions with respect to a ensemble average.Furthermore, the remaining fluctuations are expected not to be associated with a variation of a coherent motion but rather being smaller scale turbulent fluctuations.To achieve this, our method combines POD, which has been used previously for the identification and investigation of velocity field fluctuations based both on experimental and numerical data [3][4][5][6][7][8][9], and conditional averaging [10].
In the next sections, an overview of the general workflow to identify CCV using the SRS data and the modeling approach is given.Afterwards, the presented workflow will be tested for two engine applications, a simplified engine configuration [11] and the well-established benchmark Transparent Combustion Chamber (TCC) engine [12], respectively.

IDENTIFICATION OF CYCLIC VARIATIONS
In the following, an introduction of the proposed workflow to identify CCV will be given (Sect.1.4).The workflow is based on conditional averaging (presented in Sect.1.3), with the cycle-dependent POD coefficients (Sect.1.2) as conditioning variables.The subsequent analysis uses the velocity field.

Definition of Cyclic Variations
CCV in IC engines are a well-known but not fully understood phenomenon.Compared to the effects of CCV, like the variation of the peak pressure during combustion being directly accessible both in experiments and simulations, the relevant underlying physics are not yet fully identified despite a significant number of investigations, [13][14][15][16][17] and references therein as well review articles [18,19].
Although the term CCV is used frequently, to the authors' knowledge there is no generally accepted and unique definition.Here, we define CCV of the velocity field based on the following considerations: each cycle is influenced by intrinsically (e.g.turbulence) and extrinsically (e.g.acoustic effects) induced fluctuations and represents an individual process; at any instant in time within a cycle, the specific flow field directly results from these fluctuations; each of these specific flow fields is based on a less complex fundamental flow topology, which is unique and should be clearly identifiable.The fundamental flow topology is associated with a large-scale motion; cyclic variations of these fundamental flow topologies are considered to be CCV; additional fluctuations, e.g. from non-coherent structures, do exist but are not considered as CCV here; In summary, this means that for IC engines, CCV of the velocity field are characterized here by the cyclic variation of large-scale coherent structures such as e.g. the intake jet or the tumbling motion.Of course, the smaller scales do fluctuate from cycle to cycle as well but our focus in this investigation is on the coherent structures.The aim is to extract these underlying flow topologies using a suitable mathematical procedure and quantify their variations denoting this as large-scale fluctuation.Unfortunately, there is no straightforward and unique way to extract an underlying flow topology from a single flow field realization.To identify these flow topologies, a new method based on POD and conditional averaging is proposed and evaluated.It should be mentioned that recently a cluster-based analysis [20] was proposed for IC engine flows, which seems to open up a very promising route but is not the scope of this study.

Proper Orthogonal Decomposition
Proper Orthogonal Decomposition (POD) is a powerful and well-established method to analyse complex processes.By means of a linearisation, a low-dimensional approximation of the problem can be obtained, which is often easier to handle.POD has been successfully applied in model reduction and data analysis in several scientific fields.In fluid dynamics, POD was introduced by Lumley [21] to identify coherent structures in turbulent flows followed by a number of investigations of complex flow fields in several applications [22,23].
Here, we focus on the snapshot POD method [6,7,24].A detailed description of the underlying theory and the application for experimental and simulation data can be found in [21,24,25].Starting point for the POD is an approximation of an arbitrary function f by a sum of weighted linear basic functions: To obtain such an approximation, several approaches such as the Fourier decomposition exist.For POD, the base functions W k ðxÞ have to fulfill two requirements.In addition to orthogonality of the basis function W k ðxÞ, the approximation of f has to be as good as possible in a least square sense for any arbitrary k.The singular value decomposition corresponding to Equation (1) reads as follows: Matrix A represents the discrete form of function f(x,t) and includes the data set under investigation.In general, there is no explicit rule how the data in matrix A in terms of rows and columns should be arranged.However, depending on the data set structure, a particular arrangement can reduce the computational costs to solve Equation (2).For a typical numerical investigation, the number of points in space is usually large compared to the number of points in time and this also applies for the data sets discussed later.Thus, the following structure is used here: each row of A represents a single variable over the considered time period; each column of A represents a single time step and contains all variables of interest.This results in an M Â N matrix A, while the number of rows directly results from the absolute number of considered values in the complete domain of interest and N represents the number of time steps ("snapshots'').The columns of the orthogonal M Â M matrix U , which represents W k ðxÞ in Equation (1), are the proper orthogonal modes.The diagonal elements R ii , also called singular values, represent the energy content of the individual modes and contain only non-negative values.All off-diagonal elements of the M Â N matrix R are zero.Each row of the orthogonal N Â N matrix V corresponds to a POD mode and acts as a time-dependent weighting factor describing the dynamics of the mode.The factor a k ðtÞ in Equation ( 1) is obtained from the product of the row and column of R and V T .
In this work, an eigenvalue decomposition of a temporal correlation matrix using the data set under investigation is used to solve the singular value decomposition.For this reason, Equation ( 2) is multiplied with its transpose, see Equation (3): For a data set only consisting of real values, as it is the case here, it follows for the orthogonal matrix U T ¼ U À1 , leading to: Equation ( 4) can be solved by an eigenvalue decomposition and the eigenvalues k are the squared singular values R ii .Afterwards, the matrix U , consisting of the proper orthogonal modes, can be calculated from Equation (2).Finally, the rows and columns of the matrices are sorted in descending order (in terms of the energy content i.e. singular values).This means that the first mode is the highest energetic and the last mode is the lowest energetic proper orthogonal mode.
It is interesting to note that POD has been applied regularly in the analysis of IC engine flows to obtain a better understanding of cyclic variations [3,5,7,8,14,17,26].Based on these investigations, there is a broad consensus that the first (highest energetic) mode represents the ensembleaveraged flow field.Furthermore, it has become common practice to associate the other modes to fluctuations, while the modes with more energy are often associated with large-scale fluctuations.The scope of this study is to extend this idea by finding within these higher modes a number of flow topologies, which can be used for conditional averaging.

Principles of Conditional Averaging
Statistical analysis is a fundamental part for the investigation of turbulent and complex fluid flows.The use of this methodology is discussed in the context of IC engine analysis.Instead of averaging all available flow fields, the calculation of mean values and statistics is intuitively limited to flow fields, which are identical in terms of the crank angle.The advantage of this methodology is that the resulting average is based on instantaneous flow fields, which are similar in terms of the geometry and its boundary conditions (piston speed) and accordingly to the fundamental flow topology.The ensemble average is calculated as: while the components of velocity fluctuation u 0 i ðx; u; nÞ are determined by the difference of the instantaneous and the ensemble averaged flow field.For simplicity, in the following the dependency on space and crank angle is dropped unless necessary.However, a straightforward use of this averaging procedure and decomposition can lead to less meaningful results in terms of large-scale fluctuations, which can be linked to CCV as described in Section 1.1.
To remedy this problem, the average can be calculated based on a reduced set of instantaneous flow fields, which are identical in terms of a conditioning vector /, whereby this conditioning vector must be able to detect the cycle-dependent dynamics of the large-scale structures.This method of conditional averaging [10] is defined according to: As can be seen, Equation ( 8) is very similar to Equation ( 5), while the value M ðUÞ denotes the number of instantaneous flow fields for which / ¼ U is valid.
The high-energetic POD modes are linked to large-scale structures [3,7,26], therefore the cycle-dependent coefficients of these modes (matrix V in Sect.1.2) offer a starting point to define the conditioning vector /.For this specific application here (and this will be shown later), it is possible to reduce the conditioning vector to a single conditioning variable.This conditioning variable can be represented by a set of four combinations of POD coefficients.These combinations, called subsets below, are denoted by the simple variable sðnÞ, whereby the dependency of the individual cycle n is not explicitly shown in the following for notational convenience.Furthermore, the conditioning variable s can take the discrete values 1, 2, 3 and 4 (see details below).Accordingly, the subset-dependent averaged flow field is denoted as ûðsÞ.
Based on this conditional average, the following decomposition of the instantaneous flow field becomes feasible: Compared to a Reynolds decomposition (Eq.6), which results into a cycle-independent averaged flow field u and cycle-dependent fluctuations u 0 ðnÞ, a subset-dependent averaged flow field ûðsÞ and cycle-dependent fluctuations u Ã ðnÞ appear in Equation (9).Averaging the resulting subsetdependent flow fields yields the ensemble average: where aðsÞ, defined as: is the ratio of independent samples within each subset and the overall number of independent samples.The quantification of the subset specific small-scale fluctuations is performed similarly to Equation ( 7): whereby the subset specific small-scale fluctuations u Ã i ðnÞ are calculated by Equation ( 9).Please note that each cycle n is associated with one subset s.The root mean square of the coherent structure fluctuations is: and this can be seen as an indicator for large-scale CCV.
The fluctuations are defined by the difference of the globally averaged flow field (Eq.5) and the subset average (Eq.8) with an overall number of S subset averages (for the applications here, S ¼ 4).

Workflow to Identify Cyclic Variations
In this section, based on the above ideas, the procedure for the identification of cyclic variations in IC engines is described.As mentioned before the starting point of this method is the established capability of POD to identify coherent structures.A significant portion of an engine cycle is characterized by high-energetic flow structures (e.g.intake jet and large vortex structures), which are associated with high energetic POD modes.Cyclic variations of these characteristic high-energetic flow structures are reflected by a variation of the high-energetic POD coefficients from cycle to cycle.For that reason, these coefficients, which act as a "cycle-dependent weighting factor" for the proper orthogonal modes, provide a good base to cluster cycles into subsets.Following the discussion above and the definitions introduced in Section 1.1, each subset is defined in such a way that it includes all independent samples with an almost identical fundamental flow structure (i.e. the structure of the most energetic flow patterns).Averaging within each of these subsets yields the fundamental flow topology of each subset.The deviation of the subset average ûðsÞ to the ensemble average u does represent large-scale cyclic variations.The general workflow is shown in Figure 1.
It is interesting to note that this approach, in contrast to previously published methods, does not require an interpretation of the resulting POD modes.This is an important point because POD is a purely mathematical method and the resulting modes do not necessarily allow a physical interpretation.However, the conditionally averaged flow fields are intuitive in terms of the flow physics and always allow an interpretation.
In a previous investigation [27], we looked at the vortex shedding behind a triangular obstacle.We revisited this configuration (not shown here) as a first step to evaluate and test our workflow according to Figure 1 for the identification of coherent structures representative for the vortex formation at the edge of intake valve during the intake stroke, which we focus on below.Starting with two subsets (based on the first POD mode) the number of subsets was increased by adding two subsets (based on the highest energetic mode not used up to this point).The conclusions from this simple test case were: the coefficient of the first POD mode (RANS-like field) is almost constant with respect to the different samples.Thus, it is not suitable in terms of a subset classification; the coefficients of the other POD modes vary in time and can be used for subset classification; vortex shedding can be captured with even a small number of conditional averages/subsets based on the high energetic POD modes; although the use of more subsets leads to an improved and more detailed information about the large-scale  structures, 4 subsets based on mode 2 and 3 provide a reasonable representation of the vortex shedding.It has to be mentioned that this finding might differ for other, especially more complex, flow situations.Of course, in practical applications, there is only a limited number of independent samples available.This leads to a limited number of subsets, which is directly linked to the number of considered POD modes.Thus, every subset can only be an approximation of the underlying flow topology of an individual cycle.For this reason, a compromise has to be found between: the statistical convergence within each subset; the quality of the fundamental flow topology approximation.
Based on the discussion above, for the cases considered in the following, we restrict ourselves to four subsets as a compromise between statistical convergence and representation of the underlying flow topology of the subsets.

NUMERICAL METHOD AND TURBULENCE MODELING
In this section, an overview of the CFD solver, the numerical method and the used turbulence model will be given.

Numerical Method
All simulations are performed using ANSYS CFX Release 15.0, which is an implicit and node centered code.The governing equations are solved based on a finite volume method.For fully compressible flow (e.g.Sect.3.2), the governing Equations ( 14) to (16) for mass, momentum and enthalpy are Favre-filtered: where q i is the molecular heat flux.The quantity H, defined as: describes the total enthalpy, while s ij is the molecular stress tensor.The quantities q mod j and s mod ij are associated to turbulent fluctuations and have to be modelled.Furthermore, the system is closed by the equation of state for perfect gases: For an incompressible fluid, which is the case for the simplified engine setup described in Section 3.1, only a reduced set (no enthalpy equation) is considered and no Favre-filtering is applied.For all simulations, we use a central differencing scheme in space and a second order backward Euler scheme in time.
A crucial topic for the simulation of IC engines is the rather complex mesh generation taking into account the moving piston and valves, which require special treatment.
Here we use the so-called key grid approach.Starting from an initial mesh, the grid is deformed according to the moving boundaries until a critical mesh quality is reached.At this point, a new mesh is generated and the previous result is interpolated onto it.Both, the implications of the mesh morphing and the interpolation in ANSYS CFX was investigated in detail by Hasse et al. [16].

SAS-SST Turbulence Model
The Scale Adaptive Simulation Shear Stress Transport (SAS-SST) turbulence model was introduced by Menter and Egorov [28].The fundamental idea is derived from the exact transport equation of the two-point correlation developed by Rotta [29].This transport equation is linked to the turbulent kinetic energy k and the integral length scale L and includes the von Kármán length scale.This length scale, including the second derivative of the velocity field, was finally introduced in the x-equation of the well known SST turbulence model [30].The modified x-equation for the SAS-SST model reads: The last term on the right-hand side defines the transformation of e to x and is a fundamental part of the SST model.The additional source term: includes the von Kármán length scale, which is defined as: Due to this additional source term, the SAS-SST turbulence model has the capability, unlike URANS models, to resolve fluctuations smaller than the integral length scale [31,32] in case of sufficient spatial and temporal resolution.For an insufficient resolution in space and time the model reverts back to standard SST behavior.

ENGINE APPLICATIONS
In this section, the methodology presented above to identify and quantify CCV will be tested on two different engine setups.The first case, based on a simplified engine geometry, is meant to investigate the workflow in general in a welldefined and geometrically simple setup, which however still retains a significant number of engine specific characteristics.Numerical studies of this engine have been published by several authors [33][34][35].
In a second step, the method will be evaluated using data from a realistic engine case.The analysis will be performed for two specified crank angles.
Both engines are motored and therefore no injection, ignition and combustion are considered.It is important to note that we focus on fluctuations of the intake jet, sometimes referred to as "jet flapping", characterized by the flow separation at the intake valve and large-scale vortical flow structures.As explained below, this has a significant impact on the choice of the conditioning variables.

Simplified Two-Stroke Engine
This case is based on the engine setup experimentally investigated by Morse et al. [11].Due to the simplicity of the geometry, this setup is suitable for the first analysis.The direct numerical simulations by Schmitt et al. [34] provide a very detailed insight into the complex flow structure inside the combustion chamber and confirmed an axisymmetric average flow field as well as axial symmetry for the coherent large-scale and small-scale structures.The data sets are generated from slices rotated around the cylinder axis with a non-dimensional radius of 0:2 to 1 (Fig. 2).Although the slices and thus the POD snapshots are only 2D, for a sufficiently large number of samples, are relevant flow structures should be captured since the azimuthal direction is homogeneous.The average distance between the adjacent samples is chosen such that it is more than one integral length scale, see Figure 2.This yields 26 samples for each cycle.A total number of nine consecutive cycles for the analysis, leaving out the first cycle, leads to a total number of 208 samples considered below.

Numerical Setup
The geometry features a straight intake pipe, a flat dome and piston and a single non-moving valve.Although the experiments were performed in an horizontal arrangement, Figure 2 and all subsequent figures use the more common vertical arrangement.
All parts are concentric with respect to the cylinder liner with a bore diameter of 75 mm.The engine runs at 200 rpm with a sinus-shaped piston motion zðuÞ, which results in a mean and maximum piston speed of 0.40 m/s and 0.63 m/s, respectively.A compression ratio of 3 results from a clearance height of 30 mm and a piston stroke of 60 mm.Because of the non-moving valve, there is a constant valve gap size of 4 mm and a full engine cycle corresponds to 360 °CA with an intake and exhaust stroke of 180 °CA each.Corresponding to the experimental setup, the simulations were performed with air at an ambient pressure and temperature of 1.03 bar and 300 K. Resolved turbulent fluctuations were superimposed on the filtered inlet boundary condition.
The simulations are carried out on a highly-resolved hexahedral mesh with about 8 million grid points and a mesh refinement at the wall leading to a minimum edge length of 40 lm.Due to the non-moving valve and the simple geometry, the complete cycle can be simulated on a single mesh preventing interpolation errors.The necessary mesh morphing is limited to the combustion chamber.The time step width was set to 0.025 °CA, which leads to a CFL number smaller than 1 for the largest part of the simulated domain.

Results
Figure 3a shows the instantaneous velocity field on the mid plane at 90 °CA with a maximum velocity magnitude of about 12 m/s in the intake jet region.In addition an isosurface for a relative pressure of p rel = À1 855 Pa is shown.
Figure 3b shows the ensemble-averaged velocity field (left half plane) as well as the ensemble-averaged rms values of the velocity fluctuations (right half plane) extracted from slices as shown in Figure 2. The averaged velocity field shows a well-defined intake jet as well as a vortex structure corresponding to the isosurface illustrated in Figure 3a.As expected, the velocity fluctuations show a maximum in the shear layer region of the intake jet.It is interesting to note that the rms values also have a local maximum close to the vortex core region of the averaged velocity field (Fig. 3b) although no significantly increased level of fluctuations can be recognized in the independent samples (not shown here).This is a first indication that this local maximum, which occurs for wide ranges of the intake stroke, does not result from turbulent small-scale fluctuations but rather from a large-scale structure fluctuation, which is represented by the deformation and displacement of the large vortex structure with respect to the cylinder axis.Thus, the high rms values directly result from the ensemble-averaging procedure.Borée et al. [5] identified experimentally such a large scale jet fluctuation in another simplified engine setup.Hasse et al. [15] performed scale-resolving simulations of this setup and concluded in agreement with the experimental observations, that fluctuations of the large scale structure are falsely identified as turbulent kinetic energy when standard ensemble-averaging is used.These observations are another indicator that the standard averaging procedure might lead to results difficult to interpret and in fact, our method aims to reduce this issue.Figures 3c and d show the averaged and the instantaneous axial velocities (normalized by mean piston speed) on lines (normalized by cylinder radius) placed at 10 mm and 30 mm below the cylinder head.Furthermore, the standard deviation of the instantaneous velocities (related to the averaged velocity) is shown as error bars.The largest deviations for the instantaneous velocity profiles at -10 mm can be identified at a nondimensional radius of 0.6 to 0.8, which represents the unstable jet position during the intake stroke.At -30 mm a broad region between 0.45 and 0.95 with large fluctuations can be identified.

Identification of Cyclic Variations
To identify large-scale fluctuations, the workflow described in Section 1.4 is applied at 50 °CA.The first step is a POD analysis using all 208 snapshots.Figure 4 shows the POD modes (each normalized by its maximum value) resulting from the velocity field at 50 °CA.
As mentioned before, there is an obvious similarity between the first (i.e.highest energetic) POD mode in Figure 4a and the ensemble-averaged velocity field in Figure 5a.The interpretation of the 2nd, 3rd and 4th mode is less clear at this point due to the more complex flow topology.It seems that especially the 2nd and 3rd mode are strongly connected to the intake jet and the vortex structure due to the dominating influence of the jet shear layer and the vortex core close to the cylinder axis.In connection with the POD coefficients, it seems that these modes are linked to the jet direction and its width as well as the position and strength of the vortex structure of the flow field in physical space.An interpretation of the 4th mode appears more difficult at this crank angle due to missing dominant structures and associated field topologies.Most likely, this mode can be associated to the vortex core.
As a next step, the instantaneous flow fields are grouped into different subsets followed by conditional averaging.As described in Section 1.4, this classification is based on the POD coefficients.The relevance of the individual POD modes can quantified by their singular values, which are shown in Figure 6.
As expected the singular value of the first mode is high compared to the other modes.For wide ranges over the intake stroke, the singular value of the second mode is also comparatively large, which is an indicator for a large-scale structure in physical space and the higher modes seem to contain only a smaller fraction of the energy.This suggests that these modes are connected to less dominant flow structures in physical space.
For the further discussion, we use the following observations, strictly valid only for these crank angle: due to the structure of mode 1, 2 and 3, see Figure 4, these modes can be associated with the intake jet; for most of the intake stroke, the singular values of the first and second mode are high compared to the singular values of higher modes (Fig. 6); the POD coeficient of the first mode is almost constant from sample to sample and thus less useful for a classification, see Figure 7.This is similar to the findings described in Section 1.4; the POD coeficients of the other modes show strong variations from sample to sample (Fig. 7; only coefficients of 1st, 2nd and 3rd POD modes are shown for readability).
Based on this observation and the necessity of a compromise between the number of subsets and the samples within each subset, as described in Section 1.4, in the following we focus on a small number of high energetic modes for the analysis of large-scale fluctuations.This results into four subsets, defined by the coefficient of 2nd and 3rd POD mode, which ensures a significant number of individual samples necessary for statistical convergence within each subset, Table 1.It is important to note that this is a restriction of the current analysis and we cannot claim that the general applicability reaches beyond this jet intake flow configuration.Certainly further investigations are necessary to reach more generally applicable procedures.Nevertheless, as will  be shown below, this very straightforward approach (shown in Tab. 1 and Fig. 7) chosen here can provide an interesting physical insight and offers a different way of looking at large scale fluctuations.
Beside the classification criterion, Table 1 also shows the number of instantaneous velocity fields clustered into each subset.In a final step, an averaging within the subsets leads to the subset averages conditioned by the POD coefficients.Figure 5 shows the ensemble-average as well as the subset averages at 50 °CA.
As expected, the general flow topology is characterized by an intake jet and a vortex structure, which can be identified for the ensemble average as well as for the subset averages.An explicit comparison of the averages is illustrated in Figure 5f by means of velocity isolines of 6 m/s (intake jet) and 1 m/s (vortex structure), respectively.The intake jet is identical for all averages in the region close to the valve gap.This is the direct result of the strong wall influence on the flow.However, a variation of the jet direction and the penetration depth can be clearly seen in the combustion chamber.The intake jet of the 1st and 3rd subset average in Figures 5b and d show a increased penetration depth compared to the ensemble average in Figure 5a.In the direct comparison the averaged velocity fields of the 2nd and 4th subset in Figures 5c and e show a larger deviation from the ensemble-averaged velocity field with respect to the jet direction.
A variation can also be identified for the vortex structure.With the exception of the 3rd subset, each subset average shows a clear variation of the vortex core position compared to the ensemble average.It can be clearly seen in Figure 5f that the shift of the vortex core position is directly linked to the position of the intake jet.For example, in the second subset the most significant bending of the intake jet towards the cylinder axis can be identified.The corresponding vortex core is shifted significantly towards the cylinder axis.Since the vortex core and the jet are strongly coupled, their observed fluctuations are also well correlated.
Figure 8a shows the rms of the velocity fluctuations calculated by the ensemble-averaged and the total number of instantaneous velocity fields according to Equation (7).Close to the valve gap, the velocity fluctuations are observed in a well-defined and thin shear layer generated by the intake jet.
Figures 8b-e show the rms values calculated from the subset averages and the corresponding instantaneous velocity fields as described in Equation (12).Most notable is that the velocity fluctuations are very well-defined with respect to the subset averages in Figures 5b-e.As in Figure 8a, each subset shows very thin regions of velocity fluctuations in the vicinity the valve gap, corresponding to the thin shear layer.In contrast to the ensemble average, the fluctuations of the subset averages are constrained to the shear layer region near to the intake jet.An increased level of velocity fluctuations can be seen near the piston or the cylinder liner, which causes the final jet breakup.The well-defined structure of velocity fluctuations indicates that the subset averages are a good approximation for the fundamental flow topology representative for the corresponding instantaneous velocity fields.Figure 8f shows the large-scale fluctuations, calculated by Equation ( 13), which represents the as defined in Section 1.1.The structures of the large-scale fluctuations look very similar to the 2nd POD mode although the conditional averaging was done based on the POD coefficients of the 2nd the 3rd mode.This quite remarkable effect can be explained by the energy contents of the modes.As can be seen in Figure 6 the second POD mode shows an increased energy level of about 75% compared to the third POD mode and represents a more dominant flow field topology (i.e.large-scale structures).Another interesting observation is the occurrence of the thin region, bounded by the high rms values, which marks an overlap between all the subset and the ensemble average.The ensemble and subset-averaged velocity fields at 90 °CA, which corresponds to the crank angle of maximum piston speed, can be seen in Figures 9a-e, respectively.In general, the same trends as for 50 °CA can be observed.The downwards moving piston induces an intake jet.Due to the constant height of the valve gap, the velocity of the intake jet and the piston speed are directly coupled and also the maximum jet velocity is found at this crank angle.The jet direction remains almost unchanged until it is finally deflected by the cylinder liner.The direct comparison to 50 °CA shows that this jet bending is pushed further away from the axis.It is interesting to note that the influence of the liner seems to reduce the cyclic variations of the intake jet, see Figure 9f illustrating the conditionally and ensembleaveraged positions of the jet and the vortex, respectively.This liner influence will also be observed in the second example below.
The vortex structure appears more regular when compared to 50 °CA.The subset averages of the vortex core position seem to indicate a precession around the ensemble-averaged vortex core.Looking at the ensembleaveraged velocity field, this precessing vortex core results in an increased level of turbulent fluctuations depicted in Figure 10a and this finding agrees well with the investigations of Borée et al. [5].In contrast to the ensemble-averaged data field, the subset averages yield more homogeneous fluctuating fields within the combustion chamber and no local maximum at the vortex core, Figure 10b, is found.Thus, the potentially false identification of a turbulent fluctuation due to the precessing vortex core can be avoided using the proposed subset averaging.Figure 10c shows the large-scale fluctuations calculated by Equation (13) and Figure 10d shows the 2nd POD mode.Interestingly, the structures of the large-scale fluctuations and the 2nd POD mode appear to be almost identical, which can be explained by the higher energy content (more than 50% compared to mode 3) of the 2nd mode.

TCC Engine
The second application considered here is a realistic engine setup.The considered case is motored and a full cycle takes 720 °CA, while the firing TDC is at 0 °CA.Due to the fact that this engine has an straight vertical intake port without machined tumble edge, the resulting tumble motion is relatively small.For this reason, we focus on the intake phase (358 °CA-593 °CA) in the following.120 cycles/ ensembles of the TCC engine [12], being an internationally well-established benchmark case, were simulated.The optically accessible combustion chamber provides a highquality experimental data base for comparison to numerical investigations.

Numerical Setup
The domain for the numerical simulation is shown in Figure 11.
The setup features a flat piston and dome and a single intake and exhaust valve.The bore and the stroke are 92 mm and 86 mm, respectively.An effective compression ratio of 8 is based on the clearance height of 9.44 mm.All simulations were performed at the same atmospheric conditions as in the experiment.The engine runs at 800 rpm in motored operation.
As can be seen in Figure 11, the simulations use a hybrid tetrahedra/prism mesh with a maximum number of 2.5 million grid points during the phase when both valves are open.The time step varies during the cycle and reaches a minimum value of 0.25 °CA during the early phase of the intake stroke.Resolved turbulent fluctuations (similar to Sect.3.1.1)were imposed additionally as inlet boundary condition in combination with the SAS-SST model (Sect.2.2).

Results
Figures 12 shows the instantaneous velocity fields in the middle plane of cycle 6, 51 and 81 at 450 °CA.
Due to the orientation and shape of the intake port, the engine does not produce a significant tumbling motion.Indeed, the fluid motion shows no large-scale structures other than the intake jet with a velocity magnitude of about 50 m/s.After deflection, the intake jet is mainly following the direction of the cylinder liner and only small variations from cycle-to-cycle can be seen, which is consistent with the finding above for the simplified engine at 90 °CA.Furthermore, no other dominant structures are found in the rest of the flow field, which mostly features small-scale structures.Following up on the discussion above for the first test case, it is interesting to see whether the proposed method is also capable of describing this different situation with mostly small-scale and only limited large-scale fluctuations.

Identification of Cyclic Variations
The following investigation is carried out at 450 °CA during the intake stroke.The conditioning variables for this analysis are chosen in a fully analogous manner as in Section 3.1.3.A POD analysis based on the full 3D velocity field in the combustion chamber is performed.Identically to Section 3.1.3the POD coefficient of the first mode is almost constant with respect to the cycle count.Thus, the 2nd and 3rd mode are used for a flow field classification.Due to the relatively small number of 120 independent samples/cycles, the averages generated by the four subsets are potentially less reliable in a statistical sense than in Section 3.1.3.After classification, the averaging procedure, see Equation ( 8), within the different subsets is applied and the rms of the large-scale fluctuations according to Equation ( 13) is calculated.
Figure 13 compares the ensemble-averaged velocity field with the rms of the large scale fluctuations and the conditional averages of subsets 1 and 4. At a first glance, the ensemble average and the subset averages 1 and 4 appear almost identical.However, there is a notable deviation from subset to subset, which is found in close proximity of the intake jet region, which is not directly influenced by the cylinder liner and correspondingly the largest values for the rms of the large-scale fluctuations are found in this specific region.In the wall-guided part of the intake jet, the differences are very small.
In Figure 14, the velocity magnitude of selected averaged and instantaneous velocity fields along a line in the middle plane 30 mm below the cylinder head (Fig. 11) is illustrated.A direct comparison of the ensemble and subset-averaged velocity fields is shown in Figure 14a.In the region of the wall-guided intake jet, all averages are almost identical.A larger deviation can be identified for the position of the maximum velocity magnitude and subset 4 shows the largest deviation with respect to the ensemble average.
Figures 14b and c show the averaged as well as the instantaneous velocity profiles of subset 1 and 4. In general, both subsets confirm that the instantaneous velocity fields are classified correctly in terms of the dominant structure and thus the instantaneous profiles can be associated with the corresponding subset average, e.g.looking in detail at the difference in velocity magnitude near the cylinder axis or the shift of the maximum velocity in terms of the dimensionless radius.Instantaneous velocity magnitude in the TCC engine in the middle plane at 450 °CA (i.e.intake stroke).

SUMMARY
A new methodology to identify and quantify large-scale CCV in internal combustion engines is proposed.The method is applicable both for experimental and numerical data.
Here, it is used and evaluated for 2D and 3D velocity data sets from multi-cycle scale resolving simulations of two different engines.The general idea behind the method is to identify the large-scale coherent motion and its deviation from the ensemble average.This is achieved by first decomposing the instantaneous velocity field into a mean and a large (energetic)-as well as a smaller-scale contribution to the fluctuations, respectively.While the mean field is determined by a standard ensemble average, the identification of the two fluctuating quantities is based on a combination of POD and conditional averaging.POD is known to have the capability to identify coherent structures, and we use combinations of the cycle-dependent POD coefficients of the first modes to identify suitable subsets, which are then the conditioning variable.
The instantaneous flow fields are then grouped into different subsets with the aim of averaging and calculating fluctuations within these subsets.Since the number of realizations within each subset is smaller than the overall number of realizations, this averaging potentially requires a larger number of engine cycles.The difference of the subset average to the ensemble average is then denoted as large-scale fluctuation.On the other hand, the difference between the instantaneous velocity field with respect to its specific subset average is denoted a small-scale fluctuation.
This method was tested using two numerical data sets from scale resolving engine simulations.The focus was to assess the ability to detect large scale flow topologies and their cyclic fluctuations in addition to small-scale fluctuations.For the first engine test case, 208 independent samples were available and 4 subsets could be identified based on the POD coefficients of the first modes.The results show a clear difference of the position and penetration depth of the intake jet for each of the 4 subsets and this is directly identified by the method.Furthermore, a variation of the vortex core position, similar to a precessing vortex core, can be detected.Small-scale structures are found in regions, where high turbulent fluctuations are expected, e.g. in the high shear region of the jet.
The second case is the well-established TCC engine benchmark, for which previously computed 120 cycles were used for the analysis.Again, 4 different subsets could be identified and the same averaging procedure was performed.Similar conclusions concerning the decomposition of the velocity field as for the first case could be drawn and detailed results were shown.
Finally, it can be concluded that the newly proposed method allows to separate fluctuations of the velocity field into a large-scale coherent and a small-scale turbulent contribution.Thus, using this method allows to distinguish the large-scale fluctuation, such as a precessing vortex core, from the classical turbulent contributions.

Figure 1 General
Figure 1 General workflow to calculate conditional averages and fluctuations in context of IC engines at a specified crank angle.

Figure 2 a
Figure 2 a) Side and b) top view of hexahedral engine mesh and engine geometry (measurements in millimeter).Red lines illustrates the location of slices for data analysis.

Figure 3 a
Figure 3 a) Instantaneous velocity magnitude in the mid plane and isosurface at relative pressure of À1 855 Pa in cycle 2; b) shows the ensemble-averaged velocity field (left half plane; legend as a)) and the root mean square of velocity fluctuations (right half plane), c,d) show the averaged (black line) and instantaneous (grey lines) axial velocities (normalized by average piston speed) as well as the standard deviation of velocity (bars) at À10 mm and À30 mm.All results are shown for 90 °CA.

Figure 6 Figure 5 a
Figure 6 Singular values of the first 4 proper orthogonal modes as a function of crank angle.0 °CA to 180 °CA for intake and 180 °CA to 360 °CA for exhaust stroke.The two crank angles, 50 °CA and 90 °CA, respectively, discussed in Section 3.1 are marked.

Figure 4 Proper
Figure 4Proper orthogonal modes 1 to 4 (each normalized by its maximum value) resulting from the velocity field at 50 °CA in the simplified engine setup.

Figure 7
Figure 7 Sample dependent POD coefficients of 1st, 2nd and 3rd mode at 50 °CA.Furthermore, exemplarily a subset classification (corresponding to Tab. 1 is shown).

Figure 9 aFigure 8
Figure 9 a) Ensemble-averaged and b-e) conditionally-averaged velocity fields at 90 °CA in the simplified engine setup.f) Comparison of the conditionally-averaged (thin lines) and ensembleaveraged (thick line) using velocity isolines (6 m/s for intake jet and 1 m/s for vortex core).Corresponding structures are indicated by the same color.

Figure 10 a)
Figure 11 a) Geometry and b) numerical mesh of the TCC.

Figure 13 a
Figure 13 a) Ensemble-averaged velocity fields and subset-averaged velocity fields b,c) in the TCC engine at 450 °CA; d) shows the root mean square of the large-scale fluctuation according to Equation (13).In addition two isosurfaces for the velocity magnitude (50 m/s) and for ũrms (1.5 m/s) are included in all plots.

Figure 14 a)
Figure 14 a) The velocity magnitude of the ensemble and subset-averaged flow fields on line z = 30 mm below cylinder head at 450 °CA; b) and c) the averaged as well as the instantaneous velocity magnitudes of subset 1 and subset 4.