Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 47, 1, pp. 55-67, Warsaw 2009 INDICATED MEAN EFFECTIVE PRESSURE OSCILLATIONS IN A NATURAL GAS COMBUSTION ENGINE BY RECURRENCE PLOTS Grzegorz Litak Technical University of Lublin, Department of Applied Mechanics, Lublin, Poland; UniversitaPolitecnica delleMarche, Dipartimento di Architettura, Costruzioni e Strutture, Ancona, Italy; e-mail: g.litak@pollub.pl Arkadiusz Syta Technical University of Lublin, Department of Applied Mathematics, Lublin, Poland Bao-Feng Yao Guo-Xiu Li Beijing Jiaotong University, School of Mechanical, Electronic and Control Engineering, Beijing, China We investigate the nonlinear time series of Indicated Mean Effective Pressure (IMEP) of the spark ignition engine cyclic combustion process of a natural gas.By applying the embedding theoremand the recurrence plots technique, we show changes in the engine dynamics for different equivalence ratios. Especially, we provide arguments for intermittency behaviour. Key words: engine, combustion, recurrence plot, nonlinear oscillations 1. Introduction Theproblemof harmful cycle-to-cycle fluctuations in combustionwas the sub- ject of intensive research in the last fewdecades (Daily, 1988; Daw et al., 1998; Green et al., 1999;Hu, 1996;Kantor, 1984; Li andYao, 2008; Litak et al., 2007; Sen et al., 2008; Wagner et al., 2001; Wendeker et al., 2004). Identification of main factors influencing cycle-to-cycle combustion variations by Heywood (1988) shed somemore lights on its dynamics. These factors included aerody- namics in the cylinder during combustion, the amount of fuel, air and recycled exhaust gases supplied to the cylinder and the mixture composition near the 56 G. Litak et al. spark plug. It became clear that fluctuations can appear by a non-periodic character of cycle-to-cycle dynamics (Daily, 1988; Green et al., 1999; Kantor, 1984). Basing on the balance of recycled exhaust gases, Daw and collaborators (Daw et al., 1998) derived a simple theoretical model of a lean spark ignition internal combustion process. They supported the hypothesis that the com- bustion instability develops as a noisy period-doubling bifurcation using the equivalence ratio which as a control parameter by numerical simulations and experiments. 2. Experimental studies and obtained results In this note, we continue experimental investigations in this direction. By de- creasing the equivalence ratio frommore stoichiometric to very lean conditions in the engine fueled by a natural gas, we observe how combustion fluctuations appear. In spite of the different fuel engine to that examined by Daw et al. (1998), Green et al. (1999), Wagner et al. (2001), the combustion instabilities have the similar nature. These fluctuations seem to depend considerably on the equivalence ratio of the air-fuel mixture. To measure the pressure inside the engine, we used piezoelectric sensors. The schematic picture of our experi- mental stand is presented in Fig.1. The engine specifications can be found in Table 1. Further details on our experimental stand can be found in Li andYao (2008). After measuring pressure, we have estimated the Indicated Mean Ef- fective Pressure (IMEP) which is defined as the equivalent constant pressure in a given combustion cycle. Fig. 1. Scheme of the experimental setup: 1 – engine; 2 – dynamometer; 3 – dynamometer controller; 4 – high-speed data acquisition board; 5 – pressure transducer; 6 – optical encoder; 7 – ECU; 8 – computer Indicated mean effective pressure oscillations... 57 Table 1.Engine specifications Cylinder number 6 Bore× store 126mm× 130mm Displaced volume 9.726 L Compression ratio 10.5 Intake valve opens 2◦BTDC Intake valve closes 208◦ATDC Exhaust valve opens 227◦BTDC Exhaust valve closes 5◦ATDC The corresponding plot of IMEP against consecutive cycles i is presented inFig.2.Cases 1-4 differ by thedecreasing equivalence ratio: Φ=0.781, 0.677, 0.595, and 0.588. Fig. 2. Cycle-to-cycle changes of IMEP(i) [MPa], i enumerates the successive cycle engine. The equivalence ratio was chosen Φ=0.781, 0.677, 0.595, and 0.588 for cases 1-4, respectively Our strategy in analysing these nonlinear time series would be to recon- struct amultidimensional phase space IMEP basing on the scalar (IMEP(i)) time series. Thus for i > (M−1)∆i IMEP = (2.1) ={IMEP(i−(M−1)∆i),IMEP(i−(M−2)∆i), . . . ,IMEP(i−∆i),IMEP(i)} where ∆i and M are the time delay and the characteristic embedding dimen- sion, respectively (Boguś andMerkisz, 2005; Takens, 1981). These quantities are to be found by a detailed examination of the average mutual information (AMI) (Fraser and Swinney, 1986; Hegger et al., 1999; 58 G. Litak et al. Kantz and Schreiber, 1997) and the false nearest neighbour fraction (FNNF) (Abarbanel, 1996; Hegger et al., 1999; Kennel et al., 1992; Sen et al., 2008). AMI is defined via conditional probabilities of event sequences AMI(δi) =− ∑ kl pkl(δi)ln pkl(δi) pkpl (2.2) where, for some partition (16 equal parts) of the cyclic effective pressure, IMEP ∈ [IMEPmin,IMEPmax]. In the above formula pk is the probability to find a time series value in the k-th interval, while pkl is the joint probability that an observation falls later into the l-th interval. and the observation time is given by δi. The optimal time delay ∆i = δi is to be determined by the first AMI minimum for which the examined events are independent enough to define a new coordinate. Note, that AMI is positively defined and its smallest value (theoretically AMI =0) can be reached when pkl can be factorised to indivi- dual probabilities pk and pl (pkl ≈ pkpl) for any k and l far from each other by δi. On the other hand, to get a proper FNNF, one has to choose the point indicated by IMEPi and calculate the distance to its nearest neighbour point IMEPj in the m-dimensional space. For an Euclidean distance, which is typically used here, it is |IMEPi−IMEPj|m. By iterating both points along the time series, we compute the control parameter Qi,m defined as Qi,m = |IMEPi−IMEPj|m+1 |IMEPi−IMEPj|m (2.3) By comparing the above value to a chosen threshold Qc, we calculate the frac- tion of cases for which Qi,m exceeds the threshold value Qc. Finally, FNNF can then be estimated from the following expression FNNF(m)= 1 N ∑ i Θ(Qi,m−Qc) (2.4) where N is the number of vector elements in the vector time series, Θ(x) is the Heaviside step function. This so called fraction analysis is repeated by choosing different values of the dimension m. The optimal value M = m is defined when the fraction of false nearest neighbours tends to zero (note, in some cases, depending on the Qc value with respect to the standard square Indicated mean effective pressure oscillations... 59 deviation of the examined time series IMEP(i), some points are omitted and FNNF reaches the minimum value for the optimal dimension m=M) lim m→M FNNF(m)→ 0 (2.5) Using the above definitions for AMI (Eq. (2.2)) and FNNF (Eqs. (2.3)- (2.5)), we have estimated the embedding for the IMEP time series. Consequ- ently, in Figs.3a,b, one can easily find that the optimal values are ∆i=1 and M =5 for the optimal embedding in all of the considered cases (Fig.2). Fig. 3. Averagemutual information (AMI) and the false nearest neighbour fraction (FNNF) for all four examined cases 1-4. Basing on Eq. (2.2) and Eqs. (2.3)-(2.5), the embedding parameters (M,δi)= (5,1) These initial results on the embedding space presented above would be a natural frame for further recurrence studies: recurrence plots (RP) and recur- rence quantification analysis (RQA) techniques. Here, we assumed that these quantities are slowly evaluating during rock cutting. The recurrenceplot is usuallydefinedby the followingdistancematrix form R m,ǫ with the corresponding elements Rm,ǫij (Casdagli, 1997; Eckmann et al., 1987; Marwan, 2003, 2006; Marwan et al., 2007; Thiel et al., 2004; Webber and Zbilut, 1994; Wendeker et al., 2004) R m,ǫ ij =Θ(ǫ−‖IMEPi−IMEPj‖) (2.6) having 0 and 1 elements to be translated into the recurrence plot as an empty place and a black dot, respectively. In other words, Rm,ǫij = 1 measure the recurrence of the physical state IMEP with the tolerance ǫ. In this method (RP), we examine patterns showing diagonal and vertical or horizontal struc- tures of the lines. After obtaining such a structure, one can easily classify the dynamics of the studied system (Litak et al., 2008a,b; Marwan et al., 2007). 60 G. Litak et al. Fig. 4. Recurrence plots of IMEP for four values of the equivalence ratio Φ=0.781, 0.677, 0.595, and 0.588, (a)-(d) respectively. RR was fixed to 0.15 for all cases In Figs.4a-d, we mapped the corresponding matrix elements Rij of the investigated cases into recurrence plot graphs. For quantitative analysis, we define the recurrence rate RR RR= 1 N2 N∑ i,j=1 R m,ǫ ij for |i− j| ­w (2.7) which determines the black dots fraction in the RP graph. w denotes the Theiler window used to exclude identical and neighbour points of the time series IMEP (seeEq. (2.1) andFig.2) fromthe above summation (Eq. (2.7)). In our case, w was equal 1. Furthermore, the RQA can be used to identify vertical or diagonal lines through their lengths up to Lmax, Vmax for diagonal andvertical lines, respec- tively. In its frame, the RQA enables one to perform probability p(l) or p(v) distribution analysis of the lines according to their length l or v (for diagonal and vertical lines). In practice, they are calculated as p(y)= Pǫ(y) ∑N y=ymin Pǫ(y) (2.8) Indicated mean effective pressure oscillations... 61 where y = l or v depending on the diagonal or vertical structures in the specific recurrence plot, Pǫ(y) is a histogram of the lengths y of the diago- nal or vertical lines with the tolerance of recurrence ǫ (Eq. (2.6)). For various collections of thediagonal andvertical lineswith respect to their lengths distri- butions, Shannon information entropies (LENTR and VENTR) can be defined via (Marwan, 2003) LENTR =− N∑ l=lmin p(l)lnp(l) VENTR =− N∑ v=vmin p(v)lnp(v) (2.9) Other properties of RP as determinism DET and laminarity LAM as well as the trapping time TT are also based on the probabilities Pǫ(x) DET = ∑N l=lmin lPǫ(l) ∑ N i,j=1R m,ǫ i,j LAM = ∑N v=vmin vPǫ(v) ∑ N v=1vP ǫ(v) (2.10) TT = ∑N v=vmin vPǫ(v) ∑N v=vmin Pǫ(v) In the above equations, lmin and vmin (lmin = vmin = 2 in our case) denote minimal lengths of the diagonal and vertical lineswhich should be chosen for a specific dynamical system. The determinism quantity DET is the measure of the predictability of the examined time series and gives the ratio of recurrent points formed in the diagonals to all recurrent points. Note that in a perio- dic system all points would be included in the lines. On the other hand, the laminarity LAM is a similar measure which corresponds to points formed in the vertical lines. For small point-to-point changes (laminar), the consecutive recurrence points form a vertical line, while turbulent or chaotic changes pro- duce singular points or short lines in the vertical direction. These measures tell about dynamics behind sampling points changes and are strictly connec- ted to the points fraction spanning the diagonal (DET) and vertical (LAM) patterns, respectively. These diagonal and vertical line patterns form the base of deterministic features, while any singular point corresponds to randomness in the examined system.Note, for randomnumbers the recurence plot is filled uniformly without any patterns. Finally, the trapping time TT refers to the average length of vertical linesmeasuring the time scale (in terms of sampling intervals) of these small changes in the examined time history. We performed calculations (using the numerical code provided inMarwan (2006)) of all the specified quantities for our time series (Fig.2) and inclu- ded them into Table 2. For better clarity, we plotted the tendencies of DET 62 G. Litak et al. and LAM in Figs.5a and 5b, respectively. Interestingly, DET reaches the maximum for the smallest equivalence ratio (case 4). It is associated with the minimum in LAM and a considerable local increase in TT (Table 2). Table 2. Summary of statistical properties and the recurrence quantifica- tion analysis (RQA) of IMEP for different equivalence ratios Φ. The engine speed of 1600r/min was fixed for all cases. Note, the embedding parameters (M,δi) = (5,1). In all examined cases, ǫwas appropriately chosen to give the same recurrence rate RR=0.15 Case Φ 〈IMEP〉 σIMEP DET LAM Lmax Vmax LENTR VENTR TTNo. 1 0.781 0.829 0.0234 0.851 0.530 16 21 1.831 1.246 2.853 2 0.677 0.447 0.0276 0.864 0.564 22 36 1.929 1.590 3.482 3 0.595 0.342 0.0497 0.841 0.430 17 21 1.814 1.373 3.009 4 0.588 0.376 0.1502 0.915 0.032 21 10 2.136 1.545 3.344 Fig. 5.DET (a) and LAM (b) calculated for measured cases 1-4 (see Fig.2). In calculations, the RR parameter was fixed to the same value RR=0.15 In Figs.6a and 6b, we plot the corresponding entropies obtained from the diagonal and vertical lines statistics. The considerable increase of LENTR (ca- se 4) confirms the increasing complexity of the system response for lean com- bustion. Note that VENTR also reaches its local maximum for case 4. On the other hand, basing on the entropy results (Fig.6), case 3 could be classified as the more periodic as we observe the local minima in LENTR and VENTR. Surprisingly, Lmax =17 is relatively short for this case (case 3) comparing the neighbour cases (case 2: Lmax =22, case 4: Lmax =21). Finally, the change in Vmax is monotonous through cases 2-3-4 with decreasing tendency. Note, the similar tendency was observed for the parameter LAM . These symptoms together can be associated with the intermittency bifurcation associated with Indicated mean effective pressure oscillations... 63 the interrupted systemfluctuations of twodifferent types. In theRPdiagrams, we should be able to observe characteristic vertical patterns of line collections (Wendeker et al., 2004). Fig. 6.LENTR (a) and VENTR (b) calculated for measured cases 1-4 (see Fig.2). In calculations, the RR parameter was fixed to the same value RR=0.15 In summary,wewould like to add that by changing the equivalence ratio Φ from more stoichiometric to very lean conditions in the engine fueled by a natural gasweobservedadramaticdecrease of the laminarityparameter LAM (Fig.5a) and simultaneous increase of the determinism DET (Fig.5b). These results could indicate that going to a more lean combustion condition drives the engine to a less stable combustion. Note that this limit was also a subject of investigation by Daw and collaborators (Daw et al., 1998) who discovered chaotic oscillations of heat release in an engine fueled by the petrol fuel. The additional indication supporting this coincidencewas a relatively large value of diagonal (Fig.6a) and vertical (Fig.6b) lines entropy. Finally, we observe clear features of intermittency. Payingmore attention to Figs.4a-d, one can observe an interesting evolution of the vertical lines. Starting from the most thick, basically square structures in Fig.4a, we could see a more thin line structure in Fig.4b and Fig.4c. The results presented in Fig.4d look quantitatively different. Here one can distinguish delicate skeletons of the vertical diagonal lines and narrow vertical stripes. 3. Concluding remarks The obtained results from RQA analysis have been related to the traditional statistical measure of square deviation σIMEP (Table 2). One can see that 64 G. Litak et al. σIMEP increases monotonically with decreasing Φ but the transition from ca- se 3 to 4 is associated with an extremely large increase (about three times) of σIMEP. In the same time, the average pressure 〈IMEP〉 decreases with de- creasing Φ up to case 3, and then slightly increases in case 4. This is also a signature of intermittency. Much broader distribution of IMEP in lean com- bustion conditions must be caused by the effect of alternate less and more efficient combustion cycles. After each relatively bad combustion or misfire in the preceding combustion cycle, the fresh intake charge is mixed with residu- al gases producing a richer mixture. The richer mixture causes more efficient combustion in the current combustion cycle but their residual gases influence worse the mixture in the next cycle. This effect can be investigated by recur- rence plots by considering every second cycle. We have done such an analysis to compare cases 3 and 4. The relating plots are presented in Fig.7a and 7b. Comparing these two figures, one can observe the appearance of characteristic ’checkerboard’ patterns created by changing Φ from 0.595 to 0.588. The clear square-like shape of coloured regions along the diagonal line implies the exi- stence of type I intermittency (Klimaszewska and Żebrowski, 2007; Pomeau andManneville, 1980). Fig. 7. Recurrence plots IMEP plotted for even (every second) i and j and Φ=0.595 (a), and 0.588 (b). RR was fixed to 0.15 as in Fig.4̇, while diagonal points i= j were added here additionally Note that RP and RQA analyses provide strong arguments that relative- ly short time series can be investigated by these tools. There are of course a several drawbacks of the method which include the lack of more detailed in- formation about trends and periodics. For instance, larger empty places inRP (Figs.4a,b) informabout non-stationarities. This effect should be investigated using longer time series or/and using the other measures like the multisca- le entropy (Costa et al., 2005). Our results provide important indications to the nature of the combustion process and may be used to the improvement Indicated mean effective pressure oscillations... 65 of combustion control (Matsumoto et al., 2007). However, to tell more about system efficiency in a particular case one needs to perform more systematic studies including fuel consumption rates. Acknowledgements G.L. hasbeenpartially supportedby the 6thFrameworkProgramme,MarieCurie Actions,Transfer ofKnowledge,GrantNo.MTKD-CT-2004-014058andby thePolish Ministry of Science and Higher Education, Grant No. 9008/B/T02/2007/33. References 1. AbarbanelH.D.I., 1996,Analysis of ObservedChaotic Data, Springer,Berlin 2. Boguś P., Merkisz J., 2005,Misfire detection of locomotive diesel engine by non-linear analysis,Mechanical Systems and Signal Processing, 19, 881-899 3. Casdagli M.C., 1997, Recurrence plots revisited,Physica D, 108, 12-44 4. Costa M., Goldberger A.L., Peng C.-K., 2005, Multiscale analysis of biological signals,Phys Rev E, 89, 021906 5. Daily J.W., 1988, Cycle-to-cycle variations: a chaotic process?, Combustion Science and Technology, 57, 149-162 6. DawC.S.,KennelM.B.,FinneyC.E.A.,ConnollyF.T., 1998,Observing andmodelling dynamics in an internal combustion engine,Physical Review E, 57, 2811-2819 7. Eckmann J.-P., Kamphorst S.O., Ruelle D., 1987, Recurrence plots of dynamical systems,Europhys. Lett., 5, 973-977 8. Fraser A.M., Swinney H.L., 1986, Independent coordinates for strange at- tractors frommutual information,Phys. Rev. A, 33, 1134-1140 9. Green J.B. Jr., Daw C.S., Armfield J.S., Finney C.E.A., Wagner R.M., Drallmeier J.A., Kennel M.B., Durbetaki P., 1999, Time irre- versibility and comparison of cyclic-variability models, SAE Paper No. 1999- 01-0221 10. Hegger R., Kantz H., Schreiber T., 1999, Practical implementation of nonlinear time series methods: The TISEAN package,Chaos, 9, 413-435 11. Heywood J.B., 1988, Internal Combustion Engine Fundamentals, McGraw- Hill, NewYork 12. Hu Z., 1996, Nonlinear instabilities of combustion processes and cycle-to-cycle variations in spark-ignition engines, SAE paper No. 961197 66 G. Litak et al. 13. Kantor J.C., 1984,A dynamical instability of spark-ignited engines, Science, 224, 1233-1235 14. KantzH., SchreiberT., 1997,Non-linear Time Series Analysis, Cambridge University Press, Cambridge 15. Kennel M.B., Brown R., Abarbanel H.D.I., 1992, Determining embed- dingdimension for phase-space reconstructionusing ageometrical construction, Physical Review A, 45, 3403-3411 16. KlimaszewskaK., Żebrowski J.J., 2007,Detection of the type of intermit- tency and the associated bifurcation using characteristic patterns in recurrence plots, preprint 17. Li G.-X., Yao B.-F., 2008, Nonlinear dynamics of cycle-to-cycle combustion variations in a lean-burn natural gas engine,Applied Thermal Engineering, 28, 611-620 18. Litak G., Gajewski J., Syta A., Jonak J., 2008a, Quantitative estima- tion of the tools wear effect in a ripping head by recurrence plots, Journal of Theoretical and Applied Mechanics, 46, 521-530 19. LitakG.,Kamiński T., Czarnigowski J., ŻukowskiD.,WendekerM., 2007, Cycle-to-cycle oscillations of heat release in a spark ignition engine,Mec- canica, 42, 423-433 20. Litak G., Sawicki J.T., Kasperek R., 2008b, Cracked rotor detec- tion by recurrence plots, Nondestructive Testing and Evaluation, DOI: 10.1080/10589750802570836, in press 21. MarwanN., 2003,Encounters with Neighbours: Current Development of Con- cepts Based on Recurrence Plots and their Applications, PhD Thesis, Univer- sität Potsdam, Potsdam 22. Marwan N., 2006, Recurrence Plots Code, http://www.agnld.uni- potsdam.de/∼marwan/6.download/rp.php 23. Marwan N., Romano M.C., Thiel M., Kurths J., 2007, Recurrence plots for the analysis of complex systems,Physics Reports, 438, 237-329 24. MatsumotoK.,Tsuda I.,HosoiY., 2007,Controlling engine system: a low- dimensional dynamics in a spark ignition engine of amotorcycle,Zeitschrift für Naturforschung, 62A, 587-595 25. Pomeau Y., Manneville P., 1980, Intermittent transition to turbulence in dyssipative dynamical-systems,Communications in Mathematical Physics, 74, 189-197 26. SenA.K., LongwicR., LitakG.,GórskiK., 2008,Cycle-to-cycle pressure oscillations in a Diesel engine, Mechanical Systems and Signal Processing, 22, 362-373 Indicated mean effective pressure oscillations... 67 27. Takens F., 1981, Detecting strange attractors in turbulence,Lecture Notes in Mathematics, 898, Springer, Heidelberg, 366-381 28. Thiel M., Romano M.C., Read P.L., Kurths J., 2004, Estimation of dy- namical invariants without embedding by recurrence plots,Chaos, 14, 234-243 29. WagnerR.M., Drallmeier J.A., DawC.S., 2001,Characterization of lean combustion instability inpre-mixed charge spark ignition engines, International Journal of Engine Research, 1, 301-320 30. Webber C.L. Jr., Zbilut J.P., 1994,Dynamical assessment of physiological systems and states using recurrence plot strategies, J. App. Physiol., 76, 965- 973 31. Wendeker M., Litak G., Czarnigowski J., Szabelski K., 2004, Nonpe- riodic oscillations in a spark ignition engine, Int. J. Bifurcation and Chaos, 14, 1801-1806 Analiza wykresów rekurencyjnych dla średniego efektywnego ciśnienia indykowanego w silniku spalinowym zasilanym gazem naturalnym Streszczenie Zbadano nieliniowe przebiegi czasowe Średniego Ciśnienia Indykowanego cyklicz- nego procesu spalania w silniku spalinowym zasilanym gazem naturalnym. Stosując twierdzenie o zanurzeniu oraz technikę wykresów rekurencyjnych, pokazano istotne różnicewdynamice silnika przy zmianie składumieszanki paliwawyrażonejwartością współczynnika równowagowego.Wszczególnościpodanoargumentyprzemawiająceza występowaniem zjawiska intermitencji. Manuscript received July 9, 2008; accepted for print October 13, 2008