This work is based on two different types of numerical analysis. Firstly, a basic 0D analysis based on the concept of the mixing line is carried out to achieve some general conclusions on the combustion characteristics and the influence of the fuel temperature on the evolution of the in-cylinder properties. This type of analysis is also useful to define the fuel properties that could improve this specific technology. It is followed by CFD simulations to extend the conclusions drawn with the 0D model and to capture the influence of the injection parameters on the fuel stratification and the combustion evolution. Finally, the CFD model is coupled to the evolutionary algorithm and processing tools used for the optimization.
3.1. Mixing Line
The goal of this preliminary study is to model, in the simplest way possible, the evolution of the fuel jet temperature during the mixing process, in order to understand its influence on the Ignition Delay (ID) of the mixture. In this way, it is possible to determine the parameters that have more influence on the combustion evolution. This analysis is performed using Cantera, an open-source suite of tools for problems involving chemical kinetics, thermodynamics and transport processes, and the adopted model is described here. The temperature evolution as a consequence of the mixing process is calculated from the first law of thermodynamics:
where
is the fuel temperature at the exit of the injector hole,
is the average in-cylinder temperature, and
is the specific heat. The fuel is always considered in the gaseous state; thus, the heat of vaporization is never considered in the model, as it is suitable only for supercritical injection.
The fuel and air mass are related by the equivalence ratio:
The total mass is set from the state equation for ideal gas, considering the in-cylinder thermodynamic conditions at the end of the direct injection. The local temperature equation is simply a function of the equivalence ratio (
), and it models the temperature evolution in an ideal instantaneous spray where the cells do not interact between each other; the diffusion phenomena are therefore neglected. For each point defined in the
space, a constant pressure reactor equation is used to set the local ignition delay. This step requires the use of a chemical–kinetic mechanism. For this purpose, the ID is defined as the time when the global normalized progress variable reaches a value equal to 0.1 [
7]. Finally, a relation such as
is obtained. This analysis aims to model, in the simplest possible way, an abstract injection by defining the temperature stratification and the resultant stratification of ID. It is firstly necessary to set the temperature evolution in the range from Z = 0—a condition not affected by the injection—to Z = 1—the condition at the exit of the nozzle hole. It is important to notice that in real conditions, at the start of the DI, the port injected fuel is already homogeneously mixed in the charge. The presence of this fuel has been neglected since it is a constant and does not affect the trend of the evolution. Furthermore, this analysis is qualitative and not quantitative; therefore, this omission does not influence the general conclusion achieved. The first relevant consideration is that, to have a significant fuel stratification, the injection timing for these type of LTC strategies is close to the TDC (between −40 and −20 CAD), and thus the in-cylinder temperature is close to 800/900 K. Moreover, the fuel still cools down the charge temperature in the vast majority of the possible operating conditions.
Figure 3 represents the temperature evolution of the mixing line considering a charge temperature equal to 900 K and a sweep of fuel temperature between 550 and 750 K; different properties can be discussed here. At first, the range of mixture fractions interested by the combustion is characterized by a small fuel mass. Therefore, the variation of the charge temperature is much more influential with respect to the fuel temperature. The sensibility of the LTC technologies to variation of the charge conditions is well known, and it is responsible for their limited control.
Once the temperature evolution has been set, it is possible, for each temperature-equivalence ratio point, to perform a constant pressure reactor calculation to define the corresponding ignition delay. It is known that one of the most important fuel properties for stratified LTC technologies is the so-called phi sensitivity [
8], defined as follows:
This property reflects the effectiveness of the stratification, since an increase of the equivalence ratio follows a decrease of the ignition delay. This property promotes a slower combustion, avoids a sudden pressure increase and allows the maximum load achievable by the engine to be increased, decreasing the knock tendency. Note that it is normalized to properly compare conditions that have very different ignition delays, and it has been chosen to use the negative sign to use a positive value for correct behavior.
Figure 4 depicts the behavior of iso-octane under the following conditions: the in-cylinder temperature and pressure are 900 K and 65 bar, while the fuel ranges from 550 to 700 K. These values are representative of the thermodynamic state of the medium load case of the engine (defined as case b in the further section of validation) at −15 CAD. They have been chosen to be relatively close to the start of the main combustion to guarantee a good fuel reactivity. Since these conditions come after the end of the injection, the difference between the fuel and the in-cylinder temperature is higher, and thus the effect of the fuel temperature is probably quite over-estimated. However, this issue is limited since the diffusion of the direct injected fuel is not instantaneous, and the mixing continues well after the injection. The first important result is that iso-octane shows a negative phi sensitivity in the whole range of interest. It is also interesting that the ID of this fuel is almost not influenced by the fuel temperature—this characteristic is probably the result of the high stoichiometric air–fuel ratio of hydrocarbons.
It must be noted that this analysis was carried out in gaseous conditions. The total ignition delay of a supercritical injection is, however, always smaller than a cold injection, since the fuel evaporation is much faster. The results of the analysis suggest that, for iso-octane, the benefits resulting from a further increase of the injection temperature above the supercritical one are minimum. These results are useful for the determination of the initial conditions for CFD simulations and to have an idea of the expected combustion evolution.
3.2. CFD
The CFD analysis is performed in LibICE, a set of applications and libraries for Internal Combustion Engine (ICE) simulations, with pre and post-processing implemented in the OpenFOAM framework, developed by the Politecnico di Milano–ICE group [
9,
10,
11,
12]. It provides tools and solvers to model physical and chemical phenomena, such as liquid spray dynamics and evolution, combustion processes and exhaust gases after treatment. The simulations performed here are defined as a closed cycle. They range from Inlet Valve Closing (IVC) to Exhaust Valve Opening (EVO) and allow the combustion characteristics to be captured without modeling the air exchange process; thus, the required computational time is significantly reduced. The models used for combustion, turbulence and injection are described below.
3.2.1. Combustion Model
This technology is based on a stratified spontaneous ignition premixed combustion. The engine operates with commercial pump gasoline, while the simulations are performed using iso-octane. The combustion model used is based on the flamelet approach, and it is the well-known Tabulated Well Mixed (TWM) [
13]. The main advantage of the flamelet approach is the decoupling of the flame structure from the flow dynamics [
14]. This property is obtained via the introduction of the conserved scalar: mixture fraction Z, commonly defined based on the elements conservation. In the well-mixed model, each computational cell of the reactive mixture is treated as a closed homogeneous system, neglecting turbulence–chemistry interaction for the local flow conditions. Furthermore, this model allows the chemical computation to be decoupled from the CFD simulations, since the reaction equations are previously solved and the relevant quantities stored in a table as a function of the state variables of the system and the combustion progress variable. The adopted reduced chemical kinetic mechanism is a gasoline surrogate reduced mechanism developed by the Creck Group [
15].
3.2.2. Injection Model
The injection phenomenon in this particular technology is challenging to model and investigate since the physics and the modeling techniques change with the passage from liquid to super-critical injection. Commercial pump gasoline is a complex, multi-component fluid with a boiling point that ranges from 533 to 618 K [
16], while the reference fuel used in the simulation (iso-octane) has a critical point equal to 544 K and 25.6 bar [
17]. The super-critical injection of gasoline has not been deeply investigated up to now; however, a larger bibliography is present in the aerospace field for different liquids such as LN
. It has been experimentally proved that super-critical injection shows similar characteristics with respect to turbulent gas jets when the injection pressure is above the critical pressure [
18], as in the analyzed case. Conversely, the behavior is close to liquid jets for sub-critical pressure. An internally developed virtual gaseous injection model is therefore used. This model directly injects inside the cylinder gaseous fuel excluding from the computational domain the injector nozzle. This allows to avoid to model the complex phenomena that happens inside the injector, like the phase change from liquid to super-critical. The model approach is based on the work of Baratta et al. [
19,
20].
A set of source cells in correspondence of the injector hole is defined to directly introduce the gaseous fuel inside the cylinder domain based on an imposed injection profile. The complex thermo-physical phenomena that occur in the injector pipe are also neglected. For the injector cell set, an additional source term proportional to the mass flow rate is added for the finite-volume equations of density, turbulent kinetic energy, velocity and mixture fraction to take into account the gaseous injection:
where
is the injection mass flow rate and
is a generic transported property [
19,
20]. The injected fuel is considered as an ideal gas. This choice has been made since this approximation has a small influence on the mixture characteristics far from the injector, where the combustion reactions take place. A better description of the near nozzle processes will be a matter of investigation in a future work. So far, the adopted approach seems to provide a good estimation of experimental data.
3.2.3. Turbulence Model
Turbulence plays an important role in the combustion process of ICE. The model used in this work is the well known
[
21]. It is the most common turbulence model and can provide a stable and accurate description of the effects of turbulence on the average fluid motion, since it solves RANS equations.
3.2.4. Mesh Generator
The mesh is fully automatically computed by Python-based software developed by the Politecnico di Milano–ICE group that is able to provide a spray-oriented and fully hexahedral mesh. The procedure is based on a angular extrusion from a base 2D mesh composed of three regions: (1) the spray region, which is the region where the spray mainly evolves; (2) the layer region, which is the one involved in the layer addition process in the moving mesh; and (3) the piston region [
10]. The regions are automatically computed based on the piston geometry, the spray angle and the injector hole position. This meshing procedure relies on the definition of some constants, which allows the meshing to be tuned to different geometries. However, when the piston geometry and the spray angle change in a wide range, it is difficult to achieve a good mesh in the extremes of the domain.
Table 2 summarizes the information regarding the models used.
3.2.5. Validation of the CFD Model
Two relevant engine operating conditions were selected to validate the the numerical model against experimental data.
Table 3 details the operating conditions of the three cases: the case named
a represents the low load point; meanwhile, the case
b is a medium load case. Regarding the fuel temperature
, in the CFD simulations, it is related to the temperature at the exit of the injector nozzle. Meanwhile, in the experiments, it corresponds to the value measured with the thermocouple after the heating resistance, immediately before the injector.
The comparisons between the simulated and experimental in-cylinder pressure and heat release rate are shown in
Figure 5. Computed traces are in good agreement with the experimental ones. It is of particular interest how the simulations are able to correctly predict the ignition delay and the peak of HRR; therefore, the model captures the evolution of the combustion reactions and the fuel distribution. The chemical mechanism is characterized by cool flames around −25 CAD, and this intensity is overestimated with respect to the experimental results.
The main effect of the direct injection is to create a rich region in the piston bowl. Here, the richer mixture is the first to ignite, and the combustion diffuses towards the rest of the combustion chamber, as represented in
Figure 6.
The direct injection does not reach the squish region. Here, the resultant homogeneity of the mixture could promote a knocking-like combustion in this region when the load increases, since this area tend to ignite simultaneously. On the other hand, close to TDC, the influence of the thermal boundary layer and the heat exchange with the walls are relevant. These properties entail a lower unburned gas temperature that hinders the combustion promoting
and
emissions close to the cylinder wall (
Figure 6).
It is now clear that the diesel-like piston geometry is not optimized for this kind of premixed combustion, since the direct injection only reaches the volume inside the piston bowl.
An optimized piston geometry for RCCI technology has been proposed in [
22,
23]. The authors observed that hydrocarbon emissions are reduced by a piston design with a small squish distance and large bowl volume. Then, the next step of this study is to find a suitable combustion chamber architecture that maximizes the performance of the TCRCI concept.
3.2.6. Set-Up Difference between Validation and Optimization Cases
The set-up used for the optimization simulations differs from the validation case due to the necessity to speed up the computational time and to automatically mesh a large variety of piston geometries. It follows that a fully automatized mesh and coarser time steps are used.
In particular, the time step was doubled during the injection, and it is increased by five times between the end of injection and the start of combustion (
Table 4).
Figure 7 shows the changes in the mesh configuration. The image at the top of
Figure 7 represents a 60 degree sector, and it is manually adjusted to achieve the best possible configuration. The grid size was selected after a preliminary mesh independence analysis based on the previous experience of the authors [
5,
6]. The domain has around 65,160 cells at TDC position [
24,
25]. The bottom image shows the automatic mesh generated to perform the optimization process in order to reduce the time required for each simulation.
These changes do not significantly affect the simulation results, as represented by the comparison of the traces shown in
Figure 8. The only small discrepancy is in the ID. This difference has been considered acceptable since the combustion is chemically driven, and therefore the combustion is sensitive to small local differences. On the other hand, the simulation configuration adopted for the PSO reduces almost half of the computational time required, which was reduced from 40 h to 24 h when decomposing the simulation in eight cores.
3.3. Particle Swarm Optimization
PSO is a population-based optimization technique that iteratively optimizes a candidate based on their previous information and that of its neighbors with regard to a given measure of quality. This algorithm is defined as meta-heuristic since it does not require assumptions of the problem and it is able to search very large spaces of candidate solutions. The mathematical formulation is presented here:
The minimization problem consists of
The D-dimensional domain of
f is the search space. The search space is mapped by a fixed number of position vectors called particles. The set of particles takes the name of swarm
X:
For each t-iteration, the i-particle position is updated based on the velocity vector
:
where the velocity vector is computed based on the relative position of the particle with respect to the current best position of
:
, and the global best positions among all the swarm particles
g:
where
w is the inertia weight,
is the individual weight, and
is the social weight.
and
are random vectors uniformly distributed between 0 and 1. The used algorithm differs from this basic formulation since it divides the swarm into two subgroups that differ for the definition of the velocity: conquerors and explorers particles. The conquerors follow the formulation described above and tend to consistently converge to a minimum. The explorers look into the unexplored or less attractive region of the search space thanks to the application of the Novelty Search concept. For a deeper explanation of the implementation of the algorithm and its performance analysis, please refer to the work of Rodrıguez et al. [
26]; from which the small explanation reported above is taken. The same algorithm has been successfully applied to the optimization of a CI engine operating with OME fuel [
24].
3.5. Definition of the Parameters to Optimize
The parameters to optimize are presented in this section. There are nine parameters in total, and they have been selected due to the high influence on the engine performance:
Hardware: two geometrical parameters to define the bowl profile, and one that specifies the spray included angle.
Air management: pressure and temperature at IVC and Exhaust Gas Recirculation (EGR) level.
Fuel management: the distribution between direct and port injected fuel, and the direct injected fuel temperature.
Some observations regarding the parameters are proposed below:
The pressure and temperature are changed independently; therefore, the total in-cylinder mass and the equivalence ratio vary. Then, it is necessary to update the mixture fraction at IVC for each case to maintain the premixed fuel mass constant. This procedure is made according to the ideal gas equation.
The distribution of direct and port-injected fuel affects different fields. It influences the mixture fraction at IVC:
It also affects the injected mass; since the injector is gaseous and is chocked in the ideal gas hypothesis, the steady mass flow rate is kept constant, and the injection duration is adapted to the total mass. As a consequence, the time-step division is adapted to the new injection duration.
The EGR variation is controlled by changing the chemical table. The reason is that the EGR level affects the environment composition, which is initialized in the table inputs. The EGR level discretization step is 2% to reduce the number of total tables.
The search space for each parameter is summarized in
Table 5:
The flexibility of the geometric parameters and the spray angle are limited by the meshing process because it is difficult to achieve an acceptable mesh quality for wide piston bowls and narrow spray inclusive angles. The thermodynamic conditions at IVC are centered on the baseline values, and the range of temperature might seem small, but the sensitivity of this type of LTC technology is large, as the combustion is chemically driven. The upper limit of is close to the baseline since the injected fuel mass is already very small and the injector is already working almost in the ballistic region. It is thus technically difficult to further decrease the direct injected mass. The upper limit of SOI has been chosen to produce an effective fuel stratification. The bottom limit of the fuel temperature has been chosen to be 20 K above the critical point to be sure to have a super-critical injection.
3.6. Output and Merit Function
This optimization is based on maximizing the efficiency and restricting the maximum pressure and pressure rate. The merit function is configured to favor an efficiency increase; meanwhile, it limits the maximum pressure rate. The merit function is formulated considering a different weight for each parameter. In this case, the efficiency is most important, while the pressure indicators are used only as restrictions.
3.6.1. Efficiency
Since the boosting pressure is different for each tested particle, it is necessary to consider the pumping losses. A loss parameter is introduced to account for the energy required to heat the fuel, since both the direct injected fuel mass and temperature change.
The Gross Indicated Efficiency (GIE) is computed directly from the pressure integration in the volume using the rectangle rule:
The pumping losses are computed using a model based on the turbocharging equation [
1] as follows:
This coefficient ranges almost linearly with the pressure at IVC within three and five percentage points. The constants of the model have been tuned based on the experimental pressure traces of an operating condition very close to the baseline.
Moreover, the amount of energy required to heat the fuel is also considered in a coefficient. This parameter is a function of the fraction of direct injected fuel, the thermo-physical properties and the final fuel temperature. Note that this formulation does not take into account the losses in the heating process. The equation for the efficiency losses due to the fuel heating is the following:
Figure 10 presents a map of this parameter for iso-octane in the range of interest for this application. The heat is computed with CoolProp [
28] considering a constant pressure equal to the injection pressure (300 bar) and a base fuel temperature equal to 380 K. It is possible to notice that the maximum loss in the considered range is 1.24%. It is worth noticing that the energy required could be reduced by exploiting the residual thermal energy of the exhaust gasses.
It must be noticed that pumping losses are one order of magnitude bigger than the fuel heating losses.
emissions are not considered in the merit function since its evolution is strictly correlated to the efficiency, as is possible to notice in
Figure 11.
3.6.2. Maximum Pressure Rate
Knock is the effect of an improper combustion phenomenon. It produces pressure waves that travel inside the cylinder and are reflected by the walls, which result in an oscillatory behavior of the pressure trace. The combustion model used in this work is not able to directly model the detonation, and an index proportional to the knock propensity must be therefore introduced. One of the most often used is the Ringing Index proposed by Eng et al. [
29], which is the result of a theoretical analysis on the pressure waves.
where
is a scale factor determined from the experimental data. In [
30],
is directly used as limiting parameter, since the scaling factors
,
are used to include engines with different characteristics. These parameters have a small influence for similar operating points, since they have the same trend. It follows that the pressure rate is directly used as an index. Furthermore, it has been chosen to move from the time to the angular coordinate, since the engine speed is constant.
3.6.3. Merit Function
The performance of the new combustion system obtained from each case of the optimization process is evaluated by means of the merit function that considers the efficiency, maximum pressure and maximum pressure rate simultaneously. The merit function equation is described as follows:
where
It can be noticed that the efficiency function has a different formulation with respect to the other function: the logarithm is used to increase the relevance of this parameter. For the same reason, the efficiency has a coefficient much greater than the other ones. The limit efficiency is the baseline one; it has been chosen to give an upper tolerance of with respect to the baseline condition to the limit of pressure and pressure rate. The reason is that the baseline condition is far from the limits of the engine (24 bar/CAD and 240 bar). As a consequence, even if it is better to control the maximum pressure rate, a value slightly larger than the baseline is acceptable.
3.7. Baseline Condition
This optimization is made in a relatively high load and high speed condition, since the high load conditions are the most critical for LTC engines, with the charge being premixed and prone to knock at a high in-cylinder temperature and with richer mixtures. The peculiarity of this optimization is that the baseline condition is not directly the case with the diesel hardware but is the result of the previous intermediate optimization. This choice has been imposed by the original spray angle (120°), which limits the mesh quality when the bowl diameter increases. The purpose of the intermediate optimization is to find an intermediate hardware that allows a larger flexibility. This optimization involves only the piston geometry, the spray angle and a small flexibility to adjust the combustion phasing with the variation of the mixture distribution.