A one-dimensional model for the motor oil-fuel dilution under gasoline engine boundary conditions

Nowadays the climate change caused by the green-house effect increasing is a world-wide issue with which scientists have to face. Being the passenger vehicles filled with fossil fuels, the emission of the carbon dioxide green-house gas is unavoidable. In order to slowdown both the global warming and the fossil fuels wasting, lower fuel consumptions, thus high efficiency engines, are required. Currently, the coupled use of downsizing and direct injection in spark ignited engines meets both high efficiency and power-demand requirements. Being downsized engines more compact, the distance between the injector and the engine walls are shorter. Thus, depending on the engine operating condition, the fuel spray wall impingement may occur leading to the formation of a liquid fuel film. If the wall impingement occurs against the cylinder liner wall, which is wetted by a thin oil layer, the motor oil and the landed fuel dilute until the piston arrival. At this time, the mass transport by diffusion have promoted the creation of an oil-fuel mixture. Thus, while thicker part of the mixture layer is scraped into the piston top land crevice, the thinner, whose properties are degraded by the fuel contamination, remains on the cylinder liner. The mixture accumulated inside the crevice may be scattered due to the piston inertia leading to the oil detachment and transport in the combustion chamber. The latter may cause both abnormal combustions named preignitions at high loads and increased particulate emissions at low loads, cold states, and cat-heating. This paper reports the implementation of a onedimensional numerical model for the dilution between fuel and motor oil in a stroke of a gasoline engine. The model aims to give the composition of the oil-fuel mixture both on the cylinder liner after the piston arrival and in the piston top land crevice.


Introduction
Currently, the increasing public awareness about the impact of the passenger vehicles on global warming and air pollution is leading to ever tightening emission limits of green-house gases (e.g. carbon dioxide (CO2)) and pollutants (e.g. Nitrogen Oxides (NOx), soot). Downsized Turbocharged Gasoline Direct Injection (TGDI) engines have shown their potential to comply with both the most recent emissions regulation and the need of fuel economy against the classical Port Fuel Injection (PFI) technology. Since the smaller displacement of downsized engines impacts on power demand and driveability requirements, the support of other strategies such as Direct Injection (DI), turbocharging and Variable Valve Timing (VVT, e.g. Miller cycle) is mandatory.
One of the main drawbacks of coupling DI with downsizing is the small distance between the injected liquid fuel spray and the engine walls with respect to the time needed by the fuel to fully evaporate. This promotes the liquid fuel wall impingement, likely resulting in liquid fuel film formation. When the liquid film is formed on the cylinder wall, which is wetted by a thin oil layer, the oil-fuel overlap leads to the local film thickening and to the oil-fuel dilution until the piston arrival. Due to the large difference between the typical viscosities of oil and gasoline, oil-fuel dilution leads to the oil viscosity degradation, which may cause higher frictions, mechanical losses, and parts wear [1]. At the piston arrival during compression stroke, the thicker part of the non-evaporated film is scraped by the piston top ring into the piston top land crevice [2]. Then, the liquid accumulated into the top land crevice may be scattered into the combustion chamber under the form of flying liquid droplets due to the piston inertia at the end of compression stroke [3]. Furthermore, oil droplets or particles may be carried from the top land crevice to the combustion chamber by the gas upward flow during the blowdown phase [4]. These droplets may burn according to diffusive flames and, since the motor oil is mainly composed by hydrocarbons, they increased the emission of organic Particulate Mass (PM) and Particulate Number (PN) at low load [5]. Furthermore, since the motor oil is characterized by a higher chemical reactivity (lower ignition delay) than that of the fuel, the flying liquid droplets may self-ignite before the spark ignition triggering mega knock events at high load and low speed known as Low Speed Pre-Ignition (LSPI). As a consequence, severe engine damages such as the piston head crack and the increased spark plug wear may result with concerns on the engine lifetime and reliability.
In [6] Luo et al. discussed the difficulties to clean the combustion chamber from the flying oil droplets and proposed the control and the dump of the LSPI intensity as a strategy to operate the engine mega knock free without eliminating the source of these abnormal events. The Authors conducted experiments on an in-line 4-cylinder TGDI engine with particular focus on the test of the mixture enrichment and the coolant heat up that showed to be promising methods to mitigate mega knocks. In [7] Kassai et al. conducted an extensive experimental campaign on a 4-cylinder TGDI prototype engine by implementing a Design Of Experiments (DOE) that included the injection timing, the coolant temperature, the motor oil components and additives, the fuel components. As a result, the Authors reported that lower coolant temperatures, retarded Start Of Injection (SOI), calcium and sodium additives, heavier fuel, promoted the LSPI frequency increase. Higher coolant temperatures and earlier SOI were found to enhance the fuel evaporation leading to a reduced fuel mass available for oil-fuel dilution whilst lighter fuels contributed to slowdown the diffusion between oil and fuel. Zinc and molybdenum additives were found to reduce the LSPI frequency.
At now, the motor oil-induced increase of LSPI frequency and particulate emission is one of the main causes that slowdown the further downsizing of TGDI engines. Therefore, the investigation of oil-fuel dilution under TGDI conditions plays a key role to allow the next advances in the development of high efficiency engines.
In order to avoid the high time and cost spent of experimental campaigns involving LSPI, numerical simulation is an attractive and affordable method to investigate oil-fuel dilution under engine conditions. In the current literature a few authors have reported significant contributions on this topic. Zöbinger and Lauer [8] studied the influence of the oil presence into a fuel droplet under engine environmental conditions by coupling a zerodimensional model of a distillation facility step with a three-dimensional Computational Fluid Dynamics (CFD) simulation. Both the oil and the fuel were modelled as multicomponent hydrocarbon liquids, whose composition was defined respectively according to experimental GC-analysis (oil) and the zero-dimensional distillation model (fuel). The Authors provided the modelled oil and fuel to a CFD code as components of a single droplet to simulate their evaporation under different droplet detachment timing and oil amount. Yu and Min [9] presented a mathematical model for the vapor fuel absorption and desorption and the liquid fuel dilution in the oil film with the aim to model the contribution of these phenomena on the emission of Unburnt Hydrocarbons (UHCs). In the dilution model, the fuel and the oil were modelled as pseudo-pure liquids at uniform temperature available for mass transport by diffusion and surface evaporation. Zhang et al. [10], based on the work of Yu and Min, developed a model to assess the oil available to be scraped and the fuel available to dilute the oil film at the piston arrival after the wall impingement. The Authors included the solid wall modelling, the heat transfer and modelled the fuel and the oil as multicomponent liquids. In [10], based on accuracy-steadiness-computing cost compromises, the Authors proposed a fixed-space grid with adaptive time-step and space-fixed mixture properties calculated with the mass-weighted averaging.
This paper deals with the development in Python 2.7 of a one-dimensional numerical code that models the oil-fuel dilution process and capable to predict the oil viscosity degradation and the maximum oil amount that may be transported into the combustion chamber. The model was implemented by the present Authors with the aim to provide a fast tool to support the very early stage of the downsized TGDI engine development. The code can accomplish this target by performing parametric analysis to compare different engine configurations under review (e.g. spray pattern, injection timing, injection temperature, air-fuel ratio, wall temperature, oil/fuel species, oxygenates blending, compression ratio etc.) and ranking them based on the viscosity reduction and the maximum oil available to be scraped at the end of the simulations. Implementing the code, the mathematical model of Zhang et al. was used as benchmark and then improved with the use of more accurate sub-models, mixing rules and a varying grid size approach. After the description of the main engineering considerations, equations and sets, the code is tested with literature engine data under LSPI conditions to show and discuss the main results achievable with this implementation.

Methodology
The present paper is focused on the modeling and the numerical implementation of oil-fuel dilution for engine application. As far as the oil-fuel film is exposed to the engine cylinder wall on the one side and to the in-cylinder gas on the other side, the model accounts for: a) heat exchange with the solid wall; b) mass diffusion between liquids (oil and fuel); c) heat and mass exchange (evaporation only) with the gaseous phase. The solution domain comprises the liquid film and the solid wall whilst the liquid phase change due to evaporation is approached by using a moving boundary at the liquid end side according to [10].

Figure 1. Scheme of the computation domain
A one-dimensional modelling along the thickness direction is performed to discretize the 2-layers (solid-liquid) solution domain (Fig. (1)). For sake of simplicity the dispersed phase (droplets) was not modeled neither considered, thus, the model is valid until both oil and fuel are in their continuum phase.

Cylinder wall modelling
The cylinder wall is represented by a solid layer hS thick, which is thinner than the cylinder liner height according to different computation tests performed by the present Authors that highlighted the limited heat penetration depth in comparison with the cylinder liner full height. It must be underlined that reducing the solid layer size to the non-varying temperature height (hS) allows faster computations. The cylinder wall modelling includes the conductive heat transfer through the solid wall and the Boundary Conditions (BCs) at the coolant side and at the Soli-Liquid Interface (SLI), which is assumed as a flat contact surface. Conduction is modelled according to the classical Fourier's equation considering no heat source and constant solid thermo-physical properties in both space and time. Due to both the high thermal capacity of the coolant and the high thermal inertia of the solid in comparison with the computation time (dozens ms considering low engine speeds) a Dirichlet BC on the temperature (Eq. (1)) is adopted at the coolant side. Since gravity and shear forces in the liquid phase are not significant due to the typical wall film mass and area, a zero liquid velocity is considered at the SLI neglecting convection. Therefore, being conduction the leading heat transfer mechanism at the SLI, the heat flux and temperature continuity (Eq. (2)) is adopted as BC.

Liquid film modelling
The analysis of the liquid film is conducted under the following assumptions: a) thin film; b) no viscous dissipations; c) oil-fuel dilute mixture; d) pseudo-pure liquids. The mass transfer between oil and fuel in the liquid film is modelled according to the classical 2 nd Fick's diffusion law. The original formulation, which returns the molar concentration mol/m 3 , is rearranged by dividing by the molecular weight both the right and the left hand side of the equation to return the mass concentration kg/m 3 . Once the fuel concentration is returned by the Fick's law, the oil concentration is automatically determined by Eq. (3). Eq. (3) is achieved by writing the volume of the cell as the sum of the oil and the fuel volumes and then expressing each volume as function of the mass concentration, e.g.
After a deep literature review and several tests on solute-solvent combinations close to the oil-fuel mixture (e.g. hydrocarbons with carbons number between C14-C16 as representative of the oil and C6-C8 representative of the gasoline), the estimation of the liquid phase diffusion coefficient has been approached by using the empirical correlation of Siddiqi and Lucas [11] for organic solutions (Eq. (4)). In Eq. (4) the molar volume at the normal boiling point (@NBP) is estimated as the ratio between the molecular weight and the density at the normal boiling temperature with zero vapor fraction.  As shown in Tab. (1), Eq. (4) shows better agreements, i.e. errors around 5%, with respect to the experimental data from [12] against the Wilke and Chang [13] correlation (adopted in [10]), i.e. errors > 10%.
Concerning the mass transfer BCs, a no slip and no penetration (Eq. (5)) condition is adopted at the SLI whilst a Robin condition is adopted at the Liquid-Gas Interface (LGI) to account for the convective mass transfer (Eq. (6)). The driving force of the mass flux is the fuel vapor concentration difference between the LGI and the gas free stream that is estimated by considering thermodynamic equilibrium at the LGI in the perfect gas law (Eq. (7)).
The Mass Transfer Coefficient (MTC) in Eq. (6) is estimated according to the Chilton-Colburn j-factor analogy [14] (Eq. (8)) for heat, mass and momentum transfer that is a powerful method to estimate one of the transfer coefficients by knowing one of the other two transfer mechanisms and the fluid properties. Thus, considering the so called Chilton-Colburn factors for heat (jH) and mass (jM) transfer and assuming to know the heat transfer mechanism, the MTC is given by Eq. (9), where the Lewis number (Le) is defined as the ratio between the gas thermal diffusivity and the gas phase diffusion coefficient between fuel vapor and air. The gas thermal diffusivity is calculated for air as the ratio between thermal conductivity, density and specific heat whilst the gas phase diffusion coefficient is estimated with the Fuller et al. [15] correlation (Eq. (10)) where v is the so called diffusion volume that is calculated with the tables provided in [15]. 2 Since the convective mass transfer is assumed to be the leading mechanism for the liquid film consumption, once the MTC is determined, the liquid film shrinking is calculated as in Eq. (11) and then it is subtracted by the instantaneous liquid layer height hL. L mt    = (11) The heat transfer through the liquid layer is modelled with the classical Fourier's equation considering no heat source and constant liquid thermo-physical properties in time. The SLI BC is consistent with the one adopted at the solid (Eq. (2)). In order to account for the convective heat transfer, a Robin BC is adopted at the LGI (Eq. (12)), where the free stream temperature (T∞) is the gas temperature reduced by the fuel latent heat of vaporization induced cooling i.e. ̇• ⁄ .

( )
Due to diffusion, each point in the liquid layer may contain both oil and fuel with different concentrations, thus, the calculation of the mixture properties in each point is needed. The mixing rules adopted for the oil-fuel binary liquid mixture are the followings: a) mass averaging for the density and the specific heat (Eq. (13)); b) the Filipov law [16] for thermal conductivity (Eq. (14)); c) the simplified Grunberg-Nissan law [17] for the dynamic viscosity (Eq. (15)). Then, the thermal diffusivity in each point is calculated as the ratio between the thermal conductivity (Eq. (14)), the density and the specific heat (Eq. (13)). ,,

Time and space discretization
In order to capture the liquid phase diffusion phenomena occurring in heights of few dozens of microns (typical for oil-fuel film under TGDI conditions) while optimizing the simulation time, a base grid-size (Δx) of 0.1 µm was achieved. Once the base grid-size was determined, the code was tested with different time-steps (Δt) with the aim to determine the best compromise between computing cost and numerical stability. Since several tests highlighted temperature divergence in the liquid layer close to the SLI and to the initial oilfuel interface, the Courant-Friedrichs-Lewy (CFL) stability condition for the heat equation is used to automatically calculate the fixed maximum stable time-step (Eq. (16)). As mentioned before, due to liquids diffusion each cell may contain a binary oil-fuel mixture, thus, the thermal diffusivity varies cell by cell. In order to ensure the tightest stability condition, Eq. (16) is written considering the liquid with the lower thermal diffusivity (the oil in this case). Comparing the solid height (order of millimeters) with the liquid height (order dozens of microns) and considering that in the solid wall do not occur diffusion phenomena of interest, the grid-size along the solid layer was coarsened up to 0.1 mm in order to speed up the simulations while the liquid layer time-step (Eq. (16)) was maintained constant for sake of simplicity.

Liquid-Gas interface evaporation
As mentioned before, the evaporation at the LGI is approached with a moving boundary that is implemented by removing the liquid end cell from the solution domain while retracting the LGI boundary to the face of the neighbour cell. Once the shrinking is calculated (Eq. (11)) and subtracted by the cell size, the liquid end cell is removed or not depending on the remaining liquid height in comparison with the full cell size. Following the fixed-space grid approach proposed in [10], given the mass loss due to evaporation while the solution volume (i.e. Δx) remains fixed, the mass concentration in the liquid end cell would be underestimated until the cell removal.
Applying varying-cell size at the liquid end cell, as the shrinking is calculated, the liquid end cell shrinks in its turn until the removal. Nevertheless, it must be underlined that the higher is the shrinking before the cell removal, the thinner becomes the liquid end cell. Thus, as the evaporation proceeds, the liquid end cell may collapse with concerns on the simulation steadiness.
In order to accomplish both the requirements of accurate and steady simulations, the present Authors proposes an evaporation threshold criterion together with the cell resizing approach (Fig. (2)). Named Φ the ratio between the shrinking and the cell size for a given time-step, the evaporation threshold is defined as the limit of Φ to be exceeded to assume the liquid end cell as fully evaporated. At this point, before reaching the evaporation threshold, the liquid end cell is maintained with reduced size every time-step according to the film shrinking (Δ − ) and its concentration is adjusted as • Δ /(Δ − ). Once the evaporation threshold is exceeded, the liquid end cell is removed. It may happen that the calculated shrinking has exceeded the threshold without exceeding the cell size. In order to comply with both the mass conservation and the evaporation rate calculation, in this case the liquid end cell is removed and the liquid height that is out of the threshold is merged in the neighbour cell (the new liquid end cell) by increasing its size (2Δ − ) and by adjusting its concentration as • Δ /(2Δ − ). It must be underlined that an early cell removal criterion is implemented for the cells where the solution has led to temperatures higher than the fuel critical temperature.

Validation test case
Since the validation of the implemented code through quantitative data is not viable because of the lack of experimental measures of the oil mass scraped by the piston, a qualitative validation is performed with the attempt to capture the main experimental evidences observed in previous works on LSPI.
The TGDI engine data in Tab. (2) by Luo et al. [6] are considered as the baseline to setup the simulations. Then, the same engine is tested in two different conditions that have shown to reduce the LSPI frequency: a) the use of a higher oil viscosity grade [18], which is simulated by replacing the SAE 5W-30 properties with the SAE 10W-30 properties; b) the wall heating-up, which is simulated by increasing the temperature of the solid layer from 90 °C to 105 °C as tested in [6]. Focusing on the baseline, the 2-layer solution domain is composed by: i) the cylinder wall, which is modelled with 1 mm (maximum heat penetration depth) of cast iron; ii) the oil layer, which is modelled with 4 µm of the SAE 5W-30 synthetic oil, whose pseudo-pure liquid properties are calculated with the Simscape TM Thermal Liquid block [19]; iii) the RON93 gasoline liquid film, which is modelled with 20 µm of the pseudo-pure liquid resulted from the averaging of a four-component RON93 surrogate (40.7 mol% i-octane, 13 mol% n-heptane, 40.7 mol% toluene, 5.7 mol% 1-pentene) from [20], whose properties are calculated according to the temperature-dependent-correlations in [21] and with REFPROP 10.0 [22] by giving the name and the mass fraction of each of the four components (saturation pressure).
The HTC is estimated according to the classical Woschni's correlation ( [23]) in Eq. (17), where the recommended values of the coefficients C1 and C2 are provided in [23] depending on the engine stroke (i.e. intake/exhaust, compression or combustion).
According to these fluid data, the time step is about 4E-8 s and the computation grid is made by 250 cells.
An in-house developed zero-dimensional engine thermodynamic two-zone model provides the engine in-cylinder time-varying conditions volume, pressure and temperature needed by both the Woschni's correlation and the BCs at the LGI. The engine data and the fuel data were put inside the zero-dimensional model, which requires the tuning of the cylinder filling coefficient, the Wiebe function and the non-dimensional turbulence (u'/up) at the ignition. The cylinder filling coefficient was calculated by rearranging the Brake Mean Effective Pressure (BMEP) equation to meet the baseline experimental target of 22 bar. The Wiebe function and the non-dimensional turbulence were tuned to capture the experimental fired pressure signal. Once the tuning was performed, the match shown in Fig.  (3) was achieved.

Discussion
For sake of illustration, Fig. (4) shows the oil-fuel dilution (top) and the liquid heating up (bottom) predicted by the dilution model for the baseline case. In Fig. (4) (up) is visible that as the time goes on, the penetration through thickness of both the fuel in the oil layer (x from 0 to 4 µm) and the oil in the non-evaporated fuel film (x from 4 µm to the maximum film height) increases. The profile traced by the iso-CA curves progressive stops represents the time-evolution of the liquid height due to the evaporation induced shrinking. According to Fig. (4) (up), at the piston arrival the very first microns of the remaining oil film are diluted with 30% of fuel while around the same contamination can be observed in the very first microns of the non-evaporated fuel film. Furthermore, at the piston arrival the model predicts an oil presence under the 10% at the LGI and an oil dilution around the 10% at the cylinder wall. Fig. (4) (bottom) shows that, due to the heat transfer with the heated-up liquid fuel film, the mean temperature of the oil layer increased from 90 °C (simulation start, oil considered as fully warm) to 112 °C at the piston arrival. This, together with the oil dilution with the fuel, leads to the further oil viscosity reduction, which enhance the oil scrap and squeeze caused by the piston upward.

Fuel film
Oil layer and the mixing between the two fluids due to the lower viscosities at higher temperatures, the oil scrap is significantly lower because of the faster evaporation time with respect to the diffusion time-scale, as reported in [6]. Moreover, consistently with the small reduction of the average LSPI events observed in [18] when the oil grade turns from 0W-30 to 10W-30, the predicted oil scraps with the SAE 5W-30 (black curve) oil and the SAE 10W-30 oil (green curve, properties from [24]) are comparable.  As visible in the red curves in Fig. (6), under the wall heating up solution the SAE 5W-30 oil suffers from a significant viscosity reduction due to the higher contamination by fuel (18% vol) with concerns on the increase of the oil scrap and squeeze and the parts wear. Indeed, on the one hand the oil viscosity reduction to values around 10 cP may reduce the friction between the piston and the cylinder wall leading to lower Friction Mean Effective Pressures (FMEPs) [25]. On the other hand, it was proven that the lower is the oil viscosity, the higher is the oil consumption and the piston ring wear with increments of dozens of µg/h [26].
baseline: SAE 5W-30 @90 C set-A: SAE 5W-30 @105 C set-B: SAE 10W-30 @90 C set-C: SAE 10W-30 @105 C According to Fig. (5), while the use of higher oil grades is not effective as the wall heating-up in terms of LSPI frequency, it proven to be a powerful method to slow down the oil layer degradation (Fig. (6), green curve). In Fig. (6) the viscosity decreasing trend of the SAE 10W-30 oil (green curve, reduction around the 2%) is both delayed and weaker than that of the SAE 5W-30 oil (black curve, reduction around the 15%) due to the higher frictions between the molecules in oil with higher grades. Thus, the wall heating-up (90 °C to 105 °C) was implemented on the SAE 10W-30 oil with the aim to explore a compromise solution. As visible from the blue bar (Fig. (5)), the combination of the same wall heating-up with the higher grade SAE 10W-30 oil leads to the significant LSPI risk reduction (lower oil mass scraped) followed by around the same baseline viscosity (Fig. (6), blue curve), resulting in a promising solution to comply with both combustion stability and engine parts saving over multiple cycles.

Conclusion
The present paper deals with the numerical one-dimensional modelling of oil-fuel dilution under TGDI engine conditions. The model is based on two literature findings whose several aspects such as the mixture properties calculation, the diffusion sub-models, the numerical implementations are improved and discussed. A baseline conditions based on engine experimental data from literature was simulated and considered as the LSPI reference case. Then, two different solutions that are commonly used to reduce the LSPI events, namely the wall heating-up and the oil grade increase, are implemented and compared. The model is validated in a qualitative manner since it has proven to be capable to capture some key experimental observations on the reduction and control of the LSPI frequency.