Based on the results of the previously discussed experimental activities, a modeling approach is employed here consisting of the integration of a complete 1D fluid-dynamic model of the whole engine and of properly developed correlations capable of quantifying the cycle-to-cycle variation of Indicated mean effective pressure (IMEP) and In-cylinder pressure peak (*pmax*). Although the adopted 1D modeling can be considered a simple numerical method, it shows, on one hand, the advantage of a reduced computational effort and, on the other hand, it is enhanced with advanced 0D sub-models to refine the description of in-cylinder processes also including the cyclic dispersion (CCV). The CCV sub-model is capable of reproducing on a physical basis the measured sequence of consecutive pressure traces in each cylinder. In a first step, 1D engine model was validated against the experimental data considering both the overall performance of the engine and the ensemble-average in-cylinder pressure cycles. The CoV correlations of IMEP and *pmax* with average combustion parameters were developed by numerical fitting of the outcomes deriving from the post-processing of experimental spark-sweep tests. Then the correlations were validated with the measured data set already adopted for the 1D model validation (Table 3). Once validated, the combination of the 1D engine model and of the experimental CoV correlations allows the provision of a combustion optimization of the individual engine cylinder. In addition, this numerical methodology provides benefits in terms of engine fuel economy and cyclic dispersion level, as well. A detailed analysis of the employed numerical models, of their potentialities and limitations, are reported in the following sub-sections.

#### *3.1. D Engine Modeling*

The geometry of the analyzed engine is fully schematized in GT-Power 1D commercial code, through the combination of 0D and 1D modeling approaches. In particular, the intake and exhaust sub-systems described through a 1D approach to properly reproduce the flow and the pressure waves propagation. The typical flow equations (Navier–Stokes) are solved in each pipe of intake and exhaust manifolds. This means that the solutions of continuity, momentum, enthalpy and energy conservation equations are provided. The conservation equations solved by GT-Power are shown below:

Continuity:

$$\frac{dm}{dt} = \sum\_{\text{boundary}} \dot{m}\_{\text{\textdegree}} \tag{1}$$

Energy:

$$\frac{d(me)}{dt} = -p\frac{dV}{dt} + \sum\_{\text{boundaries}} \left(\dot{m}H\right) - hA\_s \left(T\_{f\text{fluid}} - T\_{\text{wall}}\right) \tag{2}$$

Enthalpy:

$$\frac{d(\rho HV)}{dt} = \sum\_{b \text{undaries}} \left( \dot{m}H \right) + V \frac{dp}{dt} - hA\_s \left( T\_{fluid} - T\_{wall} \right) \tag{3}$$

Momentum:

$$\frac{d\dot{m}}{dt} = \frac{A d\mu + \sum\_{\text{boundaries}} \left(\dot{m}u\right) - 4\mathcal{C}\_f \frac{\rho u|u| A dx}{2D} - \mathcal{K}\_p A \left(\frac{1}{2}\rho u|u|\right)}{d\mathbf{x}}\tag{4}$$

where *m˙* , *m*, *V*, *p*, ρ, *e*, H and u are mass flux, mass, volume, pressure, density, specific internal energy, enthalpy and velocity at boundary, respectively. *A* and *As* are the flow area and the heat transfer surface area, *Tfluid* and *Twall* the fluid and wall temperatures, *h* the heat transfer coefficient and *D* the equivalent diameter. *Cf* and *Kp* are the friction and pressure loss coefficients, while *dx* and *dp* represent the discretization length and the pressure differential across *dx*.

The in-cylinder processes, including fuel evaporation, air/fuel mixture formation, turbulent combustion, knock and heat transfer are reproduced based on a 0D modeling. Concerning the port fuel injection, it is actuated by a typical "injector object" imposing the air-fuel ratio.

A value equal to 30% of the injected gasoline is assumed to instantaneously evaporate after the injection event, as results from the advised setting of the adopted 1D code. Furthermore, PFI injection is modeled by neglecting the liquid wall film generation along the intake pipes and the spray dynamics. The cylinder head flow permeability is described through the definition of the steady flow coefficients measured at the engine test bench, both in forward and reverse flow conditions. Turbocharger operation is characterized by means of the measured performance maps of compressor and turbine components. The mechanical friction losses of the engine are reproduced by an empirical correlation according to the Chen-Flynn model. Friction correlation includes the variability of the engine speed and on the in-cylinder peak pressure. This correlation is validated through the available experimental data, demonstrating a high accuracy level. The following friction correlation for FMEP (friction mean effective pressure) is employed here:

$$\text{FMEP} = 0.876 + 2.66 \times 10^{-3} \cdot p\_{\text{max}, \text{cyl}} - 6.9 \times 10^{-2} \cdot v\_{\text{pm}} + 8.1 \times 10^{-3} \cdot v\_{\text{pm}}^2 \tag{5}$$

where *pmax,cyl* is the cylinder peak pressure and *vpm* the mean piston speed.

As detailed above, the in-cylinder phenomena are simulated by advanced sub-models, implemented into the code under user routines. The combustion event is reproduced by a 0D two zone (burned and unburned gas) sub-model, where the burning rate term is computed through the well-known 'fractal' approach. According to the fractal model, the burning rate term is written as:

$$\frac{dm\_b}{dt} = \rho\_{\rm u} A\_T S\_L = \rho\_{\rm u} A\_L \sum S\_{L\prime} \tag{6}$$

where ρ*u* is the unburned gas density, *AL* and *AT* the area of the laminar and turbulent flame fronts, respectively, *SL* the laminar flame speed (LFS) and Σ the wrinkling factor.

It was demonstrated to present a high level of accuracy in describing the turbulent flame front propagation inside the conventional combustion chamber of spark–ignition engines [25]. As emerged from the fractal theory available in the scientific literature [26], some relevant turbulence parameters are included into the wrinkling factor (Σ) of the fractal-derived burn rate expression. This requires a proper evaluation of the above mentioned in-cylinder turbulence variables by means of a phenomenological turbulence sub-model. To this aim, the combustion model is coupled to a turbulence one. The latter represents a K-k-T model, which solves three balance equations to furnish the engine cycle evolution of mean flow kinetic energy (K), turbulent kinetic energy (k) and tumble vortex momentum (T) [27]. In previous authors' works [27,28], the adopted turbulence model was proved to adequately estimate the in-cylinder turbulence parameters at varying the engine speeds and the valve strategies, without a model tuning variation. A kinetically-derived LFS correlation, included within the fractal burning rate, is adopted because of its better prediction capabilities with respect to the experimentally-derived LFS formulations [29]. LFS correlation is obtained by fitting 1D LFS computations via a chemical kinetic solver and it accounts for the in-cylinder thermodynamic state, equivalence ratio and charge dilution.

Finally, AL is estimated by an automatic procedure implemented into CAD software processing the actual 3D geometry of the engine combustion chamber.

Referring to the cycle-to-cycle variation model, starting from the input data of CoV of the in-cylinder pressure peak (CoVpmax), representative *faster-than-average* and *slower-than-average* cycles are computed basing on a proper algorithm. This algorithm, extensively discussed in [30], derives *faster* and *slower* burning profiles by the average fractal-predicted burn profile through two stretch factors. The obtained burnt profiles (high, average and low) correspond to three pressure traces. For each in-cylinder pressure cycle, the related maximum value (*pmax*) is computed. The developed CCV model automatically modifies the stretch factors until the target levels of pmax for high and low cycles are realized. In this way, the measured CoVpmax imposed as input data in the model is numerically reproduced. It is worth emphasizing that the adopted CCV model was also demonstrated to be very efficient in terms of computational effort [30].

Knock occurrence is detected in the model by the auto-ignition (AI) calculation of the unburned air/fuel mixture. AI is computed by using a tabulated approach [31]. The AI table is obtained by chemical kinetics simulations of auto-ignition, performed in a homogeneous reactor at constant pressure (CP) conditions. In particular, AI table is generated through chemical kinetics simulations carried out by considering the kinetic scheme of Andrae [32]. The AI time is stored in the corresponding table as a value depending on the in-cylinder thermodynamic condition (pressure and temperature), mixture quality (equivalence ratio), and residual content (internal EGR). Knock occurrence is recognized when the AI Integral exceeds a prescribed threshold level [31]. Concerning the heat transfer modeling, the in-cylinder heat transfer (gas-to-wall) is described by the modified Hohenberg correlation [33], while convective, conductive and radiative heat transfer models are considered for the exhaust pipes to improve the prediction of temperature at Turbine Inlet. A standard finite element (FEM) approach, implemented into the 1D code, allows the computation of both piston, cylinder liner and head temperatures by assigning reasonable initial temperature levels. The Hohenberg correlation includes the dependency on pressure and volume of gas and also accounts for the cylinder geometry (bore) and mean piston speed. The overall convection multiplier of the Hohenberg model was selected in order to reach a good experimental/numerical agreement for the in-cylinder pressure cycles at different speed/load points, especially in the expansion phase. *Energies* **2020**, *13*, x FOR PEER REVIEW 9 of 22 both piston, cylinder liner and head temperatures by assigning reasonable initial temperature levels. The Hohenberg correlation includes the dependency on pressure and volume of gas and also accounts for the cylinder geometry (bore) and mean piston speed. The overall convection multiplier of the Hohenberg model was selected in order to reach a good experimental/numerical agreement for the in-cylinder pressure cycles at different speed/load points, especially in the expansion phase. The convective heat transfer coefficients for engine coolants (water and oil) are estimated by

The convective heat transfer coefficients for engine coolants (water and oil) are estimated by proper simulations of coolant circuits and subsequently implemented into the 1D code, assuming a dependency on the engine speed according to an assigned power law [34]. Cooling boundary conditions (temperatures for oil and water) are imposed in the model according to the levels recommended by the engine manufacturer. proper simulations of coolant circuits and subsequently implemented into the 1D code, assuming a dependency on the engine speed according to an assigned power law [34]. Cooling boundary conditions (temperatures for oil and water) are imposed in the model according to the levels recommended by the engine manufacturer. The 0D quasi-dimensional combustion model was tuned by a trial-and-error procedure to

The 0D quasi-dimensional combustion model was tuned by a trial-and-error procedure to achieve the best agreement with the in-cylinder pressure cycles. It should be stated that a unique set of tuning constants was selected for all the tested operating conditions. achieve the best agreement with the in-cylinder pressure cycles. It should be stated that a unique set of tuning constants was selected for all the tested operating conditions.

#### *3.2. Model Validation 3.2. Model Validation*

The developed 1D model is utilized to simulate the measured operating points on the engine domain, as reported in Table 3. Regarding the setup of the model, the experimental values of combustion phasing (MFB50), air/fuel (A/F) ratio and intake valve strategy (inlet valve closing, IVC) are imposed in the 1D model simulations. In addition, the measured IMEP level is matched by a proportional integral derivative (PID) controller acting on the throttle valve opening at low loads or on the waste-gate valve opening at high loads. EGR valve was not regulated during the simulation and it was fixed to a substantially closed setting since the experimental points to be reproduced were investigated at a null external EGR rate. The model reliability is proved here in terms of overall engine performance variables, combustion characteristics, in-cylinder pressure traces and cycle-to-cycle variation. The numerical/experimental assessments of global performance parameters such as Air flow rate, ISFC, and In-cylinder pressure peak for two cylinders are reported in Figures 2 and 3. The developed 1D model is utilized to simulate the measured operating points on the engine domain, as reported in Table 3. Regarding the setup of the model, the experimental values of combustion phasing (MFB50), air/fuel (A/F) ratio and intake valve strategy (inlet valve closing, IVC) are imposed in the 1D model simulations. In addition, the measured IMEP level is matched by a proportional integral derivative (PID) controller acting on the throttle valve opening at low loads or on the waste-gate valve opening at high loads. EGR valve was not regulated during the simulation and it was fixed to a substantially closed setting since the experimental points to be reproduced were investigated at a null external EGR rate. The model reliability is proved here in terms of overall engine performance variables, combustion characteristics, in-cylinder pressure traces and cycle-to-cycle variation. The numerical/experimental assessments of global performance parameters such as Air flow rate, ISFC, and In-cylinder pressure peak for two cylinders are reported in Figures 2 and 3.

**Figure 2.** Experimental/Numerical comparison of Air Flow rate (**a**) and indicated specific fuel consumption (ISFC) (**b**) at various engine operating points. **Figure 2.** Experimental/Numerical comparison of Air Flow rate (**a**) and indicated specific fuel consumption (ISFC) (**b**) at various engine operating points.

**Figure 3.** Experimental/Numerical comparison of In-cylinder pressure peak of Cyl#1 (**a**) and Cyl #2

**2000@5**

**2000@7**

**2000@9**

**2000@11**

**2000@13**

**3000@5**

**3000@7**

**3000@12**

**3000@16**

**3000@18**

**4000@7**

**4000@9**

**4000@11**

**4000@13**

**4000@16**

(**b**) at various engine operating points.

**2000@5**

**In-cylinder pressure peak, bar**

**2000@7**

**2000@9**

**2000@11**

**2000@13**

**3000@5**

**3000@7**

**3000@12**

**3000@16**

**3000@18**

**4000@7**

**4000@9**

**4000@11**

**4000@13**

**4000@16**

**2000@72000@92000@112000@133000@53000@73000@123000@163000@184000@74000@94000@114000@134000@162000@52000@72000@92000@112000@133000@53000@73000@123000@163000@184000@74000@94000@114000@134000@16**

**ISFC, g/kWh**

*Energies* **2020**, *13*, x FOR PEER REVIEW 9 of 22

both piston, cylinder liner and head temperatures by assigning reasonable initial temperature levels. The Hohenberg correlation includes the dependency on pressure and volume of gas and also accounts for the cylinder geometry (bore) and mean piston speed. The overall convection multiplier of the Hohenberg model was selected in order to reach a good experimental/numerical agreement for the in-cylinder pressure cycles at different speed/load points, especially in the expansion phase. The convective heat transfer coefficients for engine coolants (water and oil) are estimated by proper simulations of coolant circuits and subsequently implemented into the 1D code, assuming a dependency on the engine speed according to an assigned power law [34]. Cooling boundary conditions (temperatures for oil and water) are imposed in the model according to the levels

The 0D quasi-dimensional combustion model was tuned by a trial-and-error procedure to achieve the best agreement with the in-cylinder pressure cycles. It should be stated that a unique set

The developed 1D model is utilized to simulate the measured operating points on the engine domain, as reported in Table 3. Regarding the setup of the model, the experimental values of combustion phasing (MFB50), air/fuel (A/F) ratio and intake valve strategy (inlet valve closing, IVC) are imposed in the 1D model simulations. In addition, the measured IMEP level is matched by a proportional integral derivative (PID) controller acting on the throttle valve opening at low loads or on the waste-gate valve opening at high loads. EGR valve was not regulated during the simulation and it was fixed to a substantially closed setting since the experimental points to be reproduced were investigated at a null external EGR rate. The model reliability is proved here in terms of overall engine performance variables, combustion characteristics, in-cylinder pressure traces and cycle-to-cycle variation. The numerical/experimental assessments of global performance parameters such as Air flow rate, ISFC, and In-cylinder pressure peak for two cylinders are reported in Figures 2 and 3.

*Energies* **2020**, *13*, 5548 consumption (ISFC) (**b**) at various engine operating points.

**Air Flow rate, kg/h**

**2000@5**

recommended by the engine manufacturer.

*3.2. Model Validation*

of tuning constants was selected for all the tested operating conditions.

**Figure 2.** Experimental/Numerical comparison of Air Flow rate (**a**) and indicated specific fuel

**Figure 3.** Experimental/Numerical comparison of In-cylinder pressure peak of Cyl#1 (**a**) and Cyl #2 (**b**) at various engine operating points. **Figure 3.** Experimental/Numerical comparison of In-cylinder pressure peak of Cyl#1 (**a**) and Cyl #2 (**b**) at various engine operating points.

The air flow rate (Figure 2a) is predicted with an acceptable accuracy (average absolute percent error of 1.9%), denoting an adequate schematization of the engine geometry. A good experimental/numerical correlation is also realized for the ISFC (average absolute percent error of 2.4% in Figure 2b), hence demonstrating the model capability to properly take into account the combined effects of flow, combustion, and in-cylinder heat transfer phenomena.

The numerical in-cylinder pressure peaks of two cylinders are in a satisfactory agreement with the experimental counterparts. Referring to the in-cylinder pressure peak predictions, larger errors arise for Cyl #1 compared to Cyl #2. However, the average absolute percent error for the Cyl #1 case is below the maximum allowable percent level of 5%, which is considered by the authors to preserve a valuable model accuracy.

The reliability of the combustion modeling is demonstrated by the outcomes reported in Figure 4a,b that shows the experimental/numerical assessments of the combustion core duration (MFB 10–50%) for Cyl #1 (Figure 4a) and Cyl #2 (Figure 4b). The combustion sub-model correctly reproduces the combustion development for the selected engine operating conditions (speed/load point, A/F ratio), mainly including the effects of IMEP variation on combustion evolution. It is should be highlighted that a reduced error in the reproduction of the combustion durations has to be considered a fundamental pre-requisite to ensure a proper application of the CCV correlations, as discussed in the next section. The predictivity of the combustion model is also proved by the experimental/numerical assessments of in-cylinder pressure traces for two cylinders shown in Figures 5–7. They refer to three different engine speeds and loads, labelled as 2000@13, 3000@18 and 4000@7, without the activation of the external EGR circuit (i.e., e-EGR% = 0%). In the analyzed operating points, some disagreements arise along the expansion phase of pressure cycles for Cyl #1, probably due to minor inaccuracies in forecasting the combustion tail and/or the heat transfer. However, from a global point of view, numerical in-cylinder pressure cycles present a good agreement with the experimental traces. A similar accuracy is detected for the other investigated speed/load points (Table 3), even if the related plots are not shown here for brevity. In addition, the CCV model is applied to the experimental operating points. To this aim, the experimentally derived coefficient of variation of the in-cylinder peak pressure (CoVpmax) for two cylinders is imposed in the simulations and the computed cycles (high and low ones), representative of the cycle-to-cycle variation, and these are compared to the measured sequence of 270 consecutive pressure cycles. As an example, Figure 8 shows the experimental/numerical assessments for pressure cycles in a medium/low load point and medium speed, i.e., 3000@7. The figure demonstrates that the numerical high and low cycles are in a satisfactory agreement with the extreme *faster-than-average* and *slower-than-average* experimental cycles, respectively. The consistency of the CCV modeling approach is also tested in other speed/load points listed in Table 3. Although not reported here for brevity, in many cases the shape of numerical high and low cycles is close to the measured ones. Of course, inaccuracies in the computation of the ensemble-average in-cylinder pressure trace may reflect CCV prediction, too. Once validated, the CCV model can be used in a predictive way to reproduce the cycle-to-cycle variation of operating points different from those considered in the validation phase. To this aim, the CCV model requires the CoV equation capable of quantifying the variability on the in-cylinder

a valuable model accuracy.

pressure peak. The latter correlation together with the one of IMEP variability will be extensively discussed in terms of generation and validation in the next paragraph. of the external EGR circuit (i.e., e-EGR% = 0%). In the analyzed operating points, some disagreements arise along the expansion phase of pressure cycles for Cyl #1, probably due to minor inaccuracies in forecasting the combustion tail and/or the heat transfer. However, from a global point of view,

different engine speeds and loads, labelled as 2000@13, 3000@18 and 4000@7, without the activation

*Energies* **2020**, *13*, x FOR PEER REVIEW 10 of 22

combined effects of flow, combustion, and in-cylinder heat transfer phenomena.

The air flow rate (Figure 2a) is predicted with an acceptable accuracy (average absolute percent error of 1.9%), denoting an adequate schematization of the engine geometry. A good experimental/numerical correlation is also realized for the ISFC (average absolute percent error of 2.4% in Figure 2b), hence demonstrating the model capability to properly take into account the

The numerical in-cylinder pressure peaks of two cylinders are in a satisfactory agreement with the experimental counterparts. Referring to the in-cylinder pressure peak predictions, larger errors arise for Cyl #1 compared to Cyl #2. However, the average absolute percent error for the Cyl #1 case is below the maximum allowable percent level of 5%, which is considered by the authors to preserve

The reliability of the combustion modeling is demonstrated by the outcomes reported in Figure 4a,b that shows the experimental/numerical assessments of the combustion core duration (MFB 10–50%) for Cyl #1 (Figure 4a) and Cyl #2 (Figure 4b). The combustion sub-model correctly reproduces the combustion development for the selected engine operating conditions (speed/load point, A/F ratio), mainly including the effects of IMEP variation on combustion evolution. It is should be highlighted that a reduced error in the reproduction of the combustion durations has to be considered a fundamental pre-requisite to ensure a proper application of the CCV correlations, as discussed in the next section. The predictivity of the combustion model is also proved by the experimental/numerical

#### *3.3. Experimentally-Derived Correlations for Cyclic Dispersion and Validation* numerical in-cylinder pressure cycles present a good agreement with the experimental traces. A

In the present section, the experimental CCV is extensively analyzed with main aim to define its relationship with the representative engine combustion variables. For this purpose, it is worth underlining that CCV is usually quantified through the Coefficient of Variation (CoV) of the pressure-related parameters, such as the indicated mean effective pressure (IMEP) and the in-cylinder pressure peak (*pmax*). In particular, the first one plays a certain role in the drivability of the engine, while the second one has an influence on the knock occurrence. Based on the available experimental results reported in Table 2, the cycle-to-cycle variation was preliminarily analyzed for each engine cylinder at varying spark timing. As an example, the following Figure 9a,b and Figure 10a,b show the spark timing influence on the dispersions of IMEP and *pmax*, considering operations at low load and at speed of 3000 rpm. In these figures, the average values of IMEP and *pmax* and the related CoVs are reported as continuous lines while the corresponding cycle-derived levels as vertical dots. similar accuracy is detected for the other investigated speed/load points (Table 3), even if the related plots are not shown here for brevity. In addition, the CCV model is applied to the experimental operating points. To this aim, the experimentally derived coefficient of variation of the in-cylinder peak pressure (CoVpmax) for two cylinders is imposed in the simulations and the computed cycles (high and low ones), representative of the cycle-to-cycle variation, and these are compared to the measured sequence of 270 consecutive pressure cycles. As an example, Figure 8 shows the experimental/numerical assessments for pressure cycles in a medium/low load point and medium speed, i.e., 3000@7. The figure demonstrates that the numerical high and low cycles are in a satisfactory agreement with the extreme *faster-than-average* and *slower-than-average* experimental cycles, respectively. The consistency of the CCV modeling approach is also tested in other speed/load points listed in Table 3. Although not reported here for brevity, in many cases the shape of numerical high and low cycles is close to the measured ones. Of course, inaccuracies in the computation of the

As expected, an increasing cycle-to-cycle variation is observed moving the spark timing towards the expansion phase, due to the worsening in the combustion development. For each cylinder, a greater *pmax* variation is realized with respect to the IMEP variation. A comparison between the outcomes of two cylinders (Figures 9 and 10) highlights that Cyl #2 (rich A/F mixture) shows lower CoV levels than Cyl#1 (lean A/F mixture) when varying the spark timing. Figures 9 and 10 also confirm the experimental evidence of higher IMEP and pmax values attained by the richer Cyl #2 if compared to those realized by Cyl #1. ensemble-average in-cylinder pressure trace may reflect CCV prediction, too. Once validated, the CCV model can be used in a predictive way to reproduce the cycle-to-cycle variation of operating points different from those considered in the validation phase. To this aim, the CCV model requires the CoV equation capable of quantifying the variability on the in-cylinder pressure peak. The latter correlation together with the one of IMEP variability will be extensively discussed in terms of generation and validation in the next paragraph.

**Figure 4.** Experimental/Numerical comparison of combustion core duration mass fraction burned (MFB) 10–50% for Cyl#1 (**a**) and Cyl #2 (**b**) at various engine operating points. (MFB) 10–50% for Cyl#1 (**a**) and Cyl #2 (**b**) at various engine operating points.

**Figure 5.** Experimental/Numerical comparison of in-cylinder pressure traces for Cyl#1 (**a**) and Cyl #2 (**b**) at 2000@13. **Figure 5.** Experimental/Numerical comparison of in-cylinder pressure traces for Cyl#1 (**a**) and Cyl #2 (**b**) at 2000@13.

**Figure 6.** Experimental/Numerical comparison of in-cylinder pressure traces for Cyl#1 (**a**) and Cyl#2

**Figure 7.** Experimental/Numerical comparison of in-cylinder pressure traces for Cyl#1 (**a**) and Cyl#2

(**b**) at 3000@18.

(**b**) at 4000@7.

*Energies* **2020**, *13*, 5548 (**b**) at 2000@13.

(**b**) at 2000@13.

**Figure 5.** Experimental/Numerical comparison of in-cylinder pressure traces for Cyl#1 (**a**) and Cyl #2

**Figure 5.** Experimental/Numerical comparison of in-cylinder pressure traces for Cyl#1 (**a**) and Cyl #2

*Energies* **2020**, *13*, x FOR PEER REVIEW 11 of 22

**Figure 4.** Experimental/Numerical comparison of combustion core duration mass fraction burned

*Energies* **2020**, *13*, x FOR PEER REVIEW 11 of 22

**Figure 4.** Experimental/Numerical comparison of combustion core duration mass fraction burned

(MFB) 10–50% for Cyl#1 (**a**) and Cyl #2 (**b**) at various engine operating points.

(MFB) 10–50% for Cyl#1 (**a**) and Cyl #2 (**b**) at various engine operating points.

**Figure 6.** Experimental/Numerical comparison of in-cylinder pressure traces for Cyl#1 (**a**) and Cyl#2 (**b**) at 3000@18. **Figure 6.** Experimental/Numerical comparison of in-cylinder pressure traces for Cyl#1 (**a**) and Cyl#2 (**b**) at 3000@18. (**b**) at 3000@18.

**Figure 7.** Experimental/Numerical comparison of in-cylinder pressure traces for Cyl#1 (**a**) and Cyl#2 (**b**) at 4000@7. **Figure 7.** Experimental/Numerical comparison of in-cylinder pressure traces for Cyl#1 (**a**) and Cyl#2 (**b**) at 4000@7. *Energies* **2020**, *13*, x FOR PEER REVIEW 12 of 22

**Figure 7.** Experimental/Numerical comparison of in-cylinder pressure traces for Cyl#1 (**a**) and Cyl#2

**Figure 8.** Experimental/Numerical comparison of cycle-to-cycle variation for Cyl#1 (**a**) and Cyl #2 (**b**) at 3000@7. **Figure 8.** Experimental/Numerical comparison of cycle-to-cycle variation for Cyl#1 (**a**) and Cyl #2 (**b**) at 3000@7.

*3.3. Experimentally-Derived Correlations for Cyclic Dispersion and Validation*  In the present section, the experimental CCV is extensively analyzed with main aim to define its relationship with the representative engine combustion variables. For this purpose, it is worth underlining that CCV is usually quantified through the Coefficient of Variation (CoV) of the pressurerelated parameters, such as the indicated mean effective pressure (IMEP) and the in-cylinder pressure peak (*pmax*). In particular, the first one plays a certain role in the drivability of the engine, while the second one has an influence on the knock occurrence. Based on the available experimental results reported in Table 2, the cycle-to-cycle variation was preliminarily analyzed for each engine cylinder Similar considerations can be applied to the other analyzed points of the spark timing sweep (Table 2), but the related experimental outcomes are not reported here for sake of brevity. In addition, as already demonstrated in a previous authors' work [35], the CoV trend does not substantially depend on the engine speed. On the other hand, the cycle-to-cycle variation shows a certain dependency on the combustion-related parameters, i.e., combustion phasing and characteristic duration. In order to optimize the engine's overall performance, the optimization of each cylinder operation was realized to eliminate the cylinder-to-cylinder non uniformities. In this perspective, the relationship between CCV characteristics and the combustion-related variables is analyzed here by averaging the experimental data of the two engine cylinders. In this way, experimental-derived CCV correlations can be defined

as continuous lines while the corresponding cycle-derived levels as vertical dots.

3000 rpm and low load.

at varying spark timing. As an example, the following Figure 9a,b and Figure 10a,b show the spark

**Figure 9.** IMEP and coefficient of variance (CoV)IMEP vs. Spark timing for Cyl#1 (**a**) and Cyl#2 (**b**) at

at 3000@7.

and referred to the examined engine. Once validated, engine CCV correlations can be easily employed to forecast the cyclic dispersion under optimized operating conditions for both cylinders. Based on the above considerations, the following Figures 11 and 12 show the dependency of the engine CoVIMEP and CoVpmax on the averaged combustion parameters, including 50 operating points (Table 2) and the data measured in all cylinders. reported in Table 2, the cycle-to-cycle variation was preliminarily analyzed for each engine cylinder at varying spark timing. As an example, the following Figure 9a,b and Figure 10a,b show the spark timing influence on the dispersions of IMEP and *pmax*, considering operations at low load and at speed of 3000 rpm. In these figures, the average values of IMEP and *pmax* and the related CoVs are reported as continuous lines while the corresponding cycle-derived levels as vertical dots.

second one has an influence on the knock occurrence. Based on the available experimental results

*Energies* **2020**, *13*, x FOR PEER REVIEW 12 of 22

**Figure 8.** Experimental/Numerical comparison of cycle-to-cycle variation for Cyl#1 (**a**) and Cyl #2 (**b**)

*Energies* **2020**, *13*, x FOR PEER REVIEW 13 of 22

In the present section, the experimental CCV is extensively analyzed with main aim to define its relationship with the representative engine combustion variables. For this purpose, it is worth underlining that CCV is usually quantified through the Coefficient of Variation (CoV) of the pressurerelated parameters, such as the indicated mean effective pressure (IMEP) and the in-cylinder pressure

*3.3. Experimentally-Derived Correlations for Cyclic Dispersion and Validation*

**Figure 9.** IMEP and coefficient of variance (CoV)IMEP vs. Spark timing for Cyl#1 (**a**) and Cyl#2 (**b**) at 3000 rpm and low load. **Figure 9.** IMEP and coefficient of variance (CoV)IMEP vs. Spark timing for Cyl#1 (**a**) and Cyl#2 (**b**) at 3000 rpm and low load. *Energies* **2020**, *13*, x FOR PEER REVIEW 13 of 22 Similar considerations can be applied to the other analyzed points of the spark timing sweep (Table 2), but the related experimental outcomes are not reported here for sake of brevity. In addition,

as already demonstrated in a previous authors' work [35], the CoV trend does not substantially

**Figure 10.** In-cylinder pressure peak and CoV pressure peak (*pmax*) vs. Spark timing for Cyl#1 (**a**) and Cyl#2 (**b**) at 3000 rpm and low load. **Figure 10.** In-cylinder pressure peak and CoV pressure peak (*pmax*) vs. Spark timing for Cyl#1 (**a**) and Cyl#2 (**b**) at 3000 rpm and low load. including 50 operating points (Table 2) and the data measured in all cylinders.

dependency on the combustion-related parameters, i.e., combustion phasing and characteristic **Figure 11.** CoVIMEP vs. averaged MFB50 (**a**) and MFB10–50 (**b**). **Figure 11.** CoVIMEP vs. averaged MFB50 (**a**) and MFB10–50 (**b**).

depend on the engine speed. On the other hand, the cycle-to-cycle variation shows a certain

duration. In order to optimize the engine's overall performance, the optimization of each cylinder operation was realized to eliminate the cylinder-to-cylinder non uniformities. In this perspective, the relationship between CCV characteristics and the combustion-related variables is analyzed here by averaging the experimental data of the two engine cylinders. In this way, experimental-derived CCV correlations can be defined and referred to the examined engine. Once validated, engine CCV correlations can be easily employed to forecast the cyclic dispersion under optimized operating conditions for both cylinders. Based on the above considerations, the following Figures 11 and 12 show the dependency of the engine CoVIMEP and CoVpmax on the averaged combustion parameters,

**Figure 11.** CoVIMEP vs. averaged MFB50 (**a**) and MFB10–50 (**b**).

**Figure 12.** CoVpmax vs. averaged MFB50 (**a**) and MFB10–50 (**b**). **Figure 12.** CoVpmax vs. averaged MFB50 (**a**) and MFB10–50 (**b**).

Referring to the CoVIMEP, higher values for delayed combustions are obtained (Figure 11a). A minimum level is also detected at a combustion phasing substantially corresponding to the one realizing the MBT condition. CoVIMEP also shows a clear increasing trend with the combustion core duration MFB10–50 (Figure 11b). Concerning the CoVpmax, it increases with a more delayed combustion process (Figure 12a), while a certain dispersed trend between *pmax* variation and combustion core duration is established in Figure 12b. However, a longer combustion phenomenon involves a greater CoVpmax. The discussed experimental data on CCV in Figures 11 and 12 were utilized to define proper correlations with averaged combustion parameters through the 'curve fitting' module available in the Matlab code. In particular, for the examined engine, two polynomial correlations were identified by linking the CoV levels with the averaged combustion phasing and core duration: Referring to the CoVIMEP, higher values for delayed combustions are obtained (Figure 11a). A minimum level is also detected at a combustion phasing substantially corresponding to the one realizing the MBT condition. CoVIMEP also shows a clear increasing trend with the combustion core duration MFB10–50 (Figure 11b). Concerning the CoVpmax, it increases with a more delayed combustion process (Figure 12a), while a certain dispersed trend between *pmax* variation and combustion core duration is established in Figure 12b. However, a longer combustion phenomenon involves a greater CoVpmax. The discussed experimental data on CCV in Figures 11 and 12 were utilized to define proper correlations with averaged combustion parameters through the 'curve fitting' module available in the Matlab code. In particular, for the examined engine, two polynomial correlations were identified by linking the CoV levels with the averaged combustion phasing and core duration:

$$\mathbf{^0CoV}\_{\mathbf{M}\mathbf{H}\mathbf{B}\mathbf{F}} = \mathbf{a}\mathbf{g}\_0 + \mathbf{a}\mathbf{a}\mathbf{M}\mathbf{B}\mathbf{B}\mathbf{g}\_0 + \mathbf{a}\mathbf{a}\mathbf{b}\mathbf{M}\mathbf{B}\mathbf{B}\mathbf{g}\_{0\mathbf{b}\mathbf{c}\mathbf{0}} + \mathbf{a}\mathbf{a}\mathbf{b}\mathbf{M}\mathbf{B}\mathbf{B}\mathbf{g}\_0^2 + \mathbf{a}\mathbf{a}\mathbf{M}\mathbf{B}\mathbf{H}\mathbf{B}\mathbf{g}\_0\mathbf{M}\mathbf{B}\mathbf{H}\mathbf{B}\mathbf{g}\_{0\mathbf{b}\mathbf{c}\mathbf{0}}\tag{7}$$

$$\mathbf{CoV}\_{\text{flmmec}} = \mathbf{tb}\_0 + \mathbf{tb}\_1 \mathbf{MBB} \mathbf{\hat{s}}\_{\text{f0}0} + \mathbf{b}\_2 \mathbf{b} \mathbf{MBB} \mathbf{\hat{s}}\_{\text{f0}\text{-}6\text{-}6} + \mathbf{b}\_3 \mathbf{b} \mathbf{MBB} \mathbf{\hat{s}}\_{\text{f0}\text{-}5} + \mathbf{b}\_4 \mathbf{b} \mathbf{BHB} \mathbf{\hat{s}}\_{\text{f0}\text{-}50} \tag{8}$$

where ai and bi represent the constants for the best fit with the experimental outcomes. The proposed correlations (7) and (8) show a good coefficient of determination (Table 4). where ai and bi represent the constants for the best fit with the experimental outcomes. The proposed correlations (7) and (8) show a good coefficient of determination (Table 4).

**Table 4.** CoV correlation parameters**. Table 4.** CoV correlation parameters.

The reliability of the experimental-based CoV correlations (7) and (8) was tested against the measured CoV values. The numerical/experimental comparisons for CoV of IMEP and *pmax* are reported in Figure 13: satisfactory agreements are observed for both correlations changing the speed/load points and at a null external EGR rate. A slightly better prediction is realized for the CoVIMEP. In agreement with the experimental literature [35], Figure 13 confirms that IMEP cyclic dispersion is lower than the corresponding in-cylinder pressure peak cyclic dispersion. The reliability of the experimental-based CoV correlations (7) and (8) was tested against the measured CoV values. The numerical/experimental comparisons for CoV of IMEP and *pmax* are reported in Figure 13: satisfactory agreements are observed for both correlations changing the speed/load points and at a null external EGR rate. A slightly better prediction is realized for the CoVIMEP. In agreement with the experimental literature [35], Figure 13 confirms that IMEP cyclic dispersion is lower than the corresponding in-cylinder pressure peak cyclic dispersion.

Once validated, the developed CoV correlations were linked to the 1D computed combustion parameters (average cycle) in order to provide a valuable prediction of the cycle-to-cycle variation in the operating conditions not explored during the experimental activity. To this aim, the adoption of refined combustion and turbulence sub-models is mandatory. As extensively discussed in the following section, the main use of the developed correlations will be the prediction of the CCV level in the numerically optimized engine points, where the individual cylinders show the same nominal operation.

**Figure 13.** Experimental/Numerical comparison of Coefficient of Variation (CoV) of IMEP (**a**) and Coefficient of variation of in-cylinder pressure peak (**b**) at various engine operating points. **Figure 13.** Experimental/Numerical comparison of Coefficient of Variation (CoV) of IMEP (**a**) and Coefficient of variation of in-cylinder pressure peak (**b**) at various engine operating points.

Once validated, the developed CoV correlations were linked to the 1D computed combustion

#### **4. Numerical Optimization of Combustion: Results and Discussion**

parameters (average cycle) in order to provide a valuable prediction of the cycle-to-cycle variation in the operating conditions not explored during the experimental activity. To this aim, the adoption of refined combustion and turbulence sub-models is mandatory. As extensively discussed in the following section, the main use of the developed correlations will be the prediction of the CCV level in the numerically optimized engine points, where the individual cylinders show the same nominal operation. **4. Numerical Optimization of Combustion: Results and Discussion**  As already reported in previous sections, the operating points measured in the entire engine domain (Table 3) were considered to perform a combustion optimization, assuming the same behavior for each cylinder. This means that the numerical optimization process takes into account the absence of non-uniformities among cylinders, which practically requires only partial modifications to the existing fuel injection system of the original engine. The numerical optimization is performed by employing the validated 1D model and the validated CoV correlations referring to the examined engine. It consists of a virtual engine re-calibration, mainly acting on the fuel injection and on the spark timing. In particular, at each measured point in Table 3, a stoichiometric A/F ratio is imposed in each cylinder by a refined control of injected gasoline mass. This condition eliminates the experimental cylinder-to-cylinder variation. At low loads, a PID controller acting on the throttle valve opening is used to match the prescribed IMEP while the waste-gate valve is fully opened. At high loads, a PID controller acts on the waste-gate valve opening to match the target IMEP, while the throttle valve is considered fully opened. The spark timing is modified with the aim to identify the MBT condition and in the case of knock occurrence it is automatically delayed to maintain the numerical knock index below a prescribed threshold level. In addition, at low engine speeds a safety margin is applied to the maximum boost pressure to avoid the compressor surge [36]. At high loads, both the in-cylinder pressure peak and the temperature at turbine inlet are kept below maximum allowable levels to preserve the engine safety. The numerical results obtained by the model-related optimization (labelled as 'Opt') are compared to the model prediction based on the standard ECU calibration (labelled as 'Base'). The above comparisons are realized with the aim to highlight the potential benefits arising from the optimal operation of the engine cylinders. These numericalderived advantages are here presented in terms of combustion phasing (MFB50) and core duration (MFB10–50), ISFC, and cyclic dispersion levels (CoVIMEP and CoVpmax). The following Figures 14–16 show the comparison of 'Base' and 'Opt' numerical solutions for three different engine speeds (2000, 3000 and 4000 rpm) and at various IMEP values. For the knock-free part load points at 2000 rpm (below 9 bar IMEP), the model does not substantially improve the combustion process with respect to the Base combustion parameters (Figure 14a,b). The latter are taken as the average of combustion variables between the two cylinders. The ISFC improvements at low loads (Figure 14c) are exclusively to be attributed to the suppression of the A/F ratio unbalance between cylinders and a maximum As already reported in previous sections, the operating points measured in the entire engine domain (Table 3) were considered to perform a combustion optimization, assuming the same behavior for each cylinder. This means that the numerical optimization process takes into account the absence of non-uniformities among cylinders, which practically requires only partial modifications to the existing fuel injection system of the original engine. The numerical optimization is performed by employing the validated 1D model and the validated CoV correlations referring to the examined engine. It consists of a virtual engine re-calibration, mainly acting on the fuel injection and on the spark timing. In particular, at each measured point in Table 3, a stoichiometric A/F ratio is imposed in each cylinder by a refined control of injected gasoline mass. This condition eliminates the experimental cylinder-to-cylinder variation. At low loads, a PID controller acting on the throttle valve opening is used to match the prescribed IMEP while the waste-gate valve is fully opened. At high loads, a PID controller acts on the waste-gate valve opening to match the target IMEP, while the throttle valve is considered fully opened. The spark timing is modified with the aim to identify the MBT condition and in the case of knock occurrence it is automatically delayed to maintain the numerical knock index below a prescribed threshold level. In addition, at low engine speeds a safety margin is applied to the maximum boost pressure to avoid the compressor surge [36]. At high loads, both the in-cylinder pressure peak and the temperature at turbine inlet are kept below maximum allowable levels to preserve the engine safety. The numerical results obtained by the model-related optimization (labelled as 'Opt') are compared to the model prediction based on the standard ECU calibration (labelled as 'Base'). The above comparisons are realized with the aim to highlight the potential benefits arising from the optimal operation of the engine cylinders. These numerical-derived advantages are here presented in terms of combustion phasing (MFB50) and core duration (MFB10–50), ISFC, and cyclic dispersion levels (CoVIMEP and CoVpmax). The following Figures 14–16 show the comparison of 'Base' and 'Opt' numerical solutions for three different engine speeds (2000, 3000 and 4000 rpm) and at various IMEP values. For the knock-free part load points at 2000 rpm (below 9 bar IMEP), the model does not substantially improve the combustion process with respect to the Base combustion parameters (Figure 14a,b). The latter are taken as the average of combustion variables between the two cylinders. The ISFC improvements at low loads (Figure 14c) are exclusively to be attributed to the suppression of the A/F ratio unbalance between cylinders and a maximum percent gain of 5.6% at 2000@5 is reached. Conversely, at higher loads under knock-limited operation (2000@11 and 2000@13), a combustion phasing advance is realized (Figure 14a), allowing a shortening of combustion core duration (Figure 14b). ISFC benefits at high loads (2.2% for both 2000@11 and 2000@13 in Figure 14c) have to be ascribed to an improved combustion process under in-cylinder stoichiometric mixture conditions. Obviously, the combustion optimization at high load also reflects on the reduction of CoV levels (Figure 14d), while at knock-free loads the CCV remains unchanged. The numerical optimization outcomes at speeds of 3000 and 4000 rpm show similarities. Indeed, for both speeds and under knock-limited operations (above 12 bar IMEP at 3000 rpm and 11 bar IMEP at 4000 rpm), the identified

percent gain of 5.6% at 2000@5 is reached. Conversely, at higher loads under knock-limited operation

optimal combustion phasing is advanced (Figure 15a/Figure 16a) with respect to the Base level and, consequently, the combustion core duration is also reduced (Figure 15b/Figure 16b). Lower CoV levels are predicted by the developed CCV correlations for the optimal solutions, mainly thanks to the improved combustion process (Figure 15d/Figure 16d). Concerning the ISFC advantages, the previous considerations for the case at low speed can be also applied for the medium and high speeds. It is worth highlighting that greater ISFC percent gains are realized in all the high load points where the ECU standard calibration required a rich combustion (Table 3). As confirmed in Figure 15c, higher ISFC benefits are achieved at 3000@12, 3000@16 and 3000@18 with gains of 9.7%, 13.3% and 12.2%, respectively. In a similar way, Figure 16c shows a very high ISFC advantage (equal to 14.7%) only at 4000@16. Finally, at knock-free loads and medium/high speeds (3000@5, 3000@7, 4000@7 and 4000@9) limited improvements in the combustion process and CoV indexes are achieved through the optimization. In these cases, the fuel consumption is still improved (ISFC gain of 4.4% at 3000@5 in Figure 15c and 3.5% at 4000@7 in Figure 16c), due to the elimination of an uneven A/F ratio between cylinders. Summarizing, the suppression of cylinder-to-cylinder variation through the optimization of single cylinder operation demonstrated the possibility of achieving improvements in terms of combustion evolution and stability also contributing to reduced fuel consumption while preserving engine safety both in terms of mechanical and thermal stresses. at 4000 rpm), the identified optimal combustion phasing is advanced (Figure 15a/Figure 16a) with respect to the Base level and, consequently, the combustion core duration is also reduced (Figure 15b/Figure 16b). Lower CoV levels are predicted by the developed CCV correlations for the optimal solutions, mainly thanks to the improved combustion process (Figure 15d/Figure 16d). Concerning the ISFC advantages, the previous considerations for the case at low speed can be also applied for the medium and high speeds. It is worth highlighting that greater ISFC percent gains are realized in all the high load points where the ECU standard calibration required a rich combustion (Table 3). As confirmed in Figure 15c, higher ISFC benefits are achieved at 3000@12, 3000@16 and 3000@18 with gains of 9.7%, 13.3% and 12.2%, respectively. In a similar way, Figure 16c shows a very high ISFC advantage (equal to 14.7%) only at 4000@16. Finally, at knock-free loads and medium/high speeds (3000@5, 3000@7, 4000@7 and 4000@9) limited improvements in the combustion process and CoV indexes are achieved through the optimization. In these cases, the fuel consumption is still improved (ISFC gain of 4.4% at 3000@5 in Figure 15c and 3.5% at 4000@7 in Figure 16c), due to the elimination of an uneven A/F ratio between cylinders. Summarizing, the suppression of cylinder-to-cylinder variation through the optimization of single cylinder operation demonstrated the possibility of achieving improvements in terms of combustion evolution and stability also contributing to reduced fuel consumption while preserving engine safety both in terms of mechanical and thermal stresses.

*Energies* **2020**, *13*, x FOR PEER REVIEW 16 of 22

(2000@11 and 2000@13), a combustion phasing advance is realized (Figure 14a), allowing a shortening of combustion core duration (Figure 14b). ISFC benefits at high loads (2.2% for both 2000@11 and 2000@13 in Figure 14c) have to be ascribed to an improved combustion process under in-cylinder stoichiometric mixture conditions. Obviously, the combustion optimization at high load also reflects

both speeds and under knock-limited operations (above 12 bar IMEP at 3000 rpm and 11 bar IMEP

**Figure 14.** Comparison of Base and Optimized solutions for different engine operating points in a IMEP sweep at constant speed of 2000 rpm: Ave. MFB50 (**a**), Ave. MFB10–50 (**b**), ISFC (**c**), CoVIMEP and CoVpmax (**d**). **Figure 14.** Comparison of Base and Optimized solutions for different engine operating points in a IMEP sweep at constant speed of 2000 rpm: Ave. MFB50 (**a**), Ave. MFB10–50 (**b**), ISFC (**c**), CoVIMEP and CoVpmax (**d**).

CoVpmax (**d**).

*Energies* **2020**, *13*, x FOR PEER REVIEW 17 of 22

**Figure 15.** Comparison of Base and Optimized solutions for different engine operating points in a IMEP sweep at constant speed of 3000 rpm: Ave. Ave. MFB50 (**a**), Ave. MFB10–50 (**b**), ISFC (**c**), CoVIMEP and CoVpmax (**d**). **Figure 15.** Comparison of Base and Optimized solutions for different engine operating points in a IMEP sweep at constant speed of 3000 rpm: Ave. Ave. MFB50 (**a**), Ave. MFB10–50 (**b**), ISFC (**c**), CoVIMEP and CoVpmax (**d**). IMEP sweep at constant speed of 3000 rpm: Ave. Ave. MFB50 (**a**), Ave. MFB10–50 (**b**), ISFC (**c**), CoVIMEP and CoVpmax (**d**).

**Figure 16.** Comparison of Base and Optimized solutions for different engine operating points in a IMEP sweep at constant speed of 4000 rpm: Ave. MFB50 (**a**), Ave. MFB10–50 (**b**), ISFC (**c**), CoVIMEP and **Figure 16.** Comparison of Base and Optimized solutions for different engine operating points in a IMEP sweep at constant speed of 4000 rpm: Ave. MFB50 (**a**), Ave. MFB10–50 (**b**), ISFC (**c**), CoVIMEP and CoVpmax (**d**). **Figure 16.** Comparison of Base and Optimized solutions for different engine operating points in a IMEP sweep at constant speed of 4000 rpm: Ave. MFB50 (**a**), Ave. MFB10–50 (**b**), ISFC (**c**), CoVIMEP and CoVpmax (**d**).

#### **5. Conclusions**

In this work, a small, turbocharged Spark Ignition engine is numerically and experimentally studied with the aim to optimize the individual cylinder combustion and, consequently, to reduce the engine fuel consumption and the cycle-to-cycle variation. In a first phase two different experimental steady analyses were carried out. In particular, a spark advance sweep was performed both at high and low loads, including various engine speeds and air/fuel mixture qualities. The acquired data were post-processed to find proper engine CoV correlations of IMEP and in-cylinder pressure peak with the average combustion phasing and core duration. A second experimental investigation was realized, mainly consisting of a load sweep at three different speeds (namely 2000, 3000 and 4000 rpm), covering both knock-free and knock-limited operations in the engine domain. The operating points considered in this phase of experiments correspond to the conditions imposed by the ECU in the commercial vehicle. The engine was schematized in a 1D model which was validated in terms of global performance variables, combustion evolution and average pressure cycles obtained from the above experimental IMEP sweep. Then, the developed CoV correlations were validated against the experimentally derived CoV values of measured IMEP-sweep points. Once validated, 1D model and the experimentally-related CoV correlations were coupled with the aim to realize a numerical optimization of the engine performance under knock-free and knock-limited conditions. The optimization process reveals that the individual cylinder operation can be considerably improved from ECU-defined settings, especially at high loads. The indicated specific fuel consumption is reduced if compared to the experimental calibration and an ISFC percent benefit ranging from 2.2% up to 14.7% is attained at varying the speed/load points. At low loads, the optimized ISFC level is reached mainly thanks to the possibility of stoichiometric A/F ratio operation in each cylinder. At high loads, ISFC is numerically optimized through the definition of a stoichiometric A/F ratio in each cylinder and an improvement in the combustion process. The optimization at high loads allows the definition of a new setting for the knock-limited spark advance (KLSA) which implies an advanced combustion phasing (MFB50) and a reduced combustion duration if compared to the experimental ECU-defined counterparts. Furthermore, the validated CoV correlations, starting from the numerically optimized combustion parameters, furnish an empirical-based prediction of the cycle-to-cycle variation improvement. Consistently with the combustion optimization, significant reduction in the cycle-to-cycle variation is achieved at high loads for each investigated speed. Summarizing, the adopted numerical procedure allows the identification of the optimal combustion operation, contributing to support the control of each cylinder operation. In this way, the proposed integrated numerical models represent a useful and valuable tool to limit the fuel consumption, the cycle-to-cycle and the cylinder-to-cylinder variations during the engine development phase.

**Author Contributions:** Conceptualization, L.M.; methodology, C.T., L.M. and L.T.; software, L.T.; validation, L.T.; formal analysis, C.T., L.M. and L.T.; investigation, C.T., L.M. and L.T.; data curation, C.T., L.M. and L.T.; writing—original draft preparation, C.T. and L.T.; writing—review and editing, C.T., L.M. and L.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors are grateful to Alfredo Mazzei (Istituto Motori, CNR) for the technical support that has been fundamental for the preparation of the engine test bench.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

#### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Arc-Phase Spark Plug Energy Deposition Characteristics Measured Using a Spark Plug Calorimeter Based on Di**ff**erential Pressure Measurement**

#### **Kyeongmin Kim, Matthew J. Hall \*, Preston S. Wilson and Ronald D. Matthews**

Department of Mechanical Engineering, The University of Texas at Austin, Austin, TX 78712, USA<EMAIL_ADDRESS>(K.K<EMAIL_ADDRESS>(P.S.W<EMAIL_ADDRESS>(R.D.M.) **\*** Correspondence<EMAIL_ADDRESS>Tel.: +1-512-758-1361

Received: 3 June 2020; Accepted: 7 July 2020; Published: 10 July 2020

**Abstract:** A spark plug calorimeter is introduced for quantifying the thermal energy delivered to unreactive gas surrounding the spark gap during spark ignition. Unlike other calorimeters, which measure the small pressure rise of the gas above the relatively high gauge pressure or relative to an internal reference, the present calorimeter measured the differential rise in pressure relative to the initial pressure in the calorimeter chamber. By using a large portion of the dynamic range of the chip-based pressure sensor, a high signal to noise ratio is possible; this can be advantageous, particularly for high initial pressures. Using this calorimeter, a parametric study was carried out, measuring the thermal energy deposition in the gas and the electrical-to-thermal energy conversion efficiency over a larger range of initial pressures than has been carried out previously (1–24 bar absolute at 298 K). The spark plug and inductive ignition circuit used gave arc-type rather than glow-type discharges. A standard resistor-type automotive spark plug was tested. The effects of spark gap distance (0.3–1.5 mm) and ignition dwell time (2–6 ms) were studied for an inductive-type ignition system. It was found that energy deposition to the gas (nitrogen) and the electrical-to-thermal energy conversion efficiency increased strongly with increasing gas pressure and spark gap distance. For the same ignition hardware and operating conditions, the thermal energy delivered to the gap varied from less than 1 mJ at 1 atm pressure and a gap distance of 0.3 mm to over 25 mJ at a pressure of 24 bar and a gap distance of 1.5 mm. For gas densities that might be representative of those in an engine at the time of ignition, the electrical-to-thermal energy conversion efficiencies ranged from approximately 3% at low pressures (4 bar) and small gap (0.3 mm) to as much as 40% at the highest pressure of 24 bar and with a gap of 1.5 mm.

**Keywords:** spark ignition; calorimeter; thermal energy; spark plug; natural gas engine

#### **1. Introduction**

Spark ignition engines rely on the thermal energy deposited by the spark plasma into the in-cylinder gases to initiate combustion each engine cycle. The thermal energy density must be sufficient to initiate and sustain the nascent flame kernel until the chemical heat release from the flame kernel is enough to create a self-sustaining pre-mixed turbulent flame that propagates throughout the combustion chamber on a millisecond time scale. The reliability of the ignition process depends on many factors including fuel type, fuel–air mixture ratio, mixture homogeneity, dilution (e.g., exhaust gas recirculation (EGR) level), residual fraction, in-cylinder bulk flows and turbulence, spark plug geometry, spark timing and duration, and the characteristics of the electrical energy delivery [1]. These factors can impact the amount of thermal energy deposition to the gas from the plasma and its distribution in the vicinity of the spark gap. High dilution rates and lean mixtures pose especially challenging conditions for achieving reliable combustion initiation that minimizes cycle-to-cycle variation.

Ever more detailed ignition submodels that couple to engine computational fluid dynamics (CFD) simulation codes such as Converge CFD® are being developed [2–4]. They are increasingly based on spatially-resolved spark gap physics and can include coupling to dynamic ignition circuit models. Validation of the electrical-to-thermal energy conversion and gas-to-electrode heat transfer rates is needed for these models.

For the above reasons, it is important to understand the characteristics of the thermal energy deposition from the engine spark plasma to the gases. Many fundamental studies of spark breakdown behavior have been conducted. Some of these studies are summarized in Meeks [5]. A wealth of information and data characteristic of the behavior of spark ignition systems of the type used in engines can be found in Maly [6]. By measuring the ignition voltage and current characteristics at the spark plug, it is possible to determine the electrical energy delivered to the spark gap [7], but that does not quantify the thermal energy deposited in the gas in the gap that is responsible for flame initiation. For this purpose, spark plug calorimeters have been developed that will measure the thermal energy deposition. These calorimeters measure the pressure rise of the gas in a small chamber in which the spark plug is located to determine the thermal energy input to the gas.

This work reports on the development of a spark plug calorimeter that uses very sensitive interchangeable pressure sensors to measure the differential pressure rise of the chamber gases relative to the initial chamber pressure and is designed to function at very high pressures (>20 bar). The calorimeter was then used to investigate the thermal energy deposition and efficiency of electrical-to-thermal energy conversion of arc-type spark plug discharges for an inductive ignition system for different pressures, gap distances and dwell times.

High pressure measurements are of increasing interest due to the current push by many engine manufacturers to develop ever higher brake mean effect pressure (BMEP) engines that rely on increasingly high boost pressures. Most previous studies that have used spark calorimeters have relied upon pressure sensors that measured the spark plug chamber gauge pressure relative to the ambient atmospheric pressure or are designed for use at high pressures and measure relative to an internal reference. By measuring the differential chamber pressure relative to the initial pressure, more sensitive measurements with a higher dynamic range are possible.

Several studies that used spark plug calorimeters are found in the literature. Roth et al. [8] first introduced the measurement of spark energy delivered to the gas with both constant-volume and constant-pressure calorimeters. They performed a fundamental study of the effects of the electrode diameter, gap distance, and thermal diffusivity of the gas on the electrical-to-thermal energy conversion efficiency with monatomic gases. They found that more of the electrical energy from the spark was lost to the electrodes with larger electrode diameter, shorter gap distance, and gases with larger thermal diffusivity. Merritt [9] was first to create a spark calorimeter consisting of two chambers and a differential pressure transducer. It was constructed of plastic and was used to measure spark energy thermal deposition at atmospheric pressure. The extra reference volume was said to eliminate the short-term fluctuations in atmospheric pressure. This differential pressure sensing concept for a spark calorimeter was adopted by Franke and Reinmann [10]. They measured the energy delivered to the gas comparing five different ignition systems that included both capacitive and inductive ignition types at pressures up to 16 bar; all of the systems delivered breakdown followed by glow discharges, without a discernable arc phase. The mean electrical-to-thermal energy conversion efficiency increased from 5% to 50% as the pressure varied from 1 to 16 bar, and for a given level of delivered electrical energy, the efficiency increased as the duration of the discharge increased.

While the concept of a spark plug calorimeter based on differential pressure measurement is the same here as that presented by Franke and Reinmann, the present paper goes beyond their work in several significant ways, both in terms of the experimental technique and the conditions analyzed. With respect to the calorimeter design, one difference is that in this study we used inexpensive and

easily changeable solid-state semiconductor pressure sensors with different sensitivities such that we could exchange them to adapt to different expected pressure rises, allowing us to maximize the dynamic range of the measurements. Another novel feature of this calorimeter is the system that allowed the pressure sensor to be calibrated at high pressure for very small pressure differentials. The relatively small calorimeter volume resulted in greater pressure rises for a given thermal energy deposition, further increasing the potential accuracy of the measurements.

With respect to the differences in experimental conditions, the present study extended the pressure range of the measurements by 50% relative to that of the study by Franke and Reinmann. This is especially significant for application to advanced natural gas engines under development by various engine manufacturers. These engines are designed for ever higher boost conditions and high BMEPs that result in ever higher gas densities at the time of ignition. In addition, this study is primarily focused on arc-type discharges. Franke and Reinmann stated that "In this experiment, no transition to arc mode has been observed", implying that they believed that the plasma discharges that they observed were in the glow regime. Further, in this study, we present empirically-based analytical expressions for how differential or incremental changes in thermal energy deposition scale with changes in either gap distance or gas pressure.

Other spark calorimeter studies found in the literature used calorimeters that measured the chamber pressure relative to the outside ambient pressure, a simpler method, but one that compromises sensitivity and dynamic range, particularly at high initial pressures. Teets and Sell [11] used such a calorimeter to study the thermal energy deposition characteristics of three different ignition systems that included an inductive system, a plasma jet ignitor and an ultra-short pulse ignitor. Measurements were to pressures up to 7 bar. They found conversion efficiencies varied from 5% to 65%, increasing with pressure and gap distance. They found that electrical-to-thermal energy conversion efficiency decreased as the delivered electrical energy increased, and also found the efficiency of the ultra-short pulse system higher than the inductive system.

Verhoeven [12] took a different approach to measure thermal energy deposition to the gas from a spark plasma, using an optical holographic technique. Uncertainties associated with the technique were greater than for calorimeters, but it had the advantage of allowing measurements in flowing gases. However, the flow measurements were only made at a pressure of 1 atm. They found that with a cross-flow velocity of 5 m/s, energy deposition efficiency was approximately three-fold that for a quiescent gas.

In a limited study, Abidin et al. [13] investigated the effect of dwell time, gap distance, and pressure on the breakdown voltage and the electrical-to-thermal energy conversion efficiency using a standard spark calorimeter. The electrical-to-thermal efficiency increased with larger gap distance, shorter dwell time, and higher pressure. However, the highest pressure investigated was 9 bar. Alger et al. [14] used a spark calorimeter to measure the energy delivered to the gas while investigating the effect of spark plug design on the initial flame kernel development. They found typical electrical-to-thermal energy conversion efficiencies of approximately 20% and found that spark plugs with high internal resistances had higher thermal energy depositions, which resulted in faster flame kernel development.

In this study, a spark plug calorimeter was used to measure the thermal energy delivered to the nitrogen gas surrounding the spark gap during spark ignition. Unlike most other spark plug calorimeters, which measure the small pressure rise of the gas above the relatively high gauge pressure or relative to an internal reference, the present calorimeter measured the differential rise in pressure relative to the initial pressure in the calorimeter chamber. Using this calorimeter, a study was carried out, measuring the thermal energy deposition in the gas and the electrical-to-thermal energy conversion efficiency over a larger range of initial pressures than has been carried out previously. The measurements were made at pressures up to 24 bar. A pressure of 24 bar at the ambient temperature at which the measurements were made (approximately 22 ◦C) results in gas densities equivalent to approximately 50 bar for a gas temperature at the time of ignition of 600 K or 70 bar at 900 K. The spark plug and inductive ignition circuit used gave arc-type rather than glow-type discharges. For example, for a reasonable intercooler exit temperature of 350 K, with this also taken as the intake temperature, and assuming a compression ratio of 10 while assuming polytropic compression with a typical polytropic exponent of 1.33, yields a TDC temperature of approximately 750 K. An ignition timing at or near TDC is often utilized in modern SI engines to enhance catalyst heating after startup.

#### **2. Calorimeter Design**

A photo of the calorimeter is shown in Figure 1a, along with a cross-sectional drawing shown in Figure 1b. It was machined from stainless steel in two pieces and accommodates a 14 mm spark plug. The calorimeter has two chambers, a small cylindrical chamber into which the spark plug is inserted and a second chamber in which the chip-based pressure sensor is located. The pressure sensor is mounted through a small hole that connects the two chambers. Two valves are used. One valve is used to fill the two chambers with gas (with both valves open). That valve is then closed to stop the flow of gas into the chambers. The second valve is then closed, which isolates the two chambers except for the small hole between them into which the pressure sensor entrance passage is sealed. The spark plug chamber is cylindrical in shape, with a diameter of 12.8 mm, a height of 9.5 mm, and a volume of 1.8 cm3 . The pressure sensor is sealed within the second cavity of the calorimeter. A high pressure electrical feedthrough is passed the necessary wires for signal and power source to the outside.

(**b**)

**Figure 1.** (**a**) Photo of spark calorimeter; (**b**) section drawing of the calorimeter.

The pressure sensors used were board mounted integrated silicon pressure sensors. Two types were used for the measurements presented here—one with a 3.9 kPa pressure range and another with a 10 kPa pressure range—depending on the expected pressure levels. The model numbers of the two sensors available from Freescale (USA) were NXP MP3V5004GC6U and MP3V5010GC6U for the 3.9 and 10 kPa sensors respectively. The 1/e response time of the sensors was approximately 0.2 ms. This was fast enough to quantify the trends in the cumulative energy deposition of an ignition event, but was not fast enough to resolve the thermal energy deposition associated with the breakdown process, which occurs on a submicrosecond time scale.

Figure 2 shows a schematic of the calorimeter experimental setup. For the measurements, the calorimeter was pressurized with nitrogen gas to the desired level. The charging pressure was measured using precision Bourdon tube pressure gauges. Time-resolved measurements of spark

plug voltage and current were taken to determine the electrical energy delivered to the spark plug. A Tektronix Model P6015A high voltage probe measured the breakdown and follow-on voltages at the top of the spark plug. The current-dependent resistance of the plug was determined and the IR voltage subtracted from the voltage measured at the top of the plug to obtain the gap voltage. A Pearson Model 110 current sensor was used to measure the discharge current. The voltage, current, and pressure sensor signals were recorded using a 100 MHz 4-channel Tektronix oscilloscope. Breakdown voltages were recorded separately since a faster time-base setting was needed to resolve these very short duration events.

**Figure 2.** Schematic of calorimeter setup.

The chip-based sensors used to measure the change in chamber pressure that resulted from the spark discharge were intended, by the manufacturer, to measure small differential pressures relative to the ambient barometric pressure. Since the intention was to make measurements at high pressure, the sensors needed to be tested for their response at high pressure. To do this, a calibration procedure was developed that tested the pressure sensors in situ. The calibration device is depicted in Figure 3. A glass pipette was inserted through the spark plug hole and coupled to the passage tube of the pressure sensor with a short piece of plastic tubing. The pipette was then filled with water to a water column height that gave the desired pressure differential, the pipette thus acting as a manometer. A closed-end stainless steel tube was then fit over the pipette and secured with tube fittings. The two chambers of the calorimeter were when charged with nitrogen to the same gas pressure; the difference in pressure between the two chambers at the pressure sensor was the pressure imposed by the vertical column of water. This method allowed the pressure sensors to be calibrated for different high initial pressures and provided pressure measurements relative to the initial pressure.

The calibration data for the 10 kPa range sensor is shown in Figure 4, which shows the measured voltage output of the pressure sensor versus the differential pressure imposed by the water column height for initial total gas pressures from 1 atm to 24 bar. The response of the pressure sensor is shown over the measured range; there was no systematic differential response to the initial pressure. The scatter/reproducibility of the calibration data provides an estimate of the uncertainty in the accuracy of the pressure measurements of less than ± 0.3 kPa.

**Figure 3.** Photo of calorimeter with calibration tube attached.

**Figure 4.** Calibration graph for the differential pressure sensor.

#### **3. Results and Discussion**

A standard representative 14 mm automotive spark plug (Champion RS12YC) was used for all of the tests. The center electrode (cathode) had a diameter of 2.5 mm and the width of the J-type ground strap was also 2.5 mm. It had a nominal internal resistance of 6.1 kohms. The inductive-type ignition coil had a secondary inductance value of 30 H. The default dwell time used was 4 ms. However, the effect of dwell times of 2 and 6 ms were studied as well. Figure 5a presents an example of the measured current delivered to the spark plug and a voltage trace corrected for the voltage drop across the internal resistor, and thus, represents the gap voltage. The integrated product of the corrected gap voltage and the current was used to calculate the electrical energy delivered to the gap. These representative traces were for conditions of 8 bar absolute gas pressure, with a 0.9 mm gap and a dwell time of 4 ms. Figure 5b shows a typical voltage trace having a short time scale as used to measure the breakdown voltage. The breakdown voltage for this trial was 8.67 kV and the spark duration was 3.45 ms.

**Figure 5.** (**a**) Current and voltage vs. time for entire discharge duration and (**b**) voltage vs. time for breakdown.

It is interesting to note that for all gas pressures above 1 atm, the ignition system produced an arc-type discharge, whereas for 1 atm pressure and the smallest spark gap of 0.3 mm, an arc to glow transition occurred immediately following breakdown. This was readily distinguished by the level of the relatively flat voltage profile following breakdown. The arc discharges were in a range of approximately 100–300 volts, whereas glow discharges were in the range of approximately 400–500 volts. For all of the measurements, the peak current following breakdown was approximately 100 mA and decreased nearly linearly with time, but with different spark durations.

Figure 6a shows an example of recorded pressure traces for different initial gas pressures ranging from 1 to 24 bar. In general, the pressure rise was very rapid during the early part of the discharge, consistent with the high current levels during that period and representative of the large energy deposition during the breakdown process and shortly thereafter.

**Figure 6.** (**a**) Recorded pressure rise vs. time for different initial gas pressures and (**b**) smoothed pressure traces used for analysis. Gap distance 0.9 mm; dwell 4 ms.

As seen in Figure 6a, there is a periodic ringing superimposed upon the pressure signal. For the 1 atmosphere case, this sinusoidal pressure variation superimposed on the signal has a frequency of approximately 2.4 kHz. Considerable effort was made to understand the cause of this. It was initially attributed to an acoustic mode excited in the short (4 mm in length) passage tube that was part of the pressure sensor. To test for this, the passage tube was cut off and a new holder for the pressure sensor was machined to fit and seal the new sensor geometry. Subsequent tests revealed that the ringing

remained but that it occurred at an increased frequency of approximately 6.1 kHz. At this point, a more detailed acoustic analysis was performed that considered possible acoustic modes associated with the calorimeter chamber and also considered the system response of the pressure transducer.

The calorimeter chamber is to first order a right circular cylinder with rigid acoustic boundary conditions and is the largest acoustic cavity in the system, and hence must exhibit the lowest acoustic cavity resonance frequency. Given the aforementioned dimensions, and considering the 1 atmosphere case and room temperature, the lowest-order resonance frequency is for approximately one-half wavelength acoustic fluctuations standing in the largest spatial dimension (12.8 mm), which yields a resonance frequency of approximately 13.4 kHz. All other acoustic standing wave modes must be higher in frequency. For example, the next largest dimension is the length of the calorimeter chamber, which at 9.5 mm yields approximately 18 kHz for the lowest-order resonance frequency. These frequencies are significantly higher than the fluctuations observed in the experimental data, and hence the observed ringing should not be attributed to acoustic standing wave modes in the system. To verify this, a spark plug and the pressure sensor were placed close together outside of the calorimeter, and operated in free space. A spark discharge was measured and the ringing in the pressure signal was still observed.

Below the lowest-order standing wave modal frequency, acoustic lumped-parameter behavior can occur and the pressure sensor contains elements that form a Helmholtz resonator. The air inside the passage tube acts as an acoustic mass, and the air inside the cavity within the sensor acts as an acoustic compliance. Together, these form a simple harmonic oscillator, which can be excited by a broadband source of energy, such as the shock generated by the spark [15]. Examination of both the sensor's specification sheet and the sensor itself yielded estimates of the dimensions of the structures required for the calculation of the Helmholtz resonance frequencies using

$$f\_0 = \frac{c}{2\pi} \left(\frac{S}{L'V}\right)^{\frac{1}{2}}\tag{1}$$

where *c* = 343 m/s is the speed of sound at room conditions, *S* is the inner area of the passage tube, and *V* is the volume of the cavity. The effective length of the passage tube is *L* ′ = *L* + 1.4*a*, where *L* is the actual length of the passage tube and *a* is its inner radius [16]. For the case with the passage tube attached, the predicted Helmholtz resonance frequency is 2.4 kHz and the observed frequency was 2.1 kHz. With the passage tube removed, *L* = 0, and the appropriate effective length is *L* ′ = 1.7*a* [16]. The predicted Helmholtz resonance frequency is 6.4 kHz and the observed frequency was 6.1 kHz. Both observations match the predictions within the uncertainty of the estimate of the dimensions, hence the observed fluctuations appear to be an artifact due to Helmholtz resonator behavior of the sensor.

The tabulated sensor design specifications include a first-order system time constant of 0.45 ms, which in turn yields a low-pass system response with an upper band limit of approximately 350 Hz, hence the sensor is not designed to accurately measure the effects of the initial shock, nor dynamic pressure fluctuations above approximately 350 Hz. Before performing the analyses relating the pressure rise to the thermal energy deposition, the pressure curves were smoothed to remove the Helmholtz resonator artifacts, as shown in Figure 6b.

The thermal energy deposition to the chamber gas was derived from the measured pressure rise using Equation (2) [13]. Of the electrical energy delivered to the gap, a portion is retained as thermal energy within the gas with the remainder lost as heat transfer, primarily to the electrodes. Equation (2) expresses the effect of the retained thermal energy of the gas on the rise in pressure assuming chemical and thermal equilibrium and using the simplifying assumptions of constant specific heat and ideal gas behavior.

$$E\_{therm} = \frac{V}{\gamma - 1} \Delta P \tag{2}$$

In Equation (2), *V* is the chamber volume, ∆*P* is the maximum pressure rise and γ is the ratio of specific heats of nitrogen. The relative contributions of possible heat loss mechanisms are not clear. The difference in temperature between the plasma/gas and the spark plug surfaces will contribute to thermal heat transfer, while there may also be Ohmic losses along the surfaces of the electrodes. To be representative of an engine, the heat losses should reflect the sources of loss in an engine, and not for example, due to heat loss to the walls of the calorimeter. A small decrease in chamber pressure was sometimes observed toward the end of the spark discharge process (e.g., Figure 6) indicating continued heat transfer. It is likely that this thermal heat transfer was to the spark plug electrodes and not to the chamber walls over this period since an estimate of molecular thermal diffusion to the chamber walls yields a characteristic time of 100 s of milliseconds. In this study, we defined the amount of thermal energy delivered to the gas to be that reflected by the maximum of the measured pressure rise. Over the short duration of the spark, the gas is not in thermal equilibrium, as temperature gradients exist. As a result, the application of Equation (2) implies a spatial mass average of the thermal contributions to the pressure of the different temperature layers. This should be valid because of the linear relationship between mass, temperature and the pressure as given by the ideal gas law and the linear relationship between thermal energy and temperature inherent in the constant specific heat assumption. In addition, during the spark discharge, not all of the gas exists as N2 since there is dissociation of the nitrogen in the spark plasma; however, only a small fraction of the gas is so affected and equilibrium processes will have largely returned this gas to N2 by the time the spark event has ended.

The measured spark breakdown voltages as functions of the initial gas pressure and spark gap are shown in Figures 7 and 8 for a dwell time of 4 ms. Each data point represents an average of five measurements.

**Figure 7.** Breakdown voltage vs. initial pressure for different gaps, for a dwell time of 4 ms.

**Figure 8.** Breakdown voltage vs. gap for different pressures, for a dwell time of 4 ms.

As expected, the breakdown voltage increases with both the gas pressure and the gap distance, but it is interesting to note that the breakdown voltage was not directly proportional to either of these variables even though linear relationships are predicted by Paschen's law [5]. In the case of gas pressure, the breakdown voltage increased at a rate that diminished as the gas pressure increased. While the trend was closer to a linear increase for a changing spark gap, there was a considerable zero offset for the breakdown voltage at all of the pressures for the data extrapolated back to zero gap. In Figures 9 and 10, the of effect dwell time on breakdown voltage is shown. There was no appreciable effect of either gap distance or gas pressure on the breakdown voltage within the uncertainty of the measurements over the dwell times ranging from 2 to 6 ms; however, a very slight increase in breakdown voltage with increasing dwell time is suggested from Figure 9. Since a longer dwell time results in more energy stored in the secondary coil and since the energy stored in an inductor follows Equation (3),

$$E = \frac{1}{2}LI^2\tag{3}$$

where *L* is the coil inductance and *I* is the current, it follows that the current will be greater in the inductor at higher energies. It may be the case that for the higher secondary currents associated with greater stored energies, the electric field in the gap rises slightly faster prior to breakdown, resulting in a higher voltage at the gap before the actual breakdown occurs.

**Figure 9.** Breakdown voltage vs. dwell for different gaps, for a pressure of 12 bar.

**Figure 10.** Breakdown voltage vs. dwell for different pressures, for a 0.9 mm gap.

The discharge time or spark duration as a function of gas pressure, gap distance, and dwell time is given in Figures 11–14. As seen in Figure 11, the spark duration is a weak function of gas pressure, but generally decreases with increasing pressure. The exceptions to this were the measurements at

1 atm, which show a shorter spark duration. This was associated with the glow-type discharge that was mentioned previously. The discharge times were in the range of approximately 3–4 ms. Figure 12 shows that the discharge times tended to shorten slightly as the gap distance increased; however, there is an exception for the 1 atm measurements. Figures 13 and 14 show that spark duration increased with dwell time as more electrical energy was delivered to the gap, with the increase weaker from 4 to 6 ms than from 2 to 4 ms.

**Figure 11.** Discharge time vs. pressure for different gaps, for a dwell time of 4 ms.

**Figure 12.** Discharge time vs. gap for different pressures, for a dwell time of 4 ms.

**Figure 13.** Discharge time vs. dwell for different gaps, for a pressure of 12 bar.

**Figure 14.** Discharge time vs. dwell for different pressures, for a 0.9 mm gap.

Figures 15–18 show the calculated electrical energy delivered to the gap based on the measured current traces and the corrected traces of gap voltage. From Figures 15 and 16, it can be seen that the electrical energy delivered to the gap generally increased, both with increasing gas pressure and with gap distance. As both the gas pressure and gap distance increase, the breakdown voltage increases as well; increasing the peak voltage in the secondary side of the ignition circuit will lead to greater electrical energy storage according to the secondary circuit and spark plug capacitance (typically on the order of 15 pF [10]). Both higher gas pressure and larger gap result in a greater gap resistance and in a more thermally isolated gap region since molecular diffusivities scale as the square of the characteristic diffusion distance. It is interesting that electrical energy delivered to the gap increased as these two parameters increased even though spark duration decreased. The relatively high electrical energy delivery for the case of 1 atm gas pressure and the smallest gap of 0.3 mm was a consequence of the high glow voltage relative to the arc voltages; the current traces were similar for arc and glow.

**Figure 15.** Electrical energy input vs. pressure for different gaps, for a dwell time of 4 ms.

**Figure 16.** Electrical energy input vs. gap for different pressures, for a dwell time of 4 ms.

**Figure 17.** Electrical energy input vs. dwell for different gaps, for a pressure of 12 bar.

**Figure 18.** Electrical energy input vs. dwell for different pressures, for a 0.9 mm gap.

The electrical energy delivered to the gap as a function of dwell time is shown in Figures 17 and 18. As expected, a longer dwell time leads to greater electrical energy storage in the secondary side of the coil and, consequently, more energy delivered to the gap. It is interesting to note the range of delivered electrical energy to the gap as parameters changed; delivered energies ranged from approximately 15

to 70 mJ. The literature shows minimum ignition energies for stoichiometric hydrocarbon-air mixtures that are in the range of 0.3–2 mJ [17].

The thermal energy deposition to the gas is depicted in Figures 19–22. The quantity of thermal energy delivered to the gas and its distribution in the vicinity of the spark gap is important for flame initiation, particularly for conditions for which the development of a strong flame kernel is challenging, such as in lean, highly dilute and inhomogeneous mixtures and where fluid motion near the gap can strongly convect both the spark plasma and the developing flame kernel.

Figure 19 shows the thermal energy deposition versus gas pressure as a function of gap distance and for a dwell time of 4 ms. The thermal energy delivered to the gas is seen to increase monotonically with increasing pressure. If one considers the results for the different gap distances, the normalized change in thermal energy deposition with pressure is similar for all of them; that is

$$\frac{E(P\_2) - E(P\_1)}{P\_2 - P\_1} \approx \text{constant} \tag{4}$$

where *E* is the thermal energy to the gas in mJ and *P* is the gas pressure in bar. For the pressures shown, the value of the constant in Equation (4) is in the range of 0.2–0.4 mJ/bar, with an average close to 0.3 mJ/bar. Using this average value of 0.3 mJ/bar and rearranging Equation (4) as a differential, one obtains the following relationship for the differential change in the normalized thermal energy deposition with gas pressure that is independent of gap distance. ܲଶ − ܲଵ

$$\frac{dE}{E} = \left(0.3 \text{bar}^{-1} dP\right) - 1\tag{5}$$

It would require more investigation to determine the validity of Equation (5) for different spark plug geometries and internal resistances, the type of ignition circuit, dwell times, and of course, gap fluid motion, but it may provide an approximate scaling rule in the absence of other data. It may be more appropriate to express this relationship in terms of gas density rather than pressure through their relationship via the ideal gas law.

Figure 20 shows the thermal energy deposition versus gap distance as a function of initial absolute gas pressure. Thermal energy deposition increased monotonically with gap distance for the entire range of gas pressures. Differences in gap electrical resistance and differences in heat transfer to the electrodes are probably responsible for these trends. As the gap distance increases, the positive column portion of the spark plasma increases proportionally, along with its electrical resistance. Molecular thermal diffusivities increase as the square of the characteristic heat transfer length, which one might take as one-half the gap distance. Despite the possible non-linearity, it is interesting to note that increases in thermal energy deposition with gap distance follow an approximately linear trend.

**Figure 19.** Thermal energy input vs. pressure for different gaps, for a dwell time of 4 ms.

− − experimental data is in the range of approximately 4–6 mm−1 , with an average value close to 5 mm−1 . Equation (6) is the resulting relationship for the effect of gap distance "*g*" on the change in thermal energy deposition that is independent of gas pressure

$$\frac{dE}{E} = \left(5\text{mm}^{-1}d\text{g}\right) - 1\tag{6}$$

Again, the sensitivity to other parameters needs to be explored further.

Figures 21 and 22 show the thermal energy deposition as a function of dwell. A longer dwell time results in more electrical energy delivered to the gap, so the observed trend of increasing thermal energy deposition with increasing dwell time would be expected. ܧ݀ ܧ = ሺ5 ିଵ݀݃ሻ − 1

Finally, the conversion efficiencies of electrical energy to thermal energy delivered to the gap are shown in Figures 23–26. Figure 23 shows the conversion efficiency as a function of initial gas pressure for the different gap distances. Conversion efficiency increased monotonically with gas pressure, but the dependency was not linear.

**Figure 21.** Thermal energy input vs. dwell for different gaps, for a pressure of 12 bar.

**Figure 22.** Thermal energy input vs. dwell for different pressures, for a 0.9 mm gap.

The degree to which the conversion efficiency increased with gas pressure was quite striking. At a pressure of 1 atm and for the smallest gap distance of 0.3 mm, the conversion efficiency was less than 1%, whereas with the largest gap of 1.5 mm and at the highest gas pressure of 24 bar, the conversion efficiency was over 40%. The rate of increase in the conversion efficiency with pressure tended to diminish with increasing pressure, in contrast to the thermal energy deposition for which the rate of increase showed an increase with increasing gas pressure.

Of course, the gas density associated with the room temperature measurements at 1 atm are not representative of the gas density in an engine at the time of ignition. However, the laboratory gas densities at pressures starting at 4 bar become relevant to engine ignition. At 4 bar pressure, conversion efficiencies were less than 15% for gap distances less than 1 mm.

The dependence of conversion efficiency on spark gap distance was also very strong, as depicted in Figure 24. A smaller gap distance effectively increases the ratio of the cold metal surface area subject to heat transfer relative to the volume of the spark plasma, leading to greater heat losses from the arc.

The conversion efficiencies versus dwell time are shown in Figures 25 and 26. A long dwell time results in greater electrical energy supplied to the gap. However, it is interesting that, in general, the efficiency of electrical-to-thermal energy conversion diminished with increasing dwell time. The reason for this is not clear, but longer dwell time resulted in longer spark durations and, consequently, more time for heat transfer to occur.

**Figure 23.** Conversion efficiency vs. pressure for different gaps, for a dwell time of 4 ms.

**Figure 24.** Conversion efficiency vs. gap for different pressures, for a dwell time of 4 ms.

**Figure 25.** Conversion efficiency vs. dwell for different gaps, for a pressure of 12 bar.

**Figure 26.** Conversion efficiency vs. dwell for different pressures, for a 0.9 mm gap.

#### **4. Summary and Conclusions**

A spark plug calorimeter was designed and built that utilizes a high sensitively chip-based differential pressure sensor to detect the rise in gas pressure associated with the thermal energy deposition to the gas from the spark discharge. It differs from most spark plug calorimeters in that the differential pressure rise is measured relative to the initial pressure of the calorimeter chamber rather

than relative to the outside ambient pressure or internal reference. This feature creates the possibility of higher precision measurements from a wider dynamic range of the sensor and its high sensitivity. The strength of its stainless steel construction allowed tests to be safely performed at gas pressures at least as high as 24 bar.

The present calorimeter was used to measure the thermal energy deposition and the efficiency of electrical energy to thermal energy conversion for arc discharges for a standard 14 mm automotive resistor-type spark plug for gas pressures between 1 atm and 24 bar, gap distances between 0.3 and 1.5 mm, and inductive ignition circuit dwell times between 2 and 6 ms.

Spark breakdown voltages were measured as well. They were in the expected range, based on the literature, but they did not closely follow the expected linear trends relative to changes in gas density and gap distance suggested by Paschen's law (e.g., Figures 7 and 8), with the rate of increase in breakdown voltages with increasing gap and increasing pressure less than linear.

Arc-type, as opposed to glow-type spark discharges, were investigated, as distinguished by the sustaining arc voltages during discharge in the range of 100 V–200 V (e.g., Figure 5). The measured electrical energies delivered to the gap increased with both increasing gas pressure and gap distance, but the thermal energy delivered to the gap increased much more strongly. For example, delivered electrical energies, for the larger gaps, increased by approximately 50–90% between 4 and 24 bar pressures, whereas thermal energy deposition increased approximately 200% over the same pressure range. For the same ignition hardware and operating conditions, the thermal energy delivered to the gap varied from less than 1 mJ at 1 atm pressure and a gap distance of 0.3 mm to over 25 mJ at a pressure of 24 bar and a gap distance of 1.5 mm (e.g., Figures 19 and 20).

These trends were reflected in the electrical-to-thermal energy conversion efficiency, which increased monotonically with gas density. For gas densities representative of those in an engine at the time of ignition, the conversion efficiencies ranged from approximately 3% at low pressures (4 bar) and small gap (0.3 mm) to as much as 40% at the highest pressure of 24 bar and with a gap of 1.5 mm (e.g., Figures 23 and 24).

For the particular spark plug and ignition circuit used, empirically calibrated expressions were derived for the incremental change in thermal energy deposition with changes in either gas pressure or gap distance.

**Author Contributions:** K.K. was responsible for data curation, validation and formal analysis. P.S.W. was responsible for the acoustic analysis of the calorimeter/sensor system. R.D.M. and M.J.H. were responsible for the conceptualization and supervision of the project. All authors participated in the writing, review, and editing of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project was made possible through the funding provided by Cummins Inc, through the University of Texas at Austin's site of the National Science Foundation Center for Efficient Vehicles and Sustainable Transportation Systems (EVSTS), NSF Award 1650483.

**Acknowledgments:** The authors wish to express their gratitude to Sachin Joshi, Daniel J. O'Connor and Douglas L. Sprunger of Cummins Inc. for many helpful discussions and we wish to thank Cummins Inc. for their financial support of this project.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Experimental Study of Premixed Gasoline Surrogates Burning Velocities in a Spherical Combustion Bomb at Engine Like Conditions**

#### **Miriam Reyes \* , Francisco V. Tinaut and Alexandra Camaño**

Department of Energy and Fluid-Mechanics Engineering, University of Valladolid, Paseo del Cauce, 59, 47011 Valladolid, Spain<EMAIL_ADDRESS>(F.V.T<EMAIL_ADDRESS>(A.C.)

**\*** Correspondence<EMAIL_ADDRESS>
Received: 10 June 2020; Accepted: 1 July 2020; Published: 3 July 2020

**Abstract:** In this work are presented experimental values of the burning velocity of iso-octane/air, n-heptane/air and n-heptane/toluene/air mixtures, gasoline surrogates valid over a range of pressures and temperatures similar to those obtained in internal combustion engines. The present work is based on a method to determine the burning velocities of liquid fuels in a spherical constant volume combustion bomb, in which the initial conditions of pressure, temperature and fuel/air equivalence ratios can be accurately established. A two-zone thermodynamic diagnostic model was used to analyze the combustion pressure trace and calculate thermodynamic variables that cannot be directly measured: the burning velocity and mass burning rate. This experimental facility has been used and validated before for the determination of the burning velocity of gaseous fuels and it is validated in this work for liquid fuels. The values obtained for the burning velocity are expressed as power laws of the pressure, temperature and equivalence ratio. Iso-octane, n-heptane and mixtures of n-heptane/toluene have been used as surrogates, with toluene accounting for the aromatic part of the fuel. Initially, the method is validated for liquid fuels by determining the burning velocity of iso-octane and then comparing the results with those corresponding in the literature. Following, the burning velocity of n-heptane and a blend of 50% n-heptane and 50% toluene are determined. Results of the burning velocities of iso-octane have been obtained for pressures between 0.1 and 0.5 MPa and temperatures between 360 and 450 K, for n-heptane 0.1–1.2 MPa and 370–650 K, and for the mixture of 50% n-heptane/50% toluene 0.2–1.0 MPa and 360–700 K. The power law correlations obtained with the results for the three different fuels show a positive dependence with the initial temperature and the equivalence ratio, and an inverse dependence with the initial pressure. Finally, the comparison of the burning velocity results of iso-octane and n-heptane with those obtained in the literature show a good agreement, validating the method used. Analytical expressions of burning velocity as power laws of pressure and unburned temperature are presented for each fuel and equivalence ratio.

**Keywords:** iso-octane; n-heptane; toluene; surrogate fuels; burning velocity; combustion bomb

#### **1. Introduction**

The laminar burning velocity is a fuel property of fundamental importance for predicting and studying the performance of internal combustion engines, which can be extremely useful in the analysis of fundamental processes and serve as a design utility during the engine design stage. The laminar burning velocity of premixed flames has been the focus of comprehensively experimental and numerical investigations. This property is essential in the analysis, design and implementation of internal combustion engines (ICE), because the burning velocity has a direct influence on the

efficiency, emissions and burn rate in the combustion engine [1]. Burning velocities are implemented in combustion models to validate kinetic mechanisms. The laminar burning velocity of fuel mixtures can be obtained using several techniques, typically by combining experimental methods (pressure register, optical) with thermodynamic models and instabilities studies to consider the effects of flame stretch.

Precise measurements of values of the burning velocity are necessary to characterize fuels in premixed combustions and to validate combustion models [2]. A technique widely used in the literature, to present the results of the combustion rates of different fuels, is a correlation (power law) based on the initial conditions of the pressure, temperature and equivalence ratio of a given fuel mixture for a determined range of the pressure and temperature. These correlations provide very important information for evaluating the effect of fuels in spark ignition (SI) engines.

Due to the complexity of the gasoline composition and its variability, gasoline surrogates are used both experimentally and numerically to simplify calculations. Gasoline surrogates refers to fuels with a simpler representation of a fully mixed fuel, which can be primary reference fuels or binary mixtures based on the research octane number. Some of these gasoline surrogates are used as fuels for advanced combustion engines [3]. Surrogate fuels have a simpler composition which can facilitate simulations and evaluate the property effects and fuel-composition of the in-cylinder processes: vaporization, mixing and combustion, determining the processes in engine efficiency, emissions, performance, and requirements about the after-treatment systems [4]. Therefore, surrogate fuels have a great value as reference fuels which can be used to determine different parameters of interest in engine combustion without the effect of changes in the fuel-composition. The most common gasoline surrogates are iso-octane and n-heptane, the primary reference fuels (usually named PRF's) or binary mixtures of them for determining the research octane numbers (RON) in SI fuels.

Iso-octane and n-heptane are standard gasoline surrogates for modeling combustion in SI engines because the processes of oxidation of n-heptane and iso-octane represent the ignition process and combustion characteristics of a gasoline fuel. For that reason, the effect of composition changes of multicomponent fuels can be simulated using binary mixtures of them as a first approximation [5].

Some investigations have been developed to determine the laminar burning velocity of iso-octane, n-heptane and binary mixtures of these two fuels, and, for that reason, iso-octane and n-heptane are used in this work to validate the methodology presented. Gülder [6], Bradley et al. [7], Metghalchi et al. [8], Galmiche et al. [9] and Müller et al. [10] obtained different correlations for the burning velocity of iso-octane as a function of pressure and temperature for diverse fuel/air equivalence ratios. Galmiche et al. [9] measured flame velocities in spherically expanding flames, from which the corresponding laminar burning velocities at a stretch rate of zero are derived. Iso-octane/air mixtures at initial temperatures between 323 and 473 K, and pressures between 0.1 and 1.0 MPa are studied over an extensive range of equivalence ratios, using a high-speed shadowgraph system. Varea et al. [11] used a high pressure and temperature combustion chamber to obtain the laminar burning velocity of iso-octane with a new method, from the difference between the flame speed and the intake gas velocity. They also obtained correlations as a function of pressure and temperature to express their effect. Marshal et al. [12] obtained the laminar burning velocity of iso-octane and n-heptane at elevated pressures and temperatures and they included combustion residuals in their experiments. They used a constant volume combustion vessel equipped with a Schlieren technique and obtained values for different ranges of pressure, temperature and equivalence ratios. The burning velocity of n-heptane has been obtained by different researchers: Davis and Law (1998, [13]) used a counterflow twin flame configuration, Huang et al. (2004, [14]) obtained the burning velocity of n-heptane and iso-octane (PRFs) in a counterflow configuration using a digital particle image velocimetry, Van Lipzig et al. (2011, [15]) investigated the adiabatic laminar burning velocity of n-heptane and iso-octane and their mixtures in a perforated plate burner, Sileghem et al. (2013, [16]) used the heat flux method to obtain the laminar burning velocities of iso-octane, n-heptane and toluene, Kwon et al. (2000, [17]) measured burning velocities in a spherical windowed chamber, Kumar et al. (2007, [18]) obtained the burning velocity of iso-octane and n-heptane with the counterflow flame technique (ambient pressure

and different temperatures and equivalence ratios), Chong et al. (2011, [19]) obtained the burning velocity using the jet-wall stagnation flame configuration and the particle imaging velocimetry and Dirremberger et al. [20] used a perforated plate burner at 1 atm and 358 K, for iso-octane, h-heptane and the toluene mixture.

In this work, the combustion process of iso-octane is studied and characterized in a constant volume combustion bomb with spherical geometry. After that, the burning velocities of n-heptane and a mixture of n-heptane and toluene are experimentally obtained at engine like conditions. This experimental facility has been used and validated before for the determination of the burning velocity of gaseous fuels and is validated in this work for liquid fuels. The values obtained for the burning velocity are expressed as power laws of the pressure, temperature and equivalence ratio. The experimentally registered pressure is the input for a two-zone thermodynamic model used to determine the burning velocity of different fuel mixtures.

This paper presents results of the burning velocities for gasoline surrogates at a high pressure and temperature, i.e., engine like conditions. In the literature there are not many results for gasoline surrogate fuels, and present results expand the burning velocity database to different conditions and mixtures including toluene for a surrogate fuel as a binary mixture with n-heptane. In addition, the burning velocity values are important for the studies developing explosion protection. Additionally, surrogate fuels are currently used as fuels in simulations for advanced combustion engines (for example HCCI, homogeneous charge compression ignition engines, or PCCI, premixed charge compression ignition engines), which means it is necessary to obtain results for the burning velocity of these fuels.

#### **2. Experimental Facility and Combustion Model**

The experimental installation used in this work consisted of a test facility designed for the characterization and investigation of the combustion process of gaseous and liquid fuels. The main part of the facility is a spherical constant volume combustion bomb (CVCB) with an acquisition system to register the parameters during the combustion and obtain information about the flame development, and supply lines for the introduction of fuels in the combustion bomb. Fuel evaporation is required for liquid fuels in order to obtain a homogeneous mixture at the beginning of the combustion process. In Figure 1, a scheme of the experimental facility can be seen, where the CVCB is a spherical space made of stainless-steel 200 mm in diameter, pressure and temperature transducers and two optical accesses which are radial and horizontal to detect and study the chemiluminescence emitted by the flame. The CVCB has been designed to withstand pressures up to 40 MPa and temperatures up to 1073 K during the course of the combustion process. There are two electrodes inside the CVCB between which the spark is discharged to start the combustion at the geometrical center of the sphere, for more details see [21].

The initial conditions of pressure, temperature and fuel/air ratio were set up at the beginning of each combustion test. Liquid fuels were directly introduced in the combustion chamber and the air was added later with the corresponding line supply. A mixture time is necessary to get a homogeneous mixture inside the CVCB. Once the combustion was initiated, a spherical flame front propagates inside the combustion bomb compressing adiabatically and burning the fresh mixture. During the combustion process development, a piezoelectric transducer, Kistler 7063 type (maximum calibration error of 0.06%), registered the evolution of the pressure. This transducer was connected to a KISTLER 5018A1000 charge amplifier (maximum calibration error of 0.3%). The output signal of the charge amplifier was recorded on a Yokogawa DL750 Scopecorder (16 bits AD converter). The estimated error of the pressure acquisition is 0.36% over the measuring range. With the acquisition system it is possible to obtain values of up to 30,000 Hz. Additional details of the experimental facility and of the use of OH and CH chemiluminescence for flame characterization can be seen in Tinaut et al. [22].

The burning velocity was determined by the two-zone combustion analysis model, in which the main input was the temporal evolution of the pressure registered during the combustion, in addition to the initial values of the fuel-composition and mass of the fuel blend.

**Figure 1.** Experimental facility.

݉ ˙ *ρ* The thermodynamic combustion model considers the division of the combustion chamber in two different zones: the burned zone (designated with a b subscript) and unburned zone (designated with a ub subscript). In each zone, there were applied conservation and ideal gas equations [23–25]. Once the combustion has been initiated by the spark plug, a flame front appears and the spherical flame moves concentrically with the vessel walls, as the unburned or fresh mixture transforms into burned mixture. The outputs of this model are the temperatures of both zones, the burned mass fraction, the flame front surface and the burning velocity, among other variables [23]. The burning velocity, Cc, was determined from the mass burning rate ( . *mb* ), the flame front surface (*Af* ) and from the unburned density (ρ*ub*) according to the next expression:

$$\mathbf{C}\_{\mathbf{c}} = \frac{\dot{m}\_b}{\rho\_{ub} A\_f} \mathbf{'} \tag{1}$$

As indicated in [26], some parts of the burning velocity plots were discarded to calculate the laminar burning velocity: some part of the initial points (due to the effects of the spark plug on the combustion onset, numerical oscillations and stretch rate) and some part of the final points (the autoignition processes and cellularity effects due to instabilities). A cellular combustion is obtained under certain conditions in which a cellular flame front ends in an apparent burning velocity higher than the laminar one when both velocities are referred to a smooth flame front—they can also be referred to the spherical one. The last part of the combustion process also had to be discarded because the flame front reached the CVCB wall and it could disturb the free flame development. The unburned temperature could also be locally perturbed due to heat transmission, and the hypothesis of adiabatic compression for the unburned gas temperature calculation will stop being fulfilled and some buoyancy effects may appear.

#### **3. Results**

Initially, the methodology for the determination of the burning velocity was presented and validated for liquid fuels. Iso-octane is the fuel used because of the existence of the extensive literature studying combustion. After that, the burning velocity of n-heptane and a mixture of 50% n-heptane/50% toluene are obtained for different conditions of the pressure, temperature (engine like conditions) and equivalence ratio, and are compared with literature data.

#### *3.1. Validation of the Methodology for Liquid Fuels: Iso-Octane*

In this section, the methodology for the determination of the burning velocity of liquid fuels is validated using iso-octane as a fuel. This methodology has been previously validated with gaseous fuels.

The influence of pressure and temperature is studied in a stoichiometric mixture of iso-octane and air for the conditions shown in Figure 2, with an initial pressure between 0.1 and 0.5 MPa and initial temperature between 360 and 450 K. In Figure 3, the pressure evolution versus the flame front radius in a combustion process of stoichiometric iso-octane is plotted. In Figure 3, it is possible to see that the pressure reaches 50% of its maximum value when the flame front is close to the wall of the CVCB (94.2 mm radius). In Figure 4, the temporal evolution of the pressure (Figure 4a, experimental values) and burning velocity (Figure 4b, obtained by means of the 2-Z model) is plotted versus unburned temperature, for stoichiometric combustions of iso-octane at different initial conditions. It can be seen in both Figures that the burning velocity reduces as the initial pressure is increased.

**Figure 2.** Pressure–temperature lines of adiabatic compression of iso-octane used for the calculation of the burning velocity for stoichiometric conditions.

**Figure 3.** Pressure increase versus the flame front radius in a combustion of iso-octane, with a stoichiometric fuel/air equivalence ratio, 0.1 MPa and 353 K initial conditions.

The range of valid data used for the burning velocity determination excludes initial data affected by the ignition process (a value of the mass fraction burned, equal to 0.05, is used as the initial validity point for the burning velocity). It also excludes the final points and data affected by instabilities and the autoignition process. Burning velocities are obtained with the two-zones thermodynamic model. The burning velocity shows a smooth increment, except for the lower pressure combustion where it is possible to appreciate a bump in the curve, caused by the instabilities of the flame front. Additionally, autoignition phenomena may arise, generating oscillations in the pressure line and a characteristic sharp bump in the burning velocity plots (as can be seen in some experiments presented in following sections).

**Figure 4.** (**a**) Pressure; (**b**) burning velocity. The temporal evolution of the pressure and burning velocity of iso-octane versus the unburned temperature for stoichiometric mixtures with varying initial pressures.

α β The burning velocity values shown in Figure 4b can be fitted to a power law correlation as a function of pressure and temperature, of the type proposed by Metghalchi–Keck [8], as seen in Equation (2), where *Cc* is the burning velocity in m/s, *Cco* is the reference burning velocity at the *To* and *po*, α and β are the exponents of temperature and pressure to be determined. The value of the correlation coefficient R2 is 0.96 and the standard error of estimation is 4%.

$$\mathcal{C}\_{\mathfrak{C}} = \mathcal{C}\_{\mathcal{O}} \left( \frac{T}{T\_o} \right)^{\alpha} \left( \frac{p}{p\_o} \right)^{\beta} \,. \tag{2}$$

$$\mathcal{C}\_{\varepsilon}(iso-octane) = 0.32 \left(\frac{T}{300}\right)^{1.91} \left(\frac{p}{0.1}\right)^{-0.25} \text{.}\tag{3}$$

ܥ ܶ ൬ 0.32) = ܽ݊݁ݐܿ െ ݏ݅) 300൰ ଵ.ଽଵ ቀ 0.1ቁ ି.ଶହ , In Table 1 the correlation obtained in the present work is presented together with the expressions obtained by Gülder [9], Bradley et al. [7], Metghalchi et al. [10], Galmiche et al. [15] and Müller et al. [11] for the stoichiometric fuel/air equivalence ratio. In Müller's expression, *To* is the temperature of the interior zone, which marks the transition between the reaction zone and the unburned zone (the temperature that limits chain brain and chain break reactions), *Yf* is the unburned fuel mass fraction and Tad is the adiabatic temperature.

**Table 1.** Correlations for the burning velocity of stoichiometric iso-octane.

ܶ ൬ 0.40 = ܥ 300൰ ଵ.ଵ ቀ 0.1ቁ ି.ଶ଼ ܶ ൬ 0.56 = ܥ 423൰ ଵ.଼ଽ ቀ 0.1ቁ ି.ଶ ܥ = 2.93 ∙ 10ଷ ቀ 3.80 ∙ 10 ቁ ൫ଵ.ଶଶ∙ଵషయ൯ ܻ .ହ ൬ ܶௗ െ ܶ ܶௗ െ ܶ ൰ ଶ.ହଶ ൬ ܶ ܶ ൰ In Figure 5a, the burning velocity of the stoichiometric mixtures of iso-octane and air is plotted versus the unburned temperature for an initial pressure of 0.5 MPa. The line obtained in this work is placed at the center of the plot, with a slope similar to the one obtained by Müller et al., in spite of the fact that the Müller expression is different and they obtain it from a program with an extensive reaction mechanism. Galmiche et al. [9] measured the flame velocities of spherically expanding flames using a combustion chamber equipped with a high-speed shadowgraph system. The differences between the expression of Bradley et al. [7] and Metghalchi et al. [8] are due to the fact that they considered all the data obtained during the combustion processes without any distinction between the laminar and cellular zone, which explains the higher slope of the Metghalchi correlation. Gülder [6] obtained the expression for the burning velocity in a combustion bomb using ionization sensors in the combustion zone, leading to higher values.

**Figure 5.** The burning velocity for stoichiometric iso-octane versus temperature (**a**) and pressure (**b**).

In Figure 4b, the burning velocity versus the pressure is represented for an initial temperature of 400 K. The tendencies of the present work results are similar to those explained before: the slope is similar to the one obtained by Müller et al. [10]. The differences with the correlations obtained by Gülder [6] and Bradley et al. [7] are due to a different measurement method being used, while the differences with the Metghalchi and Keck [8] values are due to the range of considered data.

The stoichiometric iso-octane burning velocity is represented in Figure 6 versus the temperature and pressure, for the expressions shown in Table 1. The highest values of the burning velocity are obtained for elevated temperatures and low pressures, in all the cases. In Figure 6, all the correlations are represented together in order to compare them more clearly. In addition, some experiments with varying equivalence ratios were developed to compare with the literature data, as seen in Figure 7. The values obtained in the present work (black points) are compared with the values of Galmiche et al. (2012, [9]), Metghalchi and Keck (1982, [8]), Varea et al. (2013, [11]), Marshall et al. (2011, [12]) and Gülder (1982, [6]). It can be seen that the values obtained in the present work agree with the rest of the data.

**Figure 6.** The burning velocity of iso-octane as a function of temperature and pressure for the stoichiometric equivalence ratio. Comparison between the expressions obtained by Gülder [6], Metghalchi and Keck [8], Bradley et al. [7], Galmiche et al. [9] and those obtained in the present work.

As a conclusion of this section, it is possible to say that the methodology for the determination of the burning velocity of liquid fuels in the CVCB has been validated and the results obtained for other fuels will be presented in the following sections.

**Figure 7.** Comparison of the burning velocity of iso-octane versus the equivalence ratio, for 0.8 MPa and 400 K of the initial conditions.

#### *3.2. Burning Velocity of n-Heptane*

In this section, the results for the burning velocity of n-heptane/air mixtures in the CVCB are presented for the different initial conditions of the pressure and temperature. The initial conditions are chosen to be able to reach the range of pressure and temperature obtained by the internal combustion in engines during the progress of the combustion, as can be seen in Figure 8a. The objective of these experiments is to characterize the burning velocity of n-heptane through correlations as a function of pressure and temperature and the fuel/air equivalence ratio, Equation (2).

#### 3.2.1. Results of n-Heptane for a Stoichiometric Mixture

A pressure–temperature diagram of all the data used to determine the burning velocity of n-heptane in stoichiometric conditions is shown in Figure 8a. In this figure, points represent the adiabatic evolution of the pressure and temperature during the combustion process inside the CVCB, starting from the initial conditions of pressure and temperature. Figure 8a delimits the validity range of the obtained correlation, with an initial pressure between 0.1 and 1.2 MPa and initial temperature between 370 and 650 K.

(**a**)

(**b**)

**Figure 8.** (**a**) Pressure–temperature diagram; (**b**) pressure; (**c**) burning velocity. The temporal evolution of the pressure and burning velocity of n-heptane versus the unburned temperature for stoichiometric combustions.