Study of droplet splashing on a liquid film with a tunable surface tension pseudopotential lattice Boltzmann method

A tunable surface tension pseudopotential lattice Boltzmann method (LBM) with an axisymmetric boundary is applied to study a droplet splashing on a thin film. Our work focused on the crown behavior influenced by the different parameters, including the Reynolds number, Weber number, liquid film thickness, and gravity acceleration. In addition, the total kinetic energy of the crown is proposed to explain the evolution of the crown shape from the perspective of energy. Based on the achieved results, it is found that the crown radius changes with time are unaffected by these parameters and consistent with the power-law in previous studies. However, the height of the crown and the satellite droplet formation process are affected by the referred parameters. Besides, it is found that the energy consumption during the collapse process and splash angle are two key factors affecting the height of the crown. There is a threshold liquid film thickness with the maximum height of the crown, which is between 0.2 and 0.25 times the droplet radius. This study shows that the proposed LBM pseudopotential model is a robust and effective tool for the study of the droplet splashing on the thin film and has the potential to predict the droplet splashing phenomena in the presence of complex boundaries.A tunable surface tension pseudopotential lattice Boltzmann method (LBM) with an axisymmetric boundary is applied to study a droplet splashing on a thin film. Our work focused on the crown behavior influenced by the different parameters, including the Reynolds number, Weber number, liquid film thickness, and gravity acceleration. In addition, the total kinetic energy of the crown is proposed to explain the evolution of the crown shape from the perspective of energy. Based on the achieved results, it is found that the crown radius changes with time are unaffected by these parameters and consistent with the power-law in previous studies. However, the height of the crown and the satellite droplet formation process are affected by the referred parameters. Besides, it is found that the energy consumption during the collapse process and splash angle are two key factors affecting the height of the crown. There is a threshold liquid film thickness with the maximum height of the crown, which is between 0.2 and 0.2...


I. INTRODUCTION
The splashing phenomenon of liquid droplet impinging on the liquid film is widespread in industrial and engineering processes, such as the injection in the turbo engine to improve the combustion efficiency, the atomization in the hydraulic engineering during the flood discharge process, and the spraying process of the ink jet printer. The splashing shape of this complex multiphase phenomenon is affected by many factors, including the gas-liquid density ratio, viscosity coefficient ratio, impact velocity, surface tension coefficient, etc. [1][2][3] Due to the complexity of the gas-liquid interface and the diversity of the evolution process of the gas-liquid interface in the process of droplet impact, the dynamics of this complex multiphase phenomenon is still not clear.
The flow and heat transfer of a liquid film on different surfaces can be obtained by theoretical analysis, 4-7 but the phenomenon of droplet splashing on a liquid film with complex boundary conditions needs many simplified assumptions, which are difficult to be studied by theoretical analysis. Earlier researchers have made a lot of efforts to understand the physics of the interaction between the drop and liquid film at different length scales. [1][2][3]8,9 Experimental results of the evolution process have been derived with a droplet splashing on the liquid film under the influence of different parameters, which is being captured by a high-speed camera. These experiments mainly focus on how dimensionless parameters affect the shape of the liquid film, such as the Weber number, Reynolds number, and film thickness. 10,11 Because of the impact of the different parameters, the evolution process can be divided into fusion, rebound, partial rebound, crown, and splash. 12 With the development of high-speed camera and image processing technology, the experiment of liquid crown evolution under different experiment conditions [13][14][15] and the process of liquid film oscillation can further be accurately captured. 16,17 Although the experiment can show the evolution process of droplet splashing on the liquid film directly, it cannot provide physical insight.
With the development of modern computers, computational fluid dynamics (CFD) has become an important tool to study the ARTICLE scitation.org/journal/adv evolution process of droplet splashing on the liquid film. Numerical simulation can reveal the evolution of pressure and velocity fields in the surrounding gas during this splashing process. Because of the dramatic changes in the gas-liquid interface and the large density gradient at this interface, numerical simulation is still a challenge. Macroscopic numerical simulation methods are mainly based on the Navier-Stokes and interface tracking methods, [18][19][20] such as the volume of fluid (VOF) method 18,19 and level-set method. 20 The flow field evolution during the impact process can be obtained by a numerical study, and the deformation of the crown and satellite droplet formation are analyzed from the perspective of velocity and vorticity field. These results indicated that the crown instability is the main reason for the crown deformation and the satellite droplet formation, but there are some discrepancies in these studies. [21][22][23][24][25][26] Allen 21 indicated that the Rayleigh-Taylor instability is the main reason for satellite droplet formation, but Fullan, 22 Roisman, 23 and Nikolopoulos 24,25 proposed that the formation of the satellite drop is due to the Plateau-Rayleigh instability, while Krechtnikov 26 reported that the Richtmyer-Meshkov instability is the main reason. As a mesoscopic method, the lattice Boltzmann method (LBM) has become an efficient and powerful tool to study a complex multiphase phenomenon based on Boltzmann kinetic molecular dynamics, such as cavitation 27 and crystallization. 28 Compared with conventional macroscope numerical methods based on Navier-Stokes equations, the LBM has many advantages. The LBM is based on a series of linear equations, but the convective terms in Navier-Stokes equations are nonlinear. In addition, the pressure field can be calculated using the equation of state (EOS), which does not need to solve the Poisson equation. Furthermore, complex boundaries can be dealt with the LBM easily. 29 After 30 years of development, a series of LBM multiphase flow models have been established; all of them can be categorized as the color gradient model, 30 the pseudopotential model, 31,32 the free energy model, 33,34 and the

Reynolds
Weber Dimensionless Gravity number (Re) number (We) film thickness (h * ) acceleration (g) interface tracking model. 35 In recent years, the LBM gradually occupies a competitive place in the field of numerical studies of a droplet splashing on the thin film. [36][37][38][39][40][41][42][43] By improving the phasefield model proposed by He, 35 Lee 36 simulated the droplet splashing process with a density ratio of 1000. Mukherjee 37  extended this model to three dimensions to simulate the droplet impacting the moving liquid film. However, all the above simulations are based on the large density phase-field model proposed by Lee. 36 In this model, the collision step is divided into a pre-streaming collision step and post-streaming collision step. In addition, the gradient and Laplacian need to be solved in the calculation process, which requires a lot of computing resources. After that, Liu 39 used this model to study the process of droplet splashing on a liquid film with horizontal initial velocity. Recently, Wang et al. 40 and Li et al. 41 have built a multiphase Lattice Boltzmann flux solver, which can be used to simulate compressible flow with a high density ratio with a body-fitted grid, and the droplet splashing on the thin film is also simulated successfully. The pseudopotential model was proposed by Shan and Chen 31 in 1993. This model is widely used because of its simplicity and computational efficiency. The model uses the interaction between particles to form an interface automatically, avoiding the use of the interface tracking method. Li 42 proposed a high density ratio pseudopotential model, which satisfies the thermodynamic consistency and can adjust the surface tension independently, and this model was used to simulate the process of droplet splashing on the liquid film on different Reynolds numbers and Weber numbers. Recently, Kharmiani 43 studied the effects of the gas-liquid viscosity ratio, surface tension, and density ratio on the process of droplet impacting on a liquid film with this modified pseudopotential model. In our study, this modified pseudopotential LBM is extended to an axisymmetric coordinate system and a droplet splashing on a thin film with a different Reynolds number, Weber number, gravity, and liquid film is studied. Furthermore, the influence of the total kinetic energy (TKE) of the liquid crown on the crown height and the satellite droplet formation is studied. This paper is organized as follows: The large-density-ratio multi-relaxation LBM pseudopotential model is described in Sec. II. The evolution process of droplet splashing on the liquid film under different conditions is studied, and the results are presented in Sec. III. Finally, our conclusions are given in Sec. IV.

II. MODEL DESCRIPTION
The pseudopotential model with a BGK collision operator has poor numerical stability and large spurious currents when simulating the multiphase flow phenomenon with a large density ratio. To achieve the numerical stability and small spurious currents with a high-density ratio multiphase phenomenon, an MRT collision operator is applied in the present study, and the corresponding particle distribution equation with extra terms is 42 where fi is the particle distribution function, f eq j 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 D2Q9 lattice applied in the present study with the discrete velocities is given as 29 where c = Δx/Δt is the lattice constant, and the diagonal matrix Λ can be expressed as 42 The relaxation time τν is related to the kinematic viscosity ν = 1/[c 2 s (τν−0.5)Δt], and cs is the lattice sound speed. M is an orthogonal transformation matrix, and M −1 is the inverse matrix of M. The equilibrium moment m eq can be obtained by projecting the equilibrium distribution f eq j onto the moment space with m eq = Mf eq , which can be described as 44 The modified external forcing scheme proposed by Li 42 is applied in the present study, which is expressed as Here, ε is a parameter that can be adjusted to achieve thermodynamic consistency. In addition, ε is chosen as 0.112 in the present study. Fm is the interaction force between particles, which can be obtained by 42 where G is the interaction strength between two particles, ψ is the interparticle potential, and ωi are the weights. For the D2Q9 LBM model, ω 1−4 = 1/3 and ω5− 8 = 1/12. 45 The ψ calculated after a nonideal-gas EOS is introduced as 46 where p EOS is the pressure calculated using a non-ideal-gas EOS and the Carnahan-Starling (C-S) EOS is applied in our study, 46 where a = 0.4963R 2 g T 2 c /pc and b = 0.1873RgTc/pc, Tc is the critical temperature, and pc is the critical pressure, and the parameters are chosen as a = 0.25, b = 4, and Rg = 1.
A source term C is applied in the present study to achieve the tunable surface tension, which is not dependent on the density ratio. Here, the source term C is expressed as 42 and Q can be obtained by 42 where κ is a parameter to tune the surface tension and its value is between 0 and 1. Finally, the collision process can be expressed as and the streaming process can be expressed as Finally, the macroscopic density ρ and macroscopic velocity u are obtained by where F = Fm + G + ⋅ ⋅ ⋅ is the total force acting on the fluid particles, and the gravity force G can be obtained by where g is the gravity acceleration and ρg is the initial gas density.

III. MODEL VALIDATION
The evaluation of the Laplace law, Δp = σ/R, where σ is the surface tension and R is the droplet radius after the fluid system reaches an equilibrium state, is used to validate the proposed method. The computational domain is a 401lu × 401lu square region, periodic boundary conditions are applied at all the boundaries of the computational domain, and a stationary droplet of radius R 0 is located in the center of the computational domain. A series of initial radii of droplets R 0 = 40, 50, 60, 70, and 80lu are chosen to validate Laplace's law. In the C-S EOS, the corresponding temperature is chosen as Te = 0.5Tc, with density ratio as ρ l /ρg = 720. The relaxation time is chosen as τν = 0.5375, and a series of κ = 0.0, 0.2, 0.5, and 0.9 are selected to tune the surface tension. Figure 1 shows the density and velocity distributions after the fluid system reaches equilibrium with ρ l /ρg = 720, R 0 = 50lu, and κ = 0.5. The spurious currents are mainly distributed in the gas phase, with the maximum spurious current in the flow field being only 0.016lu ⋅ tu −1 . As shown in Fig. 2, a linear relationship is observed between the pressure difference Δp and the reciprocal of the droplet radius, 1/R, when the fluid system reaches equilibrium. The numerical results agree well with the theoretical analysis, with density ratios ρ l /ρg = 720 and different κ's, and the surface tension decreases with the increase in κ.

IV. NUMERICAL RESULTS AND DISCUSSION
The proposed pseudopotential LBM is used to study a droplet splashing on the thin film. Several dimensionless parameters governing the droplet splashing on the liquid film are introduced in our study, including 9 where Re is the Reynolds number, We is the Weber number, h * is the dimensionless film thickness, and t * is the dimensionless time.
The TKE of the crown and satellite is analyzed from the energy perspective, and the TKE is defined as 47 If there is no special statement, the unit of mass discussed in this paper is mu, the unit of length is lattice length lu, the unit of time is ts, and the unit of velocity is lu ⋅ ts −1 .
The computational domain is a 501lu × 301lu rectangular region, as shown in Fig. 3. An axisymmetric boundary condition is applied at the left boundary, nonequilibrium extrapolation conditions are applied at the right boundary, and the no-slip boundary conditions are applied at both the top and bottom boundaries. A spherical droplet with radius R 0 = 50lu is initialized in the computational domain, H is the thickness of the liquid film, and the initial velocity of the droplet is U 0 .
The density field around the droplet is initialized as follows: where (x 0 , y 0 ) is the center of the droplet. In addition, 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, the initial interface width w in Eqs. (15) and (16)  Several dimensionless parameters are introduced in our study to analyze the influence of governing parameters, including the dimensionless radius of the crown R * , the splash angle θ, and the dimensionless crown height H * , which can be expressed as follows:

ARTICLE scitation.org/journal/adv
where Rc is the radius of the crown and Hc is the height of the crown, as shown in Fig. 3(b). In the present study, the initial velocity U 0 is set as 0.125lu ⋅ tu −1 and different Reynolds numbers and Weber numbers are obtained by adjusting τν and κ; the simulation cases are shown in Table I.

A. Effects of the Reynolds number
The snapshots of the evolution of a droplet splashing on the liquid surface with different Re are shown in Fig. 4. When the droplet is splashing on the liquid film, a jet will form at the neck between the droplet and the liquid film, and then, the jet forms a crown (t * = 0.75). Then, the height of the crown increases and its thickness decreases (t * = 1.25). With the development of the crown, under the condition of Re = 500 and Re = 2000, the crown will break up and satellite drops will be formed under the influence of the Plateau-Rayleigh instability with Re = 500 and Re = 1000 (t * = 2.0). With the increase in Re, which means that the viscosity of the liquid decreases, the crown thickness decreases and the crown height increases, and the satellite droplets are easier to form. However, there is no obvious crown formed with Re = 40 because the liquid viscosity is too large.
The evolution of dimensionless height H * and radius R * with t * is shown in Fig. 5. With the increase in the Re number, the height of crown increases, but the radius of the crown does not change with the Re number. According to the studies of Yarin, 3 there is a linear relationship between the dimensionless crown radius R * and dimensionless time √ t * . This linear relationship between the dimensionless crown, 43 radius R * , and dimensionless time √ t * is also proved by other previous studies, and the best curve fitted on our results is R * = 1.39 √ t * . Figure 6 shows the change in the TKE of the crown with dimensionless time t * under different Re, and the maximum TKE of the crown increases with the increase in Re, which means that the initial TKE of the droplet consumed in the impact process of a high Re number is less than that of the low Re number. In addition, the difference between the TKE of the crown with Re = 500 and Re = 1000 is small, and the satellite drops formed under both conditions indicate that the less TKE is consumed in the impact process and the crown is easier to break up with the same surface tension.

B. Effects of the Weber number
The different We can be obtained by adjusting the surface tension σ through the parameter κ. The influence of We is discussed in this part. High We number means low surface tension, and the crown is more prone to deformation and break up. The snapshots of the evolution of a droplet splashing on the liquid surface with different We are shown in Fig. 7. When dimensionless time t * = 0.75, the satellite droplet has formed with We = 1165.5. With the evolution of the crown, the crown height increases, while the crown thickness decreases and the Plateau-Rayleigh instability of the crown also increases. When dimensionless time t * = 2.0, the satellite droplet is formed with We = 87.8, 139.4, and 1165.5, while the crown remains stable with We = 69.0.
The evolution of dimensionless height H * and radius R * with t * under different We numbers is shown in Fig. 8. With the increase in the We number, the height of the crown increases and the radius of the liquid crown does not change with the We number. The linear relationship between the dimensionless crown radius R * and dimensionless time √ t * is R * = 1.33 √ t * . Figure 9 shows the change in the TKE of the crown with dimensionless time t * under different We; the maximum TKEs under different cases are almost identical while We = 67.0, 87.8, and 139.4. However, the maximum TKE with We = 1165.5 is smaller than that in the other three cases. Because

ARTICLE
scitation.org/journal/adv of the non-physical phenomena of the pseudopotential model, the smaller droplet is "eaten". 29

C. Effects of the liquid film thickness
The influence of the liquid film thickness has been studied in previous studies. The liquid film dimensionless thickness varies from 0.1 to 0.5 in our work. The snapshots of crown evolution with different dimensionless liquid film thicknesses h * = 0.1, 0.25, and 0.5 are shown in Fig. 10; with the decrease in the liquid film, the crown is more likely to break up and the breaking time decreases. The crown breaks up at t * = 1.125 while h * = 0.1, while the crown breaks at t * = 1.5 while h * = 0.25. However, when h * increases to 0.5, the satellites do not form, and the numerical results are consistent with the numerical results of Mukherjee. 37 Figure 11 shows the relationship between the crown height and radius with liquid film thickness. When the liquid film thickness increases from 0.1 to 0.25, the crown height increases rapidly, but when the liquid film thickness increases from 0.25 to 0.5, the crown height decreases. The crown radius is consistent with different dimensionless film thicknesses except for h * = 0.1. The possible explanation for the increase in the dimensionless height of the crown H * with the dimension film thickness h * from 0.1 to 0.25 is that the growth of the crown needs the liquid film with enough thickness. Furthermore, with the increase in liquid film dimensionless thickness from 0.1 to 0.25, the liquid film can prevent the vertical momentum from changing to the radical momentum and the diffusion process of the radical momentum. The linear relationship between the dimensionless crown radius R * and dimensionless time √ t * is R * = 1.39 √ t * with h * = 0.2, 0.25, 0.3, 0.4, and 0.5. Figure 12 shows that the TKE of the crown changes with time under different liquid film thicknesses. The TKE of the crown decreases with the increase in liquid film thickness, and the more the thickness of the liquid film, the more the consumption of the TKE when transforming the vertical momentum into radical momentum during collision. With the decrease in the maximum TKE of the liquid crown, the lower the Plateau-Rayleigh instability of the crown is, the less likely the crown is to break up.
The maximum dimensionless crown height and splash angle that vary with time are shown in Fig. 13. When the dimensionless liquid film thickness increases from 0.1 to 0.2, the splash angle increases from 32 ○ to 50 ○ and the droplet splash dimensionless height also increases rapidly from 0.99 to 1.51. However, with the increase in the film thickness from 0.2 to 0.5, the droplet splashing angle continues to increase from 50 ○ to 63 ○ , but the increasing speed becomes slower and the dimensionless film height decreases from 1.52 to 1.31.
Furthermore, the influence of film thickness on crown height is analyzed through the velocity field of the jet during the initial stage of crown formation. The velocity field with different dimensionless liquid film heights h * = 0.1, 0.25, and 0.5 is shown in Fig. 14. As mentioned earlier, the TKE for crown formation increases with the decrease in dimensionless liquid film thickness, but the splash angle, as shown in Fig. 14(a), is smaller with h * = 0.25 and 0.5, which leads to smaller crown height. However, when the thickness of the liquid film h * increases from 0.25 to 0.5, the cushion of the liquid film is more obvious, and the energy absorbed by the liquid film during collision is more. Compared with Figs. 14(a) and 14(b), it can be seen that the jet velocity decreases with the increase in the thickness of the liquid film even when the splash angle is very close.

D. Effects of the gravity
Furthermore, the influence of gravity is also discussed in our study and four different gravity accelerations are selected for simulation. The snapshots of crown evolution with dimensionless time t * = 0.75, 1.25, and 2.0 are shown in Fig. 15.
In the initial stage, the drop height decreases with the increase in g (t * = 0.75), which means that with a larger g, the vertical momentum changes into the radical momentum more rapidly and the vertical and horizontal velocity of the neck jet is larger. However, under the action of larger g, the vertical velocity decreases. The crown height decreases with the increase in g (t * = 1.25), but the crown radius still increases with the increase in g (t * = 1.25). When t * increases to 2.0, satellite droplets are formed in all four cases. With the decrease in g, the Plateau-Rayleigh instability is stronger and more satellite droplets are generated under the conditions of g = 0 and 10 −4 lu ⋅ tu −2 .
The relationship between crown height and radius with the gravity acceleration is shown in Fig. 16. The results indicate that the height of the crown decreases with the increase in g, but the crown radius increases with the increase in g. With different gravity acceleration, the crown radius has a linear relationship with dimensionless time √ t * , but the coefficients are different. Due to gravity acceleration g, the TKE consumed in the collision process becomes smaller and more TKEs of droplets are converted into the TKE of the crown, as shown in Fig. 17.

V. CONCLUSIONS
A tunable surface tension LBM pseudopotential model has been used to study a droplet splashing on the thin film. The results demonstrate the robustness and accuracy of the proposed method. The pseudopotential model has been used to study the influence of the Re number, We number, interface thickness, and gravity acceleration on the liquid crown. The following conclusions can be drawn: 1. The evolution of crown radius is not affected by the Re number, but the crown height is affected. The crown height increases with the increase in the Re number. When Re is smaller, the energy consumption caused by the viscosity of the liquid in the collision process becomes greater, which leads to a smaller crown. No satellite droplets are formed until Re ≥ 500. 2. The height of the crown increases with the increase in the We number, but the radius of the crown does not change. When the We number changes from 69 to 139.4, the TKE of the liquid crown does not change. When the We number increases to 1165.5, the TKE of the liquid crown decreases, which is due to the defect of the LBM pseudopotential model; the small broken satellite droplets are "eaten." 3. The splash angle and the energy consumption of the liquid film during the collision are two key factors that influence the crown radius and height. When the dimensionless thickness of the liquid film is between 0.2 and 0.5, the difference between crown radii is small, but the splash angle decreases significantly when the dimensionless thickness of the liquid film is 0.1, which leads the radius of the crown to become larger. The height of the crown increases first and then decreases with the increase in the thickness of the liquid film from 0.1 to 0.5, and the maximum height of the crown can be obtained when the thickness of the liquid film changes from 0.2 to 0.25. 4. The height of the crown decreases with the increase in the gravity acceleration, but the radius of the crown increases with the increase in the gravity. The smaller the gravity acceleration, the easier the formation of the satellite.
This study shows that the proposed LBM pseudopotential model is a robust and effective tool for the study of the droplet splashing on the thin film and has the potential to predict droplet splashing phenomena under complex conditions such as coating, painting, fuel injection and atomization of flood discharge of the high dam, and so on.