LATTICE BOLTZMANN PSEUDO-POTENTIAL MODELLING OF MULTIPHASE DROPLET PHENOMENA

The multiphase modeling of a droplet in a multiphase system is considered becoming a fundamental problem in fluid dynamics. A complex understanding of droplet behavior is critical to reveal a deeper insight into a more complex multiphase system. Droplet behavior studies are necessary to obtain a better understanding of solving multiphase problems in both the science and industrial aspect. The droplet behavior is characterized by a non-dimensional number such as the Eötvös number. In this study, numerical simulation was performed using the Lattice Boltzmann method. Parametric studies of Eötvös number was done. The parametric study of the Eo number is obtained using LBM. Based on the results obtained, it is concluded that the gravitational force influences the downwards rate of the droplet. Furthermore, the shape of the droplet during falling was depended on the Eo number as well. The higher Eo number means higher gravitational force, hence the velocity of the droplet is increasing as well as the reaction force of surface tension. This study is beneficial to give a deeper explanation of multiphase phenomena as well as contribute to the modeling of multiphase phenomena using an alternative numerical method of LBM. Keyword: Lattice Boltzmann method, Numerical, Droplet, Eötvös number


Introduction
The multiphase modeling of a droplet in a multiphase system is considered becoming a fundamental problem in fluid dynamics [1]. A complex understanding of droplet behavior is critical to reveal a deeper insight into a more complex multiphase system such as spray cooling [2,3], nuclear waste management [4,5], soil-vapor extraction [6][7][8], ink-jet printing [9,10], and droplet spreading related to covid-19 phenomena [11]. The droplet behavior is characterized by a non-dimensional number such as the Eötvös number which is defined as the ratio of gravitational force towards the interfacial tension. This dimensionless number represents the shape of a droplet moving in both certain surrounding fluid and gravitational force effects. Droplet behavior studies are necessary to obtain a better understanding of solving multiphase problems in both the science and industrial aspect.
In general, the parametric study based on the Eötvös number can be done by analytical approach, experimental apparatus as well as numerical simulation. However, the first approximation of an analytical solution is only limited by a simple problem. Furthermore, the experimental procedure provides actual data and conditions. Previous experimental studies have been conducted by Roisman et al [12], Bird et al [13], and Farhangi et al [14]. Those studies were clarified that droplet behavior plays an important role in the multiphase dynamics analysis. Unfortunately, due to its complex setup and costly, some experimental approach is considered to be difficult to overcome. Various numerical investigations have also been conducted by previous studies. Bussmann et al [15] applied the volume-of-fluid (VOF) method to simulate the splashing droplet impact phenomena. It was obtained that the contact angle was an important parameter in the development of splash phenomena. Tasoglu et al [16] simulated the impact of the droplet by using the front-tracking method. The demonstration revealed that droplet inner deformation was increasing in accordance with the Reynolds number and surface tension ratio of air to the droplet. Zheng et al [17] recently found that the finite-volume-based method was able to simulate the droplet impact as well as temperature parameter. Mehdi-Nejad et al [18] used the VOF model to observe the effect of air bubble entrapment on the droplet impact.
A lot of work has been performed on the droplet impact investigation, but some issues are not addressed yet, for instance, the effect of impact velocity, gravitational force, and droplet impact behavior are not considered. In this study, the numerical investigation of droplet impact as well as the gravitational force is performed by using a parametric study of Eötvös number. Furthermore, the Lattice Boltzmann method (LBM) as an alternative method considerably more efficient numerical method than Navier-Stokes (N-S) equation solver has been used in several studies [19]. The numerical simulations in this paper were carried out using the pseudo-potential lattice Boltzmann method (LBM) that has emerged as a useful approach for computational fluid dynamics analysis. The LBM uses a discretization based on the mesoscopic kinetic equation. It relatively is easier to handle multiphase flow which is difficult to be solved in the conventional method based on the N-S equation. The rest of the paper is organized as follows: Sec. 2, the brief of LBM methodology is presented. In section 3, we give the numerical results of the droplet simulation as well as the parametric study of the Eötvös number. A brief of trademark conclusion is given in section 4.

Method
The lattice Boltzmann method emerges as a powerful tool for numerical simulations for this decade [20]. One of the most efficient LBM methods is the pseudopotential Shan-Chen method [21]. The Shan-Chen LBM considers a behavior between particles and their neighbor. This method uses pseudopotential forces as a form of particle interaction. The fluid motion is represented by a distribution of particles stated by using distribution functions. Those particles are experiencing a collision and streaming that is discretized based on Bhatnagar-Gross-Krook (BGK) collision operator as shown in Eq (1). The left operator indicates the streaming terms, besides the right one is the collision operator Where, , is the function of density distribution, , is the equilibrium distribution function, is the position, is the time, is the lattice speed ( x/∆ ) and denotes the lattice direction. The fluid density ( ) and velocity ( ) are calculated by using Eqs. (3) and (4). After the collision, the particles will reach an equilibrium state within a certain time which is noted as relaxation time ( ). The equilibrium distribution function is shown in Eq (2). is related to the viscosity as = 2 − 0.5 ∆ .
is the lattice sound speed where sets as 2 3 . The grid spacing ( x) and timesteps (∆ ) are set to 1. The simulation uses the 3D model of D3Q19 with 19 velocity directions at a given point (Fig. 1). The discrete velocity is shown below: In the Shan-Chen model, the fluid interaction force causes the phase separation and is stated in Eq (5). The terms of denotes the component of fluid. Where G is the parameter that control the strength of inter-particle force. is the mean-field potential that is shown in Eq 6. Considering the Shan-Chen forces, the value of fluid velocity (Eq 4) is modified shown in Eq 7. In this study, the relaxation time of both fluid is assumed to be equal. 0 is the normalization constant which is adjusted as 1.
In this model, the half bounce-back boundary condition is chosen at the top and the bottom walls. The right, left, front and back walls are using periodic boundary condition. The simulation results are represented in dimensionless quantities In this paper, the dimensionless parameters that will be analyzed is Eötvös number (Eo). Here, the non-dimensional parameters related to the development of droplet are defined, as follows: Where, is the droplet radius, body force of gravitational acceleration is represented by , 0 is the impact velocity of droplet, relates to the droplet viscosity, is the droplet density and denotes the surface tension between droplet and surrounding fluid.
The numerical simulation is performed by using C programming language. The code is calculated using open-source compiler, Codeblock. Meanwhile, the post-processing process was done by using open source software, Paraview. The typical numerical calculation was about 4 hours of computational time for each case.

Result and Discussion
The Eötvös number is measuring the effect of gravitational forces compared to the droplet surface tension. It is used to determine the shape of droplets moving in a certain surrounding fluid. In this study, the radius of droplet sets to 12 with fixed Ohnesorge Number (Oh) = 0.1. Dimensionless Ohnesorge (Oh) number relates to the ratio of internal viscous dissipation to the surface tension. The lower Oh means the weaker friction force due to the viscous force. The Oh can be written in terms of the square root of the Weber number to the Reynolds number. Hence, by assuming a fixed value of Oh, the effect of the droplet is only influenced by gravitational force in terms of Eo. Simulation domain sets to 151x301x80 within the 3Dlattice domain. The numerical variations are shown in Table 1. The simulation results for each case are shown in Fig.1. Based on those results, can be seen that the shape of the droplet for case A is perfectly symmetrically rounded. However, for case B and C, the shape of droplets begin to be deformed as an oval form. This phenomenon relates to the previous studies [1,19]. The droplet shape changed is caused by the interaction between gravitational force and surface tension. Hence, for each Eo number is obtained distinction results of droplet form. The gravitational force which is perpendicular with surface tension produces the action-reaction force. Further, the gravitational force leads the droplet to move downwards at a certain velocity and acceleration. In linear, the higher gravitational force yields a higher value of force action-reaction towards surface tension as stated by the 3rd law of Newton. The droplet that is pressed by a couple of perpendicular forces exerts the water into the radial axis. Thus, it produces the oval form of a droplet. It also can be found that weaker Eo will delay the droplet falling and remains its form. Here, it can be seen that Case C will rapidly go downwards than both cases A and B because of the higher value of gravitational force as well as Eo number. The droplet will collide to the surface in terms of droplet impact. The droplet will oscillate up to attaining the equilibrium state. After impact, the droplet surface will be flattened then back to its original rounded form. The oscillation phenomena occur due to the presence of both gravitational force and surface tension in the form of action-reaction force interaction. In the case of a smaller value of gravitational force, the action-reaction force will be decreased as well. This phenomenon causes oscillation of droplets will take a longer duration before reaching the equilibrium state. In contrast, the equilibrium state is quickly attained for a higher value of gravitational force. The simulation result shows that cases A, B, and C will reach the equilibrium state at t = 7.42t, 6.79t, 6.35t, respectively.

Conclusion
In this study, the numerical investigation of droplet behavior, as well as its impact on the surface, is performed. The parametric study of the Eo number is obtained using LBM. Based on the results obtained, it is concluded that the gravitational force influences the downwards rate of the droplet. Furthermore, the shape of the droplet during falling was depended on the Eo number as well. The higher Eo number means higher gravitational force, hence the velocity of the droplet is increasing as well as the reaction force of surface tension. During impact, the droplet will fluctuate to attain the equilibrium state that is also influenced by the Eo number.