High-speed impact of micron-sized diesel drop trains—Splashing dynamics, secondary droplet formation, and effects of pre-existing film thickness

Surface impingement of drops is prevalent in many industrial applications such as spray cooling, inkjet printing, and fuel injection in internal combustion engines (ICEs). Despite efforts to fully understand drop impingement, there exists a fundamental lack of knowledge with respect to the underlying mechanisms that govern impingement and splashing dynamics. This knowledge gap is especially large for simultaneous and successive impingements of high-speed micrometer-sized fuel drops under ICE relevant conditions.¹ Such extreme impingement conditions are most common in diesel fuel injection, during which micrometer-sized drops impinge at high frequencies with velocities approaching 100 m/s.^2–4 The impact of fuel drops on the cylinder wall and piston head may result in full deposition, splashing (partial deposition), or rebound. The ability to predict the outcomes of fuel drop impingement is vitally important to accurately estimate combustion efficiency and hydrocarbon/particulate matter emission, as deposited fuel leads to increased soot and pollutant emission.^5–8 

Computational fluid dynamics (CFD) tools are commonly used to study fuel injection and predict air–fuel mixing and combustion efficiency. The CFD simulations of fuel spray are generally performed in a Lagrangian–Eulerian (LE) framework where a statistical representation of spray is used, e.g., the discrete droplet modeling (DDM) approach,⁹ due to the disparate length and temporal scales. Since individual drop impingements are not captured, a spray–wall interaction (SWI) sub-model is necessary to predict drop impingement outcomes; for example, how much of the drop would be deposited on the surface. The SWI sub-models require a comprehensive understanding of the drop impingement process to predict the transition from deposition to splashing, splashed mass ratio, and secondary droplet size and velocity distributions. The need for accurate SWI sub-models has motivated many studies on single and multidrop impingement, many of which are summarized by Yarin,¹⁰ Moreira et al.¹¹ and Liang and Mudawar.¹² 

Due to the simplified dynamics, many studies have focused on single drop impingement. One study which has been widely used in SWI sub-models is the work of Mundo et al.¹³ Single drop impingement experiments were performed, where drops impinged onto smooth and rough surfaces under multiple impingement conditions and angles.¹³ The transition to splashing was identified and detailed secondary droplet size and velocity profiles were reported. In more recent works, Zhang et al.¹⁴ investigated the effects of liquid viscosity on splashing for single drop impingement on a dry surface and concluded that increasing viscosity suppresses splashing. A model was proposed to predict the splashing threshold for impacts on smooth surfaces which agrees well with the experiments. Guildenbecher et al.¹⁵ used digital in-line holography to analyze the secondary droplets formed during drop impingement onto a thin liquid film. Detailed secondary droplet size, number, and velocity distributions were reported at various times throughout the splashing process. Studying drop impingement numerically, Nikolopoulos et al.¹⁶ performed simulations of single drop impingement onto a liquid film. Two distinct splashing dynamics were highlighted and secondary droplet distributions were presented. Although many studies have focused on quantifying the secondary droplets formed during splashing, the underlying physics which govern secondary droplet formation are not well understood.

There exists much debate on the instability mechanism that drives lamella breakup during splashing events. One common view attributes secondary droplet formation to Rayleigh–Plateau instability^12,17 which describes the breakup of a cylindrical jet into droplets, where surface tension is the dominant force. This theory has been supported by various numerical and experimental studies of, for example, Rieber and Frohn¹⁸ and Zhang et al.¹⁹ among others. Other mechanisms which have been attributed to lamella breakup are Rayleigh–Taylor instability, Richtmyer–Meshkov instability,¹² and a nonlinear amplification mechanism.²⁰ Villermaux and Bossa²¹ performed an extensive study on a liquid drop impacting a solid target of the same size. The impact process, formation of a liquid sheet, and secondary droplet formation were described in detail. From their analysis,²¹ it was determined a global Rayleigh–Taylor instability mechanism is unlikely during the entire lamella/sheet formation process, while a Rayleigh–Taylor like mechanism governs the rim instability leading to ligament, and eventually secondary droplet formation. In the numerical study by Agbaglah et al.,²² the breakup of liquid rims was simulated, and the contribution of both Rayleigh–Taylor and Rayleigh–Plateau instabilities were quantified. While the study of single drop impact and liquid rim breakup has provided a better understanding of drop impingement dynamics and secondary droplet formation, it has led to single-drop-based SWI sub-models which may not apply to complex sprays.

The knowledge gap pertaining to drop–drop interaction has motivated many studies on multidrop impingement. To investigate the impingement of two drops within close vicinity both Raman et al.²³ and Li et al.²⁴ performed simulations using the lattice Boltzmann method. Both studies characterized the central jet formation due to the impact of neighboring drops onto a liquid film, a phenomenon not observed during single drop impingement. Important factors in jet formation were concluded to be drop spacing, film thickness, and liquid viscosity. Raman et al.²³ also studied the effect of liquid film velocity on the symmetry of jet formation. In the work of Li et al.,²⁴ asymmetry of the central jet was also investigated by varying the phase at which neighboring drops impinge. Ersoy and Eslamian²⁵ performed an impressive array of experiments in which two dyed drops impinged in close vicinity on dry and wet surfaces. By varying the delay between the two drops, three distinct splashing modes were identified. Furthermore, it was determined impingement onto a liquid film creates a larger liquid sheet with increased instability. In an effort to approach the chaotic spray impingement process, Liang et al.²⁶ performed numerical simulations of three drops simultaneously impinging onto a liquid film. This complex impingement was found to generate liquid sheets between impinging drops, leading to splashing below the threshold proposed for single drop impingement. While such studies provide insight into the shortcomings of single-drop-based SWI sub-models, they neglect the effect of subsequent drop impingements.

To understand the effects of repeated drop impingement, numerous studies have focused on successive drop impingement and drop train impingement. Yarin and Weiss²⁰ studied the normal impingement of monodispersed trains of drops using ethanol and an ethanol–water–glycerol mixture. From the experiments, the transition to splashing was determined and quantities such as splashed mass ratio and secondary droplet distributions were quantified as functions of an impingement frequency-based parameter. In the work of Zhang et al.,²⁷ drop train impingement was studied experimentally and numerically, focusing on crater formation and crown propagation, where a new crown propagation model was proposed. More recently, Markt et al.²⁸ performed ethanol drop train simulations following the work of Yarin and Weiss.²⁰ The simulations captured the transition from deposition to splashing and agreed well with the experimental splashed mass ratio.²⁰ Extending upon this work, Markt et al.²⁹ performed further investigations on splashing cases and the formation of secondary ethanol droplets. From the simulations, detailed characterization of secondary droplet size and velocity was reported. In the numerical work of Liang et al.,³⁰ the successive impingement of two drops onto a liquid film was studied. Consistent with the results of Yarin and Weiss,²⁰ splashing was suppressed for the trailing drop at high impingement frequencies due to the film flow established by the impingement of the previous drop. Li et al.³¹ performed drop train experiments under low Weber number conditions focusing on the film fluctuations in splashing and non-splashing regimes.

In ICEs, fuel drops generally impinge on heated surfaces which has motivated multiple studies on the effects of heat transfer during drop train impingement. Richter et al.³² investigated the effect of surface temperature and impingement frequency on secondary droplet formation within isooctane drop trains. In the experimental work of Castanet et al.,³³ secondary droplet size and velocity were characterized for the impingement of water drops onto a heated surface at various impingement angles. In a similar study by Qui et al.,³⁴ the splashing angle and spreading diameter were investigated for the impingement of a water drop train onto heated surfaces. More recently, Zhang et al.³⁵ performed experiments with different arrangements of multiple impinging drop trains. The effect of flow rate and drop train spacing on heat transfer characteristics was quantified from the experiments. Furthermore, the interaction of neighboring drop trains was investigated with respect to drop train spacing.

While undoubtedly important to the understanding of SWI, all aforementioned studies focus on impingement conditions that are vastly different from diesel fuel injection. Still, most SWI sub-models are based on millimeter-sized drops impinging at relatively low velocity, which warrants investigation of their performance for high-speed impingement of micrometer-sized drops. This is supported by the work of Aboud et al.³⁶ who studied the oblique impingement of millimeter-sized drops. They found there was a distinct change in the asymmetry of prompt splash for water drops when the diameter was above 1.7 mm, suggesting studies of millimeter-sized drops may not reveal dynamics of micrometer-sized drops. Zhang et al.³⁷ compared the impingement of millimeter and micrometer-sized drops onto a liquid film using numerical simulations. It was determined that the liquid properties, mainly the dynamic viscosity, play a significant role in splashing for micrometer-sized drops. The splashing threshold was found to increase when compared to millimeter-sized drops, stressing the need to study micrometer-sized drops. Focusing on micrometer-sized drops with high impingement velocity, Guo and Lian³⁸ performed simulations of the oblique impingement of a single water drop onto a liquid film. The effect of impingement angle, film thickness, and drop to film density ratio was studied. Extending upon this work, Guo and Lian³⁹ performed numerical simulations in which two water drops obliquely impinged onto a thin liquid film. Multiple cases were presented in which the drop spacing, impingement angle, and impingement velocity were varied. Splashing asymmetry was noted and qualitative analysis of lamella and secondary droplet characteristics was presented. In the work of Markt et al.,⁴ the impingement of a single train of high-velocity micrometer-sized diesel drops was simulated. Focusing on SWI, the numerical results were compared to the model of O'Rourke and Amsden,⁴⁰ which was found to significantly overpredict the splashed mass ratio. Although CFD tools have been used to investigate specific cases of high-speed micrometer-sized drop impingement, few experimental studies exist. This problem is difficult to study experimentally due to the short length and temporal scales. While the numerical studies provide insight into the splashing process, experimental results are for engine-relevant impingement conditions to provide validation datasets. Since few studies have focused on drop impingement under engine-relevant conditions, further investigations are necessary to develop SWI sub-models which are applicable to fuel spray impingement.

This study focuses on the high-speed impingement of diesel drop trains under engine-relevant impingement conditions for drop size, velocity, and temperature. A train of monodispersed drops is a simplified analog to the complex spray impingement process; however, its investigation can be used to approach spray impingement in a systematic manner. First, experimental and numerical results are presented from which a dynamic contact angle model is proposed for diesel drops impinging upon a stainless steel substrate. Then using the newly developed contact angle model, numerical simulations of diesel drop trains are presented for a highly splashing case. Introducing a pre-existing film, the effect of film thickness on splashing characteristics of drop trains is investigated. Using a robust algorithm, the temporal evolution of secondary droplets is studied, revealing distinct phases of secondary droplet formation. Time-averaged secondary droplet size and velocity distributions are then presented for various impingement frequencies and the secondary droplet velocity angle is reported as a function of time. The detailed analysis of high-speed micrometer-sized drop trains presented here is beneficial to the efforts of developing engine-relevant SWI sub-models for hydrocarbon fuel injection.

To perform numerical simulations an in-house multiphase flow solver⁴¹ was used, which solves the full Navier–Stokes equations in a fully Eulerian framework. Both the liquid drop and surrounding gas are included in the simulations. The fluids are assumed to be Newtonain and the flow immiscible and incompressible, which results in the conservation of mass and momentum, Eqs. (1) and (2), respectively,

where ρ is the density, $U→$ is the velocity vector, P is the pressure, $F→B$ is the body force, e.g., gravity, $F→ST$ is the surface tension force and τ is the stress tensor. The volume-of-fluid (VOF) method is used to track the phase (i.e., liquid and gas) volumes in a sharp manner.⁴² The fluid volumes are denoted using a scalar function, α, which is defined as unity within the liquid and zero outside of the liquid. The fluid volume is then tracked by solving the following transport equation:

Equation (3) is solved using the geometric approach of Youngs.⁴³ This method employs piecewise linear interface calculation (PLIC), in which, the interface is represented by a planar polygon in each computational cell.

The equations are discretized following the finite-volume method, and a staggered grid Marker and Cell (MAC) methodology is employed. The two-step projection method^41,44 is used to solve the discretized Eqs. (1) and (2), and a fast multigrid pre-conditioned Bi-CGSTAB Poisson solver is employed which has been found to significantly reduce simulation runtime.⁴⁵ To avoid numerical deformation of the liquid interface in high-density ratio flows, a consistent scheme is used to perform mass and momentum transport.⁴⁶ The surface tension is implemented numerically following the balanced force algorithm,⁴⁷ in which the pressure discretely balances the surface tension force, reducing spurious currents. The flow solver follows a message passing interface (MPI)-parallel framework, which was shown to have great scalability.⁴⁸ A full description of the governing equations and numerical methods is detailed in Pathak and Raessi,⁴¹ while the accuracy of the interface reconstruction is demonstrated in Pathak and Raessi.⁴² 

In previous studies,^41,42,45 the flow solver has been extensively tested and validated using various multiphase flows. More recently, the flow solver's ability to capture the splashing dynamics of fuel drops has been validated.^4,28,29 In Markt et al.,²⁸ the passive scalar routine of Bothe and Fleckenstien⁴⁹ was adopted, providing the capability to track a specific portion of liquid. The routine was extended to include an unbounded number of passive scalars, allowing each impinging drop to be distinctly tracked.²⁸ Each scalar passively tracks a portion of the same liquid, which allows many scalars to be used without significantly increasing the computational cost. In the current simulations, a passive scalar color function, c_n, is used to “tag” the nth drop in the drop train. The function is advected through the following transport equation:

The passive scalar transport equation is solved using the same geometric VOF methodology^42,43 ensuring that the liquid from each drop is sharply tracked with no mixing. As a result, the passive scalar is always associated with the assigned liquid without diffusing. The passive scalar has no effect on the breakup dynamics, and it is possible for secondary droplets to be composed of multiple passive scalars. In the current study, passive scalars are used when simulating drop train impingement to provide accurate splashed mass calculations and gain further insight into the splashed mass composition, which is extremely difficult to perform experimentally.

Accurate contact angle implementation is necessary to capture the impingement dynamics of a train of drops on an initially dry surface; although drops succeeding the first one in a drop train land on a liquid film formed by preceding drops, the contact angle plays an important role in spreading dynamics of the film. Also, imprecise predictions of film spreading/thickness may affect the splashing dynamics, which will be discussed in Sec. IV B 5. The numerical contact angle implementation builds upon Markt et al.,⁵⁰ where the dynamic contact angle model of Yokoi et al.⁵¹ was validated. In this work, the model of Yokoi et al.⁵¹ is used to capture the spreading of diesel drops after empirically determining the model parameters described below. This model⁵¹ contains capillary-dominated and inertia-dominated regimes for low and high capillary number flows, respectively. In the capillary-dominated regime the instantaneous contact angle is determined from Tanner's law,⁵² shown in Eq. (5),

Here, θ_d is the dynamic contact angle, θ_e is the equilibrium contact angle, k is a material-related constant, and Ca is the capillary number defined as

where μ is the dynamic viscosity, σ is the surface tension and U_CL is the contact line velocity. In the inertia-dominated regime, a maximum or minimum contact angle is prescribed. During the advancing phase, the maximum advancing contact angle (θ_mda) is applied, while in the receding phase the minimum receding contact angle (θ_mdr) is imposed. Therefore, the inertia-dominated regime may be expressed as follows:

By integrating both contact angle regimes, Yokoi et al.⁵¹ obtained the following dynamic contact angle model:

where k_a and k_r are the material-related constants for the advancing and receding phases, respectively. Both parameters are determined empirically for specific liquid–substrate pairs along with θ_mda, θ_mdr, and θ_e. In this study, the model parameters are determined based on experimental results of diesel on stainless steel, presented in the  Appendix.

To quantify the secondary droplet characteristics a robust algorithm29 was used to determine the size, number, velocity, and velocity angle for all secondary droplets in the domain at a given instant. The algorithm, presented and validated in Markt et al.,²⁹ was previously used to characterize the secondary droplets formed from the impingement of ethanol drop trains.²⁹ Scrutinizing the volume fraction field, the algorithm identifies all isolated liquid structures, or secondary droplets. The majority of the secondary droplets are non-spherical; therefore, the volume of each secondary droplet is determined and an equivalent spherical diameter is calculated using the droplet volume. For a detailed description of the secondary droplet characterization algorithm refer to Markt et al.²⁹ 

Experiments of single diesel drop impingement were performed to aid in the development of a dynamic contact angle model for diesel on stainless steel. To repeatedly generate a single drop with a specific diameter, a precision syringe pump was used with a volumetric flow rate of 0.2 ml/min. The drops accelerate due to gravity before impacting a heated stainless steel impingement plate. To capture the impingement process, a high-speed camera was used along with a high intensity light source. To obtain collimated cylindrical light, an analog mode LED lamp was focused by an iris and passed through a plano–convex lense. The light is then captured by a Photron Fastcam SA 1.1 high-speed camera with the use of a Nikon Nikkor lens with f-stop of 4, positioned opposite of the LED lamp. Further details on the experimental setup can be found in Zhao et al.,⁵³ Zhao.⁵⁴ 

To achieve thermal conditions representative of drop impingement in ICEs, the diesel fuel and impingement plate were heated to 150 °C. Two temperature probes were used to ensure isothermal conditions were present during impingement. Each temperature probe consists of two “J” type thermocouples which provide a fast response time. One probe was located on the surface of the impingement plate at the impingement location; a second probe was embedded 2 mm below the surface of the impingement plate at the impingement location. This setup provided the ability to measure the heat flux through the plate to ensure isothermal impingement conditions were achieved. More information on the heat flux measurement apparatus and methodology can be found in Zhao et al.⁵⁵ 

Analyses of the contact line velocity, dynamic contact angle, spreading factor, and height ratio were performed using an in-house MATLAB image processing code. The code analyzes each frame captured by the high-speed camera and subtracts the background from the image. The image is then converted to binary using Otsu's method.⁵⁶ From the binary image, the liquid–gas and liquid–solid interfaces are then identified and used to calculate quantities of interest. By locating the contact line and measuring the angle between the tangent to the liquid–gas interface and the solid–liquid interface, the dynamic contact angle was measured as a function of time. For a detailed description of the image processing algorithm and post-processing method see Zhao et al.,⁵³ Zhao,⁵⁴ Zhao et al.⁵⁵ 

This section presents our computational simulation results. We first begin with the impact of a single diesel drop onto a stainless steel surface, where a dynamic contact angle model is proposed. Then, simulations of the surface impingement of a train of diesel drops are presented.

To acquire the necessary contact angle model parameters and obtain an experimental validation case, the impingement of a single diesel drop onto a stainless steel plate was investigated. Note that heat transfer effects were eliminated as both the drop and the impingement plate were heated to 150 °C representative of ICE fuel spray impingement conditions. The diesel properties, taken from the CONVERGE fuel property database⁵⁷ for diesel #2 at 150 °C, are shown in Table I. In both the experiments and simulations, a drop diameter of D = 2.44 mm and an impingement velocity of V₀ = 1.4 m/s were used. Relying on the experimental results and following the procedure detailed in the  Appendix, the contact angle model parameters relevant to Eq. (8) were obtained and are presented in Table II.

The proposed dynamic contact angle model, presented in Table II, was used to perform a single diesel drop impingement simulation under the impingement conditions mentioned above. The simulation was performed at a resolution of 80 cells per diameter (CPD) of the drop, deemed sufficient based on results from preliminary simulations.⁵⁰ 

In Fig. 1, the simulation (left) and experimental (right) results are presented which agree well qualitatively. A short time after impingement of the drop a thin film is formed while the drop retains its spherical shape as seen in Figs. 1(a) and 1(b). At 1.25 ms after impingement, the thin film continues to spread while the top of the drop becomes an elongated hemisphere. Advancing forward to 4.45 ms after impingement, the drop is nearing its maximum spreading diameter but continues to spread. The spreading regime is clearly evident due to the relatively large apparent contact angle seen in both the simulation and experiment. Finally, at 10.0 ms, the droplet has reached its maximum spreading diameter and relaxation of the contact angle is observed, where the equilibrium contact angle, measured as 10°, is present. Both the simulation and experiments show contact angle relaxation as seen in Figs. 1(g) and 1(h).

To quantitatively assess the accuracy of the proposed dynamic contact angle model, the instantaneous diameter d_spread and height h_drop of the impinged drop were measured to compare the spreading factor (⁠$dspread/D$⁠) and height ratio (⁠$hdrop/D$⁠) between the experiments and simulations, which are presented in Figs. 2(a) and 2(b), respectively. For both quantities, an excellent agreement is seen between the computational and experimental results. Beginning with the spreading factor, the simulation results show a maximum relative error of 2.5%, where numerically the drop spreads slightly faster around 4 ms. At t = 10 ms, the spreading phase is completed and the maximum spreading factor of ∼4.6 is observed for both the simulation and experiment. In general, the dynamic contact angle model is able to accurately capture the drop spreading and maximum spreading diameter very well.

The height ratio from the simulation is measured using two methods. The first one follows the method used in the experiments, where the film height is measured from a side view of the drop as seen in Fig. 1. The simulation result using this method is shown with the solid red line and labeled Measurement 1 in Fig. 2(b). Comparing the simulation and experiment, the height ratio is captured well throughout the impingement with a slight underprediction between 5 and 10 ms. After 5 ms, a saturation in the height ratio is seen both in the experiments and simulation due to the measurement method. Once the rim at the advancing liquid front is formed, the height of the rim remains fairly stable. Therefore, the height ratio remains constant during this time, motivating an alternative approach to measuring drop height. The second method measures the average film thickness across a square of sides equal to the drop diameter (D), which is centered at the impingement point [see inlay in Fig. 2(b)]. The result using the second method, labeled measurement 2, is shown with the dashed line in Fig. 2(b). This method eliminates the saturation of the height ratio by excluding the liquid at the film front. Measuring the height ratio in this way, which may not be straightforward experimentally, provides insights into the instantaneous film thickness for drop train impingement. The film thickness during the impingement of a subsequent drop is an important parameter that affects the splashing of the impinging drop, which will be discussed later.

The agreement between computational and experimental results in the spreading factor and height ratio provides confidence in the ability of the dynamic contact angle model to capture the spreading of diesel drops on stainless steel. It should be noted that in the computational simulations, and by extension the proposed contact angle model, the surface is considered to be ideally smooth. Capturing the effects of surface roughness requires revisions to the contact angle model. Accurately capturing the spreading dynamics of the liquid film is vital to predict quantities such as the splashed liquid mass and the transition from deposition to splashing when simulating the impingement of a train of drops.

We will now present simulations in which a train of monodispersed diesel drops impinge on a stainless steel substrate. Although a train of drops is a simplified representation of fuel spray, the interaction between subsequent drop impingements provides insights into the spray impingement phenomenon. Drop size and impact velocity were adopted from Markt et al.,⁴ in which LE fuel injection simulation results were analyzed to extract a Sauter mean diameter and normal impingement velocity. A representative drop diameter of D = 5.97 μm and a normal impingement velocity of V₀ = 77 m/s were selected for all simulations in the current study. Additional details on the injection conditions and droplet selection can be found in Markt et al.⁴ 

Using the diesel properties presented in Table I, the above drop impact condition corresponds to a Reynolds (Re) and Weber (We) number of Re = 736 and We = 1443. The gas properties are also held constant for all simulations and are ρ = 1.226 kg/m³ and ν = 1.45 × 10⁻⁵ m²/s. Following the formulations of Yarin and Weiss,²⁰ a nondimensional velocity, u, defined in Eq. (9), is used to characterize the impact conditions,

where f is the impingement frequency. Yarin and Weiss²⁰ experimentally determined the transition from non-splashing to splashing, or the splashing threshold, to be between $u=16−18$ (at u = 18 the drop train becomes splashing). While varying impingement frequency has a significant effect on the impingement outcomes, such an investigation is outside the scope of this work and will be addressed in future studies. The current study is only focused on investigating the splashing dynamics for the highly splashing case of u = 24, corresponding to $f=$ 3.103 MHz for the diesel drop at the above representative D and V₀.

The three-dimensional (3D) computational domain is of size $4.19D×4.19D×3.22D$⁠. On the impingement surface at the bottom of the domain a no-slip boundary condition is imposed. The side walls and top surface of the domain allow the secondary droplets to exit the domain. Unless it is mentioned otherwise, all simulations were performed at a resolution of 80 cells per diameter (CPD) of the drop, which will be shown to provide sufficient accuracy. Although the length scales are very small, the continuum flow assumption is still valid based on the estimated mean free path of the liquid. In most simulations, the train consists of 17 drops.

We will begin with a drop train impingement onto an initially dry surface at u = 24. Simulation results are presented in Fig. 3, where subsequent drops are tagged using four repeating passive scalar color functions starting with the second drop: clear, red, yellow, and green, respectively. The first drop is clear too. Figures 3(a)–3(f) show the impingement of the second–seventh drops in the drop train, respectively. A dimensionless time, $T=tf$ is used to characterize the impingements, where t is the physical time measured after the impingement of the first drop. The dimensionless time allows for easy identification of the impinging drop at any time, as a dimensionless time T corresponds to the impingement of the $T+1th$ drop. So, for example, T = 1.055 corresponds to a short time after the impingement of the second drop within the drop train.

For impingement on an initially dry surface, the first drop within the drop train is deposited on the surface creating a liquid film as depicted in Fig. 3(a). All subsequent drops which impinge on the film consistently lead to highly splashing dynamics, characterized by lamella instability and secondary droplet formation as seen in Figs. 3(b)–3(f). The splashing lamellae are from the previously impinged drops (second–sixth drops, respectively). The highly splashing dynamics are representative of the impingement of all subsequent drops in the drop train, suggesting the number of impinged drops does not significantly affect the splashing behavior. Focusing on Fig. 3(d), the lamella present in the image was formed by the impingement of the fourth drop (yellow). Here, defined cusps are seen in the lamella which eject secondary droplets. During this particular impingement, the secondary droplets formed upon splashing are composed of red and yellow liquid, from the third and fourth impinging drops, respectively. The large number of red secondary drops suggests that the splashed mass is composed of liquid from both the impinging drop and film formed from the previous impingement. The same observation holds for subsequent drop impacts, indicating that during splashing the upper portion of the lamella, which produces the first secondary droplets, is composed of film liquid from the previous drop, while the lower portion of the lamella is composed of liquid from the impinging drop itself. With the use of passive scalars, it is possible to quantify the composition of the splashed mass. Once pseudo-steady state impingement dynamics are present (⁠$T≈$ 6), 58% of the splashed liquid from each impingement is composed of the impinging drop. This means 42% of the splashed liquid mass is composed of film liquid deposited from the previous impingement. This unique insight provided by the passive scalars can be used to better understand drop train impingements that do not occur under thermal equilibrium. For example, under the conditions present in ICEs, the film liquid may have significantly different thermo-physical properties, due to heat transfer from the impingement surface. Therefore, the assumption that all secondary droplets are composed of liquid from the impinging drop may lead to inaccuracies in the prediction of secondary droplet vaporization and combustion efficiency. The composition ratio reported above can be used to better approximate the properties of splashed mass, using the film and impinging drop properties. Further studies that include heat transfer and evaporation are needed to understand their effects on splashing.

Since liquid film formation is an expected phenomenon during fuel injection, the effect of a pre-existing liquid film must be better understood. Due to the nature of drop train impingement, a transient film formation period is present, in which the drop train impinges on an initially dry surface. As the first few drops impinge on the surface, they create a liquid film that eventually reaches a pseudo-steady state in terms of thickness and momentum at the center of impact and its vicinity. Therefore, the splashing dynamics experience a transient period before reaching pseudo-steady state. We will use the ratio between the splashed liquid mass (M) and impinged liquid mass (M₀), termed “splashed mass ratio” (⁠$RM=M/M0$⁠), to quantitatively distinguish the transient and pseudo-steady state periods. As will be shown below, when R_M asymptotes, it is an indication that effects of film formation during the transient period have been diminished and that the pseudo-steady state is reached.

Here, we will examine the effects of a pre-existing liquid film on the duration of the transient period leading to the asymptotic splashed mass ratio. The simulation at u = 24 was repeated with a pre-existing film thickness of 0.31 μm. The thickness was determined from the initially dry surface simulation and was measured when the pseudo-steady state dynamics were established. The measurement was done over the footprint of impinging drops in a square of size D, illustrated on the inlay of Fig. 2(b). The film was initialized uniformly and at rest; a naturally developing film would have inertia and may vary in thickness further away from the impingement point. Therefore, a transient period is still expected due to the development of film inertia. To hold the liquid film in the domain, no-penetration boundary conditions were imposed at the bottom of the sidewalls of the domain up to a height that is twice the film thickness.

Figure 4 depicts the results from both the initially dry surface and pre-existing film simulations at a dimensionless time of T = 5.213. The walls holding the liquid film are also visualized in Fig. 4(b). In both cases, shortly after the impingement of the sixth drop, the lamella from the fifth drop (green) breaks into secondary droplets as seen in Fig. 4. The splashing dynamics seen in both simulations are very similar, suggesting the effects of the pre-existing liquid film are minimal; however, we will examine the effects of the initial film thickness later in Sec. IV B 5. In both cases, secondary droplets are formed from the previously impinged third and fourth drops (red and green, respectively).

The splashed mass ratio (R_M) was calculated using a “splashing threshold height,” which differentiates between impinged, pre-impinged and splashed liquid at any instant throughout the simulation. The procedure to calculate the splashed mass and ultimately R_M has been detailed in Markt et al.;²⁸ briefly, it is as follows: First, a passive scalar is used to tag all drops before impingement. As a drop then crosses the splashing threshold height moving toward the surface (before impingement), the liquid below the splashing threshold height is no longer tagged and is considered impinged. The liquid is only considered splashed after crossing the splashing threshold height with an upward velocity after impingement. In the current study, a splashing threshold height of 4.0 μm was adopted from Markt et al.,⁴ where it was verified to be suitable for calculations of the splashed mass ratio for the impingement of micrometer-sized diesel drops. The splashing threshold height is shown in Fig. 5, where all secondary droplets have crossed the splashing threshold height for the highly splashing case at u = 24.

The splashed mass ratio for the initially dry surface and pre-existing film simulations are shown in Fig. 6 as a function of dimensionless time. To calculate the splashed mass ratio, the total splashed mass for all drops (M) and the total impinged mass (M₀) are used, therefore, the splashed mass ratio is the average of all drop impingements. When the impingement occurs on a dry surface the first drop in the drop train is completely deposited, while impingement onto a pre-existing liquid film results in splashing of the first drop. Since the first drop in the initially dry surface simulation is deposited, the splashed mass ratio for the first drop is zero. Therefore, to directly compare the splashed mass ratio between the two cases, the mass of the first drop must be excluded for the dry surface case. In Fig. 6, periodic undulations are observed in the splashed mass ratio, which is expected due to the calculation method. Once the drop crosses the splashing threshold height it is considered impinged. Since the splashed mass ratio is calculated at short temporal scales, there exist times when the drop has impinged but not yet splashed, temporarily reducing the splashed mass ratio. It is not until the impinging drop begins to splash that the splashed mass ratio begins to increase again.

As seen in Fig. 6, the asymptotic splashed mass ratio for the two simulations is within one percent of each other, 58% for the initially dry surface and 57% for the pre-existing film case. The initially dry surface simulation reaches the asymptotic splashed mass ratio at a similar time when compared to the pre-existing film simulation. The pre-existing film case has a greater splashed mass ratio only for the initial drop impingements, which may be attributed to the transient splashing dynamics. The results suggest that the pre-existing film does not affect the asymptotic R_M and the pseudo-steady state dynamics when compared to the initially dry case with a naturally developing film. As we will discuss in Sec. IV B 5, the pre-existing film thickness, however, does affect splashing dynamics during the transient period.

To ensure a resolution of 80 CPD sufficiently resolves the splashing phenomena, additional simulations were performed at u = 24 with a pre-existing film thickness of 0.31 μm at resolutions of 40 and 160 CPD, which correspond to spatial resolutions of $Δx$ = 0.15 μm and $Δx$ = 0.0375 μm, respectively. The splashed mass ratios at the three resolutions are presented in Fig. 7 as a function of dimensionless time, which shows that the results from 80 CPD and 160 CPD simulations agree well, while at 40 CPD the splashed mass is underpredicted. The liquid film at 40 CPD is only resolved by two computational cells, insufficient to capture the film dynamics which directly affect splashing. As a result, the splashed mass ratio at 40 CPD is 7% less than that at 80 CPD after the impingement of the seventh drop. Comparing the splashed mass ratio from the 80 and 160 CPD simulations, less than a 2% difference exists after the impingement of the seventh drop, indicating a second order convergence. The splashed mass ratio is compared after the impingement of the seventh drop due to the extreme computational cost of the 160 CPD simulation. At a resolution of 160 CPD, 400 CPU (central processing unit) cores were required with a runtime of 4 months to obtain results through the impingement of the seventh drop. Based on the convergence seen in the splashed mass ratio, a resolution of 80 CPD was deemed sufficient to capture the splashing dynamics.

Next, we will use the 160 CPD simulation for a detailed view of the lamella instability and breakup process. In Fig. 8, results of the 160 CPD simulation at u = 24 are presented for the impingement of the first through fifth drops of the drop train. Similar to the 80 CPD results the first impinging drop and the pre-existing liquid film are clear, and four passive scalars are used in a repeating scheme for subsequent impingements: clear, red, yellow, and green, respectively.

In Fig. 8(a) the first drop has just impacted the liquid film forming a thin lamella. The first phase of secondary droplet ejection can be seen, as a ring of small secondary droplets that have been ejected from the lamella rim, shown in detail in the figure inlay. This phase of secondary droplet formation produces the smallest droplets when compared to the subsequent phases. As the lamella continues to stretch, instabilities in the lamella rim begin to propagate toward the lower portions of the lamella, where holes in the lamella form a web of ligaments which appear to occur at a specific wavelength. This is depicted in Figs. 8(b) and 8(c) for the impingement of the second and third drops, respectively. A similar formation of holes in the lamella was observed experimentally by Burzynski and Bansmer⁵⁸ who studied drop splashing on thin moving films at high Weber number. As the lamella continues to breakup, the second phase of secondary droplet ejection begins as seen in Fig. 8(c). Here instabilities cause the upper portion of the lamella to form cusps from which secondary droplets pinch off and eject. Many such secondary droplets can be seen in Fig. 8(c). As secondary droplets are ejected from cusps, smaller satellite droplets pinch off from the cusps as reported by Wang and Bourouiba.⁵⁹ The final phase of secondary droplet formation is due to the breakup of the lower lamella rim. As the lamella expands outward from the impingement point it eventually separates from the liquid film, forming a lower lamella rim; this detachment is shown in Fig. 8(e). The lower rim continues to stretch until it breaks into a ring of larger secondary droplets. Just before breakup, the intact lower lamella rim is clearly seen in Fig. 8(f) after the impingement of the fourth drop of the drop train (yellow). The breakup of the lower rim is shown in Fig. 8(d) for the impingement of the third drop (red), marking the final phase of secondary droplet formation which produces the largest secondary droplets. Leveraging the capabilities of the passive scalars, additional details of secondary droplet formation are observed. One interesting dynamic is the composition of the secondary droplets. In Fig. 8(c), the upper portion of the lamella and the first ejected secondary droplets are composed of the film liquid (clear), while the lower portion of the lamella is composed of the third drop (red) which has just impinged. This suggests the film liquid significantly contributes to the secondary droplet mass. The lamella breakup dynamics described above are representative of all subsequent drop impingements.

Although an instability analysis of lamella breakup is outside the scope of this study, a brief comparison with previous studies can be made. One quantity that has been used to identify the splashing mechanisms is the instability frequency and the corresponding number of cusp sites. In the 160 CPD simulation presented in Fig. 8, there are approximately 48 cusps that regularly form during the splashing of each drop. In the work of Marmanis and Thoroddsen,⁶⁰ the impact of liquid drops on a solid surface was studied at high Reynolds and Weber numbers. The number of fingers was determined from the experiments and the following relationship was proposed:

Applying expression 10 to the diesel drop train simulation, one obtains a prediction of 47 cusps, which agrees very well with the 48 cusps observed in the simulation. Villermaux and Bossa²¹ studied single drop impact on a similar sized solid surface and determined a Rayleigh–Taylor like mechanism is responsible for the formation of corrugations in the lamella rim which corresponds to cusp sites. As described by Villermaux and Bossa,²¹ due to deceleration of the rim, there exists a body force, proportional to $−ρR¨$⁠, where $R¨$ is the rim acceleration, pushing the fluid in the rim radially outward forming corrugations. The growth rate, ω, of the rim perturbations was derived within the planar approximation and infinite depth limit (see Bremond and Villermaux⁶¹) for wavenumber k and reads²¹ 

Villermaux and Bossa²¹ determined the optimal wavenumber, k_c, corresponding to the maximum growth rate, can be expressed as

Using k_c and the diameter of the rim at the onset of instability (⁠$≈$26 μm), the predicted number of cusps for the diesel drop train case is 84, which is fairly close (order of magnitude) to the number of cusps in the simulation. The agreement with the above models,^21,60 suggests that a Rayleigh–Taylor like mechanism is a possible cause of cusp formation at the lamella rim of the diesel drop train.

It was previously determined that the introduction of a pre-existing liquid film does not affect the asymptotic splashed mass ratio for the specific case where the pre-existing film thickness is equivalent to the naturally developed film thickness from impingement onto an initially dry surface. To quantify the effect of pre-existing film thickness, two additional drop train simulations were performed at u = 24 where the naturally developed film thickness (0.31 μm) was halved (0.155 μm) and doubled (0.62 μm). In Fig. 9, the splashed mass ratio is presented for the initially dry, 0.155, 0.31, and 0.62 μm pre-existing film cases. After the impingement of the first drop, which is representative of single drop impact on a liquid film, the splashed mass ratio is 15%, 34%, and 95% for the 0.155, 0.31 and 0.62 μm pre-existing film cases, respectively. While at the extremes a difference of 80% exists after the impingement of the first drop, the effect of initial film thickness diminishes as subsequent drops impinge. At the asymptotic values (T > 10), the splashed mass ratios are 53%, 57%, and 62% for the 0.155, 0.31, and 0.62 μm pre-existing film cases, respectively. On the initially dry surface, R_M asymptotes to 57% too.

The near convergence of the asymptotic splashed mass ratio for the three cases suggests the effect of pre-existing film thickness is less significant in drop trains than in single drop impingement. Since the splashed mass ratio is the average of all drop impingements, the later drops must have similar splashing dynamics to reduce the effect of changes in the initial film thickness. The results indicate a similarity in the evolution of the liquid film which will be verified next.

To quantify the evolution of the liquid film, the average film thickness was measured in a cube of length D centered upon the impingement point in all four cases (dry surface and pre-existing films). This control volume is the same shown in the inlay of Fig. 2. In Fig. 10, the average film thickness is presented as a function of dimensionless time for all cases. Before the impingement of the first drop (⁠$T<$ 0), the prescribed pre-existing film thickness can be seen for all cases. Upon impingement of the first drop, the film thickness increases to a maximum value due to the impinging drop itself. After impingement, as the lamella extends radially and the impinging drop is spreading, the film thickness decreases until the impingement of the next drop. Therefore, the film thickness at the time of drop impingement corresponds to the low points in Fig. 10. For all cases, the film thickness becomes nearly identical at the impingement of the second drop (T = 1). By the impingement of the third drop, all cases have the same film thickness under the drop footprint. This result suggests that regardless of the initial film thickness (below $0.62 μ$m), the impinging drop has enough inertia to make the film assume its natural thickness that would arise from impingement onto an initially dry surface. Therefore, by the impingement of the third drop, there is no dependence on the initial film thickness, yielding the same splashing dynamics for all cases. We can then conclude that the initial difference in splashed mass ratio, seen in Fig. 9, is due to the first drop as well as film inertia, which was discussed earlier. Consequently, SWI sub-models which are based on single drop impingements may lead to an overprediction of the splashed mass ratio for dense fuel sprays. A more accurate and robust approach is studying the impingement of multiple (e.g., a train of) drops to capture the interaction between the drops and film.

By analyzing the secondary droplet characteristics at specific times during the impingement of a single drop within the drop train, insights into the lamella breakup process and phases of secondary droplet formation are gained. We will now present the temporal evolution of the secondary droplets formed for the drop train impingement at u = 24 with a pre-existing film thickness of 0.31 μm at resolutions of 80 and 160 CPD; note that for brevity, we did not include figures from the 80 CPD simulation, showing the time history. In Fig. 11 instantaneous secondary droplet distributions are presented for the impingement of the fourth drop, at dimensionless times which correspond to the simulation results shown in Figs. 8(d)–8(f). The histograms in Figs. 11(a), 11(c), and 11(e) show probability density functions (PDFs) of secondary droplet size and the histograms in Figs. 11(b), 11(d), and 11(f) show the percent of secondary droplet volume each bin represents. The secondary droplet diameter, d, is normalized by the impinging drop diameter, D, providing the diameter ratio d/D. In all secondary droplet histograms, the bin size is equal to D/80. Note that the incoming drop is excluded from the secondary droplet characterization as is any liquid structure that resides or extends below the splashing threshold height.

The PDFs of secondary droplet size show little variation through time for both simulation resolutions. The results indicate the majority of secondary droplets are quite small, $d/D=0−0.05$⁠, and the droplets appear to follow an exponential distribution. Comparing the evolution of secondary droplets over time, the number of small secondary droplets increases throughout the splashing process. This trend is seen in the simulations at both 80 and 160 CPD. In general, the PDFs at both resolutions agree quite well, the only significant difference is seen in the smallest bin of secondary droplets. At 160 CPD the largest number of secondary droplets are of size d/D = 0.0125, while at 80 CPD the majority of secondary droplets are of size d/D = 0.025. In the 80 CPD simulation, any drops of size d/D = 0.0125 are only resolved by one computational cell or less. Since the drops are significantly under resolved, their formation may be suppressed leading to the difference in the PDFs between 80 and 160 CPD. In all of the PDFs, the smaller secondary droplets outnumber the larger ones. However, the secondary droplet volume distributions show that the smaller droplets only represent a small percentage of the splashed liquid and the majority of splashed volume consists of the larger droplets.

As described in Sec. IV B 4, three distinct phases of secondary droplet ejection are seen. The histograms in Figs. 11(a) and 11(b) correspond to a time just after the first stage of secondary droplet ejection [see Fig. 8(d)], where a ring of small droplets has detached from the lamella rim. At this time, the lower lamella rim from the previous drop impingement [red drops in Fig. 8(d)] is breaking into secondary droplets which are still connected by small liquid filaments. The signature of the lower lamella rim is seen in Fig. 11(b). Due to its large volume, the rim is registered as a large droplet of size $d/D≈$ 0.4 at both simulation resolutions. At this instant, a large portion of the secondary droplet volume is composed of droplets from the breakup of the previous (red) drop. The size of these drops is between $d/D=0.05−0.18$ as seen in Fig. 11(b).

A short time later at T = 4.065, the second phase of secondary droplet ejection has occurred. Here the upper portion of the lamella has completely broken up as seen in Fig. 8(e) and reflected in Fig. 11(d) where the volume of medium sized secondary droplets (⁠$d/D=0.15−0.2$⁠) is increased. For the simulation at 80 CPD, the presence of a large secondary droplet is still observed. This is because at 160CPD the rim quickly breaks up into smaller secondary droplets as seen in Fig. 8(e), while at 80 CPD the lamella rim first breaks into large ligament structures as seen in Fig. 3(f) for the case on an initially dry surface.

Finally, Fig. 11(f) depicts the distribution of secondary droplets at the end of the repeating breakup cycle, where the absence of larger secondary droplets shows the lamella has fully broken into secondary droplets. At this time secondary droplet formation has concluded, yielding a final size range of $d/D=0.05−0.18$⁠. The secondary droplet distributions shown in Fig. 11 are repeated for subsequent drop impingements. The temporal distribution of secondary droplets in the current work (We = 1443) is comparable to those presented in Markt et al.,²⁹ where a similar analysis was performed on the impingement of ethanol drop trains (We = 382; u = 22). Although the impingement conditions are vastly different in the two studies, the same phases of secondary droplet formation are seen. However, while the nondimensional velocities are close, larger secondary droplets (⁠$d/D>0.2$⁠) were observed in Markt et al.²⁹ That is consistent with the experimental results of Mundo et al.¹³ (⁠$We=580−1182$⁠), where the secondary droplet size was found to decrease with increasing Weber number.

Comparing the 80 and 160 CPD results in Fig. 11, the agreement in the droplet size range is very good. Focusing on the secondary droplet volume, while the distributions vary slightly throughout time [Figs. 11(b) and 11(d)], the distributions at the completion of splashing [Figs. 11(e) and 11(f)] are very similar between the two resolutions. Note that SWI sub-model development is in particular, concerned with the distributions at the end of splashing events, and the good agreement at such time further establishes that a resolution of 80 CPD is sufficient. The time-averaged secondary droplet distributions of 80 CPD simulations are presented next.

While the temporal evolution of secondary droplets, presented in Sec. IV B 6, provides details of the splashing process, the time-averaged analysis describes the outcomes of the drop train impingement as a whole, which are more helpful for SWI sub-model development. Figures 12(a) and 12(b) illustrate time-averaged secondary droplet distributions from the impingement of the above drop train on the initially dry surface and the surface covered with a pre-existing film of $0.31 μ$m thickness, respectively. The secondary droplet distributions are averaged over the impingement of 16 splashing drops within the train. Well before the sixteenth drop impact, the pseudo-steady state impingement dynamics are present and the influence of the transient film formation period is negligible. For all drop train impingements onto a dry surface, the first drop is deposited, therefore, the results in Fig. 12(a) represent 17 total drop impingements. For impingements onto a pre-existing film, splashing occurs at the impingement of the first drop. For both simulations, histograms show a probability density function (PDF) of secondary droplet size and the percent of secondary droplet volume. The bin size for the histograms is $Δx$⁠, or $D/80$⁠.

Comparing the secondary droplet distributions between the initially dry and pre-existing film simulations, it is clear the effect of the pre-existing film is insignificant. For both cases, the PDFs of secondary droplet size are very similar, showing that the majority of secondary droplets fall within the range $d/D=0.0125−0.1$⁠. In both cases, the largest number of secondary droplets are of size d/D = 0.025, although there are less secondary droplets of size $d/D=0.0125$ for the initially dry case. While the PDF of secondary droplet size gives quantitative evidence of the most likely secondary droplet sizes, the percent of secondary droplet volume provides further insight into the most dominant secondary droplets. As seen in Fig. 12, although the majority of secondary droplets are quite small, as evident in the PDF, they represent a negligible portion of the splashed liquid mass. By quantifying the percent of secondary droplet volume the dominance of larger secondary droplets becomes clear. In both cases, the majority of secondary droplet volume is composed of droplets of size $d/D=0.025−0.2$⁠. There is also a peak in the percent of secondary droplet volume centered at d/D = 0.4 due to the lamella rim, which registers as a large secondary droplet, before it breaks into smaller droplets.

The analysis reveals that neither the total splashed mass ratio (Fig. 6) nor the time-averaged distribution of secondary droplets show significant changes when a pre-existing film is introduced. That is mainly because the splashing and formation of secondary droplets in both surface conditions (initially dry and pre-existing film) have the same characteristics at pseudo-steady state, although they differ in the transient regime. Finally, by only focusing on the PDF of secondary droplet size a large portion of the secondary droplet volume may be overlooked (as high as 62%), leading to significant inaccuracies in predicting fuel atomization and combustion.

Equally as important for SWI sub-model development is the secondary droplet velocity. Inaccurate predictions of secondary droplet velocity may lead to discrepancies in spray rebound after impingement, which affects the fuel-air mixing and combustion efficiency. To characterize the secondary droplet velocity, the velocity in each cell that comprised the secondary droplet is averaged, excluding all interfacial cells, to ensure the surrounding gas velocity does not adversely affect the velocity calculation. In Fig. 13(a), the secondary droplet velocity is presented for the case at u = 24, with a pre-existing film thickness of 0.31 μm. Here the time-averaged velocity, averaged over the impingement of 16 drops, is normalized by the impingement velocity. The time-averaged velocity is most representative of the secondary droplets throughout the drop train impingement process. The average secondary droplet velocity increases with increasing secondary droplet size until reaching an asymptotic value of 0.8V₀ for droplets larger than $d/D=0.38$⁠. This behavior is attributed to the velocity difference between the impinging drop and the liquid film. At u = 24, the strong interaction with the liquid film causes a reduction in velocity of the lower portion of the lamella, while the lamella rim (larger drops) have already detached. Detailed velocity distributions such as those presented in Fig. 13(a) may be used to develop more accurate SWI sub-models for use in Lagrangian–Eulerian solvers.

Figure 13(b) illustrates the average secondary droplet velocity angle. The inlay in Fig. 13(b) depicts the velocity angle, θ, measured with respect to the surface. As seen in Fig. 13(b), the velocity angle exhibits a transient period in which the angle is initially larger than the pseudo-steady state value. Similar to the splashed mass ratio analysis, presented earlier, this transient period is attributed to film velocity and thickness, which undergoes a development phase. Once the film has reached pseudo-steady state the velocity angle becomes stable with slight periodic variation due to the various stages of secondary droplet formation for each impinging drop. At pseudo-steady state the average secondary droplet velocity angle is 32.5°.

Since the secondary droplets are ejected at a similar angle to the crown, the velocity angle can be compared to the crown angle model of Fredorchenko and Wang.⁶² The crown angle model⁶² is applicable to drops impinging on a thin liquid film; the model is based solely on the liquid film thickness in the form $cos(θ)=1−4H*$⁠, where $H*$ is the dimensionless film thickness defined as $H*=h0/D$⁠, where h₀ is the film thickness during impingement. In the diesel drop train simulations the pseudo-steady state film thickness is 0.31 μm and, based on the above model,⁶² the crown angle or secondary droplet angle predictions is 37.6°. Excellent agreement is seen with the model prediction for the impingement of the first two drops. Once pseudo-steady state is reached the model predicts a slightly larger crown angle than was observed in the simulations. The reason for the difference can be the interaction with the non-uniform film in both thickness and velocity, leading to crown propagation which cannot be captured by the above model.⁶² 

In this work, an in-house multiphase flow solver was used to perform high resolution simulations of high-speed micrometer-sized diesel drop train impingement. A dynamic contact angle model was developed for diesel on stainless steel using single drop impingement experiments. Following the formulations of Yarin and Weiss,²⁰ multiple drop train impingement cases were then simulated. Utilizing a passive scalar routine, each drop within the drop train was distinctly tracked providing detailed insights into the splashing of individual drops. The secondary droplet characterization routine of Markt et al.²⁹ was used to acquire time-averaged secondary droplet distributions in terms of secondary droplet size, volume, velocity, and average velocity angle. With the ability to resolve small length and time scales, the temporal evolution of secondary droplets was presented and distinct phases of secondary droplet formation were identified.

The introduction of a pre-existing liquid film provided insights into the effects of liquid films on splashing. The analysis shows that approximately 58% of the splashed liquid volume from each drop impingement originates from the impinging drop itself, while the remainder of the splashed liquid is from the already deposited film. Varying the initial film thickness has an insignificant effect on the asymptotic splashed mass ratio in a drop train impingement, contrary to what is observed in a single drop impingement. A comparison of the secondary droplet distributions between cases with an initially dry surface and a pre-existing liquid film shows that there is no significant effect on the formation of secondary droplets. Distinct phases of secondary droplet formation and lamella breakup were identified. A comparison with two instability models suggests the breakup of the lamella is due to a Rayleigh–Taylor like instability. Although the simulation results presented in the current work can facilitate the development of engine-relevant SWI sub-models by providing insights into the fuel drop impingement process, there is much work to be done in understanding the dynamics of engine-relevant drop impingement. The authors stress the need for further investigation of the effects of increased air pressure and temperature relevant to ICEs, which undoubtedly have an impact on splashing dynamics.

This material is based upon work supported by the Department of Energy, Office of Energy Efficiency and Renewable Energy (EERE), and the Department of Defense, Tank and Automotive Research, Development, and Engineering Center (TARDEC), under Award No. DE-EE007292 and the Massachusetts Clean Energy Center. The computations were performed on the HPC cluster of the UMass-Dartmouth's Center for Scientific Computing and Visualization Research supported by the ONR DURIP Grant No. N00014-18-1-2255 “A Heterogeneous Terascale Computing Cluster for the Development and Efficient Implementation of High-Order Numerical Methods,” MGHPCC and on resources provided by NSF XSEDE through the Grant Nos. ENG170004 and ENG170008.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Here we present the experimental dynamic contact angle results for single diesel drop impingement, and the selection of the dynamic contact angle model parameters. Using the image processing algorithm, the experimental results were analyzed to determine the contact line velocity (U_CL) and dynamic contact angle. In Fig. 14, the experimentally measured dynamic contact angle is shown as a function of contact line velocity. The results are the average contact angle from three repeats of the experiment to ensure accuracy of the results; an example is shown in Fig. 1 (right column). Due to the limitation of pixel resolution, there is a minimum increment of contact line velocity (one pixel per frame rate). This limitation leads to discrete values for contact line velocity; therefore, multiple contact angles may be obtained for the same contact line velocity. To show the variation of contact angle for each distinct contact line velocity, error bars are included in Fig. 14, which represents the standard deviation of the average contact angle for each distinct contact line velocity.

The maximum advancing contact angle (θ_mda) was found to be 67° as seen in Fig. 14. Although this is the maximum advancing contact angle measured from the experiment, a reduction of dynamic contact angle was observed for contact line velocities above 0.8 m/s. High contact line velocity was observed during the early stages of impingement, at which point the film front is very thin as seen in Fig. 1(b). With a very small number of pixels resolving the thin liquid film front, it becomes increasingly difficult to accurately measure the instantaneous contact angle during this time. Therefore, the expected asymptotic behavior of the advancing contact angle was not observed and the largest measured contact angle was selected for θ_mda. For the equilibrium contact angle (θ_e) a value of 10° was selected based on the last image captured from the experiments. Since no receding was observed in the experiments, the droplet continues to slowly spread until equilibrium is reached. Using the average contact angle at near zero contact line velocity would therefore overestimate the equilibrium contact angle. In Fig. 14, the apparent equilibrium contact angle is approximately 20°. Although the measurements shown are at near zero contact line velocity, they are dominated by spreading, contact angle relaxation, and drop oscillations. Therefore, such contact angle measurements are not representative of the equilibrium contact angle, and the final measured contact angle must be used. As the drop reaches its maximum spreading diameter the liquid front experiences little motion due only to drop oscillation. With no receding contact angle measurements, a receding contact angle of 1° was selected to minimize the numerical receding of the drop based on preliminary single diesel drop impingement simulations.⁵⁰ Finally, the material-related parameters (k_a and k_r) were obtained through the best fit to the experimental contact angle measurements following Tanners law [Eq. (5)]. The resulting dynamic contact angle model is presented in Fig. 14 where it is compared to the experimentally measured dynamic contact angle.