#### Abstract

The pseudopotential lattice Boltzmann method (LBM) with a tunable surface tension term is applied to study a droplet impact on a moving thin film. The effects of dimensionless parameters on the upstream and downstream crown evolution are studied, including Reynolds number (), Weber number (), liquid film thickness, and horizontal velocity of the liquid film. The movement of the liquid film causes the asymmetry development of the upstream and downstream crown. Both the instability of upstream and downstream crowns increases with the increase of and , and the upstream crown becomes more prone to break up. And a critical value of film thickness exists with the height of the upstream and downstream liquid crowns reaches the maximum value. And the velocity of liquid film restrains the development of the height of the upstream and downstream crowns, but it promotes the growth of the crown radius.

#### 1. Introduction

Droplet impact on the liquid film is a widespread phenomenon in the natural and industrial process, such as fuel injection, pesticide spraying, and ink-jet printing. Previous studies had shown that this complex multiphase phenomenon of a droplet impact on the liquid film is affected by many factors [1–4], including gas-liquid density ratio, gas-liquid viscosity ratio, and liquid film thickness. Because of the complex evolution process of the gas-liquid interface, the dynamics of this complex multiphase phenomenon is still not clear.

To date, most research studies have focused on the impact of droplets on the stationary liquid film [1, 2, 5, 6]. A lot of experiments and numerical efforts have been carried out to explain the formation and evolution process of the liquid crown at different scales. This complex phenomenon was investigated by a lot of researchers through the high-speed camera, and most of them are focused on how the dimensionless parameters affect impact behaviors and regime evolution process, including the Weber number, Reynolds number, and film thickness. Later, Huang and Zhang [7] divided the evolution process into five categories through a series of experiment results, including fusion, rebound, partial rebound, crown, and splash. However, in some cases, the liquid film moves together with the solid boundary, and the moving speed of the liquid film is even close to the same order of magnitude as the droplet impact speed. Burzynski and Bansmer [3] indicated that the droplet impacts on moving ﬁlm are different from that of on a stationary ﬁlm. Later, the researchers made different experimental studies on the droplet impact on a moving liquid film. Castrejón-Pita et al. [8] studied the process of a droplet impact on the moving liquid ﬁlm, and four different regimes are obtained, including the tranquil coalescence, violent splashing, partial and complete bouncing, and surﬁng. The results indicated the splash regimes are determined by the ratio of the velocity of the droplet and liquid film. The results of Mitchell et al. [9] showed the asymmetry of the liquid crown in the upstream and downstream when there is droplet impact on the moving liquid film. The moving liquid film can enhance the splash of the upstream crown but restrain the break of the downstream crown. Later, Gao and Li studied a single droplet impact on a flowing liquid film, and a theoretical model is developed to predict the evolution of the crown radius with the energy loss factor and derive an equation to evaluate the stretching rate of the crown jet. They also proposed a function composed by Weber and Reynolds numbers to predict splash and nonsplash phenomenon when there is droplet impact on the thin liquid film [10]. However, these research studies cannot provide a deep physical insight of the flow field of the liquid part and the surrounding gas part, and the influence of different parameters on the liquid crown evolution process cannot be given systematically, while the evolution of flow fields during the impact process can obtain the numerical simulation.

Macroscopic numerical methods are mainly based on the Navier–Stokes equations and interface tracking methods [11–13], including the volume of fluid method (VOF) and level-set method. The evolution of the crown and the satellite droplets can be analyzed from the perspective of the velocity and vorticity field obtained by numerical study. And the researchers indicated that the Rayleigh–Taylor instability is the main reason for the liquid crown fragmentation [14, 15]. As a mesoscopic method, the lattice Boltzmann method (LBM) has become a robust and effective tool to study the droplet impact on the liquid film [16–21]. By improving the phase-field model proposed by He et al. [16], Lee and Lin [22] simulated the droplet splashing process with a density ratio of 1000. It is worth mentioning that, in this model, the collision step is divided into prestreaming collision and poststreaming collision step, and several gradients and Laplacians must be solved in the calculation process, which needs a huge amount of computation resources. Later, this model is optimized by Mukherjee and Abraham [23]. The multiple-relaxation-time collision operator is applied to study the influence of different parameters, including the gas/liquid viscosity ratio and liquid film thickness. Wang et al. [24] and Li et al. [25] have built an interfacial lattice Boltzmann flux solver for high-density ratio multiphase phenomena, and the droplet splashing on the stationary thin film is also simulated successfully. Due to the existence of the shear stress in the liquid ﬁlm, the dynamic behaviors of the droplet impact on a moving liquid ﬁlm are more complex. Cheng and Luo [26] also extended the model of Lee to three dimensions to simulate the droplet impacting the moving liquid film, but the satellite droplets are not formed in his simulation. After that, Liu et al. [27] used the model of Lee and Lin [22] to study the droplet with horizontal initial velocity impact on liquid film, and the results indicated that the jet size is influence by the Bond number and the crown structure is mainly affected by the Weber number.

The pseudopotential model was firstly introduced by Shan and Chen [17, 18]. Because of its simplicity and computational efficiency, this model is widely used to simulate the complex multiphase phenomenon. Compared with conventional macroscope methods, the pseudopotential LBM has many advantages [28, 29]. Firstly, the pseudopotential LBM method is based on a series of linear equations instead of solving the nonlinear convective terms in Navier–Stokes equations. Furthermore, the Poisson equation is not solved in the pseudopotential model, and the pressure field is obtained by the equation of state (EOS). At last, the model can form the gas-liquid interface automatically, avoiding the use of the interface tracking method. However, the shortages of the multiphase model of lattice Boltzmann method is also oblivious and can be described as follows: the simulation of multiphase flow at high Reynolds number; and the numerical stability with high liquid/gas density ratio is still a challenge. Furthermore, the standard lattice Boltzmann method is not well suited to body-fitted coordinates and adaptive time-stepping. At last, the LB method is not a method choice for steady-state computations. Recently, the pseudopotential model of Li and Luo [30] is also applied by Kharmiani et al. [31] and Yuan et al. [32] to study the influence of different collision parameters on the droplet splashing on a thin film.

To the best of our knowledge, there are few studies of a droplet impact on moving liquid film by pseudopotential model. Thus, the pseudopotential lattice Boltzmann method (LBM) with a tunable surface tension term is applied to study a droplet impact on a moving thin film. The purpose of the present study is to explore the influence of different parameters on the evolution of the upstream jet and downstream jet, including the Reynolds number, Weber number, liquid film thickness, and liquid film velocity. And the results are also validated by earlier experimental and theoretical results.

This paper is organized as follows. The lattice Boltzmann method with a tunable surface tension term is described in Section 2. In Section 3, the model is validated by several multiphase benchmarks. The evolution of a droplet splashing on the moving liquid film under different conditions is studied, and the results are presented in Section 4. Finally, our conclusions are given in Section 5.

#### 2. Lattice Boltzmann Method with a Tunable Surface Tension Term

Compared with the pseudopotential model with MRT collision operator, the pseudopotential model with BGK collision operator has poor numerical stability and large spurious currents [33]. To achieve the numerical stability and small spurious currents with a high-density ratio multiphase phenomenon, the pseudopotential model with MRT collision operator is applied in the present study, and the corresponding particle distribution equation with extra terms is [30]where is the particle distribution function, is the equilibrium distribution, is the spatial position, is the discrete velocity of the *i*th direction, is the time step, is the unit matrix, and is the collision matrix.

The discrete velocities of D2Q9 lattice applied in present study are expressed as follows [28]:where is the lattice constant and the diagonal matrix is composed of different relaxation coefficients and can be expressed as follows [28]: is the relaxation time, and the corresponding fluid kinematic viscosity can be solved from , where is the lattice sound speed. The relaxation parameters in equation (3) are chosen as and in our study. is the orthogonal transformation matrix, and is the inverse matrix of . The equilibrium moment can be obtained through [28]:

In the present study, the modified external forcing scheme proposed by Li and Luo [30] is applied to achieve better numerical stability, where can be expressed asand is a parameter to achieve thermodynamic consistency. is the interaction force between fluid particles, which is shown as follows [33]:where is the interaction strength between two particles and is set as in the present study, and is the interparticle potential. are the weights in different directions, for the D2Q9 LBM model, and [34]. While a non-ideal-gas EOS is introduced in the pseudopotential model, can be calculated as follows [35]:where is the pressure calculated by the Carnahan–Starling (C-S) EOS in our study and can be expressed as follows [31]:where and , is the critical temperature, and is the critical pressure. The parameters of EOS are chosen as in our study. Furthermore, a source term is introduced in this model to achieve the tunable surface tension, where the source term is expressed as follows [30]:and can be solved from [30]where is a parameter to tune surface tension and has a value between 0 and 1.

Finally, the macroscopic density and macroscopic velocity are obtained bywhere is the total force acting on the fluid particles and is the gravity force.

The current pseudopotential method has been successfully applied to simulate the droplet impact on a stationary liquid film in our earlier studies [32]. More details about this method can be found in our previous work. The proposed method is also used to study complex multiphase flows with heat and mass transfer, such as boiling [28] and cavitation [36].

#### 3. Model Validation

##### 3.1. Laplace Law and Maxwell Construction

The Laplace law is used to validate the proposed method in the present study: , where is the surface tension and is the droplet radius after the fluid system reaches an equilibrium state, and is the pressure difference between the internal pressure of the bubble and the external liquid. The computational domain is set as a 401*lu* × 401*lu* square region, periodic boundary conditions are applied at all four boundaries of the computational domain, and a stationary droplet of radius is initialized at the computational domain center. Four different initial radii of droplets are chosen to validate Laplace’s law. In the C-S EOS, the corresponding temperature is chosen as . The relaxation time is chosen as , and a series of are selected to tune the surface tension.

Figure 1 shows the comparison of the density coexistence curve obtained by LBM with different parameters and the Maxwell coexistence curve. The results indicate that when and 0.12573 and the density of the liquid phase agrees well with the theoretical value, while a mild difference was observed at the gas density when . However, when , both the liquid and gas density agrees well with the theoretical value.

Figure 2 shows there is a linear relationship between the pressure difference and the droplet radius at , with density ratio after the fluid system reaches equilibrium. And the numerical results of different agree well with the theoretical analysis.

##### 3.2. Droplet Impact on a Stationary Liquid Film

The computational domain is set as a 1001 × 301 rectangle region as shown in Figure 3, where is the radius of the spherical droplet, is the thickness of the liquid film, is the initial droplet velocity, and is the liquid film velocity. Nonequilibrium extrapolation boundary is applied at both left and right boundaries, and the no-slip boundary is applied at both the top and bottom boundaries [37].

Combined with the liquid density , liquid viscosity , and surface tension , all the above parameters can be reconstituted into several dimensionless parameters which are governing the droplet splashing on the liquid film, includingwhere is the Reynolds number, is the Weber number, is the dimensionless film thickness, is the dimensionless time, and is the dimensionless velocity.

The density field around the droplet is initialized as follows:where is the center of the droplet. And the initial density field of the liquid film is initialized as follows:where is the interface of the liquid film and the gas. According to earlier researches, the gas-liquid initial interface width in equations (13) and (14) is set as 5.

The dimensionless radius of the crown , the upstream splash angle , the downstream splash angle , the dimensionless upstream jet height , and the dimensionless upstream jet height are also introduced in the present study. All aforementioned parameters can be defined aswhere is the radius of the crown and and are the height of the upstream and downstream jet as shown in Figure 4.

A droplet impact on a stationary liquid film is introduced to verify the reliability of the proposed pseudopotential method. The radius of the spherical droplet is , and the thickness of the liquid film is . The temperature is chosen as with corresponding liquid/gas density ratio . The initial velocity is set as 0. 125, and the liquid film velocity is set as . is chosen as 0.5375 and is chosen as 0.5, so the corresponding Reynolds number and Weber number are 500 and 87.8, respectively. The snapshots of a droplet impact on the stationary liquid surface are shown in Figure 5. And the relationship between the dimensionless crown radius and dimensionless time is shown in Figure 6.

The relationship between the dimensionless crown radius and dimensionless is presented in Figure 6. Earlier researches indicated that a linear relationship exists between the dimensionless crown radius and dimensionless time , and the best curve fitted on our results is , which is known as the power law. And the best fit curve in the present study also agrees well with the experimental and theoretical results reported in earlier pieces of literature [1, 2, 11].

Later, Gao and Li [10] presented a detail semiempirical model to predict the crown radius evolution process; the model can be described aswhere is a coefficient related to the energy loss factor and is calculated as

The energy loss factor of the droplet impact on the stationary liquid film is determined by three-dimensionless numbers, including , , and :

The comparison between the LBM numerical result and the theoretical analysis is shown in Figure 7. The numerical results are located between two lines with and , and the corresponding is 0.82 and 1.13. The discrepancies between the numerical result and theoretical analysis may be caused by the assumption of the energy loss factor; the splash jet is formed in our simulation results before the deformation phase occurs; the liquid film in splash point is larger than . The same phenomenon is also observed in some experimental results [5] and numerical results [23, 24, 26, 31].

Furthermore, the stretch rate of the crown obtained by LBM is compared with the theoretical analysis of Gao and Li [10] and is shown in Figure 8. According to the theoretical analysis of Gao [10], the stretch rate of the crown jet when there is a droplet impact on the stationary liquid film is calculated as follows:where is the ending time of the deformation phase and can be calculated asand the corresponding radius is a function of and is expressed as

The stretch rate obtained by LBM agrees well with the theoretical analyze of Gao with and suddenly increases at . The reason may be that the satellite droplet is forming and leads to the increase of the velocity of the crown rim. After the satellite droplet formed, the stretch rate decreased suddenly at .

#### 4. Numerical Results and Discussion

In this part, the influences of the parameters are studied on the moving liquid film impacted by a droplet, including the Reynolds number, Weber number, liquid film thickness, and liquid film velocity. The corresponding simulation conditions are the same as described in 3.2 except the liquid film velocity which is not 0. And the results are analyzed and validated by the model of Gao and Li [10]. The parameters of different simulation cases are shown in Table 1.

The semiempirical model of Gao and Li [10] is also applied to validate our simulation results with flow liquid film in this section; the energy loss factor of the droplet impact on the flow liquid film is determined by four dimensionless numbers, including , , , and , and is expressed as follows:

Combined with equations (16) and (17), the crown radius evolution process can be predicted.

##### 4.1. Reynolds Number

The snapshots of a droplet impact on the moving liquid surface are presented in Figure 9. Both the upstream and downstream crown thickness become thinner and higher with the increasing of , and the satellite droplets are easier to form. Affected by the moving liquid film, the upstream and downstream crowns appear asymmetry. Due to the large viscosity, no crown is formed with . On the contrary, when , both the upstream and downstream crown are formed. The upstream crown breaks at with , and satellite droplets are formed, as shown in Figure 9(b). Compared with our earlier research [28], these satellite droplets do not form when the droplet impacts the stationary liquid film, which means the moving liquid film accelerates the breaking process of the upstream crown. When , satellite droplets are formed with , as shown in Figure 9(c). However, the crown only breaks at with when there is droplet impact on the stationary liquid film. And the downstream crown does not break in the whole impact process in all cases.

**(a)**

**(b)**

**(c)**

Figure 10 plots the time evolution of the upstream and downstream crown height and crown radius with different . Both the upstream and downstream height increase with the increasing of number. The upstream crown height is higher than the downstream with the same number as shown in Figures 10(a) and 10(b). The reason is that the moving liquid is squeezed by the droplet which has reverse velocity and leads the upstream crown to be more likely to break with the increase of . Because of the break of the upstream crowns, the height difference between the upstream crowns is smaller than that of the downstream crowns with and 1000. The crown radius increases with the increase of , as shown in Figure 10(c), while the crown radius does not change with when there is droplets impact on the stationary liquid film except . The reason is that, with the increase of , the liquid viscosity decreases and the energy consumed in the collision also decreases. More kinetic energy is used to overcome the liquid film velocity and excite the liquid crown. Our simulation results also agree well with the prediction of Gao and Li [10] except ; this may be caused by the significance of the Reynolds number is underestimated. And equation (22) may be not suitable for predicting the energy loss rate with low Reynolds number, and most kinetic energy is lost in the collision process with high liquid viscosity.

**(a)**

**(b)**

**(c)**

##### 4.2. Weber Number

The parameter can be used to adjust the surface tension to obtain different numbers. High number means low surface tension and leads crown to be more prone to break up. And the influence of on the crown deformation is discussed in this part. When dimensionless time , the satellite droplet has been formed at the upstream crown with . When dimensionless time , the satellite droplets are formed with all four cases, while the downstream crown remains stable as shown in Figure 11(b). And the secondary breakup happens at the upstream crown at , and downstream crown remains stable.

**(a)**

**(b)**

**(c)**

The evolution of dimensionless upstream height , downstream height , and radius with under different numbers is shown in Figure 12. Both the height of the upstream and downstream crown height increases with the increase of number. But the crown radius also does not change with the number. Because the downstream crown does not break and all its kinetic energy is used for the growth of the downstream crown, the growth rate of the downstream crown is higher than the upstream crown as shown in Figures 12(a) and 12(b). The crown radius evolution processes of different cases are shown in Figure 12(c), the evolution processes of different numbers almost overlap with each other, and the simulation results agree well with the prediction of Gao [10]. The energy loss factor ranges from 0.49 to 0.52, which means the number has little effect on the energy loss in the impact process.

**(a)**

**(b)**

**(c)**

##### 4.3. Film Thickness

As mentioned above, both and influence the crown regime. Now we turn our attention to the influence of the dimensionless liquid film thickness , and the dimensionless film thickness varies from 0.1 to 0.5. The snapshots of crown evolution with different dimensionless liquid film thickness are shown in Figure 13. When the liquid film thickness is quite small, the upstream becomes stable during all impact process while the downstream crown break up at with , which is quite different from other cases. When increases to 0.25, the upstream crown breaks up at , while the downstream crown stays stable. Compared with , both the upstream and downstream splash angles become larger at and lead the crown becomes higher. When increases to 0.5, both the spreading and deformation tendency become weak, and the upstream and downstream crowns stay stable.

**(a)**

**(b)**

**(c)**

Figure 14 shows the evolution of the upstream crown height, downstream crown height, and crown radius. Both the upstream and downstream crown heights increase rapidly with the liquid film thickness increases from 0.1 to 0.25. However, when the liquid film thickness increases from 0.25 to 0.5, both the upstream and downstream crown heights decrease. When , there is a sudden change with the downstream height; this is caused by the satellite droplet formation. The crown radius is shown in Figure 14(c), the radius at increases with the increase of , and the results agree well with the prediction of Gao [10]. The energy loss factor increases from 0.48 to 0.58 with the dimensionless liquid film height decreasing from 0.5 to 0.1, which means the energy loss in impact increases with the increase of liquid film height. The results also agree well with the studies of Mukherjee and Abraham [23] and Kharmiani et al. [31].

**(a)**

**(b)**

**(c)**

##### 4.4. Horizontal Velocity of the Liquid Film

The horizontal velocity of the liquid film is discussed in this part. Figure 15 shows the crown evolution process of a drop impact on the moving liquid film with . When there is droplet impact on the moving liquid film, the asymmetric crown is found in Figures 15(a) and 15(b), and the instability of the upstream crown is enhanced and the downstream crown is attenuated. Satellite droplets are formed in both cases. With the increase of , the suppression effect becomes more obvious. Comparing three cases, it can be observed that , the downstream crown becomes thicker and shorter, and this increases the stability of the downstream crown. The present results agree well with the numerical results of Cheng and Luo [26].

**(a)**

**(b)**

Figure 16 shows the evolution of the upstream crown height, downstream crown height, and crown radius with different liquid film velocity. Both the upstream and downstream crown heights decrease with the dimensionless liquid film velocity increasing from 0.5 to 1.0. However, the crown radius in different cases almost overlaps with each other with the increases of the dimensionless velocity as shown in Figure 16(c). It is also found that the energy loss factor increases from 0.52 to 0.62 with the dimensionless velocity decreasing from 0.5 to 1.0.

**(a)**

**(b)**

**(c)**

Figure 17 shows a comparison of the liquid velocity fields of the upstream crown and downstream crown between and 0.5. When , the spreading of the crown is dominated by the radial vectors as shown in Figures 17(a) and 17(c). However, the crown formation is dominated by both the liquid film velocity and the radial velocity. The upstream crown is squeezed by the film velocity and the radial velocity and leads to the upstream crown becoming unstable, as shown in Figure 17(b). However, the radial velocity and the film velocity have the same direction at the downstream crown, which makes the downstream crown become more stable.

**(a)**

**(b)**

**(c)**

**(d)**

#### 5. Conclusions

A droplet impact on the moving thin film is studied by the tunable surface tension LBM pseudopotential model. The influence of number, number, interface thickness, and liquid film velocity on the upstream and downstream crown is proposed. The following conclusions can be drawn:(1)The movement of liquid film causes the asymmetry development of upstream and downstream crown. The instability of upstream and downstream crowns increases with the increase of and , and the upstream crown becomes more prone to break up.(2)The height of the upstream and downstream crown varies with the thickness of the liquid film. When , the height of the upstream and downstream liquid crowns reaches the maximum value.(3)The velocity of liquid film restrains the development of the height of the upstream and downstream crowns, but it promotes the growth of the crown radius. With the increase of the velocity of the liquid film, the upstream liquid crowns become unstable and prone to break up and produce satellite droplets, but the downstream liquid crown becomes more stable.

#### Data Availability

The processed data required to reproduce these findings cannot be shared at this time as the data also form part of an ongoing study.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 51979183) and the China National Key R&D Plan (Grant no. 2016YFC0402006).