Marcellin Perceau
email: marcellin.perceau@renault.com

Philippe Guibert

Stéphane Guilain

Analysis of a

Keywords: Dual PIV, Clustering, k-means algorithm, Miller Cycle

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

A large-scale tumbling motion is generated during the intake phase of most gasoline internal combustion engines. It enables the storage of kinetic energy which will be transformed into small scale movements during the combustion. These movements, usually called turbulence, homogenize the air-fuel mixture and wrinkle the flame front. This accelerates the combustion process and makes the engine more efficient (John L. [START_REF] Lumley | Engines: An Introduction[END_REF][START_REF] Ferguson | Internal Combustion Engines: Applied Thermosciences[END_REF].

Recent years have seen the development of the Miller cycle [START_REF] Miller | High-pressure supercharging system[END_REF]1956;[START_REF]Supercharged engine[END_REF][START_REF] Sporleder | The Efficient Gasoline Engines in the New Opel Astra K[END_REF][START_REF] Budack | The New Audi 2.0-l I4 TFSI Engine[END_REF][START_REF] Demmelbauer-Ebner | The New 1.5-l Four-Cylinder TSI Engine from Volkswagen[END_REF], which increases the efficiency of thermal engines (Perceau et al. 2020;Perceau, Guibert, and Guilain 2020;2021b). The improvement comes from its extended expansion stroke, which allows more work to be recovered for the same amount of fuel. This is achieved by decoupling the compression and expansion phases by closing the intake valves before the piston reaches the bottom dead center (BDC). The engine thus has a longer expansion than the effective compression stroke. However, unlike a classic Otto cycle, which closes the intake valves shortly after BDC, the premature closing of the intake valves here disturbs the tumble motion. Indeed, the tumble movement is no longer maintained by the valve jet before the BDC.

The question then arises as to whether the design phase of a new Miller engine can be done in the same way as in the case of an Otto engine. Most of the time, an intense movement in the middle of the intake phase or close to the BDC gives an intense movement at the end of the compression stroke. For this reason, curves characterizing the intensity level of the tumble rotation as a function of the valve lift are often used to compare different cylinder head geometries. The tests are done on stationary test benches or numerically and the selected geometries are those generating the most intensity with the valves open. It is then necessary to check, by an experimental approach closer to a regular operation, that this method can still be done in the case of a Miller cycle. Indeed, generating a strong tumble is costly in terms of pressure losses at the inlet. If a strong tumble and a tumble of medium intensity tend to produce the same level of turbulence, it will be necessary to reconsider the validation process of the geometry of a new cylinder head. The aerodynamics will be studied here experimentally on a transparent cylinder engine associated with two particle image velocimetry (PIV) systems temporally linked. These systems are usually called dual PIV. They can measure two velocity fields within the same engine cycle and thus follow the temporal evolution of the in-cylinder flow structures. This type of system has been used for various applications [START_REF] Wernet | Dual PIV Systems for Space-Time Correlations in Hot Jets[END_REF][START_REF] Bridges | Effect of Heat on Space-Time Correlations in Jets[END_REF][START_REF] Souverein | Application of a Dual-Plane Particle Image Velocimetry (Dual-PIV) Technique for the Unsteadiness Characterization of a Shock Wave Turbulent Boundary Layer Interaction[END_REF][START_REF] Schreyer | Development of a Dual-PIV System for High-Speed Flow Applications[END_REF] and rarely for the study of the internal flow of combustion engines [START_REF] Guibert | Dual Particle Image Velocimetry for Transient Flow Field Measurements[END_REF]. PIV systems using fast lasers, in the kHz range and usually called HS-PIV or TR-PIV [START_REF] Druault | Use of Proper Orthogonal Decomposition for Time Interpolation from PIV Data[END_REF][START_REF] Roudnitzky | Proper Orthogonal Decomposition of In-Cylinder Engine Flow into Mean Component, Coherent Structures and Random Gaussian Fluctuations[END_REF][START_REF] Müller | Flow Field Measurements in an Optically Accessible, Direct-Injection Spray-Guided Internal Combustion Engine Using High-Speed PIV[END_REF][START_REF] Sick | High-Speed Imaging for Direct-Injection Gasoline Engine Research and Development[END_REF][START_REF] Voisine | Spatio-Temporal Structure and Cycle to Cycle Variations of an in-Cylinder Tumbling Flow[END_REF][START_REF] Buschbeck | Dependence of Combustion Dynamics in a Gasoline Engine upon the In-Cylinder Flow Field, Determined by High-Speed PIV[END_REF][START_REF] Zegers | Pre-and Post-Injection Flow Characterization in a Heavy-Duty Diesel Engine Using High-Speed PIV[END_REF], are generally preferred because of the higher number of fields that can be measured during a single engine cycle. They can thus be used to monitor the flow more closely. These systems also require only one laser, which means lower maintenance costs and simpler adjustments than dual systems. However, this type of system can't measure as many cycles as the dual PIV systems. Indeed, the high number of velocity fields measured for each cycle quickly fills the memory available on the acquisition system. The dual PIV system allows to monitor 300 engine cycles, which represents 300 couples of velocity fields in one experiment. This large number of instantaneous flow fields is therefore well suited for statistical studies and a clustering method, that classifies the dual-PIV data into different groups with similar properties, is chosen in this paper to describe the flow. A first part will focus on the center of rotation of the motion and will define a criterion to obtain strong rotations. A second part will classify the velocity fields in three groups representing strong, medium and weak flows. It will thus reveal whether strong rotation rates at the beginning of the Miller cycle generate strong turbulence rates at the end of the compression stroke. The classification is performed here using a k-means algorithm [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. It is beginning to be applied to the analysis of phenomena involved in internal combustion engines such as combustion chemistry [START_REF] Perini | High-Dimensional, Unsupervised Cell Clustering for Computationally Efficient Engine Simulations with Detailed Combustion Chemistry[END_REF], vortex pattern detection [START_REF] Zhao | K-Means Clustering-Driven Detection of Time-Resolved Vortex Patterns and Cyclic Variations inside a Direct Injection Engine[END_REF] or cyclic variations [START_REF] Cao | Cluster-Based Analysis of Cycle-to-Cycle Variations: Application to Internal Combustion Engines[END_REF] of the internal flow.

The paper is organized as follows. At first in Sect. 2 the experimental setup is presented. Sect. 3 presents the results of the clustering on the rotational centers and the velocity fields. The article is then ended by several conclusions and perspectives in Sect. 4.

EXPERIMENTAL SETUP

Model engine

The model engine was representative of a 0.9L three-cylinder spark ignition car engine. It was made of a single transparent cylinder. This offered wide optical access to study the aerodynamic movements inside the combustion chamber. Note that it did not allow for firing and that the rotation is provided by an electrical engine. The experiments are carried out at 1200 rpm in order to be tuned to the 10Hz frequency of the lasers and allow two observations per cycle by the dual-PIV system. This speed also avoids damaging the transparent engine while being close to the idle speed and providing good visualization of the aerodynamic movement. The cylinder head had four valves, two for the intake and two for the exhaust, in a pentroof geometry. All the measurements were made with an intake pressure of 1 bar. The engine had a 12.3 compression ratio, an 81.2mm stroke and a 72.2mm bore. A flat piston head was used. Indeed, a real piston head shape can induce distortions in the laser sheet that could lead to inaccurate positioning and reflections. This modification was made to maintain the value of the compression ratio. The engine was running in a Miller configuration, with an early intake valve closing at 170°CA (Crank Angle) and short opening duration of 140°CA.

Dual PIV measuring system

The measuring system consists of two synchronized PIV systems, Figure 1. They work together to measure two consecutive velocity fields, during the same engine cycle, which allows a temporal tracking of the flow structures.

The in-cylinder flow was seeded with oil particles (Rhodorsil 47V100). A generator (LaVision Aerosol Generator) produced particles of approximately 1µm at a concentration of more than 10 8 particles per cm 3 . These particles were then illuminated using a light sheet with a thickness of less than 1 mm. The light was emitted from two double pulsed Nd-YAG lasers (Specta Physics PIV 200) at a wavelength of 532 nm. The maximum output energy was around 200 mJ per pulse, with a pulse duration in the range of 8-10 ns. These lasers provided a stable cycle-tocycle illumination of the cylinder at a 10 ± 2 Hz frequency and their beams are combined by a beamsplitter optic.

After going through two Nikon AF Nikkor 50 mm f/1.8D, the scattered light was captured on two monochrome CCD cameras (Jai CV-M2 CL) with a resolution of 1600×1200 pixels. The spatial resolution was about 68µm/pixel. The first camera is equipped with an LCD shutter. This avoids the shots of the second laser to illuminate its second frame, Figure 2. It has an opening time of 1600µs and takes 600µs to close, so it must be triggered before the first laser shot. The time required to close it implies that the minimum measurable angular difference, for two consecutive fields at 1200rpm, is 4.3°CA. The synchronization of these devices was performed using a controller (Dantec Synchronizer) that took the signal of the engine encoder as input. The measurement started after the rotation speed of the engine was stable. It then took 300 pairs of images corresponding to 300 engine cycles at the two desired crank angles. An adaptive crosscorrelation PIV algorithm from Dantec DynamicStudio software is used to process the pairs of images. The interrogation area went from 32×32 to 64×64 pixels according to the particle density and a grid step size of 16×16 pixels is used. For each instantaneous velocity field, 77×78 instantaneous velocity vectors spaced 1mm apart were finally provided.

PIV measurements

Measurements were performed at four crank angles. The first point chosen was 146°CA, when the intake valves are open. The second one was 200°CA, just after BDC and the third one 270°CA at half of the compression stroke. These three points enable us to follow the evolution of the structured tumble movement. The last point is at 320°CA and corresponds to the maximum angle observable by the PIV system. To connect each of these four angles,

( 4 2 ) = 4! 2!2!
= 6 dual PIV measurements were performed, Figure 3. The temporal evolution of the flow structures can then be followed between each of these angles. Each angle is observed three times, so the measurements give 3×300 = 900 instantaneous velocity fields. This large number of velocity fields will allow defining a criterion on the center of rotation to obtain a high rotation rate, part 3.2. Figure 4 presents the mean flow evolution for the four crank angles measured. It is calculated as the average of 900 instantaneous fields and to improve readability, one vector out of three is displayed. Throughout this paper, the intake valves are on the right side of the figures.

The time-resolved measurements, made possible by the dual PIV, will then allow classifying the velocity fields by similarity throughout the compression phase, part 3.3. A classification will first be performed by a k-means algorithm on the instantaneous velocity fields of the first angle (146°CA). The temporal link between the structures at the different angles will then transmit the information of the classification to the three other angles (200, 270 and 320°CA). Part 3.3.1 details the approach. 3 CLUSTERING ANALYSIS

K-means algorithm

The classification of the data is done here using a k-means algorithm. This algorithm was developed in the '60s, it was introduced by Lloyd in 1957 [START_REF] Lloyd | Least Squares Quantization in PCM[END_REF] who was working on pulse code modulation and used the term "Least Square Quantization". Later, MacQueen was the first to use the term "kmeans" in 1967 [START_REF] Macqueen | Some Methods for Classification and Analysis of Multivariate Observations[END_REF]. This method consists in partitioning the data into K groups, called clusters, by trying to minimize the distance of each element to the mean of the elements of the cluster associated with it [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]).

The algorithm starts by initializing K average elements representative of the K clusters in an optimal way [START_REF] Arthur | K-Means++: The Advantages of Careful Seeding[END_REF]. These average elements are called centroids. The centroid term can refer to a center (physical coordinates) but also to a velocity field or other quantities depending on the data to be processed. This initialization step is followed by an iterative process composed of two steps. The first step classifies each element in its cluster, by calculating the distance between this element and the centroid of each cluster and then associating it to the closest cluster. The second step consists in recalculating the centroids. The process ends when the elements do not change clusters anymore and the centroids do not vary anymore. Figure 5 shows the flowchart. In this study, the Matlab function "kmeans" is used ("K-Means Clustering -MATLAB Kmeans -MathWorks France," n.d.). It has an algorithm to quickly initialize the average centers of the clusters. The use of this pre-programmed function makes the classification very fast and the best classification out of 1000, using different initial values of the cluster centroids, is retained.

The distance must be adapted to the data to be processed. In this study, we will treat a set of rotation centers 𝑥 ⃗, represented by their position in space, and a set of velocity fields 𝑢 ⃗⃗, represented by their components in the measurement plane.

The first case will thus use the euclidean distance which is classically used to quantify the proximity of points in space:

𝑑 1,𝑘 (𝑥 ⃗, 𝑐 𝑘 ⃗⃗⃗⃗) = (𝑥 ⃗ -𝑐 𝑘 ⃗⃗⃗⃗)(𝑥 ⃗ -𝑐 𝑘 ⃗⃗⃗⃗) 𝑡
Where 𝑐 𝑘 ⃗⃗⃗⃗ (𝑘 ∈ [1: 𝐾]) are the centroids and 𝑌 ⃗⃗ 𝑡 is the transpose vector of 𝑌 ⃗⃗ .

𝑢 ⃗⃗ = (𝑈 1 , … , 𝑈 𝑋 𝑚 , 𝑉 1 , … , 𝑉 𝑋 𝑚 )
The distance used will then be a correlation. It will be done on centered and normalized data, in order to study the whole range of flow variations. It is calculated with the following formula:

𝑑 2,𝑘 (𝑢 ⃗⃗, 𝑐 𝑘 ⃗⃗⃗⃗) = 1 - (𝑢 ⃗⃗ -𝑢 ̅ ⃗⃗ )(𝑐 𝑘 ⃗⃗⃗⃗ -𝑐 𝑘 ̅ ⃗⃗⃗⃗ ) 𝑡 √ (𝑢 ⃗⃗ -𝑢 ̅ ⃗⃗ )(𝑢 ⃗⃗ -𝑢 ̅ ⃗⃗ ) 𝑡 √ (𝑐 𝑘 ⃗⃗⃗⃗ -𝑐 𝑘 ̅ ⃗⃗⃗⃗ )(𝑐 𝑘 ⃗⃗⃗⃗ -𝑐 𝑘 ̅ ⃗⃗⃗⃗ ) 𝑡 Where 𝑢 ̅ ⃗⃗ = 1 2𝑋 𝑚 ∑ 𝑢 𝑗 2𝑋 𝑚 𝑗=1
. 1 2𝑋 𝑚 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ and 1 2𝑋 𝑚 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ is a row vector of 2𝑋 𝑚 ones. The vector 𝑐 ⃗ 𝑘 (𝑘 ∈ [1: 𝐾]) still represents the centroids and is written in the same manner as 𝑢 ⃗⃗.

Criterion on the center of rotation

This part seeks to give a criterion on the instantaneous center of rotation of the tumble to obtain high rotation rates after the BDC.

It uses the 900 instantaneous velocity fields at 200°CA from which 900 centers of rotation are extracted. An algorithm based on two detection criteria is used to detect the instantaneous centers of rotation [START_REF] Graftieaux | Combining PIV, POD and Vortex Identification Algorithms for the Study of Unsteady Turbulent Swirling Flows[END_REF][START_REF] Chaillou | Développement d'outils de Description et de Caractérisation de l'aérodynamique: Application Aux Écoulements Internes Issus Des Cylindres Des Moteurs Automobiles[END_REF]). The first one computes at any point P of the search area S the following criterion:

Γ 1 (𝑃, 𝑆) = 1 𝑁 𝑆 ∑ 𝑠𝑖𝑛(𝜃 𝑀 𝑖 ) 𝑁 𝑠 𝑖=1
Where the points 𝑀 𝑖 are in the search area 𝑆 and 𝑁 𝑆 is the number of points in 𝑆. During rotational motion, the angle 𝜃 𝑀 𝑖 between 𝑃𝑀 𝑖 ⃗⃗⃗⃗⃗⃗⃗⃗ (𝑃 ∈ 𝑆) and 𝑢 𝑀 𝑖 ⃗⃗⃗⃗⃗⃗⃗ is tending toward 90°. It is then necessary that 𝛤 1 tends to 1. Some areas of the flow may resemble vortex structures but be unclosed. They may be recirculation or shear zones. The following criterion can be used to differentiate them from centers of rotation:

Γ 2 (𝑃, 𝑆) = 1 𝑁 𝑆 √ ∑(𝑠𝑖𝑛(𝜃 𝑀 ) ̅̅̅̅̅̅̅̅̅̅ -𝑠𝑖𝑛(𝜃 𝑀 𝑖 )) 2 𝑁 𝑆 𝑖=1
This second criterion quantifies the deviation between the angle of each grid point and the average angle, too large a deviation indicates that the rotational motion is not "complete" or closed, 𝛤 2 will then need to be minimized.

The intensity of the rotational movements can then be calculated. We chose to use the criterion 𝑅 𝑡 = 𝐿/𝜔𝐼, which is the dimensionless ratio of the angular momentum 𝐿 [𝑘𝑔. 𝑚 2 . 𝑠 -1 ] to the moment of inertia 𝐼 [𝑘𝑔. 𝑚 2 ] of the motion divided by the engine speed 𝜔 [𝑟𝑎𝑑. 𝑠 -1 ]. It can be written in the form:

𝑅 𝑡 = ∑ (𝑥 𝑚 -𝑥 𝑐 )𝑉 -(𝑦 𝑚 -𝑦 𝑐 )𝑈 𝑁 𝑚 𝑚=1 2𝜋𝑁 60 ∑ (𝑥 𝑚 -𝑥 𝑐 ) 2 + (𝑦 𝑚 -𝑦 𝑐 ) 2 𝑁 𝑚 𝑚=1
Where 𝑁 𝑚 is the number of points (𝑥 𝑚 ; 𝑦 𝑚 ) in the measurement plane, (𝑥 𝑐 ; 𝑦 𝑐 ) are the coordinates of the center of rotation, 𝑈 and 𝑉 are the velocity components and 𝑁 = 1200 rpm is the engine rotation speed.

Figure 6 shows the position of the instantaneous centers of rotation on the mean field at 200°CA. The intensity of their rotation is also indicated by a colorbar. It seems that the points located on the left side have a higher intensity. It would therefore be interesting to give a criterion on the position of the center of rotation allowing to obtain a more intense movement, but also to quantify the gain in rotational intensity obtained by placing the center at this position. The question then arises of how to partition this set of points into different parts of the cylinder. This is where clustering comes in. It will optimally distribute the cloud of points in space into different zones. We have chosen to divide it into three clusters, which will represent high, medium and low rotational movements.

Figure 7 presents the result of the clustering process. The upper part shows the distribution of all the rotation centers in the 3 clusters and the lower part presents the average intensity of the rotation associated with each cluster. As indicated in the legend of the upper part, each cluster has about 300 points, so the distribution is homogeneous for a total of 900 points. We can observe that the points of cluster 1, in green, have an average rotation intensity that is 7.8% higher than that of cluster 3, in blue, with the lowest intensity used as a reference. Similarly, points in cluster 2, in red, have an average 𝑅 𝑡 that is 17% higher than those of cluster 3. It is therefore very interesting to place the centers in this area.

Figure 7 : Distribution of centers in 3 clusters of high, medium and low rotational motion intensity

Figure 8 shows the average field associated with each cluster. We can observe that the cluster with the strongest intensity is the one whose centers of rotation are farthest from the high speed zone. In the same way, the cluster with the weakest intensity has centers of rotation close to the high speed region. The momentum generated is therefore weaker, resulting in weaker rotation intensity. The return area is also more intense for the red cluster. This can be explained by the fact that the centers are further away from the right side of the cylinder. The upward motion of the tumble is thus more intense to compensate for the motion generated by the rear jet of the intake valves located on the upper right side of the cylinder.

There is then an area of the cylinder where placing the center of rotation of the tumble motion increases its intensity by about 17%. This was demonstrated close to the BDC, when the valve jet no longer alters the dynamics of the motion. The question now arises as to whether high intensity motions preserve this property during the compression phase.

Evolution of the tumbling flow

Clustering method

In order to study the evolution of the tumble flow, the velocity fields are classified by similarity into different groups. These families of fields will then represent flows of different intensities. We have therefore chosen to classify the fields into three groups representing strong, medium and weak flows. Note that the clusters are defined by similarity calculated using the distance 𝑑 2,𝑘 between a velocity field 𝑢 ⃗⃗ and the cluster 𝑘 represented by its centroid 𝑐 𝑘 ⃗⃗⃗⃗. Also note that the initial centroids are optimally chosen by the algorithm (part 3.1).

The clustering on the six velocity fields couples starts with the application of the k-means algorithm on the first three points at 146°CA. It is thus applied on a large dataset composed of 900 observations and to remove similarity between the data, the mean part of the set of three measurements is subtracted. The algorithm is therefore run on the turbulent field.

After running the algorithm, all 146°CA velocity fields are classified. The second angle of the three pairs 146/200, 146/270 and 146/320 are then directly classified by the temporal link provided by the dual PIV measurements. Figure 9 then shows how to associate each field of the other measured angles with these three groups.

The 200°CA velocity fields of the 146/200 pair are used to calculate the 200°CA centroids. Note that the centroid 𝑐 𝑘 ⃗⃗⃗⃗ of the cluster 𝑘, at a precise crank angle, is the mean value of its elements, at this crank angle. These centroids are then used to classify the 200/320 and 200/270 pairs. To do this, we compute the distance between each 200°CA velocity fields, of the two pairs, and the 200°CA centroids. The Matlab function "pdist2" is used with the correlation distance 𝑑 2 . A velocity field is associated with the group to which it is the closest. All 200°CA velocity fields are then classified and the second angles of these two pairs are directly classified by the temporal link provided by the dual PIV measurements.

The 270°CA velocity fields, of the 146/270 pair, are used to calculate the 270°CA centroids in order to classify the two angles of the last 270/320 pair, following the same approach as for the fields at 200°CA. All velocity fields are now classified into three clusters. These clusters were defined at 146°CA by the k-means algorithm, which classified the fields by similarity in an optimal way. The temporal link between the measurements was then used to pass the cluster information to the following angles.

Clustering visualization

Clustering can be visualized using a method of Proper Orthogonal Decomposition (POD) of velocity fields. This method developed by [START_REF] Lumley | The Structure of Inhomogeneous Turbulent Flows[END_REF][START_REF] Lumley | The Structure of Inhomogeneous Turbulent Flows[END_REF][START_REF] Chen | On the Use and Interpretation of Proper Orthogonal Decomposition of In-Cylinder Engine Flows[END_REF][START_REF] Fogleman | Application of the Proper Orthogonal Decomposition to Datasets of Internal Combustion Engine Flows[END_REF][START_REF] Vu | Proper Orthogonal Decomposition Analysis for Cycle-to-Cycle Variations of Engine Flow. Effect of a Control Device in an Inlet Pipe[END_REF]) allows decomposing the velocity vector into 𝑁 𝑘 components, classified by decreasing order of energy, in a new space. The POD is interesting here because the first components are the most energetic and thus contain a large part of the information. It therefore makes it possible to visualize the distribution of the clusters on these first components during compression. This method will be done here on the turbulent part of the velocity fields. The mean part, common for all fields of the same angle, is thus subtracted from the instantaneous velocity fields. The turbulent velocity, at cycle 𝑛, is written as follows:

𝑢 ′ ⃗⃗⃗⃗ 𝑛 (𝑥 𝑚 , 𝑦 𝑚 , 𝜃) = ∑ 𝑎 𝑛 (𝑘) Φ ⃗⃗⃗⃗ (𝑘) (𝑥 𝑚 , 𝑦 𝑚 , 𝜃) 𝑁 𝑘

𝑘=1

Where 𝑎 𝑛 (𝑘) is the component number 𝑘 of the turbulent velocity field at cycle 𝑛 in the basis composed of the Φ ⃗⃗⃗⃗ (𝑘) vectors.

Energy in clusters

The classification of the flow into different groups having been achieved, it is then possible to follow the evolution of the energy quantities representative of each cluster. It is chosen to follow the mean kinetic energy (𝑀𝐾𝐸), the turbulent kinetic energy (𝑇𝐾𝐸) and the rotation rate (𝑅 𝑡 ). To this end, the instantaneous velocity vector 𝑈 𝑛 (or 𝑉 𝑛 ), of cycle 𝑛 in cluster 𝑗, is divided into an average part

𝑛=1

The rotation rate formula is described in part 3.2.

Figure 11 presents the evolution of the flow energies quantities representative of each cluster, according to the crank angle. Note that at 320°CA the rotation rate has not been calculated. Indeed, no rotational motion is observable, so the calculation of 𝑅 𝑡 is not possible and has no physical meaning. It can be observed that the hierarchy between clusters remains the same over all the measurement points. The red curve, associated with cluster 3, has a level of mean kinetic energy and rotation rate that is the highest for all observed angles. Cluster 3 is therefore representative of the fields creating the most intense tumble motions. It is followed by cluster 1, green curve, representative of the fields of average intensity and then cluster 2, blue curve, associated with the least intense movements. This observation is important because it demonstrates that the most intense movements after BDC remain the most intense during the compression phase. Taking a closer look at the values at 320°CA, Figure 12, we can see that cluster 3 has the highest total kinetic energy value and MKE. Its total kinetic energy is 18.3% higher than cluster 1, so this cluster has the highest potential to create the highest turbulence rate at the TDC. Then, although the total kinetic energy levels are equal for clusters 1 and 2, we can see that cluster 1 has a higher MKE level, this MKE which will later be transformed into TKE, will allow this cluster to have the second highest turbulence level at TDC. Furthermore, the curves in Figure 11 indicate that the MKE level drops between 270 and 320°CA. This suggests that the large-scale motions turn into turbulence. This observation is verified by the stagnation of the TKE curves between 270 and 320°CA. Although the turbulence dissipates over time, the TKE level remains stable thanks to a transfer of energy from the large-scale motions. The MKE therefore tends to transform into TKE and the hierarchy between the groups should be maintained at the TDC. The intense movements before compression should thus generate the highest turbulence rates for the combustion phase of the engine.

CONCLUSION

This paper presented a clustering analysis method of the internal aerodynamics generated in a Miller cycle gasoline engine. The tumble motion, classically used to create the level of turbulence for efficient combustion, can be perturbed by the premature closure of the intake valves. The design approach for a new engine may therefore be different from that of an Otto cycle engine. In general, the performance of a cylinder head is compared by observing curves of rotational intensity versus valve lift. It is then necessary to ensure that an intense movement, before compression, gives a higher turbulence rate for the combustion phase.

The flow study was performed on a motored single-cylinder transparent engine following a Miller lifting law of 140°CA. It was associated with a dual PIV system that measures two velocity fields per engine cycle. The evolution of the flow structures, between two crankshaft angles, could thus be followed. The dual PIV systems also make it possible to record a large number of cycles per test. They are therefore well suited to statistical studies, unlike other time-resolved PIV systems, such as HS-PIV or TR-PIV.

The clustering study, using a k-means algorithm, was performed in two steps:

 First, a cloud of 900 instantaneous tumble rotation centers was partitioned into three clusters at 200°CA. This showed that there is an area in the cylinder space that achieves approximately 17% higher rotation rates. This observation may guide the development of future Miller cycle engines.

An effort could be made to better position the center of rotation of the movement in this zone by working, for example, on the inclination of the intake duct or the shape of the combustion chamber roof. One could also imagine flow control systems that would position the tumble movement in the right zone according to, for example, the offset generated by a lift law variation system. Indeed, most of the new engines have a Variable Valve Timing (VVT) system that angularly shifts the lift law of the intake valves. The maximum opening of the valves is therefore never at the same position in relation to the position of the piston, which generates movements of different intensities and centers (Perceau, Guibert, and Guilain 2021a).

The second part classified the velocity fields by resemblance to each other. It has been seen that these three groups represent strong, medium and weak flows. All the data of the dual PIV experimental design were linked in order to follow the characteristic quantities of each group (total kinetic energy, turbulent kinetic energy and rotation rate). It has shown that:  Working on the creation of intense movements before the compression phase can generate a total kinetic energy that is 18.3% higher close to the end of the compression stroke, which can lead to higher turbulence rates close to the combustion top dead center.

Future work could measure the velocity fields at angles closer to the TDC combustion in order to measure the difference in TKE between the different groups. A CFD numerical approach could also be adopted. The gains in terms of thermal efficiency of the engine could then be estimated using, for example, complete 0D models.

Figure 1 :

 1 Figure 1 : Dual PIV measuring system

Figure 2 :

 2 Figure 2 : Dual PIV system chronogram

Figure 3 :

 3 Figure 3 : Dual PIV measurements

Figure 4 :

 4 Figure 4 : Mean flow evolution -From top to bottom: 146, 200, 270 and 320°CA

Figure 5 :

 5 Figure 5 : K-means algorithm flowchart

Figure 6 :

 6 Figure 6 : Position and intensity of the 900 instantaneous centers of rotation at 200°CA.

Figure 8 :

 8 Figure 8 : Average rotational movement in each cluster -200°CA

Figure 9 :

 9 Figure 9 : Classification of the velocity fields couples

Figure 10

 10 Figure 10 presents the evolution of the clusters according to the first three components, 𝑎 𝑛 (1) , 𝑎 𝑛 (2) , and 𝑎 𝑛 (3) of the velocity vectors, in the new space using a POD. A clear distinction between the different groups can be observed for the first two angles 146 and 200°CA. The distinction then remains visible at 270°CA and becomes almost imperceptible at 320°CA. We can then say that the similarity information between 3 types of structures first moves in groups without deconstructing the initial clusters. There is then a memory effect of these clusters throughout 124°CA (between 146 and 270°CA), which can be related to the notion of the integral time scale of the flow. Note that the POD is performed on the turbulent part of the velocity fields at each crank angle. The loss of coherence at 320°CA then shows that the information about the different clusters shifts to the more distant components. The compression stroke leads to a movement of the clusters' information to the more distant components. The flow starts to lose its coherence at 270°CA due to cyclic fluctuations or turbulence within each cluster, shear, out-of-plane movements etc. At 320°CA the distinction between clusters is totally lost, which indicates that they can no longer be differentiated by the first 3 modes. The next part then focuses on the energy quantities of the total velocity field, including all components. It is then possible to define a specific behavior for each cluster.

Figure 10 :

 10 Figure 10 : Representation of the clusters according to the first three components of the POD

Figure 11 :

 11 Figure 11 : Evolution of the flow energies quantities representative of each cluster

Figure 12 :

 12 Figure 12 : Total kinetic energy at 320°CA of each cluster

  𝑗 (𝜃, 𝑚) 2 + 𝑉 ̅ 𝑗 (𝜃, 𝑚) 2 )

	𝑀𝐾𝐸(𝜃, 𝑗) =	1 2𝑁 𝑚	∑ (𝑈 ̅ 𝑁 𝑚 𝑚=1
	Where 𝑁 𝑚 is the number of points in the
	measurement plane. The turbulent kinetic energy at
	crank angle 𝜃 for cluster 𝑗 is then calculated as the
	mean value of the TKE value of cycles inside cluster
	𝑗:
	𝑇𝐾𝐸(𝜃, 𝑗) =	2𝑁 𝑗 𝑁 𝑚 1	𝑚=1 ∑ ∑ (𝑢 𝑛 ′ (𝜃, 𝑚) 2 + 𝑣 𝑛 ′ (𝜃, 𝑚) 2 ) 𝑁 𝑚 𝑁 𝑗
					𝑈 ̅ 𝑗 =	1 𝑁 𝑗	∑ 𝑁 𝑗 𝑛=1	𝑈 𝑛	(or 𝑉 ̅ 𝑗 ) over all 𝑁 𝑗 cycles of cluster 𝑗
					and a turbulent part 𝑢 𝑛 ′ (or 𝑣 𝑛 ′ ) . The mean kinetic
					energy at crank angle 𝜃 for cluster 𝑗 is calculated as
					follows:

Acknowledgements

This work was supported by Renault SA, Guyancourt, France.

DECLARATIONS

Conflict of interest

The authors declare that they have no conflict of interest.