Double droplet splashing on a thin liquid film with a pseudopotential lattice Boltzmann method

This paper studies the interaction of two droplets splashing on a stationary film. A source term is included in the large-density-ratio pseudopotential lattice Boltzmann method to achieve tuneable surface tension. This model offers excellent numerical accuracy and stability for droplet impacts on liquid films. The influence of the Reynolds number, Weber number, film thickness, and horizontal/vertical distance between the droplets on the crown geometry evolution is investigated. The energy loss during the impact process and the velocity discontinuity in the liquid film are the two key factors affecting the stability and evolution process of the crown. A smaller Reynolds number or thicker liquid film enhances the energy loss and decreases the velocity discontinuity, leading to more stable side and central jets. An increase in the horizontal distance between the droplets reduces the velocity discontinuity, causing the central jet height to decrease. An increase in the Weber number does not affect the energy loss or velocity discontinuity, but the lower surface tension leads to a dramatic deformation in both the central and side jets. A vertical distance between the two droplets causes an asymmetrical evolution of the crown geometry, and postpones the breakup time of the central jet .


Introduction
After droplets impact a solid wall, they gather and form a liquid film on the surface; later droplets then impact the liquid film. This phenomenon is widespread in nature and industrial applications, and the whole process includes a series of complex multiphase flows. Understanding how droplets interact with the liquid film and each other upon impact is important when studying spray coating, inkjet printing, pesticide spraying, and soil erosion by raindrops.
Earlier studies have found that the evolution of a droplet impacting a liquid film is dominated by several non-dimensional parameters, such as the film thickness, geometry, and other properties of the wall. A large number of studies have examined the underlying complex dynamics of a single droplet splashing a stationary thin film (Alghoul et al., 2011;Burzynski & Bansmer, 2018;Cossali et al., 2004;Hung et al., 2011;Josserand & Zaleski, 2003;Rioboo et al., 2003;Vander Wal et al., 2006;Wang, Dandekar, et al., 2018 ). According to the comprehensive CONTACT Xiaolong He xiaolong − hescu@yahoo.com.hk experimental studies of Huang and Zhang (2008), the impact mode of droplet splashes on a thin film can be classified as one of five types: fusion, rebound, partial rebound, crown, and splash. An important study was carried out by Yarin and Weiss (1995), in which a quasione-dimensional model revealed the discontinuity of the kinematic velocity in the liquid film. The discontinuity of the kinematic velocity corresponds to the formation of the crown jet and exerts a great influence on the evolution process of the crown. Yarin and Weiss (1995) also indicated that the non-dimensional crown radius is a function of the non-dimensional time. However, the relationship between the two parameters was not clearly defined. Gao and Li (2015) built an empirical model through theoretical analysis and experimental studies, and predicted the variation of the crown root radius and jet stretch rate with respect to time. Compared with experimental studies, numerical simulations can provide deep insights into the flow field during the evolution process. Macroscopic multiphase models have typically used the Navier-Stokes equations and an interface tracking model (Guo et al., 2014;Mosavi et al., 2019;Motzkus et al., 2011;Shamshirband et al., 2020;Shetabivash et al., 2014;Wang et al., 2020). Several studies have analyzed the evolution process of the flow field and found that the instability of the crown rim causes the jet to break up (Allen, 1975;Nikolopoulos et al., 2007;Roisman et al., 2006). However, some discrepancies still exist, with Roisman et al. (2006) and Nikolopoulos et al. (2007) stating that the breakup is caused by the Plateau-Rayleigh instability, while Allen (1975) suggests that the breakup is caused by the Rayleigh-Taylor instability.
The lattice Boltzmann method (LBM) has been applied to simulations of multiphase phenomena in which there is a dramatic change in the interface. The most widely used LBM multiphase models fall into four categories, namely the color model (Grunau et al., 1993), pseudopotential model (Shan & Chen, 1993, 1994, freeenergy model (Swift et al., 1995;Swift et al., 1996), and phase-field model (He et al., 1999). Compared with the traditional macroscopic numerical method, pseudopotential LBM approaches have several attractive characteristics He, Zhang, Yang, et al., 2020;Li, Luo, et al., 2016): (i) the convection terms are composed of linear equations, which avoids the complex treatment of the nonlinear convection term; (ii) the pressure is solved from the ideal or non-ideal equation of state (EOS), removing the need to solve the Poisson equation; (iii) the gas-liquid interface is automatically formed by the interaction force between particles, without using the interface tracking model; (iv) the boundary treatment is relatively simple, especially for walls with different contact angles; (v) the pseudopotential model gives better numerical results than traditional methods (Ezzatneshan, 2017); and (vi) the governing equation can be divided into stream and collision components, which is conducive to parallel computing. Lee and Lin (2005) first proposed the LBM phasefield model to simulate a droplet impact on a film with a high density ratio. This model must determine the directional gradient, which requires considerable computing resources. Based on this model, the multiple-relaxationtime (MRT) collision model was introduced by Mukherjee and Abraham (2007) and the effects of the density and viscosity ratios were studied. Recently, improved pseudopotential models have been used by Fallah Kharmiani et al. (2016) and Yuan et al. (2020) to study the impact process, allowing the breakup of the crown rim to be clarified.
In most natural phenomena, large numbers of droplets splashing on a liquid film mean that the droplets not only interact with the liquid film, but also with each other. However, most studies based on LBM have focused on a single droplet impacting the liquid film. To the best of our knowledge, few studies have used LBM to examine the underlying evolution mechanics of multiple droplets impacting a liquid film (Li, Jia, et al., 2016;Raman et al., 2015;. Raman et al. (2015) first investigated the evolution of the crown following the impact of two droplets on a liquid film using a phase-field model, and found that a critical non-dimensional film thickness produces the highest side and central jets. The authors also indicated that the growth rate of the crown height and radius decreases with increasing gas density. In their study, the crown did not break up, although the governing Reynolds (Re) and Weber (We) numbers were set as 500 and 8000, respectively. Additionally, the side jets rolled inward, which was inconsistent with experimental results (Wang, Dandekar, et al., 2018). Furthermore, Li, Jia, et al. (2016) simplified the model of Lee and Lin (2005) and studied the same phenomenon with different initial impact parameters. However, Re and We in their study were limited to a narrow range, and no satellite droplets were formed.
Note that the aforementioned studies (Li, Jia, et al., 2016;Raman et al., 2015; related to two droplets splashing were based on the LBM phase-field model -the present paper describes the first systematic study of double-droplet impacts on liquid films using the LBM pseudopotential model. In this model, a source term for tuneable surface tension is applied, thus enlarging the range of We. This strategy has been successfully applied in previous research (Fallah Kharmiani et al., 2016;Yuan et al., 2020). We analyze the effects of We, Re, and the non-dimensional film thickness on the crown structure. The kinetic energy of the liquid phase and the velocity discontinuity in the film are introduced to investigate the instability of the side and central jets. The effects of the non-dimensional horizontal and vertical distance between two droplets on the evolution process of the side and center jet height and crown radius are also studied. The remainder of this paper is organized as follows. The introduction of a source term to achieve tuneable surface tension in the pseudopotential LBM is described in section 2. Mesh convergence is analyzed and the model is validated by a droplet splashing on a film in section 3, with the results quantitatively compared with those from a semi-empirical model (Gao & Li, 2015). We then study the influence of different parameters on the evolution processes in section 4. Finally, the conclusions to this study are presented in section 5.

Model description
The original pseudopotential model (Shan & Chen, 1993, 1994 has several evident drawbacks, including: (1) poor numerical stability for high-density-ratio multiphase phenomena; (2) high spurious currents at the gas-liquid interface; (3) a non-tuneable surface tension; and (4) thermodynamic inconsistency. Significant efforts have been made to overcome these drawbacks. In this section, a large-density-ratio pseudopotential model with tuneable surface tension and thermodynamic consistency is proposed. Yu and Fan (2010) indicated that the LBM with an MRT model exhibits higher numerical stability and smaller spurious currents than with the Bhatnagar-Gross-Krook collision operator. The corresponding particle distribution equation with source terms is expressed as (Yu & Fan, 2010):

Pseudopotential lattice Boltzmann method
where x is the node position, e α is the discrete velocity in the αth direction, f α and f eq α are the particle distribution functions and equilibrium distribution functions, t is the time step, I is the unit matrix,Λ is the relaxation matrix, S α is the forcing term, and C α is the source term that is used to achieve different values of the surface tension. For incompressible flows, the equilibrium distribution f eq α is given as : where ω α denotes the weights: for a two-dimensional model with nine discrete velocities (D2Q9),ω 0 = 4/9, ω 1∼4 = 1/9, and ω 5∼8 = 1/36. c s = c/ √ 3 is the lattice sound speed, where c = x/ t is the lattice constant. ρ is the macroscopic density and v is the equilibrium velocity, set to be equal to the macroscopic velocity in the present study. In the D2Q9 lattice model, e α is given by (Yu & Fan, 2010): The diagonal relaxation matrixΛ with different relaxation coefficients is defined as : where ρ is the density, e is the energy, ζ is the energy square, j is the momentum, q is the energy flux, and υ is the kinematic viscosity. Note that τ ν is related to the kinematic viscosity by υ = 1/[c 2 s (τ υ − 0.5) t]. The relaxation parameters in Equation (4) are designated as τ −1 ρ = τ −1 j = 1.0, τ −1 e = τ −1 ζ = 1.1, and τ −1 ρ = 1.1 in the present study. The transformation matrix M can be written as (Yu & Fan, 2010): where M −1 is the inverse of M. The equilibrium moment m eq can be obtained by transforming f eq into the momentum space with m eq = Mf eq (Yu & Fan, 2010):

Source terms
The external force term S proposed by  is introduced here, as this provides better numerical stability and enables larger liquid/gas density ratios to be achieved: where ε is a tuneable parameter for thermodynamic consistency. F m is the force between particles, which is given by : where G = −1 is the interaction strength and ψ is the pseudopotential. w α denotes the weights for the interaction force, and w 0 = 0, w 1∼4 = 1/3 and w 5∼8 = 1/12 (Shan, 2008). The non-ideal EOS is applied here, where ψ can be calculated as (Yuan & Schaefer, 2006): Yuan and Schaefer (2006) found that the application of the Carnahan-Starling (C-S) EOS can achieve higher numerical stability in the pseudopotential model. Thus, p EOS is calculated as: where T c is the critical temperature, and p c is the critical pressure. The parameters in Equation (10) are a = 0.25, b = 4, and R g = 1, as used by Yuan and Schaefer (2006).
Considering that the surface tension is related to the gas/liquid density ratio in the original LBM, and cannot be adjusted independently, the parameters are limited to a very small range. Therefore, the source term proposed by  is introduced to achieve tuneable surface tension: and Q is solved from : where the parameter κ is used to obtain tuneable surface tension. In the present study, the corresponding surface tension is obtained from Laplace's law; however, the process of determining the surface tension is tedious , and is not presented here.
The macroscopic density ρ and velocity u are calculated through (Yuan & Schaefer, 2006): where F = F m + G + . . . is the resultant force and G is the gravitational force, which is calculated as where g is the gravitational acceleration and ρ g is the initial density of the gas.

Mesh analysis and model validation
The pseudopotential model proposed in this section was validated in our earlier study through Laplace's law Yuan et al., 2020), and the corresponding surface tension values were obtained. This section presents simulation results for a single droplet splashing on a film to validate our model. The results are quantitatively compared with the analytical theory in an earlier study (Gao & Li, 2015). Figure 1 shows a schematic diagram of a single droplet impacting on a liquid film. The left, right, and top boundaries are set as non-equilibrium boundaries, and the bottom is set as a no-slip boundary, which is expressed as (Li et al., 2014): (15) r 0 is the droplet radius, U 0 is the initial impact velocity, and H is the liquid film height. r c is the crown root radius, defined as the horizontal distance between the inner crown root and the droplet center, as shown in Figure 1. h is the crown height, defined as the distance between the highest point of the crown and the gas-liquid interface of the liquid film.
Combined with the liquid density ρ l , viscosity of the liquid phase υ l , and surface tension σ , all of the aforementioned parameters determine several important nondimensional parameters related to the crown evolution process, including where h * is the non-dimensional film thickness and t * is the non-dimensional time. The non-dimensional radius of the crown r * and the non-dimensional crown height h * are also introduced, and are normalized as In our study, the density field of the droplet is initialized as: where (x 0 , y 0 ) is the droplet center. The density field of the liquid film is initialized as: where y 1 is the interface location. Considering the numerical stability and physical reality , the interface width w in Equations (18) and (19) is chosen as 5 lattice units (lu). The temperature is chosen as T = 0.5T c and the liquid/gas density ratio is 720, where T c = 0.0249 in the present study. The viscosity relaxation time τ ν is chosen as 0.575, with the corresponding Re and We set to 500 and 87.8, respectively. Yarin and Weiss (1995) indicated that the main parameters affecting the crown evolution are Re and We, and gravity has little influence on the crown evolution process when the droplet is small. Thus, the gravitational acceleration in the present study is set to g = 0.
Four grid densities are considered, which are 500lu × 150lu, 800lu × 240lu, 1000lu × 300lu, and 2000lu × 600lu. The corresponding droplet radii are 25lu, 40lu, 50lu, and 100lu. To compare the simulation results quantitatively, the initial Reynolds number is held at 1000 and the Weber number is fixed to 87.8. The variation in the root radius versus non-dimensional time is then examined.
The jet shape of different mesh densities at nondimensional time t * under different grid densities is shown in Figure 2. However, the shape of the crown does not vary with the grid density. The main difference is in the size of the satellite droplets, because the LBM pseudopotential model is a diffuse method and there exists a gas-liquid interface thickness with several grids, which is related to the parameter a in the C-S EOS. According to the earlier research of , the thickness is 5lu with a = 0.25. When the initial droplet radius is reduced to 25lu, the interface thickness is 1/5 of the radius, which leads to the simulation results becoming unphysical and unable to fully capture the shape of the satellite droplets. Gao and Li (2015) developed a model for predicting the evolution process of the crown root radius and jet stretch rate: where β is used to determine the energy loss coefficient, and is solved from (Gao & Li, 2015): The energy loss factor λ is obtained from a large number of experimental results (Gao & Li, 2015): In our study, the energy loss factor is λ = 0.416, and the corresponding β is 0.824. The simulation results are compared with the theoretical analysis in Figure 3(a). The numerical results for all four grid densities are consistent with the semi-empirical theoretical analysis at t * = 0.4 − 1.4. There are small discrepancies between the numerical results and the theoretical analysis for t * = 1.4 − 2.0. This may be caused by the breakup of the crown, which is not considered in the theoretical analysis. For the grid density 500lu × 150lu, oscillations occur for t * = 1.4 − 2.0, whereas the other grid densities retain smooth curves. Furthermore, the curves for 1000lu × 300lu and 2000lu × 600lu almost overlap with each other. Taking computational efficiency, precision, and stability into consideration, a grid density of 1000lu × 300lu is suitable.
Furthermore, the stretch rate of the jet obtained by the LBM with a grid density of 1000lu × 300lu is compared with the theoretical results in Figure 3(b). The vertical jet stretch rate is described as (Gao & Li, 2015): where t i is the end time of the deformation phase, given by (Gao & Li, 2015): and r i is calculated from the non-dimensional height h * as (Gao & Li, 2015): τ 0 is the time shift, obtained as (Roisman & Tropea, 2002): The stretch rate obtained by the LBM is consistent with the theoretical analysis during the no-breakup period for t * = 0.5 − 1.5, although some discrepancies appear when t * > 1.5 after the satellites have formed. This may be caused by the semi-empirical model not considering the influence of crown breakup. Some oscillations can also be observed in the stretch rate curve, and these exist even when the grid density is increased to 2000lu ×  600lu. We can identify the grid nodes that produce the maximum stretch rate changes in each time step, from which we infer that the oscillations in the stretch rate may be caused by the shape of the crown rim changing frequently. Obviously, the semi-empirical model of Gao and Li (2015) can only accurately predict the evolution process during the no-breakup period.
To determine which parameters have the greatest influence on the evolution rate of the crown radius, many researchers have conducted systematic experimental and numerical studies. For the whole evolution process, including the no-breakup period and the breakup period, the crown root radius evolution is a function of time. This was represented by Yarin and Weiss (1995) as follows: According to Shetabivash et al. (2014), for a single droplet impacting on a thin liquid film, the parameters in Equation (27)  . The relationship between the non-dimensional time and the nondimensional crown root radius is presented in Figure  4. The simulation results show that there is a linear relationship between the non-dimensional radius and the non-dimensional time. The slope of the line of best fit is 0.893, which agrees well with the value of 0.895 given by Equation (27) for both the no-breakup and breakup periods.

Numerical results and discussion
The computational and boundary conditions are shown in Figure 5(a). The size of the computational domain is set as 1201lu × 651lu. The boundaries are as indicated in Figure 1. Two droplets of radius r 0 = 50lu are located in the computational domain. The impact velocity of the droplets is set as U 0 = 0.125lu · tu −1 . The horizontal distance between the two droplets is S x , and the vertical distance between them is S y . H is the liquid film height.
The parameters used to quantify the crown structure are presented in Figure 5(b). r is the crown rim radius, h is the crown height, and S l and S r are the left and right spread lengths, respectively. The subscripts l and r indicate the left crown and the right crown, respectively. h c is the central jet height. The density ratio is set as ρ l /ρ g = 720.
To provide a deeper insight into the crown dynamics, we introduce the total kinetic energy of the liquid phase and investigate the crown deformation and breakup from an energy perspective.
where ρ f is the density of the liquid grid nodes, distinguished by the discriminant ρ f ≥ 0.5(ρ l + ρ g ). The discontinuity of the velocity in the liquid film is also discussed.
The influence of several non-dimensional parameters on the impact dynamics is investigated systematically in the following subsections: (1) Reynolds number (Re).

Reynolds number
The influence of Re on the crown evolution is studied with three different Reynolds values. Raman et al. (2015) and Li, Jia, et al. (2016) studied the influence of Re by adjusting the impact velocity and liquid viscosity. Their results show that the crown height increases as Re rises, but the evolution of the non-dimensional extension length overlaps with different Re numbers (the non-dimensional extension length is defined as S * a = S * l + S * r , where S * l = S l /r and S * r = S r /r are the left and right non-dimensional spread lengths). The impact velocity of the droplets is kept constant, and different Re numbers are obtained by changing the liquid viscosity υ l . The Re numbers considered here are Re = 200, 500, and 1000. The other parameters are as follows: We = 87.8, non-dimensional liquid film height h * = 0.25, non-dimensional horizontal distance S * x = 2.0, and non-dimensional vertical distance S * y = 0. The crown geometry at t * = 2.5 is shown in Figure 6. The central jet is formed by the collision of two crown jets between two droplets. The central jet heights are the same for Re = 500 and 1000, and are greater than for Re = 200, which is consistent with the study of Li, Jia, et al. (2016). At t * = 2.5, both the left and right crowns break up, and satellite droplets are formed with Re = 500 and 1000. In a numerical study by Raman et al. (2015), inward bending of the left and right crowns was observed at high density ratios, a phenomenon that did not concur with earlier studies (Wang, Dandekar, et al., 2018). In the present study, the left and right crown jets bend outwards during the evolution process. Furthermore, a gas-liquid interface was observed in the central jet of the simulation results reported by Raman et al. (2015), which is not consistent with the actual physical phenomenon. The colliding crowns merge into a central jet in our study, which is qualitatively consistent with earlier simulation results (Li, Jia, et al., 2016), and no gas-liquid interface appeared in the central jet.
Due to the symmetrical geometry of the crown, only the crown rim radius and height of the left crown are analyzed. The relationship between the left crown rim radius and non-dimensional time is shown in Figure 7(a) for the different Reynolds numbers. Before the breakup of the left crown, the radius increases as Re rises, because the energy lost in the collision process decreases as Re increases, as shown in Figure 7(d). At the point at which the left crown breaks, the rim radius is greatest in the case of Re = 200. However, the rim radius with Re = 1000 is still greater than for Re = 500. The left crown height varies with time under all Re values, as shown in Figure  7(b). Before the crown breaks, the crown height increases with increasing Re. This is because the deformation of the crown is dominated by inertia, and the attenuation of the inertial force is slowed by the decrease in liquid viscosity. As the left crown breaks, the height is greatest with Re = 1000. The left crown height with Re = 500 almost overlaps that with Re = 200. Figure 7(c) depicts the evolution of the central jet height with respect to nondimensional time. The central jet breakup occurs from t * = 0.75 − 1.0 for Re = 1000. Due to some defects in the pseudopotential model, the droplets are 'eaten', and no satellite droplets are observed on the central jet in Figure 6. The total kinetic energy of the liquid phase varies with non-dimensional time during the evolution process, as shown in Figure 7(d). For the three different Reynolds numbers, E u decreases over time. During the evolution of E u , three steep drops are observed, corresponding to the adjustment of the initial flow field, the droplets impacting the liquid film, and the crown colliding to form a  central jet. During the second steep drop, the decrease in E u is smaller at higher values of Re, which indicates that more kinetic energy is consumed in the collision process at lower Re. With more residual E u , the side and central jets become more unstable, leading to a significant deformation and a greater likelihood of breaking up. Figure 8 shows the flow field of the central jet and the liquid film for Re = 200 and 1000, indicating that the velocity in the liquid film below the central jet is moving in the opposite direction. With the increase in Re, the thickness of the liquid film under the central jet decreases, and the velocity discontinuity inside the liquid film increases. This leads to an increase in the velocity of the central jet and greater deformation.
The horizontal velocity discontinuity in the liquid film is now analyzed. According to the width of the central jet, the velocity 20lu from the x-axis central line is selected, and the horizontal velocity difference at h * y = 0.2h * , 0.4h * , and0.8h * is presented in Figure 9. The horizontal velocity discontinuity at t * = 1.25 for the different Re numbers is presented in Figure 9(a). The horizontal velocity discontinuity increases by about 19.3% at h * y = 0.8h * as the Re number increases from 200 to 1000, whereas the equivalent change is approximately −2.7% at h * y = 0.2h * . However, at t * = 2.5, the differences in the horizontal velocity discontinuity decrease to 5.1% and 1.3% for h * y = 0.8h * and 0.2h * , respectively. The reason for the enhanced central jet as Re increases can therefore be largely attributed to the discontinuity in the horizontal velocity on the surface.

Weber number
In this section, the effect of We on the crown evolution is studied. Fallah Kharmiani et al. (2016) showed that, as We increases, the kinetic energy consumed by surface tension decreases and the crown becomes more unstable. The crown height increases, but its thickness decreases as We rises. The crown root radius is independent of We. In the present study, different We numbers are obtained by changing the surface tension σ . The Weber numbers are set as We = 69.0, 87.8, 139.4, and 1165.5. Figure 10 shows the crown geometry at t * = 2.5. The central jet height increases as We increases from 69.0 to 139.4, and no satellite droplets form until We reaches 1165.5. As the Weber number increases, the rim of the crown becomes thinner and more unstable on both sides, and the satellite droplets decrease in size. Although the crown rims exhibit different geometries, the crown roots remain consistent, and the radius of the jet root spreading in the radical direction is not affected by We. The results agree well with those from earlier studies. Figure 11(a) shows that the left crown rim radius varies with non-dimensional time. It can be observed that, for We = 1165.5, the left crown breaks four times, compared with only one breakup for the other three Weber numbers. As We increases, the first breakup of the crown occurs earlier, and in the early stages, the crown rim radius evolves in almost the same way under the different We numbers. When t * = 2.5, the radius of the crown rim is the largest with We = 139.4, while the radii of the other three conditions are very similar. The variation in the left crown height with non-dimensional time is shown in Figure 11(b). The interface stability decreases with increasing surface tension, and more kinetic energy is consumed. Although the crown jet breaks many times at We = 1165.5, the left crown height is still higher than for the other three We numbers. Both the left crown height and the central jet increase with increasing We. The central jet does not break until We = 1165.5, for which the central jet breaks at t * = 2.5. Before breakup occurs, the central jet height increases as We rises. Figure  11(d) shows the evolution of the kinetic energy of the liquid phase under different We numbers. The E u value of the impact loss is almost the same under all four We numbers, indicating that We has little effect on the energy loss during the impact process. The side jets break up many times with We = 1165.5 (due to some unphysical phenomena that occur within the pseudopotential model), and the satellite droplets are 'eaten', which leads to an increase in the kinetic energy loss during the crown evolution process. Figure 12 shows the flow field of the central jet and the liquid film below when We = 69 and 1165.5. The change in We makes little difference to the velocity discontinuity inside the liquid film. However, the surface tension decreases as We increases, leading to the central jet becoming less stable. The flow velocity at the end of the central jet is larger with higher We values, and the central jet is more prone to deforming and breaking up.
Furthermore, the horizontal velocity discontinuity with different We numbers is analyzed. The horizontal velocity discontinuity at t * = 1.25 with different We numbers is presented in Figure 13(a). For all We numbers, the horizontal velocity discontinuity is almost the same at different depths in the liquid film. However, for t * = 2.5, the horizontal velocity discontinuity becomes increasingly different at each depth of the liquid film, and the difference reaches about 40% at the surface. A stronger surface tension will consume more kinetic energy, leading to a stable central and lower central jet.

Liquid film thickness
Earlier studies indicated that the liquid film thickness exerts a great influence on the crown evolution process. Gao and Li (2015) obtained a semi-empirical formula from many experimental data, and noted that the film thickness has a greater effect on the crown radius and vertical stretch rate evolution process than We and Re. A previous numerical study (Yuan et al., 2020) indicated that an increase in film thickness will lead to an increase in energy consumption during the collision process and an increase in the splash angle, because more kinetic energy is used to squeeze the liquid film for the crown. Snapshots of two droplets splashing on different non-dimensional film thicknesses are presented in Figure 14. As the liquid film thickness decreases, the left liquid crown becomes more prone to break up with a smaller splash angle. During the impact process, the energy loss decreases as the film becomes thinner. The vertical interaction area between the droplet and the liquid film decreases as the non-dimensional liquid film thickness decreases, which weakens the restriction of the liquid film to the crown at both sides and reduces the splash angle. Note that there are different liquid densities in the film because the LBM is a diffuse model, and so the entrapped gas in the liquid film cannot maintain a bubble form and gradually dissolves in the liquid film, before being discharged to the gas phase. Figure 15 shows the variation of the left crown radius, left crown height, and central jet height with time. Before the liquid crown breaks, the splash angle increases with increasing liquid film thickness, leading to a decrease in the left crown radius. The left liquid crown also breaks earlier with a thinner liquid film, as shown in Figure  15(a). Before the breakup of the left crown, its height reaches a maximum at h * = 0.25, which is consistent with previous results. Figure 15(c) presents the relationship between the central jet height and non-dimensional    time, and the evolution processes of h * = 0.1 and h * = 0.25 are consistent with each other. The maximum central height occurs when the non-dimensional liquid film thickness h * = 0.5. As the non-dimensional liquid film thickness increases from 0.5-0.75, the central jet height decreases. Figure 15(d) shows the evolution of the total kinetic energy of the liquid phase with different film thicknesses h * . The E u loss increases as the liquid film thickness increases, which is consistent with experimental results (Gao & Li, 2015). When h * = 0.5, although more kinetic energy is consumed in the impact process than for h * = 0.1 and 0.25, the splash angle is also greater. The combined effects of energy loss and splash angle cause the central jet to reach its maximum height with h * = 0.5. Comparing the velocity field in the central jet and the lower liquid film with h * = 0.25 and 0.75, the greater liquid film thickness produces a smaller difference in velocity discontinuity, which makes little difference to the central jet height, as shown in Figure 16.

Horizontal distance between two droplets
The relationship between the evolution of the liquid crown and the horizontal distance between the two droplets is now investigated. Note that the left liquid crown radius and height do not vary with the horizontal distance between the two droplets, so we only study the evolution of the central jet. Some discrepancies exist in earlier studies. Both Li, Jia, et al. (2016) and Raman et al. (2015) studied the central jet evolution process for non-dimensional horizontal distances S * x from 1.5-2.1. Raman et al. (2015) found that the central jet height increases with increasing horizontal distance S * x , whereas Li, Jia, et al. (2016) reported that the central jet height reaches a maximum when S * x = 1.8. In the present study, horizontal distances of S * x = 1.5, 1.8, 2.0, and 2.5 are considered. The other parameters are fixed as follows: Re = 500, We = 87.8, h * = 0.25, and S * y = 0. Snapshots of two droplets splashing with S * x = 1.5 and 2.0 are presented in Figure 17. An earlier study noted that the crown root thickness increases with time during the early stage of the evolution process when a single droplet impacts the liquid film. An increase in the distance between the two droplets leads to an increase in the evolution time before impact, which means that the central jet root thickness also increases. Yarin and Weiss (1995) noted that the jet stability is related to the kinematic discontinuity in the liquid film. The kinematic discontinuity decreases with increasing horizontal distance. The simulation results show that when S * x = 1.5, the central jet breaks. However, for S * x = 2.0, the central jet remains stable.  The time evolution of the non-dimensional height of the central jet with different horizontal distances is presented in Figure 18(a). Before the central liquid jet breaks, the central jet height decreases as the horizontal distance increases, as does the central jet vertical stretch rate. However, the only case in which two breaks are observed is when S * x = 1.5. After that, the central jet reaches its maximum height at S * x = 1.8. The kinetic energy evolution process is depicted in Figure 18(b). The curves showing the evolution of the energy almost overlap for t * = 0 − 0.4. However, the time of jet impact is delayed as the horizontal distance rises, and the kinetic energy loss becomes faster as S * x decreases. Comparing the flow field inside the central jet and the liquid film, the velocity discontinuity is greater when S * x = 1.5 than when S * x = 3.0, which leads to a higher central jet height that is more prone to breaking up ( Figure 19). Furthermore, the horizontal velocity discontinuity with different horizontal distances is presented in Figure  20. At both t * = 1.5 and 2.5, a decrease in S * x enhances the horizontal velocity discontinuity of the liquid film surface, although the horizontal velocity discontinuity remains almost the same in the middle of the liquid film and at the bottom. The horizontal velocity discontinuity with S * x = 1.5 is 1.59 times that with S * x = 2.5 at t * = 1.5, which destabilizes the central jet and makes satellite droplets more prone to form. The main reason can be ascribed to the larger horizontal distance, which leads to a longer propagation time before the crown collapse, consumes more kinetic energy, and leads to a smaller vertical stretch rate, as shown in Figure 18(a). The horizontal velocity discontinuity with S * x = 1.5 is 0.75 times that with S * x = 2.5 when t * = 2.5, which means that not only the height of the central jet, but also its width and stability, are affected.

Vertical distance between two droplets
In the previous subsections, we considered the case in which two droplets hit the liquid surface at the same time. Roisman and Tropea (2002) noted that the time interval between the droplets impacting the liquid film plays an important role in the crown evolution process. Therefore, Figure 18. Influence of horizontal distance between two droplets on time evolution of (a) central jet height and (b) total kinetic energy of the liquid phase with Re = 500, We = 87.8, h * = 0.25, and S * y = 0.  the influence of the vertical distance between droplets on the central jet evolution process is now studied. The vertical distances are set as S * y = 0.13, 0.25, and 0.38, corresponding to non-dimensional time intervals of 0.13, 0.25, and 0.38 between the two droplets impacting the surface.
Due to the time interval between the two droplets impacting the liquid surface, the evolution of the right crown is slower than that of the left crown. The vertical impact point of the two crowns decreases with increasing vertical distance between the two droplets. After the collision occurs, the higher crown bends outward and a satellite droplet is formed. The droplet size increases with increasing vertical distance, as shown in Figure 21. After impact, the central jet exhibits asymmetric evolution.
Special attention should be paid to the asymmetry of the central jet. We, therefore, study the central jet height and the relative horizontal position of the highest point in the evolution process. The relative horizontal position x * c of the highest point is introduced to determine the horizontal deviation of the central jet. This is defined as: where x max is the horizontal coordinate of the highest point and lx is the width of the computational domain. As shown in Figure 22(a), the formation time of the central jet increases as the vertical distance between the droplets increases; however, the crown heights in the later stage almost overlap. Figure 22(b) shows that the relative horizontal coordinate of the highest position varies with nondimensional time, the central jet is inclined to the left with all three vertical distances, and the offset distance decreases with increasing vertical distance. Figure 22(c) shows the evolution process of the total kinetic energy of the liquid phase with the different vertical distances. As we are using the same Re, We, and film thickness, the energy loss during the impact is very similar, resulting in the same central jet height. However, due to the different impact times and positions, x * c is different. The flow field in the center jet and the lower liquid film with S * y = 0.13 and 0.38 at t * = 2.5 are shown in Figure 23. The difference in the velocity field in the liquid film is small. However, the flow field is quite different within the central jet. Due to the later collision time in the case of S * y = 0.38, the velocity of the central jet cannot adjust to become normal to the wall, which causes the difference in x * c .

Conclusions
A high-density-ratio pseudopotential LBM with tuneable surface tension has been developed to investigate Figure 21. Influence of the non-dimensional vertical distance between two droplets: (a) S * y = 0.13, (b) S * y = 0.25, (c) S * y = 0.38 on the crown geometry of the density contour line ρ = 0.2(ρ g + ρ l ) at different non-dimensional times with Re = 500, We = 87.8, h * = 0.25, and S * x = 2.0.
the dynamics of two droplets impacting a liquid film. The model was validated by the benchmark problem of a single droplet splashing on the film. Quantitative verification of the model was performed. The crown and central jet evolution processes were investigated with different Reynolds numbers, Weber numbers, liquid film thicknesses, horizontal distances, and vertical distances. The following conclusions can be drawn:  (1) The breakup of the liquid crown occurs earlier and the central jet height increases as the Reynolds number increases. Once the central jet height reaches a threshold, it no longer varies with the Reynolds number.
(2) As the Weber number increases, the crown breaks up more easily and the central jet height increases.
(3) The increase in liquid film thickness not only consumes more energy in the impact process, but also affects the splash angle of the jet and changes the normal velocity of the central jet. Initially, the central jet height increases with increasing liquid film thickness. Once the thickness reaches a threshold, the central jet begins to decrease with further increases in liquid film thickness.
(4) The kinematic discontinuity decreases with delayed collision time, and so the central jet height decreases with increasing horizontal distance between the droplets.
(5) A vertical distance between the two droplets leads to an asymmetric evolution of the crown geometry. The offset distance between the highest point and the satellite droplet size of the central jet increases with increasing vertical distance between droplets.
Note that the basic features of two droplets impacting on the liquid film of a central axis have been captured in the present study. However, Tanaka et al. (2011) reported that the impact point and breakup mode of the central jet should be examined in three dimensions. Further investigations using a three-dimensional pseudopotential LBM are currently underway to study the asymmetric evolution of the crown geometry. Additionally, the splashing of multiple droplets on stationary and moving liquid films will be the focus of future work.