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

Authors: Jinchao He; Hao Yuan; Xiaolong He; Chunhang Xie; Haonan Peng; Ruichang Hu
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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. TheRe effects of dimensionless parameters on the upstream and downstream crown evolution are studied, including Reynolds number (Re), Weber number (We), 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 Re and We, 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.


---

## Body

## 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 Section2. 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](1)fαx+eαΔt,t+Δt=fαx,t−∑αΛ¯αβfβ−fβeqx,t−∑jIαβ−Λ¯αβS¯βΔtx,t+Cα,where fi is the particle distribution function, fjeq is the equilibrium distribution, x is the spatial position, ei 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.The discrete velocities of D2Q9 lattice applied in present study are expressed as follows [28]:(2)eα=0,0,α=0,ccosα−1π2,sinα−1π2,α=1−4,2ccos2α−9π4,sin2α−9π4,α=5−8,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]:(3)Λ=diagτρ−1,τe−1,τζ−1,τj−1,τq−1,τj−1,τq−1,τν−1,τν−1.τν is the relaxation time, and the corresponding fluid kinematic viscosity can be solved from ν=1/cs2τν−0.5Δt, where cs is the lattice sound speed. The relaxation parameters in equation (3) are chosen as τρ−1=τj−1=1.0,τe−1=τζ−1=1.1, and τρ−1=1.1 in our study. M is the orthogonal transformation matrix, and M−1 is the inverse matrix of M. The equilibrium moment meq can be obtained through meq=Mfeq [28]:(4)meq=ρ1,−2+3ux2+uy2,1−3ux2+uy2,ux,−ux,uy,−uy,ux−uy,uxuyT.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(5)S=ρ0,6uxFx+uyFy+12εFm2ψ2Δtτe−0.5,−6uxFx+uyFy−12εFm2ψ2Δtτe−0.5,Fx,−Fx,Fy,−Fy,2uxFx−uyFy,uxFx+uyFyT,and ε is a parameter to achieve thermodynamic consistency. Fm is the interaction force between fluid particles, which is shown as follows [33]:(6)Fm=−Gψx∑iωiψx+eiΔt,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]:(7)ψρ=2pEOS−ρcs2Gc2,where pEOS is the pressure calculated by the Carnahan–Starling (C-S) EOS in our study and can be expressed as follows [31]:(8)pEOS=ρRgTe1+bρ/4+bρ/42−bρ/431−bρ/43−aρ2,where a=0.4963Rg2Tc2/pc and b=0.1873RgTc/pc, Tc is the critical temperature, and pc is the critical pressure. The parameters of EOS are chosen as a=0.25,b=4,andRg=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]:(9)C=0,1.5τe−1Qxx+Qyy,−1.5τζ−1Qxx+Qyy,0,0,0,0,−τν−1Qxx−Qyy,−τν−1Qxy,and Q can be solved from [30](10)Q=κG2ψx∑iωiψx+eiΔt−ψxeiei,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(11)ρ=∑ifi,ρu=∑ifiei+12FΔt,where F=Fm+G+⋯ is the total force acting on the fluid particles and G 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:Δ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. The 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 R0 is initialized at the computational domain center. Four different initial radii of droplets R0=40,50,60,70lu are chosen to validate Laplace’s law. In the C-S EOS, the corresponding temperature is chosen as Te=0.5,0.55,0.6,0.7,0.8,0.9,0.95Tc. The 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.Figure1 shows the comparison of the density coexistence curve obtained by LBM with different parameters ε and the Maxwell coexistence curve. The 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 Te<0.7Tc. However, when ε=0.11895, both the liquid and gas density agrees well with the theoretical value.Figure 1
Comparison of the LBM coexistence density with differentε with the Maxwell coexistence curve.Figure2 shows there is a linear relationship between the pressure difference Δp and the droplet radius 1/R at Te=0.5Tc, 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.Figure 2
The relationship of pressure difference and drop radius withTe=0.5Tc and different κs.


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

The computational domain is set as a 1001lu × 301lu rectangle region as shown in Figure 3, where R0 is the radius of the spherical droplet, H is the thickness of the liquid film, U0 is the initial droplet velocity, and Uf 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].Figure 3
Schematic of the computational domain for the droplet splashing on the moving film.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(12)Re=2R0U0νl,We=2ρlR0U02σ,h∗=H2R0,t∗=U0t2R0,u∗=U0Uf,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.The density field around the droplet is initialized as follows:(13)ρx,y=ρl+ρg2+ρl−ρg2tanh2x−x02+y−y02−R0w,where x0,y0 is the center of the droplet. And the initial density field of the liquid film is initialized as follows:(14)ρx,y=ρl+ρg2+ρl−ρg2tanh2y1−Hw,where y1 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.The dimensionless radius of the crownR∗, the upstream splash angle θu, the downstream splash angle θd, the dimensionless upstream jet height Hu∗, and the dimensionless upstream jet height Hd∗ are also introduced in the present study. All aforementioned parameters can be defined as(15)R∗=Rc2R0,Hu∗=Hu2R0,Hd∗=Hd2R0,where Rc is the radius of the crown and Hu and Hd are the height of the upstream and downstream jet as shown in Figure 4.Figure 4
Crown parameters for the droplet splashing on the thin film.A droplet impact on a stationary liquid film is introduced to verify the reliability of the proposed pseudopotential method. The radius of the spherical dropletR0 is 50lu, and the thickness of the liquid film is H=25lu. The temperature is chosen as Te=0.5Tc with corresponding liquid/gas density ratio ρl/ρg=720. The initial velocity U0 is set as 0. 125lu⋅tu−1, and the liquid film velocity is set as Uf=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. The 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.Figure 5
Snapshots of the droplet impact on the thin film at different time.Figure 6
Relationship between the dimensionless crown radiusR∗ and dimensionless time t∗.The relationship between the dimensionless crown radiusR∗ 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.13t∗, 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(16)R∗=βt∗+16h∗−13h∗−16h∗1/2,where β is a coefficient related to the energy loss factor λ and is calculated as(17)β=2λ23h∗1/4.The energy loss factorλs of the droplet impact on the stationary liquid film is determined by three-dimensionless numbers, including Re, We, and h∗:(18)λs=0.26Re−0.05We0.07h∗0.34.,The comparison between the LBM numerical result and the theoretical analysis is shown in Figure7. The numerical results are located between two lines with λ=0.42 and λ=0.78, 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 h∗. The same phenomenon is also observed in some experimental results [5] and numerical results [23, 24, 26, 31].Figure 7
Comparison of the crown radius evolution between numerical result and the semiempirical model of Gao.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:(19)St=t∗+ti+2τ0−1,where ti is the ending time of the deformation phase and can be calculated as(20)ti=ri−0.5λ,and the corresponding radius ri is a function of h∗ and is expressed as(21)ri=16h∗1/2.Figure 8
Comparison of the stretch rate evolution process between numerical result and the semiempirical model of Gao and Li [10].The stretch rate obtained by LBM agrees well with the theoretical analyze of Gao witht∗=0.5∼1.2 and suddenly increases at t∗=1.2∼1.5. 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 t∗=1.5∼2.0.


## 3.1. Laplace Law and Maxwell Construction

The 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. The 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 R0 is initialized at the computational domain center. Four different initial radii of droplets R0=40,50,60,70lu are chosen to validate Laplace’s law. In the C-S EOS, the corresponding temperature is chosen as Te=0.5,0.55,0.6,0.7,0.8,0.9,0.95Tc. The 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.Figure1 shows the comparison of the density coexistence curve obtained by LBM with different parameters ε and the Maxwell coexistence curve. The 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 Te<0.7Tc. However, when ε=0.11895, both the liquid and gas density agrees well with the theoretical value.Figure 1
Comparison of the LBM coexistence density with differentε with the Maxwell coexistence curve.Figure2 shows there is a linear relationship between the pressure difference Δp and the droplet radius 1/R at Te=0.5Tc, 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.Figure 2
The relationship of pressure difference and drop radius withTe=0.5Tc and different κs.


## 3.2. Droplet Impact on a Stationary Liquid Film

The computational domain is set as a 1001lu × 301lu rectangle region as shown in Figure 3, where R0 is the radius of the spherical droplet, H is the thickness of the liquid film, U0 is the initial droplet velocity, and Uf 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].Figure 3
Schematic of the computational domain for the droplet splashing on the moving film.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(12)Re=2R0U0νl,We=2ρlR0U02σ,h∗=H2R0,t∗=U0t2R0,u∗=U0Uf,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.The density field around the droplet is initialized as follows:(13)ρx,y=ρl+ρg2+ρl−ρg2tanh2x−x02+y−y02−R0w,where x0,y0 is the center of the droplet. And the initial density field of the liquid film is initialized as follows:(14)ρx,y=ρl+ρg2+ρl−ρg2tanh2y1−Hw,where y1 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.The dimensionless radius of the crownR∗, the upstream splash angle θu, the downstream splash angle θd, the dimensionless upstream jet height Hu∗, and the dimensionless upstream jet height Hd∗ are also introduced in the present study. All aforementioned parameters can be defined as(15)R∗=Rc2R0,Hu∗=Hu2R0,Hd∗=Hd2R0,where Rc is the radius of the crown and Hu and Hd are the height of the upstream and downstream jet as shown in Figure 4.Figure 4
Crown parameters for the droplet splashing on the thin film.A droplet impact on a stationary liquid film is introduced to verify the reliability of the proposed pseudopotential method. The radius of the spherical dropletR0 is 50lu, and the thickness of the liquid film is H=25lu. The temperature is chosen as Te=0.5Tc with corresponding liquid/gas density ratio ρl/ρg=720. The initial velocity U0 is set as 0. 125lu⋅tu−1, and the liquid film velocity is set as Uf=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. The 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.Figure 5
Snapshots of the droplet impact on the thin film at different time.Figure 6
Relationship between the dimensionless crown radiusR∗ and dimensionless time t∗.The relationship between the dimensionless crown radiusR∗ 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.13t∗, 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(16)R∗=βt∗+16h∗−13h∗−16h∗1/2,where β is a coefficient related to the energy loss factor λ and is calculated as(17)β=2λ23h∗1/4.The energy loss factorλs of the droplet impact on the stationary liquid film is determined by three-dimensionless numbers, including Re, We, and h∗:(18)λs=0.26Re−0.05We0.07h∗0.34.,The comparison between the LBM numerical result and the theoretical analysis is shown in Figure7. The numerical results are located between two lines with λ=0.42 and λ=0.78, 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 h∗. The same phenomenon is also observed in some experimental results [5] and numerical results [23, 24, 26, 31].Figure 7
Comparison of the crown radius evolution between numerical result and the semiempirical model of Gao.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:(19)St=t∗+ti+2τ0−1,where ti is the ending time of the deformation phase and can be calculated as(20)ti=ri−0.5λ,and the corresponding radius ri is a function of h∗ and is expressed as(21)ri=16h∗1/2.Figure 8
Comparison of the stretch rate evolution process between numerical result and the semiempirical model of Gao and Li [10].The stretch rate obtained by LBM agrees well with the theoretical analyze of Gao witht∗=0.5∼1.2 and suddenly increases at t∗=1.2∼1.5. 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 t∗=1.5∼2.0.


## 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.Table. 1
Simulation cases of a droplet impact on the moving liquid film.
Reynolds number (Re)Weber number (We)Dimensionless film thickness (h∗)Dimensionless liquid film velocity (u∗)100069.00.250.550010040100069.00.25139.41165.587.80.10.20.40.50.250.81.0The 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:(22)λf=0.53Re−0.02We0.03u∗−0.26h∗0.12.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 Figure9. 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. The 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 9
Droplet impact the moving thin film with differentRe at (a) t∗=0.75, (b) t∗=1.25, and (c) t∗=2.0 with We=87.8,h∗=0.25, and u∗=0.5.
(a)(b)(c)Figure10 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. The upstream crown height is higher than the downstream with the same Re 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 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. The 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. The 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 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.Figure 10
Evolution of the (a) upstream crown height, (b) downstream crown height, (c) crown radius with dimensionless timet∗ under different Re.
(a)(b)(c)


### 4.2. Weber Number

The 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.Figure 11
Droplet impact on the moving thin film with differentWe at (a) t∗=0.75, (b) t∗=1.25, and (c) t∗=2.0 with Re=500, h∗=0.25, and u∗=0.5.
(a)(b)(c)The evolution of dimensionless upstream heightHu∗, downstream height Hd∗, 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). The 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]. The 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.Figure 12
Evolution of the (a) upstream crown height, (b) downstream crown height, and (c) crown radius with dimensionless timet∗ under different We.
(a)(b)(c)


### 4.3. Film Thickness

As mentioned above, bothRe 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. The snapshots of crown evolution with different dimensionless liquid film thickness h∗=0.1,0.2,and0.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 13
Droplet impact on the moving thin film with differenth∗ at (a)h∗=0.1, (b) h∗=0.25, and (c) h∗=0.5 with Re=500, We=87.8, and u∗=0.5.
(a)(b)(c)Figure14 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. The 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]. The 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. The results also agree well with the studies of Mukherjee and Abraham [23] and Kharmiani et al. [31].Figure 14
Evolution of the (a) upstream crown height, (b) downstream crown height, and (c) crown radius with dimensionless timet∗ under different h∗.
(a)(b)(c)


### 4.4. Horizontal Velocity of the Liquid Film

The horizontal velocity of the liquid film is discussed in this part. Figure15 shows the crown evolution process of a drop impact on the moving liquid film with u∗=0.5and1. 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.5and1.0, 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].Figure 15
Droplet impact on the moving thin film with differentu∗ at (a) u∗=0.5 and (b) u∗=1.0 with Re=500, We=87.8, and h∗=0.25.
(a)(b)Figure16 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 16
Evolution of the (a) upstream crown height, (b) downstream crown height, and (c) crown radius with dimensionless timet∗ under different u∗.
(a)(b)(c)Figure17 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. 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.Figure 17
The liquid velocity field near the upstream crown with (a)u∗=0.0 and (b) u∗=0.5 and downstream crown with (c) u∗=0.0 and (d) u∗=0.5 at t∗=1.5.
(a)(b)(c)(d)


## 4.1. Reynolds Number

The snapshots of a droplet impact on the moving liquid surface are presented in Figure9. 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. The 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 9
Droplet impact the moving thin film with differentRe at (a) t∗=0.75, (b) t∗=1.25, and (c) t∗=2.0 with We=87.8,h∗=0.25, and u∗=0.5.
(a)(b)(c)Figure10 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. The upstream crown height is higher than the downstream with the same Re 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 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. The 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. The 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 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.Figure 10
Evolution of the (a) upstream crown height, (b) downstream crown height, (c) crown radius with dimensionless timet∗ under different Re.
(a)(b)(c)


## 4.2. Weber Number

The 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.Figure 11
Droplet impact on the moving thin film with differentWe at (a) t∗=0.75, (b) t∗=1.25, and (c) t∗=2.0 with Re=500, h∗=0.25, and u∗=0.5.
(a)(b)(c)The evolution of dimensionless upstream heightHu∗, downstream height Hd∗, 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). The 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]. The 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.Figure 12
Evolution of the (a) upstream crown height, (b) downstream crown height, and (c) crown radius with dimensionless timet∗ under different We.
(a)(b)(c)


## 4.3. Film Thickness

As mentioned above, bothRe 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. The snapshots of crown evolution with different dimensionless liquid film thickness h∗=0.1,0.2,and0.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 13
Droplet impact on the moving thin film with differenth∗ at (a)h∗=0.1, (b) h∗=0.25, and (c) h∗=0.5 with Re=500, We=87.8, and u∗=0.5.
(a)(b)(c)Figure14 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. The 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]. The 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. The results also agree well with the studies of Mukherjee and Abraham [23] and Kharmiani et al. [31].Figure 14
Evolution of the (a) upstream crown height, (b) downstream crown height, and (c) crown radius with dimensionless timet∗ under different h∗.
(a)(b)(c)


## 4.4. Horizontal Velocity of the Liquid Film

The horizontal velocity of the liquid film is discussed in this part. Figure15 shows the crown evolution process of a drop impact on the moving liquid film with u∗=0.5and1. 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.5and1.0, 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].Figure 15
Droplet impact on the moving thin film with differentu∗ at (a) u∗=0.5 and (b) u∗=1.0 with Re=500, We=87.8, and h∗=0.25.
(a)(b)Figure16 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 16
Evolution of the (a) upstream crown height, (b) downstream crown height, and (c) crown radius with dimensionless timet∗ under different u∗.
(a)(b)(c)Figure17 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. 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.Figure 17
The liquid velocity field near the upstream crown with (a)u∗=0.0 and (b) u∗=0.5 and downstream crown with (c) u∗=0.0 and (d) u∗=0.5 at t∗=1.5.
(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 ofRe number, We 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 ofRe and We, 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. Whenh∗=0.25, 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.



---

*Source: 1801639-2020-07-15.xml*