© R. Vihar et al., published by IFP Energies nouvelles, 2018

1 Introduction

The increasing share of diesel engine powered automotive fleet in recent years, being largely due to high thermal efficiency of diesel engines, has become an environmental issue since it represents a major source of NO_x emissions [1]. One of the ways to comply with forthcoming increasingly stringent diesel emission limits for NO_x, which will comprise also real driving emission, is adoption of engine control strategies on the real-time (RT) basis. Accurate estimation of in-cylinder NO_x formation on the cycle-to-cycle basis can namely be used as an input parameter for active control methods (like Closed Loop Combustion Control) and also for the control of exhaust after-treatment systems.

Some modern diesel engines are equipped with in-cylinder pressure transducers, which can be used as a source of combustion related data [2, 3]. Currently in-cylinder pressure data is mainly used for obtaining information related to aging of engine components as well as for adaptation of the injection strategy with respect to the variation in the fuel quality [4]. Information on the in-cylinder pressure opens many other possibilities for more precise control of the engine especially in transient operating conditions, as it can be used for calculation of the thermodynamic parameters inside the cylinder being an enabler for the virtual NO_x sensor with higher level of predictiveness.

In general, methods for simulations of pollutants in the RT are usually divided into semi-empirical and phenomenological approaches [5]. In the semi-empirical models where NO_x formation process is based on a single-equation correlation which takes into account measured engine variables, the accuracy largely depends on the tuning of the model. Therefore their prediction quality outside the calibration area is very poor [6]. Hence, phenomenological models, e.g. [5, 7], can be considered as a preferred option in terms of general applicability of the virtual NO_x sensors. Thus the advanced virtual NO_x sensor presented in this paper is also based on the phenomenological modeling approach.

Main thermodynamic parameters that are required for prediction of the NO_x formation are rate of heat release (ROHR) and consequently the temperatures and concentrations in the combustion chamber or selected zones of the combustion chamber. These properties can be assessed with thermodynamic models whose accuracies depend on the physical depth, which is directly reflected in the computational complexity of the models [8]. Since the whole control loop of the engine (signal processing, thermodynamic model and potential control of actuators) should be in the automotive application carried out by the engine control unit (ECU), design of a computationally efficient thermodynamic model is of utmost importance.

Combustion of non-homogeneous mixture of fuel, air and recirculated exhaust gases is characterized by a high heterogeneity of the temperature and concentration field, which have high impact on NO_x generation. The most accurate methodology to simulate such phenomena is based on the use of 3D CFD models [9]. Due to their high spatial resolution these models feature a high computational complexity and thus long computing times. On the system level, 0-D combustion models, e.g. [10, 11], are often used as models that effectively combine the requirements on multi-zone spatial resolution and much shorter computational times compared to the 3D CFD models. However, multi zone 0-D models featuring many zones to resolve spray propagation are still too complex and thus computationally demanding for RT applications [12]. On the other hand, in the simplest models the entire combustion area is treated as a single zone (single-zone models) [13] and the temperature could be treated as a global variable that is uniform throughout the combustion chamber. As formation of thermal NO_x is very sensitive on temperature, the global in-cylinder temperature provided by the single zone models is not a suitable input for virtual NO_x sensors [14].

Two-zone models can, on the system level, efficiently comply with the requirements on heterogeneity of the temperature and concentration field as well as short computational times that are in the range of single zone models. Two-zone models thus separate the combustion chamber into burned and unburned zone while considering zone specific temperatures and concentrations as well as homogeneity of the pressure field [15, 16]. Three-zone models are further variants of two-zone models, where by adding complexity, computational complexity is also increased [17]. The model proposed in this paper was developed with the objective on very low computational demand and thus the two-zone combustion modeling approach offered a proper basis for achieving good trade-off between accuracy and computational complexity at low computational demand.

The two-zone models mostly differ in terms of complexity of the gas property treatment. Besides the inputs of gas properties this also influences the complexity of governing equations through the inclusion of the partial derivatives of internal energy and gas constant. Some of the RT two-zone models are using constant gas properties in the calculation of the thermodynamic parameters [18, 19]. This modeling depth might not provide sufficiently profound basis for achieving high accuracy of the results as will be shown in the Result section of this paper. Therefore, variable gas properties need to be considered to comply with the requirement on sufficient accuracy of the results. This calls for an efficient approach in terms of complexity of the gas property database and of the underlying governing equations, where in addition the computationally efficient treatment of gas property database is preferred. The modeling approach proposed in this paper thus considers dependency of the gas properties on temperature and relative air-fuel ratio. In addition, the gas properties are not given as polynomial functions, as for example in [17], which are computationally demanding and might not allow for achieving high accuracy especially if lower order polynomials are applied. In the proposed model, main thermodynamic properties, their derivatives and species equilibrium fractions are therefore stored in the equidistant look-up tables being a computation efficient approach since commercial engine ECUs mostly use the look-up tables for modeling operating point dependent parameters [20]. Moreover, the thermodynamic modeling framework was developed in a way that the thermodynamic equations could be solved in a closed form, without requiring an iterative procedure. This means that governing equations are explicit in mass of burned fuel, while considering two-zone combustion model and variable gas properties as well as using an in-cylinder pressure as an input. This thermodynamic modeling framework is crucial to effectively combine computational efficiency and sufficient modeling depth to achieve high level of predictiveness and it is more detailed compared to other models that were applied for real-time virtual NO_x sensors [18, 21].

The two-zone models also differ in the approaches for determining the evolution of relative air-fuel ratio during combustion (RAFC), which influences the mass, species and enthalpy transport from the unburned to the burned zone. The first group of the models is based on the assumption that at the end of combustion process the entire unburned zone is consumed and thus transferred to the burned zone [22], while the second group of models is based on the assumption that fuel and air are burning at stoichiometric conditions [8, 15, 23]. Both assumptions are not necessarily relevant for the entire operating range of the diesel engine. This can, besides omitting contribution of the prompt NO mechanism whose relative contribution is more pronounced at low engine loads [6, 21, 24, 25], be related to the inability of predicting NO_x emissions using extended Zeldovich mechanism with sufficient accuracy without using operating point specific NO_x multipliers [8]. As these NO_x multipliers to a large extend correct also for the temperature trace in the burned zone − being significantly influenced by the previously mentioned approaches to determine the RAFC − this paper introduces a physically motivated degree of freedom in the calibration of combustion models being the RAFC. It is also important to mention that in some studies other thermodynamic calibration parameters, e.g. heat transfer multipliers, were used to calibrate the NO_x mass [21, 24]. Heat transfer multiplier inherently influences the burned zone temperature and thus also NO_x mass, however it also has a direct impact on the energy balance and thus a significant impact on other thermodynamic results as MFB. This might lead to inconsistencies, if heat transfer multiplier that best suits for calibration of the NO_x mass violates the energy balance. Therefore, the proposed modeling framework is derived in a way that NO_x calibration does not violate the energy balance. This is made possible by the application of the RAFC as the main calibration parameter as it influences the transport from the unburned to the burned zone and thus the burned zone temperature and NO_x formation. It also has a nearly negligible influence on the overall internal energy of the burned and the unburned zone and thus on the energy balance. Later is true for lean mixtures as overall internal energy of the burned and the unburned zone is only very insignificantly influenced by the composition [26]. RAFC is thus the only independent calibration parameter of the NO_x sub-model which has minor impact on energy balance. In addition, unlike the approach proposed in [22], the approach proposed in this paper does not require that entire unburned zone is converted to the burned zone at the end of combustion. This is very important to obtain plausible results of NO_x especially at low loads, where indeed only a part of the in-cylinder charged is mixed with the fuel and thus being subjected to combustion.

The proposed two-zone thermodynamic modeling framework thus considers an additional thermodynamic NO_x calibration parameter, which differs from the two currently used approaches for determining the RAFC that either rely on constant RAFC or lead to full conversion of the unburned zone. This allows applying extended Zeldovich mechanism with fixed model parameters over the entire engine operating range. Moreover, it was shown that it is also possible to obtain plausible NO_x results using only a single RAFC input for the entire engine operating regime thus omitting any operating point specific tuning of the air-fuel ratio evolution during combustion.

The entire modeling framework, i.e. thermodynamic and NO_x model, was validated in multiple operating points of two significantly different diesel engines demonstrating its general applicability. In addition, it was validated using pressure signals from a laboratory pressure sensor and a large series low-cost glow plug pressure sensor. The model was run both in off-line analyses on PC and in on-line analyses utilizing field programmable gate array (FPGA). The proposed modeling framework can thus be considered as an approach offering low parametrization effort and optimized trade-off between accuracy and calculation time thus opening ways towards its use in the next generation large series commercial hardware.

2 Experimental setup

The experiments were carried out on a two different diesel engines covering light duty (LD) and heavy duty (HD) segment. Both engines were coupled with a Zöllner B-350AC eddy-current dynamometer controlled by Kristel, Seibt & Co control system KS ADAC. A Kistler CAM UNIT Type 2613B shaft encoder provided an external trigger and an external clock (0.1–6 deg CA) for data acquisition system. In-cylinder pressure was measured with two piezo-electric pressure transducers: calibrated AVL GH12D transducer in combination with charge amplifier AVL MICROIFEM and pre-series Hidria glow plug pressure transducer (PSG) with integrated charge amplifier [27], both connected to 16 bit, 4 channel National Instruments data-acquisition system with maximum sampling frequency 1 MS/s/ch. Top dead center (TDC) was determined by capacitive sensor COM Type 2653. An AVL 730 gravimetric balance fuel meter was employed to measure fuel consumption. Intake air flow was measured with Meriam laminar flow meter 50MC2-6F. Emissions of NO_x were measured with Horiba OBS-2200 with the use of CLD analyzer. Sampling attachments were welded on both engines in the exhaust system after the turbine and before after-treatment devices. In the development phase of the model the laboratory equipment was used for evolution of the proposed modeling approach, while validation was done only with the use of pressure signal from PSG in combination with data provided by the ECU which demonstrates the applicability of the proposed methods for RT use.

NO_x emissions were measured continuously while in-cylinder pressure was recorded over 100 successive cycles at sampling resolution of 0.2 CA for the HD and 0.1 CA for the LD engine. For off-line application, a representative pressure trace for the model development was generated by averaging 100 cycles of the individual operating point measured with AVL pressure transducer in order to partially suppress measurement noise from the measured pressure trace. Furthermore, additional filtering of the pressure trace was applied, as suppression of oscillations of the pressure signal is very important in this particular study due to the high impact of pressure oscillations on the burn rate, which strongly influences temperature in the burned zone and consequently NO_x concentration. Therefore, pressure traces were further processed by applying low pass FIR filter [28], whereby filter parameters (e.g. filter order, cut-off frequency) were carefully selected in a way to suppress unwanted contributions of vibrational eigen-modes in the combustion chamber and of the measurement noise while preserving contributions of piston kinematics and combustion [29, 30]. For validation of on-line application, 1 cycle was measured with PSG pressure transducer to validate applicability of the method in most severe real-time applications. In on-line application the same filter settings as in off-line application were applied.

Besides different engine application areas, i.e. light and heavy duty, analyzed engines feature also different emission standards and different fuel injection systems, which confirms validity of the proposed in-cylinder NO_x formation modeling approach. The LD engine was a 1.6 L, 4-cylinder, PSA turbocharged diesel engine equipped with a common rail fuel system in combination with solenoid injectors with injection pressure up to 1600 bar. The main characteristics of the LD engine are given in Table 1. The HD engine was a 6-cylinder, 4 stroke, turbocharged, 6.87 L MAN diesel engine. The main characteristics of the second engine are given in Table 2.

The proposed modeling approach was implemented in the LabVIEW environment. The model was designed in a modular way and allows for variations in modeling complexity to support analyses of the trade-off between computational speed and accuracy of the model. Off-line validation of the proposed modeling approach was performed on a standard PC, whereas on-line validation was performed on a Xilinx Kintex-7 325T FPGA integrated on a NI cRIO Controller. The matrix of operating points that was used for model validation is shown in Tables 3 and 4 for LD and HD engine, respectively.

3 Description of the modeling framework

Until the start of combustion (SOC), the in-cylinder gas mixture is modeled as a single zone model described in Ref. [31] and mean in-cylinder temperature is determined considering ideal and non-perfect gas. After the SOC, combustion chamber is virtually divided into two zones assuming homogeneous in-cylinder pressure distribution. SOC was identified by comparing the motored and the measured in-cylinder pressure after the start of injection (SOI) [6]. The crank angle resolution is 1 deg CA with the exception of first 30 deg CA after the SOC, where it is set to 0.1 deg CA resolution in order to prevent potential instability of the model in the initial phase of the two-zone calculation.

3.1 Two-zone combustion modeling framework

The two-zone modeling framework is derived based on the mass and energy balances applied to each of the zones. The burned zone is denoted with subscript B and the unburned zone with subscript U. In general, n-zone modeling frameworks are derived in a way to predict pressure derivative based on the known ROHR, e.g. [32, 33], whereas in the RT modeling ROHR is usually predicted based on the know pressure trace [8, 17, 23]. To avoid iterations, n-zone modeling framework, as for example given in [32, 33], was rearranged in a way to be an explicit function of derivative of the mass of burned fuel while considering only two zones. Governing equations are derived for a high pressure phase only while considering mass and enthalpy variation due to the injected fuel and enthalpy variation due to the heat transfer. In the analyzed case, the blow-by was not considered, whereas modeling framework can easily be extended to consider this effect. Conservation of mass in the particular zone is expressed as :     $$\frac{d{m}_B}{d\phi }=\frac{d{m}_{U\to B}}{d\phi }\text{,}$$(1) $$\frac{d{m}_U}{d\phi }=\frac{d{m}_{U,\text{inj}}}{d\phi }-\frac{d{m}_{U\to B}}{d\phi }\text{,}$$(2) where the equation for mass transport from the unburned to the burned zone is derived on the basis of equation reported in Ref. [16], considering change of fuel mass (dm_fuel) from Equation (8) and relative air-fuel ratio of the unburned zone (λ_U) from the previous calculation step as: $$\frac{d{m}_{U\to B}}{d\varphi }=\frac{d{m}_{\text{fuel}}}{d\varphi }\cdot \frac{\left({\lambda }_{\text{comb}}·L_{\text{st}}+1\right)\cdot \left(1+{\lambda }_U\cdot L_{\text{st}}\right)}{L_{\text{st}}\cdot \left({\lambda }_U-1\right)}\text{.}$$(3)

The conservation of total in-cylinder mass is expressed as: $$\frac{dm}{d\phi }=\frac{d{m}_B}{d\phi }+\frac{d{m}_U}{d\phi }\text{.}$$(4)

Conservation of energy in the particular zone is expressed as: $$\frac{d{U}_B}{d\phi }=\frac{\mathrm{d}{Q}_{B,B}}{d\phi }+\frac{\mathrm{d}{Q}_{\text{ht},B}}{d\phi }+\frac{\mathrm{d}{H}_B}{d\phi }-p\cdot \frac{\mathrm{d}{V}_B}{d\phi }\text{,}$$(5) $$\frac{d{U}_U}{d\phi }=\frac{\mathrm{d}{Q}_{\text{ht,}U}}{d\phi }+\frac{\mathrm{d}{Q}_{\text{evp,}U}}{d\phi }+\frac{\mathrm{d}{H}_U}{d\phi }-p\cdot \frac{\mathrm{d}{V}_U}{d\phi }\text{,}$$(6) where dH_B, dH_U represent change of enthalpy and dV_B, dV_U volume change for burned and unburned zone. Equations were derived under the assumption that combustion only takes place in the burned zone yielding the following expression for the ROHR: $$\frac{\mathrm{d}{Q}_{B,B}}{d\phi }=\frac{d{m}_{\text{fuel}}}{d\phi }\cdot H_{\text{LHV}}\text{,}$$(7)where H_LHV represents lower heating value of the fuel. After rearranging the equations for mass (1–4) and energy (5–6) the equation for the derivative of the mass of burned fuel reads:

See Equation (8) below: $$X=\frac{1}{\left(V-m_B\cdot \left(\frac{A_B\cdot T_B\cdot R_B}{p}\right)-m_U\cdot \left(\frac{A_U\cdot T_U\cdot R_U}{p}\right)\right)}\text{,}$$ $$Y=\frac{\left({\lambda }_{\text{comb}}·L_{\text{st}}+1\right)\cdot \left(1+{\lambda }_U\cdot L_{\text{st}}\right)}{L_{\text{st}}\cdot \left({\lambda }_U-1\right)}\text{,}$$ $$A_B=\frac{R_B+T_B\cdot \frac{\partial R_B}{\partial T_B}}{R_B+T_B\cdot \frac{\partial R_B}{\partial T_B}+\frac{\partial u_B}{\partial T_B}}\text{,}$$ $$A_U=\frac{R_U+T_U\cdot \frac{\partial R_U}{\partial T_U}}{R_U+T_U\cdot \frac{\partial R_U}{\partial T_U}+\frac{\partial u_U}{\partial T_U}},$$ $$\frac{d{m}_{\text{fuel}}}{d\phi }=\left(-\frac{dp}{d\phi }-X·p\frac{dV}{d\phi }+X·A_B·\frac{d{Q}_{\text{ht},B}}{d\phi }+X·A_U·\frac{d{Q}_{\text{ht},U}}{d\phi }+X·A_U·h_{\text{inj}}·\frac{d{m}_{U,\text{inj}}}{d\phi }+X·A_U·q_{\text{evp}}·\frac{d{m}_{U,\text{inj}}}{d\phi }+X·A_U·\left(\frac{\frac{\partial u_U}{\partial T_U}{T}_U}{1+\frac{\partial R_U}{\partial T_U}\frac{T_U}{R_U}}-u_U\right)\frac{d{m}_{U,\text{inj}}}{d\phi }\right)/\left[-X·A_B·H_{\text{LHV}}+Y·X·A_B·h_B-Y·X·A_U·h_U+Y·X·A_B·\left(\frac{\frac{\partial u_B}{\partial T_B}{T}_B}{1+\frac{\partial R_B}{\partial T_B}\frac{T_B}{R_B}}-u_B\right)-Y·X·A_U·\left(\frac{\frac{\partial u_U}{\partial T_U}{T}_U}{1+\frac{\partial R_U}{\partial T_U}\frac{T_U}{R_U}}-u_U\right)+2Y·X·A_U·\left(\frac{\left(\frac{\partial u_U}{\partial T_U}{T}_U\frac{\partial R_U}{\partial {\lambda }_U}\right)}{\left(R_U+\frac{\partial R_U}{\partial T_U}{T}_U\right)}-\frac{\partial u_U}{\partial {\lambda }_U}\right)\frac{1}{(1+{\lambda }_U{L}_{\text{st}})}-X·A_B·\left(\frac{\left(\frac{\partial u_B}{\partial T_B}{T}_B\frac{\partial R_B}{\partial {\lambda }_B}\right)}{\left(R_B+\frac{\partial R_B}{\partial T_B}{T}_B\right)}-\frac{\partial u_B}{\partial {\lambda }_B}\right)\left(1+Y·\left(\frac{1}{(1+{\lambda }_U{L}_{\text{st}})}-\frac{1}{(1+{\lambda }_B{L}_{\text{st}})}\right)\right)\right],$$(8)

and equations for temperature derivative of the burned and unburned zone read $$\frac{d{T}_B}{d\varphi }=\left[\frac{d{Q}_{B,B}}{d\varphi }+\frac{d{Q}_{\text{ht},B}}{d\varphi }+h_U\cdot \frac{d{m}_B}{d\varphi }-h_B\cdot \frac{d{m}_B}{d\varphi }+\frac{m_B\cdot R_B\cdot T_B}{p}\cdot \frac{\mathrm{d}p}{d\varphi }-m_B\cdot \left(\frac{\partial u_B}{\partial {\lambda }_B}+T_B\cdot \frac{\partial R_B}{\partial {\lambda }_B}\right)\cdot \frac{d{\lambda }_B}{d\varphi }\right]/\left(m_B\cdot \left(\frac{\partial u_B}{\partial T_B}+\left(1+\frac{T_B}{R_B}\cdot \frac{\partial R_B}{\partial T_B}\right)\cdot R_B\right)\right)\text{,}$$(9) and $$\frac{d{T}_U}{d\phi }=\left[\frac{d{Q}_{\text{ht},U}}{d\phi }+\frac{d{Q}_{\text{evp},U}}{d\phi }+h_{\text{inj}}\cdot \frac{d{m}_{U,\text{inj}}}{d\phi }+\frac{m_U\cdot R_U\cdot T_U}{p}\cdot \frac{\mathrm{d}p}{d\phi }-m_U\cdot \left(\frac{\partial u_U}{\partial {\lambda }_U}+T_U\cdot \frac{\partial R_U}{\partial {\lambda }_U}\right)\cdot \frac{d{\lambda }_U}{d\phi }\right]/\left(m_U\cdot \left(\frac{\partial u_U}{\partial T_U}+\left(1+\frac{T_U}{R_U}\cdot \frac{\partial R_U}{\partial T_U}\right)\cdot R_U\right)\right),$$(10)

m_UBZ represents the sum of the mass of the trapped air, recirculated exhaust gases and residual gases and thus determines the oxygen concentration of the unburned zone. From these three contributions it is most demanding to determine the amount of residual gases. The values of the residual gases were determined in the off-line analyses and stored in the look-up tables for on-line applications, which resembles calibration process of control functionalities. Both tested engines were running without the exhaust gas recirculation (EGR) and thus the EGR rate was set to zero. Otherwise this parameter is available in the ECU, while, if required, modeling framework proposed in this section can also be used to determine the EGR rate. As the evaporation of fuel, which occurs only in the unburned zone (Eq. 8), has minor impact on the burned zone temperature, it was assumed that there is no delay between injection rate (dm_inj/dφ) and the rate of fuel burned. (EGR)

3.2 Heat transfer sub-model

The total heat transfer between the in-cylinder charge and the combustion chamber walls is estimated by means of convective formulation [34]: $$\frac{d{Q}_{\text{ht}}}{d\varphi }=\frac{1}{6n}\cdot {\alpha }_{\text{HT}}\cdot \left(A_p\cdot \left(T_p-T_{\text{cyl}}\right)+A_l\cdot \left(T_l-T_{\text{cyl}}\right)+A_h\cdot \left(T_h-T_{\text{cyl}}\right)\right)\text{,}$$(11) where A_p, A_l and A_h are surface areas (indexes represent piston, liner and cylinder head respectively), T_p, T_l and T_h are temperatures and α_HT is heat transfer coefficient. The mean temperature in the cylinder (T_cyl) is for the purpose of heat transfer sub-model obtained from the ideal gas law. In the analyzed case, the Woschni correlation [35, 36] was applied to determine the heat transfer coefficient: $${\alpha }_{\text{HT}}=130\cdot D^{-0.2}\cdot p^{0.8}\cdot T_{\text{cyl}}^{-0.53}\cdot {\left[\left(2.28+0.308\cdot \frac{v_u}{v_m}\right)\cdot v_m+0.00324\cdot \frac{V\cdot T_{\text{IVC}}}{p_{\text{IVC}}\cdot V_{\text{IVC}}}\cdot \left(p-{\mathit{p}}_{\text{IVC}}\right)\right]}^{0.8},$$(12)where D represents cylinder bore, p is in-cylinder pressure, V represents actual cylinder volume, v_u is circumferential velocity and v_m is mean piston speed. T_IVC, p_IVC and V_IVC are temperature, pressure and volume in the cylinder at intake valve closing (IVC).

3.3 Thermodynamic data

Thermodynamic properties such as R_B, R_U, u_B, u_U and their partial derivatives ∂R_B/∂ T_B, ∂R_U/∂ T_U, ∂R_B/∂ λ_B, ∂R_U/∂ λ_U, ∂u_B/∂ T_B, ∂u_U/∂ T_U, ∂u_B/∂ λ_B, ∂u_U/∂ λ_U were pre-calculated and stored in the look-up tables as a function of temperature (300–3500 K) and relative air-fuel ratio (1–1000). The specific enthalpies of both zones must also be evaluated in order to solve the Equation (8) and are expressed as: $$h_B=u_B+R_B\cdot T_B\text{,}$$(13) $$h_U=u_U+R_U\cdot T_U\text{.}$$

3.4 Determination of RAFC

The evolution of the RAFC gives the amount of air that is consumed while burning a unit of fuel. It thus determines the mass, species and enthalpy transport from the unburned to the burned zone. Therefore, it significantly influences temperature and composition of the burned zone and thus formation of NO_x. In contrast to other published approaches, where it is considered that combustion takes place at stoichiometric conditions or that the air-fuel ratio is determined in a way that the whole unburned zone is consumed during combustion [22], in the proposed method the RAFC is a physically motivated calibration parameter. As the amount of the burned fuel is not known during the combustion, the variation of the RAFC is not given as function of the relative amount of the burned fuel. Instead, it is given as a function of the crank angle degree to ensure general applicability and robustness of the model.

At the SOC point (Fig. 1), the RAFC is initialized and it is afterwards kept constant until the end of injection (EOI). In the analyzed cases, RAFC was initialized with the values between 1.0 and 1.1. In contrast to some approaches where air-fuel ratio in the flame front is assumed to be <1 [37], in this model RAFC was kept above 1.0 during the initial combustion phase. This is aligned with the main focus of the model being NO_x formation, as at least some surplus of oxygen in the burned zone is required to allow for NO_x formation during the initial and thus also premixed combustion phase, which is indeed characterized by NO_x formation [38]. After EOI, RAFC is linearly increasing thereby leaning burned gases. This trend is motivated by the fact that during that period spray detaches from the nozzle, which provokes more pronounced mixing with air as fuel occupies larger volume. The RAFC variation as presented in Figure 1 is given by Equation (14), $${\lambda }_{\text{com}}=\{\begin{array}{c}\hfill {\lambda }_{\text{com}\_\text{init}};\mathrm{i}\mathrm{f}\phi <\text{EOI}\hfill \\ \hfill k_{\text{lam}}\cdot \phi +n_{\text{lam}};\mathrm{i}\mathrm{f}\text{EOI}\le \phi <\text{EVO}\hfill \end{array},$$(14) $$n_{\text{lam}}={\lambda }_{\text{com}\_\text{init}}-\left(\mathrm{S}\mathrm{O}\mathrm{I}\mathrm{+}\mathrm{I}\mathrm{D}\right)\cdot k_{\text{lam}}\text{,}$$ where k_lam and n_lam are constants of the linear function and λ_{com_init} is initial value of RAFC.

3.5 NO_x sub-model

The calculation of species equilibrium fractions is usually performed by the time-consuming minimization of Gibbs free energy [39]. For minimization of computation effort of complex equilibrium calculations in the present work all species equilibrium fractions were given in pre-calculated look-up tables determined on the basis of JANAF data tables [26]. The equilibrium calculation is performed for a set of a six relevant species (N, N₂, O, O₂, NO, OH) as a function of temperature (1000–3500 K) and RAFC (1–4), which covers the relevant range for the burned zone. The formation of NO_x was modeled according to the extended Zeldovich mechanism described by the following six reactions [40]: $$\mathrm{R}\mathrm{1}\mathrm{:}{\mathrm{N}}_{\mathrm{2}}\mathrm{+}\mathrm{O}\mathrm{\leftrightarrow }\mathrm{N}\mathrm{O}\mathrm{+}\mathrm{N},$$(15) $$\mathrm{R}\mathrm{2}\mathrm{:}{\mathrm{O}}_{\mathrm{2}}\mathrm{+}\mathrm{N}\mathrm{\leftrightarrow }\mathrm{N}\mathrm{O}\mathrm{+}\mathrm{O}\mathrm{\text{,}}$$ $$\mathrm{R}\mathrm{3}\mathrm{:}\mathrm{O}\mathrm{H}\mathrm{+}\mathrm{N}\mathrm{\leftrightarrow }\mathrm{N}\mathrm{O}\mathrm{+}\mathrm{H}\mathrm{\text{,}}$$ $$\mathrm{R}\mathrm{4}\mathrm{:}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\mathrm{+}\mathrm{O}\mathrm{\leftrightarrow }\mathrm{N}\mathrm{O}\mathrm{+}\mathrm{N}\mathrm{O}\mathrm{\text{,}}$$ $$\mathrm{R}\mathrm{5}\mathrm{:}{\mathrm{O}}_{\mathrm{2}}\mathrm{+}{\mathrm{N}}_{\mathrm{2}}\mathrm{\leftrightarrow }{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\mathrm{+}\mathrm{O}\mathrm{\text{,}}$$ $$\mathrm{R}\mathrm{6}\mathrm{:}\mathrm{O}\mathrm{H}\mathrm{+}\mathrm{O}\mathrm{\leftrightarrow }\mathrm{N}\mathrm{O}\mathrm{+}\mathrm{N}\mathrm{\text{.}}$$

Even though the NO₂ formation in compression ignition engines cannot be neglected [41], in the described model approximation on NO_x formation was made only on calculation of NO formation. Similar as in, e.g. [8, 42], NO formation and NO_x emissions are assumed to be correlated and one constant calibrating factor was applied on the final cycle NO_x mass to match the experimental tailpipe NO_x emissions. The rate of formation of NO in the burned zone is therefore derived by the following equation [11, 16, 40]: $$\frac{d{c}_{\text{NO}}}{d\phi }=\frac{2}{6n}\cdot C_{\text{NO,}\text{multi}}\cdot \left(1-{{\alpha }_{\text{NO}}}^2\right)\cdot \left[\frac{k_1\cdot c_{e,{\mathrm{N}}_2}\cdot c_{e,O}}{1+{\alpha }_{\text{NO}}\cdot K_2}+\frac{k_4\cdot c_{e,{\mathrm{N}}_2\mathrm{O}}\cdot c_{e,O}}{1+K_4}\right],$$(16) where α_NO is the ratio of actual NO concentration and its equilibrium concentration calculated according to Equation (17): $${\alpha }_{\text{NO}}=\frac{c_{\text{NO}}}{c_{e,\text{NO}}},$$(17)

_n is engine speed and C_NO,multi is multiplicative parameter that also takes into account difference between predicted NO formation and actual NO_x formation. Furthermore, c_e represents equilibrium species concentration read from the look-up table. K₂ and K₄ are help variables calculated according to Equation (18): $$K_2=\frac{k_1\cdot c_{e,N_2}\cdot c_{e,O}}{k_2\cdot c_{e,O_2}\cdot c_{e,N}+k_3\cdot c_{e,\text{OH}}\cdot c_{e,N}}\text{,}$$(18) $$K_4=\frac{k_4\cdot c_{e,N_{2}O}\cdot c_{e,O}}{k_5\cdot c_{e,O_2}\cdot c_{e,N_2}+k_6\cdot c_{e,OH}\cdot c_{e,N_2}}\text{.}$$

Variables denoted with k are determined by the following equations [16]: $$k_1=4.93\cdot {10}^{13}\cdot T_B^{0.047}\cdot e^{\left(\frac{-38048.01}{T_B}\right)}\text{,}$$(19) $$k_2=1.48\cdot {10}^8\cdot T_B^{1.5}\cdot e^{\left(\frac{-2859.01}{T_B}\right)}\text{,}$$ $$k_3=4.22\cdot {10}^{13}\text{,}$$ $$k_4=4.58\cdot {10}^{13}\cdot e^{\left(\frac{-12130.6}{T_B}\right)}\text{,}$$ $$k_5=2.25\cdot {10}^{10}\cdot T_B^{0.825}\cdot e^{\left(\frac{-50569.7}{T_B}\right)}\text{,}$$ $$k_6=9.14\cdot {10}^{07}\cdot T_B^{1.148}\cdot e^{\left(\frac{-36190.66}{T_B}\right)}$$

The adjustment of the model to the measurements was done with the scaling parameter C_NO,multi, which was calibrated in the model calibration phase and was then kept constant for the entire operating range of an engine.

4 Results

4.1 Impact of code execution source

The applied modeling framework of NO_x model was developed and calibrated in Lab View environment on a PC, while execution of the code was performed on the FPGA integrated circuit. Due to the different architecture of the chips, the code implemented on the FPGA integrated circuit was restructured in a way to avoid all nonlinear dependencies. Furthermore, on the FPGA integrated circuit it is necessary to specify word length and integer bites for every mathematical operator in order to improve computational times and reduce memory consumption. Therefore, it is important to verify consistency of the results calculated on a PC and on a FPGA integrated circuit. This was done by using the same pressure trace with all the related thermodynamic data as an input for models on both hardware platforms and by comparing the result of calculated temperature in the burned zone and the total NO_x mass. As discernable from the results shown in Figure 2, it can be concluded that the difference in the result is within the line thickness and the error between the execution sources is less than 1% regarding total NO_x mass. These results thus confirm consistency of the models on both hardware platforms.

4.2 Impact of the gas property treatment

Before analyzing NO_x emissions, it is meaningful to assess plausibility of the thermodynamic modeling framework and gas property treatment presented in section 3. As some of the RT two-zone NO_x models are using constant gas properties in the calculation of the thermodynamic parameters [18, 19], in this analysis importance to consider variable gas properties to establish a sound basis for modeling NO_x emissions is demonstrated. To demonstrate this, two different cases were analyzed (Fig. 3): 1. the modeling depth as considered in Equations (8)–(10), i.e. variable gas properties of R, u and their derivatives (denoted with: T_base) and 2. a simplified approach with constant values of R, u (denoted with T_u const_R const). The impact of gas property treatment is demonstrated on the case of temperature in the burned zone (Fig. 3) being one of the most influential parameters of NO_x kinetics. It can be seen in Figure 3 that the simplified gas property treatment results in unrealistically high temperatures. This can mainly be attributed to neglecting the dissociation. Consequently, simplified gas property treatment yields too low values of specific heats at high temperatures, whereas neglecting associated partial derivatives additionally increases the error. These results indicate that variable gas properties need to be considered to establish a sound basis for modeling NO_x emissions.

4.3 Thermodynamic model results

In this section, the traces of thermodynamic parameters and NO_x mass are shown in Figures 4 and 5 for the two operating points of the LD engine at 3000 min⁻¹ and at 20 Nm and 160 Nm to provide an insight into underlying phenomena and to form the basis for further analyses where only data at exhaust valve opening (EVO) will be shown. Figure 4 shows measured in-cylinder pressure trace being an input to the model (Fig. 4a)), ROHR (Fig. 4b) calculated by Equations (7) and (8), masses of the burned and unburned zone (Fig. 4c) calculated by Equations (1)–(3), the RAFC (Fig. 4d) given by Equation (14), temperatures of the burned zone, unburned zone and mean in-cylinder temperature (Fig. 4e) calculated by Equations (9) and (10) as well as actual NO_x mass calculated by Equation (16) compared to the equilibrium NO mass at instant temperature and composition in burned zone (Fig. 4f) at 3000 min⁻¹ and 20 Nm.

Figure 4b clearly indicates combustion of fuel injected in the pilot and the main injection, which is a direct consequence of the shape of the pressure trace in Figure 4a. ROHR (Fig. 4b) is also reflected in mass trace of the unburned to burned zone (Fig. 4c), where there is first a mass transfer due to combustion of a pilot injection followed by a short period of no mass transfer and then a more pronounced mass transfer due to combustion of the fuel injected during the main injection. As the RAFC (Fig. 4d) increases after the EOI, the rate of mass transfer from the unburned to the burned zone per mass of burned fuel increases towards the EOC. It can be seen in Figure 4c that a large amount of mass is not transported from the unburned to the burned zone despite increasing RAFC during the end of combustion. This can be attributed to high relative air-fuel ratio in the combustion chamber due to small mass of injected fuel during this low load operating point. Therefore, only a part of the in-cylinder charge is mixed with the fuel and thus being subjected to combustion at low loads. ROHR (Fig. 4b) trace is also reflected in the traces of mean in-cylinder temperature (Fig. 4e), whereas it also dictates the shape of burned zone temperature and NO_x traces (Fig. 4f).

Within the proposed modeling framework, equilibrium NO_x mass and NO_x concentration in the burned zone feature increased values during the combustion of the fuel injected during pilot injection. Afterwards, there is a slight drop in both traces due to no heat addition by combustion and presence of heat extraction through the heat transfer to the combustion chamber walls. During the combustion of the fuel injected during main injection there is first a slight increase in equilibrium concentration, which is followed by its drop as the temperature of the burned zone is decreasing as the piston moves towards the BDC. Unlike equilibrium concertation in the burned zone, equilibrium NO_x mass is increasing until approx. 17 deg. CA, as mass and volume of the burned zone are still increasing. As NO kinetics according to Zeldovich mechanism is slow, actual NO_x values significantly deviates from the equilibrium ones. It can thus be seen that within the proposed modeling framework actual NO_x concentration in the burned zone peaks in the decreasing phase of the equilibrium concentration as at that instant there are favorable conditions for NO_x formation. During the combustion of the fuel injected in the main injection actual NO_x concentration in the burned zone is thus decreasing, however actual NO_x mass still increases, which is again related to increasing mass of the burned zone. Although, NO_x mass can in the present study be experimentally validated only as engine-out emissions, i.e. at the EVC, the predicted trend qualitatively coincides well with the published results of in-cylinder NO_x traces, e.g. [38, 41], which also indicate non-negligible NO_x formation during the premixed burning period. This confirms meaningfulness of keeping RAFC above value 1.0 during the initial combustion phase, which ensures surplus of oxygen in the burned zone being required for NO_x formation.

In Figure 5, the same thermodynamic properties are shown for operating point with higher engine load. Consequently, this operating point is characterized by higher in-cylinder pressure (Fig. 5a) and higher ROHR values (Fig. 5b) corresponding to higher injected mass of fuel. Therefore, higher portion of mass of unburned zone is transferred to burned zone. In comparison with Figure 4c where less than 50% mass of unburned zone was used during combustion, at this higher load around 80% of mass of unburned zone took part in the combustion process within the proposed modeling framework. This is due to low relative air-fuel ratio in the combustion chamber and also due to more intense combustion in the later combustion phases, which are characterized by larger values of the RAFC (Fig. 5d). Relatively long period of high burned zone temperature (Fig. 5e) is the main reason for approximately 8 times higher NO_x mass (Fig. 5f) in this high load operating point compared to the low load point analyzed in Figure 4.

4.4 Model validation

As presented in the Introduction, RAFC represents the main calibration parameter of the virtual NO_x sensor. To analyze the influence of this calibration parameter two approaches were analyzed: 1. RAFC function (Eq. (13) featuring a single set of parameters was used for the entire engine operating range denoted as “const_PAR”, and 2. parameters of the RAFC function (Eq. (13) are operating point dependent denoted as “var_PAR”.

4.4.1 LD engine

In the model calibration phase, sensitivity analysis on the parameters of the RAFC function Equation (13) was performed. Based on this knowledge, parameters of the RAFC function were determined manually as this study is aimed to demonstrate capabilities of the modeling framework and thus automated model parametrization was not the scope of the research.

For the case with constant RAFC parameters, initial value of RAFC was thus set to 1.01 and constant of linear function k_lam was set to 0.02 (Tab. 5). Good agreement of the NO_x model results with the measurement even without calibration (Tab. 5) was demonstrated in the low and middle load operating points (e.g. 20 Nm and 100 Nm) at all engine speeds with the most pronounced deviation in the operating point 1200_20. In this operating point, constant RAFC parameters lead to higher burned zone temperature compared to the case with variable RAFC (Tab. 6) thus yielding too high NO_x mass. At high engine load (e.g. 160 Nm) deviation between the simulated and measured values is slightly higher. Compared to the case with variable RAFC (Tab. 6) this can be attributed to the more intense transport of mass from unburned zone to burned zone in the late combustion phase as a result of a relatively high k_lam parameter.

This case represents a very limiting case, where entire engine operating range is covered by a single set of model parameters of the emission formation model within a mechanistically based modeling framework, which is to the best of author's knowledge unique in this modeling depth. By additionally considering the fact that extremely complex physicochemical phenomena in diesel engines are treated with a very simplistic approach, it can be concluded that the modeling framework provides satisfactory results of the NO_x emissions.

As most of the studies related to virtual NO_x sensors use operating point specific calibration parameters [17, 43], an operating point specific input of the RAFC function was also applied in this study. As presented in the Introduction it is expected that RAFC, as a physically motivated calibration parameter, allows for high NO_x emission prediction capability while using only RAFC as a calibration parameter. In Table 6, the results for the LD engine are shown for this case. It can be seen from Table 6 that only relatively small spread of parameters of the RAFC function was required for the parameter λ_{com_init} and larger spread is characteristic for the parameter k_lam. At certain engine speed the highest value of k_lam was set at lower load and was decreased with increasing load. That is in agreement with the change of the relative air-fuel ratio in the unburned zone and confirms plausibility of the applied approach. Additional graphical interpretation of the NO_x results presented in Tables 5 and 6 is given in Appendix A.

4.4.2 HD engine

The main difference in the hardware configuration between the LD and the HD engine is in fuel injection systems influencing the preparation of the air and fuel mixture. Fuel injection pressure is lower in the HD engine and HD engine does not feature a pilot injection. In addition, HD engine generally operates with leaner mixture than the LD engine. These differences have impact on the concentration and the temperature at the SOC and also later during combustion phase. This leads to different engine specific set of calibration parameters used for the case with constant RAFC parameters. As can be seen from Table 7, for the case with constant RAFC parameters, parameters of the RAFC function feature larger values for the HD engine compared to the LD one (Tab. 5). Initial value of RAFC was set to 1.1 and constant of linear function k_lam was set to 0.04.

On the HD engine, differences between predicted and measured values of mass of the NO_x in the case with constant parameters are even smaller (Tab. 7) than for the LD engine. One of possible explanations for this trend is related to the absence of the pilot injection in the HD engine. Current approach, due to the sake of generality and simplicity, assumes equal RAFC values for all injections, which is in general not the case in real engines. As dwell time between injections varies this simplification impacts accuracy of the NO_x prediction for the engines with multiple injections as in the analyzed case presented for the LD engine.

At lower engine speed, the matching of the predicted values with measurements is very good at low to middle load where even without additional tuning results were within the desired accuracy (Tab. 7). At high engine load were predicted values of the NO_x mass quite higher than the measured in the case with constant parameters which is similar pattern that was already observed at LD engine at middle and high engine speed where the RAFC at late combustion period were decreased to match the measured values. In this case, where the RAFC is already higher due to initialization value it should be further increased in order to increase mass transport from unburned zone to burned zone and consequently decrease temperature in burned zone. This difference between the engine specific tuning procedures is the result of much higher relative air-fuel ratio in combustion chamber in HD engine. At higher engine speed a constant shift in the predicted values compared to measured data can be observed which can be attributed to selected initial value of RAFC and was effectively solved by increasing initial RAFC value as can be seen from Table 8. Additional graphical interpretation of the NO_x results presented in Tables 7 and 8 is given in Appendix B.

4.5 Computational times and memory consumption

The aim of this work was to develop a model of such a format and with sufficiently low demands of computational power that it could be calculated on an embedded system and is potentially compliant with next generation large series commercial hardware. The use of look-up tables was crucial to obtain a fast and computationally efficient algorithm since the FPGA architecture does not allow calculation of nonlinear dependencies. Despite that, commercial engine ECUs are also using look-up tables to model nonlinear and operating point dependent behaviors [20].

Precision for every mathematical operating on the FPGA was manually determined in order to improve computational times and reduce memory consumption. Total flash memory used for implementation of the virtual NO_x sensor on FPGA was 2.34 MBs (14.6% of total flash memory available on FPGA), which proves that the modeling framework could be implemented on a commercial ECU of which the most powerful version has 8 MB of flash memory [44].

Calculation time for one cycle was determined by measuring time for calculation of 100 consecutive cycles. The calculation time of one engine operating cycle is 21.8 ms, which means that the FPGA could perform the complete thermodynamic and NO_x calculation of one engine cycle up to engine speed at 5500 min⁻¹, which is higher than the maximum speed of most of the modern automotive diesel engines. Such computational time thus ensures ample margin to adapt the control parameters in the subsequent cycle. Since the commercial ECUs consist of a microcontroller with integrated flash memory and Application Specific Integrated Circuit (ASIC) chip architecture [45] and it is generally excepted that the FPGA is not as fast as ASIC [23], it is therefore concluded that the calculation could also be executed on an commercial ECU with similar or even shorter computational times.

5 Conclusion

A real-time capable virtual NO_x sensor for diesel engines was presented in the paper. The main characteristics of the modeling framework can be summarized as: (1) the thermodynamic modeling framework was developed in a way that the thermodynamic equations could be solved in a closed form while applying a detailed two-zone modeling framework and variable gas properties, (2) the NO_x sub-model relies only on one innovative physically motivated calibration parameter being relative air-fuel ratio during combustion (RAFC). Both features are crucial to effectively combine computational efficiency and sufficient modeling depth to achieve high level of predictiveness. Besides computational efficiency and high level of predictiveness, an important advantage of the proposed modeling approach is also relatively low calibration effort of the RAFC function. It was shown in the paper that the model features good prediction capability if only unique RAFC input is used for the entire engine operating regime. It was also shown that prediction capability can be further improved with additional operating point specific calibration of the RAFC input. This parameter can be later stored in the look-up table and programed on the ECU. It was shown that demands on computational resources are within the specifications of the current ECUs. High prediction capability of the virtual NO_x sensor was demonstrated on two significantly different engines featuring also two different fuel injection systems, which proves general applicability of the method.

Abbreviations

A: surface area [m²]

BDC: bottom dead center

CA: crank angle [°]

CFD: computational fluid dynamics

CLD: chemiluminescence detector

CNO, multi: NO scaling parameter

D: cylinder bore [m]

ECU: engine control unit

EGR: exhaust gas recirculation

EOI: end of injection

EVO: exhaust valve opening

FIR: finite impulse response

FPGA: field programmable gate array

H: enthalpy [J]

HD: heavy duty

H_LHV: lower heating value [J/kg]

LD: light duty

L_st: stoichiometric ratio

MFB: mass fraction burned

N: nitrogen

N₂O: nitrous oxide

NO: nitric monoxide

NO₂: nitric dioxide

NOx: nitric oxides

O: oxygen

OH: hydroxide

PSG: glow plug pressure transducer

Q: heat energy [J]

R: gas constant [J/(kg K)]

RAFC: relative air-fuel ratio during combustion

ROHR: rate of heat release

RT: real time

SOC: start of combustion

SOI: start of injection

T: temperature [K]

TDC: top dead center

U: internal energy [J]

V: cylinder volume [m³]

c: species concentration [mol/cm³]

c_e: equilibrium species concentration [mol/cm³]

h: specific enthalpy [J/kg]

k_lam: constant of the linear function of RAFC

m: mass [kg]

n: engine speed [min⁻¹]

n_lam: constant of the linear function of RAFC

p: pressure [Pa]

q_evp: latent heat [J/kg]

u: internal energy [J/kg]

v_m: mean piston speed [m/s]

v_u: circumferential velocity [m/s]

α_HT: heat transfer coefficient

α_NO: ratio of actual NO concentration and its equilibrium concentration

λ: relative air-fuel ratio

${\lambda }_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{\_}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}$: initial value of RAFC

φ: angle [°]

Acknowledgments

The research is co-funded by the Republic of Slovenia and the European Union from the European Regional Development Fund within the framework of the EVA4GREEN program. The authors also acknowledge the financial support from the Slovenian Research Agency (research core funding No. P2-0401 and young researcher program with evidence No. 35395).

Appendix A

Appendix B

References

All Tables

All Figures