Development of On-Board In-Cylinder Gas Flow Model for Heat Loss Estimation of Diesel Engines

Diesel engines have been demanded to further increase the thermal efficiency through precise engine control under transient driving conditions, especially, it is essential to optimize the fuel injection timing and quantity cycle-by-cycle. Conventionally, fuel injection have been controlled by control maps, which resulted in large numbers of experiments and increase in cost. In order to overcome these problems, the present study focused on the model-based control and developed the onboard gas flow model, because the heat loss is affected by the turbulent intensity. Firstly, a validation of the CFD simulation is evaluated. The CFD simulation was used to validate the developed models and to determine unknown parameters used in the model. Secondly, modeling of in-cylinder gas flow is presented. To estimate the injection timing within 0.5 deg. against the target value, the heat loss must be estimated within the error range of 7.6%. Finally, as results, the error of heat loss obtained from gas flow model was 1.6%, and gas flow model fully met the requirement of tolerance range. From the viewpoint of calculation time, the calculation time of the model was 50.6 s per cycle, and thus the model is capable of the use of on-board applications.


Introduction
Diesel engines have been complying with strict emission regulations by employing post-combustion treatment systems such as selective catalytic reduction (SCR) and diesel particulate filters (DPF) [1][2].As the regulations become stricter to ensure low emissions under real driving conditions, diesel engines must increase the thermal efficiency and reduce the emissions through precise combustion control without relying on the post-combustion treatment systems.More precise combustion control is possible by optimizing the fuel injection timing for each cycle with feed-forward control [3][4][5].Conventionally, fuel injection has been controlled by control map which is developed from many laboratory experiments that are costly and time consuming [6].Alternatively, a new control method called model-based control directly installs on-board models to ECU for calculations each cycle.On-board model can predict state of in-cylinder at each cycle using various information from engine sensor without an engine control map, and are applicable even if engine operation condition change [7].Traditionally, for model-based control, on-board experimental model have been used.On the other hand, in theory, on-board physical model which is based on fundamental description of physics are applicable to ECU even if kind of engine change.Thus, model-based control can prevent thermal efficiency from reducing and reduce the number of experiments in the development process [8][9].Since calculation time of model to be applied on-board must be faster than the actual engine cycle, for instance, it is required to be faster than 1 ms which is time spent in one cycle at maximum revolution speed of general diesel engines.Even a 1-D simulation and a 0-D simulation as well as a computational fluid dynamics (CFD) simulation cannot satisfy this demand [10].Furthermore, accurate estimation of heat loss is a crucial factor when estimating the gas temperature for optimization of fuel injection timing.To estimate the injection timing within 0.5 deg.against the target value, the heat loss must be estimated within the error range of 7.6 %.
Although the models for model-based control of diesel engines [3] and gasoline engines [11] have been developed previously, they were either empirically estimated or low in robustness.The previous physicsbased heat loss models [12][13] could not reduce the computational load, because it is necessary to calculate flow velocities in boundary layers.To achieve the requirement of calculation load, such a model has been previously presented by Ravi et al. and Yamasaki et al. dealing with homogeneous charge compression ignition (HCCI) engines [14][15].In both studies, the in-cylinder state of a single cycle was discretized into several representative points.Although Inagaki et al. and Fujikawa et al. have developed heat loss models for diesel engines with high robustness, however, they did not consider gas flow effects [16][17].The turbulence is caused by numerous factors including axial direction, squish, swirl, injection flows, which greatly affects heat loss [18].Various gas flow models with different calculation speeds and based on different concepts have been developed.The CFD simulation can calculate the local distribution of flow velocity and analyze turbulent flow field in detail [19][20][21].However, the CFD simulation imposes a heavy computational load and is not necessarily useful for engine control.Hiroyasu et al. have developed a spray tip penetration model from experimental information and theoretical analysis of the intact length of high speed liquid jet [22].Inagaki et al. improved the model for combustion simulation soft for diesel engines [23].Murakami et al. assumed that the charge in each of the two partial volumes rotates as a solid-body swirl and was confirmed experimentally by LDA measurements in combustion chamber [24].Schubert et al. have developed an axial flow model and a squish flow model for heat transfer model [25].However, these four models are not targeted by the onboard application.
The objective of this study is to develop an on-board gas flow model which estimates the heat loss and gas temperature before the timing of fuel injection.Gas flow of diesel engine is mainly divided into four different flows, such as axial flow, squish flow, swirl flow and injection flow.There are few swirl models and heat loss models considering four kinds of flow.Therefore, in this study, swirl model was developed and applied to a heat loss model developed in author's previous study [26].This model is derived from the continuity equation and the energy equation with a low computational load considering convection as well as conduction with formulation of thermal boundary layers.Other gas flow models were referred to other studies [22,[24][25] described above and modified to be applicable to ECU.Since the heat loss due to convection varies vastly depending on chamber regions, the engine cylinder was divided into six different regions.The gas flow velocity and its fluctuation caused by turbulence at each region are considered with influences of axial direction, squish, swirl, and injection flows based on simple physical theories.In addition, to satisfy a requirement of calculation time, the in-cylinder state of a single cycle was discretized into several representative points, such as the valve timing (IVO, IVC, EVO and EVC) and the points to maintain the calculation accuracy of the gas flow model.This paper presents in-cylinder gas flow models of diesel engines with a low computational load for model-based control and identification of undetermined constants, which are required to improve on-board heat loss model developed in author's study.In Section 2, validation of the CFD simulation was evaluated by comparing in-cylinder pressure and heat release rate obtained from the CFD simulation and those from experiments.The CFD simulation was used to validate the developed models and to determine unknown parameters used in the model.Modeling of gas flow and calculation procedure of determined constants are presented in Section 3. In Section 4, results of identification of determined constants are shown, and to evaluate the validation of the gas flow model, results of velocity of gas flows and heat loss obtain from gas flow model are compared with those of the CFD simulation.Moreover, the computational load was evaluated at last to verify its capability of on-board application.The paper finishes with conclusions in Section 5.

Condition of CFD Simulation
Figure 1 is a schematic of the computational model used in the CFD simulation.The engine specifications are listed in Table 1, and the calculation conditions are shown in Table 2.The CFD simulation was used to validate the developed models and to determine the unknown parameters used in the model under the steady driving condition.The fuel was injected by three stages, whose timings were listed in Table 2.

Validation of CFD Simulation
Unsteady simulation for identifying determined constants and evaluating gas flow models were performed using a commercial code (Converge, Convergent Science Inc.), whose models are listed in Table 3.Another method called FVM, adaptive mesh refinement code, which automatically refines the grid based on fluctuating and moving conditions such as temperature and velocity, was used for creating 3-D meshes.Spray breakup was simulated by N.A. and KH-RT models [27][28][29], and SAGE model [30] was used to solve combustion.Wall heat flux was calculated using the law of the wall, which was shown good agreement with experimentally-measured wall heat flux [31][32][33][34] and given below.The flow was computed as compressible, since in the cylinder, volume of gas changes in the closed space.The boundary conditions of inlet are set for the transient mass flow rate and temperature, on the other hand, those of wall of cylinder are set for the steady temperature.Velocity-pressure coupling is addressed with the well-known PISO strategy.RNG k- model [35], which is one of the RANS model and widely used in engine simulations, was employed.F. Perini et al. [36] and P. C. Ma et al. [37] guaranteed the accuracy of RANS model by comparing results of RANS model with results of PIV.The initial time step was set to 1 × 10 −5 s, and the maximum mesh number was 2,127,510 and the minimum one was 112,243.
  where q w,k is the heat flux from the gas to the cylinder   For the validation of the CFD simulation, the incylinder pressure and heat release rate (HRR) obtained from the CFD simulation were compared to experimental results as shown in Figures 2 and 3.The average error of in-cylinder pressure was evaluated 4.0%, and the simulated pressure and HRR histories agreed with the experimental data.

Reduction of Computational Model
The in-cylinder states during a single cycle were discretized into 22 points for reducing the computational load for on-board applications.These points are specific points of a common diesel cycle, such as the valve timing (IVO, IVC, EVO and EVC) and the points to maintain the calculation accuracy of the gas flow model.22 discretized points and the corresponding crank angles are shown in Figure 4 and Table 4, and each discretized point will be referred as "point".

Heat Loss
Heat loss is the total amount of heat transfer between the in-cylinder gas and cylinder wall, and is varied depending on the regions of the cylinder.Therefore, the cylinder was divided into six regions such as piston top (k = 1), cavity side (k = 2), cavity bottom (k = 3), liner (k = 4), inner head (k = 5), outer head (k = 6) as shown in Figure 5.The wall heat flux was estimated for each region, and the total heat flux is written as the summation of heat fluxes due to the conduction (the first term on the right-hand side in Eq. ( 3)) and convection (the second term on the right-hand side in Eq. ( 3)) [26].

ˆ' 4
where q w,i,k is the heat flux from the gas to the cylinder wall of region k at discretized point i [W/m 2 ], C λ is the ratio of thermal conductivity of air to the gas temperature [W/(m•K 2 )], P 0 is the in-cylinder pressure at TDC during the intake stroke (point 1) [Pa], τ is dimensionless time [-], T w,k is the wall temperature of region k in the previous cycle [K], ψ (=0.4) is the Karman constant [-], , ˆik u is the turbulent intensity at region k at point i [m/s], whose detail was described in Section 3.3, and t is the elapsed time from TDC during the intake stroke [s].Additionally, the rate of heat transfer from the in-cylinder gas to the in-cylinder wall surface can be obtained by multiplying the heat flux by the area of each region.The heat loss is estimated by time-integrating the rate of heat transfer.

Turbulent Intensity
For the calculation of the heat flux with Eq. ( 3), the turbulent intensity, , ˆik u , is required, which is defined as the vertical turbulence component of the gas flow at the wall surface and whose value greatly differs depending on the cylinder region.Thus, the present study modeled the turbulent intensities at the six regions as the combinations of four different kinds of flows, such as axial direction, squish, swirl and injection flows, which were calculated as:

The Axial Direction Flow Model
The axial direction flow velocity is changed by the vertical movement of piston.Assuming the velocity near the cylinder head is negligibly small and that near the piston is equal to the piston speed, u p , the instantaneous axial direction flow velocity, u axial , is modeled as:

The Squish Flow Model
A schematic of the squish flow inside the cylinder is shown in Figure 6.The radial direction flow velocity heading from upper piston top (Region 2 in Figure 6) to upper cavity (Region 1) was defined as squish 1 velocity, u sq1 , and the axial direction flow velocity heading from region 1 to cavity (Region 3) was defined as squish 2 velocity, u sq2 .During the compression stroke, assuming that the mass of in-cylinder gas is preserved, u sq1 is modeled as: where A g is the side surface area of region 1 [m 2 ], A sq is the surface area of piston top [m 2 ], V bowl is the volume of cavity [m 3 ], and V is the volume of the entire combustion chamber [m 3 ].Similarly, u sq2 is modeled as:

The Swirl Flow Model
Assuming the swirl flow consists of two rigid body flows at cavity region (r = 1) and squish region (r = 2), u θ is modeled as: , where x k is the distance from the center of cylinder to the wall (r = 1 when k = 2, 3 and 5, and r = 2 when k = 1, 4 and 6) [m], and ω r is the angular velocity [rad/s], which is obtained using Eq. ( 13).Applying the law of conservation of momentum to cavity region and squish region, the total momentum inside the cylinder is modeled as:  12) and ( 13), we employed the previous models for M f , M sq , and M ν [7] and developed the new models for M in and M out .When the area of intake valve opening is assumed to be A in [m 2 ], the inflow angular momentum by intake M in is modeled as: where r is the radius of the intake port [m], v is the gas flow velocity passing through the intake valve [m/s], n is the vertical unit vector toward the interface [-], and ρ is the density of intake gas [kg/m 3 ].The inflow angular momentum of helical port and tangential port was estimated based on Eq. ( 16) as the basic equation.
Figure 7 shows a schematic of intake and exhaust ports of diesel engine, and Figure 8 depicts schematic of tangential port and helical port viewed from the axial direction.In case of tangential port, assuming that fresh air flows in from one side of intake valve, M in is modeled as: In case of helical port, assuming that fresh air flows at an angle of 45 deg.from the whole area of intake valve to the vertical axis, M in is modeled as: where a is the distance from the center of the cylinder to the center of the intake port [m], v Tan is the intake velocity of tangential port [m], V Hel is the intake velocity of helical port [m], φ [rad], θ [rad], and l [m] are angle or distance shown in Figure 8.
A schematic of exhaust port from axial direction is shown in Figure 9. Assuming that in-cylinder gas flows out from one side of exhaust port entrance, similar to inflow angular momentum, M out is modeled as:

The Injection Flow Model
Hiroyasu and Arai model [22] was used to calculate the fuel spray penetration during the injection period.The fuel spray penetration at the tip of injection was assumed to be attenuated by air resistance.Then, as the injection flow collides with the wall, it assumed to flow along the wall.Assuming the spray angle does not change after collision, the tip velocity of injection, u inj , at time t is modeled as: where t inj is the injection period [s], u ie is the tip velocity of injection [m/s], ΔP is the pressure difference between atmosphere and nozzle sac [Pa], d 0 is the diameter of nozzle hole [m], φ is the velocity damping coefficient considering air resistance [-].

Calculation Procedure of Turbulence Intensity Coefficient
The turbulent intensity coefficients, C a , are expressed as the ratio of mainstream velocity of gas flow and its fluctuating velocity component at the wall surface, and their values vary depending on gas flows.In the present study, the turbulent intensity coefficients of the axial direction, squish, swirl, and injection flows were identified by using the CFD simulation, which is calculated as: where u' is the fluctuating vertical velocity at the wall surface [m/s], and u  is the mainstream velocity of gas flow [m/s], which was determined using the flow velocity parallel to each wall surface of cylinder by the CFD simulation.On the other hand, RNG k-ε model used in this analysis is considered turbulence as isotropic, and thus, u' was calculated using turbulent kinetic energy [m 2 /s 2 ], k, obtained from the CFD simulation as: For the axial direction, squish, and swirl flows, the area-averaged value of u' at the wall surface was used.On the other hand, for the injection flow, the maximum value of u' at the wall surface was used.Furthermore, for the calculation of u' and u  of the axial direction and the squish flows, the CFD simulation was carried out without initial swirl flow and fuel injection from IVC to the end of injection period, because both flows are induced by only the movement of piston.For the swirl flow, the simulation was carried out without fuel injection during the main injection period, and for the injection flow, the simulation was carried out without initial swirl flow during the main injection period.

Evaluation of In-Cylinder Gas Flow Models
The validation of the gas flow model was evaluated by comparing the CFD simulation and the model results.Under the condition as listed in Table 2, the axial direction flow velocity, squish 1 and 2 velocities, ω 1 and ω 2 , were compared to the CFD simulation results.As representative results, the comparisons of squish 1 and squish 2 velocities, and ω 1 and ω 2 obtained from the CFD simulation and those from the model are shown in Figures 10 and 11.  10, before TDC (CA 360 deg.), the squish 1 velocities of the CFD simulation have the negative values because the gas in the squish area flows into the cavity.On the other hand, after TDC, the squish 1 velocities have the positive values because the gas in the cavity returns to the squish area.Also, the model reproduced the tendency of the CFD simulation well.As shown in Figure 11, it is indicated that ω 1 increases when the piston approaches toward TDC.Near TDC, since the gas enters from the squish area to the cavity, the mass of gas and angular momentum in the cavity increases.Then, ω 1 and ω 2 fall slightly due to the increase in the moment of inertia and the effect of wall friction.Figure 11  is slightly larger than that of ω 2 .This can be attributed to one of the assumptions of the swirl model, which is a rigid body flow.In reality, the flow is not a rigid body flow and the swirl center might be deviated.

Identification of the Turbulent Intensity Coefficient
With the procedure described in Section 3.8, the turbulent intensity coefficients for each flow were calculated.The turbulent intensity coefficients, C α , of each flow are shown in Table 5.As an example, the time evolution of the turbulent intensity coefficient of axial direction flow, C axial , is shown in Figure 12, and C axial was evaluated to be 0.028 by averaging the results during the main injection period (CA 365 -374 deg.).

Evaluation of Heat Flux and Heat Loss Using Heat Loss and Gas Flow Models
Under the conditions listed in Table 2, the wall flux and heat loss estimated by the model were compared to those by the CFD simulation.The wall heat flux at each region was calculated by Eq. ( 3) with the gas flow model.As an example, the comparison of the wall heat flux at the piston top and cavity side (k = 1 and 2) between the CFD simulation and the model are shown in Figure 13.The heat flux was increased during main injection period (CA 365 -374 deg.) at both regions.The heat fluxes at the cavity side were larger than those at the piston top, due to the diffusion flame directly colliding with the wall.The model was able to reproduce the heat flux obtained by the CFD simulation at all regions reasonably well.However, during combustion, the performance of the model deteriorated, because this heat loss model does not include a combustion model [38].
The heat losses at the piston, head, liner and total area from CA 0 deg. to 482 deg.were calculated with the gas flow model.The results are compared to the CFD simulation as shown in Figure 14.The RMS errors of heat loss at the piston, head, liner and total area between the model and the CFD simulation were summarized in Table 6.Although the errors at the head and liner are slightly larger than those at the piston and total area, their influence on the heat loss of total area is small, because the heat losses at the head and liner are much smaller than those at the piston as shown in Figure 14.Also, this result verified the validity of the gas flow model for the accurate estimation of heat loss.

Evaluation of the Computational Load
The computational load of the model was evaluated to verify the applicability to the model-based control.Calculations are performed on a typical personal computer (OS: Windows10 Home 64bit, CPU: Intel Core i5-4200U@1.6GHz,Memory: 4GB), whose time for a cycle was approximately 50.6 μs.By contrast, one engine cycle at an engine speed of 2000 rpm takes MHz) runs slower than the PC used in this study, it is apparent that the calculation time of the model is faster than the time of a single engine cycle.The model is capable of being implemented in a feed-forward controller for real-time prediction of heat loss.

Conclusions
The present study developed an on-board incylinder gas flow model, and identified undetermined constant, the turbulent intensity coefficient.Firstly, a validation of the CFD simulation was evaluated by comparing in-cylinder pressure and heat release rate obtained from the CFD simulation and those from experiments.Secondly, modeling of gas flow and calculation procedure of determined constants were presented.Thirdly, results of of determined constant was shown.Moreover, to evaluate the validation of the gas flow model, results of velocity of gas flows and heat loss obtained from gas flow model were compared with those of the CFD simulation.Finally, the computational load was evaluated to verify its capability of on-board application.The following conclusions were obtained: 1.The gas flow model was consisted of axial direction flow, squish flow, swirl flow, and injection flow, which were calculated at 22 points for a cycle to reduce the computational load.The gas flow velocities of the developed model reproduced the tendency of the CFD simulation results.The turbulent intensity coefficients for axial direction flow, the squish flow, the swirl flow, and the injection flow were estimated to be 0.028, 0.039, 0.041 and 0.102, respectively.2. The heat fluxes and the heat losses were evaluated using the gas flow model and the turbulent intensity coefficients.The RMS error of heat loss at total area between the model and the CFD simulation was estimated 1.6%, and fully met the requirement of tolerance range within 7.6 %. 3. The calculation time of the model was 50.6 μs per cycle.It is faster than the time of a single engine cycle, considering the CPU clocks of present massproduced ECUs typically run 10 % slower than the PC used in this evaluation.Thus, the model is capable of the use of on-board applications.

Figure 1 .
Figure 1.Schematic of computational model k is the wall temperature [K], P rm is the molecular Prandtl number [-], y is the distance between wall surface and center of gravity of first cell layer [m], P rt is the turbulent Prandtl number [-], y + is the dimensionless distance [-], B is the parameter of wall function [-].

Figure 2 .
Figure 2. Comparison of in-cylinder pressure between experiment and 3-D CFD simulation

Figure 3 .
Figure 3.Comparison heat release rate between experiment and 3-D CFD simulation

Figure 4 .
Figure 4. Relationship between the in-cylinder pressure of a typical diesel engine and the dots corresponding discretized points

Figure 5 .
Figure 5. Schematic of in-cylinder regions

)
where u θ,k is the in-cylinder circumferential direction flow velocity due to swirl flow [m/s], u sq1 is the radial direction flow velocity due to squish flow [m/s], u sq2 is the axial direction flow velocity due to squish flow [m/s], u axial is the axial direction flow velocity due to piston movement [m/s], u inj is the tip velocity of injection [m/s], R 1 is the radius of cavity [m], R 2 is the radius of cylinder [m], and C α (α = axial, squish, swirl and inj) is the turbulent intensity coefficient[-].

Figure 6 .
Figure 6.Schematic of squish flow model ) where I r is the moment of inertia [kg•m 2 ], M in is the inflow angular momentum by intake [N•m], M f is the decrease in angular momentum due to wall friction [N•m], M sq is the angular momentum due to squish flow [N•m], M ν is the decrease in angular momentum due to fluid friction [N•m], and M out is the outflow angular momentum by exhaust [N•m].To calculate Eqs. ( where ρ is the density of in-cylinder gas [kg/m 3 ], v out is the outflow velocity of exhaust port [m/s], A out is the area of exhaust valve opening [m 2 ], and , φ [rad], θ [rad], and l [m] are angle or distance shown in Figure 9.

Figure 7 .Figure 8 .Figure 9 .
Figure 7. Schematic of intake and exhaust ports of diesel engine

Figure 10 .Figure 11 .
Figure 10.Comparison of squish velocity obtained from CFD simulation with the model shows that ω 2 of the model can generally track the tendency of the CFD simulation.On the other hand, the estimation error of the model of ω 1

Figure 13 .Figure 14 .
Figure 13.Comparison of the wall heat flux at the piston top and cavity side obtained from CFD simulation and the heat loss model about 60.0 ms.Although a typical ECU (about 240

Table 1 .
Engine Specifications

Table 3 .
Simulation model used in 3-D CFD simulation

Table 5 .
The turbulence intensity coefficient of each flow

Table 6 .
Root mean square error of total heat loss from CA 0 deg. to 482 deg. between CFD simulation and the model