DF Individual torque estimation utilizing one in-cylinder presssure sensor and flywheel acceleration Master’s thesis in Systems, Control and Mechatronics ERIK WILLIAMSSON UGO MAURIS Department of Electrical Engineering CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2019 Master’s thesis 2019:EENX30 Individual torque estimation utilizing one in-cylinder pressure sensor and flywheel acceleration ERIK WILLIAMSSON UGO MAURIS DF Department of Electrical Engineering Chalmers University of Technology Gothenburg, Sweden 2019 Individual torque estimation utilizing one in-cylinder pressure sensor and flywheel acceleration ERIK WILLIAMSSON UGO MAURIS © ERIK WILLIAMSSON & UGO MAURIS, 2019. Supervisors: Johan Engbom, Anton Kjellin, Gabriella Lygnestrand, Volvo Technol- ogy AB Examiner: Torsten Wik, Division of Signals and Systems, Chalmers University of Technology Master’s Thesis 2019:EENX30 Department of Electrical Engineering Division of Systems and Control Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Individual combustion torques calculated from pressure sensor data, coming from six individual sensors, during a full combustion cycle (720◦) Typeset in LATEX, template by David Frisk Printed by Chalmers Reproservice Gothenburg, Sweden 2019 iii Individual torque estimation utilizing one in-cylinder pressure sensor and flywheel acceleration ERIK WILLIAMSSON UGO MAURIS Department of Electrical Engineering Chalmers University of Technology Abstract The diesel engine, although studied for years, is still subject to modifications and improvements to limit its emissions and adapt to new regulations. An important axis of development consists in estimating the torque contribution for each cylin- der using a virtual sensor based solely on the pressure data of one cylinder and the angular acceleration of the crankshaft. Access to such information could improve en- gine performance by optimizing the amount of diesel injected into each combustion chamber, leveling out the well-known differences that can occur between cylinders. To do so, this thesis proposes three methods to estimate combustion torque parame- ters, one based on estimating the Indicated Mean Effective Pressure (IMEP) giving some knowledge about mean torque over each cylinder cycle, the other one based on identifying transfer functions between the flywheel acceleration and torques at stake, catching the dynamics of the torque curves. The third method employs a discrete state space system with a Kalman filter. It is demonstrated that the first method, if well calibrated, can achieve an IMEP accuracy of more than 95 percent with a low computational load, but lacks in robustness. Whereas the second method, more robust, can potentially catch peak cylinder torque with a high accuracy but is com- plex to implement in an EMS and requires more computations. The third method, utilizing a Kalman filter, can estimate total combustion torque in simulation, but also, individual torques by employing a separation approach. Keywords: Combustion torque, IMEP, Fourier series, Transfer functions, Flywheel acceleration, Torque estimation, ICPS, Individual cylinder torque, In-cylinder pres- sure, Flywheel position sensor, Separation algorithm, Kalman filter. iv Acknowledgements We are very grateful for all the people that has helped us during this master’s thesis project. A special appreciation goes to Anton Kjellin and Grabriella Lygnestrand as well as Johan Engbom at Volvo Technology AB, who have supported us on a regular basis. A very big thank you goes to Torsten Wik who has been the examiner at Chalmers, and has been a great support during this thesis. Also, a thank you goes to Michael Di-Loreto, teacher at INSA Lyon, for giving us good input and advice on some solutions and literature. Erik Williamsson & Ugo Mauris, Gothenburg, June 2019 vi Contents List of Figures x List of Tables xiii Abbreviations xiv Acronyms xvii 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.3 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.4 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.5 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 System overview 4 2.1 Combustion cycle in a four-stroke engine . . . . . . . . . . . . . . . . 4 2.2 Pressure sensor for in-cylinder measurements . . . . . . . . . . . . . . 5 2.3 Angular sensor at the flywheel . . . . . . . . . . . . . . . . . . . . . . 6 2.4 Cylinder ID . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.5 Signal processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3 Physical models 8 3.1 Torque model with crankshaft . . . . . . . . . . . . . . . . . . . . . . 8 3.1.1 Combustion torque Tcomb . . . . . . . . . . . . . . . . . . . . . 8 3.1.2 Mass torque Tm . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.1.3 Friction torque Tf . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.1.4 Driveline torque TD . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 Powertrain model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3 Linear IMEP model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.4 Torque model utilizing a Kalman filter . . . . . . . . . . . . . . . . . 14 4 Model identification 17 4.1 Identification of Fourier series . . . . . . . . . . . . . . . . . . . . . . 17 4.2 Powertrain model identification . . . . . . . . . . . . . . . . . . . . . 23 4.3 Linear IMEP model identification . . . . . . . . . . . . . . . . . . . . 26 viii Contents 5 Torque estimation 30 5.1 Torque estimation using Powertrain model . . . . . . . . . . . . . . . 30 5.1.1 Sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . . 34 5.2 Estimation using Linear IMEP model . . . . . . . . . . . . . . . . . . 35 5.3 Torque estimation utilizing a Kalman filter . . . . . . . . . . . . . . . 38 6 Implementation on engine 40 6.1 Experimental set-up and implementation . . . . . . . . . . . . . . . . 40 6.2 Online estimation using the IMEP model . . . . . . . . . . . . . . . . 42 7 Discussion 44 7.1 IMEP estimator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 7.2 Estimation using Fourier series . . . . . . . . . . . . . . . . . . . . . 45 7.3 Kalman filter approach . . . . . . . . . . . . . . . . . . . . . . . . . . 45 8 Conclusion 46 8.1 The thesis outcome . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 8.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Bibliography 48 A Appendix 1 I A.1 Mass torque derivation . . . . . . . . . . . . . . . . . . . . . . . . . . II A.2 Engine data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III A.3 Correlation coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . IV ix List of Figures 1.1 Framework for this thesis work. . . . . . . . . . . . . . . . . . . . . . 3 2.1 The four different combustion phases in a four-stroke engine [1]. . . . 4 2.2 A dynamic pressure sensor capable of doing high precision measure- ments during a combustion cycle [2]. . . . . . . . . . . . . . . . . . . 5 2.3 Interpolated pressure signal in the power phase of the combustion. The dotted line marks where the TDC-angle is. . . . . . . . . . . . . 5 2.4 Representation of the experimental setup along the flywheel. . . . . . 6 2.5 Low-pass filtered flywheel velocity (solid line), with ωL = 36ω0, and raw signal (dotted line) versus crank angle. The window is for one full combustion cycle [720◦,1440◦]. . . . . . . . . . . . . . . . . . . . . 7 3.1 The crankshaft piston model: s is the piston’s displacement from the TDC, β is the connecting rod angle, l the connecting rod length, θ the crank angle and r is the crankshaft web length. . . . . . . . . . . 9 3.2 Powertrain with cylinder torques, crankshaft, flywheel and the driv- eline torque TD [3]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3 PV diagram of a complete combustion cycle for an engine running at 1000 rpm and a nominal torque of 1800 Nm. . . . . . . . . . . . . . 14 4.1 Fitted least squares θ0 = 2ω0t+K with tooth time data θ(tj) for each crank angle sample. The engine frequency is the slope in the fit. . . 18 4.2 Fitted Fourier series of ψ(t) in a steady state engine condition. . . . . 19 4.3 Fitted Fourier series with N = 6 to the acceleration signal θ̈(t). The root mean square error is erms = 262.52 rad/s2. . . . . . . . . . . . . 20 4.4 Fitted Fourier series with N = 24 to the acceleration signal θ̈(t). The root mean square error is erms = 97.35 rad/s2. . . . . . . . . . . . . . 21 4.5 Fitted Fourier series with N = 6 of T = Tcomb + Tm for the first cylinder (cylinder 1). The crank angle range is [720◦,1440◦], one full combustion cycle. Root mean square error erms = 530.68 Nm. . . . . 21 4.6 Fitted Fourier series with N = 24 of T = Tcomb + Tm for the first cylinder (cylinder 1). The crank angle range is [720◦,1440◦], one full combustion cycle. Root mean square error erms = 33.36 Nm. . . . . . 22 4.7 The driveline torque TD with its corresponding fitted Fourier series (dotted line), for N = 6, with an error erms ≈ 64 Nm. . . . . . . . . . 22 4.8 The driveline torque TD with its corresponding fitted Fourier series (dotted line), for N = 24, with an error erms ≈ 41 Nm. . . . . . . . . 23 x List of Figures 4.9 The fitted acceleration output from the powertrain model and the measured acceleration, for order N = 14. . . . . . . . . . . . . . . . . 25 4.10 Flywheel velocity and torque produced before the flywheel on the same window. Cylinder 1 is indicated with a cross at its maximum torque, firing order follows. . . . . . . . . . . . . . . . . . . . . . . . 26 4.11 The flywheel velocity during a smaller range, sampled as in the real engine (every 6 degrees). The window here is (8,1). . . . . . . . . . . 27 4.12 IMEP for four operating points at the same engine speed (1000 RPM) but different nominal torques. The x-axis is expressed in velocity change (velocity difference between 2 points) and not acceleration, time difference being almost constant. . . . . . . . . . . . . . . . . . . 28 4.13 Correlation coefficient along the different operating points for each configuration. Matlab correlation function was used for this curve. . . 29 5.1 Schematic overview on how the separation strategy is made. This example is estimating the coefficients in the second cylinder, Ĉ∗n,2. . . 32 5.2 Cylinder separated torques, with no faults, in firing order T1 (green), T5, T3, T6, T2 and T4. The dotted lines are the estimated torques from the separation algorithm. This estimate has an error of erms ≈ 138 Nm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 5.3 Cylinder separated torques, with −10% fuel in cylinder 3, in firing order T1 (green), T5, T3, T6, T2 and T4. The dotted lines are the estimated torques from the separation algorithm. As one can see, cylinder 3 has less peak torque that the estimator catches. This esti- mate has an error of erms ≈ 131 Nm. . . . . . . . . . . . . . . . . . . 33 5.4 Cylinder separated torques, with −10% fuel in cylinder 3, in firing order T1 (green), T5, T3, T6, T2 and T4. The dotted lines are the estimated torques from the separation algorithm. As one can see, cylinder 3 has less peak torque that the estimator catches. This esti- mate has an error of erms ≈ 139 Nm. . . . . . . . . . . . . . . . . . . 34 5.5 Sensitivity to offsets in gain. . . . . . . . . . . . . . . . . . . . . . . . 34 5.6 Sensitivity to offsets in phase. . . . . . . . . . . . . . . . . . . . . . . 35 5.7 Real and estimated averaged IMEP for cylinder 2 during a simulated run at OP 9 to 12. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5.8 Flowchart showing the estimator principle that will be used for Tar- getLink implementation. . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.9 The estimated state x̂2(n) by the Kalman state observer (dotted) and the measured state x2,measured coming from engine data (solid). . . . . 38 5.10 The estimated torque by the Kalman state observer (dotted) and the measured torque coming from the pressure sensor (solid), for combus- tion torque in cylinder 1. . . . . . . . . . . . . . . . . . . . . . . . . . 39 6.1 Calibration set-up for tuning the ADC-gain. . . . . . . . . . . . . . . 41 6.2 Experimental set-up for doing online estimation with the IMEP esti- mator running on the EMS. . . . . . . . . . . . . . . . . . . . . . . . 41 xi List of Figures 6.3 Reference combustion torque and estimated one thanks to the IMEP estimator during a run from OP 13 to 16 with an EMA coefficient of 0.05. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 6.4 Estimated IMEP for all cylinders from OP 13 to 16. . . . . . . . . . . 43 xii List of Tables 5.1 Average relative error [%] between real and estimated averaged IMEP EMA-averaged for a 4 seconds simulation time. . . . . . . . . . . . . 36 A.1 Engine data from Lamella 3. . . . . . . . . . . . . . . . . . . . . . . . III A.2 Engine data from Lamella 3 with faulty injections. . . . . . . . . . . . IV A.3 Correlation coefficients to choose the windows for the IMEP model. . V xiii Abbrevations Small letters b Least squares system matrix cf Friction coefficient erms Root mean square error h0 Tooth time sampling interval [ s ] j Index for time or Complex number l Connecting rod length [ m ] mosc Mass of the oscillating part (piston) [ kg ] mrot Total mass from the rotational part (both piston and crankshaft) [ kg ] n Number of teeth ncyl Number of cylinders (6) nr Number of revolutions per cylinder during a cycle p Pressure [ Pa ] pi Pressure in the cylinder i [ Pa ] p0 Atmospheric pressure [ Pa ] r Crankshaft web length [ m ] s Piston’s displacement from the TDC position [ m ] t Time [ s ] tj Time stamps [ s ] ttimes Tooth time vector [ s ] ∆txy Time between the x and y positions [ s ] v Derivative of s w.r.t θ [m/rad] v̄ Number of pulses of the flywheel sensor per second [ Hz ] xiv Capital letters A Least squares system matrix Apiston Area of the piston surface [ m2 ] Ax,n Real Fourier coefficient Bx,n Real Fourier coefficient Cconv Proportionality coefficient between average torque and IMEP Cx,n Complex Fourier coefficient Fcomb,i Force given by the combustion of cylinder i [ N ] FO Firing order GD Driveline’s transfer function [ kg.m2 ] Gi Individual cylinder’s transfer function [ kg.m2 ] IMEPi IMEP for cylinder i [ Pa ] J Lumped inertia from the flywheel and connections [ kg.m2 ] Ji Equivalent inertia of cylinder i [ kg.m2 ] JI Inertia of the input-shaft and clutch [ kg.m2 ] K Constant [ rad ] L Crank lever function [ m3 ] N Order of the Fourier series P Power during a cycle [ W ] Pcomb,i Power given by the combustion of cylinder i [ W ] P720 Pressure vector containing 120 values (720 degrees) [ Pa ] T Resulting torque on the crankshaft [ Nm ] Tcomb,i Combustion torque given by the combustion of cylinder i [ Nm ] Tcomb,tot Combustion torque from all cylinders [ Nm ] T̄comb,cycle Averaged total combustion torque per cycle [ Nm ] TD Driveline torque [ Nm ] Tf Friction torque [ Nm ] Tf,tot Total friction torque for all the cylinders [ Nm ] Ti Torque signal for each individual cylinder [ Nm ] TI Input-shaft torque [ Nm ] Tm Mass torque [ Nm ] Tm,tot Mass torque from all cylinders [ Nm ] Ttot Total resulting torque before the flywheel [ Nm ] Ttot,cont Total counteracting torque [ Nm ] Ttot,cont,i Equivalent total counteracting torque from cylinder i [ Nm ] T720 Averaged total combustion torque per cycle [ Nm ] V Volume of the chamber [ m3 ] Vd Volume swept by the piston [ m3 ] Wcycle Work done during a whole cycle (720 degrees) [ J ] xv Greek letters α Exponential moving average coefficient αi Identified combustion torque coefficient [ Nm ] β Connecting rod angle [ rad ] βi Connecting rod angle of cylinder i [ rad ] γi Identified combustion torque coefficient [ Nm ] ζ Shift phase [ rad ] θ Crankshaft angle relative to TDC [ rad ] |θ| Total length of the whole engine cycle [ rad ] ∆θ Sample interval (6 degrees) [ rad ] θF Flywheel angular position [ rad ] θi Crankshaft angle relative to cylinder’s i TDC [ rad ] θI Input-shaft angular position [ rad ] θ̇x Crankshaft angular velocity at point x [ rad/s ] θ0 Angular function [ rad ] χ Polynomial function ρX,Y Pearson’s correlation coefficient τi Constant shift time [ s ] φ Real coefficient ψ Zero mean angular position (oscillation function) [ rad ] ω Engine angular velocity [ rad/s ] ω̄ Average engine angular velocity [ rad/s ] ω0 Fundamental frequency of oscillations [ rad/s ] ωL Passband frequency filtering flywheel velocity [ rad/s ] xvi Acronyms BDC Bottom Dead Center. 4 CA50 Crank Angle position where 50% of the heat is released. 12 CAD Crank-angle Degree. 6 ECU Engine Control Unit. 2 EMA Exponential Moving Average. 35 EMS Engine Management System. 40 ICPS In-Cylinder Pressure Sensor. 1 IMEP Indicated Mean Effective Pressure. iv, 1 PCP Peak Cylinder Pressure. 1 TDC Top Dead Center. 4 xvii 1 Introduction In this chapter the background of the thesis is presented, as well as the purpose. Then a brief definition of the problem formulation is given, together with scope of the work. 1.1 Background Today’s heavy-duty diesel engines have stringent requirements on emission levels, such as the most recent regulation, Euro-VI [4]. This creates more restrictive oper- ating conditions for the combustion process. Therefore, new methods are being de- veloped to make the combustion process more efficient, with the help of In-Cylinder Pressure Sensor (ICPS). One way to both reduce energy and emission levels, is to use precise engine control [5]. Especially, Indicated Mean Effective Pressure (IMEP), which is proportional to the mean combustion torque, is one way to control the combustion process. This is realized by installing one or multiple pressure sensors in the cylinders of the engine. To improve a closed-loop engine control, one would want to extract precise information about torque and pressure in all the cylinders. However, the computational load could be high [6], and the cost for the pressure sensors is high, which makes the task challenging. It is therefore of high interest to develop a method to extract torque and pressure data in multiple cylinders by only using the output of one sensor combined with a flywheel position sensor. 1.2 Purpose The purpose of this thesis work is to estimate the torque contribution of all six cylinders in a 13 liter, Euro-VI diesel engine, by only utilizing the output of one in-cylinder pressure sensor and a crankshaft position sensor. This is to potentially lower the cost by not having to install six pressure sensors in trucks for commercial use, and also reduce emissions. Both estimation of average combustion torque and torque dynamics, during the combustion process, are of interest. They could for future research and work, be converted into pressure curves. Information about Peak Cylinder Pressure (PCP) and heat release rate could then be extracted to be able to control the self-ignition process to reduce certain types of emissions [5]. 1 1. Introduction 1.3 Problem formulation The problem to solve will rely on a two-fold approach. The first one is to make an estimation of the mean combustion torque using the relation between IMEP for the cylinder pressures and crankshaft acceleration, mainly inspired by [6]. This, in order to have a low computational method for an Engine Control Unit (ECU) that is suited for implementation. The second problem is to use a grey-box model to extract the torque functions, by using cylinder pressures and crankshaft acceleration, inspired by the model in [7]. The two approaches complement each other well, since in the first problem it is only possible to extract mean torque, but in the second, it is only possible to extract dynamical behaviours. Also, a third method utilizing a Kalman filter could be investigated by simulation. This could be a suitable case, since the system and measurements are expected to carry noise. 1.4 Methodology To solve the problem raised in Section 1.3 the framework in Figure 1.1 is used. First, a physical model with its assumptions is developed. Secondly, the model is identified on previous engine data as in [3]. Thirdly, the model is verified and validated by testing it on different operating points for the same engine and with the same parameters. Lastly, if simulation results are promising and implementation barriers reasonable, the method will be tested on a real engine in a test bench. Simulations are done using Matlab and Simulink. For implementation the software is written using a library called Targetlink (part of Simulink), that has an embedded coder for compiling C-code, which can be run on the EMS. 2 1. Introduction Figure 1.1: Framework for this thesis work. 1.5 Scope Only stationary operating points for the engine will be considered. This means that the average engine velocity and indicated torque are held constant. No estimation is being considered at transient points, but it is not to be ruled out that the derived model includes the possibility of working in transient operations as well. 3 2 System overview In the following sections the main system overview is explained. A brief explanation of the four-stroke engine is given together with the two sensors that are being used: the pressure sensor and the angular sensor. Also, an explanation of how the signals are processed is given. This is crucial in order to filter out noise, but also convert from angle to time domain. 2.1 Combustion cycle in a four-stroke engine A four-stroke diesel engine cycle consists of four different combustion phases, as can be seen in Figure 2.1 below. At the beginning, the inlet valve opens, which allows air to flow into the cylinder. The next step compresses the air, which makes the pressure rise rapidly. At this point, the fuel is being injected into the cylinder. Because of the increase in pressure, temperature rises rapidly, which self-ignites the gas mixture (a diesel engine is a compression ignited engine). The total pressure from the ignition in the gas mixture and the compression push the piston downwards. This generates the output power from the cylinder and happens right after the lever passes the Top Dead Center (TDC) position. When the lever has passed the Bottom Dead Center (BDC) position, the outlet valve opens and the exhaust gas is pushed out of the cylinder. Lastly, when the lever reaches TDC again, the inlet valve opens and the process is repeated. Figure 2.1: The four different combustion phases in a four-stroke engine [1]. Since every four-stroke cycle takes two engine revolutions, the total crankangle swept during one cycle is 4π rad = 720◦. 4 2. System overview 2.2 Pressure sensor for in-cylinder measurements An ICPS will in this thesis be used to measure the pressure signal in one of the 6 cylinders in a 13 liter diesel engine. To measure the pressure during combustion, a piezo-electric pressure sensor (Figure 2.2) is used. The sensor generates a charge from a quartz crystal when exposed to pressure. This sensor can only be used for measuring dynamic pressures (when the pressure changes), the reason being that the charge quickly leaks to zero. We are therefore counting this sensor to drift during longer samples. A way to compensate for this drift is to set the cylinder pressure at BDC before the compression phase equal to the intake manifold pressure [7]. Figure 2.2: A dynamic pressure sensor capable of doing high precision measure- ments during a combustion cycle [2]. The output signal from the sensor is sampled, and is hence discrete. To construct an estimated continuous pressure function p(θ) (where θ = θ(t) is the crank angle), interpolation is used. In Figure 2.3 below, one can see the pressure during the combustion phase of one cylinder. 600 650 700 750 800 850 0 2 4 6 8 10 12 14 16 18 10 6 Figure 2.3: Interpolated pressure signal in the power phase of the combustion. The dotted line marks where the TDC-angle is. 5 2. System overview 2.3 Angular sensor at the flywheel On the flywheel there are 18 teeth and 2 holes (locations where the teeth are missing) every 120◦. They are separated equidistantly with an angle of 6◦. An analog voltage sensor is measuring the periphery of the flywheel and gives a discrete time signal sampled every 6◦. A 32 values toothtime vector is fed into the ECU every 120 Crank- angle Degree (CAD). The first values are the latest ones measured by the sensor (as opposed to the pressure sensor). The missing-teeth time values are computed thanks to an interpolation. 2.4 Cylinder ID Another value we have access to is the cylinder ID. It represents the cylinder in which the combustion TDC is occurring. The whole setup is represented along the flywheel in the Figure 2.4 below, where n is the number of the teeth, p the pressure vector coming from the sensor and ttimes the tooth time vector. The direction of rotation is anticlockwise. TDC𝑐𝑦𝑙1&6 TDC𝑐𝑦𝑙5&2 TDC𝑐𝑦𝑙3&4 𝑝 = 𝑝20, … , 𝑝1 𝑡𝑡𝑖𝑚𝑒𝑠 = [𝑡1, … , 𝑡32] Outputs from the sensors Figure 2.4: Representation of the experimental setup along the flywheel. 6 2. System overview 2.5 Signal processing The angular sensor in Section 2.3 samples in angle domain, with a sample interval of ∆θ = 6◦. Some signals, such as acceleration, that is a time derivative of θ(t), needs to be expressed in time domain in order to approximate their derivatives [3]. This is especially crucial in Section 4.1 in order to find the oscillation functions (oscillations in θ(t)) and acceleration coefficients. It is therefore necessary to interpolate the angular signal into a time signal. This is done with linear interpolation using the function interp1 in Matlab. The tooth times are then sampled with a constant sampling interval h0 = tj+1 − tj. All raw signals that are to be analyzed in Matlab (for the most part identification), are then interpolated with the same method. The number of pulses for the flywheel sensor with 60 teethes is 60 v̄ 60 = v̄ per second, where v̄ is the engine speed in RPM. Hence, the sampling interval becomes h0 = 1 v̄ = 1 60ω̄/(2π) = 2π 60ω̄ , (2.1) where ω̄ is the average engine angular velocity, or more commonly, engine frequency. In [3], it is concluded that no frequencies higher than 12ω̄ (12th engine order) are analyzed. The engine order is the number of the engine natural frequency’s har- monics. Therefore, when doing identification of Fourier series in Section 4.1 the maximum order is N = 24, or 24ω0 half engine orders, where ω0 = ω̄ 2 is the funda- mental frequency of oscillations. Above this frequency a lot of noise occurs and the signal is therefore also low-pass filtered before further analysis. In Figure 2.5 the flywheel velocity is filtered with a bandwidth of ωL = 36ω0. 800 900 1000 1100 1200 1300 1400 950 960 970 980 990 1000 1010 1020 1030 1040 1050 Raw signal Filtered signal Figure 2.5: Low-pass filtered flywheel velocity (solid line), with ωL = 36ω0, and raw signal (dotted line) versus crank angle. The window is for one full combustion cycle [720◦,1440◦]. 7 3 Physical models In this chapter the three models are presented. The models allow to describe the relation between different physical quantities. In this thesis, the focus is on the combustion torque being generated from the cylinders, and its contribution to the flywheel acceleration. This is a key in order to make an estimator that uses angular data combined with pressure sensor data for estimation. 3.1 Torque model with crankshaft The 6 cylinders in the diesel engine are all linked to the crankshaft via a piston and connecting rods. In each combustion cycle, the force from the cylinder pressure performs a work on the piston. The pistons push the lever down on the crankshaft which results in a torque. The main components that are involved are, the combus- tion torque Tcomb, which is the torque created from the pressure, the mass torque Tm that is generated from the inertia, and finally the friction torque Tf , which is all lumped losses from heat between the moving parts. Together, these three torque contributions creates a resulting torque on the crankshaft [7], T = Tcomb − Tm − Tf . (3.1) At the phase where positive work is done on the crankshaft (Tcomb > 0) the mass torque counteracts the motion in Equation (3.1). The net work is positive for the combustion torque, but zero for the mass torque. Friction is however breaking the system at all times, and always generates negative work. 3.1.1 Combustion torque Tcomb In order to describe the relationship between the force from the piston translates to crankshaft torque, a geometrical relationship will be derived [8]. The physical system for one connected piston is shown in Figure 3.1 below. As one can see the angle β is seen from the cylinder, and θ is the crank angle. The model presented below in Figure 3.1 will be used to characterize the displacement of one piston with regards to the combustion top dead center position (TDC). 8 3. Physical models 𝑠 𝜃 𝛽 𝑙 𝑟 𝒄𝒐𝒈 Figure 3.1: The crankshaft piston model: s is the piston’s displacement from the TDC, β is the connecting rod angle, l the connecting rod length, θ the crank angle and r is the crankshaft web length. In this position, the torque given to the crankshaft is equal to zero. The piston stroke s can be expressed as a function of the two angles θ and β as in (3.2). We can note that θi = θ1− 4π 6 ([1, 5, 3, 6, 2, 4]−1) depends on the firing order. The index i stands for the cylinder number i = {1, 2, ..., 6}. In the following sections, we take the rotation of one cylinder as a reference (cylinder one for instance: θref = θ = θ1), s(θi, βi) = (r + l)− (lcos(βi) + rcos(θi)) = l(1− cos(βi)) + r(1− cos(θi)), (3.2) an expression that links θ and β is, r sin(θi) = l sin(βi) . (3.3) Therefore, sin(βi)2 = r2 l2 sin2(θi) cos(βi)2 = 1− r2 l2 sin2(θi) cos(βi) = √ 1− r2 l2 sin2(θi) , (3.4) 9 3. Physical models which inserted into (3.2) gives, s(θi) = r 1− cos(θi) + l r 1− √ 1− r2 l2 sin2(θi)  . (3.5) The power given by the combustion can be expressed as, Pcomb,i(θi) = Fcomb,i(θi) ds(θi) dt = dθi dt Tcomb,i(θi) , (3.6) where Fcomb stands for the force given by the combustion to the piston mass. Thanks to this expression, Tcomb,i(θi) = Fcomb,i(θi) ds(θi) dt dt dθ = Fcomb,i(θi) ds(θi) dθ = (pi(θi)− p0)Apiston ds(θi) dθ︸ ︷︷ ︸ L(θi) = (pi(θi)− p0)L(θi), (3.7) where pi is the pressure in the cylinder chamber i, p0 the atmospheric pressure (we approximate the crankcase pressure by the atmospheric pressure [8]), Apiston is the area of the piston surface in the chamber and L(θi) is called the crank lever function and is given by differentiation of (3.5) L(θi) = Apistonr sin(θi) ( 1 + r l cos(θi)√ 1− r2 l2 sin2(θi) ) . (3.8) Finally, the total combustion torque is given by Tcomb,tot(θ) = 6∑ i=1 (pi(θi)− p0)L(θi). (3.9) 3.1.2 Mass torque Tm All connecting rods and cylinders in the piston assembly have a mass, which gives the system inertia, the purpose of the mass torque model is to collect all inertia contributions into one torque contribution. The mass torque can be calculated as in [8], Tm = (mrotr 2)θ̈ + ( mosc ncyl∑ i=1 v2(θi) ) θ̈ + ( mosc ncyl∑ i=1 v(θi) dv(θi) dθ ) θ̇2, (3.10) where ncyl = 6, v(θi) = ds(θi) dθ , mrot is the total mass from the rotating part (both piston and crankshaft), mosc is the mass of the oscillating part (piston), defined as mosc = mrod lcog l +mpiston , (3.11) mrot = 6mrod ( 1− lcog l ) +mcrank , (3.12) 10 3. Physical models where mrod is the mass of the connecting rod, mpiston the mass of the piston and mcrank is the mass of the crankshaft, approximated as a point mass with mcrank ≈ Jcrank/r 2. For the full derivation of mass torque please see Appendix A. On a complete combustion cycle [0,720◦) the mean value of the mass torque is zero. Hence, when only considering IMEP and mean torque calculations at steady state operating conditions, this contribution can be neglected. However, when identifying torque or pressure functions, this needs to be accounted for. 3.1.3 Friction torque Tf The friction torque can be modeled with the help of Coulomb’s law as [8], Tf,tot = cf 6∑ i=1 ( ds(θi) dθ ) · θ̇, (3.13) where cf is the friction coefficient. A more detailed friction model is given in [3]. Friction is however being neglected when doing model identification with Fourier series, as it is in [3] for instance. When it comes to the linear IMEP model in Section 3.3, the friction contribution is included in the load term. This is however the average friction contribution which is constant for each cylinder for a given operating point. 3.1.4 Driveline torque TD The torque acting after the flywheel is affecting the whole powertrain, since the total torque balance of the system is, Tcomb,tot − Tm,tot − Tf,tot︸ ︷︷ ︸ Ttot −TD − Jθ̈F = 0 (3.14) or simply, Ttot − TD − Jθ̈F = 0 (3.15) where J is the lumped inertia from the flywheel together with the connections in the crankshaft and θ̈F is the angular acceleration of the flywheel. What the flywheel experiences in terms of torque is Jθ̈F = Ttot − TD, (3.16) hence the driveline torque TD is acting as a load. In a test-cell environment this can be a break with a spring-system, whereas in a real powertrain, it is the consequence of the inertia and torque acting from the input shaft (see Figure 3.2). This load torque can be modelled as in [3], TD = JI θ̈I + TI , (3.17) where JI is the inertia for the input-shaft and clutch parts, TI is the input shaft torque and θ̈I the input shaft angular acceleration. 11 3. Physical models 3.2 Powertrain model In [3] the following powertrain model is presented θ̈F = GD(d)TD + 6∑ i=1 Gi(d)Ti , (3.18) where Gi and GD are transfer functions in the derivative domain d, TD the driveline torque and Ti is the individual cylinder torque in Equation (3.1), but with Tf = 0. Each Gi is the physical channel between the cylinder torques and their contributions to the flywheel and GD is the physical channel from the driveline to the flywheel. The model is assumed to be linear in the dynamics between each section. However, in reality there are some non-linearities, but most of the dynamical components are still captured in the transfer functions. This is a combination of a physical model which captures the cylinder torques Ti, but with a black-box approach, capturing dynamics in the transfer functions GD, which makes it a grey-box model. In Figure 3.2 the powertrain is shown with the physical channels. Figure 3.2: Powertrain with cylinder torques, crankshaft, flywheel and the driveline torque TD [3]. 3.3 Linear IMEP model Several values such as CA50, PCP and IMEP can be used to characterize the com- bustion process of any reciprocating engine. In this section, IMEP is explained. IMEP is defined by the following equation [8]: IMEP = Wcycle Vd = 1 Vd ∮ cycle (p(V )− p0)dV , (3.19) which represents the work done, Wcycle, or the integral of the pressure cycle in the PV diagram, during one combustion cycle divided by the volume swept, Vd, by the piston. Figure 3.3 shows a whole pressure combustion cycle against the volume of the chamber. One of the advantages of this value is that it does not depend on the volume of the combustion chamber and therefore the size of the engine [9]. Also, the IMEP represents an average over a whole cycle and is suspected to be less computationally heavy than working on a curve reconstruction. Another advantage of IMEP is that it can very easily be calculated from the mean combustion torque 12 3. Physical models over a cycle T̄comb,cycle. This can be shown by considering the work done during a cycle, Wcycle = Pnr ω̄/2π , (3.20) where P stands for the power during a cycle and is linked to the work by Equation (3.19), nr is the number of revolutions per cylinder (for a four-stroke engine nr = 2) during a cycle and ω̄ is the average engine angular velocity. The following relation also connects torque to power: P = ω̄T̄comb,cycle, (3.21) where T̄comb,cycle is computed as follows [2]: T̄comb,cycle = 1 |θ| ∫ θ Tcomb,tot(θ)dθ = 1 |θ| ∫ θ ncyl∑ i=1 (pi(θi)− p0)L(θi)dθi = ncyl∑ i=1 1 |θ| ∫ θ (pi(θi)− p0)L(θi)dθi = ncyl∑ i=1  1 |θ| ∫ θ pi(θi)L(θi)dθi − 1 |θ| p0 ∫ θ L(θi)dθi︸ ︷︷ ︸ =0  = ncyl∑ i=1 1 |θ| ∫ θ pi(θi)L(θi)dθi . (3.22) Therefore, combining Equation (3.19), (3.20) and (3.21) yields IMEP = 2πnr Vd T̄comb,cycle . (3.23) Therefore, IMEP can be computed as, IMEPi = 2πnr Vd 1 |θ| ∫ θ pi(θi)L(θi)dθi . (3.24) where pi stands for the pressure in cylinder i. 13 3. Physical models Figure 3.3: PV diagram of a complete combustion cycle for an engine running at 1000 rpm and a nominal torque of 1800 Nm. As IMEP is proportional to the mean combustion torque, as shown in Equation (3.23), a simplified, empirical linear model of the overall crankshaft can be expressed with IMEP [6], Jθ̈F = Tcomb,tot − Ttot,cont (3.25) IMEP ∝ Jθ̈F + Ttot,cont, (3.26) where Ttot,cont is the total torque counteracting the combustion torque so Tm,tot − Tf,tot − TD. In other words, IMEP is an affine function of the acceleration of the flywheel under the assumption that the load torque variations are negligible. 3.4 Torque model utilizing a Kalman filter It is possible to model the flywheel acceleration change, coming from combustion torque, by using the torque balance model in [8], i.e. J(θ)θ̈ = Tcomb(θ)− f(θ) · θ̇2 − T ∗load(θ) , (3.27) where Tcomb is the total combustion torque (sum of all 6 cylinders) and T ∗load(θ) = Tf (θ)+Tload(θ) is the extended load torque, caused by both friction and load coming 14 3. Physical models from the driveline or the brake. J(θ) and f(θ) are defined as, J(θ) = mrotr 2 +mosc 6∑ j=1 (dsj(θ) dθ )2 , (3.28) f(θ) = 6∑ j=1 ( mosc dsj(θ) dθ · d 2sj(θ) dθ2 ) , (3.29) resulting from a two mass model approach for the crankshaft. The reason why the inertia is varying is because of the piston’s motion when moving up and down. This puts the center of gravity closer to the crankshaft, and hence the inertia is changing, with the other oscillating parts. The second derivative θ̈ may be formulated as θ̈ = dθ̇ dθ · dθ dt = dθ̇ dθ · θ̇ . Substituting this into Equation (3.27) yields θ̇dθ̇ = 1 J(θ) ( Tcomb(θ)− f(θ) · θ̇2 − T ∗load(θ) ) dθ . (3.30) By integrating both sides we get an equation that is only dependent on the square of the crankshaft velocity, θ̇2(n+ 1)− θ̇2(n) = 2 J(θ) ∫ θ(n+1) θ(n) ( Tcomb(θ)− f(θ) · θ̇2 − T ∗load(θ) ) dθ . (3.31) During the discrete angular step ∆θ = θ(n+ 1)− θ(n), the integral may be approx- imated as θ̇2(n+ 1)− θ̇2(n) ≈ 2∆θ J(θ) ( Tcomb(n)− f(n)θ̇2(n)− T ∗load(n) ) . (3.32) By regarding the square of the crankshaft rotational velocity θ̇2 as a state, we can linearize the model. The linear discrete state space model becomes x1(n+ 1) = ( 1− 2∆θ J(n) · f(n) ) x1(n) + 2∆θ J(n)x2(n) (3.33) with the state variables, x1(n) = θ̇2 (3.34) x2(n) = Tcomb(n)− T ∗load(n) (3.35) x3(n) = x2(n+ 1) (3.36) The Kalman filter employ a Markovian system model i.e. a first order discrete linear system excited by white noise. This behavior is not fully representing the oscillations of combustion and load torque since it is decreasing at higher engine orders [8]. To model this behavior, the second state is put to be the output of a second order low pass system H(z) which is excited by white noise U(z), X2(z) = H(z)U(z) . (3.37) 15 3. Physical models This output is called colored noise, which decreases in power with increasing engine order. H(z) has a double pole on the real axis, and is chosen as H(z) = (1− exp(−δ∆θ))2 (z − exp(−δ∆θ))2 . (3.38) where δ is a time step. In the discrete time angle domain, the state x2(n) can be modelled as, x2(n+ 2)− 2x2(n+ 1) · exp(−δ∆θ) + x2(n) · exp(−2δ∆θ) = (1− exp(−δ∆θ))2u(n) . (3.39) The discrete state space representation now becomes x(n+ 1) = A(n)x(n) +B(n)u(n) y(n) = C(n)x(n) + w(n) where u(n) and w(n) is white noise, and the matrices and states are x(n) = [ θ̇2(n) Tcomb(n)− T ∗load(n) Tcomb(n+ 1)− T ∗load(n+ 1) ]T (3.40) y(n) = θ̇2(n) (3.41) A(n) = 1− 2f(n)∆θ J(n) 2∆θ J(n) 0 0 0 1 0 −exp(−2δ∆θ) 2 · exp(−δ∆θ)  (3.42) B(n) = 0 0 0 0 0 0 0 0 (1− exp(−δ∆θ))2  (3.43) C(n) = [ 1 0 0 ] (3.44) All dicrete steps are expressed as the angle step ∆θ = 2π 60 , since the crankshaft position sensor detects 60 teeths. The Kalman filter state updates can be expressed as [10], x̂n+1 = Anx̂n + Ln(yn − Cnx̂n) (3.45) Ln = AnPnC T n (CnPnCT n +Rn)−1 (3.46) Pn+1 = AnPnA T n +BnQnB T n − AnPnCT n (CnPnCT n +Rn)−1CnPnA T n (3.47) with initial conditions P0 = cov(x0) and x0 = K0, typically we set x0 to the zero vector and the covariance to 10I, where I is a 3x3 unit matrix. Qn and Rn are positive definite matrices that are chosen (tuned), δ is also a parameter that needs to be chosen (tuned). The Kalman state observer can now estimate the difference between combustion torque and load x2(n) = Tcomb(n)− T ∗load(n) by measuring the square of the crankshaft velocity y(n) = θ̇2(n). Appropriate tuning for model noise (errors coming from imperfect model) Qn = σ2 u and measurement noise (coming from the crankshaft position sensor) Rn = σ2 w, will ensure that x̂n, Ln and Pn converges to the reference model, regardless of combustion model [8]. The power of the Kalman filter is based on this proposition, since the model doesn’t have to be perfect, but still very valuable results can be obtained. The simulation results on data for this model can be found in Section 5.3. 16 4 Model identification In this chapter, the different identification processes are explained: the one for identifying Fourier series as well as the one for identifying the models derived in Chapter 3. The data set used is a very large set of data available at Volvo. This data comes from the engine of interest and the same type that an online estimation algorithm would operate on, as in Section 6. The data is sampled with a higher frequency than what is available for this thesis work. A 153 teeth crankshaft sensor is used and hence the sample step is 360/153 ≈ 2.35 degrees, instead of 6 degrees. To simulate a 60 teeth crankshaft sensor from this data, linear interpolation was used to perform down-sampling. In Section 4.3 the data is interpolated every 6 degrees, but in Section 4.1 the maximum sample rate is used (2.35 degrees) when identifying Fourier series. A typical table of the data can be seen in Appendix A.2 and the data is processed according to Section 2.5 before doing analysis in Matlab. 4.1 Identification of Fourier series When the engine is at steady state operating conditions, all raw signals can be approximated by Fourier series [11]. This is possible since each raw signal is periodic in time. Since this is the case, the engine model presented in (3.2) is also valid when transformed to Fourier series. Any signal x(t) (pressure, torque and angle) can be expressed by Fourier series, i.e. x(t) = x{N,θ0} ≡ A0,x + N∑ n=1 An,xcos(nθ0(t)/2) + Bn,xsin(nθ0(t)/2) , (4.1) where θ0 ≡ K + ω̄t ≡ K + 2ω̄0t is the angular function, ω̄ is the average angular velocity (frequency) of the flywheel, also called "engine order" and K is a constant. Ax,n and Bx,n are called Fourier coefficients, and need to be identified for each signal. The order of the series is in theory N → ∞ and will, if the series converges, give an exact solution to x(t). For practical use N is finite, and the number of terms in the series becomes 2N + 1. Since each combustion cycle consists of two crank angle cycles, ω0 ≡ ω̄/2 is the normal engine frequency. To obtain the normal engine frequency in the engine, each crank angle θ(t) is plotted versus the time stamps tj, j = 0, 1, 2...ns, and fitted with least squares. This can be seen in Figure 4.1. 17 4. Model identification 0.2 0.202 0.204 0.206 0.208 0.21 0.212 0.214 0.216 0.218 0.22 790 800 810 820 830 840 850 860 870 880 890 Figure 4.1: Fitted least squares θ0 = 2ω0t+K with tooth time data θ(tj) for each crank angle sample. The engine frequency is the slope in the fit. Once the engine frequency is determined, it is possible to construct the Fourier series using this data. For each signal x(t) we calculate the Fourier coefficients Ax,n and Bx,n by solving the least squares problem below (it is also possible to calculate these with traditional Fourier integrals).  x(tj=1) x(t2) x(t3) ... x(tns)  ︸ ︷︷ ︸ b =  1 cos(θ0,11/2) sin(θ0,11/2) . . . cos(θ0,1n/2) sin(θ0,1n/2) 1 cos(θ0,21/2) sin(θ0,21/2) . . . cos(θ0,2n/2) sin(θ0,2n/2) 1 cos(θ0,31/2) sin(θ0,31/2) . . . cos(θ0,3n/2) sin(θ0,3n/2) ... ... ... ... ... ... 1 cos(θ0,ns1/2) sin(θ0,ns1/2) . . . cos(θ0,nsn/2) sin(θ0,nsn/2)  ︸ ︷︷ ︸ A  A0 A1 B1 ... An Bn  ︸ ︷︷ ︸ x where n = 1, 2, 3, ..., N is the order of the Fourier series, and t = 1, 2, 3, ..., ns is the number of time samples. The system of equations gives an over-determined system Ax = b which can be solved by the regular solution to the least squares problem, ATAx = ATb (4.2) where x is the sought solution of the Fourier coefficients. To identify the engine dynamics, one needs to consider θ0 with zero mean, i.e. deviations from its average value. We define this function as the oscillation function ψ(t) ≡ θ(t)− θ0(t) = θ(t)− (2ω0t−K). (4.3) 18 4. Model identification Since ψ(t) has zero mean we can express it as ψ(t){N,θ0} = A0︸︷︷︸ =0 + N∑ n=1 An,ψcos(nθ0/2) + Bn,ψsin(nθ0/2) = N∑ n=1 Cn,ψe jnθ0/2 (4.4) where Cn,ψ is the complex Fourier coefficients of ψ(t) given by Cn,ψ = (An,ψ − jBn,ψ)/2, ∀ n = {1, 2, 3, ..., N}, (4.5) only expressed in the dynamical coefficients n = 1, 2, 3, ..., N . An important property of ψ(t) is that d2 dt2 ψ(t) = θ̈ − d2 dt2 (2ω0t−K)︸ ︷︷ ︸ =0 = θ̈ which means the θ̈ also is zero mean. This lets us relate the Fourier coefficients in ψ(t) directly to acceleration θ̈ on the flywheel. To see this, we express both signals in complex Fourier series and compare the coefficients: θ̈ = d2 dt2 N∑ n=1 Cn,ψe jnθ0/2 = N∑ n=1 −(nω0)2Cn,ψ︸ ︷︷ ︸ Cn,θ̈ ejnθ0/2. This means that it is possible to retrieve the coefficients of the accelerations signal by only considering the flywheel angle, tooth time and subtraction of the average. The relation then becomes Cn,θ̈ = −(nω0)2Cn,ψ, ∀ n = {1, 2, 3, ..., N}. (4.6) Figure 4.2 below shows the fitted Fourier series to ψ(t) with N = 6. 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 -0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02 Figure 4.2: Fitted Fourier series of ψ(t) in a steady state engine condition. 19 4. Model identification The interpolated flywheel acceleration signal θ̈ is fitted with Fourier series of degree N = 6 and N = 24 in Figure 4.3 and 4.4 below. Cn,θ̈ was determined by solving (4.2) and using (4.5). The higher the order of the Fourier series, the better the fit. This is because the signal consists of higher frequency components that are periodic. How well a given Fourier series identifies the raw signal is measured in the root mean square level as an error vector ej = x(tj) − x̂(tj), where x̂(tj) is the fitted Fourier signal for time sample j, erms = √√√√ 1 ns ns∑ j=1 e2 j . (4.7) All the calculated Fourier coefficients Cn,T , Cn,θ̈ and Cn,D are identified and saved for a given N . These are used in Section 4.2 to identify the grey box model between the flywheel and torque contribution. 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 -1000 -500 0 500 1000 1500 Figure 4.3: Fitted Fourier series with N = 6 to the acceleration signal θ̈(t). The root mean square error is erms = 262.52 rad/s2. 20 4. Model identification 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 -1000 -500 0 500 1000 1500 Figure 4.4: Fitted Fourier series with N = 24 to the acceleration signal θ̈(t). The root mean square error is erms = 97.35 rad/s2. In the same way it is possible to determine the torque signal Ti for each individual cylinder i, and extract the set of coefficients Cn,Ti , i = {1, 2, .., 6}. Since friction is not identified, the torque signal is T = Tcomb + Tm, for each piston assembly. The fit for N = 6 and N = 24 can be seen in Figure 4.5 and 4.6 below. 800 900 1000 1100 1200 1300 1400 -2000 -1000 0 1000 2000 3000 4000 5000 Figure 4.5: Fitted Fourier series with N = 6 of T = Tcomb+Tm for the first cylinder (cylinder 1). The crank angle range is [720◦,1440◦], one full combustion cycle. Root mean square error erms = 530.68 Nm. 21 4. Model identification 800 900 1000 1100 1200 1300 1400 -2000 -1000 0 1000 2000 3000 4000 5000 6000 Figure 4.6: Fitted Fourier series with N = 24 of T = Tcomb + Tm for the first cylinder (cylinder 1). The crank angle range is [720◦,1440◦], one full combustion cycle. Root mean square error erms = 33.36 Nm. In Figure 4.7 and 4.8 below, the fitted Fourier series is plotted versus the driveline torque model TD, from Equation (3.17). The coefficients Cn,D are identified in the same way as before. 800 900 1000 1100 1200 1300 1400 -2500 -2400 -2300 -2200 -2100 -2000 -1900 -1800 -1700 Figure 4.7: The driveline torque TD with its corresponding fitted Fourier series (dotted line), for N = 6, with an error erms ≈ 64 Nm. 22 4. Model identification 800 900 1000 1100 1200 1300 1400 -2500 -2400 -2300 -2200 -2100 -2000 -1900 -1800 -1700 Figure 4.8: The driveline torque TD with its corresponding fitted Fourier series (dotted line), for N = 24, with an error erms ≈ 41 Nm. 4.2 Powertrain model identification To identify the model in Section 3.2 we consider the general form of a transfer function G with input u and output y [12], N∑ n=1 Cy,ne jnω0t = N∑ n=1 G(jnω0)Cu,nejnω0t. (4.8) The same equality can be expressed as Cy,n = G(jnω0)Cu,n ∀n = 1, 2, 3, ..., N (4.9) Since the real coefficients Ay,n and By,n can be expressed as Ay,n = 2Re(G(jnω0)Cu,n) (4.10) By,n = −2Im(G(jnω0)Cu,n) (4.11) we can find the real and imaginary part of G by using G(jnω0)Cu,n = Re(G(jnω0))Re(Cu,n)− Im(G(jnω0))Im(Cu,n)+ (4.12) j(Im(G(jnω0))Re(Cu,n) + Re(G(jnω0))Im(Cu,n)) (4.13) which can be expressed as G(jnω0)Cu,n = 1/2 (Re(G(jnω0))Au,n + Im(G(jnω0)Bu,n)︸ ︷︷ ︸ Ay,n − (4.14) 1/2j (−Im(G(jnω0))Au,n + Re(G(jnω0))Bu,n)︸ ︷︷ ︸ By,n . (4.15) 23 4. Model identification This yields the following system of equations[ Ay,n By,n ] = [ Au,n Bu,n Bu,n −Au,n ] [ Re(G(jnω0)) Im(G(jnω0)) ] . (4.16) For each torque contribution in Equation (3.18) the coefficients are derived. Super- position of all contributions then gives the total system of equations (when applied to the powertrain model) [ Aθ̈,n Bθ̈,n ] = [ ATD,n BTD,n AT1,n . . . BT6,n BTD,n −ATD,n BT1,n . . . −AT6,n ]  GD,Re,n GD,Im,n G1,Re,n ... G6,Im,n  (4.17) where the parameters GRe,n = Re(G(jnω0)) and GIm,n = Im(G(jnω0)) is a represen- tation with a shorter notation. This system is under-determined, and needs to have at least 7 operating points for a constant engine frequency ω0. However, in [12] it is concluded that this is not a sufficient condition. In order to identify the transfer function, a parameterization on the following form is suggested instead GD(jω) = χ(ω)φD,Re + jχ(ω)φD,Im, (4.18) Gi(jω) = χ(ω)φi,Re + jχ(ω)φi,Im, i = 1, . . . , 6 (4.19) where φD,Re, φD,Im, φi,Re and φi,Im are parameter vectors of length nP + 1, which are the coefficients in the polynomial χ(ω) = [ 1 ω ω2 . . . ωnP ] . (4.20) Expressing Equation (4.17) with the help of this parameterization now provides the vector GD,Re,n GD,Im,n G1,Re,n ... G6,Im,n  =  χ(nω0)φD,Re χ(nω0)φD,Im χ(nω0)φ1,Re ... χ(nω0)φ6,Im  =  χ(nω0) 0 0 . . . 0 0 χ(nω0) 0 . . . 0 0 0 χ(nω0) . . . 0 ... ... ... ... ... 0 0 0 0 χ(nω0)   φD,Re φD,Im φ1,Re ... φ6,Im  (4.21) Note that the zeros in this matrix represents the zero vector. For nP = 0 the matrix χ(nω0)I becomes the identity matrix I and creates 14 unknowns with 2 equations, therefore it can be solved uniquely with n = 7 at a single operating point. With nP = 1 there are 28 unknowns and can be solved with n = 7 for two operating points, or n = 14 for one operating point. As an example we set up the equation systems for nP = 1 using (4.21) [ Aθ̈,n Bθ̈,n ] = [ ATD,n BTD,n AT1,n . . . BT6,n BTD,n −ATD,n BT1,n . . . −AT6,n ] 1 nω0 0 0 0 0 . . . 0 0 0 0 1 nω0 0 0 . . . 0 0 0 0 0 0 1 nω0 . . . 0 0 ... ... ... ... ... ... ... ... ... 0 0 0 0 0 0 0 1 nω0   φD,Re1 φD,Re2 φD,Im1 φD,Im2 φ1,Re1 φ1,Re2 ... φ6,Im1 φ6,Im2  (4.22) 24 4. Model identification where e.g. Im1 and Im2 denotes the vector elements and n = 1, 2, 3, ..., N . The system is solved using least squares, just as (4.2). Then for each n the imaginary and real part of G(jnω0) is extracted. This yields a transfer function G that is a matrix of 7 rows and N columns, when converted to complex Fourier coefficients. It is now possible to describe the powertrain model (3.18) with the use of Fourier coefficients, Cn,θ̈F = GD(jnω0)Cn,D + 6∑ i=1 Gi(jnω0)Cn,i, ∀n = 1, 2, 3, . . . , N. (4.23) In Figure 4.9 below, the calculated Cn,θ̈F converted to time domain, compared to the filtered acceleration signal in the engine is displayed. This confirms that the model has been identified correctly in a certain operating region. The conversion from coefficients to time domain, is done by converting Cn,θ̈F to real coefficients and then multiply A in (4.2) by this vector. The model was identified on operating point 7 and 8 (see Appendix A). Figure 4.9: The fitted acceleration output from the powertrain model and the measured acceleration, for order N = 14. 25 4. Model identification 4.3 Linear IMEP model identification To identify the model in Equation (3.25), the first step is to extend it for all the individual cylinders. Indeed, each cylinder stroke brings a specific contribution to the crankshaft rotational speed that can be seen in Figure 4.10. Therefore, it is possible to compute the flywheel acceleration provided by each cylinder and link it to the corresponding average combustion torque. This is necessary in order to identify each individual cylinder contribution. Figure 4.10: Flywheel velocity and torque produced before the flywheel on the same window. Cylinder 1 is indicated with a cross at its maximum torque, firing order follows. One can note that the minimum angular speed occurs just after the TDC (6 degrees after TDC to be precise) due to the contributions of the other torques at stake during the process. To identify the linear IMEP model, a small window of several degrees (defined by the positions x and y after TDC) where the flywheel angular velocity is almost linear is chosen in order to compute the flywheel acceleration. On the other hand, average torque (or IMEP) is computed for a specific cylinder. The overall equation system for identification is, 26 4. Model identification T̄comb,1 = CconvIMEP1 = J1θ̈[θx,θy ],1 + Ttot,cont,1 T̄comb,2 = CconvIMEP2 = J2θ̈[θx,θy ],2 + Ttot,cont,2 ... T̄comb,6 = CconvIMEP6 = J6θ̈[θx,θy ],6 + Ttot,cont,6 (4.24) where T̄comb,i is the average combustion torque provided by the cylinder i = {1, 2, 3, . . . , 6} during an entire cycle (720 degrees), Cconv = Vd 4π is the proportionality coefficient between average torque and IMEP, IMEPi is the IMEP computed for cylinder i, Ji and Ttot,cont,i the coefficients to determine, assumed to be constant, re- spectively representing the total inertia and torques - other than combustion torque - linked to cylinder i. θ̈[θx,θy ],i is the acceleration during a windows where cylinder i is contributing on the overall acceleration. The window (noted (y, x)) is defined by two different positions on the velocity curve located after TDC as illustrated in Figure 4.11. For this plot, the window is (8, 1). Figure 4.11: The flywheel velocity during a smaller range, sampled as in the real engine (every 6 degrees). The window here is (8,1). To compute the acceleration, the velocity difference between the two points is com- puted and divided by the amount of time between them (since the curve is almost linear after TDC). The acceleration θ̈(t) is computed as follows θ̈(t) = θ̇y − θ̇x ∆txy , (4.25) where θ̇x (respectively θ̇y) is the velocity at point x (respectively y) and ∆txy the amount of time between the two points. 27 4. Model identification Then, IMEP is plotted against the flywheel acceleration (both averaged after a 15 second run, or half the entire run). An approximately linear relationship is then obtained (Figure 4.12) and a least squares method is used to determine the one degree polynomial for optimal fit. 1.5 2 2.5 3 3.5 1.8 2 2.2 2.4 2.6 2.8 x 10 6 Crankshaft velocity change, [rad/s] IM E P a v er a g ed o v er th e 3 0 se co n d s ru n , [P a ] Cyl1 Cyl2 Cyl3 Cyl4 Cyl5 Cyl6 Figure 4.12: IMEP for four operating points at the same engine speed (1000 RPM) but different nominal torques. The x-axis is expressed in velocity change (velocity difference between 2 points) and not acceleration, time difference being almost constant. The least squares method has been iterated for 11 window configurations and a correlation coefficient has been computed for each. The results are displayed in Figure 4.13. 28 4. Model identification Figure 4.13: Correlation coefficient along the different operating points for each configuration. Matlab correlation function was used for this curve. The black dashed line represents the average over all the different OPs. As shown, the two best correlation coefficients stand for configurations (8,1) and (9,2). The spread of the correlation coefficients is also small for those configurations. With the current experimental setup, the latest tooth-time value available in every cy- cle is TDC+9, therefore, the latest velocity that can be computed is for TDC+8. Hence, for implementation, the only possible configuration is (8,1). The calculated correlation coefficients for this window are given in Table A in Appendix. All are calculated with Pearson’s correlation coefficient, ρX,Y = cov(X, Y ) σXσY = E[(X − µX)(Y − µY )] σXσY , (4.26) where µX or µY is the mean value and σX or σY the standard deviation for each random variable X or Y . The coefficient is a ratio ρX,Y ∈ [−1, 1]. After choosing the window with the best fit, the coefficients α1,2,...,6 and γ1,2,...,6 are used in the following equation system: T̄comb,1 = α1 ¨̄θ[6◦,48◦],1 + γ1 (4.27) T̄comb,2 = α2 ¨̄θ[6◦,48◦],2 + γ2 (4.28) ... (4.29) T̄comb,6 = α6 ¨̄θ[6◦,48◦],6 + γ6, (4.30) (4.31) where the coefficients are determined by a linear least squares fit to the points in Figure 4.12. 29 5 Torque estimation In this chapter both the estimation from Fourier series and the IMEP model is explained, as well as torque estimation from the Kalman filter. When using IMEP the average combustion torque is related by a constant as shown in Section 3.3, however the estimations are displayed in this unit. All estimations in this chapter are simulated in Matlab on operating points (see Appendix A.2), coming from engine data by Volvo. 5.1 Torque estimation using Powertrain model In Section 4.2 the representation of the powertrain model was done in terms of complex Fourier coefficients. Equation (4.23) is the key to make an estimation of the cylinder torques. The goal is to make an estimation of five cylinder torques using one pressure sensor, not having the information of all six. This lets us only know the torque contribution of one torque function, which in this case is T1. The Equation (4.23) now becomes a multiple input and single output system (MISO) and hence it has no unique solution. The new equation becomes Cn,θ̈F = Gn,DCn,D +Gn,1Cn,1 + 6∑ i=2 Gn,iC ∗ n,i, ∀n = {1, 2, . . . , N}, (5.1) where the coefficients C∗n,2 . . . C∗n,6 denotes the unknown coefficients for cylinder torques T2 . . . T6. The idea is now, based on inspiration of a separation algorithm in [13], to separate each cylinder torque contribution from the flywheel acceleration based on an initial estimate. This is however not done in time or angle domain, as in [13], but in frequency domain. It is possible to know the phase-shift in the coef- ficients, since each cylinder stroke is shifted with 120◦ in crank angle, or 2π 3 radians. The firing order (FO) is also known: FO = [ 1 5 3 6 2 4 ] , (5.2) and indicates which cylinder-ID is fired. This means that cylinder 1 is phase-shifted 2π 3 w.r.t cylinder 5, and cylinder 3 is phase-shifted 22π 3 w.r.t cylinder 1, and so on. One cylinder is also at TDC in the combustion phase, when θF = n4π, n ∈ Z, which in this case is the first cylinder (cylinder 1). This lets us relate the crankshaft angle θi for each cylinder i = 1, 2, . . . , 6 , by the equation θi = θF − ζi, (5.3) 30 5. Torque estimation where ζ is given by ζ = 2π 3 [ 0 4 2 5 1 3 ] , (5.4) and i is the index for the elements in ζ. The phase-shift ζi is approximately equal to a constant shift in time τi (if angle is interpolated into time), which is τi ≈ ζi ω0 . (5.5) In frequency domain, such a time delay corresponds to a phase-shift, since for the general case T (t− τ) = N∑ n=1 Cne jnω0(t−τ)/2 = N∑ n=1 e−jnω0τ/2Cne jnω0t/2. (5.6) The phase shift for each cylinder torque coefficient Cn,i is thus Cn,ζi = e−jnω0τi/2 ≈ e−jnζi/2, ∀n = {1, 2, . . . , N}, (5.7) where i = {1, 2, . . . , 6}. Now it is possible to make an initial estimation Ĉn,i, i = {2, 3, . . . , 6}, based on Cn,1, by using this phase-shift, Ĉn,2 = Cn,1e −jnζ2/2 = Cn,1e −jn 2π 3 4/2 = Cn,1e −jn 4π 3 Ĉn,3 = Cn,1e −jnζ3/2 = Cn,1e −jn 2π 3 2/2 = Cn,1e −jn 2π 3 Ĉn,4 = Cn,1e −jnζ4/2 = Cn,1e −jn 2π 3 5/2 = Cn,1e −jn 5π 3 Ĉn,5 = Cn,1e −jnζ5/2 = Cn,1e −jn 2π 3 1/2 = Cn,1e −jnπ3 Ĉn,6 = Cn,1e −jnζ6/2 = Cn,1e −jn 2π 3 3/2 = Cn,1e −jnπ. By inserting these initial estimates into Equation (5.1), it is possible to estimate each torque contribution Ĉ∗n,2, Ĉ∗n,3,. . . ,Ĉ∗n,6 for all n = {1, 2, 3, . . . , N}. In Figure 5.1 it is shown how the separation is done for retrieving the complex Fourier coefficients for the second cylinder. It is only necessary to switch the index in order to solve this for cylinders 3 to 6. 31 5. Torque estimation Figure 5.1: Schematic overview on how the separation strategy is made. This example is estimating the coefficients in the second cylinder, Ĉ∗n,2. The results from the separation strategy can be seen in Figure 5.2 below, where one full combustion cycle, 720◦ is plotted. Each individual torque got reconstructed with its corresponding dynamics and phase-shift. Cylinder 1 is always the reference torque (which is known) and therefore has no estimate. 800 900 1000 1100 1200 1300 1400 -3000 -2000 -1000 0 1000 2000 3000 4000 5000 6000 Figure 5.2: Cylinder separated torques, with no faults, in firing order T1 (green), T5, T3, T6, T2 and T4. The dotted lines are the estimated torques from the separation algorithm. This estimate has an error of erms ≈ 138 Nm. Since the model was identified on two operating points, one with faults and one without faults, it is of interest to see how the separation strategy can catch faults. The results from this can be seen in Figure 5.3 below. 32 5. Torque estimation 800 900 1000 1100 1200 1300 1400 -3000 -2000 -1000 0 1000 2000 3000 4000 5000 6000 Figure 5.3: Cylinder separated torques, with −10% fuel in cylinder 3, in firing order T1 (green), T5, T3, T6, T2 and T4. The dotted lines are the estimated torques from the separation algorithm. As one can see, cylinder 3 has less peak torque that the estimator catches. This estimate has an error of erms ≈ 131 Nm. It is clear that the separation catches cylinder 3 with −10% less fuel injection. What is important to notice is that this comes from the changes in dynamics at the flywheel Cn,θ̈F and driveline torque Cn,D as a consequence of this dynamical change. In general, most of the correction comes from the flywheel, since it is more sensitive than the driveline torque, TD. To see this we make the driveline Cn,D fix to the identified value at operating point 7 and then switch to estimate on the faulty operating point 8. This means that the new faulty driveline signal torque TD is not in use and Equation (5.1) becomes Cn,θ̈F = Gn,DCn,fix +Gn,1Cn,1 + 6∑ i=2 Gn,iC ∗ n,i, ∀n = {1, 2, . . . , N}, (5.8) where Cn,fix are the coefficients of the driveline from operating point 7. The estima- tion results for the faulty operating point 8, with a constant driveline torque Cn,fix can be seen in Figure 5.4 below. The error has gone up to erms ≈ 139 Nm from previously erms ≈ 131 Nm, by only calculating the fault with the flywheel and not the driveline. 33 5. Torque estimation 800 900 1000 1100 1200 1300 1400 -3000 -2000 -1000 0 1000 2000 3000 4000 5000 6000 Figure 5.4: Cylinder separated torques, with −10% fuel in cylinder 3, in firing order T1 (green), T5, T3, T6, T2 and T4. The dotted lines are the estimated torques from the separation algorithm. As one can see, cylinder 3 has less peak torque that the estimator catches. This estimate has an error of erms ≈ 139 Nm. 5.1.1 Sensitivity analysis By changing the input torque T1 with corresponding coefficients Cn,1 by an offset in gain, it is possible to see how sensitive the estimator is to errors in torques. In Figure 5.5 below the offset is between 0.8Cn,1 and 1.2Cn,1 for the coefficients. -20 -15 -10 -5 0 5 10 15 20 100 150 200 250 300 350 400 450 500 Figure 5.5: Sensitivity to offsets in gain. Offsets in incorrect phase in the coefficients are plotted in Figure 5.6 below. As one can see the model is quite sensitive to imperfections in phase. 34 5. Torque estimation -15 -10 -5 0 5 10 15 100 200 300 400 500 600 700 Figure 5.6: Sensitivity to offsets in phase. 5.2 Estimation using Linear IMEP model The estimation of the IMEP in each cylinder is based on the four following steps: 1. As explained in Section 4.3, the calibration of the engine that consists of the calculation of mean IMEP and mean (crankshaft velocity) amplitudes, or acceleration, during a run of several seconds for as many operating points as required at the same angular speed for all six cylinders. Then, a basic fitting between amplitudes and IMEP with a degree one polynomial function using least squares method is made for every cylinder. The 12 coefficients (two coefficients per cylinder) need to be saved in the truck ECU. 2. While the engine is running at a steady state operating point at the desired speed, the ECU computes the amplitude of the crankshaft velocity every 120 degrees for all the cylinders. 3. Thanks to the linear fitting, the ECU can estimate the corresponding IMEP to the amplitudes made during step 2. 4. An Exponential Moving Average (EMA) with parameter α = 0.4 is applied to filter the IMEP values. It is worth noting that the calibration, or identification process, requires to have a pressure sensor in all six cylinders. The EMA is defined as y(t) = αycalc(t) + (1− α)y(t− 1) , (5.9) where y(t) is the averaged value at time t, ycalc(t) is the calculated value at time t and α is the EMA coefficient that determines the contribution of the calculated values to the average. 35 5. Torque estimation 0 20 40 60 80 100 120 1.6 1.8 2 2.2 2.4 2.6 2.8 3 x 10 6 subset number IM E P va lu e, [P a ] IMEP calculated IMEP estimated Figure 5.7: Real and estimated averaged IMEP for cylinder 2 during a simulated run at OP 9 to 12. Table 5.1: Average relative error [%] between real and estimated averaged IMEP EMA-averaged for a 4 seconds simulation time. Average relative error [%] OP Cylinder 1 Cylinder 2 Cylinder 3 Cylinder 4 Cylinder 5 Cylinder 6 9 0.2066 0.0110 0.8162 0.0208 0.3982 0.3510 10 0.7757 0.5679 0.1727 1.1289 0.4996 0.1233 11 0.6471 0.1438 0.0274 0.2745 0.0845 0.2509 12 4.4062 6.9195 2.6772 3.7561 6.8469 2.4450 13 1.9033 1.5011 1.3460 0.7669 2.4705 1.1012 14 1.5719 2.9090 2.1576 0.8912 2.6017 1.1181 15 0.7892 0.5089 0.5231 0.7704 0.1932 0.3310 16 0.2737 0.3862 1.1064 0.0527 0.3879 0.4008 As shown in Table 5.1, one can estimate IMEP with an error below 3% for each cylinder taking an average of the amplitudes over ten values for operating points 9 to 11. For operating point 12, the error is greater and reaches more than 6% for some cylinders. This is because the relation between IMEP and amplitudes is not linear anymore for these OPs when reaching high torques (2650 Nm), and the estimator therefore has a tendency to overestimate IMEP. As stated in Section 3.3, the main goal of this work is to be able to determine if a cylinder is weaker or acting differently from the others during a certain period of time (more than 20 seconds for instance). As a consequence, we want our system to catch a trend and correct it rather than being able to control every combustion cycle. Figure 5.8 is a flowchart showing how the IMEP estimator works on a real engine. Note that the computation of T720 from P720 is achieved thanks to the Riemann sum of Equation (3.22) [2]: 36 5. Torque estimation T720 = T̄comb,cycle = ncyl 1 n n−1∑ k=0 L(θi(k))p(θi(k)) , (5.10) where n is the number of samples during a full cycle of 720 degrees, that is to say 720/6 = 120. Retreive P120, Toothtimes120, CylinderID Compute the velocity difference (from toothtimes) between 2 points of the velocity rise Compute IMEPi from this value and the saved coefficients EVERY 120° CylinderID = 1 Concatenate P120 to have a full cycle and obtain P720 CylinderID = i Pressure sensor location = j EVERY 720° Compute T720 from the P720 vector of the whole cycle IMEP > 0 Compute the EMA of IMEPi END IMEP values corrected with the calculated offset when CylinderID = 1 Convert T720 to IMEPj (this value is the real computed IMEP value) Compute the relative error between the real and estimated IMEPj value Error < 7% END (Warning) Correction is saved and fed on all IMEP estimates YES NO NO YES NO YES Conversion to mean torque: Mean torque estimate for cylinder i SENSOR DATA Figure 5.8: Flowchart showing the estimator principle that will be used for Tar- getLink implementation. 37 5. Torque estimation 5.3 Torque estimation utilizing a Kalman filter The estimation of x2(n) (presented in Section 3.4) is displayed Figure 5.9 below. This simulation is done on data coming from a 13 L 6 cylinder heavy duty engine at Volvo. The table of data used is operating point 9 in Table A.1, and is interpolated to be sampled at ∆θ = 6◦. The Kalman filter was coded in Matlab and tuned with the parameters δ, Rn and Qn. The solid line comes from combustion torque estimation using 6 ICPS’s in each cylinder and subtracting the load torque, which in this case is the driveline torque TD. Hence calculated as, x2,measured = 6∑ i=1 pi(θi(n))L(θi(n))− TD(n) . (5.11) 800 900 1000 1100 1200 1300 1400 -5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000 6000 Figure 5.9: The estimated state x̂2(n) by the Kalman state observer (dotted) and the measured state x2,measured coming from engine data (solid). By neglecting friction we can make the estimate T̄ ∗load(n) ≈ T̄D (5.12) Hence, we neglect the slowly varying oscillations in the load torque and assume it is constant during a combustion cycle. Therefore, it is possible to estimate Tcomb by adding back this value, i.e. T̂comb(n) = x̂2(n) + T̄D . (5.13) We can then use similar methodology as the separation strategy in Section 5.1, but subtracting the other torque contributions, assuming they have the same pressure 38 5. Torque estimation p2(θ) (assuming the ICPS is located in cylinder 2) coming from the ICPS. For an example, estimating the combustion torque in cylinder 1, yields the equation T̂comb,1(n) = x̂2(n) + T̄D︸ ︷︷ ︸ T̂comb − 6∑ i=2 p2(θi(n))L(θi(n)) , (5.14) the result of which is displayed in Figure 5.10 below. 800 900 1000 1100 1200 1300 1400 -3000 -2000 -1000 0 1000 2000 3000 4000 5000 6000 7000 Figure 5.10: The estimated torque by the Kalman state observer (dotted) and the measured torque coming from the pressure sensor (solid), for combustion torque in cylinder 1. 39 6 Implementation on engine In this chapter the implementation part of the thesis work is explained. To validate the IMEP estimator explained in Section 5.2 an online estimation is performed using a test bench set-up. The tests were performed on a Volvo 13 L diesel engine. 6.1 Experimental set-up and implementation There are two layers when running a program in an engine [2]: the platform layer, which is retrieving the raw data from the sensors and converting them into usable data (Vector of in-cylinder pressures, angular position or tooth-time vectors for instance) that is saved during 120 degrees and then, overwritten. The second layer includes estimators and controllers and is named ’application Layer’. This is the layer in which the TargetLink estimator is operating, as shown in Figure 5.8. The engine in which we are performing the validation of our IMEP estimator is a Euro-VI 13 liter diesel engine, mounted in a test cell with a dynamo (brake) to simulate the driveline. The ECU in which the TargetLink model is flashed is an Engine Management System (EMS) v2.3 and the data retrieval is made with Vision v5.2 software. The device used to convert piezo-electric charge to measurable voltage is an AVL amplifier (Microifem). It outputs values from 0 to 10V [2]. Since the EMS needs a voltage range lying between 0 and 5V, a divider is used with a 2.407 coefficient. The ADC-converter needs tuning in order for the gain to be calibrated for our set up. The calibration set-up is shown in Figure 6.1. Since the sensor has a given conversion from pC to Bar, one can use the emulator in order to tune the conversion from analog signal to the correct digital output in Vision. 40 6. Implementation on engine Emulator AVL Voltage- divider EMS ComputerVision Figure 6.1: Calibration set-up for tuning the ADC-gain. The test bench set-up for testing the IMEP estimator is displayed in Figure 6.2 below. All sensors are connected to the EMS that is sampling the pressure from the ICPS and crankshaft position at the same rate. The EMS is then connected to the engine which is running at a specific operating point. Everything is controlled via the test bench platform, which controls the fuel injection, and produce an indicated torque (which we use as a reference in Section 6.2). This indicated torque would simulate what a driver in an actual truck demands as torque when pushing the pedal. ICPS AVL Voltage- divider EMS ComputerVision Figure 6.2: Experimental set-up for doing online estimation with the IMEP esti- mator running on the EMS. 41 6. Implementation on engine Only operating points 13 to 16 were studied in the engine and the chosen EMA coefficient was 0.05. 6.2 Online estimation using the IMEP model Figure 6.3 displays the estimated global torque by the IMEP-based estimator in a run ranging from operating point 13 up to operating point 16. To achieve this result, the estimated Tcomb was computed as follows T̄comb = Vd 4π 6∑ i=1 IMEPi . (6.1) Figure 6.3: Reference combustion torque and estimated one thanks to the IMEP estimator during a run from OP 13 to 16 with an EMA coefficient of 0.05. 42 6. Implementation on engine Figure 6.4: Estimated IMEP for all cylinders from OP 13 to 16. Figure 6.4 shows the estimated IMEP in every cylinder during the run. The spread between cylinders is around 20%. Since the experimental setup does not include a pressure sensor in all of the 6 cylinders (except cylinder 2), it is not possible to plot and compare with the real IMEP values. As shown on Figure 6.3, the estimator follows the combustion torque trend with an offset for lowest and highest torque values. 43 7 Discussion In this chapter the discussion part of the report is presented. The IMEP estimator, estimation using Fourier series and the torque estimation with a Kalman filter are discussed in separate sections. 7.1 IMEP estimator The first thing to underline is that the Tcomb estimator (c.f. Figure 6.3) provides accurate results. The relative error between the reference and the estimation is always less than 3% and less than 1% for OP 14, 15 and 16. In spite of a very low EMA coefficient (0.05), the combustion torque value does not have a significant delay compared to the reference. This is because the cycle is updated at a very high rate, roughly every 18 milliseconds at this engine speed. One disadvantage of this filter is that it could smoothen some high-frequency dynamics; but since the aim is not to develop a cycle to cycle controller, a more responsive system is not needed. Figure 6.4 represents the contribution of each cylinder individually. One impor- tant thing to note is that the spread of values is significant between cylinders. This result is perplexing. The maximum difference in combustion torque between cylin- ders computed during simulation was not exceeding 5%, whereas this plot shows a maximum of 20% difference between the cylinders. This leads us to think that a calibration of the slope coefficients computed in Chapter 4.3 on the specific test- cell engine would be needed. Additionally, the correction added every 720 degrees (right-hand part updated every 720 degrees on the Flowchart 5.8) has not been able to work since the estimator that relies on a former thesis [2] has a tendency to over-estimate the combustion torque. Therefore, the relative error between the real and estimated IMEP value was greater than 7%, which terminated this part of the system, see Figure 5.8. Considering the data signals used to compute the slope coefficients, one can no- tice that in cylinder 1 the combustion torque is slightly higher than the average combustion torque of all the 6 cylinders, whereas cylinder 2 is slightly lower. Figure 6.4 does indeed show cylinder 2’s torque as the lowest one but cylinder one is not the highest. The only way to know if this is a normal behaviour would be to have a pressure sensor in each cylinder and compute the IMEP. Finally, one noteworthy thing is that the system response during transient opera- 44 7. Discussion tion is very satisfactory. The estimation of the torque remains close to the reference torque and not as noisy as one could have feared. 7.2 Estimation using Fourier series Using the powertrain model with Fourier series and the separation algorithm in Sec- tion 5.1 requires quite a lot of computational heavy steps. In Matlab this method was suitable to try and simulate on old engine data, but was never implemented in the actual EMS, since it was suspected to require too long computational time, and also long time to implement. It was concluded that doing estimation of average combustion torque using the IMEP-model was more interesting for this thesis time frame and goal. Also, average combustion torque (or indicated torque) is one of the most common parameter for controlling the applied torque (requested from the driver). The method itself is however interesting if one would like to estimate the individual torque functions for each cylinder. Interesting values of the combustion process can be estimated by developing a method to convert the torque into e.g. pressure, in order to extract estimates of the PCP. However, this needs further investigation since the crank lever function L(θ) has a singularity at TDC. Also, since the model is relating acceleration of the flywheel (θ̈F ) to torque, the Fourier series have zero mean and their index is n ≥ 1. Therefore, it is not possible to retrieve any informa- tion about the average torque or IMEP. There are also limitations in the model identification conducted in this thesis, since it was only identified on a certain region of engine velocity and torque (operating point 7 and 8). The model is local, which means that it is valid for a certain range of engine speeds and torque (in a window). In [12] the powertrain model is also defined by taking the differences in torque and acceleration, with respect to a nomi- nal torque and acceleration. This gave higher accuracy than achieved in this thesis, but was never implemented. Since the model is local one needs to save parts of the model for different operating regions, which is also a limitation. 7.3 Kalman filter approach The estimations coming from the Kalman filter model looks promising. It is also a robust model since tuning of Rn, Qn, and δ is possible. Since the measured state y(n) = θ̇2 contains noise, it is of big advantage to tune Rn. The noise is caused by different vibrations coming from the engine. Also, the model is not expected to be perfect, since it assumes a stiff crankshaft and an approximated two mass model, which is why tuning Qn is a big advantage. To estimate combustion torque in a real engine, one would need to measure T ∗load(n) with high accuracy, which would improve estimations. In [10] this information was available, and the combustion torque was estimated with very promising results. 45 8 Conclusion In this chapter conclusions are made from the outcome of this thesis. Also, some interesting topics for future work are presented. 8.1 The thesis outcome Estimating the pressure or equivalent combustion torque in all the cylinders of a combustion engine is a recurring problem that has been addressed many times. This thesis presented three different concepts that allow to deduce, from the pressure in one of the cylinders and the angular position of the inertia wheel, the combustion torque in the 6 cylinders. The first method, based on transfer function identification between flywheel acceleration and torque, gave very good results in simulation for catching dynamics and retrieving the whole torque curve for each cylinder. How- ever, this method has not been implemented as a TargerLink model and tested on a real engine due to its complexity and the Master’s Thesis time frame. The sec- ond method derived a very simplified engine model between load torque, flywheel acceleration and combustion torque in order to collect average combustion torques over an important number of cycles. Since the computational cost and complex- ity of this solution is low and compatible with TargetLink as well as the EMS, an implementation was possible. Thus, this method was tested on a real engine and gave good results to determine the average combustion torque coming from the 6 cylinders. The results for the individual torque estimation were not as good as ex- pected and predicted by the simulation. A better calibration would lead to better results. A third method has also been analyzed during this thesis: The estimation of the difference between combustion and load torque curve with a Kalman filter. It is possible to estimate the combustion torque by approximating the load torque to be constant. Also, individual torque contributions was possible to estimate in simulation, by using separation with the help of the ICPS. A more accurate load torque signal would improve these estimations. However, the observer still gave some promising results, and may be suitable for this problem since both noise in model and measurements exists. Due to the time frame of this thesis, it never got implemented in the real engine, but was tested on engine data. 46 8. Conclusion 8.2 Future work There are some directions and recommendations for future work. It could be in- teresting to investigate the IMEP method by correctly calibrating the slope values, checking other operating points and introducing injection faults. To reduce the num- ber of saved slope values, one possible idea is to fit a polynomial between different slopes values at different engine speeds for the same cylinder. This would decrease the number of saved values and potentially allow transient operations in speed. For the transfer function method, the main development axis would be to directly write a C code that is compatible with embedded systems instead of using TargetLink. The Kalman state observer could possibly be implemented by doing approximations to the model, and simplifying the filter updates. Also, one could try to simplify the state space model, assuming non-varying inertia such that An becomes constant. In a real time EMS one would want to tune Q, R and δ in real time and have ac- cess to these parameters. A load torque sensor capable of measuring T ∗load(n) would also be interesting to investigate, since this is needed in order to extract the total combustion torque. 47 Bibliography [1] Saber Fallah. Electric and Hybrid Vehicles - Technologies, Modeling and Con- trol: A Mechatronic Approach. 04 2014. [2] Anton Kjellin and Per-Sebastian Pettersson. Torque estimation from in-cylinder pressure sensor for closed loop torque control. Master’s thesis, Chalmers Uni- versity of Technology, 2017. [3] Marcus Hedegärd. Torque estimation project, Report 1: Initial signal processing. Technical report, 2018. [4] Martin Williams and Ray Minjares. A technical summary of Euro 6/VI vehicle emission standards. The international council of clean transportation, 2016. [5] Maennel Ralf Schiefer Dirk and Nardoni Wesley. Advantages of Diesel Engine Control Using In-Cylinder Pressure Information for Closed Loop Control. SAE Technical Paper, 2003. [6] Eunhwan Kang Kyunghan Min, Jaesung Chung and Myoungho Sunwoo. Indi- vidual Cylinder IMEP Estimation using a Single Cylinder Pressure Sensor for Light-duty Diesel Engines. SAE International, 2014. [7] Marcus Hedegärd. Estimation of Torque in Heavy Duty Vehicles with Focus on Sensor Hysteresis. PhD thesis, Chalmers University of Technology, 2016. [8] Nielsen Lars Kiencke Uwe. Automotive Control Systems For Engine, Driveline, and Vehicle. Springer, 2005. [9] John B. Heywood. Internal combustion engine fundamentals. McGraw-Hill, 1988. [10] Philippe Moulin Michel Castagne Nicolas Petit Jonathan Chauvin, Gilles Corde and Pierre Rouchon. Real-Time Combustion Torque Estimation on a Diesel En- gine. Test Bench Using Time-Varying Kalman Filtering. 43rd IEEE Conference on Decision and Control, 2004. [11] Marcus Hedegärd. Torque estimation project, Report 2: Identification of Fourier Series. Technical report, 2018. [12] Marcus Hedegärd. Torque estimation project, Report 3: Engine model identifi- cation relying on identified Fourier series. Technical report, 2018. [13] Albert Nistor. Individual cylinder torque estimation from angular speed using a dynamic model. Master’s thesis, Chalmers University of Technology, 2007. 48 A Appendix 1 I A. Appendix 1 A.1 Mass torque derivation According to [8], mass torque can be computed with the derivative of the kinetic energy, i.e. Kmass(θ) = 1 2J(θ)θ̇2 , (A.1) Tmass = dKmass dθ . (A.2) Before going any further, it is important to define the considered masses in this model and their action. The assumption made in [8] is that the connecting rod can be divided in two parts: An upper part that is oscillating and a lower part connected to the crankshaft that is in rotation. The center of gravity is the part where we divide the masses. Let us call the distance between the top of the lever (connected to the piston) and the center of gravity lcog. Then, when adding the mass mpiston of the piston to the oscillating one and the mass of the crankshaft mcrank to the rotating one, mosc = mrod lcog l +mpiston , (A.3) mrot = 6mrod ( 1− lcog l ) +mcrank . (A.4) Note that mosc is defined for one cylinder and mrot is for the whole crankshaft. An expression of the inertia J is derived in [8], and gives J(θ1, ..., θncyl) = mrotr 2 +mosc ncyl∑ i=1 v2(θi), (A.5) where v(θi) = ds(θi) dθ . Finally, using (A.1), (A.2) and (A.5), Tmass can be expressed as follows: Tmass(θ) = 1 2 ( θ̇2dJ(θ1, ..., θncyl) dθ + J(θ1, ..., θncyl) dθ̇2 dθ ) = J(θ1, ..., θncyl)θ̈ + 1 2 θ̇ 2dJ(θ1, ..., θncyl) dθ = ( mrotr 2 +mosc ncyl∑ i=1 v2(θi) ) θ̈ + 1 2 θ̇ 2 d dθ ( mrotr 2 +mosc ncyl∑ i=1 v2(θi) ) = (mrotr 2)θ̈ + ( mosc ncyl∑ i=1 v2(θi) ) θ̈ + 1 2 ( mosc ncyl∑ i=1 dv2(θi) dθ ) θ̇2 = (mrotr 2)θ̈ + ( mosc ncyl∑ i=1 v2(θi) ) θ̈ + ( mosc ncyl∑ i=1 v(θi) dv(θi) dθ ) θ̇2 (A.6) II A. Appendix 1 A.2 Engine data Table A.1: Engine data from Lamella 3. Lamella 3 matrix Data Operating point number Nominal velocity [RPM] Nominal Torque [Nm] 1 800 1800 2 800 2100 3 800 2350 4 800 2650 5 900 1800 6 900 2100 7 900 2350 8 900 2650 9 1000 1800 10 1000 2100 11 1000 2350 12 1000 2650 13 1100 1800 14 1100 2100 15 1100 2350 16 1100 2650 17 1200 1800 18 1200 2100 19 1200 2350 20 1200 2650 21 1300 1800 22 1300 2100 23 1300 2350 24 1300 2650 25 1400 1800 26 1400 2100 27 1400 2350 28 1400 2650 29 1500 1800 30 1500 2100 31 1500 2350 32 1500 2650 33 1600 1800 34 1600 2100 35 1600 2350 36 1600 2650 37 1700 1800 38 1700 2100 39 1700 2350 40 1700 2650 41 1800 1800 42 1800 2100 43 1800 2350 44 1800 2650 III A. Appendix 1 Table A.2: Engine data from Lamella 3 with faulty injections. Lamella 3 faulty injection (-10%) Data Operating point number Nominal veloc- ity [RPM] Nominal Torque [Nm] Faulty cylinder 1 800 2200 2 800 2200 1 3 850 2200 4 900 2200 5 900 2200 2 6 950 2200 7 1000 2200 8 1000 2200 3 9 1050 2200 10 1100 2200 11 1100 2200 4 12 1150 2200 13 1200 2200 14 1200 2200 5 15 1250 2200 16 1300 2200 17 1300 2200 6 18 1350 2200 19 1400 2200 20 1400 2200 1 21 1450 2200 22 1500 2200 23 1500 2200 2 24 1550 2200 25 1600 2200 26 1600 2200 3 27 1650 2200 28 1700 2200 29 1700 2200 4 30 1750 2200 31 1800 2200 32 1800 2200 5 A.3 Correlation coefficients Table for the calculated correlation coefficients for the linear IMEP model in Section 3.3. IV A. Appendix 1 Table A.3: Correlation coefficients to choose the windows for the IMEP model. Correlation coefficients Windows OP9 to OP12 OP13 to OP16 OP17 to OP20 OP21 to OP25 OP25 to OP28 Average TDC+6+0 0.9927 0.9683 0.9936 0.9697 0.9925 0.98336 TDC+7+0 0.9949 0.9762 0.9956 0.9752 0.9977 0.98792 TDC+7+1 0.9963 0.9935 0.9977 0.9974 0.9949 0.99596 TDC+7+2 0.9967 0.9941 0.9975 0.9937 0.9921 0.99482 TDC+8+0 0.9973 0.988 0.9968 0.9892 0.9986 0.99398 TDC+8+1 0.9963 0.9975 0.9974 0.9981 0.9945 0.99676 TDC+8+2 0.9964 0.9963 0.9984 0.9964 0.9945 0.9964 TDC+9+0 0.9976 0.9847 0.9956 0.9816 0.9978 0.99146 TDC+9+1 0.9961 0.9952 0.9965 0.9984 0.9954 0.99632 TDC+9+2 0.9963 0.998 0.998 0.9977 0.9962 0.99724 TDC+9+3 0.9847 0.9537 0.9972 0.975 0.9825 0.97862 V List of Figures List of Tables Abbreviations Acronyms Introduction Background Purpose Problem formulation Methodology Scope System overview Combustion cycle in a four-stroke engine Pressure sensor for in-cylinder measurements Angular sensor at the flywheel Cylinder ID Signal processing Physical models Torque model with crankshaft Combustion torque Tcomb Mass torque Tm Friction torque Tf Driveline torque TD Powertrain model Linear IMEP model Torque model utilizing a Kalman filter Model identification Identification of Fourier series Powertrain model identification Linear IMEP model identification Torque estimation Torque estimation using Powertrain model Sensitivity analysis Estimation using Linear IMEP model Torque estimation utilizing a Kalman filter Implementation on engine Experimental set-up and implementation Online estimation using the IMEP model Discussion IMEP estimator Estimation using Fourier series Kalman filter approach Conclusion The thesis outcome Future work Bibliography Appendix 1 Mass torque derivation Engine data Correlation coefficients