Paper received: 04.09.2009 Paper accepted: 25.09.2009 Application of Optimization Techniques to Determine Parameters of the Vibe Combustion Model Ivo Prah1, - Tomaž Katrašnik2 1 AVL-AST d.o.o., Maribor, Slovenia 2 University of Ljubljana, Faculty of Mechanical Engineering, Ljubljana, Slovenia The capability of the optimization algorithms employed for tuning parameters of the Vibe combustion model included in 1-D thermodynamic engine cycle simulation tool AVL BOOST is analyzed while simulating high pressure phase of the in-cylinder process. The objective of the analysis includes the ability check of the optimization methods to determine the parameters of the Vibe combustion model within the tolerance range needed to set up a high fidelity model and to reach accuracy threshold of the commonly applied analytic approach for determining combustion parameters. The employed method presents an influence of different merit functions and constraints on the accuracy of the results and proposes the methodology for their quality analysis. It can be concluded that the accuracy of the results calculated by combustion parameters determined by the optimization techniques reaches accuracy of the results calculated by a special analytically based software tool. © 2009 Journal of Mechanical Engineering. All rights reserved. Keywords: internal combustion engine, high pressure phase, optimization, Vibe combustion model 0 INTRODUCTION Software tools for thermodynamic modeling of internal combustion engines (ICEs) [1] to [4] have become indispensable for developing and optimizing the ICEs. Due to hardware performance constraints and due to computational time limitations, commercial simulation tools for modeling the complete internal combustion engine (ICE), including intake and exhaust manifolds, do not incorporate 3-D methods, or just take advantage of using coupling to the 3-D software to resolve phenomena in a specific selected component. Therefore, although the models are based on mechanistic models, they incorporate many tuning parameters needed to be adjusted to set up the high fidelity ICE simulation models. Tuning parameters have a clear physical interpretation and are, therefore, adjusted within a meaningful range characteristic for particular phenomena. Tuning parameters are employed due to the inability of the models to fully capture 3-D fluid flow and heat transfer phenomena, which are additionally coupled with mass transfer phenomena and chemical kinetics mechanisms during fuel injection, evaporation and combustion phase. Adjusting of a tuning parameter to meet the accuracy threshold required by the customers (agreement of measured and calculated engine performance should typically be in the range of 1 to 3%) is very time consuming and presents a considerable part of the work load of the whole project. Therefore, employment of optimization methods seems a promising approach for determining tuning parameters. It has already been shown that optimization methods are suitable for determining optimized configurations of components, e.g. manifold geometry [5] and [6], shape of injection rate [7], exhaust gas recirculation (EGR) rates and multiple injections [8], injection timings [9], valve lift profiles and timings [6] and [10] to [13], as well as constants of heat transfer model [14]. These analyses either use a calibrated ICE simulation model as a starting point for evaluating optimized configuration [5], [6], [10], [12] and [13] or tune specific sub-model by optimization methods and compare simulation results to the experimental data [14]. The number of tuning factors (parameters) required to calibrate the simulation of the conventional turbocharged compression ignited ICE with variable turbine geometry (VTG) and EGR system, can vary in a range between 15 and 30. This might be inconvenient for an attempted employment of optimization methods to calibrate simulation model since the effectiveness of optimization techniques decreases with an increase the number of optimization factors. Therefore, it is proposed to divide the simulation model of ICE to feasible sub-systems and carry *Corr. Author's Address: AVL-AST d.o.o., Trg Leona Štuklja 5, 2000 Maribor, Slovenia, ivo.prah@avl.com out the calibration process in two steps. The goal of the first step is to tune factors of sub-systems, while the goal of the second one is to adjust the reduced number of the most dominated parameters while optimizing the whole simulation model. Fig. 1 shows the proposed division of ICE simulation model into subsystems. In the proposed study, the high pressure phase (HPP) of the in-cylinder process is analyzed (subsystem SS1_6 in Fig. 1) with the objective to analyze the ability of the optimization methods to determine the parameters of the combustion model within the tolerance range needed to set up a high fidelity model capable of reaching the above accuracy constraint. Therefore, engine performance data, in-cylinder combustion parameters, rate of heat release (ROHR) and pressure traces calculated by combustion parameters derived through optimization methods are compared to the measured and analytically derived data and to the data calculated by combustion parameters derived from analytically evaluated ROHR (commonly applied approach). Combustion parameters are determined through optimization methods by coupling of the optimization techniques with thermodynamic engine simulation model [1]. Analysis is performed with the design of experiments (DoEs) technique and NLPQL optimization method based on quadratic programming. Methods are included in software package [15]. Additionally, the influence of different merit functions on the determined tuning parameters is analyzed to reveal their substantial influence on the results. In order to ensure high interpretative value of the analysis, a single Vibe combustion model was selected. It includes combustion duration (CD) and shape parameter m [16]. Moreover, the start of combustion (SOC) was also assigned as an optimization parameter, since it significantly influences engine performance. 1 HIGH PRESSURE PHASE AND SIMULATION MODEL Cylinder process of ICE thermodynamic cycle includes an air plus, optionally, recirculated exhaust gas induction, compression, and combustion together with an expansion and exhaust phase. High pressure phase (HPP) comprises the cylinder process during the time period from the intake valve closing (IVC) to the exhaust valve opening (EVO). Fig. 1 presents the proposed division of the ICE simulation model into subsystems marked with SS1_6, SS2_6, ..., SS6_6 for the case of the available measurement data labeled at positions (0,1, ..., 8). The division of the ICE model should be driven depending on the available measurement data. Quantities p1, T denote static pressure and temperature at a certain position, whereas subscript i denotes ref refers ambient, 1 compressor inlet, 2_1 compressor outlet, 2_2 intercooler outlet, im intake manifold, 3 turbine inlet and 4 conditions at turbine outlet. m and air m^, represent air and fuel mass flow rate, respectively, Pen brake power and n rotational speed of ICE. SS1_6 covers cylinder with HPP, SS2_6 inlet orifice - compressor inlet, SS3_6 compressor outlet - intake manifold inlet, SS4_6 turbine outlet - exhaust orifice, SS5_6 turbocharger with compressor/turbine inlet and outlet, SS6_6 intake manifold inlet, engine block with cylinders, exhaust manifold outlet (up to turbine inlet). The required condition enabling evaluation of the HPP in the cylinder e.g. with [1] or [17] is initial thermodynamic state (ITDNS) in the cylinder at the time of IVC. During HPP in-cylinder variables are not influenced by other engine components, since the intake and exhaust valves are closed, and thus only ITDNS is influenced by other engine components (e.g. intake/exhaust system during gas exchange process). ITDNS is determined either by the cylinder pressure at IVC ( pcyl _ JVC X mair, m fuel and residual gas concentration (RGC), or by the temperature, pressure and gas compositions at IVC. Gas compositions are defined by air to fuel ratio, fuel vapour and combustion products concentration. The application of the measured data for Pcyl _ JVC , mair, m fuel aM estimated RGC for HPP analysis, excludes the impact of the other engine components. RGC cannot be measured. It is, therefore, determined either by specialized software tools, e.g. [1] and [18] where complete ICE cycle has to be simulated, or based on the user experience. For presented work RGC values have been determined based on the user experience. Fig. 1. Conceptual division of ICE into subsystems The objective of the HPP thermodynamic analysis is to determine appropriate ROHR resulted from in-cylinder combustion process based on known (measured) cylinder pressure trace and ITDNS. The determination of the ROHR from the measured data is commonly performed by the analytical method via specialized software tools, e.g. [17] and [18]. This gives ROHR as a function of the CA rotation (table). Besides, [17] and [18] also evaluate SOC and Vibe combustion model parameters (CD, m) from ROHR by analytical approach. Optionally also, a direct approach to tune parameters of a certain combustion model included in [1] might be applied. Thereby, sufficiently good agreement in simulated and measured pressure represents the target for parameters tuning of certain combustion model. Tuning of the combustion model parameters through the direct approach can be realized by performing simulations with an activated check box start high pressure (SHP) within the cylinder element of the simulation model presented in Fig. 2 either with an application of trial-and-error method or by optimization techniques. The latter has been carried out in the presented work. The input factors of the selected combustion model being the subject of optimization tuning have to be assigned as variables (parameters). [1] includes several physical models for ROHR prediction whose number of tuning parameters depends on the complexity of certain model. A single Vibe model was applied in the analysis, since it enables a comparison between parameters evaluated with [17] and the same set of parameters determined directly with an optimization techniques. Fig. 2. Simulation model created with [1] for combustion model parameters tuning A simulation model consists of a cylinder element (C1), attached pipes 1, 2 that model intake and exhaust port, and internal boundaries (IB1, IB2) that involve static boundary conditions (e.g. pressure, temperature) at a certain position in the intake and exhaust port. ITDNS for HPP analysis is determined in SHP input dialog. The required input data for IB1, IB2 conditions and for pipes 1 and 2 attached to the C1 do not influence the simulation results, but they have to be specified due the layout structure of [1]. The measured data of MAN turbocharged compression ignited (TCI) diesel engine were used for the presented work. Basic engine geometry and valve timings are listed in Table 1. Table 1. Engine geometry and valve timings of Bore [m] 0.108 Stroke [m] 0.125 Connection rod length [m] 0.1825 Compression ratio [-] 18 Piston pin offset [m] 0.0005 IVO [CA] 17 BTDC = 343 IVC [CA] 33 ABDC = 573 EVO [CA] 57 BBDC = 123 EVC [CA] 25 ATDC = 385 Number of valves/cylinder 2 Number of cylinders 6 Cylinder displacement [l] 1.145 Engine displacement [l] 6.87 Valve timings are specified according to firing top dead center (FTDC) determined in [1] with 0 CA. 2 CYLINDER BALANCE EQUATIONS In a 0-D single zone model, balance equations of mass, energy and species concentrations form a sufficient set of dependent variables, and time or equivalently crank angle rotation represents the independent variable. Fig. 3 shows a physical model and state variables of the cylinder and intake/exhaust ports. Fig. 3. Schematic of the considered four-stroke cylinder connected to the intake and exhaust manifold including state variables and fluxes The framework of balance equations is laid out in a general way enabling a consideration of an arbitrary number of species W =[ Wl K ^ f (l) In the presented analysis species vector W, represents burned fuel (FB), combustion products (CP) and fuel vapor (FV). W = [ WFB , WCP , WFV ] (2) whereas species concentration of air is derived by: Wair = 1 - WCP - WFV • (3) In deriving governing equations the dependency of internal energy (u) u = u{T, p, w1 k wm+1) (4) and specific gas constant (R) R = R(T, p, w K ww+1) (5) on temperature (T), pressure (p), and species concentrations (wFB, wCP, wFV, wair) are considered. A revised ideal gas equation [19]: ZM pV = m-T = mRT M (6) adequately captures deviations of the real gas from the ideal one in the range of temperatures and pressures characteristic for the in-cylinder processes of ICEs; herein V is volume, m is mass, Z is compressibility factor, ^ is universal gas constant and M molar mass. The mass balance equation is: dm dm j dm dç = 1 inj dm bb dç dç dç (7) where (p is the independent variable representing CA rotation, and indexes inj injection and bb blow-by. The change of mass depends on the fluxes through the attached valves and the amount of the injected fuel. The rate of change of species conservation comprises species concentration variation due to mass transfer and species concentration variation due to: dwk dç naves dmj Ç ( • '- , * +- dw. m d ç where index d denotes I cylinder for dm < 0 d = | plenum for dm > 0 (9) The first law of thermodynamics for an open system equals where d^ = dß - pdV + dH dQ = dQc + dQHT nvalves dH = X hddmj + hbbdmbb (10) (11) (12) and index d is defined by eq. (9). Term dQC in Eq. (11) models rate of heat release (ROHR) while dQm in Eq. (11) models heat flux from the gas within the combustion chamber to the surrounding walls. ROHR according to Vibe model [16] is determined by the following equation dQC = q d^L = df F df LHVminj— (m + 1) ^^ e ^ Af n Af ^ Af where QF = minjLHV denotes fuel energy, minj mass of the fuel injected into the cylinder, LHV lower heating value of the fuel, a completeness of combustion [16], Af combustion duration [16], f0 start of combustion and m shape factor [16]. A general equation for heat transfer (to the walls) evaluation considered in [1] has a form (13) dqpht- = aia (p)htcf(t (p)- twii) dp (14) where A, represents the surface area of the surrounding walls, e.g. (cylinder head, piston liner), aw heat transfer coefficient [20], HTCF heat transfer multiplier, Twi wall temperature of the surface A;. Index i denotes cylinder head, piston or liner. The thermodynamic engine cycle simulation model [1] calculates in-cylinder temperature and pressure due to mass and enthalpy flows, combustion and piston kinematics. Therefore, Eq. (10) is algebraically reformulated to be explicit in temperature leading to (derivation was done in analogy to [21] where additionally a more general dependency of the internal energy (Eq. (4)) and the gas constant (Eq. (5)) on species concentration is taken into account) T = B dp m dQ dH . v — +-+ (Km -1) dp dp p ŠL -K + KTRm)) dp dp ( -m where dR du KmT + dwk dwk dwk dp K = du v h - P dR R dp (15) (16) and B = ■ 1 KmR 1+T dR R dT du ~8T (17) The direct approach of combustion model parameters tuning (SOC, CD, m) via the optimization technique is based on the coupling of the optimization techniques with thermodynamic engine simulation model [1] based on the Eq. (15). Combustion analysis is based on the inverse procedure where ROHR is calculated from measured pressure trace, piston kinematics, heat transfer, and ITDNS data. Therefore, Eq. (4) is algebraically reformulated to dß + dH dp dp T du u ~A~8T 1 + - 1 du RA dT dm — du dR dV >-- dp m where: du dwk du dp dT dw, k RA dw, — TC du pA dT (18) dp = 0 A = 1 + T dR R dT (19) and C = 1 - p cR R cp (20) Analytic determination of the ROHR from the measured data is performed via specialized software tool [17] based on Eq. (18). 3 COMBUSTION MODEL PARAMETERS TUNING WITH OPTIMIZATION METHODS The proposed study was investigated with the combination of DoE and NLPQL optimization method included in the optimization software package [15]. DoE method was used to determine feasible design space of tuning factors (parameters) and to determine a good starting point for a further optimization process employing NLPQL algorithm. NLPQL [15] and [22] is gradient based local optimization method, which requires a good starting point close to the global minimum or the merit function must be convex featuring only one local minimum. The purpose of DoEs is to gain a maximum amount of information about the simulation model response, while minimizing the number of simulations. Full factorial design (FFD) [15] was implemented for factors screening. Factors by which minimum merit function was attained have been used as a starting point for optimization process. Based on former analyses it was decided to perform optimization with NLPQL instead of e.g. genetic algorithm, since combining DoEs and NLPQL methods requires significantly less iterations (5 to 10 times) to attain the similar accuracy of the results. The following steps were performed in tuning the Vibe combustion model parameters using DoEs and NLPQL optimization algorithms: 1. Preprocessing measured cylinder pressure trace and analytical evaluation the ROHR. [17] includes options to filter the measured cylinder pressure trace in order to compensate for noise, errors or inaccuracies. Besides top dead center (TDC) and pressure offset error can be detected. A preprocessed measured cylinder pressure trace was later used for an analytical evaluation of the ROHR based on Eq. (18) considering the measured ITDNS as well for the optimization process. Besides, the calculated Vibe parameters (m, CD), SOC, and CA rotations at 5% (ps%mfb) and 90% (p90%mfb burnt mass fractions were extracted from [17]. 2. Geometric parameterization. The desired tuning factors of Vibe combustion model (m, CD) and SOC used in the simulation model presented in Fig. 3 were assigned to the variables enabling their variation during optimization. 3. The simulation of the response and the definition of the optimization objectives. With the initial set of tuning factors (m, CD, SOC) one simulation was performed by [1] and afterwards merit and constraints functions were determined by [23]. Single and multi-objective optimization problems were dealt with equality constraints. Besides, selected merit functions were analyzed also as an unconstrained optimization problem with the aim to reveal the influence of constraints. The optimization problem is expressed with the following equations: min F(x) Subject to: hi(x) = 0 i = 1 .. 4 x = [xi, x2, x3] = [m, CD, SOC] xk < xk < xku k = 1 .. 3 (21) whereas F(x) represents single or multi-objective merit function, hi(x) equality constraint and x vector of tuning factors. Four different single (F1,2,4,5) and two multi-objectives merit (F3,5) functions have been investigated: pevo F1(x) = i(cyl _ s (x,p Pivc - pcyl _ m ((p)YdV F2(x) = J (x,p) (22) DUR ' V dP Pivc dQP _ a dp (P) dp (23) 2 F3(x) = F1(x) + F2(x) 1 p0% MFB F4(x) = c- ¡(fey' _ -(x' p) (24) " Pcyl _ n( d -80 c ^ -60 C -40