Flow Field Parametric Interpolation Using a Proper Orthogonal Decomposition: Application to the Variable Valve Timing Effect on a Tumble In-cylinder Miller Engine Mean Flow

: The current article presents a method to reconstruct the mean velocity ﬁeld of a cyclic ﬂow for an input parameter value that has not been measured, allowing for the number of tests to be reduced. It is applied to the tumble ﬂow of a gasoline engine following a Miller cycle. New engines often include variable valve timing (VVT) systems to maximize the efﬁciency of such over-expanded cycles for different operating points. The reconstruction was thus carried out using different offset values of the intake valve lift timing. Experimental data were collected from a transparent engine in an early intake valve closing (EIVC) conﬁguration using particle image velocimetry (PIV). The mean velocity ﬁeld reconstruction was based on the interpolation of the proper orthogonal decomposition (POD) coefﬁcients. The accuracy of the method was evaluated at different points by comparing the interpolated and the measured ﬂow ﬁelds. The accuracy was estimated by calculating the error in the rotation rate of the tumble and the position of its center of rotation. The new mean velocity ﬁeld set allowed for the position of the tumble’s center of rotation to be closely tracked according to the input parameter and a rotation rate map to be made. Some results on Miller’s cycle could thus be found and the data generated could guide future developments.


Introduction
This study focused on the tumble flow inside an internal combustion spark-ignition engine. During the intake phase of an air-gasoline mixture into an engine, an aerodynamic rotational movement is created. This large-scale motion allows for the storage of kinetic energy that will be transformed into smaller-scale motions during the compression phase. These small-scale turbulent motions allow for wrinkling the flame front during the combustion phase and thus increase the flame propagation speed in order to obtain a good thermal efficiency. Several studies were undertaken to improve its stability and intensity in the classical Otto cycle configuration [1][2][3][4][5][6][7]. However, this cycle is now starting to be left behind. Indeed, the environmental constraints of the last few years have seen the application of over-expanded cycles, such as the Atkinson and Miller cycles, which further increase the efficiency of thermal engines [8][9][10][11][12][13]. Car manufacturers are beginning to apply them to their new engines [14][15][16][17][18]. One drawback is that these cycles affect the intake valves by increasing or reducing their opening time, which strongly disturbs the tumble movement. This study focused on the tumble flow generated by a Miller lift law in the EIVC (early inlet valve closing) configuration. Most of the new engines also include a system that allows for varying the valve lift laws [19][20][21] so that the engine operation can then be optimized for different loads. A system that enables the angular shift of the intake valve lift law was studied here, which is called VVT (variable valve timing), and it was on the angular offset of this system that the study was done.
To characterize the tumble motion is to obtain information on its position in the cylinder and on its intensity and thus to obtain the velocities composing the motion. The velocities can be calculated numerically or measured experimentally. The study presented in this paper was experimental. It was carried out on a transparent cylinder engine. It allowed for comparisons to be made with results that were obtained numerically and thus to check the correct operation of codes that model the flow of new Miller engines. Two main methods are used to experimentally obtain velocities on an engine, namely, LDV (laser doppler velocimetry) [22] and PIV (particle image velocimetry) [23]. PIV was chosen for its greater spatial resolution, which provides the complete velocity field inside the cylinder.
The VVT angular offset is a parameter that modifies the aerodynamics and thus changes the characteristics of the tumble motion. A high number of points must be measured to be able to characterize the tumble for different offset values and thus optimally control the new Miller engine. This is why we present in this paper a method to interpolate a velocity field according to a parameter modifying the internal aerodynamics. By knowing the velocity fields that are obtained for two values of the parameter, it is possible to reconstitute the average field generated by a value of the parameter included in the starting range. This method reduces the number of tests to be carried out during the study of the mean motion of a cyclic flow.
Parametric interpolation is made possible by the proper orthogonal decomposition (POD) method, which is generally used to study the cyclic stability of engine flows [7,[24][25][26][27][28][29]. The parametric interpolation principle is presented very briefly in [30] and is mainly used for model reduction and command control. Some analysis of data obtained using timeresolved particle image velocimetry (TR-PIV) presented the temporal interpolation of velocity fields [31,32]. Indeed, temporal tracking, which is made possible by the high acquisition frequency, allows for a good reconstruction of the fields during a cycle. Parametric interpolation is, however, limited by this type of measurement, which can only retain a limited number of cycles; the convergence of average speeds is not sufficient in this case. The interpolation can also be done between 2D velocity fields to reconstruct a 3D field [33]. In this experimental study, 2D PIV data that was measured from a transparent engine was used. The flow was observed several times at a given instant of the cycle for different VVT offsets. The 300 observations allowed for obtaining a representative mean velocity field of the flow and thus created usable reconstructed fields after parametric interpolation.
The paper is organized as follows. In Section 2, the experimental setup of the measured data is outlined. Section 3 then presents the flow field parametric interpolation method. The effect of the Miller cycle on the tumble mean flow characteristics is then described in Section 4. The article ends with several conclusions.

Model Engine
The model engine was representative of a 0.9 L 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, the 1200 rpm rotation was provided by an electrical engine. 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.2 mm stroke and a 72.2 mm bore. A flat piston head was used; 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 PIV experiments were performed in the tumble plane, which ran through the axis of the cylinder and the middle of the two intake and exhaust valves ( Figure 1). Throughout this paper, the intake valves are on the right side of the figures. The VVT offset adjustment was made with the grey part in Figure 2. This part locked the camshaft before installing the timing belt and the angle of its rectangular openwork set the VVT offset. As the camshaft made one rotation for two rotations of the crankshaft, the angle of the rectangular openwork was half the VVT offset in crank angle degrees.

PIV Measuring System
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 a double pulsed Nd-YAG laser (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. This laser provided a stable cycle-to-cycle illumination of the cylinder at a 10 ± 2 Hz frequency, Figure 4. After going through a Nikon AF Nikkor 50 mm f/1.8D, the scattered light was captured on a monochrome CCD camera (Jai CV-M2 CL) with a resolution of 1600 × 1200 pixels. The spatial resolution was about 68 µm/pixel. 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 desired crank angle. An adaptive cross-correlation 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, 71 × 59 instantaneous velocity vectors were finally provided with a spatial resolution of 1.08 mm.

PIV Measurements
As mentioned in the introduction, the studied parameter that modified the aerodynamics was the VVT offset, corresponding to an angular shift of the intake valve timing. The measurements were taken at several crank angles indicated in Table 1. The first point chosen was at 130 • CA because the tumble motion is established for any VVT offset at this angle. The last one was at 270 • CA and corresponded to the middle of the compression stroke. As previously mentioned, these experiments were performed at 1200 rpm and atmospheric pressure.
The results of a convergence check for the 300 cycles are displayed in Figure 5. It shows the absolute mean value evolution of the velocity components according to the cycle number, where the mean horizontal component U i = 1 i ∑ i n=1 U n was calculated with the instantaneous components U n over i cycles. Calculations were made in four different locations (x ; y) = (±Bore/4; 0) and (0 ; ±Stroke/4) using the θ = 200 • CA ; θ VVT = 0 • CA flow field. The convergence was estimated with the error e i,U = It was below 1 and 2% for the U and V components, respectively, after 200 cycles (i > 200). This meant that the value of the average velocity field almost did not vary when the number of cycles was higher than 200. The mean flow could, therefore, be accurately calculated from 300 cycles and could represent the global behavior of the flow.

Proper Orthogonal Decomposition
The interpolation of the mean flow field according to a parameter is made possible by the POD method thanks to two of its properties. First, it allows for defining the velocity vector using a new basis whose coordinates do not depend on the position but only on the time and the parameters influencing the fluid movements. Moreover, this new orthogonal basis is an optimal decomposition in the sense of kinetic energy, i.e., the most energetic portion of the flow, i.e., the mean field, is found in the first few components. The mean field can therefore be easily extracted by determining the number of components needed.
The velocity field at a given time in cycle n and for a parameter θ VVT is classically written as follows: where X m represents the position in space and U n (or V n ) is the component along the x-(or y-) axis along the elemental vector → e x (or → e y ). After changing the basis using the POD method, the vector is written as follows: This new formulation gives rise to N k components a (k) n that are associated with the elementary vectors → Φ (k) . The components have a precise order, the total kinetic energy they contain decreases with k and we refer to them as modes to designate the value of k ( Figure 6). We can notice that here the position information is included in the elementary vectors and not in the components. As expressed in the introduction of this part, the mean field could be recomposed by interpolating the components along the parameter θ VVT , which represents the VVT offset, and selecting the appropriate number of modes. This decomposition was first applied to turbulent flows in 1967 by Lumley [34]. It consists of looking for a function φ that follows the flow as much as possible, which tends to be parallel to the velocity vector. Therefore, the projection of the velocity vector onto φ should be maximized, with the maximum being reached when the two vectors are parallel. It will also be necessary to remove the influence of the amplitude of φ on the test of parallelism of the two vectors, which is why the norm of the desired function appears in the problem formulation. Mathematically, the problem is written as follows: is the scalar product, ||.|| is the associated norm and . is the mean operator.
The method for solving this problem uses a velocity vector correlation tensor whose eigenvalues are calculated. The associated eigenvectors are used to define the components a (k) n of the basis, which are then used to compute the elemental vectors → Φ (k) . A matrix formatting of the data allows this method to be applied quite easily. The methodology is clearly described by Vu [35] and the reader can refer to it for a detailed formulation. In this study, the POD was applied to a single angle θ and for four values of the parameter θ VVT = {0, 10, 20, 30}.

Number of Required Modes
The interpolation must allow for the reconstruction of the mean flow field; therefore, it must be ensured that it has a sufficient number of modes to capture all the information. Figure 7 shows the evolution of the POD coefficients as a function of the VVT offset. The coefficients of the individual cycles are shown in gray and the average coefficient over all the cycles a (k) is shown in blue. This figure gives a first idea of the number of modes that need to be selected to represent the average field. Indeed, the more the instantaneous coefficient of a cycle moves away from the average, the more it represents the behavior of the fluid that is specific to the cycle in question and not an average behavior over all the cycles. It can be observed that the instantaneous coefficients seemed to have the same tendencies and thus represented the average field up to mode 4. The average coefficient of modes 5 and 6 was, however, non-zero and continued to vary with θ VVT , as well as for mode 8. It was then necessary to go further in the analysis in order to allow for the choice of the cutoff mode. To ensure that enough modes were selected, it was possible to compare the mean field → u = ∑ N c n=1 → u n over N c cycles of an angle and VVT offset of the POD with its reconstruction → u K from K modes. The reconstructed field is written as follows: The comparison was made at the 130 • CA angle, to remain in the same conditions as Figure 7 by looking at the absolute value of the subtraction of the reconstructed → u K from the mean flow fields → u . Figure 8 presents the results for the different VVT offsets on flow fields that were reconstructed with 1, 2, 6 and 8 modes. It shows that the results with six and eight modes were very similar.  Figure 9 shows the total energy error between the mean and reconstructed flows as a function of the number of modes. The classical formula e = E mean f ield −E reconstructed f ield, k modes E mean f ield × 100 was used. It is presented at 130 (middle of the intake stroke), 220 (beginning of compression stroke) and 270 • CA (middle of the compression stroke) for all VVT offsets. It can be seen that the total energy error was less than 1% from K = 6. For the rest of the study, we set the number of modes to six to reconstruct the mean field. The same approach was applied to the other angles of the experimental study for any VVT shift and the value of six modes was also sufficient.

Accuracy of the Reconstruction
With the observations made in the previous section, we chose to take six modes to reconstruct the mean flow field. Figure 10 shows the average field at θ = 130, 220, 270 • CA and θ VVT = 30 • CA and the associated reconstruction. For all fields presented in this paper, the inlet valves are on the right (Figure 1), and for improved readability, one vector out of two is displayed. The direction of the vectors seems to be the same for both fields from the same crank angle presented. However, deviations on the norm of the velocity vectors were observed for the 130 • CA flow fields, mainly in the high-velocity areas. Two quantities could be calculated to quantify the accuracy of the reconstruction: the position of the tumble motion center of rotation (x c ; y c ) and the intensity of its rotation.
An algorithm based on two detection criteria was used to detect the instantaneous centers of rotation [36,37]. The first one computes the following criterion at any point P of the search area S: where the points M i are in the search area and N S is the number of points in S (Figure 11).
During rotational motion, the angle between → PM i and → u M i tends 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: Figure 11. Center detection scheme.
This second criterion quantifies the deviation between the angle of each grid point and the average angle: too large of a deviation indicates that the rotational motion is not "complete" or closed; Γ 2 will then need to be minimized. Figure 12 shows the mean velocity fields of the experimental plan crank angles for a VVT offset of 0 • CA with their centers of rotation detected using the algorithm. One out of two vectors is displayed to improve the readability.
To quantify the intensity of rotation, we chose to use the criterion R t = L/ωI, which is the dimensionless ratio of the angular momentum L to the moment of inertia I of the motion divided by the engine speed ω. It is referred to as the rotation rate throughout this paper. This quantity is particularly suitable for the comparison of VVT systems. Indeed, the VVT offset generates a variation in the air load of the engine and thus the rotating air mass. Assuming a homogeneous distribution of the mass in the chamber, the division of L by I eliminates it from the equation. It can therefore be written in the form: where N m is the number of points in the measurement plane (4189 for θ = 180 • CA) and N = 1200 rpm is the engine rotation speed.  Table 2 shows the deviations in the position of the centers of rotation and the rotation rates. The positional deviations of the center of rotation are expressed as the number of mesh points and those of the rotation rate in percent. We chose to do the calculations for three crank angles of the experimental plan. This shows that the deviations in the position of the center of rotation were rather small. They were on average less than one mesh point (two mesh points were about 1 mm apart). The maximum mean deviation in the rotation rate was about 2.3%.

Interpolation
As noted above, the first six modes of the POD could be used to reconstruct the mean velocity field of the tumble motion with good accuracy. The objective of interpolation is to generate the mean field for a VVT offset value outside the POD database decomposition for values different than {0, 10, 20, 30}. In order to validate the method, additional measurements were taken for θ VVT = 25, which was a value outside the POD decomposition. As a reminder, the new POD basis was calculated at each crank angle of the experimental plan (Table 1) for θ VVT = {0, 10, 20, 30}. Supplementary fields were measured for θ VVT = 25 • CA at three angles θ = {130, 220, 270}. Therefore, it was possible to compare the average field measured for this value of VVT to the interpolated field using the POD.
Three interpolation methods were compared: linear, polynomial (pchip in Matlab) and spline (piecewise-defined function with polynomials) interpolation. Figure 13 shows the evolution of the average coefficients of the POD that were interpolated by following these three methods. As in the previous section, the center position and rotation rate of the tumble were used to compare mean and interpolated fields. Table 3 presents the results for the interpolation at θ VVT = 25 • CA. We observed that the linear interpolation gave the best results. Differences in the center position are low with a maximum deviation of 2 points for the y value of the 270 • CA flow field. The mean error in the tumble rotation rate was 4.16%. Therefore, the interpolation provided a good estimation of the mean field. Table 3. Comparison between mean and interpolated fields of points outside the POD.

Linear
Cubic Spline  Figure 14 presents the comparison between mean and interpolated fields for a linear interpolation at 130, 220 and 270 • CA (from top to bottom). As shown in Table 3, the results at 130 • CA are very similar, where only a difference in the general shape of the rolling motion can be observed. At 220 • CA, the tumble motions seemed to be very similar. However, the interpolated flow field presented a lower intensity next to the right wall, which explained the 5% difference in the rotation rate shown in Table 3. The same observation was made for the last crank angle flow field. The increase in the difference with the angle value could be explained by an increase in the cyclic variability of the flow during the compression phase of the tumble motion. By making the motion less intense, the strong VVT offset could thus destabilize the flow during the compression, making the mean field less representative of the flow. This could be the subject of a future study. In the rest of this paper, we focus on the characterization of the mean tumble motion for all the VVT shifts studied.

Tumble Mean Flow Characteristics
It is well known that the tumble motion in gasoline engines is a complex threedimensional movement. It also has an instantaneous characteristic, in other words, it is likely to vary from one engine cycle to another. However, it is common to refer to an ensemble average field in order to describe its global behavior and its evolution during the intake and compression strokes. It can also be considered that the tumble movement is globally symmetrical with respect to the "tumble plane" (Figure 1). Two-dimensional measurements in this area can give a good idea of the intensity of the rotational movement. By selecting the first modes, the POD method can extract the mean flow field from instantaneous measurements. Interpolation of the POD coefficients then allows for reconstructing the mean flow for unmeasured VVT offsets; the previous section gives the accuracy of such a reconstruction on the studied model engine. Any average flow field of a measured crank angle with a VVT offset in the measurement interval can therefore be known. This section focuses on the effect of the VVT on the tumble mean motion, more specifically on its intensity and the position of its center of rotation. The data presented were calculated from velocity fields reconstructed with six modes using linear interpolation.

Mapping of the Rotation Rate
The mean field was interpolated for any θ VVT belonging to [0,1,2,...,30]. The center detection algorithm was then applied to each field at a fixed crank angle and for any VVT offset so that each rotation rate could be calculated with respect to its center of rotation. Figure 15 shows the mapping of the results. It can be observed on this map that the rotation rate was low after the BDC. This was due to the Miller intake valves' lift law, which had an opening time of 140 • CA. This short duration of the intake valves opening made it difficult for the tumble movement to be as strong as for an Otto cycle with the same engine configuration. This is why the values tended to be less than 1.
Comparisons were made with the same engine configuration [38] following an Otto cycle ( Figure 16). Note that the Miller intake valve lift had a VVT offset of 0 • CA. The rotation movement of the classical law (Otto cycle) was present at 100 • CA, while with the Miller law, the rotation was only visible from 110 • CA; this is due to the early opening of the intake valves of the Otto law. Between 120 and 146 • CA, the rotation rates were very close and then moved away until they differed by about 15% at 180 • CA. The Miller law did not allow for maintaining the motion as long as the Otto law; it then lost intensity. The loss of rotation intensity was about 20% at 270 • CA, which was not negligible.
A second observation can be made from Figure 15. Before the BDC (180 • CA), the VVT offset is profitable for the rotation rate. Indeed, the valves opened earlier (Figure 3). The tumble movement thus had more time to form. An optimum could also be seen around the 20 • CA VVT offset. This trend did not continue after the BDC. Indeed, after 200 • CA, the smallest VVT offsets produced the highest rotation rates. Even though the movement was more intense for high VVT offsets at the beginning of the intake stroke, it struggled to persist due to the early closing of the intake valves.   As mentioned before, the intake valves are on the right. The valve jet followed the left wall of the cylinder, then the piston, and finally moved up along the right wall to form a rolling motion. In the middle of the intake phase, the strong VVT offsets were more intense. The rotational movements were driven by the valve jet, which pushed them to the opposite side, to the right of the cylinder. During the compression phase, the center distribution followed the same pattern. Low VVT offsets tended to be centered and away from the piston, while high offsets tended to be closer to the walls. The phenomenon became more pronounced for the angles 220 and 270 • CA. The early closure of the intake valves made the movements created by the strong VVT offsets less intense. This loss of energy decentered the movements that approached the walls and thus became even less intense.

Conclusions
This paper describes a method to interpolate the mean velocity field at a given time according to the value of a parameter that modified the aerodynamics of a cyclic flow. This technique was applied to the tumble motion inside the cylinder of a Miller cycle gasoline engine. The parameter on which the interpolation was performed was the VVT offset and the data were collected experimentally in a transparent engine using PIV measurements.
The method used a proper orthogonal decomposition, which was chosen to reconstruct the average field with the first six modes of the POD. The maximum mean value error on the rotation rate between the mean and the reconstructed fields of the POD was 2.3%. The error of the position of the center of rotation was below 1 mm. The quality of the interpolation was then estimated using the measured fields for VVT offset values outside the POD.

•
The parametric interpolation error on the tumble rotation rate was about 5% and the errors on the position of the center of rotation were below 1 mm.
This method allowed for reducing the number of points that were measured experimentally to extract information from the mean velocity field. Knowing the R t at different angles for any VVT offset would help to guide ignition strategies. A low R t led to slower combustion and thus to a higher risk of knocking. The optimal ignition, which avoids knocking while having the best possible efficiency, could be evaluated. The optimization could then be done at the same time with a VVT offset and injection timing. In the same way, the knowledge of the position of the center of the tumble linked to that of the R t for any VVT would make it possible to optimize fuel injection. A better homogenization of the load could be obtained by choosing the optimal tuning of injection according to the position of the injectors, the type of injection (direct or indirect), the position of the tumble and its intensity. Knowing the R t for a large VVT offset interval would also allow for the VVT system to be predimensioned. When the offset creates a low R t , leading to a risk of knocking that cannot be avoided by changing the ignition timing, the use of VVT to regulate the load should be discontinued and intake duct throttling should be used. In this way, the maximum offset of the VVT system can be estimated.

•
A mapping that shows the tumble rotation rate for any VVT offset of the studied range and at different angles was given. The evolution of the tumble center position with the VVT offset at different angles was presented.
The offset in the timing of the intake valves in the direction of a valve crossing led to a decrease in the rotational movement. The movement was no longer sustained by the valve jet before the piston reached the bottom dead center. In addition, this VVT offset led to a decrease in the flow rate and, therefore, in the mass of air in rotation; the movement then had less inertia and had difficulty conserving its energy after the intake valves were closed. This loss of energy was associated with a decentering of the movement in the cylinder and the center got closer to the walls, which led to an even greater decrease in the intensity of the rotation. It was then necessary to adapt the engine to the Miller cycle by optimizing, for example, the intake duct and the shape of the combustion chamber roof to limit the loss of rotation intensity.

•
The Miller valve timing generated a less intense rotational motion than a conventional one. The decrease in rotation rate intensity was about 20% at the mid-compression for the same engine configuration.
Parametric interpolation was applied to the VVT offset but could have been done on other parameters related to the internal aerodynamics, such as the masking of the valves or the inclination of the intake duct. It could also be done on the external aerodynamics of a vehicle by investigating, for example, the inclination of a spoiler. More generally, it could be applied to any system using blades. Consequently, this technique can have many applications.