Giant Electroresistance in Ferroionic Tunnel Junctions

Standard Article

A real-time pressure wave model for knock prediction and control

International J of Engine Research 1–15 Ó IMechE 2019 Article reuse guidelines: sagepub.com/journals-permissions DOI: 10.1177/1468087419869161 journals.sagepub.com/home/jer

Ruixue C Li and Guoming G Zhu

Abstract This article develops a control-oriented real-time zero-dimensional knock pressure wave model for spark-ignited engines based on the reaction-based two-zone combustion model developed earlier. This knock pressure wave model is capable of predicting the in-cylinder pressure oscillations under knocking combustion in real-time and can be used for the model-based knock prediction and control. A pressure wave equation including the knock deadening behavior is proposed, simplified, and used to calculate the pressure perturbations generated by knock combustion. The boundary and initial conditions at knock onset are analyzed and the analytic solution of the pressure wave equation is obtained. The model is calibrated and validated over two engine operating conditions at knock limit. The chemical kinetic-based Arrhenius integral and the maximum amplitude of pressure oscillations are used as the evaluation methods for knock onset and intensity prediction, and the knock frequency is studied using a fast Fourier transform of the filtered incylinder pressure oscillations. The simulation results show that the proposed knock pressure wave model is able to predict the knock onset timing, frequency, and intensity accurately under two engine operating conditions. Furthermore, the knock characteristics associated with gas mixture properties at intake valve closing is analyzed based on the experimental data, and their effect to knock cycle-to-cycle variation is also studied for the proposed model.

Keywords In-cylinder pressure oscillations, knock prediction, real-time simulation, spark-ignition

Date received: 25 April 2019; accepted: 16 July 2019

Introduction Energy crisis and environmental issues are today’s biggest challenges worldwide that motivate the automotive industry to improve engine efficiency and reduce emissions. Over 60% of vehicles powered by the traditional internal combustion engines are spark-ignited (SI) and the fuel consumption of SI engines is 20%–30% higher than the diesel engines.1 As a result, improving the thermal efficiency and hence reduce the fuel consumption of SI engines is very important. However, operating engine at its knock limit turns out to be one of the main challenges in improving the thermal efficiency of SI engines, especially for turbocharged engines. Knock is an abnormal combustion due to auto-ignition of the end-gas ahead of the propagated flame front, resulting in the rapid chemical energy release with aggressive combustion in the unburned mixtures. The unevenly distributed heat release in the burned and unburned mixtures causes high-frequency shock waves and pressure oscillations in the combustion chamber, limiting further improvement of thermal efficiency and even damaging the engine. Therefore, the prediction and

control of engine knock is very important to meet the fuel economy and emission regulations. Engine knock phenomenon is extensively studied in past decades from the knock detection methods to three-dimensional (3D) modeling of the knocking combustion. Knock detection is the first step to study knock characteristics, and the pressure sensor is one of the most popular methods used to detect engine knock in the laboratory; accelerometer (also called knock sensor) is widely used in production vehicles for detecting engine knock but the mechanical vibrations caused by engine valve events and interaction among different combustion chambers add significant noises into the knock sensor, sometimes making knock detection

Department of Mechanical Engineering, Michigan State University, East Lansing, MI, USA Corresponding author: Guoming G Zhu, Department of Mechanical Engineering, Michigan State University, 1497 Engineering Research Complex, E148, East Lansing, MI 48824, USA. Email: [email protected]

2 impossible. Some new methods are proposed in the past decade. Zhu et al.2 developed an ignition coil–based ionization detection circuit for detecting the engine knock, along with ionization-based closed-loop combustion control strategy by regulating engine sparking timing to operate the engine at its knock limit and providing the best thermal efficiency. The knock onset prediction is another important topic for improving engine knock control. Two prediction methods are widely used in past decades. One is the auto-ignition delay model3 based on calculated auto-ignition delay time t defined as duration between the end of compression and knock onset time. The pressure and temperature of the unburned mixtures are required but no chemical reaction and flame dynamics are considered. The other method is the chemical kinetic method based on the Arrhenius integral (ARI) utilizing the chemical reaction rate and concentration of the species involved in the auto-ignition. Both methods have their advantages and disadvantages. The auto-ignition delay method has been improved4–6 and is widely used due to its simplicity and low computational cost. However, without considering the important chemical reaction process in the combustion chamber, the auto-ignition process cannot be predicted accurately. The ARI method is widely used in 3D computational fluid dynamic (CFD) models for predicting the knock onset but the chemical reaction steps and species involved can be fairly large (up to thousands), resulting in high computation load. It is not suitable for the model-based knock control design. Therefore, a real-time model that is able to predict knock onset accurately is highly desired. In our previous study,7 a two-zone reaction-based combustion model for SI engines was developed and this model, with a simplified two-step chemical reaction mechanism, is capable of predicting the thermodynamic status and the species concentration of unburned mixtures in real time. This article further extends the work to predict the knock onset using the chemical kinetic method for real-time simulations. Engine knock phenomenon is complex and not consistent. In-cylinder pressure oscillation information has great potential in identifying major knock characteristics8 such as knock frequency and intensity and it is widely studied. Katsumata et al.9 investigated the endgas shock waves under knocking combustion using a high-speed direct and Schlieren photography method. Kawahara and Tomita10 visualized the auto-ignition of end-gas and pressure waves in a hydrogen sparkignition engine using a high-speed camera, and furthermore, they were able to observe the cycle-to-cycle images of the auto-ignited kernel in the end-gas region. The visualization results indicate that a large amount of unburned mixtures generates strong shock waves caused by the auto-ignited kernel explosion. These optical diagnostics methods are very helpful to understand the physical process of knock combustion. Modeling in-cylinder pressure waves under knock combustion and associated numerical simulations are widely used to

International J of Engine Research 00(0) study the knock phenomenon. Yao et al.11 developed a 3D CFD model to study the propagation and reflection of pressure waves generated by the end-gas auto-ignition. They used the 3D simulation results to optimize the combustion chamber shape. Terashima and Koshi12 proposed a one-dimensional (1D) pressure wave model including detailed chemical kinetic mechanisms to study the pressure wave generated during end-gas auto-ignition to show the influence of wall temperature to the knock intensity. Richard et al.13 developed a 1D CFD model to predict the cycle-to-cycle variability in SI engines based on the analysis of large eddy simulation (LES). These 1D and 3D models are able to predict the in-cylinder pressure oscillations accurately but with the tremendous computational load. Therefore, they are good for off-line studying knock characteristics but not suitable for the model-based knock prediction and control design. To overcome the limitations of the 1D and 3D models discussed above and reduce the computational load, simplified zero-dimensional (0D) pressure wave model is required. However, there are few studies in the area of real-time pressure wave modeling. Draper14 is a pioneer of applying the acoustic wave equation to internal combustion engines to study the pressure wave characteristics in the combustion chamber but the pressure oscillation magnitude deadening behavior, caused by the energy loss and piston movement, was not considered. Note that the deadening behavior is a natural process and must be considered for accurately predicting the real-time pressure oscillations. Based on Draper’s work, Brecq and Le Corre15 proposed a pressure envelop curve to predict the pressure oscillation peaks. Although this envelop curve shows the magnitude deadening behavior, it is a curve-fitting model that is not able to predict the actual pressure oscillations. Based on the envelop curve model, an improved 3D pressure wave model is proposed by Di Gaeta et al.16 to describe the actual pressure oscillations and magnitude decay behavior. A general solution to this wave equation is provided but the 3D model is still too complex to be used for model-based knock prediction and control design. This article focuses on developing a control-oriented knock pressure wave model capable of predicting the major characteristics of knock phenomenon in real time, utilizing our earlier work7 of a two-zone reactionbased real-time combustion model. The main contribution of this article is the proposed 0D reaction-based real-time knock pressure wave model capable of predicting the knock onset timing and in-cylinder pressure wave (used to predict knock frequency and intensity) under knock combustion. Furthermore, the cycle-bycycle knock variability is also demonstrated by the proposed model. This article is organized as follows. In the next section, the reaction-based combustion model developed earlier is reviewed briefly and a 0D pressure wave equation for modeling the in-cylinder pressure oscillations is

Li and Zhu

3

derived. The boundary and initial conditions used for deriving the pressure wave equation is further discussed. Next, the ARI for predicting knock onset and knock intensity is discussed and the proposed model is calibrated using the experimental data from a fourcylinder SI engine under knocking combustion. The simulation results are presented to show the model’s capability of predicting the knock onset, knock frequency, and intensity under different engine operating conditions. The knock cycle-to-cycle variability is also analyzed using both experimental data and simulation results, and simulation results confirm its capability of predicting the cycle-by-cycle variability. Finally, the conclusions are drawn.

Reaction-based combustion model The two-zone reaction-based combustion model developed earlier can predict the engine combustion process and thermodynamics of in-cylinder mixture. The proposed knock pressure wave model is based on the outputs of the two-zone reaction-based combustion model under knocking combustion. Note that the chemical kinetics in the combustion model is very important to provide the initial conditions for the knock pressure wave model and this section provides a brief review. The reaction-based combustion model focuses on the compression, combustion, and expansion processes using a two-zone structure, see Figure 1. It is assumed that the combustion chamber is divided into two zones: the unburned zone consisting of the pre-mixed fresh air, fuel, and the trapped residual gas, and the burned (reaction) zone where the chemical reaction takes place. The two zones are separated by the flame front and both are assumed to have the hemisphere shape. It is further assumed that the mixtures in two zones are homogeneous and the reaction products stay within the reaction zone, resulting in the expansion of this zone after the combustion starts. Both zones have heat loss to the cylinder boundary and they also interact via heat and mass transfer at the interface. The end-gas auto-ignition could occur in the unburned zone after combustion is initiated, and the chemical energy released from the knock combustion and the combustion in the reaction zone results in uneven heat release distribution, which generates shock waves in the combustion chamber. This is the main cause of in-cylinder pressure oscillations. Instead of the conventional Wiebe-based combustion model, the reaction-based combustion model is based on a simplified but practical two-step chemical reaction mechanism and the specific reaction steps are shown as follows C8 H18 +

17 k1 O2 ! 8CO + 9H2 O ðStep 1Þ 2 k2

1 * CO + O2 ) k3 CO2 2

ðStep 2Þ

Figure 1. Two-zone reaction-based combustion model.

where the proportionality factors ki (i = 1, 2, 3) are the specific reaction rate constants and can be obtained by the Arrhenius function.17 Note that ki = Ai eEa, i =Ru T

ð1Þ

where Ea, i is the reaction activation energy (J/mol), Ai is the pre-exponential factor, Ea, i and Ai are constants to be calibrated, and T is the corresponding temperature of the species involved in the reaction. The first reaction step is non-reversible so the reaction rate is vs1 = k1 ½C8 H18 n1 ½O2 n2

ð2Þ

where ½C8 H18  and ½O2  are molar concentrations of species fuel and O2 , respectively; and n1 and n2 are associated reaction orders and are empirically determined. The second step chemical reaction is a reversible reaction and the net reaction rate is given by vs2 = k2 ½COn3 ½O2 n4  k3 ½CO2 n5

ð3Þ

where ½CO and ½CO2  are molar concentrations of CO and CO2 , respectively; n3 n5 are associated reaction order and are also empirically determined. In order to calculate the chemical reaction rates of these two steps, the involved molar concentrations of associated species need to be determined first. ½Xi  in this article represents the molar concentration (moles per unit volume) of species i and is defined by ½Xi  =

Ni V

ð4Þ

where i stands for different species in each zone; Ni is the associated molecular number of species i, where Ni = mi =MWi ; and V is the associated zone volume.

4

International J of Engine Research 00(0)

The molar concentration of each species ½Xi  changes during the combustion process and the associated zone volume also changes due to the piston movement. Therefore, the rate of change for the molar concentration ½Xi  is denoted by ½X_ i  (kmol/m3 s) and given by d½Xi  dðNi =VÞ = ½X_ i  = dt dt 1 dmi 1 dV =  Ni 2 VMWi dt V dt

ð5Þ

where the first and second terms in this equation describe the influence of mass and volume changes to the molar concentration, respectively. Equations (4) and (5) are applicable to both zones and the only difference is the first term for mass changes. The mass change of each species in the reaction zone is due to the chemical reaction and mass transfer over the interface of two zones. However, the mass change of the species in the unburned zone is only due to the mass transfer. Note that the mass change in the reaction zone due to the chemical reaction can be obtained based on equations (2) and (3). Based on the second law of thermodynamics, conservation of mass, conservation of energy, and the chemical reaction mechanism discussed above, the temperature rates of change for both zones can be derived and results are presented in equations (6) and (7), see detailed derivation in Men and Zhu18 T_ r =

Q_ r Vr

SðN_ i hi Þ V_ r + Ru Tr S½X_ i   S½X_ i hi  V S½Xi hi + Vr r   S½Xi  cp, i  Ru

ð6Þ T_ u =

Q_ u Vu

V_ u Vu

+ Ru Tu S½X_ i   S½X_ i hi  S½Xi hi +   S½Xi  cp, i  Ru

SðN_ i hi Þ Vu

ð7Þ

where hi and cp, i are molar enthalpy (kJ/kmol) and molar specific heat (kJ/kmol K) of species i, respectively, with respect to temperature; Q_ r (kJ/s) and Q_ u (kJ/s) are the net heat transfer rate for the reaction zone and unburned zone, respectively; and Ru (J/mol K) is the universal gas constant. Therefore, the average in-cylinder pressure can be obtained based on the ideal gas law p =

mtot Ru Tavg Vcyl MWmix

ð8Þ

where mtot is the total mass of mixtures in the chamber; Tavg and MWmix are the average temperature and molecular weight of in-cylinder mixtures, respectively. Therefore, the in-cylinder pressure under knocking combustion can be defined as pcyl = p + Dp

ð9Þ

where Dp is the pressure wave generated by the shock waves due to knocking combustion and it is zero before

Figure 2. Relationship between the two-zone reaction-based combustion model and pressure wave model.

the knock onset. A knock pressure wave model is developed in the next section. For the further development of the pressure wave model, the initial conditions at knock onset, such as the in-cylinder pressure, volumes, and their rates of the two zones, are required and they can be obtained based on the two-zone reaction-based combustion model. The in-cylinder pressure is obtained based on equation (8), and volumes and their rates of the two zones can be calculated based on the following two equations mu Ru Tu _ dVu , Vu = dt pMWu

ð10Þ

Vr = Vcyl  Vu , V_ r = V_ cyl  V_ u

ð11Þ

Vu =

where mu , Tu , and MWu are mass, temperature, and average molecular weight of the unburned zone mixture, respectively; Vr and Vcyl are the volume of the reaction zone and the cylinder, respectively. The relationship between the two-zone reactionbased combustion model and the proposed pressure wave model is shown in Figure 2.

Derivation of pressure wave model The general 3D wave equation was first applied to the internal combustion engine by Draper14 to study the incylinder pressure oscillations. Moreover, the pressure

Li and Zhu

5

wave is a physical fluid motion that can be described by the pressure wave equation in rectangular coordinates (see equation (12)) ∂2 Dp ∂2 Dp ∂2 Dp 1 ∂2 Dp + + = 2 2 2 2 2 ∂x ∂y ∂z c ∂t

ð12Þ

where Dp is the in-cylinder pressure perturbations and c is the sound of speed. Assuming that mixture in the unburned zone is ideal and adiabatic gas, the smallamplitude sound speed theory19 can be applied, where c is a constant, leading to the adiabatic sound sffiffiffiffiffiffiffiffi g p0 ð13Þ c= r0 where g is the heat capacity ratio, assumed to be a constant for the ideal gas. Note that the subscript 0 indicates that the pressure and density in equation (13) are taken at the equilibrium conditions. This 3D wave equation (12) can be written in the cylindrical coordinates to describe the pressure oscillations conveniently (see equation (14)) ∂2 Dp 1 ∂Dp 1 ∂2 Dp ∂2 Dp 1 ∂2 Dp + + + = ∂r2 r ∂r r2 ∂u2 ∂z2 c2 ∂t2

ð14Þ

Draper14 gave a general solution to this pressure wave equation to study the pressure wave frequencies and vibration modes. However, the magnitude decay of the in-cylinder pressure wave, a very important knock property, was not considered. To study the actual pressure waves under knock combustion, Di Gaeta et al.16 proposed a damped wave equation (DWE) by introducing a time-dependent dissipation term in Draper’s wave equation (see equation (15)) 2

2

2

∂ Dp 1 ∂Dp 1 ∂ Dp ∂ Dp + 2 + + 2 2 ∂r r ∂r r ∂u ∂z2  2  1 ∂ Dp ∂Dp = 2 +s c ∂t2 ∂t

ð15Þ

where s(∂Dp=∂t) is the damping term used to describe the pressure oscillations deadening behavior due to the progressive energy loss of the in-cylinder pressure wave caused by the heat transfer, friction, and piston movement. Moreover, s is the damping coefficient to be calibrated. This article focuses on developing a control-oriented real-time knock pressure wave model capable of predicting the main characteristics of engine knock, such as the knock onset timing, intensity, and frequency. A simplified pressure wave equation is proposed, assuming that the pressure wave is uniform in the circumferential and axial directions. The assumption of uniformity in the circumferential direction is due to the preassumed knock onset location at 60% of the unburned zone radius. Therefore, the in-cylinder pressure perturbation Dp becomes a function of radius

location r and time t. That is, Dp = Dp(r, t) (see equation (16))   ∂2 Dp 1 ∂Dp 1 ∂2 Dp ∂Dp = + + s ð16Þ ∂r2 r ∂r c2 ∂t2 ∂t Assuming that t and r in equation (16) are independent, the wave equation can be solved using separating variable method and Fourier series. For the method of separating variables, the solution of wave equation (16) can be written as follows ð17Þ

Dp(r, t) = R(r)T(t)

That is a product of two functions, R(r) and T(t), and each depending on only one of the two variables r and t. Differentiating equation (16) yields the following two ordinary differential equations (ODEs) T00 + sT0 + l2 T = 0

ð18Þ

and R00 +

1 0 l2 R + 2 R=0 r c

ð19Þ

where  l2 is the separation constant. Solving equations (18) and (19) leads to the general solutions of T(t) and R(r) s

TðtÞ = e 2 t ½A cos vt + B sin vt ð20Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where v = l 1  (s=2l)2 ; A and B are constants determined by the initial and boundary conditions. Equation (19) can be reduced to Bessel’s equation by setting k = l=c and s = kr. Then, the derivative of R can be obtained by the chain rule R0 = k

dR ds

and R00 = k2

d2 R ds2

ð21Þ

Substituting equation (21) into equation (19), the following Bessel’s equation yields with order n = 0 d2 R 1 dR +R=0 + ds2 s ds

ð22Þ

Solutions of equation (22) are the Bessel function J0 and Y0 of the first and second kind. However, Y0 becomes infinite at 0 and cannot be kept. Therefore, the solution of equation (22) is R(r) = J0 (s) = J0 (kr), where J0 is the Bessel function of the first kind with order 0, which will be discussed in detail in the next section.

Boundary and initial conditions Engine knock is the auto-ignition of unburned mixture in the combustion chamber, so the initial conditions for the knock event can be grouped into two parts: before and at knock onset, the pressure wave magnitude is zero and at the moment right after knock onset, the pressure wave rate is non-zero in the unburned mixture

6

International J of Engine Research 00(0)

Table 1. Root of the Bessel functions and their derivatives for order n = 0 and 1.20 m

J0 (a0, m )

J1 (a1, m )

_J0 (b0, m )

_J1 (b1, m )

1 2 3 4 5

2.4048 5.5201 8.6537 11.7915 14.9309

3.8317 7.0156 10.1735 13.3237 16.4706

3.8317 7.0156 10.1735 13.3237 16.4706

1.8412 5.3314 8.5363 11.7060 14.8636

and zero in the burned mixture, see equations (23) and (24) Dp(r, t = 0) = 0,

r 2 ½0, Rc 

ð23Þ

and ∂Dp jr, t = 0 = ∂t



0 in the reaction zone fðrÞ in the unburned zone

ð25Þ

for all t ø 0

Bessel’s equation The solution J0 to Bessel’s equation (22) has an infinite number of roots denoted by a0, m (m = 1, 2, . . . ), and its derivative J_0 also has an infinite number of roots denoted by b0, m (m = 1, 2, . . . ). The root of the Bessel functions and their derivatives for orders 0 and 1 are shown in Table 1. It indicates that the roots are not regularly spaced on the axis. Note that m is the corresponding vibration normal mode and the roots in the third and fourth columns will be used in this article. To satisfy the boundary condition (equation (25)), the pressure wave equation (17) can be derived as  ∂Dp   _ r = Rc , t = TðtÞJ0 ðkrÞ r = Rc = 0 ∂r

ð26Þ

The time function T(t) is non-zero, and therefore, the derivative of the Bessel function is J_0 ðkRc Þ = 0

ð27Þ

Since the derivative J_0 of the Bessel function with order n = 0 has an infinite number of positive roots (Table 1), the boundary condition for equation (24) can be satisfied and equation (27) leads to kRc = b0, m ,

thus,

k = km =

b0, m , Rc

Time function By satisfying the first initial condition (23), the time function T(t) provides an initial condition T(t = 0) = A = 0. Thus, the general solution of the time function T(t) is TðtÞ = Be 2 t sin vt

ð30Þ

where B is a constant and B 6¼ 0; s is the damping coefficient to be calibrated; and the wave frequency is v=(2p) (Hz), where v is computed by sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi s 2 , lm = km c v = vm = lm 1  ð31Þ 2lm Note that the adiabatic sound speed c, which is a constant, can be further derived based on the smallamplitude of sound speed theory and ideal gas law sffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi g p0 Ru  = g T ð32Þ c= r0 MWu where Ru is the universal gas constant, MWu is the average molar weight of unburned zone mixtures, and T is the average temperature of end-gas at the equilibrium condition. The sound of speed is used for calculating the pressure wave frequency. It is found experimentally that the influence of temperature T to the knock frequency is very small and can be neglected, compared with the effect of engine geometry. Furthermore, the knocking combustion occurs rapidly and the temperature change can also be neglected. Therefore, T is assumed to be constant (constant sound of speed) and becomes a calibration parameter based on the experimental knock frequency. The damping coefficient s is a bounded constant. Since parameter v in equation (31) is defined as a positive real number, the solution of the following equation should be positive 

s 1 2lm

2 ð33Þ

.0

and therefore m = 1, 2, . . . ð28Þ

where roots (b0, m ) are given in Table 1.

Note that although the knock onset location is assumed at 60% of unburned zone radius, the solution R(r) varies as a function based on radius r due to the expansion of the reaction zone during combustion.

s

ð24Þ

where Rc is the radius of the engine cylinder and f(r) is the pressure wave rate at knock onset. For the boundary conditions, it is assumed that the cylinder wall is a rigid body, so the pressure oscillation rate at the wall is zero, see equation (25) ∂Dp j =0 ∂r r = Rc , t

Therefore, the solution of ODE (equation (19)) satisfying the boundary condition (equation (25)) is   b0, m RðrÞ = J0 ðkm rÞ = J0 r , m = 1, 2, . . . ð29Þ Rc

0 \ s \ 2lm = 2km c =



2c 2c b0, m \ min b0, m Rc Rc ð34Þ

Li and Zhu

7

Based on the characteristics of the Bessel function ( min½b0, m  = b0, 1 = 3:8317), damping coefficient s is up-bounded by s \ 2cb0, 1 =Rc . Then, the corresponding general solution of ODE (equation (18)) is

R Ðc

Bm =

0

vm

R Ðc 0

Tm ðtÞ = Bm e

s2 t

ð35Þ

sin vm t

Hence, the general solution of the pressure wave equation (16) satisfying the boundary condition (25) and the first initial condition (23) is Dp(r, t) =

‘ X

Rm ðrÞTm ðtÞ

m=1 ‘ X

  b0, m s = Bm e 2 t sinðvm tÞJ0 r Rc m=1

ð36Þ

Coefficient Bm Knock is an auto-ignition of the unburned end-gas in the combustion chamber, so the pressure wave equation for knock prediction is for the unburned zone. Therefore, the first initial condition in equation (24) regarding pressure waves in the reaction zone is ignored. As a result, coefficient Bm can be determined using the second initial condition in equation (24). Differentiating equation (36) with respect to t and using initial condition (24) lead to   ‘ X b0, m ∂Dp jr, t = 0 = Bm vm J0 r = f ðrÞ ∂t Rc m=1

ð37Þ

This is the Fourier–Bessel series representing f(r) in terms of J0 (b0, m r=Rc ). And the corresponding coefficient Bm can be determined as 2   Bm = 2 2 vm Rc J1 b1, m

R ðc

rfðrÞJ0

  b0, m r dr Rc

ð38Þ

 2  R

Ðc

J0 b0, m r = rJ2 b0, m r dr 0

Rc Rc 0   Rc 2 2 = 2 J1 b1, m

ð39Þ

Therefore, term J21 (b1, m ) can be obtained by   2 b1, m = 2 Rc

R ðc

rJ20

b0, m Rc

r dr

b rJ20 R0,cm



r dr

ð41Þ

where r is in the interval of ½0, Rc .

Pressure rate at knock onset The mass conservation equation for the in-cylinder flow field can be written as ∂Dp + c2 r  ðruÞ = 0 ∂t

ð42Þ

where r is the mass density, u is the velocity vector, and Dp is the in-cylinder pressure perturbations under knock combustion. Using this equation for the unburned zone, where knock occurs, and calculating the integral over the total unburned zone volume yield ð ð ∂Dp dVu + c2 r  ðruÞdVu = 0 ð43Þ ∂t Vu

Vu

and it can be further simplified to ∂Dp c2 rV_ u = ∂t Vu

ð44Þ

Combining equations (32) and (44), the in-cylinder pressure rate at knock onset, denoted by f(r), can be obtained as follows ∂Dp 1 dVr, 0 jr, t = 0 = fðrÞ = g p0 ∂t Vu, 0 dt

ð45Þ

where g is the heat capacity ratio, p0 is the average incylinder pressure at knock onset, and Vu, 0 and Vr, 0 are the volume of the unburned and reaction zone at knock onset, respectively. All these variables can be obtained by the two-zone combustion model developed earlier.

0

where J1 (b1, m ) is the Bessel function of the first kind with order n = 1. Since the norm square of the Bessel function J0 is

J21



rfðrÞJ0

  b0, m r dr Rc

ð40Þ

0

Then, substituting equation (40) into equation (38), the coefficient Bm can be written as

Solution of pressure wave equation for knock prediction Based on these results from the previous section, the general solution of the pressure wave equation satisfying the boundary and initial conditions is in a series form, see equation (36), and each term in the equation can be calculated, where m is the vibration mode of the in-cylinder pressure wave. Therefore, equation (36) is a general solution covering all vibration modes from m = 1 to ‘. The details for different vibration modes can be found in He et al.21 In this article, the first vibration mode (m = 1) is used for modeling the knock pressure wave and the associated solution is   b0, 1 s2 t Dpðr, tÞ = B1 e sinðv1 tÞJ0 r ð46Þ Rc

8

International J of Engine Research 00(0)

where coefficient B1 is g p0 V_ r, 0 B1 = v1 Vu, 0

R Ðc 0 R Ðc 0

Table 2. Engine parameters.



b0, 1 Rc

r dr



b rJ20 R0,c1

r dr

rJ0

ð47Þ

Note that V_ r, 0 = dVr, 0 =dt, b1 is the first (minimum) root of the derivative of the Bessel function J0 , and v1 is the damped frequency. Combining equations (28) and (31), v1 can be obtained as follows sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b0, 1 2 c2 s v1 =  ð48Þ 4 Rc 2 Therefore, the in-cylinder pressure pcyl can be obtained from equation (9).

Evaluation methods of knock phenomenon Knock onset The Livengood–Wu correlation22 was widely used as an empirical method to predict the knock timing of spark-ignition engines. And based on the developed reaction-based combustion model, a chemical kineticbased method based on the ARI is developed from Livengood–Wu correlation and used in this article to predict the start of combustion (SOC) in the unburned zone, where ARI is defined as follows ðui ARI =

Ta, knk

Aknk pcyl aknk ½C8 H18 bknk ½O2 cknk e Tu ðuÞ du ð49Þ

where Aknk , aknk , bknk , cknk are auto-ignition coefficients to be calibrated based on experimental data under knocking combustion; ½C8 H18  and ½O2  are the molar concentration of fuel and oxygen, respectively; and Tu is the unburned zone temperature. Note that the molar concentration and zone temperatures are provided by the reaction-based two-zone combustion model developed earlier; Ta, knk is the activation temperature defined by Ea, knk Ru

Value

Bore Stroke Rod Compression ratio Displacement Intake valve closing (IVC) Exhaust valve opening (EVO)

86 mm 86 mm 143.6 mm 11:1 499.56 cm3 190 BTDC 156 ATDC

BTDC: before top dead center; ATDC: after top dead center.

ð50Þ

Ea, knk (J/mol) in equation (50) is the activation energy of the chemical reaction in the unburned zone, and Ru is the universal gas constant. Note that Ea, knk is a constant to be calibrated. Equation (49) can be interpreted as an integral of chemical reaction rate of the unburned mixture (endgas); and this equation also indicates that the ARI is positive and the integral is monotonically increasing. The criterion for auto-ignition is at the crank angle when ARI reaches one, that is

ð51Þ

ARI = 1

Knock intensity Engine knock properties are generally characterized in terms of knock onset timing, knock intensity, and frequency based on in-cylinder pressure waves. In reality, the knock onset timing is the most important factor and is investigated by various studies. The knock intensity and frequency are the other two important factors to describe the knock severity and physical characteristics. There are three widely used measures for knock intensity and they are maximum amplitude of pressure oscillations (MAPOs), integral of modulus of pressure gradient (IMPG), and integral of modulus of pressure oscillation (IMPO), where the MAPO method is used in the rest of this article to represent the knock intensity due to its simplicity, which is shown as follows MAPO =

uIVC

Ta, knk =

Parameter

N 1X max jDpj N 1 u0 ;u0 + z

ð52Þ

where Dp is the pressure perturbations, N is the number of pressure wave peaks, u0 is the crank angle of knock onset, and z is the knock window length.

Experiment setup and model calibration The experimental data used to calibrate and validate the proposed reaction-based 0D pressure wave model is collected from a four-cylinder, four-stroke SI engine through dynamometer experiments. The engine parameters are listed in Table 2. The test data for two typical steady-state engine operating conditions (high load at 1500 and 2000 r/min) are used for validation, as shown in Table 3, where the relative air–fuel ratio is controlled to be close to stoichiometric. Note that the ignition timing under these two operational conditions is controlled to be near the engine knock limit, making it possible to validate the pressure wave model. An A&D CAS (combustion analysis system) is used to record the in-cylinder pressure, intake manifold pressure, ignition coil dwell current, and so on. The calibrated parameters are listed in Table 4. Parameters

Li and Zhu

9

Table 3. Engine operating conditions (at knock limit). Case

1

2

IMEP (bar) Engine speed (r/min)

7.5 1500

8.23 2000

IMEP: indicated mean effective pressure.

Table 4. Calibrated parameters. Group 1 Group 2

Aknk 3e7 s 0.015

aknk 0.5e–10 T 2080 K

bknk 0.25

cknk 1.5

Ea, knk 1.8691 J/mol

listed in Table 4 can be classified into two groups. The first group consists of Aknk , aknk , bknk , cknk , and Eknk , which are related to the knock onset prediction; and  are related to the second group parameters, s and T, knock frequency. The first group parameters are calibrated based on the knock onset timing obtained from experimental data, where aknk ;Ea, knk are coefficients corresponding to the chemical reactions of knocking combustion and have low sensitivity to engine geometry and operating conditions, and they are calibrated first. Aknk has high sensitivity and is calibrated by keeping the other four parameters constant. s is used to describe the pressure wave deadening behavior and it is bounded, see equations (33) and (34), and T is the average temperature of unburned zone mixture, and was assumed to be constant. Based on equation (32), T is used to calculate the sound of speed and then, combining with equation (48), to obtain the knock frequency. As a summary, s and T are calibrated based on the knock frequency obtained from the experiment data. Note that these parameters remain unchanged under different operating conditions, which is one of the major advantages of the proposed pressure wave model and can be conveniently used for the model-based knock control.

Simulation results and discussion Generally, the knock pressure wave is dominated by its first and second vibration modes; and the first mode frequency is around 6 kHz and the second mode frequency is around 12.5 kHz. In the simulation study, the first mode frequency is investigated in this article. It is well known that the in-cylinder pressure signal contains rich information of key knock characteristics. Figure 3 shows the experimental in-cylinder pressure of the first knock cycle and its bandpass-filtered pressure signals at 1500 r/min with indicated mean effective pressure (IMEP) of 7.5 bar. The black solid line presents the unfiltered in-cylinder pressure signal under knock

Figure 3. Experimental in-cylinder pressure and filtered pressure wave of the first knock cycle at 1500 r/min with IMEP = 7.5 bar (case 1).

combustion and the red dash-dotted line presents the pressure wave obtained using a band-pass filter of 3– 10 kHz (for the first knock mode). The knock window is between knock onset and 40 crank angle degree (CAD) after top dead center (ATDC), where the data within the knock window is used to calculate MAPO. As a result, knock onset timing and intensity can be obtained using the in-cylinder pressure signal, and the knock frequency can be found by processing the filtered pressure wave using fast Fourier transform (FFT), which will be discussed later. The model’s ability of predicting the knock onset timing, knock intensity and frequency will be discussed in subsections ‘‘Knock onset and intensity’’ and ‘‘FFT analysis of incylinder pressure waves’’ by analyzing the simulated incylinder pressures. The second part of model validation is regarding the knock cycle-to-cycle variability, which is kind of random but could be related to the mixture properties at intake valve closing (IVC). For a fixed engine operating condition, the knock intensity could change cycle-bycycle and vary from low to high. Although it is difficult to predict the random knock intensity variation accurately, the correlation between knock intensity and the in-cylinder mixture properties at IVC will be studied using the simulated pressure signal from the proposed model in subsection ‘‘Cycle-to-cycle variability.’’

Knock onset and intensity In this section, the capability of estimating knock onset and intensity is demonstrated under two engine operating conditions: case 1 (1500 r/min, 7.5 bar, Figure 3) and case 2 (2000 r/min, 8.2 bar, Figure 6), where the first knock cycle pressure signals for both cases are used. The main reason to use the first knock cycle data

10

International J of Engine Research 00(0)

Figure 5. Calculated in-cylinder pressure wave at 1500 r/min with IMEP = 7.5 bar (case 1). Figure 4. Simulated ARI and its rate in the unburned zone for the first knock cycle at 1500 r/min with IMEP = 7.5 bar (case 1).

is that the in-cylinder mixture properties at IVC are relatively consistent for the first knock cycle and after knock occurs, the properties at IVC are influenced by the previous knock cycle. To validate the proposed model, first, the experimental in-cylinder pressure of case 1 is processed using a band-pass filter of 3–10 kHz, where the first knock cycle data are shown in Figure 3. The black solid-line is the recorded in-cylinder pressure of the first knock cycle, and the red dash-dotted line is the pressure wave obtained using a band-pass filter. Figure 3 shows that the knock onset timing is 21.83 CAD ATDC and the knock intensity MAPO calculated based on the filtered pressure wave within the knock window is 0.7519 bar. For predicting knock onset timing based on the developed model, equations (49) and (51) are used to simulate the unburned zone ARI and the associated rate for case 1, see Figure 4. Note that the upper plot in Figure 4 is the ARI rate in the unburned zone and the lower one is the corresponding ARI. The rate of ARI is a measure representing the knock severity. For instance, a steep increment of ARI rate indicates a heavy knock cycle. Based on the early discussion, ARI is positive and monotonically increasing (see equation (49)). The calculated ARI confirms it and increases from 0 and reaches 1 (marked with a red dot) at 21.3 CAD ATDC, indicating knock onset timing of 21.3 CAD ATDC. Note that the experimental data show a knock onset timing of 21.83 CAD ATDC, and therefore, the predicting error is about 2.43%. The initial conditions at knock onset, such as the incylinder pressure and its rate, the unburned zone volume and the volume rate of the reaction zone, are required to model the knock pressure waves. With the predicted knock onset timing, the initial conditions can be obtained via the two-zone reaction-based combustion model developed earlier, see Table 5.

Table 5. Initial conditions at knock onset for the first knocking cycle of cases 1 and 2.

Case 1 Case 2

p0 (bar)

Vu (m3)

V_ r (m3/s)

∂p ∂t jr, t = 0

34.45 32.16

1.186e–5 1.382e–5

0.03796 0.03614

1.1696e5 1.1791e5

(bar/s)

The simulated in-cylinder pressure and pressure wave are shown in Figure 5, where the predicted incylinder pressure is filtered with a band-pass filter of 3– 10 kHz (the same one used for processing the experimental data). Next, the knock intensity MAPO is calculated, and it is 0.7999 bar with a predicting error of 6.38% over the experimental data shown in Figure 3. Similarly, the proposed model is simulated at 2000 r/ min using the boundary conditions, as shown in Table 5. Case 2 experimental in-cylinder pressure of the first knock cycle is filtered using the same band-pass filter and the results are shown in Figure 6, indicating a knock onset timing of 24 CAD ATDC and the MAPO of 1.173 bar. For predicting the knock onset timing under this operating condition, the simulated ARI and its rate in the unburned zone are shown in Figure 7. The ARI increases from 0 and reaches 1 at 23.6 CAD ATDC, indicating the knock onset timing of 23.6 CAD ATDC. Compared with the experimental onset timing of 24 CAD ATDC (see Figure 6), the predicting error is only 1.67%. The simulated in-cylinder pressure and pressure wave obtained using the same band-pass filter are shown in Figure 8. The calculated MAPO is 1.1969 bar with a predicting error of 2.037%, compared with the experimental data in Figure 5. As a summary, the model’s capability of predicting the knock onset timing and intensity under two different engine operating conditions are demonstrated by

Li and Zhu

Figure 6. Experimental in-cylinder pressure and filtered pressure wave of a typical knock cycle at 2000 r/min with IMEP = 8.23 bar (case 2).

11

Figure 8. Simulated in-cylinder pressure and pressure wave with a band-pass filter of 3–10 kHz at 2000 r/min with IMEP = 8.23 bar (case 2).

Figure 9. FFT analysis of the experimental in-cylinder pressure wave at 1500 r/min, IMEP = 7.5 bar (case 1). Figure 7. Simulated ARI and its rate in the unburned zone for the first knocking cycle at 2000 r/min with IMEP = 8.23 bar (case 2).

comparing the experiment results. The maximum predicting error for knock intensity is 6.38% (case 1) and maximum prediction error for knock onset timing is 2.43% (case 1).

FFT analysis of in-cylinder pressure waves The knock pressure wave frequency is mainly determined by engine geometry and remains unchanged throughout the combustion process under different operating conditions. In general, the FFT is a wellknown method for analyzing the knock signal in the

frequency domain. To find the knock frequency, FFT is used to analyze both experimental and simulated incylinder pressure waves for cases 1 and 2. The FFT results of experimental simulated data of cases 1 and 2 are shown in Figures 9 and 10. Figures 9 and 10 indicate that the experimental knock frequency is 6.303 kHz and remains unchanged under different engine operating conditions while the predicted knock frequency for both cases (see Figures 9 and 10) is 6.310 kHz. Note that the magnitude of the pressure waves in Figures 9 and 10 is different due to the difference in knock intensity. These results match with the experimental ones, and the predicting error of knock frequency is less than 0.2%, indicating that the proposed model is able to predict knock frequency accurately.

12

International J of Engine Research 00(0)

Figure 12. MAPO of 20 consecutive engine cycles calculated based on the experimental in-cylinder pressure at 2000 r/min, IMEP = 8.23 bar (case 2, SD: standard deviation).

Figure 10. FFT analysis of the experimental in-cylinder pressure wave at 2000 r/min, IMEP = 8.23 bar (case 2).

Figure 11. Experimental in-cylinder pressure of six consecutive cycles under knocking combustion, engine operated at 1500 r/min, IMEP = 7.5 bar (case 1).

Cycle-to-cycle variability The engine knock phenomenon shows cycle-to-cycle variability even under a fixed engine operating condition. Figure 11 shows the experimental in-cylinder pressures for six consecutive engine cycles under knocking combustion when the engine is operated at 1500 r/min with IMEP of 7.5 bar (case 1 in Table 3). It is clear that

the knock intensity changes cycle-by-cycle, where cycle 1 has light knock, cycle 2 has heavy knock, cycles 3–5 have medium knock, and cycle 6 has normal combustion (without knock–baseline cycle), where the baseline cycle pressure has a small rate of rise after the ignition and the peak pressure is around 27 bar. Compared with the baseline cycle, the knock cycles show a faster pressure rate of rise and the heavy knock cycle has the most steep rate of rise. The difference in peak in-cylinder pressure between heavy knock and based-line cycles can be up to 17 bar. The experimental in-cylinder pressures under 2000 r/ min with 8.23 bar IMEP (case 2) are also analyzed. As shown in Figure 12, 20 consecutive engine cycles under knocking combustion condition are bandpass-filtered and the corresponding MAPOs are calculated. The knock intensity plot shown in Figure 12 indicates cycleto-cycle variability. Although it is kind of random but with a repeatable pattern, a high knock intensity cycle is always followed by a low knock intensity cycle for most of the cases. Note that MAPO pattern repeats but its value is different, indicating that the knock servility varies, and the MAPO can drop down to zero (no knock), after repeating the pattern a few times, see cycles 14 and 19. Furthermore, the mean value and standard deviation (SD) of knock intensity for these 20 consecutive engine cycles are 1.0009 bar and 0.5713, respectively. The cycle-to-cycle variability of knock phenomenon may be related to the in-cylinder mixture properties at IVC, which is influenced by the last combustion cycle. Therefore, the correlation between cycle-to-cycle knock variability and mixture properties at IVC will be studied in this section using the calibrated pressure wave model. In this study, the external exhaust gas recirculation (EGR) valve was closed for engine tests and hence there is no external EGR. As a result, the mixture at IVC consists of the pre-mixed fresh air, fuel, and the trapped residual gas from the last cycle. Based on our study, the mixture properties, especially the mixture temperature at IVC (TIVC ), has a significant influence on the combustion characteristics, heat transfer, and

Li and Zhu

13

Table 6. Simulation results for six engine cycles with monotonically decreasing TIVC . Cycles

TIVC (K)

TEVO (K)

MAPO (bar)

1 2 3 4 5 6

389 376 370 366 358 349

1537 1549 1558 1565 1586 1634

2.3759 1.7679 1.5773 0.9912 0.8211 0

MAPO: maximum amplitude of pressure oscillations.

gas temperature at exhaust valve opening (EVO). For a knock cycle, the knocking combustion in the unburned zone leads to high pressure rate of rise and rapid heat loss to the wall. As a result, the net heat transfer rate in the unburned zone Q_ u will be smaller than the normal combustion cycle. Based on equation (7), the unburned zone temperature at EVO (TEVO ) will also decrease. Therefore, the mixture temperature at IVC for the following cycle will decrease due to the reduced temperature of the trapped residual gas, which further affects the next combustion cycle. This could be one of the reasons for the cycle-to-cycle knock variability. To verify this assumption and study the influence of TIVC to the knock intensity and TEVO , six combustion cycles are simulated at 2000 r/min with 8.23 bar IMEP for a given set of monotonically decreasing TIVC , along with the proposed pressure wave model, to generate pressure and its wave signals. The associated results are shown in Table 6 with the highest TIVC of 389 K and lowest TIVC of 349 K. The exhaust temperate at EVO (TEVO ) is predicted by the two-zone reaction-based combustion model (see equation (7)) developed earlier. The predicted unburned zone temperature is shown in Figure 13. As shown in Table 6 and Figure 13, the first cycle starts with the highest TIVC of 389 K, and the combustion is very fast, leading to high reaction rate and unburned zone temperature, resulting in fast heat loss to the cylinder wall and low rest net heat transfer in the unburned zone (small Q_ u ). Based on equation (7), the exhaust temperature will decrease. In addition, knocking combustion could also destroy the cylinder-wall oil film and lead larger heat transfer (loss), compared with the normal combustion cycle. Based on the energy balance, increased energy (heat) lost during knock combustion results in reduced exhaust temperature. Therefore, the first cycle has the lowest exhaust temperature (TEVO ), compared with five other cycles. A little bit lower TIVC is given to the next (second) cycle, and the simulation results indicate that TEVO increases by 12 K and the corresponding knock intensity (MAPO) decreases by 25.5%. This trend continues as the temperature at IVC reduces. Therefore, the simulation results confirm the hypothesis that TIVC is highly correlated to knock intensity. That is, the high the TIVC is, the high the knock intensity is.

Figure 13. Unburned zone temperature (six cycles) with monotonically decreasing TIVC at 2000 r/min, IMEP = 8.23 bar (case 2).

Figure 14. Curve fitting for predicting the mixture temperature at IVC of next cycle based on the exhaust temperature at the current cycle.

In addition, for the engine cycles under knock combustion at a fixed operating condition, the mixture temperature at IVC of cycle k + 1 (TIVC (k + 1)) could be influenced by the residual gas (exhaust) temperature of cycle k (TEVO (k)). Generally, the knock variability could be caused by many factors, such as the mixture properties at IVC, hot spots, the concentration of fuel, air, and residual gas throughout the chamber, locations of auto-ignition. Among them, TIVC is one of the main factors. In order to study the cycle-to-cycle knock variability due to TIVC variation, the experimental incylinder pressures of 50 consecutive engine cycles and the proposed pressure wave model are used to predict the next cycle TIVC (k + 1) based on the exhaust temperature TEVO (k) of current cycle. Assuming that TIVC is only affected by the exhaust temperature of the previous cycle, 50 engine cycles are simulated to generate in-

14

International J of Engine Research 00(0)

Table 7. Simulation results of 20 consecutive cycles using predicted TIVC (case 2). Cycles

TIVC (K)

TEVO (K)

MAPO (bar)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

348 390.25 357.49 370.09 351.86 372.19 350.64 371.49 357.87 363.20 347.87 387.67 355.87 360.18 356.91 351.52 368.49 347.52 378.51 350.47

1634 1536 1588 1557 1593 1554 1604 1555 1587 1570 1635 1540 1591 1576 1590 1594 1562 1635 1546 1608

0 2.4512 0.7986 1.5429 0.4073 1.7591 0.3985 1.6881 0.8043 0.9509 0 2.3088 0.7139 0.8915 0.7453 0.4014 1.3589 0 1.9567 0.3052

Figure 15. Predicted MAPOs of 20 consecutive engine cycles.

MAPO: maximum amplitude of pressure oscillation.

cylinder pressure and its knock wave. For each cycle, TIVC is optimized by minimizing the prediction error of in-cylinder pressure. The results are shown in Figure 14, where each blue-square marker represents one cycle experimental data point and the error bar presents the influence of other factors. Note that there are many repeated temperate points TIVC (k + 1) and the number of markers in Figure 14 is less than 50. As shown in Figure 14, the experimental TIVC (k + 1) and TEVO (k) show an obvious correlation, and the fitted curve (red solid line) can be used to predict the TIVC of next cycle based on the exhaust temperature of current cycle. However, the fitted curve only reflects the influence of TIVC due to the last cycle exhaust temperature. In order to have an accurate prediction of TIVC , the influence of other factors is included by adding a random variable to this fitted curve, shown as follows TIVCðk + 1Þ = fcf ðTEVO ðkÞÞ + c

ð53Þ

where fcf is the fitted curve (red solid line) in Figure 14, and c is a zero-mean random variable with normal distribution, that is, c;N (0, 31). Therefore, the correlation between knock intensity and mixture temperature at IVC is further studied based on equation (53), where the simulated current cycle exhaust temperature is used to find the mixture temperature at IVC for the next cycle. For the simulation study, the mixture temperature at IVC for the first cycle is assigned to 348 K (relative low) and the predicted exhaust temperature (using the proposed pressure wave model) is 1634 K and note that this cycle is not knocking with MAPO of 0. Based on equation (53), the predicted TIVC of the second cycle is 390.25 K. Therefore, the second cycle simulation is based on

TIVC (k = 2) of 390.25 K. This process continues for 20 engine cycles, see Table 7 for calculated TIVC and simulated MAPO results. As shown in Table 7 and Figure 15, it is obvious that the heavy knock cycle (high knock intensity) is followed by a light knock cycle (low knock intensity) for most cases. The simulation results match the experimental results (shown in Figure 12) in terms of knock intensity mean and SD, where the mean and SD from simulations are 0.9741 and 0.7605, respectively, and these from experiments (shown in Figure 12) are 1.0009 and 0.5713, respectively. The capability of the proposed pressure wave model in predicting cycle-to-cycle variability based on in-cylinder mixture temperature at IVC is demonstrated. We believe that the difference could be caused by the un-modeled mixture characteristics at IVC. In summary, the simulation results with the given monotonically decreasing TIVC confirm that the high TIVC is, the high knock intensity is. Also, heavy knocking combustion leads to low exhaust temperature, and hence, low mixture temperature at IVC for the next cycle. The simulation results of 20 consecutive engine cycles show that current cycle TIVC is affected by the last cycle knock intensity and variation of TIVC is one of the main causes for cycle-to-cycle knock variability at a fixed engine operating condition.

Conclusion A control-oriented (0D) knock pressure wave model, based on outputs of a two-zone reaction-based combustion model, for SI engines is developed to predict knocking combustion in this article. The developed model is calibrated and validated by experimental data. The simulation results confirm model’s capability of predicting the key knock characteristics such as knock onset timing, knock frequency and intensity. The maximum prediction error under two engine operating conditions is less than 2.5% and the maximum prediction error for knock intensity (MAPO) is less than 6.38%; the knock frequency prediction is even small (less than 0.2%). In addition, the capability of predicting the cycle-to-cycle knock variability is studied by correlating the knock intensity to the in-cylinder mixture temperature at IVC. Simulation results confirm the hypothesis

Li and Zhu

15

that the in-cylinder mixture temperature at IVC is one of the key factors leading to cycle-to-cycle knock variability, and the proposed model is able to replicate this phenomenon, where a high knock intensity cycle is often followed by a low knock intensity cycle. Declaration of conflicting interests The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article. Funding The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was partially supported by Ford University Research Program under funding number 2015-8032R. ORCID iD Guoming G Zhu

https://orcid.org/0000-0002-2101-2698

References 1. Wang Z, Liu H and Reitz RD. Knocking combustion in spark-ignition engines. Prog Energ Combust 2017; 61: 78–112. 2. Zhu GG, Haskara I and Winkelman J. Closed-loop ignition timing control for SI engines using ionization current feedback. IEEE T Contr Syst T 2007; 15(3): 416– 427. 3. Pipitone E, Beccari S and Genchi G. A refined model for knock onset prediction in spark ignition engines fueled with mixtures of gasoline and propane. J Eng Gas Turbines Power 2015; 137(11): 111501. 4. Douaud A and Eyzat P. Four-octane-number method for predicting the anti-knock behavior of fuels and engines. SAE technical paper 780080, 1978. 5. Elmqvist C, Lindstro¨m F, A˚ngstro¨m H, Grandin B and Kalghatgi G. Optimizing engine concepts by using a simple model for knock prediction. SAE technical paper 2003-01-3123, 2003. 6. Hoepke B, Jannsen S, Kasseris E and Cheng W. EGR effects on boosted SI engine operation and knock integral correlation. SAE Int J Engines 2012; 5(2): 547–559. 7. Li RC, Zhu GG and Men Y. A two-zone reaction-based combustion model for a spark-ignition engine. Int J Engine Res. Epub ahead of print 10 April 2019. DOI: 10.1177/1468087419841746. 8. Shu G, Pan J and Wei H. Analysis of onset and severity of knock in SI engine based on in-cylinder pressure oscillations. Appl Therm Eng 2013; 51(1–2): 1297–1306.

9. Katsumata M, Morikawa K and Tanabe M. Behavior of shock wave and pressure wave of SI knocking with super rapid compression machine. SAE technical paper 201101-1875, 2011. 10. Kawahara N and Tomita E. Visualization of autoignition and pressure wave during knocking in a hydrogen spark-ignition engine. Int J Hydrogen Energ 2009; 34(7): 3156–3163. 11. Yao A, Xu H and Yao C. Analysis of pressure waves in the cone-type combustion chamber under SI engine knock. Energ Convers Manage 2015; 96: 146–158. 12. Terashima H and Koshi M. Mechanisms of strong pressure wave generation in end-gas autoignition during knocking combustion. Combust Flame 2015; 162(5): 1944–1956. 13. Richard S, Dulbecco A, Angelberger C and Truffin K. Invited review: development of a one-dimensional computational fluid dynamics modeling approach to predict cycle-to-cycle variability in spark-ignition engines based on physical understanding acquired from large-eddy simulation. Int J Engine Res 2014; 16(3): 379–402. 14. Draper CS. The physical effects of detonation in a closed cylindrical chamber. Report no. 493, 1 January 1935, https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/ 19930091567.pdf 15. Brecq G and Le Corre O. Modeling of in-cylinder pressure oscillations under knocking conditions: introduction to pressure envelope curve. SAE technical paper 2005-011126, 2005. 16. Di Gaeta A, Giglio V, Police G and Rispoli N. Modeling of in-cylinder pressure oscillations under knocking conditions: a general approach based on the damped wave equation. Fuel 2013; 104: 230–243. 17. Law CK. Combustion physics. Cambridge: Cambridge University Press, 2010. 18. Men Y and Zhu GG. A multi-zone reaction-based diesel combustion model for model-based control. In: Proceedings of the ASME 2017 dynamic systems and control conference, Tysons, VA, 11–13 October 2017, p.V003T27A003, https://proceedings.asmedigitalcollection.asme.org/proceeding.aspx?articleid=2663600 19. Le Me´haute´ B. Linear small amplitude wave theories. In: Le Me´haute´ B (ed.) An introduction to hydrodynamics & water waves. Berlin; Heidelberg: Springer, 1976, pp.198–212. 20. Abramowitz M, Stegun IA and Romer RH. Handbook of mathematical functions with formulas, graphs, and mathematical tables. Am J Phys 1988; 56: 958. 21. He X, Qi Y, Wang Z, Wang J, Shuai S and Tao L. Visualization of the mode shapes of pressure oscillation in a cylindrical cavity. Combust Sci Technol 2015; 187(10): 1610–1619. 22. Livengood J and Wu P. Correlation of autoignition phenomena in internal combustion engines and rapid compression machines. Symp Combust 1955; 5(1): 347–356.