Strojniški vestnik - Journal of Mechanical Engineering 59(2013)4, 223-236 © 2013 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2012.668 (1.01) Original Scientific Paper Received for review: 2012-06-19 Received revised form: 2012-12-28 Accepted for publication: 2013-01-18 On the Convergence, Stability, and Computational Speed of Numerical Schemes for 0-D IC Engine Cylinder Modelling Tomaž Katrašnik1* - Henrik Aleš Schuemie2 - Johann C. Wurzenberger2 1 University of Ljubljana, Faculty of Mechanical Engineering, Slovenia 2 AVL List GmbH, Advanced Simulation Technologies, Austria The development of real time capable 0-dimensional internal combustion engine models places high demands on convergence, stability, and computational speed of the applied numerical methods. The cylinder model represents the crucial element in attaining high computational speed and accuracy of results. A basic example comprising a single cylinder connected to two plenums is analysed with different numerical schemes in order to reveal methods effectively associating accuracy requirements with computational time constraints. The integration performance to solve a system of coupled ODEs was compared for explicit Euler and explicit fourth order Runge-Kutta schemes, as well as for multi-step methods including backward differentiation formulas and Adams-Moulton formulas. The performed analysis emphasizes two major points. First, the numerical accuracy of integration schemes differs significantly at equal computational effort revealing the necessity of selecting an adequate scheme for a specific task. Second, the comparison of integral engine parameters (e.g. indicated mean effective pressure, mean engine torque), calculated by different methods, with a numerically assumed exact solution should not be used as an estimate for the convergence and stability of the applied numerical approach, since good agreement in integral parameters does not imply good agreement in cycle resolved traces of thermodynamic variables. This paper provides clear guidelines for selecting the appropriate numerical integration methods with respect to the intended application. Analyses are also based on innovative test examples. Finally, a comparison of numerical and experimental in-cylinder pressure traces is shown for a series production engine confirming the applicability and accuracy of the cylinder model. Keywords: internal combustion engine modeling, integration schemes, convergence, stability, computational speed 0 INTRODUCTION The steadily increasing demands on vehicle fuel economy, emissions, and performance are a major driving force in developing and optimising future powertrains. In order to meet customers' and legislative requirements it is of particular importance not only to simulate and optimise different powertrain elements, but also to develop simulation models capable of modeling the performance of the powertrain system as a whole. The latter is of particular importance for developing and testing components in a virtual environment, as well as for accompanying the conventional development steps from design phases in the office to the application and validation phases on the test bed. Test bed validation, development and calibration of components in the hardware-in-the-loop (HiL) environment and modelbased engine controlling require real-time capable simulation models. Moreover, computational speed is also of considerable importance when simulating complete emission test cycles where engine and vehicle dynamics are simulated on time scales of 103 s. Recently, 0-dimensional (0-D) crank-angle (CA) resolved engine simulation models, e.g. [1] to [4], have emerged as a promising approach for real-time engine simulations of internal combustion engines (ICE) and hybrid powertrain of heavy-duty engines [5]. Unlike higher resolution 1-D, e.g. [6] to [9], and 3-D models, e.g. [10] and [11], 0-D models can comply with realtime constraints and provide a higher level of detail compared to mean-value engine models, e.g. [12] and [13]. The latter models reach their limits when future control functions need to be taken into account, which incorporate CA resolved events such as cylinder pressure based control functions, electronic control unit controlled fuel injections, and variable valve timing. Mean-value engine models are not able to model these phenomena predictably since they rely on experimental data that are used to tune or train the model as presented in [14] for camshaft variability. Besides computational speed, the convergence and stability of numerical solution algorithms are the prerequisites for the successful application of 0-D CA resolved engine simulation models. This topic has not yet been systematically analysed in the literature and therefore the presented paper aims to provide a study of the convergence, stability, and computational speed of different explicit and implicit integration schemes. The cylinder model represents the crucial element to attaining high computational speed and accuracy of results, since the time variation of the thermodynamic variables is most pronounced in the cylinder due to piston kinematics and combustion. It is common practice to validate the accuracy of simulation models by comparing numerical and experimental results. This indeed gives an indication of the consistency *Corr. Author's Address: University of Ljubljana, Faculty of Mechanical Engineering, Aškerčeva 6, SI-1000 Ljubljana, Slovenia, tomaz.katrasnik@fs.uni-lj.si 223 of the virtual model and the hardware being simulated. However, it is well known that 0-D models incorporate many parameters included in combustion and heat transfer models that can be tuned with the aim of achieving agreement between numerical and experimental results. If numerical stability and convergence of the model are not validated separately, the agreement between numerical and experimental results alone cannot be understood as sufficient validation. This is particularly true as long as model parameters are tuned with the objective of achieving this agreement, since numerical errors might thereby be tuned out of the model. Agreement between numerical and experimental results does therefore not necessarily imply convergence and stability of the solution, since the results may be non-converged and even locally non-stable. The drawbacks of such models are often revealed during transient engine simulation where significant discrepancies may also arise between the numerical and experimental results, since it is generally not possible to tune the model for the complete transient operating region. These phenomena are particularly discernible with turbocharged engines exhibiting turbocharger lag, as evident in [2] where very large integration increments are applied. Fig. 1. Schematic of the considered four-stroke cylinder connected to the intake and exhaust plenum including state variables and fluxes Moreover, integral engine parameters (e.g. indicated mean effective pressure (IMEP), brake mean effective pressure, mean engine torque, effective fuel consumption) are often used to validate numerical models. The accuracy threshold for validating high fidelity ICE models is typically 1%. The presented paper shows that the agreement of integral engine parameters (results for IMEP are shown, since other integral engine parameters are directly dependent on IMEP) calculated by the different methods and the numerically assumed exact solution might often easily comply with these accuracy requirements. However, the traces of thermodynamic parameters may be non-converged and even locally non-stable. The same consideration can also be applied to the comparison of numerical and experimental data. It is thus shown that the agreement of integral engine parameters should not be the only measure for assessing the accuracy of the model. This paper presents an analysis of convergence, stability, and computational speed of integration schemes for 0-D modelling of ICE cylinders in a simple model. A single four-stroke cylinder (Fig. 1) connected via two valves to two plenums (intake and exhaust) of infinite inertia (all thermodynamic variables are constant) is used. The model considers fluxes through the attached valves, the flux of the injected fuel, combustion of fuel, heat transfers to the combustion chamber walls and the walls of the ports, and the power given to the piston. This simplified test case is selected since it excludes other engine components, thus enabling the investigation of numerical phenomena relevant only to the low and high pressure phase of the in-cylinder process. 2 GOVERNING EQUATIONS The framework of equations was laid out in a general way enabling consideration of an arbitrary number of species: W = [[...wN (1) where wi represents the mass fraction of species i. In the present analysis, the species vector, W, represents the burned fuel (FB), combustion products (CP), and fuel vapor (FV). The species concentration of air is derived by: Wair = 1 - WCP - WFV . (2) wFB gives the concentration of the fuel that was burned, i.e. converted from FV to FB, whereas wCP gives the concentration of the corresponding combustion products, i.e. burned fuel plus air that was used in combustion. Therefore, it is possible to evaluate the air-fuel ratio of CP for the particular fuel, which is mandatory for evaluating adequate gas properties. The governing equations were derived considering the dependencies of gas properties, i.e. internal energy (u), the specific gas constant (R), and all partial derivatives of both properties, on temperature (T), pressure (p), and species concentrations (wFB , wCP , wFV . wair). Gas properties are given by input maps taking into account chemical equilibrium considerations. A database was prepared for different fuels and fuel blends. In the analysed case, the properties of Diesel fuel are used. A revised ideal gas equation [15]: pV = mT = mRT, M (3) adequately captures the deviations in the real gas from the ideal one in the range of temperatures and pressures characteristic for the in-cylinder processes of ICEs; V is the volume, m is the mass, Z is the compressibility factor, ^ is the universal gas constant, and M the molecular weight. The mass balance equation is: dm "n^e dmj dminj (4) where % is the independent variable representing either time (t) or CA rotation The change in mass depends on the fluxes through the attached valves and the amount of injected fuel. The rate of change of species concentration comprises the species concentration variation due to mass transfer and species concentration variation due to chemical kinetics, i.e. combustion. It follows that: dwk dwk, "vo^ dm ,GE - Wk )) , dm. W ^ - W „., \)- + ((J„j- Wkcy! ) ml where index GE denotes: I cylinder for dm < 0 GE = plenum for dm > 0 (5) (6) where the positive value of dm indicates that mass flows into the cylinder as given in Eq. (4). The first law of thermodynamics for an open system is: dU _ dQ _ dV dH P + ' (7) where dQ dQc , dQH Č4 Č4 dH dmj dmn h 1 (8) (9) and the index GE is defined by Eq. (6). After algebraic manipulation of Eq. (3) in differential form, and insertion together with Eq. (4) into Eq. (7), Eq. (7) was algebraically reformulated to be explicit in temperature. This derivation was done in analogy to [16] by considering a more general dependency of the internal energy and the gas constant on species concentration is taken into account leading to: dT B dQ dH , , dV — +-+ (Km -1)»--- - -(u + KTRm))-, KmT — + -9U dw. dw, dw. k y d^ (10) where: K = du and B = dpv ^ _ p dR R dp 1 (11) KmR 1 + T dR \ du (12) R dT I dT The combustion is modeled according to the single-zone Vibe model [17]. The heat transfer coefficient to the walls of the combustion chamber is modeled according to the model proposed by Hohenberg [18]. The mass flow rates through the intake and exhaust valves are calculated by the orifice equation (see e.g. Ref. [19]). 3 SOLVING THE SYSTEM OF COUPLED ODEs 3.1 Solution Vector The 0-D ICE cylinder model is described by the system of coupled ODEs described above. An overall solution vector can be written in the following form: y = [m,T, Wi... ww] (13) It is common practice to solve the system of coupled equations (i.e. Eqs. (4), (10) and N Eqs. (5)) 1 + m as an initial value problem approaching a convergent solution by advancing an independent variable. This is given by: ^ = f(l;,y), y^ ) = y0, y e RN+2, (14) where % represents either time (t) or CA rotation 3.2 Integrations Schemes The convergence, stability, and computational speed of integration schemes for 0-D modeling of the ICE cylinder (Eq. (14)) are analysed for the following integration schemes: Explicit Euler scheme [20], Explicit fourth order Runge-Kutta scheme [20], Backward Differentiation scheme [21], and Adams-Moulton scheme [21], where the first two schemes are constant step schemes and the last two are variableorder, variable-step schemes. The Euler scheme features a very high computational speed for a selected integration increment and it is therefore already applied in realtime 0-D ICE models presented in literature, e.g. [2]. However, according to [20], it is not recommended for any practical use. The Runge-Kutta fourth order scheme is known as a robust, widely applied scientific and engineering tool. The Backward Differentiation scheme and Adams-Moulton scheme are included in the CVODE solver package [21]. Variable-order, variable-step multistep CVODE methods are organized to monitor internal consistency. Numerical errors are controlled by adaptive step sizes within a global integration increment. The CVODE package offers the Adams-Moulton formulas for non-stiff problems. For stiff problems, Backward Differentiation schemes are available in the CVODE package. The governing equations for solving thermodynamic variables in the cylinder are nonstiff if physically meaningful input data are used. Thus, all four schemes could be applied in the analysis and compared for performance and robustness. The CVODE solver offers the choice of either functional iteration, suitable only for non-stiff systems, and various versions of Newton iteration applicable to both stiff and non-stiff systems. Both iteration methods were therefore applied in the present analysis and thus the performance of all combinations of non-stiff and stiff numerical schemes and iteration methods were examined. Six solvers are applied to solve cylinder balance equations. The following abbreviations are used in the text below and in the Figures: • Euler, • RK4 - Runge-Kutta fourth order, • BDF, NWT - Backward Differentiation scheme, Newton iteration, • BDF, FCN - Backward Differentiation scheme, functional iteration, • AM, NWT - Adams-Moulton scheme, Newton iteration, • AM, FCN - Adams-Moulton scheme, functional iteration. The Newton iteration requires the solution of linear systems. The selection of the method used to solve a linear system does not noticeably influence the results, since a system of only 5 equations describing a single cylinder engine is considered. In all CVODE schemes, equal methods to solve the linear system and equal iteration parameters were applied. The Euler and RK4 schemes are integrated on a CA rotation basis. With this approach, similar numerical accuracy can be achieved for different engine operating conditions. However, computational time increases with increasing engine speed. CVODE schemes were integrated on a time basis. The tradeoff between numerical accuracy and computational speed was analysed by varying integration increments of Euler and RK4 schemes, and by varying tolerances and global integration increments of CVODE solvers. The main mechanism influencing numerical accuracy and computational speed of CVODE solvers is the definition of tolerances, whereas the global integration increment does not significantly influence the results. This is because of the fact that the solver internally reduces its integration increments until local error tests are satisfied. The tolerances of the CVODE solver were given by relative and absolute tolerance values. The first were given as scalar values and the latter were provided in vector form. Components of the solution vector, Eq. (13), evaluated for cylinders of automotive engines, feature very distinct orders of magnitude when evaluated in SI units, i.e. O(m) = 10-3, O(T) = 103 and O(w,) = 1. The absolute tolerance vector was therefore defined by: TOLabs = O (m) O (T ) O (w) X TOLre (15) where Eq. (15) represents an efficient handling of tolerances, since only TOLrei is varied to analyse the trade-off between numerical accuracy and computational speed. Computational speed is inversely proportional to TOLreh Table 1. Basic cylinder data and parameters Bore [m] 0.085 Stroke [m] 0.094 Connection rod length [m] 0.1505 Compression ratio [-] 17 Piston pin offset [m] 0.002 IVO [deg. CA] 17BTDC = 703 IVC [deg. CA] 33 ABDC = 213 EVO [deg. CA] 57BBDC = 483 EVC [deg. CA] 25 ATDC = 25 Number of valves [-] 2 Vibe parameter m [ -] 0.7 Vibe parameter A [-] 6.901 Combustion duration [deg. CA] 70 4 RESULTS The assembly of the considered single cylinder engine is shown in Fig. 1. An internal mixture preparation compression ignition, i.e. Diesel, engine (see Table 1) was chosen. Besides temperature and concentration variation, the in-cylinder mass also varies during combustion due to direct fuel injection, which therefore represents a more demanding test for the accuracy of numerical schemes. Two different engine speeds are analysed to explore the characteristics of the different integration schemes. The selected speeds cover the lower and upper speed boundary of automotive engines. In order to ensure a reliable comparison between the different integration schemes, special emphasis is placed on a converged solution of cylinder traces. Thus, for Euler and RK4 schemes, the simulation was cyclically repeated until the differences in the last two cycles were below the round-off of machine arithmetic for all analysed parameters. Likewise, for the CVODE schemes, the simulation was cyclically repeated until the differences in the last two cycles are within the specified tolerance limits and it can therefore be assumed that a steady-state solution was reached. Traces of selected thermodynamic variables were analysed first to expose specific characteristics of the particular schemes. Results of these schemes Fig. 2. Pressure traces for an engine speed of 1500 rpm and an IMEP of 1.92 MPa; a) entire cylinder cycle, b) Ap for the entire cylinder cycle, c) and d) low pressure phase of the cylinder cycle were compared to the numerically assumed exact solution (denoted "exact solution" in the subsequent text) calculated using the BDF, NWT method, with TOLrei = 1012 being an unreasonably small tolerance resulting in excessive computational times. For a comparison of different numerical integration schemes applied to different engine operating points, it is necessary to condense the cycle resolved results to a clearly traceable set of lumped parameters. Therefore a non-dimensional root mean square (NDRMS) and a real time (RT) factor were introduced. NDRMS describes a relative difference in a considered variable (i.e. pressure, temperature, mass, etc.) calculated with a chosen integration scheme and the exact solution throughout the complete engine cycle. The NDRMS is defined by: NDRMS (X ) = 1 dato \ 2 ^ (Xi - Xi,exact ) X (16) where X denotes an arbitrary cylinder variable, ndata denotes number of the data within one engine cycle, exact is the exact solution, and Xm represents the maximum value of the considered variable. The NDRMS value is always calculated for the complete cylinder cycle comprising both the gas exchange and the high pressure phase. In order to define a measure for the computational effort of the different numerical schemes an RT factor is applied. RT is given by: RT = (17) physical where t computation is computational time obtained on a standard linux machine with an Intel Xeon 3.00 GHz processor, and tphysicai represents the considered physical time of the investigated engine cycle. It is obvious that efficient integration schemes simultaneously feature both a low NDRMS value and a low RT value. 4.1 Analysis of Cycle Resolved Results The accuracy of integration schemes by means of pressure traces is analysed first. Pressure is not directly calculated by the integration scheme (Eq. (13)) and therefore it represents a good indication of numerical accuracy, since it is influenced by all components of the state vector via Eq. (3). The pressure is also the driving force for mass, species, and enthalpy transfer during the gas exchange phase, and therefore is the most significant parameter for determining numerical accuracy. Additionally, the pressure trace determines the indicated work per cycle via the indicated mean effective pressure (IMEP), which is related to the engine power (P). It is given by: (18) where | pdV is the indicated work per cycle, and Vd cycle dispiacement represents the engine displacement. IMEP is numerically integrated by a trapezoidal scheme considering the pressure at the beginning and at the end of the integration step. Table 2. Boundary conditions for Test case 1: engine speed 1500 rpm and IMEP = 1.92 MPa Plenum T [K] p [Pa] Wfb [-] Wcp [-] Wfv [-] Intake 320 220000 0 0 0 Exhaust 320 200000 0 0 0 Fig. 2 shows pressure traces for an engine speed of 1500 rpm, IMEP = 1.92 MPa. The boundary conditions in the intake and exhaust plenum, which correspond to a turbocharged engine, are given in Table 2. In all graphs, 0 deg. CA corresponds to the top dead center (TDC) of the gas exchange phase. The pressure traces are plotted for different numerical schemes featuring a low RT factor (e.g. RT<0.1 for engine speed 1500 rpm) and for the exact solution. The time steps of the two explicit schemes (Euler, RK4) are chosen to keep the computational effort at a comparable level (RK4 with 1 and 4 deg. CA and Euler with 0.25 and 1 deg. CA). Additionally, the results of the Euler scheme with 2 deg. CA are shown to point out the significantly reduced accuracy. The results of the AM, FCN scheme with TOLrel =10 5 (Eq. (15)) and a global time increment corresponding to 1 deg. CA are shown, since it features both low RT and the best NDRMS vs. RT ratio out of the schemes given by the CVODE package for this example. Fig. 2a shows that integration schemes featuring low RT values and consequently lower numerical accuracy provide a reasonable agreement for the overall in-cylinder pressure trace. A more profound analysis is therefore necessary to reveal the characteristics of the particular schemes. Figs. 2c and d show the low pressure, i.e. gas exchange, phase of the cylinder cycle. It can be observed that results of the RK4/1 deg. CA scheme and of the AM, FCN n t / TOLrel = 10-5 scheme coincide very well with the exact solution (within the line thickness), whereas all other schemes produce significant deviations. These deviations strongly influence mass, species, and enthalpy transfer during gas exchange and therefore imply nonconverged solutions. This non-convergence additionally compromises the accuracy in the subsequent integrations steps. Stability issues are observed for all Euler schemes and for the RK4/4 deg. CA scheme. It can also be observed that the Euler/0.25 deg. CA scheme featuring similar computational times as RK4/1 deg. The CA scheme performs worse during the gas exchange period. This is mainly attributed to the low order of the Euler scheme. The difference is even more discernible in Fig. 2b, which shows the absolute differences between the pressure traces of the applied numerical schemes and the exact solution (Ap = | pscheme - pexact |). It can be observed that Ap reaches its maximum value in the vicinity of the firing TDC, i.e. 360 deg. CA, since deviations in thermodynamic variables at intake valve closing (IVC) result in increasingly pronounced pressure deviation during the compression stroke. A comparison of the results in Figs. 2a and b leads to the conclusion that pressure traces as shown in Fig. 2a cannot be used to derive reliable measures of the quality of different numerical schemes. Out of 6 applied schemes, the RK4/1 deg. CA scheme shows the smallest deviation from the exact solution in the vicinity of the TDC, whereas the RK4/4 deg. CA and the Euler schemes with 1 and 2 deg. CA on the other end show significant and non-negligible errors. Fig. 3. Mass of combustion products (mCP = wCPm) in the cylinder for an engine speed of 1500 rpm and IMEP = 1.92 MPa Deficiencies in convergence of particular schemes are also analysed in Fig. 3, showing the mass of combustion products (mCP = wCP m) in the cylinder for the Test case 1 (Table 2). wCP starts decreasing at intake valve opening (IVO) for the exact solution, since the in-cylinder pressure is larger than the pressure in the exhaust plenum (p > pE) during the exhaust valve opened period and the in-cylinder pressure is smaller than the pressure in the intake plenum (p < pj) during the initial sequence of the intake valve opened period (as shown in Figs. 2c and d) resulting in an inflow of fresh air at IVO. wCP decreases from IVO up to 180 deg. CA, when a flow reversal across the intake valve occurs, since p > pI as shown in Fig. 2c and d. On the other hand, mCP decreases until exhaust valve closing (EVC) and is constant thereafter until 180 deg. CA, since only fresh air flows into the cylinder due to the intake boundary condition specification. mCP, which is not directly calculated by the integration scheme and is equal to a product of two variables in the state vector (m and wCP), therefore represents a suitable criterion for testing the convergence of integration schemes. It can be seen from the figure that RK4/1 deg. CA and AM, FCN / TOLrei = 10-5 schemes provide convergent results, whereas both Euler schemes fail to provide convergent results. Fig. 4. Mass flux through the exhaust valve for a motored cylinder at an engine speed of 1500 rpm Fig. 4 shows mass flux through the exhaust valve for a motored cylinder (without combustion) at an engine speed of 1500 rpm. These operating conditions are characteristic for vehicle deceleration where the engine represents an energy sink. Modeling of such operating conditions is a demanding test for the stability of integration schemes. In order to point out the differences between the various integration schemes the mass flux through the exhaust valve is examined. Only Euler and RK4 schemes are examined. The integration schemes out of the CVODE package are not considered here, since adequate specifications of error tolerances and solver parameters prevent production of results exhibiting stability issues. The exhaust valve opens at 483 deg. CA resulting in an outflow from the cylinder. Afterwards, the cylinder approaches the bottom dead center (BDC), 540 deg. CA, and a flow reversal occurs, since no combustion occurred in the high pressure phase and some mass has flown out of the cylinder resulting in p > pE when the piston approaches the BDC. Immediately after the BDC the flow direction changes again. The pressure in the cylinder is similar to the pressure in the exhaust plenum in the vicinity of the BDC and the derivative of the cylinder volume is small. It is discernible from Fig. 4 that the first order Euler scheme produces unphysical oscillations in mass flux due to coupled mass flow and pressure fluctuations as analysed in Fig. 2. Decreasing the time step of the Euler scheme from 1 deg. CA to 0.25 deg. CA diminishes mass flow oscillations, however, it does not suppress them. Fig. 4 shows that the RK4 scheme does not produce oscillatory results. It can also be observed that the mass flow traces of the RK4/4 deg. CA and Euler/1 deg. CA schemes deviate significantly. 4.2 Test Case 1: Engine Speed 1500 rpm The trade-off between numerical accuracy and computational speed is first analysed for Test case 1 (Table 2). The integration increments of the Euler scheme were 0.1, 0.25, 0.5, and 1 deg. CA and those of the RK4 scheme were 0.25, 0.5, 1, 2, and 4 deg. CA. TOLrei values of the CVODE schemes were increased by a factor of 10 in the range from 10-7 to 10-5 for Newton iteration schemes and in the range from 10-7 to 10-3 for functional iteration schemes. Time increment of the CVODE schemes corresponds to 1 deg. CA. Additionally, for the BDF, FCN scheme results for the time increment corresponding to 4 deg. CA was analysed to investigate the influence of the global integration increment on the accuracy of results. In order to preserve readability of the figures only results for the BDF, NWT are shown. It might be noted that the results of the BDF, FCN scheme feature very similar trends as those of the BDF, NWT scheme. The TOLra of CVODE schemes thus ranges from 10-7 to 10-5 for BDF, NWT scheme, from 10-7 to 10-5 for AM, NWT scheme, and from 10-7 to 10-3 for AM, FCN scheme. In the subsequent figures, markers on a line indicating a particular integration scheme represent different integration increments or different tolerance settings. Larger integration increments or smaller tolerance settings feature lower RT factors for a particular scheme. Fig. 5 summarizes the dependency of the NDRMS on the RT factor for in-cylinder theromynamic variables (Figs. 5a to e). Additionally, Fig. 5f shows the deviation of the IMEP calculated by a particular scheme and the exact solution (AIMEP = IMEP(exact)-1) as a function of the RT factor, where IMEP is defined by Eq. (18) and represents an integral engine parameter. Efficient integration schemes simultaneously feature both a low NDRMS value and a low RT value and thus tend towards the lower left end of the diagrams as indicated in Fig. 5a. It can be observed that all analysed cylinder data feature very similar trends. It can be seen from Fig. 5, with some local deviations, that an increased RT value (increased computational time) leads to increased numerical accuracy for a particular scheme. However, it can also be observed that the difference in NDRMS values of different integration schemes could reach up to two orders of magnitude at the same RT factor. This fact underlines the necessity of selecting the appropriate numerical scheme for the specific simulation task. It can be seen that all CVODE schemes feature a pronounced increase in NDRMS values for decreasing RT values. Fig. 5 shows that CVODE schemes reach similar minimum RT values (i.e. upper computational speed), when applying tolerance and iteration parameter settings, which prevent the occurrence of non-stabilities. On the other hand, explicit schemes without error control do not include any mechanism for controlling numerical accuracy. These schemes thus allow one to attain very small RT values, however, numerical accuracy is reduced significantly. A thorough analysis of convergence and stability is therefore mandatory, especially for explicit methods, in order to avoid evaluation of non-converged simulation results. Based on the results shown in Figs. 2 and 3 and NDRMS values shown in Fig. 5, basic guidelines for determining NDRMS values that indicate converged solutions can be defined. It is shown in Figs. 2 and 3 that RK4/1 deg. CA and AM, FCN/TOLrel = 10-5 schemes provide results that coincide very well with the exact solution. These results could thus be assumed as converged considering numerical accuracy requirements in the initial stage of the engine development process or in the HiL tests. From Fig. 5 it can be concluded that both schemes, i.e. RK4/1 deg. CA and AM FCN/ TOLrel = 10-5, feature NDRMS values approximately 10-4 and lower. This value might thus be considered as a guideline for a converged solution taking into account numerical accuracy requirements needed in the engine development process. On the contrary, schemes featuring NDRMS values considerably larger than 10-4 exhibit larger Fig. 5. Dependencies of the NDRMS on the RT factor for an engine speed of 1500 rpm and an IMEP of 1.92 MPa; a) pressure, b) temperature, c) in-cylinder mass, d) mass fraction of combustion products, e) mass fraction of burned fuel, f) dependency of the AIMEP on the RT factor deviations from the exact solutions indicating non-converged solutions. Furthermore, Fig. 5f shows that the AIMEP is below 1% (a typical accuracy threshold for validating high fidelity ICE models) for all analysed integration schemes and all integration increments or specified tolerances. Considering the deficiencies of particular schemes shown in Figs. 2 to 4, it can be concluded that agreement of integral engine parameters should not be the only measure for assessing the accuracy of 231 the model, since non-converged and even locally nonstable solutions might comply with these accuracy requirements. Fig. 5 shows, with some local deviations, that RK4 schemes feature the lowest NDRMS values for the analysed RT range. It can further be concluded that the RK4/1 deg. CA scheme complies with the convergence and stability requirements and simultaneously features a very low RT value for this engine operating point. Fig. 5 also shows that the results of the RK4 scheme featuring a steep NDRMS increase in state variables for integration increments larger than 1 deg. CA. It is discernible from Fig. 5 that the Euler scheme features a worse trade-off between NDRMS and RT compared to the RK4 scheme for NDRMS values of 10-4 and lower, which is an acceptable level of numerical accuracy required in the engine development process, as discussed above. On the other hand, the Euler scheme performs better than the RK4 scheme for larger NDRMS values and thus small RT values. However it is discernible from Fig. 2 that these small RT values correspond to integration increments leading to numerical instabilities in both schemes. It can also be observed in Fig. 5 that the functional iteration applied by the AM scheme results in better trade-off between NDRMS and RT compared to the application of the Newton iteration. This can be explained by the fact that the computationally fast functional iteration brings benefits over the computationally more demanding Newton iteration for this non-stiff example. Additionally, it can be observed that the BDF scheme might even perform better than the AM, NWT scheme although it is also applicable to stiff systems. Unlike for the AM scheme, only minor differences can be observed in the trade- off between NDRMS and RT when applying different iteration methods to the BDF scheme. It is also shown for the BDF scheme that increasing the global integration increment from 1 to 4 deg. CA does not influence the trade-off between NDRMS and RT. This is in line with the above statement that the definition of tolerances represents the main mechanism influencing the numerical accuracy and computational speed of CVODE schemes. Due to these facts only the results of the BDF, NWT scheme are shown in the figures. The analysis of different engine configurations has shown that the ratio of the effective valve flow area (Mgeom) and cylinder volume (V) most significantly influences the trade-off between NDRMS and RT values of a particular engine. This is justified by the fact that the gas exchange period generally represents the numerically most demanding phase, provided simple combustion models are used (e.g. Vibe, ROHR table, and the majority of empirical models that do not incorporate chemical kinetics or droplet evaporation). This ratio together with the compression ratio strongly influences the ratio of the mass exchanged through the valves within one integration increment and the in-cylinder mass, and therefore the stability of the solution during the gas exchange period. The ratio of the mass exchanged through the valves within one integration increment and the in-cylinder mass thus most significantly determines the optimum integration increment of the explicit schemes when considering the trade-off between NDRMS and RT values. Similarly, the solvers of the CVODE schemes adapt their integration increments to consider this effect. Therefore, the ratio of RT factors of the explicit and CVODE schemes is almost unaffected by the variation of the ratio of ^Ageom to V for a selected NDRMS value. However, the RT factors of all schemes increase a) rt [-] b) RT [-] Fig. 6. Dependencies of the NDRMS on the RT factor for an engine speed of 6000 rpm and IMEP = 1.68 MPa; a) pressure, b) temperature at the same value of NDRMS if the ratio of jiAgeom to V increases. 4.3 Test Case 2: Engine Speed 6000 rpm The computational time of CA based explicit schemes without error control directly correlates with the engine speed. In contrast, time based methods with error control and adaptive integration increments might not strictly follow this rule. Fig. 6 shows results for the Test case 2 (Table 3). nE = 6000 rpm corresponds to the upper speed range of the internal mixture preparation compression ignition engines applied in racing cars. For such operating conditions it is very common for turbocharged engines to feature exhaust pressure that is higher compared to the boost pressure as shown in Table 3. Fig. 6 shows the results of NDRMS vs. RT for the cylinder pressure (Fig. 6a) and temperature (Fig. 6b), since it was shown in Fig. 5 that other state vector variables feature similar trends. It is discernible from Fig. 6 that the NDRMS vs. RT curves of the CVODE schemes move closer to the corresponding RK4 curve compared to the nE= 1500 rpm case shown in Fig. 5. This is mainly a consequence of different integration variables, i.e. time or CA rotation, as pointed out above. Although dV/dt (Eq. (10)) increases with increasing engine speed, the time needed to integrate the coupled system of equations (Eq. (14)) with time based CVODE schemes does not increase proportionally to the engine speed. Similarly to Fig. 5, it can be observed, with some local deviations, that RK4 schemes feature the lowest NDRMS values for the analysed RT range. It can be seen that the RK4/1 deg. CA scheme complies with convergence and stability requirements. The calculated RT factor of the RK4/1 deg. CA scheme is 0.19 for the applied platform. It could thus comply with real-time constraints for a very simple multi-cylinder engine configurations at nE = 6000 rpm. Fig. 6 indicates that for very low values of NDRMS, CVODE schemes tend to show a better performance than RK4 schemes. However, at this engine speed the RT factor of CVODE schemes with small values of TOLrel approaches or already exceeds 1 for this simplified test case. It can also be concluded that it is not reasonable to apply integration schemes with very low tolerances, since other errors introduced by the application of the 0-D model are much larger and therefore eliminate the benefit of very low tolerances of the integration scheme. It can be observed in Fig. 6 that the difference in the trade-off between NDRMS and RT value evaluated for the AM, FCN, and AM, NWT scheme is reduced compared to the lower engine speed (Fig. 5). This can be explained by the fact that for higher engine speeds time derivatives of the variables increase and thus the computational overload of the more sophisticated gradient based Newton method is less pronounced than that of the functional iteration method. It might be concluded that BDF, NWT scheme performs comparable to both AM schemes for this engine speed. Based on the above findings it could be assumed that time based CVODE schemes might show a better performance than RK4 at much higher engine speeds. Therefore an example with nE = 12000 rpm was analysed, although this engine speed represents a rather unrealistically high value for internal mixture preparation compression ignition engines. Again the RK4/1 deg. CA scheme featured the lowest RT values (RT = 0.38) for NDRMS values in the range of 10-4. Table 3. Boundary conditions for Test case 2: engine speed 6000 rpm and IMEP = 1.68 MPa Plenum T [K] p [Pa] Wfb [-] Wcp [-] Wfv [-] Intake 340 200000 0 0 0 Exhaust 320 280000 0 0 0 4.4 Comparison of Numerical and Experimental Results In addition to the thorough theoretical analyses presented above it is also meaningful to compare numerical and experimental in-cylinder pressure traces for a series production engine to confirm the applicability and accuracy of the cylinder model. Therefore, the crank angle resolved cylinder model presented here is embedded into a mean value based gas path model of the engine manifolds [4]. As presented in [4], the overall model is applied to a turbocharged, intercooled 4-cylinder diesel engine with a cooled external EGR and variable geometry turbine. The combustion is modelled using rate-of-heat-release curves derived from experimental in-cylinder pressure traces. This ensures a consistent comparison of numerical and experimental cylinder pressure traces [4]. Fig. 7 shows numerical and experimental cylinder pressure traces for two selected full and part load operating conditions applying a RK4/1 deg. CA scheme to integrate the cylinder governing equations. It is clear that the model is capable of accurately predicting high- and low-pressure phase of the cylinder process. A good agreement can be seen for the early combustion characteristic for the full load condition and also for the late combustion characteristic for part load If ■ 14 ■ 12-■ran) ■ - : all r\ i * : ! \ : j \ : ; : f : ß - ; bi) ; / i g i \ : t t v : / \ I.'' i i i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 500 400 ■ 300 ■ O.200 ■ 100 - 300 - ra 200 ■ 10D ■ } i t i: t t 1 i V J " i I 1 j'i / : -------experiment J - i l ; b2) r i \ \ i i i iL.-«—' i .J i i i i i i i i i 1 1 1 -eo -3d 30 60 90 120 -300 -1E0 1.BD 3B0 CA [d eg] CA [deg] Fig. 7. Comparison of numerical and experimental cylinder pressure traces for different engine load points; a) 3800 rpm and 1.54 MPa BMEP, b) 1600 rpm and 0.33 MPa BMEP conditions. Minor discrepancies in pressure traces can be observed for the exhaust stroke. The deviations are caused by the mean value gas path description of the engine manifolds neglecting gas-dynamic effects. The considered engine features nearly negligible valve overlap and a high compression ratio. Thus the discrepancies during the exhaust stroke do not notably influence the in-cylinder state at EVC. Therefore, due to the application of consistent governing equations and adequate methods for their integration, numerical and experimental results of the intake stroke, compression, combustion, and expansion show very good agreement. Good agreement between the simulation and experimental results is furthermore confirmed by a very good agreement in simulated and experimental values of the peak cylinder pressure for 3 full load and 1 part load point as given in Table 4. Table 4. Peak firing pressure In MPa Engine speed BMEP 1200 rpm 1.05 MPa 3000 rpm 1.94 MPa 3800 rpm 1.54 MPa 1600 rpm 0.33 MPa Simulation 10.29 16.33 16.37 5.25 Experiment 10.24 16.19 16.24 5.29 5 CONCLUSIONS The trade-off between numerical accuracy and computational speed was analysed for different integration schemes applied to solve governing equations of the 0-D internal combustion engine cylinder model. It was shown that the numerical accuracy of integration schemes differs significantly at equal computational effort. The analysis reveals that the RK4 scheme features the best ratio of numerical accuracy and computational speed in a realtime capable computation speed range. It was shown that explicit schemes without error control should apply CA rotation as an independent variable, since this approach ensures similar numerical accuracy for different engine operating conditions. It was also demonstrated that for a particular numerical accuracy threshold the difference in computational speed between explicit schemes applying CA rotation as an independent variable and time based methods with error control decreases with increasing engine speed. Moreover, CVODE schemes show a comparable or better performance than the RK4 scheme for high engine speeds when very high accuracy of the results is required. However, it is generally not reasonable to apply integration schemes with very low tolerances, since other errors introduced by the application of the 0-D model are much larger and therefore eliminate the benefit of very low tolerances of the integration scheme. It was shown that the global integration increment of CVODE schemes does not significantly influence the performance of these schemes, being predominantly influenced by tolerance settings. NDRMS values give valuable guidelines on the numerical accuracy of the integration schemes. However, they feature a limited capability of detecting local deficiencies produced by integration schemes. The NDRMS vs. RT approach is therefore suitable for acquiring a basic insight into the performance of integration schemes, whereas cycle resolved results need to be subsequently analysed to study the convergence and stability of the solution in detail. It was also shown that agreement of integral engine parameters should not be the only measure for assessing the accuracy of the model. Very good agreement between numerical and experimental results additionally confirms the accuracy of the cylinder model and its applicability for predictive system level simulations. 6 NOMENCLATURE A area [m2] H enthalpy [J] h specific enthalpy [m2/s2] M molar mass [kg/kmol] m mass [kg] N number of species in the species vector [- ^valves number of valves [-] nE engine speed [rpm] P power [W] P pressure [Pa] Q heat [J] R specific gas constant [J/kgK] m universal gas const. [J/kmolK] T temperature [T] t time [s] TOL tolerance [-] U internal energy [J] u specific internal energy [m2/s2] V volume [m3] wk species concentration [-] y solution vector [-] Z compressibility factor [-] 9 crank angle rotation [deg. CA] K ratio of specific heats [-] M flow coefficient [-] t independent variable: t or