MECCA 02 2017 PAGE 37 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK 10.1515/mecdc ‑2017 ‑0007 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK RASTISLAV TOMAN, JAN MACEK CTU in Prague, Faculty of Mechanical Engineering; Technická 4, Praha 6, 166 07, Czech Republic E ‑mail: rastislav.toman@fs.cvut.cz, jan.macek@fs.cvut.cz ABSTRACT The current study evaluates the predictive capabilities of a new phenomenological combustion model, available as a part of the GT ‑Suite software package. It is comprised of two main sub ‑models: 0D model of in ‑cylinder flow and turbulence, and turbulent SI combustion model. The 0D in ‑cylinder flow model (EngCylFlow) uses a combined K ‑k ‑ε kinetic energy cascade approach to predict the evolution of the in ‑cylinder charge motion and turbulence, where K and k are the mean and turbulent kinetic energies, and ε is the turbulent dissipation rate. The subsequent turbulent combustion model (EngCylCombSITurb) gives the in ‑cylinder burn rate; based on the calculation of flame speeds and flame kernel development. This phenomenological approach reduces significantly the overall computational effort compared to the 3D ‑CFD, thus allowing the computation of full engine operating map and the vehicle driving cycles. Model was calibrated using a full map measurement from a turbocharged natural gas SI engine, with swirl intake ports. Sensitivity studies on different calibration methods, and laminar flame speed sub ‑models were conducted. Validation process for both the calibration and sensitivity studies was concerning the in ‑cylinder pressure traces and burn rates for several engine operation points achieving good overall results. KEY WORDS: PREDICTIVE PHENOMENOLOGICAL MODEL; INTERNAL COMBUSTION ENGINE; SPARK ‑IGNITION; K ‑k ‑ε KINETIC ENERGY CASCADE; 0D IN ‑CYLINDER FLOW MODEL; TURBULENT SI COMBUSTION MODEL; NATURAL GAS ENGINE; GENETIC ALGORITHM; GT ‑SUITE SHRNUTÍ Predkladaný článok hodnotí prediktívne schopnosti nového fenomenologického modelu horenia, ktorý je k dispozícii ako súčasť softvérového balíka GT ‑Suite. Skladá sa z dvoch hlavných sub ‑modelov: 0D modelu prúdenia a turbulencie vo valci a zážihového turbulentného modelu horenia. 0D model prúdenia a turbulencie vo valci (EngCylFlow) používa kombinovaný prístup K ‑k ‑ε kaskády kinetickej energie na predpoveď pohybu zmesi a turbulencie vo valci, kde K a k sú stredné a turbulentné kinetické energie a  ε je turbulentná rýchlosť disipácie. Následný model turbulentného horenia (EngCylCombSITurb) určuje rýchlosť horenia vo valci na základe výpočtu rýchlosti čela plameňa a vývoja jadra plameňa. Tento fenomenologický prístup výrazne znižuje celkovú výpočtovú náročnosť v porovnaní s 3D ‑CFD, čo umožňuje výpočet úplnej charakteristiky spaľovacieho motora a jazdných cyklov vozidla. Model bol kalibrovaný pomocou meraní úplnej charakteristiky preplňovaného zážihového motora na zemný plyn so swirlovými vstupnými kanálmi. Boli vykonané štúdie citlivosti na rôzne kalibračné metódy a na rôzne sub ‑modely laminárnej rýchlosti čela plameňa. Validačný proces pre kalibrácie a štúdie citlivosti sa týkal tlaku vo valci a rýchlostí horenia pre niekoľko pracovných bodov motora, dosahujúc dobré celkové výsledky. KLÍČOVÁ SLOVA: PREDIK TÍVNY FENOMENOLOGICKÝ MODEL; ZÁŽIHOVÝ SPAĽOVACÍ MOTOR; K ‑k ‑ε KASKÁDA KINETICKEJ ENERGIE; 0D MODELU PRÚDENIA VO VALCI; TURBULENTNÝ ZÁŽIHOVÝ MODEL HORENIA; MOTOR NA ZEMNÝ PLYN; GENETICKÝ ALGORITMUS; GT ‑SUITE EVALUATION OF THE PREDICTIVE CAPABILITIES OF A PHENOMENOLOGICAL COMBUSTION MODEL FOR NATURAL GAS SI ENGINE MECCA 02 2017 PAGE 38 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK 1. INTRODUCTION Current development of the internal combustion engines (ICE) is focused on the overall efficiency improvement and emissions reduction. To fulfil these goals, downsizing of the ICE presents one of the most valuable options. But the increasing boost levels also lead to an increase in the knock likelihood, requiring spark timing retardation or mixture enrichment. Moreover, current engines use progressively also additional advanced control systems, such as advanced gas exchange systems, cylinder deactivation, or variable compression ‑ratio systems. With many iterations needed, accurate and robust modelling of the combustion process have become essential during the ICE development process, with an emphasis on the overall simulation time of one engine operating cycle. A detailed 3D Computational Fluid Dynamics (3D ‑CFD) analysis of the in ‑cylinder flows, charge motion and combustion leads to accurate prediction of burn rate (if set‑ ‑up properly), but with the obvious drawback of its high computational demands [1]. 3D ‑CFD is therefore used mostly for the analysis of separate engine operating points. Empirical combustion models usually use an approximation of a measured burn rate. The most common empirical model is a Vibe formula [2]. However, if the user wants to obtain correct burn rate values in changed ICE operating conditions, reference burn rate pattern must be adjusted by additional formulas [3], [4]. In general, empirical models are simple and work well inside the calibrated region, but their extrapolation abilities are poor [5]. Multi ‑zone models of combustion and heat transfer in SI engines present a fast, accurate, stable and above all physical‑ ‑based solution. A general theory of zone models based on the laws on conservation is described in [6], with a comparison of Lagrangian and Eulerian approaches. Recent paper of Hvezda [7] presents a specific adaptive approach to the chemical transformation. Multi ‑zone model of Hvezda models in detail the flame velocities using a turbulent coefficient, and accounts for the real geometry of the combustion chamber. Finally, phenomenological combustion models also respect the combustion chamber geometry and obtain a burn rate by the calculation of turbulent flame speed and instantaneous flame area [8]. These models need an information on in‑ ‑cylinder flow quantities as well. Several 0D turbulence models aim to reproduce the complex 3D phenomena, mainly by k ‑ε approach [9], [10], [11] or K ‑k approach recently studied in [12] and [13]. Both, the multi ‑zone and phenomenological models show very good extrapolation capabilities and low computational demands, allowing for the fast simulation of vehicle driving cycles. A combustion model evaluated in this study consists of two main sub ‑models: 0D in ‑cylinder flow model EngCylFlow (or Flow) and turbulent combustion model EngCylCombSITurb (or SITurb). Gamma Technologies is currently developing both sub ‑models as a part of GT ‑Suite software package [14]. The Flow model combines the k ‑ε approach of Morel et al. [9], [10] with the K ‑k into combined K ‑k ‑ε kinetic energy cascade model. Fogla et al. [15] describes the current Flow model, comparing with 3D ‑CFD and former model version, using two similar turbocharged gasoline ICEs, with tumble intake ports. The SITurb model – originally developed by Wahiduzzaman, Morel and Sheard [8] – uses a turbulent flame concept, directly linked to the in ‑cylinder flow and turbulence calculation. Mirzaeian et al. [16] assessed the predictive capability of the current model version adding the equation system description. They also proposed a calibration method starting with DoE calibration of the Flow model against the 3D ‑CFD data and continuing with the SITurb calibration using the Genetic Algorithms (GA) with the objective to match the burn rate against the measurement data (turbocharged gasoline ICE with tumble intake port). More detail on the combustion model follows in Section 2. 1.1 MAIN GOALS The main objectives of this paper are: • First, to calibrate the current combustion model, obtaining a single set of optimal model parameters; • Second, to test its predictive capabilities. Already mentioned papers [15] and [16] evaluated the model capabilities on turbocharged gasoline engines, with tumble intake ports. This study uses a full engine map measurement set with EGR variations and stoichiometric conditions of the ICE fueled by natural gas, with swirl intake ports. The additional goals of the paper are following: • to summarize the main features of the predictive combustion model; • to compare different calibration approaches; • to test the sensitivities of the combustion model on the laminar flame speed. 2. PREDICTIVE COMBUSTION MODEL 2.1 IN ‑CYLINDER FLOW MODEL The main equation system of the Flow model contains three differential equations (equations 1 ‑3) that govern the mean kinetic energy K=(1⁄2)U 2 (U is the mean velocity inside the cylinder), turbulent kinetic energy k=(3⁄2)u' 2 (u' is the mean fluctuating turbulent velocity inside the cylinder), and the MECCA 02 2017 PAGE 39 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK turbulent dissipation rate ε. The model assumes homogeneous and isotropic turbulent field [15].  to summarize the main features of the predictive combustion model;  to compare different calibration approaches;  to test the sensitivities of the combustion model on the laminar flame speed. 2. Predictive combustion model 2.1 In-cylinder flow model The main equation system of the Flow model contains three differential equations (equations 1-3) that govern the mean kinetic energy � � �1 2⁄ ��� (� is the mean velocity inside the cylinder), turbulent kinetic energy � � �3 2⁄ ���� (�� is the mean fluctuating turbulent velocity inside the cylinder), and the turbulent dissipation rate ϵ. The model assumes homogeneous and isotropic turbulent field [15]. ����� �� � ����1 � ������� � ��� ��� � �� (1) ����� �� � ��������� � ��� ��� � �� � ������ � �� (2) ����� �� � ������ √� �� � ��� ��� � �� � ������ √� �� � 1.�2 ��� � (3) First right-hand side terms in all three equations represent the production of each flow quantity, with the inflow energy ���. Parameter ��� indicates the fraction of inflow energy entering the cylinder directly as turbulence, although not generated by the kinetic energy cascade process. The second right-hand side term in the main equations, describes the energy out-flow through the valves, with the mass flow rate of the cylinder exit flow �� ���. Production terms �� and �� (equations 4-5) model the production of the turbulent kinetic energy and a dissipation rate from the large scale mean flows via the kinetic energy cascade process; �� represents a turbulent viscosity, � a density and �� a rate of change of density inside the cylinder. The appendix of [15] describes the evolution of these terms. �� � ���� 2�� ��� � 23�� � �� �� � 2 3��� � �� �� � (4) �� � � � ��.�6���� �� ��� � 2�� ����� � 2.64 3 ��� � �� �� � � (5) The last right-hand side term in each of the three main equations is a sink term for its respective quantity. The terms with the quantity � model the production of turbulence by the decay of the tumble macro-vortex during the compression [15]. Simple equation systems for the time rate change of the angular momentum �� ��⁄ model the rotational components of the flow – tumble and swirl – as a single macro-vortex undergoing stretching and compression during the engine intake and compression. Swirl and tumble are produced by the incoming charge, accounting for the measured swirl and tumble coefficients, and reduced by the cylinder out-flow. Equation systems for both rotational motions contain a proper decay functions. Paper [15] does not discuss the swirl decay function, but provides further information on the tumble decay function. Former Flow model ([10], [11]) accounted for the squish motion (inside the swirl model) and injection event kinetic energy and so does the current Flow model. However, the exact equation systems are not available in [15]. Since the current Flow model calculates the kinetic energy and the dissipation rate, the evolution of integral length scale �� over time is then obtained directly with the equation 6 (�� � �.�� is a standard k-ε model constant) [15]. �� � ��� �⁄ �� �⁄ � (6) (1)  to summarize the main features of the predictive combustion model;  to compare different calibration approaches;  to test the sensitivities of the combustion model on the laminar flame speed. 2. Predictive combustion model 2.1 In-cylinder flow model The main equation system of the Flow model contains three differential equations (equations 1-3) that govern the mean kinetic energy � � �1 2⁄ ��� (� is the mean velocity inside the cylinder), turbulent kinetic energy � � �3 2⁄ ���� (�� is the mean fluctuating turbulent velocity inside the cylinder), and the turbulent dissipation rate ϵ. The model assumes homogeneous and isotropic turbulent field [15]. ����� �� � ����1 � ������� � ��� ��� � �� (1) ����� �� � ��������� � ��� ��� � �� � ������ � �� (2) ����� �� � ������ √� �� � ��� ��� � �� � ������ √� �� � 1.�2 ��� � (3) First right-hand side terms in all three equations represent the production of each flow quantity, with the inflow energy ���. Parameter ��� indicates the fraction of inflow energy entering the cylinder directly as turbulence, although not generated by the kinetic energy cascade process. The second right-hand side term in the main equations, describes the energy out-flow through the valves, with the mass flow rate of the cylinder exit flow �� ���. Production terms �� and �� (equations 4-5) model the production of the turbulent kinetic energy and a dissipation rate from the large scale mean flows via the kinetic energy cascade process; �� represents a turbulent viscosity, � a density and �� a rate of change of density inside the cylinder. The appendix of [15] describes the evolution of these terms. �� � ���� 2�� ��� � 23�� � �� �� � 2 3��� � �� �� � (4) �� � � � ��.�6���� �� ��� � 2�� ����� � 2.64 3 ��� � �� �� � � (5) The last right-hand side term in each of the three main equations is a sink term for its respective quantity. The terms with the quantity � model the production of turbulence by the decay of the tumble macro-vortex during the compression [15]. Simple equation systems for the time rate change of the angular momentum �� ��⁄ model the rotational components of the flow – tumble and swirl – as a single macro-vortex undergoing stretching and compression during the engine intake and compression. Swirl and tumble are produced by the incoming charge, accounting for the measured swirl and tumble coefficients, and reduced by the cylinder out-flow. Equation systems for both rotational motions contain a proper decay functions. Paper [15] does not discuss the swirl decay function, but provides further information on the tumble decay function. Former Flow model ([10], [11]) accounted for the squish motion (inside the swirl model) and injection event kinetic energy and so does the current Flow model. However, the exact equation systems are not available in [15]. Since the current Flow model calculates the kinetic energy and the dissipation rate, the evolution of integral length scale �� over time is then obtained directly with the equation 6 (�� � �.�� is a standard k-ε model constant) [15]. �� � ��� �⁄ �� �⁄ � (6) (2)  to summarize the main features of the predictive combustion model;  to compare different calibration approaches;  to test the sensitivities of the combustion model on the laminar flame speed. 2. Predictive combustion model 2.1 In-cylinder flow model The main equation system of the Flow model contains three differential equations (equations 1-3) that govern the mean kinetic energy � � �1 2⁄ ��� (� is the mean velocity inside the cylinder), turbulent kinetic energy � � �3 2⁄ ���� (�� is the mean fluctuating turbulent velocity inside the cylinder), and the turbulent dissipation rate ϵ. The model assumes homogeneous and isotropic turbulent field [15]. ����� �� � ����1 � ������� � ��� ��� � �� (1) ����� �� � ��������� � ��� ��� � �� � ������ � �� (2) ����� �� � ������ √� �� � ��� ��� � �� � ������ √� �� � 1.�2 ��� � (3) First right-hand side terms in all three equations represent the production of each flow quantity, with the inflow energy ���. Parameter ��� indicates the fraction of inflow energy entering the cylinder directly as turbulence, although not generated by the kinetic energy cascade process. The second right-hand side term in the main equations, describes the energy out-flow through the valves, with the mass flow rate of the cylinder exit flow �� ���. Production terms �� and �� (equations 4-5) model the production of the turbulent kinetic energy and a dissipation rate from the large scale mean flows via the kinetic energy cascade process; �� represents a turbulent viscosity, � a density and �� a rate of change of density inside the cylinder. The appendix of [15] describes the evolution of these terms. �� � ���� 2�� ��� � 23�� � �� �� � 2 3��� � �� �� � (4) �� � � � ��.�6���� �� ��� � 2�� ����� � 2.64 3 ��� � �� �� � � (5) The last right-hand side term in each of the three main equations is a sink term for its respective quantity. The terms with the quantity � model the production of turbulence by the decay of the tumble macro-vortex during the compression [15]. Simple equation systems for the time rate change of the angular momentum �� ��⁄ model the rotational components of the flow – tumble and swirl – as a single macro-vortex undergoing stretching and compression during the engine intake and compression. Swirl and tumble are produced by the incoming charge, accounting for the measured swirl and tumble coefficients, and reduced by the cylinder out-flow. Equation systems for both rotational motions contain a proper decay functions. Paper [15] does not discuss the swirl decay function, but provides further information on the tumble decay function. Former Flow model ([10], [11]) accounted for the squish motion (inside the swirl model) and injection event kinetic energy and so does the current Flow model. However, the exact equation systems are not available in [15]. Since the current Flow model calculates the kinetic energy and the dissipation rate, the evolution of integral length scale �� over time is then obtained directly with the equation 6 (�� � �.�� is a standard k-ε model constant) [15]. �� � ��� �⁄ �� �⁄ � (6)  to summarize the main features of the predictive combustion model;  to compare different calibration approaches;  to test the sensitivities of the combustion model on the laminar flame speed. 2. Predictive combustion model 2.1 In-cylinder flow model The main equation system of the Flow model contains three differential equations (equations 1-3) that govern the mean kinetic energy � � �1 2⁄ ��� (� is the mean velocity inside the cylinder), turbulent kinetic energy � � �3 2⁄ ���� (�� is the mean fluctuating turbulent velocity inside the cylinder), and the turbulent dissipation rate ϵ. The model assumes homogeneous and isotropic turbulent field [15]. ����� �� � ����1 � ������� � ��� ��� � �� (1) ����� �� � ��������� � ��� ��� � �� � ������ � �� (2) ����� �� � ������ √� �� � ��� ��� � �� � ������ √� �� � 1.�2 ��� � (3) First right-hand side terms in all three equations represent the production of each flow quantity, with the inflow energy ���. Parameter ��� indicates the fraction of inflow energy entering the cylinder directly as turbulence, although not generated by the kinetic energy cascade process. The second right-hand side term in the main equations, describes the energy out-flow through the valves, with the mass flow rate of the cylinder exit flow �� ���. Production terms �� and �� (equations 4-5) model the production of the turbulent kinetic energy and a dissipation rate from the large scale mean flows via the kinetic energy cascade process; �� represents a turbulent viscosity, � a density and �� a rate of change of density inside the cylinder. The appendix of [15] describes the evolution of these terms. �� � ���� 2�� ��� � 23�� � �� �� � 2 3��� � �� �� � (4) �� � � � ��.�6���� �� ��� � 2�� ����� � 2.64 3 ��� � �� �� � � (5) The last right-hand side term in each of the three main equations is a sink term for its respective quantity. The terms with the quantity � model the production of turbulence by the decay of the tumble macro-vortex during the compression [15]. Simple equation systems for the time rate change of the angular momentum �� ��⁄ model the rotational components of the flow – tumble and swirl – as a single macro-vortex undergoing stretching and compression during the engine intake and compression. Swirl and tumble are produced by the incoming charge, accounting for the measured swirl and tumble coefficients, and reduced by the cylinder out-flow. Equation systems for both rotational motions contain a proper decay functions. Paper [15] does not discuss the swirl decay function, but provides further information on the tumble decay function. Former Flow model ([10], [11]) accounted for the squish motion (inside the swirl model) and injection event kinetic energy and so does the current Flow model. However, the exact equation systems are not available in [15]. Since the current Flow model calculates the kinetic energy and the dissipation rate, the evolution of integral length scale �� over time is then obtained directly with the equation 6 (�� � �.�� is a standard k-ε model constant) [15]. �� � ��� �⁄ �� �⁄ � (6) (3) First right ‑hand side terms in all three equations represent the production of each flow quantity, with the inflow energy Ein. Parameter αin indicates the fraction of inflow energy entering the cylinder directly as turbulence, although not generated by the kinetic energy cascade process. The second right ‑hand side term in the main equations, describes the energy out ‑flow through the valves, with the mass flow rate of the cylinder exit flow m· out. Production terms Pk and PЄ (equations 4 ‑5) model the production of the turbulent kinetic energy and a dissipation rate from the large scale mean flows via the kinetic energy cascade process; υT represents a turbulent viscosity, ρ a density and ρ· a rate of change of density inside the cylinder. The appendix of [15] describes the evolution of these terms.  to summarize the main features of the predictive combustion model;  to compare different calibration approaches;  to test the sensitivities of the combustion model on the laminar flame speed. 2. Predictive combustion model 2.1 In-cylinder flow model The main equation system of the Flow model contains three differential equations (equations 1-3) that govern the mean kinetic energy � � �1 2⁄ ��� (� is the mean velocity inside the cylinder), turbulent kinetic energy � � �3 2⁄ ���� (�� is the mean fluctuating turbulent velocity inside the cylinder), and the turbulent dissipation rate ϵ. The model assumes homogeneous and isotropic turbulent field [15]. ����� �� � ����1 � ������� � ��� ��� � �� (1) ����� �� � ��������� � ��� ��� � �� � ������ � �� (2) ����� �� � ������ √� �� � ��� ��� � �� � ������ √� �� � 1.�2 ��� � (3) First right-hand side terms in all three equations represent the production of each flow quantity, with the inflow energy ���. Parameter ��� indicates the fraction of inflow energy entering the cylinder directly as turbulence, although not generated by the kinetic energy cascade process. The second right-hand side term in the main equations, describes the energy out-flow through the valves, with the mass flow rate of the cylinder exit flow �� ���. Production terms �� and �� (equations 4-5) model the production of the turbulent kinetic energy and a dissipation rate from the large scale mean flows via the kinetic energy cascade process; �� represents a turbulent viscosity, � a density and �� a rate of change of density inside the cylinder. The appendix of [15] describes the evolution of these terms. �� � ���� 2�� ��� � 23�� � �� �� � 2 3��� � �� �� � (4) �� � � � ��.�6���� �� ��� � 2�� ����� � 2.64 3 ��� � �� �� � � (5) The last right-hand side term in each of the three main equations is a sink term for its respective quantity. The terms with the quantity � model the production of turbulence by the decay of the tumble macro-vortex during the compression [15]. Simple equation systems for the time rate change of the angular momentum �� ��⁄ model the rotational components of the flow – tumble and swirl – as a single macro-vortex undergoing stretching and compression during the engine intake and compression. Swirl and tumble are produced by the incoming charge, accounting for the measured swirl and tumble coefficients, and reduced by the cylinder out-flow. Equation systems for both rotational motions contain a proper decay functions. Paper [15] does not discuss the swirl decay function, but provides further information on the tumble decay function. Former Flow model ([10], [11]) accounted for the squish motion (inside the swirl model) and injection event kinetic energy and so does the current Flow model. However, the exact equation systems are not available in [15]. Since the current Flow model calculates the kinetic energy and the dissipation rate, the evolution of integral length scale �� over time is then obtained directly with the equation 6 (�� � �.�� is a standard k-ε model constant) [15]. �� � ��� �⁄ �� �⁄ � (6) (4)  to summarize the main features of the predictive combustion model;  to compare different calibration approaches;  to test the sensitivities of the combustion model on the laminar flame speed. 2. Predictive combustion model 2.1 In-cylinder flow model The main equation system of the Flow model contains three differential equations (equations 1-3) that govern the mean kinetic energy � � �1 2⁄ ��� (� is the mean velocity inside the cylinder), turbulent kinetic energy � � �3 2⁄ ���� (�� is the mean fluctuating turbulent velocity inside the cylinder), and the turbulent dissipation rate ϵ. The model assumes homogeneous and isotropic turbulent field [15]. ����� �� � ����1 � ������� � ��� ��� � �� (1) ����� �� � ��������� � ��� ��� � �� � ������ � �� (2) ����� �� � ������ √� �� � ��� ��� � �� � ������ √� �� � 1.�2 ��� � (3) First right-hand side terms in all three equations represent the production of each flow quantity, with the inflow energy ���. Parameter ��� indicates the fraction of inflow energy entering the cylinder directly as turbulence, although not generated by the kinetic energy cascade process. The second right-hand side term in the main equations, describes the energy out-flow through the valves, with the mass flow rate of the cylinder exit flow �� ���. Production terms �� and �� (equations 4-5) model the production of the turbulent kinetic energy and a dissipation rate from the large scale mean flows via the kinetic energy cascade process; �� represents a turbulent viscosity, � a density and �� a rate of change of density inside the cylinder. The appendix of [15] describes the evolution of these terms. �� � ���� 2�� ��� � 23�� � �� �� � 2 3��� � �� �� � (4) �� � � � ��.�6���� �� ��� � 2�� ����� � 2.64 3 ��� � �� �� � � (5) The last right-hand side term in each of the three main equations is a sink term for its respective quantity. The terms with the quantity � model the production of turbulence by the decay of the tumble macro-vortex during the compression [15]. Simple equation systems for the time rate change of the angular momentum �� ��⁄ model the rotational components of the flow – tumble and swirl – as a single macro-vortex undergoing stretching and compression during the engine intake and compression. Swirl and tumble are produced by the incoming charge, accounting for the measured swirl and tumble coefficients, and reduced by the cylinder out-flow. Equation systems for both rotational motions contain a proper decay functions. Paper [15] does not discuss the swirl decay function, but provides further information on the tumble decay function. Former Flow model ([10], [11]) accounted for the squish motion (inside the swirl model) and injection event kinetic energy and so does the current Flow model. However, the exact equation systems are not available in [15]. Since the current Flow model calculates the kinetic energy and the dissipation rate, the evolution of integral length scale �� over time is then obtained directly with the equation 6 (�� � �.�� is a standard k-ε model constant) [15]. �� � ��� �⁄ �� �⁄ � (6) (5) The last right ‑hand side term in each of the three main equations is a sink term for its respective quantity. The terms with the quantity T model the production of turbulence by the decay of the tumble macro ‑vortex during the compression [15]. Simple equation systems for the time rate change of the angular momentum dL/dt model the rotational components of the flow – tumble and swirl – as a single macro ‑vortex undergoing stretching and compression during the engine intake and compression. Swirl and tumble are produced by the incoming charge, accounting for the measured swirl and tumble coefficients, and reduced by the cylinder out ‑flow. Equation systems for both rotational motions contain a proper decay functions. Paper [15] does not discuss the swirl decay function, but provides further information on the tumble decay function. Former Flow model ([10], [11]) accounted for the squish motion (inside the swirl model) and injection event kinetic energy and so does the current Flow model. However, the exact equation systems are not available in [15]. Since the current Flow model calculates the kinetic energy and the dissipation rate, the evolution of integral length scale Lt over time is then obtained directly with the equation 6 (Cμ = 0.09 is a standard k ‑ε model constant) [15].  to summarize the main features of the predictive combustion model;  to compare different calibration approaches;  to test the sensitivities of the combustion model on the laminar flame speed. 2. Predictive combustion model 2.1 In-cylinder flow model The main equation system of the Flow model contains three differential equations (equations 1-3) that govern the mean kinetic energy � � �1 2⁄ ��� (� is the mean velocity inside the cylinder), turbulent kinetic energy � � �3 2⁄ ���� (�� is the mean fluctuating turbulent velocity inside the cylinder), and the turbulent dissipation rate ϵ. The model assumes homogeneous and isotropic turbulent field [15]. ����� �� � ����1 � ������� � ��� ��� � �� (1) ����� �� � ��������� � ��� ��� � �� � ������ � �� (2) ����� �� � ������ √� �� � ��� ��� � �� � ������ √� �� � 1.�2 ��� � (3) First right-hand side terms in all three equations represent the production of each flow quantity, with the inflow energy ���. Parameter ��� indicates the fraction of inflow energy entering the cylinder directly as turbulence, although not generated by the kinetic energy cascade process. The second right-hand side term in the main equations, describes the energy out-flow through the valves, with the mass flow rate of the cylinder exit flow �� ���. Production terms �� and �� (equations 4-5) model the production of the turbulent kinetic energy and a dissipation rate from the large scale mean flows via the kinetic energy cascade process; �� represents a turbulent viscosity, � a density and �� a rate of change of density inside the cylinder. The appendix of [15] describes the evolution of these terms. �� � ���� 2�� ��� � 23�� � �� �� � 2 3��� � �� �� � (4) �� � � � ��.�6���� �� ��� � 2�� ����� � 2.64 3 ��� � �� �� � � (5) The last right-hand side term in each of the three main equations is a sink term for its respective quantity. The terms with the quantity � model the production of turbulence by the decay of the tumble macro-vortex during the compression [15]. Simple equation systems for the time rate change of the angular momentum �� ��⁄ model the rotational components of the flow – tumble and swirl – as a single macro-vortex undergoing stretching and compression during the engine intake and compression. Swirl and tumble are produced by the incoming charge, accounting for the measured swirl and tumble coefficients, and reduced by the cylinder out-flow. Equation systems for both rotational motions contain a proper decay functions. Paper [15] does not discuss the swirl decay function, but provides further information on the tumble decay function. Former Flow model ([10], [11]) accounted for the squish motion (inside the swirl model) and injection event kinetic energy and so does the current Flow model. However, the exact equation systems are not available in [15]. Since the current Flow model calculates the kinetic energy and the dissipation rate, the evolution of integral length scale �� over time is then obtained directly with the equation 6 (�� � �.�� is a standard k-ε model constant) [15]. �� � ��� �⁄ �� �⁄ � (6) (6) Equations 1 ‑5 contain four calibration parameters that can be used to match the in ‑cylinder Flow model with 3D ‑CFD results and to enhance the predictive abilities of the consecutive turbulent combustion model. These parameters are following: • Intake term multiplier C1 multiplies the intake term Cin = 0.18C1 and thus accounts for the actual flow velocities through the valves, since these are not equal to the isentropic values assumed by the Ein; • Production term multiplier C2 adjusts the magnitudes of the production terms (equations 4 ‑5) directly through the Production term Cβ with Cβ = 0.38 C2; • Geometrical length scale multiplier C3 adjust the magnitudes of the production terms indirectly, through the scaling of Geometrical length scale Lg = Clen × min(s, 0.5B), with B being the cylinder bore, s the instantaneous piston stroke, and Clen = 0.19C3; • Tumble term multiplier Ctumb controls the intensity of the tumble decay contribution to the production of turbulence. 2.2 TURBULENT COMBUSTION MODEL Turbulent combustion model SITurb predicts the burn rate for the homogeneous charge, respecting the geometry of the combustion chamber, spark location and timing, mixture motion and fuel properties. The model simulates the development of the flame as a turbulent entrainment process followed by a burnup process in a region behind the flame front [8]. Equations 1-5 contain four calibration parameters that can be used to match the in-cylinder Flow model with 3D-CFD results and to enhance the predictive abilities of the consecutive turbulent combustion model. These parameters are following:  Intake term multiplier �� multiplies the intake term ��� � ��1��� and thus accounts for the actual flow velocities through the valves, since these are not equal to the isentropic values assumed by the ���;  Production term multiplier �� adjusts the magnitudes of the production terms (equations 4-5) directly through the Production term �� with �� � ������;  Geometrical length scale multiplier �� adjust the magnitudes of the production terms indirectly, through the scaling of Geometrical length scale �� � ���� � �����������, with � being the cylinder bore, � the instantaneous piston stroke, and ���� � ��1���;  Tumble term multiplier ����� controls the intensity of the tumble decay contribution to the production of turbulence. 2.2 Turbulent combustion model Turbulent combustion model SITurb predicts the burn rate for the homogeneous charge, respecting the geometry of the combustion chamber, spark location and timing, mixture motion and fuel properties. The model simulates the development of the flame as a turbulent entrainment process followed by a burnup process in a region behind the flame front [8]. ��� �� � �� ∙ �� ∙ ��� � ��� (7) The entrainment mass rate of unburned gas ��� ��⁄ is determined by the equation 7; with the flame front area ��, laminar and turbulent flame speeds �� and �� and finally the unburned gas density ��. Dedicated sub-model evaluates the instantaneous flame front area from the combustion chamber geometry assuming the flame front spherical in shape [8]. A model parameter ������ (initial spark size) determines the initial flame front size. For a typical spark plug, its value should be the same as the gap between the spark plug electrodes. However, the real spark size can slightly differ, especially with high-energy spark plugs. Therefore, we use the ������ as a model tuning parameter also, in a reasonable range of sizes. �� � ��� � �� ∙ �� � ����� ∙ � �� ��� � ∙ � ���� � ∙ �1 � ���� ∙ �������∙���� (8) �� � �� ∙ �� ∙ �1 � 1 1 � �� ∙ ��� ���⁄ � (9) During the initial flame kernel development, when the size of the flame kernel is still small, the unburned gas entrainment rate is limited by the laminar flame speed �� (equation 8). Then, the equation 9 accounts for the transition to the turbulent flame speed, with �� representing the mean fluctuating turbulent velocity, �� the flame radius and �� the turbulent length scale [16]. The rate of burnup ��� ��⁄ behind the flame front is proportional to the unburned mass behind the flame front, resulting in the rate equation 10 for the burned mass ��. ��� �� � �� � �� � (10) Model assumes that the burnup phase takes place in the laminar flame speed and over the Taylor microscale of turbulence �, with time constant � (equation 11). Other assumption is that the turbulence is isotropic and therefore the Taylor microscale of turbulence � can be obtained from the integral length scale �� (equations 12). � � ��� (11) � � �� ∙ ������ (12a) ��� � �� ∙ �� ∙ �� � (12b) (7) The entrainment mass rate of unburned gas dMe /dt is determined by the equation 7; with the flame front area Af, laminar and turbulent flame speeds SL and ST and finally the unburned gas density ρu. Dedicated sub ‑model evaluates the instantaneous flame front area from the combustion chamber geometry assuming MECCA 02 2017 PAGE 40 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK the flame front spherical in shape [8]. A model parameter SSinit (initial spark size) determines the initial flame front size. For a typical spark plug, its value should be the same as the gap between the spark plug electrodes. However, the real spark size can slightly differ, especially with high ‑energy spark plugs. Therefore, we use the SSinit as a model tuning parameter also, in a reasonable range of sizes. Equations 1-5 contain four calibration parameters that can be used to match the in-cylinder Flow model with 3D-CFD results and to enhance the predictive abilities of the consecutive turbulent combustion model. These parameters are following:  Intake term multiplier �� multiplies the intake term ��� � ��1��� and thus accounts for the actual flow velocities through the valves, since these are not equal to the isentropic values assumed by the ���;  Production term multiplier �� adjusts the magnitudes of the production terms (equations 4-5) directly through the Production term �� with �� � ������;  Geometrical length scale multiplier �� adjust the magnitudes of the production terms indirectly, through the scaling of Geometrical length scale �� � ���� � �����������, with � being the cylinder bore, � the instantaneous piston stroke, and ���� � ��1���;  Tumble term multiplier ����� controls the intensity of the tumble decay contribution to the production of turbulence. 2.2 Turbulent combustion model Turbulent combustion model SITurb predicts the burn rate for the homogeneous charge, respecting the geometry of the combustion chamber, spark location and timing, mixture motion and fuel properties. The model simulates the development of the flame as a turbulent entrainment process followed by a burnup process in a region behind the flame front [8]. ��� �� � �� ∙ �� ∙ ��� � ��� (7) The entrainment mass rate of unburned gas ��� ��⁄ is determined by the equation 7; with the flame front area ��, laminar and turbulent flame speeds �� and �� and finally the unburned gas density ��. Dedicated sub-model evaluates the instantaneous flame front area from the combustion chamber geometry assuming the flame front spherical in shape [8]. A model parameter ������ (initial spark size) determines the initial flame front size. For a typical spark plug, its value should be the same as the gap between the spark plug electrodes. However, the real spark size can slightly differ, especially with high-energy spark plugs. Therefore, we use the ������ as a model tuning parameter also, in a reasonable range of sizes. �� � ��� � �� ∙ �� � ����� ∙ � �� ��� � ∙ � ���� � ∙ �1 � ���� ∙ �������∙���� (8) �� � �� ∙ �� ∙ �1 � 1 1 � �� ∙ ��� ���⁄ � (9) During the initial flame kernel development, when the size of the flame kernel is still small, the unburned gas entrainment rate is limited by the laminar flame speed �� (equation 8). Then, the equation 9 accounts for the transition to the turbulent flame speed, with �� representing the mean fluctuating turbulent velocity, �� the flame radius and �� the turbulent length scale [16]. The rate of burnup ��� ��⁄ behind the flame front is proportional to the unburned mass behind the flame front, resulting in the rate equation 10 for the burned mass ��. ��� �� � �� � �� � (10) Model assumes that the burnup phase takes place in the laminar flame speed and over the Taylor microscale of turbulence �, with time constant � (equation 11). Other assumption is that the turbulence is isotropic and therefore the Taylor microscale of turbulence � can be obtained from the integral length scale �� (equations 12). � � ��� (11) � � �� ∙ ������ (12a) ��� � �� ∙ �� ∙ �� � (12b) (8) Equations 1-5 contain four calibration parameters that can be used to match the in-cylinder Flow model with 3D-CFD results and to enhance the predictive abilities of the consecutive turbulent combustion model. These parameters are following:  Intake term multiplier �� multiplies the intake term ��� � ��1��� and thus accounts for the actual flow velocities through the valves, since these are not equal to the isentropic values assumed by the ���;  Production term multiplier �� adjusts the magnitudes of the production terms (equations 4-5) directly through the Production term �� with �� � ������;  Geometrical length scale multiplier �� adjust the magnitudes of the production terms indirectly, through the scaling of Geometrical length scale �� � ���� � �����������, with � being the cylinder bore, � the instantaneous piston stroke, and ���� � ��1���;  Tumble term multiplier ����� controls the intensity of the tumble decay contribution to the production of turbulence. 2.2 Turbulent combustion model Turbulent combustion model SITurb predicts the burn rate for the homogeneous charge, respecting the geometry of the combustion chamber, spark location and timing, mixture motion and fuel properties. The model simulates the development of the flame as a turbulent entrainment process followed by a burnup process in a region behind the flame front [8]. ��� �� � �� ∙ �� ∙ ��� � ��� (7) The entrainment mass rate of unburned gas ��� ��⁄ is determined by the equation 7; with the flame front area ��, laminar and turbulent flame speeds �� and �� and finally the unburned gas density ��. Dedicated sub-model evaluates the instantaneous flame front area from the combustion chamber geometry assuming the flame front spherical in shape [8]. A model parameter ������ (initial spark size) determines the initial flame front size. For a typical spark plug, its value should be the same as the gap between the spark plug electrodes. However, the real spark size can slightly differ, especially with high-energy spark plugs. Therefore, we use the ������ as a model tuning parameter also, in a reasonable range of sizes. �� � ��� � �� ∙ �� � ����� ∙ � �� ��� � ∙ � ���� � ∙ �1 � ���� ∙ �������∙���� (8) �� � �� ∙ �� ∙ �1 � 1 1 � �� ∙ ��� ���⁄ � (9) During the initial flame kernel development, when the size of the flame kernel is still small, the unburned gas entrainment rate is limited by the laminar flame speed �� (equation 8). Then, the equation 9 accounts for the transition to the turbulent flame speed, with �� representing the mean fluctuating turbulent velocity, �� the flame radius and �� the turbulent length scale [16]. The rate of burnup ��� ��⁄ behind the flame front is proportional to the unburned mass behind the flame front, resulting in the rate equation 10 for the burned mass ��. ��� �� � �� � �� � (10) Model assumes that the burnup phase takes place in the laminar flame speed and over the Taylor microscale of turbulence �, with time constant � (equation 11). Other assumption is that the turbulence is isotropic and therefore the Taylor microscale of turbulence � can be obtained from the integral length scale �� (equations 12). � � ��� (11) � � �� ∙ ������ (12a) ��� � �� ∙ �� ∙ �� � (12b) Equations 1-5 contain four calibration parameters that can be used to match the in-cylinder Flow model with 3D-CFD results and to enhance the predictive abilities of the consecutive turbulent combustion model. These parameters are following:  Intake term multiplier �� multiplies the intake term ��� � ��1��� and thus accounts for the actual flow velocities through the valves, since these are not equal to the isentropic values assumed by the ���;  Production term multiplier �� adjusts the magnitudes of the production terms (equations 4-5) directly through the Production term �� with �� � ������;  Geometrical length scale multiplier �� adjust the magnitudes of the production terms indirectly, through the scaling of Geometrical length scale �� � ���� � �����������, with � being the cylinder bore, � the instantaneous piston stroke, and ���� � ��1���;  Tumble term multiplier ����� controls the intensity of the tumble decay contribution to the production of turbulence. 2.2 Turbulent combustion model Turbulent combustion model SITurb predicts the burn rate for the homogeneous charge, respecting the geometry of the combustion chamber, spark location and timing, mixture motion and fuel properties. The model simulates the development of the flame as a turbulent entrainment process followed by a burnup process in a region behind the flame front [8]. ��� �� � �� ∙ �� ∙ ��� � ��� (7) The entrainment mass rate of unburned gas ��� ��⁄ is determined by the equation 7; with the flame front area ��, laminar and turbulent flame speeds �� and �� and finally the unburned gas density ��. Dedicated sub-model evaluates the instantaneous flame front area from the combustion chamber geometry assuming the flame front spherical in shape [8]. A model parameter ������ (initial spark size) determines the initial flame front size. For a typical spark plug, its value should be the same as the gap between the spark plug electrodes. However, the real spark size can slightly differ, especially with high-energy spark plugs. Therefore, we use the ������ as a model tuning parameter also, in a reasonable range of sizes. �� � ��� � �� ∙ �� � ����� ∙ � �� ��� � ∙ � ���� � ∙ �1 � ���� ∙ �������∙���� (8) �� � �� ∙ �� ∙ �1 � 1 1 � �� ∙ ��� ���⁄ � (9) During the initial flame kernel development, when the size of the flame kernel is still small, the unburned gas entrainment rate is limited by the laminar flame speed �� (equation 8). Then, the equation 9 accounts for the transition to the turbulent flame speed, with �� representing the mean fluctuating turbulent velocity, �� the flame radius and �� the turbulent length scale [16]. The rate of burnup ��� ��⁄ behind the flame front is proportional to the unburned mass behind the flame front, resulting in the rate equation 10 for the burned mass ��. ��� �� � �� � �� � (10) Model assumes that the burnup phase takes place in the laminar flame speed and over the Taylor microscale of turbulence �, with time constant � (equation 11). Other assumption is that the turbulence is isotropic and therefore the Taylor microscale of turbulence � can be obtained from the integral length scale �� (equations 12). � � ��� (11) � � �� ∙ ������ (12a) ��� � �� ∙ �� ∙ �� � (12b) (9) During the initial flame kernel development, when the size of the flame kernel is still small, the unburned gas entrainment rate is limited by the laminar flame speed SL (equation 8). Then, the equation 9 accounts for the transition to the turbulent flame speed, with u' representing the mean fluctuating turbulent velocity, Rf the flame radius and Lt the turbulent length scale [16]. The rate of burnup dMb /dt behind the flame front is proportional to the unburned mass behind the flame front, resulting in the rate equation 10 for the burned mass Mb. Equations 1-5 contain four calibration parameters that can be used to match the in-cylinder Flow model with 3D-CFD results and to enhance the predictive abilities of the consecutive turbulent combustion model. These parameters are following:  Intake term multiplier �� multiplies the intake term ��� � ��1��� and thus accounts for the actual flow velocities through the valves, since these are not equal to the isentropic values assumed by the ���;  Production term multiplier �� adjusts the magnitudes of the production terms (equations 4-5) directly through the Production term �� with �� � ������;  Geometrical length scale multiplier �� adjust the magnitudes of the production terms indirectly, through the scaling of Geometrical length scale �� � ���� � �����������, with � being the cylinder bore, � the instantaneous piston stroke, and ���� � ��1���;  Tumble term multiplier ����� controls the intensity of the tumble decay contribution to the production of turbulence. 2.2 Turbulent combustion model Turbulent combustion model SITurb predicts the burn rate for the homogeneous charge, respecting the geometry of the combustion chamber, spark location and timing, mixture motion and fuel properties. The model simulates the development of the flame as a turbulent entrainment process followed by a burnup process in a region behind the flame front [8]. ��� �� � �� ∙ �� ∙ ��� � ��� (7) The entrainment mass rate of unburned gas ��� ��⁄ is determined by the equation 7; with the flame front area ��, laminar and turbulent flame speeds �� and �� and finally the unburned gas density ��. Dedicated sub-model evaluates the instantaneous flame front area from the combustion chamber geometry assuming the flame front spherical in shape [8]. A model parameter ������ (initial spark size) determines the initial flame front size. For a typical spark plug, its value should be the same as the gap between the spark plug electrodes. However, the real spark size can slightly differ, especially with high-energy spark plugs. Therefore, we use the ������ as a model tuning parameter also, in a reasonable range of sizes. �� � ��� � �� ∙ �� � ����� ∙ � �� ��� � ∙ � ���� � ∙ �1 � ���� ∙ �������∙���� (8) �� � �� ∙ �� ∙ �1 � 1 1 � �� ∙ ��� ���⁄ � (9) During the initial flame kernel development, when the size of the flame kernel is still small, the unburned gas entrainment rate is limited by the laminar flame speed �� (equation 8). Then, the equation 9 accounts for the transition to the turbulent flame speed, with �� representing the mean fluctuating turbulent velocity, �� the flame radius and �� the turbulent length scale [16]. The rate of burnup ��� ��⁄ behind the flame front is proportional to the unburned mass behind the flame front, resulting in the rate equation 10 for the burned mass ��. ��� �� � �� � �� � (10) Model assumes that the burnup phase takes place in the laminar flame speed and over the Taylor microscale of turbulence �, with time constant � (equation 11). Other assumption is that the turbulence is isotropic and therefore the Taylor microscale of turbulence � can be obtained from the integral length scale �� (equations 12). � � ��� (11) � � �� ∙ ������ (12a) ��� � �� ∙ �� ∙ �� � (12b) (10) Model assumes that the burnup phase takes place in the laminar flame speed and over the Taylor microscale of turbulence λ, with time constant τ (equation 11). Other assumption is that the turbulence is isotropic and therefore the Taylor microscale of turbulence λ can be obtained from the integral length scale Lt (equations 12). Equations 1-5 contain four calibration parameters that can be used to match the in-cylinder Flow model with 3D-CFD results and to enhance the predictive abilities of the consecutive turbulent combustion model. These parameters are following:  Intake term multiplier �� multiplies the intake term ��� � ��1��� and thus accounts for the actual flow velocities through the valves, since these are not equal to the isentropic values assumed by the ���;  Production term multiplier �� adjusts the magnitudes of the production terms (equations 4-5) directly through the Production term �� with �� � ������;  Geometrical length scale multiplier �� adjust the magnitudes of the production terms indirectly, through the scaling of Geometrical length scale �� � ���� � �����������, with � being the cylinder bore, � the instantaneous piston stroke, and ���� � ��1���;  Tumble term multiplier ����� controls the intensity of the tumble decay contribution to the production of turbulence. 2.2 Turbulent combustion model Turbulent combustion model SITurb predicts the burn rate for the homogeneous charge, respecting the geometry of the combustion chamber, spark location and timing, mixture motion and fuel properties. The model simulates the development of the flame as a turbulent entrainment process followed by a burnup process in a region behind the flame front [8]. ��� �� � �� ∙ �� ∙ ��� � ��� (7) The entrainment mass rate of unburned gas ��� ��⁄ is determined by the equation 7; with the flame front area ��, laminar and turbulent flame speeds �� and �� and finally the unburned gas density ��. Dedicated sub-model evaluates the instantaneous flame front area from the combustion chamber geometry assuming the flame front spherical in shape [8]. A model parameter ������ (initial spark size) determines the initial flame front size. For a typical spark plug, its value should be the same as the gap between the spark plug electrodes. However, the real spark size can slightly differ, especially with high-energy spark plugs. Therefore, we use the ������ as a model tuning parameter also, in a reasonable range of sizes. �� � ��� � �� ∙ �� � ����� ∙ � �� ��� � ∙ � ���� � ∙ �1 � ���� ∙ �������∙���� (8) �� � �� ∙ �� ∙ �1 � 1 1 � �� ∙ ��� ���⁄ � (9) During the initial flame kernel development, when the size of the flame kernel is still small, the unburned gas entrainment rate is limited by the laminar flame speed �� (equation 8). Then, the equation 9 accounts for the transition to the turbulent flame speed, with �� representing the mean fluctuating turbulent velocity, �� the flame radius and �� the turbulent length scale [16]. The rate of burnup ��� ��⁄ behind the flame front is proportional to the unburned mass behind the flame front, resulting in the rate equation 10 for the burned mass ��. ��� �� � �� � �� � (10) Model assumes that the burnup phase takes place in the laminar flame speed and over the Taylor microscale of turbulence �, with time constant � (equation 11). Other assumption is that the turbulence is isotropic and therefore the Taylor microscale of turbulence � can be obtained from the integral length scale �� (equations 12). � � ��� (11) � � �� ∙ ������ (12a) ��� � �� ∙ �� ∙ �� � (12b) (11) Equations 1-5 contain four calibration parameters that can be used to match the in-cylinder Flow model with 3D-CFD results and to enhance the predictive abilities of the consecutive turbulent combustion model. These parameters are following:  Intake term multiplier �� multiplies the intake term ��� � ��1��� and thus accounts for the actual flow velocities through the valves, since these are not equal to the isentropic values assumed by the ���;  Production term multiplier �� adjusts the magnitudes of the production terms (equations 4-5) directly through the Production term �� with �� � ������;  Geometrical length scale multiplier �� adjust the magnitudes of the production terms indirectly, through the scaling of Geometrical length scale �� � ���� � �����������, with � being the cylinder bore, � the instantaneous piston stroke, and ���� � ��1���;  Tumble term multiplier ����� controls the intensity of the tumble decay contribution to the production of turbulence. 2.2 Turbulent combustion model Turbulent combustion model SITurb predicts the burn rate for the homogeneous charge, respecting the geometry of the combustion chamber, spark location and timing, mixture motion and fuel properties. The model simulates the development of the flame as a turbulent entrainment process followed by a burnup process in a region behind the flame front [8]. ��� �� � �� ∙ �� ∙ ��� � ��� (7) The entrainment mass rate of unburned gas ��� ��⁄ is determined by the equation 7; with the flame front area ��, laminar and turbulent flame speeds �� and �� and finally the unburned gas density ��. Dedicated sub-model evaluates the instantaneous flame front area from the combustion chamber geometry assuming the flame front spherical in shape [8]. A model parameter ������ (initial spark size) determines the initial flame front size. For a typical spark plug, its value should be the same as the gap between the spark plug electrodes. However, the real spark size can slightly differ, especially with high-energy spark plugs. Therefore, we use the ������ as a model tuning parameter also, in a reasonable range of sizes. �� � ��� � �� ∙ �� � ����� ∙ � �� ��� � ∙ � ���� � ∙ �1 � ���� ∙ �������∙���� (8) �� � �� ∙ �� ∙ �1 � 1 1 � �� ∙ ��� ���⁄ � (9) During the initial flame kernel development, when the size of the flame kernel is still small, the unburned gas entrainment rate is limited by the laminar flame speed �� (equation 8). Then, the equation 9 accounts for the transition to the turbulent flame speed, with �� representing the mean fluctuating turbulent velocity, �� the flame radius and �� the turbulent length scale [16]. The rate of burnup ��� ��⁄ behind the flame front is proportional to the unburned mass behind the flame front, resulting in the rate equation 10 for the burned mass ��. ��� �� � �� � �� � (10) Model assumes that the burnup phase takes place in the laminar flame speed and over the Taylor microscale of turbulence �, with time constant � (equation 11). Other assumption is that the turbulence is isotropic and therefore the Taylor microscale of turbulence � can be obtained from the integral length scale �� (equations 12). � � ��� (11) � � �� ∙ ������ (12a) ��� � �� ∙ �� ∙ �� � (12b) (12a) Equations 1-5 contain four calibration parameters that can be used to match the in-cylinder Flow model with 3D-CFD results and to enhance the predictive abilities of the consecutive turbulent combustion model. These parameters are following:  Intake term multiplier �� multiplies the intake term ��� � ��1��� and thus accounts for the actual flow velocities through the valves, since these are not equal to the isentropic values assumed by the ���;  Production term multiplier �� adjusts the magnitudes of the production terms (equations 4-5) directly through the Production term �� with �� � ������;  Geometrical length scale multiplier �� adjust the magnitudes of the production terms indirectly, through the scaling of Geometrical length scale �� � ���� � �����������, with � being the cylinder bore, � the instantaneous piston stroke, and ���� � ��1���;  Tumble term multiplier ����� controls the intensity of the tumble decay contribution to the production of turbulence. 2.2 Turbulent combustion model Turbulent combustion model SITurb predicts the burn rate for the homogeneous charge, respecting the geometry of the combustion chamber, spark location and timing, mixture motion and fuel properties. The model simulates the development of the flame as a turbulent entrainment process followed by a burnup process in a region behind the flame front [8]. ��� �� � �� ∙ �� ∙ ��� � ��� (7) The entrainment mass rate of unburned gas ��� ��⁄ is determined by the equation 7; with the flame front area ��, laminar and turbulent flame speeds �� and �� and finally the unburned gas density ��. Dedicated sub-model evaluates the instantaneous flame front area from the combustion chamber geometry assuming the flame front spherical in shape [8]. A model parameter ������ (initial spark size) determines the initial flame front size. For a typical spark plug, its value should be the same as the gap between the spark plug electrodes. However, the real spark size can slightly differ, especially with high-energy spark plugs. Therefore, we use the ������ as a model tuning parameter also, in a reasonable range of sizes. �� � ��� � �� ∙ �� � ����� ∙ � �� ��� � ∙ � ���� � ∙ �1 � ���� ∙ �������∙���� (8) �� � �� ∙ �� ∙ �1 � 1 1 � �� ∙ ��� ���⁄ � (9) During the initial flame kernel development, when the size of the flame kernel is still small, the unburned gas entrainment rate is limited by the laminar flame speed �� (equation 8). Then, the equation 9 accounts for the transition to the turbulent flame speed, with �� representing the mean fluctuating turbulent velocity, �� the flame radius and �� the turbulent length scale [16]. The rate of burnup ��� ��⁄ behind the flame front is proportional to the unburned mass behind the flame front, resulting in the rate equation 10 for the burned mass ��. ��� �� � �� � �� � (10) Model assumes that the burnup phase takes place in the laminar flame speed and over the Taylor microscale of turbulence �, with time constant � (equation 11). Other assumption is that the turbulence is isotropic and therefore the Taylor microscale of turbulence � can be obtained from the integral length scale �� (equations 12). � � ��� (11) � � �� ∙ ������ (12a) ��� � �� ∙ �� ∙ �� � (12b) (12b) Parameters in the laminar flame speed equation 8 depend on the fuel type and its composition. Since the composition of the natural gas differs, GT ‑Suite offers two different parameter sets (summed ‑up in Table 1): • First set SLNG1 from Hernandez et al. [17]; • The second one SLNG2 by the work of Lio, Jiang, and Cheng [18] There are five different calibration parameters in the SITurb model that we use to match the measured burn rate. These parameters are following: • Turbulent flame speed multiplier CS scales the turbulent flame speed ST in equation 9; • Flame kernel growth multiplier CK scales the flame front evolution from the initial laminar smooth surface to a distorted turbulent flame front (equation 9); • Taylor length scale multiplier Cλ scales the Taylor microscale of turbulence λ in equation 12 • Dilution exponent multiplier DEM accounts for the dilution by the exhaust residuals and EGR, affecting the laminar flame speed in the equation 8; • Initial spark size SSinit parameter determines the size of the initial flame front. 3. EXPERIMENTAL SET ‑UP AND TEST MATRIX The set of experimental data used in this study originates from a steady state engine test bed measurements with a four‑ ‑cylinder turbocharged SI engine rebuilt from a CI variant and fueled by natural gas. The usual average composition of the natural gas is 98.39 [%vol] CH4, 0.44 [%vol] C2H6, 0.26 [%vol] higher hydrocarbons and 0.84 [%vol] N2. Table 2 summarizes the main geometrical parameters of the experimental engine. One of the necessary inputs for the Flow model is the swirl (or tumble) characteristic of the experimental ICE. Such measurements were conducted in 2005 but only the swirl TABLE 1: Laminar flame speed sub ‑model parameters TABUĽKA 1: Parametre sub ‑modelu pre výpočet laminárnej rýchlosti čela plameňa Parameter Description SLNG1 [17] SLNG2 [18] Bm [m/s] Maximum laminar speed 0.490 0.397 Bϕ [m/s] Laminar speed roll ‑of value – 0.590 – 1.649 ϕm [ ‑] Fuel/air equivalence ratio at maximum laminar flame speed 1.390 1.061 α [ ‑] Temperature exponent 0.68 × ϕ2 – 1.70 × ϕ + 3.18 5.75 × ϕ2 – 12.15 × ϕ + 7.98 β [ ‑] Pressure exponent – 0.52 × ϕ2 + 1.18 × ϕ – 1.18 – 0.925 × ϕ2 + 2 × ϕ – 1.473 MECCA 02 2017 PAGE 41 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK characteristic is available – shown in Figure 1, where CSwirl represents the swirl coefficient (definition from [14]), BS the swirl torque, and Lυ/dυ a ratio of valve lift to its diameter. The experimental ICE is equipped with a central mixer for metering and delivery of the gaseous fuel mixture downstream the compressor inlet. The fuel flow control is either manual or automatic by a closed loop lambda control. ICE features also a cooled low ‑pressure EGR system, with the EGR rate adjusted by a servo driven butterfly valve. Variable turbine geometry performs the boost pressure control and the conventional throttle, located downstream from the intercooler, controls the mixture inflow. High ‑energy ignition system ensures the sufficient spark energy, with the possibility of the spark discharge angle adjustment or closed ‑loop CA50 control. Automated data acquisition system records the engine speed and torque, fuel flow, airflow, exhaust gas composition and average temperatures in the intake and exhaust manifolds. Uncooled piezoelectric transducer installed in the glow plug hole of the first cylinder measures the in ‑cylinder pressure and two piezo resistive pressure transducers measure the intake and exhaust pressures to get a full three ‑pressure‑ ‑analysis (TPA). Details on the experimental set ‑up can be found in [5] and [19]. We used a measurement set containing 83 steady state operation points, representing the full engine map with the stoichiometric mixture and EGR ratio variations (BMEP 4.75‑ ‑19.30 bar; 1200 ‑2600 RPM). Figure 2 shows the reduced test matrix with model calibration points (in blue) and prediction points (in red). The size of a circle and the number indicate the EGR content. Most of the calibration points represent medium ICE loads; low to medium speeds; EGR rates 0 ‑5.6%. Only three ICE operating points contain high EGR rate of 17% and high ICE speeds. The prediction points then cover low load/high load parts of the map, generally with high EGR rates (except two low load points @ 1800 RPM with 0% EGR rate), to really test the predictive capabilities of the combustion model. 4. CALIBRATION PROCEDURE 4.1 BASIS TPA MODEL A proper function of the basis thermodynamic model must be ensured to allow for the calibration of the predictive combustion model. Therefore, we have calibrated the basis TPA model beforehand, correcting some model uncertainties and measurement errors, namely: effective compression ratio, convection multiplier of the heat transfer model [20], TDC positon error, intake and exhaust ports pressure shifts. Figure 2: Test matrix from the full map measurement set with stoichiometric mixture and EGR ratio variation. Obrázok 2: Testovacia matica z merania úplnej charakteristiky so stechiometrickou zmesou a variáciou pomeru EGR We used a measurement set containing 83 steady state operation points, representing the full engine map with the stoichiometric mixture and EGR ratio variations (BMEP 4.75-19.30 bar; 1200-2600 RPM). Figure 2 shows the reduced test matrix with model calibration points (in green) and prediction points (in blue). The size of a circle and the number indicate the EGR content. Most of the calibration points represent medium ICE loads; low to medium speeds; EGR rates 0-5.6%. Only three ICE operating points contain high EGR rate of 17% and high ICE speeds. The prediction points then cover low load/high load parts of the map, generally with high EGR rates (except two low load points @ 1800 RPM with 0% EGR rate), to really test the predictive capabilities of the combustion model. 4. Calibration procedure 4.1 Basis TPA model A proper function of the basis thermodynamic model must be ensured to allow for the calibration of the predictive combustion model. Therefore, we have calibrated the basis TPA model beforehand, correcting some model uncertainties and measurement errors, namely: effective compression ratio, convection multiplier of the heat transfer model [20], TDC positon error, intake and exhaust ports pressure shifts. �� � 1���� � ������ � |����� � ����| ���� ������ �� (13) The determination of the optimal set of calibrated parameters implies the formulation of the objective functions for the whole calibration. In the case of this calibration, these are derived from following model parameters:  Absolute pressure difference �� between the measured and simulated in-cylinder pressures (equation 13);  �������������� evaluated from the GT-Suite output parameter LHV Multiplier. The fuel energy LHV multiplier represents a multiplier of total fuel energy. When its value differs from unity, it indicates, that the input energy in the simulation system is different from the energy needed to follow exactly the measurement in-cylinder pressure trace. Then, �������������� � �������������� � 1�. �� � � �� � ��� �� ������ (14) Both variables are calculated for each calibration engine operation point from the calibration set (Figure 2). The average and maximum values from all calibrated operating points serve as the objective functions (leading to four objective functions �� in total). GA [21] then minimizes these objective functions. The result of a multi-parameter and multi-criterial optimization is a set of non-dominated optimal solutions on a so-called Pareto Frontier. A single optimal solution from the Pareto set is obtained by a criterial function (equation 14), whose value is calculated for each Pareto set solution. The fraction �� ������⁄ than represents a normalization, so that different objective functions �� can be combined into a single equation. ������ is a maximum value from all Pareto set solutions, for the respective objective function �� and parameter �� is a criterial function weight factor. Table 3 summarizes the objective functions �� and values of weight factors ��; Table 4 the selected optimal settings for the basis TPA model. Figure 3 displays a Pareto Frontier for this specific calibration (optimum point in red); Figure 4 the values of �� and �������������� for each calibrated ICE operating point. Objective function �� Weight factor �� Average �� 0.35 Maximum �� 0.15 (13) FIGURE 1: Intake port swirl characteristics of the experimental ICEs OBRÁZOK 1: Swirlová charakteristika sacích kanálov experimentálneho spaľovacieho motora FIGURE 2: Test matrix from the full map measurement set with stoichiometric mixture and EGR ratio variation. OBRÁZOK 2: Testovacia matica z merania úplnej charakteristiky so stechiometrickou zmesou a variáciou pomeru EGR TABLE 2: Basic experimental ICE features TABUĽKA 2: Základné charakteristiky experimentálneho spaľovacieho motora Bore 102 [mm] Stroke 120 [mm] Compression ratio 12:1 Number of Cylinders 4 Valves per Cylinder 4 IVO/IVC 342/595 [°CA aTDC] @ 0.1 mm lift EVO/EVC 123/377 [°CA aTDC] @ 0.1 mm lift Maximum Torque 600 Nm @ 1600 ‑1800 RPM Maximum Power 120 kW @ 2000 RPM MECCA 02 2017 PAGE 42 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK The determination of the optimal set of calibrated parameters implies the formulation of the objective functions for the whole calibration. In the case of this calibration, these are derived from following model parameters: • Absolute pressure difference ∆p between the measured and simulated in ‑cylinder pressures (equation 13); • ∆LHVMultiplier evaluated from the GT ‑Suite output parameter LHVMultiplier. The fuel energy LHVMultiplier represents a multiplier of total fuel energy. When its value differs from unity, it indicates, that the input energy in the simulation system is different from the energy needed to follow exactly the measurement in ‑cylinder pressure trace. Then, ∆LHVMultiplier = |LHVMultiplier – 1|. Both variables are calculated for each engine operation point from the calibration set (Figure 2). The average and maximum values from all calibration points serve as the objective functions (leading to four objective functions Xk in total). GA [21] then minimizes these objective functions. Figure 2: Test matrix from the full map measurement set with stoichiometric mixture and EGR ratio variation. Obrázok 2: Testovacia matica z merania úplnej charakteristiky so stechiometrickou zmesou a variáciou pomeru EGR We used a measurement set containing 83 steady state operation points, representing the full engine map with the stoichiometric mixture and EGR ratio variations (BMEP 4.75-19.30 bar; 1200-2600 RPM). Figure 2 shows the reduced test matrix with model calibration points (in green) and prediction points (in blue). The size of a circle and the number indicate the EGR content. Most of the calibration points represent medium ICE loads; low to medium speeds; EGR rates 0-5.6%. Only three ICE operating points contain high EGR rate of 17% and high ICE speeds. The prediction points then cover low load/high load parts of the map, generally with high EGR rates (except two low load points @ 1800 RPM with 0% EGR rate), to really test the predictive capabilities of the combustion model. 4. Calibration procedure 4.1 Basis TPA model A proper function of the basis thermodynamic model must be ensured to allow for the calibration of the predictive combustion model. Therefore, we have calibrated the basis TPA model beforehand, correcting some model uncertainties and measurement errors, namely: effective compression ratio, convection multiplier of the heat transfer model [20], TDC positon error, intake and exhaust ports pressure shifts. �� � 1���� � ������ � |����� � ����| ���� ������ �� (13) The determination of the optimal set of calibrated parameters implies the formulation of the objective functions for the whole calibration. In the case of this calibration, these are derived from following model parameters:  Absolute pressure difference �� between the measured and simulated in-cylinder pressures (equation 13);  �������������� evaluated from the GT-Suite output parameter LHV Multiplier. The fuel energy LHV multiplier represents a multiplier of total fuel energy. When its value differs from unity, it indicates, that the input energy in the simulation system is different from the energy needed to follow exactly the measurement in-cylinder pressure trace. Then, �������������� � �������������� � 1�. �� � � �� � ��� �� ������ (14) Both variables are calculated for each calibration engine operation point from the calibration set (Figure 2). The average and maximum values from all calibrated operating points serve as the objective functions (leading to four objective functions �� in total). GA [21] then minimizes these objective functions. The result of a multi-parameter and multi-criterial optimization is a set of non-dominated optimal solutions on a so-called Pareto Frontier. A single optimal solution from the Pareto set is obtained by a criterial function (equation 14), whose value is calculated for each Pareto set solution. The fraction �� ������⁄ than represents a normalization, so that different objective functions �� can be combined into a single equation. ������ is a maximum value from all Pareto set solutions, for the respective objective function �� and parameter �� is a criterial function weight factor. Table 3 summarizes the objective functions �� and values of weight factors ��; Table 4 the selected optimal settings for the basis TPA model. Figure 3 displays a Pareto Frontier for this specific calibration (optimum point in red); Figure 4 the values of �� and �������������� for each calibrated ICE operating point. Objective function �� Weight factor �� Average �� 0.35 Maximum �� 0.15 (14) The result of a multi ‑parameter and multi ‑criterial optimization is a set of non ‑dominated optimal solutions on a so ‑called Pareto Frontier. A single optimal solution from the Pareto set is obtained by a criterial function (equation 14), whose value is calculated for each Pareto set solution. The fraction Xk /Xk,max than represents a normalization, so that different objective functions Xk can be combined into a single equation. Xk,max is a maximum value from all Pareto set solutions, for the respective objective function Xk and parameter αk is a criterial function weight factor. Table 3 summarizes the objective functions Xk and values of weight factors αk; Table 4 the selected optimal settings for the basis TPA model. Figure 3 displays a Pareto Frontier for this specific calibration (optimum point in red); Figure 4 the values of ∆p and ∆LHVMultiplier for each calibrated ICE operating point. TABLE 3: Objective functions and weight factors for the basis TPA model calibration TABUĽKA 3: Objektívne funkcie a váhové faktory pre kalibráciu základného modelu TPA Objective function Xk Weight factor αk Average ∆p 0.35 Maximum ∆p 0.15 Average ∆LHVMultiplier 0.35 Maximum ∆LHVMultiplier 0.15 TABLE 4: Selected optimal settings of the basis TPA model TABUĽKA 4: Vybrané optimálne nastavenie základného modelu TPA Full Map set Effective compression ratio 12.36:1 Convection multiplier 1.33 TDC position error 0.1 [°CA] Intake port pressure shift ‑0.021 [bar] Exhaust port pressure shift 0.041 [bar] FIGURE 3: Pareto Frontiers with the optimum point for the basis TPA model calibration. OBRÁZOK 3: Pareto hranice s optimálnym bodom pre kalibráciu základného modelu TPA FIGURE 4: Values of ∆p and ∆LHVMultiplier errors for the full map calibration set and optimal settings of the basis TPA model. OBRÁZOK 4: Hodnoty ∆p a ∆LHVMultiplier základného TPA modelu pre jednotlivé pracovné body experimentálneho spaľovacieho motora a vybrané optimálne nastavenie parametrov. MECCA 02 2017 PAGE 43 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK 4.2 MAIN CALIBRATION AND PREDICTIONS Main calibration of the combustion model is conducted on the calibration part of the full map measurement set, combining in total nine Flow and SITurb parameters (Sections 2.1 and 2.2) and using the SLNG1 parameters for the laminar flame speed sub ‑model. The objective functions for Main calibration of the combustion model are derived from two output parameters: • Absolute pressure difference ∆p between the measured and simulated in ‑cylinder pressures, but now for the SITurb model (equation 13); • Burn Rate RMS Error (GT ‑Suite output parameter). The GA then minimizes four objective functions Xk: two averages and two maxima. Regarding, the values of weight factors αk in the criterial function (equation 14), ∆LHVMultiplier is exchanged for Burn Rate RMS Error. After the calibration of the combustion model, additional prediction points are also simulated. It is worth noting, that the optimal set of model parameters is universal for the whole ICE map, without any dependencies on variables such as ICE speed or ICE load. The same applies for the following optimized set for both sensitivity studies. 4.3 SENSITIVIT Y STUDIES Apart from the main calibration, we have conducted two different sensitivity studies: • Sensitivity 1 on calibration inputs, where we calibrated only the SITurb parameters, with Flow parameters fixed at default values (def = 1); • Sensitivity 2 on laminar flame speed sub ‑model, changing its settings to SLNG2. The calibration procedure, objective functions, and weight factors of the criterial function are the same for the Sensitivity 1 and Sensitivity 2 as for the Main calibration. 5. RESULTS AND DISCUSSION 5.1 MAIN CALIBRATION After the main calibration of the Flow and SITurb models, we kept constant the optimized parameters (Table 5, first column) and simulated both the 14 calibration operating points and additional 16 prediction points. Table 6 than summarizes the average and maximum values (Note: the average errors are evaluated from absolute values for the individual operating points). Figure 5 shows the IMEP percentage error between the experimental and simulated values; Figure 6 displays the CA50 error; Figure 7 the 10% ‑90% Burn Duration (MFB10 ‑90) error; Figure 8 maximum firing pressure error; Figure 9 the error of maximum firing pressure CA position. TABLE 5: Optimal values of the calibration parameters for the combustion model (Main calibration, Sensitivity 1, Sensitivity 2) TABUĽKA 5: Výsledné optimálne hodnoty kalibračných parametrov pre model horenia (Hlavná kalibrácia, Citlivosť 1, Citlivosť 2) Parameter Main calibration Sensitivity 1 Sensitivity 2 Turbulent Flame Speed Multiplier CS 1.060 0.370 1.600 Flame Kernel Growth Multiplier CK 9.040 4.210 0.080 Taylor Length Scale Multiplier Cλ 7.510 2.650 8.970 Dilution Exponent Multiplier DEM 0.830 0.830 0.710 Initial Spark Size SSinit 3.500 3.560 4.810 Intake Term Multiplier C1 1.640 def = 1 0.010 Production Term Multiplier C2 3.690 def = 1 0.946 Geometrical Length Scale Multiplier C3 0.070 def = 1 0.050 Tumble Term Multiplier Ctumb 0.300 def = 1 1.410 TABLE 6: Main calibration average and maximum errors (calibration/prediction points) TABUĽKA 6: Priemerné a maximálne odchýlky Hlavnej kalibrácie (kalibračné/predikčné body) Parameter Avg. error Max. error IMEP 1.35/0.95 [%] 4.75/2.84 [%] CA50 0.44/0.45 [°CA] 1.28/0.68 [°CA] MFB10 ‑90 0.52/0.82 [°CA] 2.18/1.14 [°CA] MFB10 ‑75 1.17/2.17 [°CA] ‑2.46/0.53 [°CA] Maximum Pressure 1.20/2.46 [bar] 3.94/5.78 [bar] CA @ Maximum Pressure 0.42/0.54 [°CA] 1.10/1.10 [°CA] FIGURE 5: IMEP percentage errors of the experimental versus simulation values, Main calibration OBRÁZOK 5: Percentuálna chyba IMEP, experiment verzus simulácia MECCA 02 2017 PAGE 44 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK FIGURE 8: Maximum firing pressure errors of the experimental versus simulation values, Main calibration OBRÁZOK 8: Chyba maximálneho tlaku, experiment verzus simulácia, Hlavná kalibrácia FIGURE 9: Maximum firing pressure crank angle position errors of the experimental versus simulation values, Main calibration OBRÁZOK 9: Chyba polohy maximálneho tlaku, experiment verzus simulácia, Hlavná kalibrácia FIGURE 6: CA50 errors of the experimental versus simulation values, Main calibration OBRÁZOK 6: Chyba CA50, experiment verzus simulácia, Hlavná kalibrácia FIGURE 7: Burn Duration 10% ‑90% errors of the experimental versus simulation values, Main calibration OBRÁZOK 7: Chyba dĺžky horenia 10 ‑90% zhoreného paliva, experiment verzus simulácia, Hlavná kalibrácia MECCA 02 2017 PAGE 45 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK Both, calibration and prediction operating points show very good overall agreement, supported by the visual comparison of the burn rate and in ‑cylinder pressure traces (e.g. Figure 10 and Figure 11). Figure 10 shows the burn rate and in ‑cylinder pressure comparison for the operating point 81 (prediction set) with the worst IMEP error (2.84%) and Figure 11 for operating point 28 (calibration set). In conclusion, Table 5 shows that values for some model parameters, such as CK, Cλ, and C1 are quite high, on the other hand the value of C3 is low. This means, that although the agreement of the overall parameters (in ‑cylinder pressure, IMEP burn rate) is very good, a comparison with 3D ‑CFD is necessary. The value for the initial spark size SSinit is reasonable. 5.2 SENSITIVIT Y 1: CALIBRATION INPUTS The first sensitivity compares the Main calibration to reduced calibration set, with Flow model parameters fixed on default values (def = 1). GA was only optimizing the SITurb values. Optimal SITurb parameters are summarized in the second column of Table 5. It is important to note that the optimal values for both sets (Main Calibration and Sensitivity 1) are comparable to each ‑other, with a possible trade ‑off effect between the SITurb parameter CS (Turbulent Flame Speed Multiplier) and Flow parameter C2 (Production Term Multiplier). Table 7 sums ‑up the average and maximum error values. The overall agreement is also very good, but compared to the Main Calibration results (Table 6), both the average and maximum errors are higher. To illustrate the effects of the Flow model parameters, Figure 12 shows the difference between the turbulent kinetic energy and turbulent flame speeds. The Flow model parameters in the Main Calibration actually dampen the in‑ ‑cylinder flow motion, which is compensated by the turbulent velocities. 5.3 SENSITIVIT Y 2: LAMINAR FLAME SPEED MODEL Sensitivity 2 deals about the effect of the laminar flame speed model, which is set to SLNG2 values. The third column of Table 5 shows the optimized values for both the Flow and SITurb sub ‑models. In this case, compared to the Main Calibration outputs, the differences are notable. Table 8 summarizes the average and maximum error values. The average and maximum error values are higher than in the Sensitivity 1. And especially those for MFB10 ‑75% show the effect of the different laminar flame speed models. Figure 13 than depicts the burn rate and in ‑cylinder pressure comparison of experimental values, Main Calibration, and Sensitivity 2. The selected low load operating point 41 is taken from the prediction set (achieves a low overall error in the Main Calibration). TABLE 7: Sensitivity 1 average and maximum errors (calibration/prediction points) TABUĽKA 7: Priemerné a maximálne odchýlky Citlivosti 1 (kalibračné/predikčné body) Parameter Avg. error Max. error IMEP 1.15/1.08 [%] 5.62/2.71 [%] CA50 0.49/0.59 [°CA] ‑1.74/ ‑1.52 [°CA] MFB10 ‑90 0.79/1.03 [°CA] 2.50/ ‑1.72 [°CA] MFB10 ‑75 3.04/3.71 [°CA] ‑4.49/ ‑5.44 [°CA] Maximum Pressure 1.30/3.13 [bar] 4.49/6.96 [bar] CA @ Maximum Pressure 0.55/0.72 [°CA] ‑1.40/ ‑1.40 [°CA] FIGURE 12: Comparison between the experimental and simulation of turbulent kinetic energy and turbulent flame speed at prediction operating point 65 (2000 RPM, BMEP 16.61 bar, 10.7% EGR) OBRÁZOK 12: Porovnanie experimentálneho a simulačného priebehu turbulentnej kinetickej energie a turbulentnej rýchlosti čela plameňa pre predikčný pracovný bod 65 (2000 RPM, BMEP 16.61 bar, 10.7% EGR) << FIGURE 10: Comparison between the experimental and simulation burn rate and in ‑cylinder pressure at 2600 RPM, BMEP 12.98 bar, 17.3% EGR (operating point 81 maximum IMEP error for the prediction set) OBRÁZOK 10: Porovnanie experimentálneho a simulačného priebehu rýchlosti horenia a tlaku vo valci pri 2600 RPM, BMEP 12.98 bar, 17.3% EGR (pracovný bod 81 s maximálnou percentuálnou chybou IMEP, z predikčného setu) < FIGURE 11: Comparison between the experimental and simulation burn rate and in ‑cylinder pressure at 1600 RPM, BMEP 15.59 bar, 2.3% EGR (operating point 28, calibration set) OBRÁZOK 11: Porovnanie experimentálneho a simulačného priebehu rýchlosti horenia a tlaku vo valci pri 1600 RPM, BMEP 15.59 bar, 2.3% EGR (pracovný bod 28, kalibračný set) MECCA 02 2017 PAGE 46 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK 6. CONCLUSION We have evaluated the predictive capabilities of a 0D phenomenological in ‑cylinder flow model, based on the K ‑k ‑ε kinetic energy cascade approach and coupled with a turbulent combustion model, using a turbocharged natural gas SI engine. First, we did a detailed model calibration using GA and then added two model sensitivity studies: regarding the calibration procedure and the laminar flame speed sub ‑model. The main detailed model calibration shows that: • Very good agreement with the experimental data can be achieved on the side of in ‑cylinder pressure traces and burn rate; • The combustion model is capable of prediction outside of its calibration range, achieving good results; • Some of the calibration parameters in the optimal set are high, which has to be further studied (e.g. comparison with 3D ‑CFD) The results from the first sensitivity on the calibration inputs (calibrating only the turbulent combustion model; the flow model set to default values) only strengthen the conviction of the necessary comparison with 3D ‑CFD: • The available results from the flow model are different than those in the main calibration; • However, if the flow model is not included in the calibration, the combustion model still shows good agreement with the experimental data and prediction abilities. The second calibration – concerning the laminar flame speed model – shows the importance of a correct model values: • We have tested two different sets of laminar flame speed sub ‑model values and both showing differences; • The overall effect on the burn rate and in ‑cylinder pressure traces is greater than in the first sensitivity. In conclusion, our work shows the importance of the in ‑cylinder flow model verification with the 3D ‑CFD, and the importance of a correct laminar flame speed model, especially for the natural gas fueled ICE. Future development concerning the natural gas ICE will focus on a comparison of the 0D phenomenological in ‑cylinder flow model with a 3D ‑CFD and further extension of the test matrix. The extended test matrix will include the air dilution and spark ‑timing sweeps. After the studies on the natural gas engine, the work will move to a gasoline SI ICE also. ACKNOWLEDGEMENTS This work was supported by the Grant Agency of the Czech Technical University in Prague, grant No. SGS16/213/OHK2/3T/12 LIST OF ABBREVIATIONS 0D/3D Zero/Three ‑Dimensional BMEP Brake Mean Effective Pressure CA Crank Angle CFD Computational Fluid Dynamics DoE Design of Experiments EGR Exhaust Gas Recirculation EVC Exhaust Valve Close EVO Exhaust Valve Open GA Genetic Algorithms ICE Internal Combustion Engine IMEP Indicated Mean Effective Pressure IVO Intake Valve Open IVC Intake Valve Close MFB Mass Fraction Burned RMS Root Mean Square SI Spark Ignition TDC Top Dead Center TPA Three ‑Pressure ‑Analysis TABLE 8: Sensitivity 2 average and maximum errors (calibration/prediction points) TABUĽKA 8: Priemerné a maximálne odchýlky Citlivosti 2 (kalibračné/predikčné body) Parameter Avg. error Max. error IMEP 0.69/0.92 [%] ‑1.97/ ‑1.93 [%] CA50 0.67/0.68 [°CA] 2.60/1.93 [°CA] MFB10 ‑90 0.37/0.41 [°CA] 0.84/0.97 [°CA] MFB10 ‑75 2.10/3.03 [°CA] ‑5.71/ ‑5.38 [°CA] Maximum Pressure 1.48/2.08 [bar] ‑3.69/ ‑5.63 [bar] CA @ Maximum Pressure 0.64/0.92 [°CA] 1.50/1.90 [°CA] FIGURE 13: Comparison between the experimental and simulation burn rate and in ‑cylinder pressure at prediction operating point 41 (1800 RPM, BMEP 8.62 bar, 0.0% EGR) OBRÁZOK 13: Porovnanie experimentálneho a simulačného priebehu rýchlosti horenia a tlaku vo valci pre predikčný pracovný bod 41 (1800 RPM, BMEP 8.62 bar, 0.0% EGR) MECCA 02 2017 PAGE 47 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK LIST OF SYMBOLS 𝐴𝑓 Flame Front Area 𝐵 Cylinder Bore 𝐵 𝑚 Maximum Laminar Speed 𝐵 𝜙 Laminar Speed Roll ‑of Value 𝐶1 Intake Term Multiplier 𝐶2 Production Term Multiplier 𝐶3 Geometrical Length Scale Multiplier 𝐶𝑖𝑛 Intake Term 𝐶𝐾 Flame Kernel Growth Multiplier 𝐶𝑆 Turbulent Flame Speed Multiplier 𝐶𝑆𝑤𝑖𝑟𝑙 Intake Port Swirl Coefficient 𝐶𝑡𝑢𝑚𝑏 Tumble Term Multiplier 𝐶𝛽 Production Term 𝐶𝜆 Taylor Length Scale Multiplier 𝐶𝜇 k ‑ε Model Constant 𝐷𝐸𝑀 Dilution Exponent Multiplier 𝐷𝑖𝑙 Mass Fraction of the Residuals in the Unburned Zone 𝐸𝑖𝑛 Intake Energy 𝐹 Criterial Function 𝑘 Turbulent Kinetic Energy 𝐾 Mean Kinetic Energy 𝐿 Angular Momentum 𝐿𝑔 Geometrical Length Scale 𝐿𝑡 Integral Length Scale 𝑚 In ‑cylinder Mass 𝑚· 𝑜𝑢𝑡 Cylinder Exit Mass Flow Rate 𝑀𝑏 Burned Mass 𝑀𝑒 Entrained Mass 𝑝 Pressure 𝑝0 Reference Pressure 𝑃𝑘 Turbulence Production Term 𝑃𝜖 Dissipation Rate Production Term 𝑅𝑓 Flame Front Radius 𝑅𝑒 Reynolds Number 𝑠 Piston Stroke 𝑆𝐿 Laminar Flame Speed 𝑆𝑇 Turbulent Flame Speed 𝑆𝑆𝑖𝑛𝑖𝑡 Initial Spark Size 𝑇𝑢 Temperature of Unburned Gas 𝑇0 Reference Temperature 𝑇 Tumble Contribution to Turbulence 𝑢' Mean Fluctuating Turbulent Velocity 𝑈 Mean Velocity inside the Cylinder 𝑋𝑘 Objective Function 𝛼 Temperature Exponent 𝛼𝑖𝑛 Intake Energy Fraction Converted Directly into Turbulence 𝛼𝑘 Weight Factor 𝛽 Pressure Exponent Δ Difference 𝜖 Turbulent Dissipation Rate 𝜆 Taylor Microscale of Turbulence 𝜈𝑇 Turbulent Viscosity 𝜌 Density 𝜌· Density Rate of Change 𝜌𝑢 Density of Unburned Gas 𝜙 Fuel/air Equivalence Ratio 𝜙𝑚 Fuel/air Equivalence Ratio at Maximum Laminar Flame Speed REFERENCES [1] Vitek, O., Macek, J., Tatschl, R., Pavlovic, Y., Priesching, P. (2012). LES Simulation of Direct Injection SI ‑Engine In ‑Cylinder Flow. SAE Technical Paper 2012 ‑01 ‑0138, doi: 10.4271/2012 ‑01 ‑0138 [2] Vibe, I.I. (1956). Semi ‑empirical expression for combustion rate in engines, In: Proc. Conference on Piston Engines, USSR Academy of Sciences, Moscow [3] Csallner, P., Woschni, G. (1982). Zur Vorausberechnung des Brennverlaufes von Ottomotoren bei geanderten Betriebsbedingungen, MTZ, No.5 [4] Vavra, J., Takats, M. (2004). Heat Release Regression Model for Gas Fuelled SI Engines, SAE Technical Paper 2004 ‑01 ‑1462, doi: 10.4271/2004 ‑01 ‑1462 [5] Skarohlid, M. (2010). Modeling of Influence of Biogas Fuel Composition on Parameters of Automotive Engines, SAE Technical Paper 2010 ‑01 ‑0542, doi: 10.4271/2010 ‑01 ‑0542 [6] Macek, J., Steiner, T. (1995). Advanced Multizone Multidimensional Models of Engine Thermodynamics. In: 21st CIMAC Congress, London, p. D10/1 ‑D10/18 [7] Hvezda, J. (2014). Multi ‑Zone Models of Combustion and Heat Transfer Processes in SI Engine. SAE Technical Paper 2014 ‑01 ‑1067, doi: 10.4271/2014 ‑01 ‑1067 [8] Wahiduzzaman, S., Morel, T., Sheard, S. (1993). Comparison of Measured and Predicted Combustion Characteristics of a Four ‑Valve S.I. Engine, SAE Technical Paper 930613, doi: 10.4271/9306 [9] Morel, T., Mansour, N. N. (1982). Modeling of turbulence in Internal Combustion Engines, SAE Technical Paper 820040, doi: 10.4271/820040 [10] Morel, T., Keribar, R. (1985). A Model for Predicting Spatially and Time Resolved Convective Heat Transfer in Bowl ‑in ‑Piston Combustion Chambers, SAE Technical Paper 850204, doi: 10.4271/850204 MECCA 02 2017 PAGE 48 Evaluation of the Predictive Capabilities of a Phenomenological Combustion Model for Natural Gas SI Engine RASTISLAV TOMAN, JAN MACEK [11] Morel, T., Rackmil, C. I., Keribar, R., Jennings, M. J. (1988). Model for Heat Transfer and Combustion in Spark Ignited Engines and Its Comparison with Experiments, SAE Technical Paper 880198, doi: 10.4271/880198 [12] Bozza, F., Gimelli, A. (2004). A Comprehensive 1D Model for the Simulation of a Small ‑Size Two‑ ‑Stroke SI Engine. SAE Technical Paper 2004 ‑01 ‑0999, doi: 10.4271/2004 ‑01 ‑0999 [13] Vitek, O., Macek, J., Poetsch, Ch., Tatschl, R. (2013). Modeling Cycle ‑to ‑Cycle Variations in 0 ‑D/1‑ ‑D Simulation by Means of Combustion Model Parameter Perturbations based on Statistics of Cycle‑ ‑Resolved Data. SAE Int. J. Engines, Vol. 6, No. 2, doi: 10.4271/2013 ‑01 ‑1314 [14] GT ‑POWER Engine Performance Application Manual. [Manual] Westmont: Gamma Technologies, Inc., 2015 [15] Fogla, N., Bybee, M., Mirzaeian, M., Millo, F. et al. (2017). Development of a K ‑k ‑ε Phenomenological Model to Predict In ‑Cylinder Turbulence, SAE Int. J. Engines, Vol. 10, No. 2, doi: 10.4271/2017 ‑01 ‑0542 [16] Mirzaeian, M., Millo, F., and Rolando, L. (2016). Assessment of the Predictive Capabilities of a Combustion Model for a Modern Downsized Turbocharged SI Engine, SAE Technical Paper 2016 ‑01 ‑0557, doi: 10.4271/2016 ‑01 ‑0557 [17] Hernandez, J., Lapuerta, M., Serrano, C. (2005) Estimation of the Laminar Flame Speed of Producer Gas from Biomass Gasification, Energy Fuels, Vol. 19, No. 5, p. 2172 ‑2178, doi: 10.1021/ef058002y [18] S.Y. Liao, D.M. Jiang, and Q. Cheng. (2004) Determination of laminar burning velocities for natural gas, Fuel, Vol. 83, Issue 9, p.1247 ‑1250. doi: 10.1016/j.fuel.2003.12.001 [19] Vavra, J., Takats, M., Klir, V., Skarohlid, M. (2012). Influence of Natural Gas Composition on Turbocharged Stoichiometric SI Engine Performance, SAE Technical Paper 2012 ‑01 ‑1647, doi: 10.4271/2012 ‑01 ‑1647 [20] Woschni, G. (1967). A Universally Applicable Equation for the Instantaneous Heat Transfer Coefficient in the Internal Combustion Engine, SAE Technical Paper 880198, doi: 10.4271/670931 [21] modeFRONTIER – Multi ‑Objective Design Environment, version 4.4.3. [CD ‑ROM], 2012 _GoBack