CHEMICAL ENGINEERING TRANSACTIONS VOL. 57, 2017 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš, Laura Piazza, Serafim Bakalis Copyright © 2017, AIDIC Servizi S.r.l. ISBN 978-88-95608- 48-8; ISSN 2283-9216 Study of the Motion of a Droplet in a Microchannel using Shan-Chen Multiphase Lattice Boltzmann Model Bagdagul Dauyeshova*a, Luis R. Rojas-Solórzanob, Ernesto Monacoc a,bNazarbayev University, School of Engineering, Kabanbay batyr ave., 53, 010000, Astana, Kazakhstan c Friedrich-Alexander Erlangen-Nurnberg University, Dept. of Mathematics, Cauerstrasse, 11, 91058, Erlangen, Germany bagdagul.dauyeshova@nu.edu.kz This paper presents a study of the motion of a deformable droplet modeled using Shan-Chen (SC) Multiphase Lattice Boltzmann Model (LBM). Firstly, a 2D static droplet in a periodic domain is modeled such that its surface tension is numerically determined and successfully validated against Young-Laplace equation. Next, simulation of a 2D droplet transported by a gas stream along a periodic microchannel is carried out. Initially, the droplet is placed close to the top wall and surrounded by a gas stream which is originally flowing in fully developed condition, as described by the classical plane Poiseuille flow. Both the droplet and the carrier fluid have same kinematic viscosities but different densities. The simulation is extended until pseudo-steady state, determined by a constant vertical position and normalized velocity, is reached. The results of this study are successfully validated using available data from previous highly-accurate numerical work. The LBM demonstrated a number of major advantages such as easiness to set up, ability to track the phase interface automatically and its inherent parallel capacity which helped reducing the computational time. However, numerical results obtained with the SC-LBM, which is presented in this study, demonstrate that the LBM also has limitations in handling large density and viscosity ratios between the phases. 1. Introduction Multiphase flow in porous media are found in many applications, such as geological sequestration of CO2, underground pollutant remedy and proton exchange membrane fuel cells (PMEFC), to name a few. Numerical techniques can greatly help when facing these problems. Each proposed scheme must be accurately validated against simple flows which are basis of more complex phenomena. A very common problem found in oil and gas, as well as in many other applications of multiphase flow transport, is the flow of droplets with the gas stream in a channel. The study of droplet distribution through the cross-section of the channel is of interest as this defines the overall flow rate of the multiphase fluid (Mortazavi and Tryggvason, 2000). Typical approaches are based on discretizing the Navier-Stokes equations, though that poses a series of challenges (Chen et al., 2014) mainly due to the nonlinear convective term. The Lattice Boltzmann Method (LBM hereafter) has arisen as a meso-scale alternative to standard CFD techniques capable of dealing with the complexity of micro-scales correctly, while reproducing the typical macro-scales characterizing a significant number of different flows. The Shan-Chen (SC) model is a multiphase LBM that has been successfully used in numerous studies (Chin et al., 2002; Kang et al., 2002; Zhang et al., 2001). Despite the fact, that LBM has shown the potential to be applied successfully for simulation of multiphase flow in porous media, the method has still issues that need to be addressed, the most notable being the challenge of treating large density and viscosity ratios between the fluid phases present in the flow (Chen et al., 2014). Indeed, all the different implementations of multiphase LBM suffer when applied to simulation of multiphase flows with high density ratios, as this leads to an increase in spurious currents (unphysical velocities at the vicinity of interfaces due to the truncation error in evaluating the numerical gradient) (Chen et al., 2014). Strong unphysical velocities are undesired since they can get as big as the typical velocity scales of the considered problem, shadowing the relevant physics to be observed. DOI: 10.3303/CET1757212 Please cite this article as: Dauyeshova B., Rojas-Solorzano L., Monaco E., 2017, Study of the motion of a droplet in a microchannel using shan-chen multiphase lattice boltzmann model, Chemical Engineering Transactions, 57, 1267-1272 DOI: 10.3303/CET1757212 1267 In this work, the droplet motion along the straight channel is studied using the SC LBM. The following sections of the paper are arranged as follows: firstly, a description of the LBM and of the SC model is provided; then the static droplet case is examined in order to find out the surface tension and to validate the interfacial results with the Young-Laplace equation. Next, the motion of the 2D droplet is simulated and its lateral position and axial velocity are analysed and compared against the data from Mortazavi and Tryggvason (2000). Finally, conclusions and directions for future work are presented. 2. Lattice Boltzmann Method (LBM) The Boltzmann equation describes the evolution of the fluid state because of streaming and collision of its molecules. Since the Navier-Stokes equations can be obtained by a multiscale expansion of the Boltzmann equation, one may think of solving a simplified form of that as an alternative to standard numerical approaches like Finite Elements or Finite Volume techniques. The LBM is obtained by a convenient discretization of the phase space of the Boltzmann equation where the infinite set of trajectories and molecular velocities are replaced by a finite set of speed vectors (the lattice) (Chen et al., 2014; Shan and Chen, 1993). From the Boltzmann equation, the LBM borrowed the idea of streaming and collision steps for a set of particle distribution functions (the fraction of the particles having velocities in the interval 𝒆𝒊 and 𝒆𝒊 + 𝑑𝒆𝒊, in the i- direction). The resulting equations are linear, explicit, and they recover the incompressible Navier-Stokes with second order accuracy. Due to its local nature, the LBM can easily be parallelised (Chen et al., 2014). The main governing equations in LBM are listed in Table 1. Table1. Main equations used in the study. Equations Descriptions 𝑓𝑖 (x + 𝑒𝑖 ∆𝑡, 𝑡 + ∆𝑡) = 𝑓𝑖 (x, 𝑡) + Ω𝑖 (𝑓(x, 𝑡), (𝑖 = 0, 1 … , M) (1) Evolution of distribution functions Ω𝑖 = 1 𝜏 (𝑓𝑖 𝑒𝑞 − 𝑓𝑖 ) (2) BGK Collision operator 𝑣 = 𝑐𝑠 2(𝜏 − 0.5∆𝑡) (3) Kinematics viscosities 𝜌 = ∑ 𝑓𝑖𝑖 𝜌u = ∑ 𝑓𝑖𝑖 𝑒𝑖 (4) Macro-density and macro-momentum equations 𝐹𝑖𝑛𝑡(𝑥) = − (𝑥) ∑ 𝐺(𝑥, �́��́� ) (�́�) (�́� − 𝑥) (5) Interaction force between particles at locations x and �́� 𝐺(𝑥, �́� ) = ∫ 𝑔, |𝑥 − �́�| ≤ 𝑐, 0, |𝑥 − �́�| > 𝑐, (6) Simplified Green’s function with single component system  = 𝜌0(1 − exp (−𝜌/𝜌0) (7) Effective mass 𝑃 = 𝑐𝑠 2𝜌 + 𝑐𝑠 2 𝑔 [ (𝜌)] 2 (8) Shan-Chen Equation of State 𝑃𝑖𝑛 − 𝑃0𝑢𝑡 = ∆𝑃 = 𝜎 𝑅 (9) Young-Laplace equation 𝑊𝑒 = 𝜌𝑖 𝑈𝑐 2𝑑/𝜎 (10) Weber number 𝑅𝑒 = 𝜌𝑖 𝑈𝑐 𝑑/𝜇𝑖 (11) Reynolds number based on droplet diameter Equation 1 describes, as mentioned above, streaming and collision steps. During collisional step, particles are allowed to exchange momentum but the sum of mass and momentum over each node must be conserved. All parameters in this study are given in lattice units (Lu). The 𝒇𝒊 in eq. (1) is the particle velocity distribution, while ∆𝒕 is the time increment, and 𝒆𝒊 is the microscopic velocity. 𝛀𝒊 is the Bhatnagar-Gross-Krook (BGK) collision operator which relaxes the generic 𝒇𝒊 towards its equilibrium value 𝒇𝒊 𝒆𝒒 which is the lattice equivalent of the Maxwell-Boltzmann distribution function. The relaxation time 𝝉 is the relaxation time that is related to 1268 kinematics viscosity through the eq (3), and 𝒄𝒔 is the lattice sound speed. Equation 4 shows that density and momentum of the fluid can easily be found through distribution functions. Interaction force between neighboring particles maintains phase separation and is expressed as discrete gradient of a local effective mass  , weighted by a Green’s function 𝑮(𝒙, �́� ). In the case of single component system, the latter parameter can be simplified to a single number, 𝒈. This number defines the interaction strength between the nearest neighbour particles, which can be repulsive or attractive. 𝒄 and 𝝆𝟎 are the lattice spacing and a normalization constant, respectively. 𝑼𝒄 is the channel centerline velocity, while 𝑷𝒊𝒏, 𝑷𝟎𝒖𝒕 are pressure inside and outside bubble/droplet, and 𝑹 is the radius of curvature. Finally, 𝒅, 𝝈, 𝝆𝒊, 𝝁𝒊 are, respectively, the diameter of the droplet/bubble, surface tension, density and viscosity of the i-th fluid. 3. Results and Discussion 3.1 2D Static Droplet If there is no interaction between particles, the Ideal Gas Equation of State (EOS) can be used. However, for two-phase systems, the ideal Gas EOS is no longer applicable and therefore, more realistic EOS is required to correctly describe the fluid behaviour under a given condition. In this study, Shan-Chen EOS (eq. 8) is used. By incorporating eqs. (7) and (8), non-monotonic relationship between pressure and density can be obtained. The parameters at the critical point can be derived by setting 1st and 2nd order pressure derivatives with respect to density (𝜕𝑃 𝜕𝜌 , 𝜕2𝑃 𝜕𝜌2 ) to be equal to zero. Hence, critical values are 𝜌𝑐 = 𝑙𝑛2 and value of 𝑔 at critical point (𝑔𝑐 ) is -4.0. Additionally, 𝑔 is related to temperature through 𝑇 = −1/𝑔. When 𝑔𝑐 < 𝑔, there is no phase separation, whereas when 𝑔𝑐 > 𝑔, the system is separated into two phases, thus two densities (for example, 𝜌𝑜 corresponds to gas density, 𝜌𝑖 is liquid density) can coexist at the same pressure. The surface tension is determined using Young-Laplace eq. (9) by measuring the pressure difference between the liquid and gas phases for a series of circular-shaped static droplets with different radii and then by finding the slope of the “∆𝑃 vs. 1/𝑅 ” relation. Figure 1 shows the expected linear dependence between the pressure difference and the inverse of the droplet radius. For the current case study, the surface tension is found to be equal to 0.1119 as depicted in Figure 1. Figure 1: Validation of the LBM results with the Laplace law. The slope of the graph is found to be 0.1119. Parameters are in lattice units (Lu) 3.2 2D Droplet Motion in Poiseuille Flow Next, the capability of the SC LBM to model the droplet flow at different Re and We numbers is tested. The results are compared against the data from Mortazavi and Tryggvason (2000), where the finite difference/front tracking method was used. Convergence criteria used in this study is based on the relative error of the velocity to be less than 1E-5. Further, grid dependency test was carried out with three grid resolutions of 64*64, 128*128, and 256*256. Figure 2 represents normalized fluid velocity along the channel length which is found using three above mentioned grid resolutions. Simulation results are taken when the fluid flow reached steady 0,493 0,494 0,495 0,496 0,497 0,498 0,499 0,500 0,028 0,030 0,034 0,040 0,044 0,050 0,059 0,064 P in -P ou t ( Lu ) 1/R (Lu) LBM simulation results Linear fit 1269 state at 199000 time-step and the position of the droplet normalized with respect to the channel height is 0.73. The points and relative position of the droplet at which the results are taken is also shown in Figure 3. The 2- Norm difference between the results of two last grid resolutions is 0.9881%. This value is considered to be low enough and, for that reason, the grid resolution of 128*128 is used in this study. Figure 2: Normalized fluid velocity along the channel length with three different grid sizes Figure 3: Points and position of the droplet in the channel at which the values of fluid velocity are taken Figure 4a shows the geometry of the domain and initial droplet position within it. Figure 4b presents the droplet position when the flow is fully developed. The droplet was initially placed close to the upper wall (Figure 4a) and its position in the channel changed slightly when pseudo-steady flow was achieved. The white line highlights this change in the final droplet vertical position when compared to the initial starting one. Periodic boundary conditions are applied at left and right side of the domain and no-slip boundary condition is implemented at the top and the bottom walls. The grid resolution is 128*128 Lu, the density ratio is 8 and the kinematic viscosity ratio is 1, the value of 𝑔 is -4.7, radius of the droplet is 16 Lu. The relevant non- dimensional numbers, the Weber and Reynolds (We and Re) numbers are 0.5 and 10, respectively. The droplet is moving together with the gas streams under the effect of the pressure gradient applied to the system. a) Initial condition b) Pseudo-steady state condition Figure 4: Droplet flow in 2D mirochannel 1270 Figure 5 shows the change in the droplet axial velocity as it passes through the entire channel length 40 times. As it can be seen from Figure 5, the droplet velocity initially increases considerably and remains steady further. Figure 5: Evolution of the droplet axial velocity at We=0.5, Red =10 Figure 6: Evolution of the droplet lateral position at We=0.5, Red =10 Figure 6 shows the change in vertical position of the droplet as it progresses towards the pseudo-steady sate condition in the channel. The droplet is initially placed close to the upper wall in both simulations, but it slightly bounces back monotonically from the wall in both cases until it reaches a constant vertical displacement. Figure 5 and 6 show an initial transient phase (about 10 times the channel length) where the behaviour is different from what was reported by the benchmark paper. That is due to the diffuse-interface nature of LBM. Initially the droplet undergoes a relevant deformation in order to accommodate the pressure distribution on the inside corresponding to the given density ratio. Once a steady droplet shape is reached and sometime elapses, the LBM correctly reproduces the droplet behaviour both in terms of vertical displacement and relative velocity with differences, respectively, in the order of 2.5% and 1% compared to the results from the benchmark paper. Compared to the front tracking method, the LBM is easier to set up and it does not require an explicit tracking of the phase interface. However, the LBM suffers from limitations in terms of density and viscosity ratios 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0 10 20 30 40 50 u /U c x/H Mortazavi and Tryggvason, 2000 LBM Results 0,5 0,6 0,7 0,8 0,9 1 0 10 20 30 40 50 z/ H x/H LBM Results Mortazavi and Tryggvason, 2000 1271 achievable. For example, the model exhibits numerical instabilities in the cases with high Re (~40) and We (> 2). Studies like Chen et al. (2014), Liu et al. (2010) or Yuan and Schaefer (2006) showed improvements in terms of maximal density ratio by using other EOS, instead of that given in eq. (8). In order to reduce the huge spurious currents arising at high density ratios, one needs to consider higher-accuracy schemes to evaluate the gradient in eq. (5), as well as alternative methods of incorporating the body force defined by eq. (5). In this study, we used the original Shan-Chen force scheme, however, there are other force schemes such as Guo’s (Hu et al., 2014) and the Exact Difference Methods (EDM) (Kupershtokh et al., 2009). These force schemes could also be tried as further improvements since they allow stable simulations at low kinematic viscosities. For example, a study of binary collisions of equally-sized liquid droplets with high density ratio (150) was performed using the SC-LBM coupled with the CS-EOS and the EDM force scheme (Monaco et al., 2014). In that study the authors also applied the Multiple-Relaxation-Time (MRT) collision operator instead of the single relaxation time BGK. The model has been proven to be effective and reliable by reproducing collision dynamics correctly. 4. Conclusions In this study, the multiphase Shan-Chen (SC) Lattice Boltzmann Method (LBM) was used to model the vertical displacement of a 2D droplet in the channel with resolution of 128*128 in comparison to a well-established method such as the Front-Tracking scheme adopted by Mortazavi and Tryggvason (2000). Except for a difference in the initial steps of the simulations, which is due to the different nature of two methods being compared, SC LBM correctly reproduces the final state of the droplet. However, the whole range of (We, Re) cases reported by Mortazavi and Tryggvason (2000) could not be reproduced due to the intrinsic limitations of the standard SC LBM adopted in this study. Future work will include three-dimensional cases, and a significant extension of the Weber and Reynolds number regimes by incorporating more effective numerical gradient schemes and new EOS together with forcing approaches, as previously discussed. Acknowledgments Authors want to acknowledge the School of Engineering at Nazarbayev University for supporting this work through the PhD scholarship to the main author Bagdagul Dauyeshova. Reference Chen, L., Q. J. Kang, Y. T. Mu, Y. L. He, and W. Q. Tao, 2014, A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications: International Journal of Heat and Mass Transfer, v. 76, p. 210-236. Chin J., Boek E.S., Coverney P.V., 2002, Lattice Boltzmann simulation of the flow of binary immiscible fluids with different viscosities using the Shan-Chen microscopic interaction model, Philos. Trans. R. Soc. A360, p. 547–58. Hu A., Longjian L., Uddin R., 2014, Force Method in a Pseudo-potential Lattice Boltzmann Model, Computational Physics, pp. 1-31. Kang Q., Zhang D., Chen S., 2002, Displacement of a two-dimensional immiscible droplet in a channel, Phys. Fluids 14, p. 3203–14. Kupershtokh A., Karpov D., and Medvedev D., 2009, On equations of state in a lattice Boltzmann method, Computers and Mathematics with Applications, 58(5), p. 965-974. Liu M., Yu Zh., Wang T., 2010, A modified pseudopotential for a lattice Boltzmann simulation of bubbly flow, Chemical Engineering Science, v. 65, pp. 5615-5623. Monaco E., Luo K. H., and Brenner G.,2014, Multiple-Relaxation-Time Lattice Boltzmann simulation of binary droplet collision, Microfluidics and Nanofluidics, 16 (1),. DOI: 10.1007/s10404-013-1202-0. Mortazavi S., Tryggvason G., 2000, A numerical study of the motion of droplets in Poiseuille flow. Part 1. Lateral migration of one droplet, Journal of Fluid Mechanics 411, 325-350. Shan, X. W., and H. D. Chen, 1993, Lattice Boltzmann Model for Simulating Flows with Multiple Phases and Components: Physical Review E, v. 47, p. 1815-1819. Yuan P., Schaefer L., 2006, Equations of state in a lattice boltzmann model, Phys. Fluids 18 (4), 042101- 042111. Zhang R., He X., Doolen G.D., Chen S., 2001, Surface tension effects on two- dimensional two- phase Kelvin- Helmholtz instabilities, Adv. Water Resour., v. 24, p.461-78. 1272