A Dynamic Simulation Strategy for PCCI Combustion Control Design

A Dynamic Simulation Strategy for PCCI Combustion Control Design — Subject of this work is a dynamic simulation strategy for PCCI combustion that can be used in closed-loop control development. A detailed multi-zone chemistry model for the high-pressure part of the engine cycle is extended by a mean value model accounting for the gas exchange losses. The resulting stationary model is capable of describing PCCI combustion sufficiently well. It is at the same time very economic with respect to computational costs. The model is further extended by identified system dynamics influencing the stationary inputs. For this, a Wiener model is set up that uses the stationary model as a nonlinear system representation. In this way, a dynamic nonlinear model for the representation of the controlled plant Diesel engine is created. This paper summarizes an important outcome of the the Collaborative Research Centre "SFB 686 Modellbasierte Regelung der homogenisierten NiedertemperaturVerbrennung" at RWTH Aachen University and Bielefeld University, Germany. 550 Oil & Gas Science and Technology – Rev. IFP Energies nouvelles, Vol. 66 (2011), No. 4 Definitions, Acronyms, Abbreviations aTDC After Top Dead Center CA Crank Angle CA50 Crank Angle of 50% burnt mass ECU Electronic Control Unit EGR Exhaust Gas Recirculated FMI Total Fuel Mass Injected IMEP Indicated Mean Effective Pressure IMEPGE Indicated Mean Effective Pressure of the Gas Exchange IMEPHP Indicated Mean Effective Pressure of the High-Pressure cycle pae Pressure in the exhaust manifold (after engine) pbe Pressure in the intake manifold (before engine) PCCI Premixed Charge Compression Ignition rpm Revolutions per minute RWTH Rheinisch Westfälische Technische Hochschule Aachen SOI Start of Injection Tae Temperature in the exhaust manifold Tbe Temperature in the intake manifold TS Timing Sweep VGT Variable Geometry Turbine X0D Zero-dimensional multi-zone chemistry solver


INTRODUCTION
There is a strong demand to explore new combustion concepts capable of meeting the stringent future emission standards, such as Tier 2 Bin 5 on the North American market or EURO VI on the European market.Premixed Charge Compression Ignition (PCCI), or Premixed Compression Ignition (PCI), as it is sometimes referred to, has emerged as an interesting alternative to conventional Diesel combustion in the part-load operating range [17,24,41].This concept involves relatively early injection timings, high external EGR rates, and cooled intake air, leading to a Low-Temperature Combustion (LTC) process with low NO x and particulate emissions.A general review of LTC processes can be found in Dec [8].However, in contrast to conventional Diesel combustion, with early direct-injection LTC it can be difficult to prevent combustion from occurring before top dead center, which often increases noise and reduces engine efficiency.Sophisticated closed-loop control of the combustion process is one means to overcome this difficulty.Against this background, this paper deals with a first step towards the development of a controller for the PCCI combustion process.
The most important characteristics for the operation of an internal combustion engine are the engine's load and its combustion efficiency.The former is directly dependent on the Indicated Mean Effective Pressure (IMEP), the latter can be characterized by the center of combustion, where 50% of the total injected fuel mass is burned (CA50).The engine's load is set by the driver in an automobile application.For CA50, a set point can be determined which is dependent on the operating point and gives the best combustion efficiency.
In the recent past, several efforts have been reported in the literature that aim at controlling engine combustion.Among these, [5,20,21,26,28,39,42] can be mentioned.The standard procedure for creating a controller includes the modeling part as the first step [1,6,27,38,40].Often these models differ in several aspects from models which include details of the combustion process, like three-dimensional Computational Fluid Dynamics (CFD) models [4,7,32,33].From the viewpoint of automatic control, the dynamics describing the dependency of the system's outputs (IMEP, CA50) on the actors (SOI, external EGR rate, and total fuel mass injected) is of highest priority.Another requirement is an acceptable calculation speed, as these models are usually applied in dynamic closed-loop simulations.This paper presents a new approach to the development of a simulation model for the use in closed-loop control development.The model is based on a recently introduced multizone model for PCCI combustion that was systematically derived from a detailed three-dimensional CFD model.This model is extended by a physically inspired description of the gas exchange part of the engine cycle.For the use in closed-loop simulations the system's dynamics have to be covered.For this reason, the stationary model is further extended by identified system dynamics influencing the stationary inputs.In this manner, a stationary exact model is extended to a Wiener-type model with a static part describing the nonlinearities and an upstream part describing the system's dynamics.This novel procedure integrates the detailed knowledge from combustion simulation tools into closed-loop control and establishes a broad field of possibilities for testing completely new controlled process variables.

COMBUSTION MODEL FORMULATION
Crucial for reacting turbulent flows is the modeling of the chemistry.Here, a multi-zone chemistry model is employed.It covers the nonlinear dependencies within the highpressure part of the engine cycle with stationary exactness.

Multi-Zone Model Derivation from a Detailed CFD Model
Recently, Barths et al. [3] developed a consistent mixing model for the two-way-coupling of CFD codes and zerodimensional multi-zone codes, yielding a hybrid CFD/zonal modeling approach for HCCI chemistry that aims at reducing the computational overhead for large chemical mechanisms.There, all aspects except chemical kinetics were computed using three-dimensional CFD, and CFD cells were binned into a relatively small number of zones for the chemistry calculations.Further on, each CFD cell was binned into one of 27 zones based on its fuel mixture fraction, dilution, and total enthalpy.Issues within the two-way-coupling included maintaining thermodynamic consistency between the two representations and interactions among zones.Following from there, an interactively coupled CFD-multi-zone approach was presented in Felsch et al. [12] that can be used to simulate HCCI combustion.Later, Felsch et al. [13] applied this interactively coupled CFD-multi-zone approach for various PCCI operating conditions in a 1.9l FIAT GM Diesel engine.In the same work, the interactively coupled CFD-multi-zone approach was systematically reduced towards a computationally efficient stand-alone multi-zone model.For this, the multi-zone model was decoupled from the CFD code.In the framework of the interactively coupled CFD-multi-zone approach, all relevant information fed from the CFD code to the multi-zone model was identified.This information was then modeled within the decoupled standalone multi-zone model.The validity of the reduction was limited to the investigated 1.9l FIAT GM Diesel engine that is also subject of the present work.The stand-alone multizone model showed a significant advantage with respect to computational costs.Due to its strict derivation from the interactively coupled CFD-multi-zone approach, the standalone multi-zone model was at the same time nearly as predictive as the detailed CFD approach.
The theoretical background of the stand-alone multi-zone model is briefly summarized in the following.For detailed information regarding the reduction methodology from the interactively coupled CFD-multi-zone approach towards the stand-alone multi-zone model, the reader is referred to Felsch et al. [13].

Multi-Zone Model
The stand-alone multi-zone model employed in the current study is X0D, a zero-dimensional chemistry solver based on multiple zero-dimensional reactors.X0D was developed internally at General Motors R&D by Hardo Barths, Tom Sloane, and Christian Hasse, and was first described in Hergart et al. [18].

Governing Equations
The equations governing species mass fractions, temperature, and pressure change in the multi-zone chemistry code are given below: and In Equation (1), Y i j denotes the mass fraction of species j in zone i, and ωij is the corresponding chemical source term.ρs i j accounts for the source term due to fuel vaporization.It is zero for all species except the fuel itself (ρ s i j = 0, j fuel).The second term on the left-hand side of Equation (1) describes the mass exchange between the zones, where ṁik is the rate at which mass is transported between zones i and k.Thereby, nz stands for the total number of zones.Similarly, in Equation (2) the second and third term on the left-hand side represent the enthalpy exchange between zones due to enthalpy stratification between zones and due to species stratification between zones, respectively.Δh v j denotes the latent heat of vaporization of species j, and Qwall,i is the wall heat transfer of zone i.In Equations ( 2) and (3), nsp accounts for the number of species employed in the underlying chemical mechanism.The equation of state is used to derive Equation (3), which solves for the pressure across the zones.Through this equation, all zones are thermodynamically coupled with each other.The mixing in the multi-zone model is accounted for by allowing the different zones to exchange mass and energy with each other, in addition to the interaction through the pressure.They exchange their scalar quantities (species composition and enthalpy) based on the rate at which they exchange their mass.The mass of each zone is kept constant and the model can accommodate any number of zones.

Heat Transfer to the Walls
Wall heat transfer is described as where n wall denotes the total number of walls in the engine.A wall,l,i is the area of wall l belonging to zone i, h wall,l,i the heat transfer coefficient of zone i to wall l, T i the temperature of zone i, and T wall,l the temperature of wall l.The heat transfer coefficient h wall,l,i is modeled with the Woschni correlation [44], which reads where b is the bore of the engine, p the mean cylinder pressure, T av. the average in-cylinder temperature, and v piston the mean piston speed.The constant C 1 is adjustable and is set to 0.5.Since the average in-cylinder temperature T av. is used within the Woschni correlation and h wall,l,i does not depend on any specific wall, the heat transfer coefficient is equal for all zones, leading to h wall,l,i = h wall .
For the specific engine investigated, the area A wall,l,i of wall l belonging to zone i has to be taken from a preceding simulation with the interactively coupled CFD-multi-zone approach mentioned above.Then, this information can be incorporated into a successive simulation with the standalone multi-zone model.

Fuel Injection and Vaporization
A preceding simulation with the interactively coupled CFDmulti-zone approach also allows for determining the amount of fuel that is injected and vaporizes in each CFD zone.This information can then be incorporated into a successive simulation with the stand-alone multi-zone model (see source term ρs i j in Eq. ( 1) and ( 2)).

Mixing Model
The mixing model within the stand-alone multi-zone model assumes that the mass exchange rate ṁik of zone i with neighbouring zone k (see Eq. ( 1) and ( 2)) is proportional to the wall heat transfer coefficient h wall and the area A both zones share [18].This is an analogy to the well known empirical Lewis law, which assumes a similarity between heat and mass transfer [2].The heat transfer to the walls and the mass exchange between the zones are both induced by the turbulence intensity so that they can be related to each other.The corresponding equation reads Average parameter r MB as extracted from the interactively coupled CFD-multi-zone approach (denoted as AC-FluX-X0D) and as incorporated into the stand-alone multi-zone (denoted as X0D) model for one selected operating condition (SOI at −26 • CA aTDC, external EGR rate of 30%, injected fuel mass of 10.2 mg/cycle).
In Equation ( 6), c pi is the average heat capacity at constant pressure of zone i. r MB is a parameter that is adjusted to yield the desired level of mixing.The area A is computed from the volumes V i and V k of the two zones i and k according to In the following, some aspects of the systematic reduction procedure from the interactively coupled CFD-multizone approach towards the stand-alone multi-zone model are discussed.It is possible to evaluate the choice of the parameter r MB in Equation ( 6) by means of a preceeding simulation with the interactively coupled CFD-multi-zone approach (Details are given in Felsch et al. [13].) Figure 1 displays the average parameter r MB around top dead center firing for a selected operating condition (SOI at −26 • CA aTDC, external EGR rate of 30%, injected fuel mass of 10.2 mg/cycle), as extracted from a simulation with the interactively coupled CFD-multi-zone approach and as incorporated into the stand-alone multi-zone model.For incorporation, r MB was set to 75.0 during compression and expansion, with a linear increase up to 500.0 during injection and a linear decrease back again to 75.0 after end of injection.Qualitatively, both evolutions are similar.Altogether, r MB decreases after EOI with decaying turbulence intensity after EOI.The evolution of r MB as calculated from the interactively coupled CFDmulti-zone approach reveals two peaks.These two peaks are due to first and second stage ignition.The corresponding integral average of the whole simulation is equal 169.7.
The SOI variation shown in Figures 3, 4, 5, and 6 was also investigated in a previous work by Felsch et al. [14], with just one constant value of 100.0 being chosen for r MB , which is of the same order of magnitude as 169.7.In Felsch et al. [14], the overall simulation results for the stand-alone multi-zone model in terms of pressure curves, IMEP HP , CA10, CA50, CO emissions, combustion efficiency, and NO/NO x emissions were not as accurate in predicting the experimental results as the ones presented in the current work, where the evolution of r MB was qualitatively adjusted to the one from the three-dimensional approach.
It should be emphasized, too, that the order of magnitude of the parameter r MB strongly depends on the level of premixing.In Hergart et al. [18], where more advanced injection timings were investigated, a constant value of 6.0 was sufficient for r MB .Therefore, this modeling parameter must be adjusted for each engine type.

Chemistry Model
Combustion chemistry in X0D is described by a detailed chemical mechanism that comprises 59 elementary reactions among 38 chemical species.This mechanism describes low-temperature auto-ignition and combustion of n-heptane, which serves as a surrogate fuel for Diesel in this work.In addition, it accounts for thermal NO formation.The chemical mechanism for n-heptane was constructed by Peters et al. [31].The NO-submechanism, which is part of the full mechanism, is the extended Zeldovich mechanism (see for example Heywood [19]).

Computational Accuracy and Efficiency
The stand-alone multi-zone model does not solve for the solution of the flow field.This yields a significant advantage in computational costs in comparison to three-dimensional CFD modeling approaches.The run time for the stand-alone multi-zone model is approximately 3.3 minutes per highpressure engine cycle for a simulation with 15 zones on a Dell PowerEdge 1950 with two Woodcrest 5160 Intel Xeon Dual-Core CPUs with 3.0 GHz.(The stand-alone multizone model applies 15 zones in order to provide a sufficient discretization and to equally achieve computational efficiency for the calculation of auto-ignition processes.)For such a simulation, the stand-alone multi-zone model written in FORTRAN 77 was compiled with the Intel Fortran Compiler 9.1.043.However, modeling the CFD information according to the procedure outlined above may reduce the accuracy of the simulation outcome.A detailed discussion of the computational accuracy and efficiency of the standalone multi-zone model can be found in [13].

STATIONARY VALIDATION
At the Institut für Technische Verbrennung at RWTH Aachen University, Germany, experiments were carried out with a 1.9l GM Fiat Diesel engine.This engine is equipped with a second-generation Bosch Common-Rail injection system and an EDC16 electronic control unit.All relevant engine data are provided in Table 1.The mounting of the engine on the test bench is shown in Figure 2.    A more detailed description regarding the engine, the test cell equipment, and the injection rate measurements can be found in Vanegas et al. [43].
The engine was operated at part-load conditions with a speed of 2000 rpm.For this study, 50 different stationary experiments were carried out with variations in Start of Injection (SOI), external EGR rate, and total fuel mass injected.
Figures 3, 4, and 5 show simulation results obtained with the stand-alone multi-zone model in terms of pressure curve, indicated mean effective pressure of the high-pressure cycle (IMEP HP ), crank angle of 10% and 50% burnt fuel mass (CA50), as well as CO emissions in g/kg fuel in the exhaust gas and combustion efficiency in percent in comparison to test bench measurements for five selected operating conditions (SOI variation from −36 until −16 • CA aTDC with an external EGR rate of 30% and an injected fuel mass of 10.2 mg/cycle).There is a very good qualitative and quantitative agreement between simulation and experiment.
Finally, Figure 6 presents a comparison between numerically and experimentally obtained NO/NO x emissions in g/kg fuel in the exhaust gas for the SOI variation.There is a good agreement.Note that the chemical mechanism used in the stand-alone multi-zone model simulations only accounts for thermal NO (see also the "Chemistry Model'' section), while both NO and NO 2 emissions are measured in the experiments.NO and NO 2 emissions are usually grouped together as NO x emissions, with NO being the predominant oxide of nitrogen produced inside the engine cylinder (see for example Heywood [19]).

GAS EXCHANGE MODELING
The stationary validation described in the previous section is restricted to the high-pressure part of the engine cycle.The usage of the stand-alone multi-zone model within a closedloop controller for PCCI combustion requires that it is capable of predicting the dependency of the controlled variables on the actuators, which are SOI, the external EGR rate, and the total fuel mass injected.The controlled variables are the IMEP and CA50.The latter is directly obtained from a simulation with the stand-alone multi-zone model.The former, however, can only be predicted for the time frame from closing of the intake until opening of the exhaust valves.For this reason, the calculation of the IMEP HP is extended to the gas exchange part of the engine cycle.
The IMEP of the whole engine cycle is composed of the high-pressure part and the part describing the losses caused by the gas exchange.From the viewpoint of automatic control, a model predicting the accumulated losses occurring from exhaust valve opening until intake valve closing is sufficient.Furthermore, as one aim of this work is model reduction, the lowest acceptable approach in terms of detail level and calculation time is desired.Thus, a mean value model for the gas exchange is presented in the following.Nevertheless, if desirable, many commercially available software packages exist for a more detailed calculation of the gas exchange part of the engine cycle (e.g., the WAVE code by Ricardo Software [34] or the GT-POWER code by Gamma Technologies [15]).
The calculation of the indicated mean effective pressure throughout the gas exchange (IMEP GE ) is physically inspired by pumping losses.This approach is derived from the one proposed in Richert [35].It was introduced in Hoffmann et al. [22].The IMEP GE is given by Equation ( 8) includes terms of the volume flow through an orifice and through a throttle, as well as the dependency on the volumetric efficiency or on the density at the intake, respectively.The latter is represented by the fraction p be /T be , which is related to the ideal gas law ρ = p/RT = const.• p/T.The engine-dependent constants x 1 , x 2 , x 3 , x 4 , x 5 , and x 6 need to be determined by a fitting against experimental data.The mean value model for the gas exchange depends on the state directly after and before the combustion process or the manifolds after and before the engine, respectively.In Equation ( 8), these states are denoted with the subscripts " ae '' and " be ''.Inputs to the model from the exhaust manifold (i.e., p ae and T ae ) can be taken from the high-pressure calculation with the stand-alone multi-zone model, while the inputs to the model from the intake manifold (i.e., p be and T be ) become influencing variables to the plant model.These influencing variables can serve as a disturbance input within the closed-loop control.
After appropriately fitting the abovementioned constants, the approach is capable of predicting the IMEP GE within 5% accuracy for most of the 50 operating conditions mentioned above (see Fig. 7).
Combining the stand-alone multi-zone model X0D with this gas exchange model therefore leads to a static model, which can be used to determine the static dependency of the controlled variables IMEP and CA50 on the actuated variables SOI, external EGR rate, and total fuel mass injected.Influencing inputs are the pressure and temperature in the intake manifold.

IDENTIFICATION OF CYCLE-TO-CYCLE ENGINE DYNAMICS
The combination of stand-alone multi-zone combustion model and mean value model for the IMEP GE of the gas exchange does so far not include any dynamics of the controlled variables.Hence, a structure has to be chosen which is suitable for adding the dynamic aspect to the static accurate model.In control theory, two types of models or a combination of both have been established for this purpose: Wiener and Hammerstein models (see for example [29,30]).As the dynamic physical behavior of the real engine shall be enforced on the whole static model, a Wienertype dynamics was chosen in Hoffmann et al. [22].With this choice, the inputs to the static model are overlaid with time attributes, which enforces the dynamic behavior on the combustion simulation.Thus, IMEP and CA50 are consequently affected by the dynamics to be identified.Figure 8 depicts the implemented approach.Here, the static nonlinear behavior is described by the combined integrated static model, or alternatively the gains K 1 through K 6 and the two offset values A 1 and A 2 , respectively.For the system identification described in this section, the latter alternative is used.Thereby, the gains K 1 through K 6 and the two offset values A 1 and A 2 are dependent on the operating point due to the system's nonlinearity.
As only load-transient conditions with a static engine speed are regarded as a first step, the identified dynamics is assumed to be given by a PT 1 .If speed-transient conditions would additionally apply, this approach would either have to be extended by time constants tabulated over engine speed or by a physical model for the mechanical dynamics of pistons, conrods, crank, and additionally gas exchange, respectively.The total model would still have to

Stand-alone model and gas exchange model SOI EGR rate
Model structure used for identifying the system's dynamic time-discrete transfer functions.The dashed line shows the part that may be substituted by the complete integrated static model.be fitted to measurements as well.With the latter approach, the unknown masses and inertias of the drivetrain, as well as parameters of the gas dynamics need to be identified.An example for such an approach is the work by Schulze et al. [37].
A description of the dynamics from the system's inputs SOI, external EGR rate, and total Fuel Mass Injected (FMI) towards the outputs IMEP and CA50 is needed.For this purpose, several approaches for arranging experiments exist in the literature.In this regard, step responses are the most common way for the identification of a system's dynamics.The experimental setup realizes this by using an ES1000 system by ETAS as the central unit for impressing the steps on the actuated variables, as well as for simultaneously collecting the data of the step responses in the controlled variables.More details regarding this Rapid Control Prototyping system by ETAS can be found in the user manual by ETAS GmbH [11].The steps in the actors were enabled by an ETK-connection of the ES1000 system to the engine's electronic control unit (ECU, type EDC16).The physical actor for the SOI and the total fuel mass injected obviously is the Common-Rail direct injector, while for the EGR rate there is no explicit acting device.For this, the EGR rate is changed by the combination of the actors EGR valve and Variable Geometry Turbine (VGT) position.It is impossible to make the engine step in EGR rate.Nevertheless, the system shows a dynamic response to a change in the EGR rate caused by a simultaneous step in EGR valve and VGT position.For the operating points to step between, the corresponding values for EGR throttle and VGT position were determined for a steady state operation.Furthermore, the test bed offers no direct possibility to measure the dynamics of the EGR rate, which is the input and actuating variable for our model.Commonly, in modern ECUs in series applications, the mass flow of fresh air aspirated by the engine is measured instead [16].The EGR rate can be calculated by EGR rate = ṁEGR ṁtrapp.mass = ṁtrapp.mass − ṁair ṁtrapp.mass (9) where ṁair denotes the mass flow of fresh air.The total trapped gas mass ṁtrapp.mass mainly depends on the revolution speed, as well as on the pressure and temperature in the intake manifold.The engine was run at a constant revolution speed.In addition, the intake manifold conditions did not vary significantly during the step response experiments.The total trapped gas mass ṁtrapp.mass is therefore considered approximately constant and, for this reason, the main dynamics affecting the EGR rate in this case is the same but negative dynamics as with the fresh air mass ṁair .This was measured by the engine's ECU, but also handed over to the measuring device ES1000 via ETK.Because the steady state value for the EGR rate is known before and after applying the step in VGT position and EGR valve, the dynamic EGR rate during the tests can be recalculated and used as an input to the model.Another prerequisite for the system identification is the measurement of the target values or outputs of the model, respectively.Both IMEP and CA50 are calculated from the pressure progression measurement.This has to be realized by an additional device, as in series applications no pressure analysis is implemented in the ECU.In this work, a FI2RE system by IAV GmbH is used [23].It analyzes the pressure signal in real time and hands it over to the measurement system ES1000 via CAN.
The step response experiments provide a dataset containing the temporally resolved signals for SOI, external EGR rate, and FMI, as well as the corresponding response in IMEP and CA50.The aim of the extraction of the system's dynamics is the abovementioned Wienertype dynamics.The advantage of this approach is also its disadvantage.By assigning the system's dynamics to the inputs, a relatively simple model structure concerning the dynamic transfer functions is achieved, as their number is reduced to the number of inputs (see Fig. 8).But this also is the main drawback of the approach, because all three inputs influence both outputs.This results in six dynamic dependencies, which have to be described by three transfer functions.The dynamic effects on IMEP and CA50 excited by the changes in SOI thus have to be represented by a transfer function G SOI , the dynamic influence on both outputs caused by changes in the external EGR rate or the total fuel mass injected by a transfer function G EGR or G FMI , respectively.
The transfer function G SOI can be determined from the evaluation of the responses of IMEP and CA50 to steps in SOI.The same sampling rate as used during the measurements with the ES1000 system was used for the timediscrete transfer function.The sampling rate T s of the calculation of IMEP and CA50 is preset by the engine's revolution speed and is given by With n eng.= 2000 rpm, a new value for IMEP and CA50 is calculated every 0.06 seconds.To ensure compliance with Shannon's theorem, the measurements were sampled with 0.01 seconds, and consequently the same sample time was chosen for the time-discrete transfer function.
In mathematics and signal processing, the Z-transform converts a discrete time-domain signal, which is a sequence of real or complex numbers, into a complex frequencydomain representation (see for example Lunze [25]).It is the discrete equivalent of the Laplace transform.The discrete time-domain function f (k * T spl ) of sample time T spl is transformed to a z-domain transfer function F(z) (T spl ) valid for the specified sample time.In this representation, z can be understood as a time shift operator: x(k * T spl − n * T spl ) c . . . . . . . . . . .s z −n X(z) (T spl ) .The measurements indicate the use of a PT 1 -transfer function.With a sample time of 0.01 seconds the identified PT 1 is given by Here, the superscript " '' indicates the temporal attribute of the signal.Selected results of this identification are shown in Figure 9a.For the system identification, the corresponding static gains K 1 through K 6 and the two offset values A 1 and A 2 were adjusted to the operating point to substitute the complete integrated static model.The sampling rate used here obviously is not consistent with the revolution speed.The static integrated model is executed every engine cycle, demanding a sample time of 0.06 seconds for 2000 rpm.To remedy this problem and to enable the implementation of the transfer function into the integrated model, the time-discrete transfer function is resampled with the new sample time of 0.06 seconds to G(z) SOI,(0.06)= SOI SOI = 0.4686 z − 0.5314 (12) In an analogue manner, the transfer functions G(z) EGR and G(z) FMI can be determined from the evaluation of the responses of IMEP and CA50 to the dynamic EGR rate or dynamic total fuel mass injected signal, respectively.Thereby, the EGR rate signal is recalculated as described above.The corresponding gains and offset values in Figure 8 are again adjusted such that the static nonlinear transfer behavior is met.The fitting of the system's dynamics is evaluated with a measuring sample time of 0.01 seconds and resampled to 0.06 seconds afterwards.The influence of the external EGR rate is given by G(z) EGR,(0.01)= EGR rate EGR rate = 0.015 z − 0.985 (13) and can be resampled with T s = 0.06 seconds to Using Equation ( 13), the results shown in Figure 9b, among others, could be achieved.
For the dynamic dependency of the model's outputs on the total fuel mass injected the transfer function was determined, which leads to the resampled transfer function G(z) FMI,(0.06)= FMI FMI = 0.3572 z − 0.6428 (16) The results for the response in IMEP and CA50 to a step in FMI from 15 mg/cycle to 25 mg/cycle and back (with an SOI of −35.7 • CA aTDC and an external EGR rate of 45%) are depicted in Figure 9c.Note that for the identification of G(z) FMI,(0.01), G(z) EGR,(0.01) has to be identified first, as it is impossible to vary FMI without changing EGR, if VGT and EGR valve are not adjusted to the new load condition.

INTEGRATED MODEL
For application within a closed-loop control simulation, the stand-alone multi-zone combustion model was transferred to an environment suitable for the conception and testing of controllers.The stand-alone multi-zone model X0D written in FORTRAN 77 was embedded into a Matlab r /Simulink r FORTRAN s-function enabling the simulation of the stand-alone multi-zone model from within Matlab r /Simulink r .Moreover, the gas exchange model was implemented into this s-function.The three identified dynamic transfer functions were added by means of timediscrete PT 1 -dynamics within appropriate function blocks.The simulation environment is realized such that the initial in-cylinder conditions reflect the engine-out conditions of the prior cycle.This establishes the cycle-to-cycle connectivity.However, as the external EGR was cooled, the measurements did not show a strong dependency of the intake conditions on the engine-out conditions, except for the external EGR composition when varying the total fuel mass injected.
Figure 10 shows IMEP and CA50 obtained from the step response experiment, the system identification with the corresponding identified discrete transfer function G(z) SOI,(0.01), and the transient simulation with the integrated model using G(z) SOI,(0.06) for an SOI step from −20.7 to −30.7 • CA aTDC and back with an external EGR rate of 30% and an injected fuel mass of 10.2 mg/cycle.The simulation results are in very good agreement with the measurements.In particular, the dynamic step responses are reproduced well.Similar simulation results (not shown here for brevity) were achieved for all other step response experiments discussed above.This validates the integrated model composed of the stand-alone multi-zone model, the gas exchange model, and the identified system dynamics.
The strict derivation of the stand-alone multi-zone model from a detailed three-dimensional CFD approach integrates the detailed knowledge from combustion simulation tools into closed-loop control.This novel procedure establishes a broad field of possibilities for testing completely new controlled process variables.In this context, Figure 11 displays further outcomes of the transient simulation for an SOI step from −20.7 to −30.7 • CA aTDC and back with an external EGR rate of 30% and an injected fuel mass of 10.2 mg/cycle.Exemplary, the NO emissions, the maximum pressure, the burning time from CA10 to CA90, and the ringing intensity are shown.
The combustion noise level is quantified applying the ringing intensity correlation developed by Eng [10].This correlation relates the ringing intensity to the maximum pressure rise rate, the maximum pressure, and the speed of sound.The ringing intensity is expressed as In Equation ( 17), dp/dt max denotes the maximum pressure rise rate, p max the maximum pressure, and T max the maximum temperature.β is a scaling factor determined from the experimental data.It relates the pressure pulsation amplitude to the maximum rate of pressure rise.In this work, it is set to 0.005 ms.In Eng [10], a value of 0.05 ms was applied.Figure 11 reveals that the NO emissions increase with advanced SOI timings due to an increasing maximum pressure or maximum temperature, respectively.Advancing the SOI also increases the burning time.This reduces the pressure rise rate, leading to a decreasing ringing intensity.These results indicate the potential of the integrated model in the framework of closed-loop control development.Beyond the controlled variables IMEP and CA50, the integrated model is capable of predicting various additional process variables, with four selected ones being briefly discussed in this subsection.These additional process variables allow for defining further process constraints or may even replace the original controlled variables, depending on the specific control task.

SUMMARY AND CONCLUSIONS
Closed-loop simulations are a necessary and common tool in the development process of controllers.A computationally efficient stand-alone multi-zone model was employed that is capable of describing the combustion characteristics for the high-pressure part of the engine cycle.The controller to be developed shall actuate SOI, external EGR rate, and total fuel mass injected to control the IMEP of the whole engine cycle and the CA50.The IMEP HP of the highpressure engine cycle and the CA50 can be extracted from simulations with the stand-alone multi-zone model.The former was combined with a mean value model for the losses of the gas exchange to calculate the IMEP of the whole engine cycle.The combination of these models leads to an accurate

FUTURE WORK
Future experimental work will have to reveal additional transient experimental results (e.g., NO emissions, maximum pressure, burning time from CA10 to CA90, etc. (see Fig. 13) in order to further assess the performance of the integrated model.The dynamic model will build the framework in which a controller for the PCCI combustion process will be developed.With the developed plant model, this controller will be laid out and tested.
The actuated variable external EGR rate is theoretical as long as there is no corresponding actor to it.Because of this, the overall goal of this research project is to develop not only a controller for the PCCI combustion process, but also an additional controller enabling the fastest possible changes in external EGR rate by simultaneously adjusting EGR valve and VGT position [9].This controller, in combination with its actuated variables EGR valve and VGT position, will serve as the virtual actor for the external EGR rate.

Figure 2
Figure 2 Engine test bench.

Figure 3 Figure 4 CA10
Figure 3Average cylinder pressure (left) and IMEP HP (right) for five selected operating conditions (SOI variation from −36 until −16• CA aTDC with an external EGR rate of 30% and an injected fuel mass of 10.2 mg/cycle).Comparison between stand-alone multi-zone model simulation and experiment.

Figure 5 CO
Figure5CO emissions in g/kg fuel in the exhaust gas (left) and combustion efficiency in percent (right) for five selected operating conditions (SOI variation from −36 until −16• CA aTDC with an external EGR rate of 30% and an injected fuel mass of 10.2 mg/cycle).Comparison between stand-alone multi-zone model simulation and experiment.

Figure 6 NO
Figure6NO/NO x emissions in g/kg fuel in the exhaust gas for five selected operating conditions (SOI variation from −36 until −16• CA aTDC with an external EGR rate of 30% and an injected fuel mass of 10.2 mg/cycle).Comparison between stand-alone multi-zone model simulation and experiment.

Figure 7 IMEP
Figure7IMEP GE of the gas exchange for all 50 experiments mentioned above.Comparison between mean value model calculation and experiment.

9
Figure 9 Comparison between experiment and system identification with different identified discrete transfer functions G(z).a) G(z) SOI,(0.01) .IMEP (left) and CA50 (right) for an SOI step from −20.7 to −30.7 • CA aTDC and back with an external EGR rate of 30% and an injected fuel mass of 10.2 mg/cycle.b) G(z) EGR,(0.01) .IMEP (left) and CA50 (right) for an EGR step from 30 to 45% and back with an SOI of −20.7 • CA aTDC and an injected fuel mass of 10.2 mg/cycle.c) G(z) FMI,(0.01) .IMEP (left) and CA50 (right) for an injected fuel mass step from 10.2 to 17.0 mg/cycle and back with an SOI of −35.7 • CA aTDC and an external EGR rate of 45%.

Figure 10 IMEP
Figure 10 IMEP (top) and CA50 (bottom) for an SOI step from −20.7 to −30.7 • CA aTDC and back with an external EGR rate of 30% and an injected fuel mass of 10.2 mg/cycle.Comparison between experiment, system identification, and integrated model.

Figure 11 Additional
Figure 11Additional process variables obtained from a transient simulation with the integrated model.NO emissions, maximum pressure, burning time from CA10 to CA90, and ringing intensity (from top to bottom) for an SOI step from −20.7 to −30.7 • CA aTDC and back with an external EGR rate of 30% and an injected fuel mass of 10.2 mg/cycle.