Droplet Impact on a Moving Thin Film with Pseudopotential Lattice Boltzmann Method

Chongqing Southwest Research Institute for Water Transport Engineering, Chongqing Jiaotong University, Chongqing, 400074, China Key Laboratory of Inland Waterway Regulation Engineering, Chongqing Jiaotong University, Chongqing, 400074, China Chongqing Xike Consultation Center for Water Transport Engineering, Chongqing Jiaotong University, Chongqing, 400074, China State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu, 610065, China


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][2][3][4], including gas-liquid density ratio, gasliquid 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. is complex phenomenon was investigated by a lot of researchers through the highspeed 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 film are different from that of on a stationary film. 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 film, and four different regimes are obtained, including the tranquil coalescence, violent splashing, partial and complete bouncing, and surfing. e results indicated the splash regimes are determined by the ratio of the velocity of the droplet and liquid film. e 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. e 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. ey 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][12][13], including the volume of fluid method (VOF) and level-set method. e 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][17][18][19][20][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]. e 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 film, the dynamic behaviors of the droplet impact on a moving liquid film 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. e 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. us, the pseudopotential lattice Boltzmann method (LBM) with a tunable surface tension term is applied to study a droplet impact on a moving thin film. e 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. is paper is organized as follows. e 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. e 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.

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] 2 Mathematical Problems in Engineering where f i is the particle distribution function, f eq j is the equilibrium distribution, x is the spatial position, e i is the discrete velocity of the ith direction, Δt is the time step, I is the unit matrix, and Λ � M − 1 ΛM is the collision matrix. e discrete velocities of D2Q9 lattice applied in present study are expressed as follows [28]: where c � Δx/Δt 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 1 in our study. M is the orthogonal transformation matrix, and M − 1 is the inverse matrix of M. e equilibrium moment m eq can be obtained through m eq � Mf eq [28]: In the present study, the modified external forcing scheme proposed by Li and Luo [30] is applied to achieve better numerical stability, where S can be expressed as and ε is a parameter to achieve thermodynamic consistency. F m is the interaction force between fluid particles, which is shown as follows [33]: where G is the interaction strength between two particles and is set as G � − 1 in the present study, and ψ is the interparticle potential. ω i are the weights in different directions, for the D2Q9 LBM model, ω 1− 4 � 1/3 and ω 5− 8 � 1/12 [34]. While a non-ideal-gas EOS is introduced in the pseudopotential model, ψ can be calculated as follows [35]: where p EOS is the pressure calculated by the Carnahan-Starling (C-S) EOS in our study and can be expressed as follows [31]: where a � 0.4963R 2 g T 2 c /p c and b � 0.1873R g T c /p c , T c is the critical temperature, and p c is the critical pressure. e parameters of EOS are chosen as a � 0.25, b � 4, andR g � 1 in our study. Furthermore, a source term C is introduced in this model to achieve the tunable surface tension, where the source term C is expressed as follows [30]: Mathematical Problems in Engineering and Q 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 u are obtained by where F � F m + G + · · · is the total force acting on the fluid particles and G is the gravity force. e 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. e proposed method is also used to study complex multiphase flows with heat and mass transfer, such as boiling [28] and cavitation [36].

Laplace Law and Maxwell Construction.
e Laplace law is used to validate the proposed method in the present study: Δp � (σ/R), where σ is the surface tension and R is the droplet radius after the fluid system reaches an equilibrium state, and Δp is the pressure difference between the internal pressure of the bubble and the external liquid. e computational domain is set as a 401lu × 401lu square region, periodic boundary conditions are applied at all four boundaries of the computational domain, and a stationary droplet of radius R 0 is initialized at the computational domain center. Four different initial radii of droplets R 0 � 40, 50, 60, 70lu are chosen to validate Laplace's law. In the C-S EOS, the corresponding temperature is chosen as T e � 0.5, 0.55, 0.6, 0.7, 0.8, 0.9, 0.95T c . e relaxation time is chosen as τ ] � 0.5375, and a series of κ � 0.0, 0.2, 0.5, 0.9 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. e results indicate that when ε � 0.12 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 T e < 0.7T c . However, when ε � 0.11895, both the liquid and gas density agrees well with the theoretical value. Figure 2 shows there is a linear relationship between the pressure difference Δp and the droplet radius (1/R) at T e � 0.5T c , with density ratio (ρ l /ρ g ) � 720 after the fluid system reaches equilibrium. And the numerical results of different κs agree well with the theoretical analysis.

Droplet Impact on a Stationary Liquid Film.
e computational domain is set as a 1001lu × 301lu rectangle region as shown in Figure 3, where R 0 is the radius of the spherical droplet, H is the thickness of the liquid film, U 0 is the initial droplet velocity, and U f 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 ρ l , liquid viscosity ] l , and surface tension σ, all the above parameters can be reconstituted into several dimensionless parameters which are governing the droplet splashing on the liquid film, including where Re is the Reynolds number, We is the Weber number, h * is the dimensionless film thickness, t * is the dimensionless time, and u * is the dimensionless velocity. e density field around the droplet is initialized as follows: where (x 0 , y 0 ) is the center of the droplet. And the initial density field of the liquid film is initialized as follows: where y 1 is the interface of the liquid film and the gas. According to earlier researches, the gas-liquid initial interface width w in equations (13) and (14) is set as 5lu. e dimensionless radius of the crown R * , the upstream splash angle θ u , the downstream splash angle θ d , the dimensionless upstream jet height H * u , and the dimensionless upstream jet height H * d are also introduced in the present study. All aforementioned parameters can be defined as where R c is the radius of the crown and H u and H d 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. e radius of the spherical droplet R 0 is 50lu, and the thickness of the liquid film is H � 25lu. e temperature is chosen as T e � 0.5T c with corresponding liquid/gas density ratio (ρ l /ρ g ) � 720. e initial velocity U 0 is set as 0. 125lu · tu − 1 , and the liquid film velocity is set as U f � 0. τ ] 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. e snapshots of a droplet impact on the stationary liquid surface are shown in Figure 5. And the relationship between the dimensionless crown radius R * and dimensionless time t * is shown in Figure 6. e relationship between the dimensionless crown radius R * and dimensionless t * is presented in Figure 6. Earlier researches indicated that a linear relationship exists between the dimensionless crown radius R * and dimensionless time �� t * √ , and the best curve fitted on our results is R * � 1.13 �� t * √ , 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 as where β is a coefficient related to the energy loss factor λ and is calculated as e energy loss factor λ s of the droplet impact on the stationary liquid film is determined by three-dimensionless numbers, including Re, We, and h * : e comparison between the LBM numerical result and the theoretical analysis is shown in Figure 7. e numerical results are located between two lines with λ � 0.42 and λ � 0.78, and the corresponding β is 0.82 and 1.13. e 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 h * . e 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 t i is the ending time of the deformation phase and can be calculated as and the corresponding radius r i is a function of h * and is expressed as e stretch rate obtained by LBM agrees well with the theoretical analyze of Gao with t * � 0.5 ∼ 1.2 and suddenly increases at t * � 1.2 ∼ 1.5. e 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 t * � 1.5 ∼ 2.0.

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. e 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]. e parameters of different simulation cases are shown in Table 1.
e 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 λ f of the droplet impact on the flow liquid film is determined by four dimensionless numbers, including Re, We, h * , and u * , and is expressed as follows: Combined with equations (16) and (17), the crown radius evolution process can be predicted.

Reynolds Number.
e 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 Re, 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 Re � 40. On the contrary, when Re � 200, 500, 1000, both the upstream and downstream crown are formed. e upstream crown breaks at t * � 1.25 with Re � 500, 1000, 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 t * � 2.0, satellite droplets are formed with  Re � 200, 500, 1000, as shown in Figure 9(c). However, the crown only breaks at t * � 2.0 with Re � 1000 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. Figure 10 plots the time evolution of the upstream and downstream crown height and crown radius with different Re. Both the upstream and downstream height increase with the increasing of Re number. e upstream crown height is higher than the downstream with the same Re number as shown in Figures 10(a) and 10(b). e 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 Re. Because of the break of the upstream crowns, the height difference between the upstream crowns is smaller than that of the downstream crowns with Re � 500 and 1000. e crown radius increases with the increase of Re, as shown in Figure 10(c), while the crown radius does not change with Re when there is droplets impact on the stationary liquid film except Re � 40. e reason is that, with the increase of Re, 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  Mathematical Problems in Engineering simulation results also agree well with the prediction of Gao and Li [10] except Re � 40; 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.

Weber Number.
e parameter κ can be used to adjust the surface tension σ to obtain different We numbers. High We number means low surface tension and leads crown to be more prone to break up. And the influence of We on the crown deformation is discussed in this part. When dimensionless time t * � 0.75, the satellite droplet has been formed at the upstream crown with We � 1165.5. When dimensionless time t * � 1.5, 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 t * � 2.0, and downstream crown remains stable. e evolution of dimensionless upstream height H * u , downstream height H * d , and radius R * with t * under different We numbers is shown in Figure 12. Both the height of the upstream and downstream crown height increases with the increase of We number. But the crown radius also does not change with the We 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). e crown radius evolution processes of different cases are shown in Figure 12(c), the evolution processes of different We numbers almost overlap with each other, and the simulation results agree well with the prediction of Gao [10]. e energy loss factor ranges  from 0.49 to 0.52, which means the We number has little effect on the energy loss in the impact process.

Film ickness.
As mentioned above, both Re and We influence the crown regime. Now we turn our attention to the influence of the dimensionless liquid film thickness h * , and the dimensionless film thickness varies from 0.1 to 0.5. e snapshots of crown evolution with different dimensionless liquid film thickness h * � 0.1, 0.2, and 0.5 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 t * � 2.0 with h * � 0.1, which is quite different from other cases. When h * increases to 0.25, the upstream crown breaks up at t * � 2.0, while the downstream crown stays stable. Compared with h * � 0.1, both the upstream and downstream splash angles become larger at h * � 0.25 and lead the crown becomes higher. When h * increases to 0.5, both the spreading and deformation tendency become weak, and the upstream and downstream crowns stay stable. 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 h * � 0.1, there is a sudden change with the downstream height; this is caused by the satellite droplet formation. e crown radius is shown in Figure 14(c), the radius at t * � 2.5 increases with the increase of h * , and the results agree well with the prediction of Gao [10]. e energy loss factor increases from 0.48 to 0.58 with the dimensionless liquid film height h * decreasing from 0.5 to 0.1, which means the energy loss in impact increases with the increase of liquid film height. e results also agree well with the studies of Mukherjee and Abraham [23] and Kharmiani et al. [31].

Horizontal Velocity of the Liquid Film.
e 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 u * � 0.5 and 1. 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 u * , the suppression effect becomes more obvious. Comparing three cases, it can be observed that u * � 0.5 and 1.0, the downstream crown becomes thicker and shorter, and this increases the stability of the downstream crown. e present results agree well with the numerical results of Cheng and Luo [26]. 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. Figure 17 shows a comparison of the liquid velocity fields of the upstream crown and downstream crown between u * � 0.0 and 0.5. When u * � 0.0, 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. e 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.

Conclusions
A droplet impact on the moving thin film is studied by the tunable surface tension LBM pseudopotential model. e influence of Re number, We number, interface thickness, and liquid film velocity on the upstream and downstream crown is proposed. e following conclusions can be drawn: (1) e movement of liquid film causes the asymmetry development of upstream and downstream crown. e instability of upstream and downstream crowns increases with the increase of Re and We, and the upstream crown becomes more prone to break up.
(2) e height of the upstream and downstream crown varies with the thickness of the liquid film. When h * � 0.25, the height of the upstream and downstream liquid crowns reaches the maximum value. (3) e 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 e 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
e authors declare that they have no conflicts of interest.