Lattice Boltzmann Methods for Industrial Applications - American

Subscriber access provided by Imperial College London | Library

Review

Lattice Boltzmann Methods for Industrial Applications Keerti Vardhan Sharma, Robert Straka, and Frederico Wanderley Tavares Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.9b02008 • Publication Date (Web): 02 Aug 2019 Downloaded from pubs.acs.org on August 5, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Lattice Boltzmann Methods for Industrial Applications Keerti Vardhan Sharma,∗,†,‡ Robert Straka,∗,¶ and Frederico Wanderley Tavares∗,†,‡ †Escola de Qui´ımica, Federal University of Rio de Janeiro, CEP:21949-900, Rio de Janeiro, Brazil ‡PEQ/COPPE, Federal University of Rio de Janeiro, CEP: 24210-240, Rio de Janeiro, Brazil ¶Department of Heat Engineering and Environment Protection, Faculty of Metals Engineering and Industrial Computer Science, AGH University of Science and Technology, Al. Mickiewicza 30, 30-059, Krakow, Poland E-mail: [email protected]; [email protected]; [email protected] Abstract We present an overview of recent advancements and potentialities for industrial applications of Lattice Boltzmann Methods (LBM). Various processes and phenomena with industrial importance such as flow in microchannels, dissolution and precipitation of solids, bubble dynamics, droplet formation and breakage, chemical reactions and reactors, mass transfer and optimization in fuel cells, flow in stirred tanks, multicomponent-multiphase flow, flow in porous media, and natural and convective heat transfer have been covered. The LBM implementation for such applications is presented and the results are analyzed here. Various LBM collision schemes and their impact on the numerical stability and accuracy in the above-mentioned processes are also discussed.

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction Lattice Boltzmann Methods (LBM) are popular and credible numerical tools to simulate a large number of physical and chemical phenomena associated with the flow of fluid and heat transfer. 1–5 It originates from the Lattice Gas Automata (LGA), which suffered from numerical instability and inaccuracy due to the Boolean nature and highly constrained simulation domain. 2,6–8 Lattice Boltzmann Equations (LBE) were then proposed, which considered the fluid as a collection of particles with continuous velocity distribution functions followed by other improvements, to remove the drawbacks associated with the LGA. 8–10 The generalized form of the LBE can be obtained by discretizing the Boltzmann Transport Equation (BTE). The BTE is a non-linear partial integro-differential equation which governs the advection of fluid and the collisions between those fluid particles by obeying the laws of conservation. The LBE is then adapted to a lattice model which fulfills the conditions to recovering the correct form of e.g. the Navier-Stokes equations, when the LBE is expanded using Chapman-Enskog (C-E) expansion or techniques such as Equivalent Partial Differential Equations (EPDE). 6,10–15 The mathematical formulation of LBM consists of two important steps, (1) streaming or propagation of fluid particles and (2) collisions among the particles in the fluid. It is also a common practice to separate the LBE into two operators, (1) streaming operator and (2) collision operator. The collision operator conserves mass, momentum, and energy. The collision operator in its original integral form in the BTE is highly non-linear. BhatnagarGross and Krook (BGK) linearized the non-linear collision operator, by representing it by simple relaxation process, which is known as the BGK collision scheme. 16 The BGK operator is also known as the single relaxation time collision scheme because it considers that the relaxation of all non-conserved quantities to their equilibrium is with identical relaxation time. Before getting into more details on collision schemes, it is necessary to mention that the evolution of fluid particles is governed by BTE and the system is described by the

2


Page 2 of 105

Page 3 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


velocity distribution functions, which are nothing but the probability of localization of such fluid particles at certain time and space interval. 1,17–19 Another important step is to choose a suitable lattice model. A lattice model is an arrangement of nodes in space which can also be defined as a cartesian grid with lattice nodes. 20 The fluid particles move from one node to another with microscopic velocities via links that connect these nodes, and the collisions occur at the nodes. The BGK collision scheme was originally defined in the velocity space. 16 When the LBE is written in discretized form for a lattice model, some restrictions on its performance are imposed as a result, but a lattice model must support the recovery of Navier-Stokes equations and maintaining its isotropic properties. Lattice models are limited velocities models, and this constraint results in the violation of the Galilean invariance. 21–23 The Galilean invariant terms can be recognized as the coefficients of the third order derivatives in the recovered hydrodynamic equations by applying the C-E expansion or EPDE techniques. Because all Lattice Boltzmann equations are discretized, therefore, it is expected to have the drawbacks associated with the limited velocity present in the numerical simulations. These drawbacks can be sorted out by increasing the number of velocities, e.g. using high order lattices or using entropic LBM approach, which we discuss later in details. The BGK collision scheme also suffers from numerical instability and inaccuracy due to a single relaxation parameter used in computations. The result is that all non-conserved hydrodynamic quantities relax to equilibrium with identical relaxation rate, which is not a very physically correct assumption. In fluid flow, different events occur at different time scale, therefore different hydrodynamic quantities must relax differently. There exist some non-conserved ghost quantities which do not appear in the second order truncated N-S equations, but they appear in the LBM scheme marring its stability. Moreover, one relaxation parameter limits the application of the LBM scheme for variable temperature problems in thermal flow because it fixes the Prandtl number. Therefore, to increase the stability of the

3



LBM numerical scheme, to get more control over the ghost moments, and to be able to solve thermal flow problems with a range of temperature variation, Multiple Relaxation Times collision scheme was developed. 22,24 In the MRT scheme, streaming is performed in velocity space and collision is performed in the moment space. Each non-conserved moment can be relaxed to equilibrium with different and desirable relaxation parameters. The MRT schemes result in variable Prandtl number which renders the MRT LBM scheme suitable for thermal flows with varying transport coefficients. Also, the ghost moments which were somehow uncontrolled in BGK LBM scheme can be controlled independently in MRT LBM. 19 The MRT LBM scheme possesses greater stability and accuracy compared to the BGK numerical scheme. It must be stated here as well that MRT schemes are not any better compared to the BGK when it comes to the Galilean invariance preservation. There exist another collision scheme called cascaded collision scheme, which was developed to further increase the stability of the LBM scheme. 25,26 Cascaded collision scheme performs collision in central moment space, which is nothing but a moment space shifted by the macroscopic velocity of the system. This helps to reduce the degree of violation of the Galilean invariance. Recently many groups have adopted the cascaded LBM scheme to solve complex flow problems. 27–30 Another family of LBM schemes is called Entropic LBM and Karlin-Bosch-Chikatamarala (KBC) LBM, which are relatively recent methods, which are robust and stable. 31–34 The Entropic LBM is currently being investigated for further improvements and have shown good performance, and reduce the degree of violation of Galilean invariance, as well. Mostly, entropic formulation of the LBM schemes have produced unconditional stability and complete Galilean invariance 35–37 . A recent trend in the literature show that the entropic implementations are gaining much popularity recently in addressing complex flow problems more prone to numerical inaccuracies and instabilities to achieve better stability and accuracy. LBM schemes are not only limited to solving fluid flow, but thermal problems can also be solved successfully using the described LBM collision schemes. The LBM frameworks to deal with thermal problems are the following; (1) Multi-Speed (M-S) method, (2) Hybrid 4


Page 4 of 105

Page 5 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


method, and (3) Double Distribution Function (DDF) method. 18,19 M-S method deals with only one distribution function set and one lattice model and considers energy as a conserved quantity obtained from the velocity distribution function. 38–45 Hybrid LBM uses two numerical techniques to solve thermal problems, where LBM is used to solve the flow field, and one of the conventional numerical methods such as finite element, finite volume, finite difference method can be used to solve the temperature or energy equation. 46–50 The DDF scheme solves two LBEs separately for flow field and energy (or temperature). Both LBEs can be solved on one lattice model, or two different lattice models can be used to solve flow field and temperature simultaneously. 51–59 The high numerical stability and accuracy of the LBM schemes made them useful for various complex athermal and thermal problems. 60–65 As a result, researchers quickly adopt LBM schemes to solve various complex phenomena and processes with industrial importance. 66–84 This paper not only describes the existing LBM models implementations in field scale applications but also explores the potential of the LBM approaches that can efficiently address various problems found in the industry. The advantages and disadvantages of LBM approaches for large-scale applications are highlighted in Tab. 1. It must be noted that there is no well-defined barrier between LBM application and LBM industrial application. Industrial scale problems are large scale-multiscale problems, which need enormous computational resources to handle such heavy simulations. Threrefore, computational cost of any method play an important role to qualify the method as an efficient industrial scale solver. In this review, we have also addressed this issue. Parallelization is one of the strongest advantages of the LBM solvers. And as a result, recent trends show that GP-GPU computing, parallel computing are being used extensively to solve large scale simulations. From the commercial aspect, LBM is also used in recent commercial solvers like Powerflow, XFlow, and ultraFluidX. The first two are general solvers, the latter is focused on the aerodynamics of cars and buildings. The XFlow was e.g. used to predict noise in aeroacoustic

5



Table 1: Main advantages and disadvantages of using LBM for industrial applications Advantages Highly stable and accurate in solving incompressible flows, easy extension to solve complex thermal flows, ease of implementation of complex boundary conditions, effective and large scale parallel implementations.

Disadvantages High memory demands, nonhydrodynamic moments values not well-defined, solving highly compressible flows and large temperature fluctuations is challenging.

test cases. 85 The advantages of LBM application in simulation of industrial, and especially large scale problems, compared to other numerical methods are as follows; the easy handling of complex geometry using e.g. immersed boundary methods (IBM), off-lattice collisions or body-fitted mesh refinement, and effective and highly scalable parallel implementations on both multi-CPU, and GPU and multi-GPU. Disadvantages connected to LBM could be summarized as; high demands on memory used in computations which can pose some problems for 3D problems on single GPU, the precise values for non-hydrodynamic moments relaxation times and their influence on stability and accuracy is not known and still under investigation, and highly compressible flows and flows with large thermal differences are still challenging tasks. Bubble dynamics is one of the most important problems with industrial importance. From microfluidics devices to oil and gas industry bubbles flow are present. Li et al. 86 recently implemented a 3D MRT LBM to study bubble behavior in nitrogen and glycerol/water mixtures system. Shu and Yang 87 presented the direct numerical simulation (DNS) of bubble dynamics using MRT LBM. Abdollahzadeh et al. 88 recently developed a BGK-LBM considering molecular collision to simulate homogeneous and heterogeneous chemical reactions. Benamour et al., 89 also very recently, combined volume penalization method with BGK-LBM to simulate fluid-structure interaction. Droplet formation, deposition, and breakage are some of the very important processes in chemistry research with huge industrial importance. Wang et al. 90 used BGK-LBM to study droplet formation in microchannels. A ternary BGK-LBM was developed by Fu et al. 91 to study mixing inside a double-emulsion droplet. Wu et al. 92 6


Page 6 of 105

Page 7 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


proposed a hybrid LBM method to simulate the impact of a falling droplet on a liquid film. Another important droplet application was discussed by Hazlett and Vaidya, 93 the authors studied contact angle hysteresis and wettability phenomena in convergent-divergent media. Although the authors do not specify the type of the collision model and information about the used LBM, the study is of good importance for various applications in water and oil and gas industry. Zhang et al. 94 recently implemented a pseudopotential MRT-LBM to investigate the deposition of multi-droplets and its line morphology. Bararnia et al. 95 used BGK-LBM to study the coalescence and breakup of falling multi-droplets. The increased number of LBM implementations to study bubble dynamics show the popularity LBM models are gaining for such complex problems. We can see from the recent trends that due to the simplicity of the BGK collision scheme, it continues to be the favorite choice of the engineers or researchers. Though MRT collision schemes increase the level of complex mathematical formulation slightly but they produce better stability and accuracy. Flow in porous media is one of the most active fields in petroleum and chemical engineering and chemistry research, of great industrial importance. LBM has been frequently used to study various lab-scale, pore-scale and large scale flows. Porous media or rocks are complex geometry arrangements of solid matrix and pore space, which impose complex boundary conditions to the fluid being transported through the complex geometry. Due to the local particulate and kinetic nature of the LBM, the complexity associated with the boundary conditions is naturally addressed by the model, using bounce-back boundary concepts. 96 This makes the LBM good choice for simulations of fluid and thermal flow through complex geometry. Yiotis et al. 97 studied the viscous coupling effects in immiscible multi-phase flow using BGK-LBM and the free energy LBM framework for multiphase flow. Zhao et al. 98 recently simulated pore-scale study of the gas flow in shale using the regular BGK-LBM. Using technology such as micro-computed tomography (m-CT), 3D reconstruction of pore-space can be done, which is then used as the flow domain for LBM simulations. 99–105 Fractured 7



porous media are of great importance when it comes to hydrocarbon production. Gunde et al. 106 presented the pore-scale interfacial dynamics and oil-water relative permeabilities of capillary driven flow in fractured porous media using a two-color framework for multicomponent LBM with BGK collision scheme. Free energy framework for multiphase flow LBM was implemented to study immiscible displacement efficiency using pore-scale approach. 107 Flow in porous media is a very general term, there might exist single-phase flow, multiphasemulticomponent flow, phase transition and many other flow characteristics through complex geometry. Therefore there is no standard approach to deal with such flows, but depending on the flow type the numerical tools are selected. Most importantly, correct implementation of the boundary conditions or interfaces for such flows is given more importance than LBM collision scheme, without undermining the effects of fluid properties such as viscosity, diffusivity. For example, for very high viscosity or very low viscosity flow in porous media, the BGK collision scheme is not an appropriate choice. For simple athermal flows, due to just one relaxation parameter, boundary conditions and flow coefficients depend on the same parameter which renders the scheme inefficient. Various other recent LBM applications in porous media are e.g. analysis of heterogeneity and permeability in complex rocks, 108 where the authors adopted digital rock physics (using micro-computed tomography, m-CT) to study carbonate rocks. The effect of particle size distribution and packing compression on fluid permeability can be found in the paper by Vidal et al. 109 This study is of particular interest to industrial application point of view, as the authors describe the massively parallel implementation of the LBM and the effect of morphometric properties on fluid transport in porous media. Another important work is the evaluation of permeability and non-darcy flow in macroporous structures such as vuggy carbonates which bear great industrial importance by Sukop et al. 110 The authors used digital rock physics to study various types of real-world rocks formations such as limestone from Florida and used LBM to predict flow behavior. This work can be considered as a reference to study fluid flow in highly heterogeneous porous media. For multiphase flow in 8


Page 8 of 105

Page 9 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


porous media, liquid-vapor interface configurations in porous media were studied by Sukop et al. 111 Tian and Wang presented studies related to reactive transport in porous media 82 for the CO2 gas. The authors studied reactive transport in fractured porous media. The readers interested in carbon capturing in fractured porous media using LBM must refer to this paper as it presents a great insight into the geological sequestration of the CO2 . Moving to LBM applications in energy storage devices, for example, Li-ion battery, the electron transport in Li-ion battery during the discharge process was studied by Jiang and Qu. 112 More interesting examples where LBM models were implemented to study porous media flows with importance to oil and gas industry are multi-seam coalbed methane production simulations, 84 particle filtration process inside porous walls, 113 bubble generation in porous media (very important for enhanced oil recovery (EOR) process), 114 and flow through membranes. 115 The magnetohydrodynamic nanofluid transport in porous media was recently developed by Sheikholeslami et al., 116 which can be an interesting problem for further LBM implementation. LBM implementation to study fluid flow and heat transfer in packed beds has been studied by various groups to address various issues regarding transport in packed beds. 70,117–123 These examples show some complex porous media flow cases, where sometimes regular or uniform grids are inefficient in handling the complexity of the flow or the geometry. To provide a solution to this problem, the pioneering work by He et al. 124 comes very handy. The authors developed an interpolation-supplemented LBM to perform efficient and accurate simulations of the curvilinear coordinate system, and turbulent isothermal and thermal flows. It is important to know that the use of non-uniform mesh grids open up paths for much larger scale applications where simulations in the curvilinear coordinate system are of importance. Energy conversion is another very important industrial application of chemistry and chemical engineering. More specifically, the conversion of chemical energy to electrical energy takes place in electrochemical devices. LBM has also found ground to simulate chemical

9



reactions and heat transfer present in such devices due to reacting fuels. Fuel cells are one example of such electrochemical devices. The flow of fluid and heat transfer through a solid oxide fuel cell is a common process and can be simulated by LBM. 125 LBM has also been used to address various issues regarding the determination of various properties of fuel cell gas diffusion layer. 126 Mass transfer in fibrous media for flow battery electrodes was simulated by Salah et al. 75 who applied LBM to simulate water drainage in a gas channel to optimize polymer electrolyte membrane fuel cell optimization. Meso-scale modeling of coal devolatilization was performed by using BGK-LBM. 127 On large scale applications, the following are some examples where LBM was used to simulate flow characteristics on larger domains, to address applications with industrial importance. Recently Shu and Yang computed flow in stirred tanks using LBM and Large Eddy Simulation (LES) techniques. 128 The LES is the numerical method used to address turbulent flow. Zhao et al. used Shan-Chen type multi-component LBM with BGK collision scheme to track the displacement interface present while cementation of the horizontal wells. 81 Also, LBM has been frequently used in simulating fluid flow and transfer process in microfluidic devices. 129–131 Mino et al. simulated the flow of oil-in-water emulsion through a coalescing filter using free energy framework of LBM for multiphase flow. 132 Buick and Cosgrove applied BGK-LBM to compute flow field in a single-screw extruder. 133 Recently, a two-phase flow LBM 134 was implemented to study gas-solid flow and soot loading in diesel particulate filters, where LBM was used to simulate flow field and the cell-automation probabilistic model was used to study the migration of solid particles. 135 Melting of thermal storage media such as icy slurry used in cold storage systems was simulated using two-fluid BGK-LBM. 136,137 Laminar solid-fluid flow in microchannels for agglomerate breakage was simulated with the LBM. 138 LBM applications in the area of nanotechnology are growing recently and showing very good performance. Nanoparticles dispersed in various base fluids are used for Enhanced

10


Page 10 of 105

Page 11 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Oil Recovery (EOR) purpose in petroleum industry. 139,140 Sun et al. recently presented a magnificent critical review of recent progress in applications of nanoparticles in EOR. 141 Hammou et al. developed an LBM scheme based on two-relaxation times (TRT) parameters to compute solute transport in unsaturated water flow 142 . Very interestingly, Afrooz et al. presented the experimental study of co-transport of gold nanospheres and carbon nanotubes in saturated porous media, 143 which is an ideal problem for numerical validation using multicomponent LBM. Song et al. coupled multiphase LBM and population balance equation (PBE) to simulate dynamics particle aggregation in flowing and heated colloidal suspensions. 144 Another very recent study regarding the deposition of monosized spherical particles in homogeneous isotropic turbulence using LBM was presented by Wang et al. 145 Heat and mass transfer phenomena are of paramount importance and occur more often in various industrial applications and chemical processes. The simplicity of LBM designed for fluid flow or isothermal problems facilitated its extension to deal with thermal problems. Various LBM frameworks for thermal flow have been discussed earlier. A wide range of problems which have been addressed by LBM are such as Rayleigh-Benard convection, 56,146,147 two-phase Rayleigh-Benard convection with a deformable interface, 148 Rayleigh-Taylor instability simulations, 149 gas-liquid interface mass transfer, 150 conjugate mass transfer, 5,32,151–155 impact of magnetic field on natural convection heat transfer, 156 heat transfer between porous walls and air in thermal energy storage, 157 evaporation of thin seawater film on horizontal tubes, 158 convection of thermosolutal fluid with Soret and Durfour effects, 159 gravity effect on phase change heat transfer, 160 natural and forced convection. 27,161–169 It is evident from the literature review presented above that LBMs have been extremely successful and credible numerical tool to solve a large range of problems covering many engineering streams. Most of the applications of LBM presented above are very recent, and it hints that LBMs are gradually gaining grounds from other conventional CFD techniques. The ease of parallelization of the LBM code e.g. for GPUs makes them very attractive and

11



Page 12 of 105

computationally effective. In the next sections, we present the mathematical formulation of the LBM and their implementations to address various applications in the industry.

Lattice Boltzmann Equation The advection and collision of fluid particles in a system can be represented by the wellknown Boltzmann Transport Equation (BTE)

∂t f (x, ξ, t) + ξ · ∇f (x, ξ, t) = Ω(f ),

(1)

where the term ξ · ∇f (x, ξ, t) is known as the advection term and Ω(f, f ) is the collision operator, f is the velocity distribution function, which is responsible for bringing the kinetic nature to the equation, and ξ is the velocity attained by the fluid particles. The discretization of above BTE in space and time leads towards the LBE, which reads

f (x + 4x, t + 4t) − f (x, t) = Ω(f ).

(2)

The above equation can be rearranged in a new form

f (x + ci 4 t, t + 4t) − f (x, t) = Ω(f ),

(3)

where ci is the discretized form of the microscopic velocity, 4t is the time step. The term f (x + ci 4 t, t + 4t) is known as the advection or the streaming step, and f (x, t) + Ω(f ) is defined as the collision step. Every LBM scheme consists of these two steps. The collision operator conserves the mass, momentum and energy. The macroscopic quantities such as

12


Page 13 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


density ρ, velocity u and energy E can be computed as follows ρ=

X

fi (x, t)

i

ρu =

X

ci fi (x, t)

(4)

i

ρE =

1X fi (x, t)(ci − u)2 2 i

Lattice Models Lattice models can be understood as the set of vectors which impose limitations on the way of fluid streams and where collisions occur. Lattice models can be defined in 1D, 2D, and 3D. Each lattice model consists of nodes which are orderly arranged in the space. Each node (white circles in Fig. 1) is connected with its nearest neighbor nodes through links (denoted as solid arrows). These links limit the directions and movement of the fluid particles, limiting the number of velocities considered in the LBE. The node that lies in the center of the unit cell (resting node) has zero velocity, and starting from this resting node particles can move to eight different directions with independent characteristic velocities. The total of nonzero and zero velocities (which is nine in the shown example) are known as nine velocity lattice model. Because the lattice model is defined for the 2D case, this lattice model is called a D2 Q9 lattice model. For this lattice model, characteristic velocity vector c~i = (cix , ciy )T can be defined for each lattice node such that 



1 1 1 0 0 −1 −1 −1 0 [~ c1 , . . . , c~9 ] =  , 0 1 0 −1 −1 −1 0 1 1 Each characteristic velocity is provided with a weight wi , which appears in the definition of the equilibrium distribution function for each node (provided later). The lattice model must obey conservation of mass and momentum for athermal flow and must preserve the rotational isotropy. 3,19 Some popular lattice models in 2D and 3D are D2 Q9 , D2 Q5 , D3 Q15 , 13



Page 14 of 105

D3 Q19 , D3 Q27 . Texts by Kruger et al. 19 and Wolf Gladrow 20 provide the properties and selection criteria of various lattice models in great details.

Figure 1: Structure of a nine-velocity lattice model in 2D. The characteristic velocity c~i is denoted for ith lattice link.

The collisions occur only at nodes, which are governed by the term Ω(f ). Actually, the collision term impose the relaxation of the particle towards the equilibrium in one time step, and collision is followed by streaming, in other words, the rate of change of velocity distribution function fi of ith velocity towards fieq is defined by the collision operator. fieq is the equilibrium distribution function, which is given by the Maxwell-Boltzmann distribution function

fieq

= ρwi

~u · ~ci (~u · ~ci )2 ~u · ~u − 2 1+ 2 + cs 2c4s 2cs

14


,

(5)

Page 15 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where wi is the weight factor and cs is the speed of sound in lattice units. The above expression comes from the Taylor expansion of the Maxwell-Boltzmann distribution truncated at second order in Mach number. This expression for the equilibrium distribution works well with various 2D and 3D lattice models. For the D2 Q9 lattice model, the weight factors wi are such that

[w1 , . . . , w9 ] =

4 , 9

1 , 36

1 , 9

1 , 36

1 , 9

1 , 36

1 , 9

1 , 36

1 9

.

The text books by Kruger et al. 19 and Gladrow 20 discuss these lattice models and their corresponding weight factors in details. It must be noted that for each lattice model, the total number of conserved and non-conserved quantities are equal to the number of velocities of the lattice model. For example, for D2 Q9 lattice model, mass, and momentum in x and y direction are the conserved quantities, therefore, six non-conserved quantities can be defined for the lattice model. These non-conserved quantities are often called non-conserved P n moments, denoted as mαβ , where mxm yn = i cm ix ciy fi . The most common non-conserved moments defined for a nine-velocity models are mxx , myy , mxy , mxxy , mxyy , mxxyy . While m00 , mx , my are the conserved moments. m00 =

X

fi = ρ

i

mx =

X

cix fi = ρux

(6)

i

my =

X

ciy fi = ρuy ,

i

Collision Operators, Stability and Accuracy As mentioned earlier, the BGK, MRT, and cascaded collision operator are the most used collision schemes in LBM simulations. The mass, momentum, and energy are conserved in the collision. Mathematically, the conservation laws for mass and momentum can be

15



Page 16 of 105

represented as follows X

Ω(fi , fi ) = 0,

i

X

ci Ω(fi , fi ) = 0.

i

The basic details regarding the streaming step and collision step of these collision schemes are provided below one by one. In particular BGK is the most popular collision scheme. 12,38,170–176 It is a single relaxation time collision scheme (SRT), where Ω(fi , fi ) = −

fi (x,t)−fieq (x,t) , τ

for

which the BGK-LBE can be written such that fi (x, t) − fieq (x, t) fi (x + ci 4 t, t + 4t) − fi (x, t) = − , τ

(7)

where τ is the relaxation time by which all non-conserved quantities approach equilibrium. For BGK, the streaming step is fi (x + ci 4 t, t + 4t) and the collision step is fi (x, t) − fi (x,t)−fi (x,t)eq . τ

It is of equal importance to show that the LBE is consistent with the fluid

flow behavior which is governed by the Navier-Stokes equations. The consistency of an LBM scheme is established by recovering of N-S equations and analyzing them. The LBE can be expanded using e.g. C-E expansion 177–179 or EPDE technique 15,59,180 to recover the N-S equations for consistency analysis purpose. Readers are suggested to refer to the text by Kruger et al. for details on discussions on recovery and its analysis. The recovered hydrodynamic equation for momentum in x direction from BGK-LBM 59,181 reads

1 = τ− ∂xx 2

∂t ρux + ∂x p + ρu2x + ∂y ρux uy 2 2 ρux + cs ∂yy ρux + 2cs ∂xy ρuy .

(8)

The standard N-S equation for the momentum in x direction 19 reads

∂t ρux + ∂y ρux uy

= −∂x p + ∂y η ∂y ux + ∂x uy .

16


(9)

Page 17 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Comparing the N-S equation recovered from LBM, Eq. (8) and Standard N-S equation, Eq. (9), some relations between transport coefficients and collision scheme can be established, which directly reflects on the stability and accuracy of the LBM scheme for that particular collision operator. The pressure can be described as

p = ρc2s ,

and viscosity is given by η η = ρc2s τ −

1 . 2

The stability criteria for the BGK-LBM can be determined from the above expressions. The necessary condition for the stability of the BGK-LBM should be that the viscosity must not be negative or zero. These two situations are non-physical and render the BGK scheme unstable. Therefore, the necessary condition for the stability reads

τ ≥ 1/2.

It must be noted that for turbulent flow modeling using BGK-LBM, without increasing grid resolution, the viscosity is set to be very low, which means τ must be assigned values pretty close to 1/2, which induces instability. Therefore, it is a common practice to increase the grid resolution, which comes at the cost of increased computational cost. On the other hand, for flow with very high viscosity, τ must be set very large which induces numerical inaccuracy. The BGK-LBM also produces inaccuracy for such cases, where boundary conditions have a role to play. The position of the wall location directly depends on the relaxation time, therefore the location of the wall in numerical simulations is different from the correct physical location of the wall, which in turn produces inaccurate results. One example of this problem is solving fluid flow in porous media using BGK-LBM. The MRT collision operator was proposed to increase the stability and accuracy of the 17



Page 18 of 105

LBM scheme, which was encountered in BGK-LBM. The streaming in MRT schemes is performed in velocity space, and collisions are performed in the moment space. 24,47,171,182,183 MRT LBM scheme increases the number of degrees of freedom compared to the BGK collision scheme, offering more control over stability and accuracy. The non-conserved quantities (moments) relax to their equilibrium with different relaxation times, increasing the degree of freedom. The LBE for MRT reads

fi (x + ci 4 t, t + 4t) = fi (x, t) − M−1 R(mi (x, t) − meq i (x, t)),

(10)

where M is the transformation matrix through which velocity space f is mapped to moment space m, such that m = Mf , R is the diagonal matrix with relaxation times. For example for a D2 Q9 lattice model, one must define three conserved moments (mass, and momentum conservation in x and y direction) and six non-conserved moments, as described earlier already.  0  0    0   0   R= 0   0   0   0   0

 0 0

0

0

0

0

0

0 0

0

0

0

0

0

0 0

0

0

0

0

0

0 0

1 τ1

0

0

0

0

0 0

0

1 τ2

0

0

0

0 0

0

0

1 τ3

0

0

0 0

0

0

0

1 τ4

0

0 0

0

0

0

0

1 τ5

0 0

0

0

0

0

0

0  0    0   0   0 .   0   0   0  

1 τ6

The conserved moments do not relax to the equilibrium as they are in the equilibrium at every time step X i

fi =

X

fieq = ρ,

i

18


Page 19 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


X

ci f i =

X

i

ci fieq = ρu,

i

From definition of moments provided earlier, the above expressions can be reformulated as follows m00 = meq 00 = ρ, mx = meq x = ρux , my = meq y = ρuy By adopting the similar C-E or EPDE procedure to recover the N-S equations, the accuracy and stability of the MRT-LBM can be analyzed.

1 ∂xx ρux ∂t ρux + + + ∂y ρux uy = τ1 − 2 1 1 1 ∂yy ρc2s ux + τ1 − + τ3 − ∂xy ρc2s uy . + τ3 − 2 2 2

∂x ρc2s

ρu2x

(11)

Comparing Eq. (11) with Eq. (9), we can determine the transport coefficients as before, the pressure p = ρc2s , the shear viscosity

η=

ρc2s

1 , τ3 − 2

the bulk viscosity νb can be written as

νb =

ρc2s

1 η τ1 − − . 2 3

The stability and accuracy analysis of the above presented MRT scheme can be done

19



by analyzing the above mentioned numerical viscosity terms. It can be seen that due to more parameters the degree of freedom of the MRT scheme is higher. In MRT, we have more control over the relaxation of the high degree moments such as mxxy , mxyy , mxxyy , which do not contribute to the physical behavior of the fluid flow as per N-S equations but contribute to the numerical stability. 19,184–187 It must be noted that the recovered N-S equations are truncated at O(u2 ) and the rest of the terms of order O(u3 ) and beyond are error terms, where appear the relaxation parameters for high order moments. By controlling these parameters the accuracy of the LBM scheme can be increased, which is impossible in the BGK case. Regarding stability, the issues generally occur at very low shear viscosity (turbulent flow). To maintain the MRT stable in turbulent flow conditions, for very small η the parameter τ3 has to be very small, therefore the bulk viscosity νb can be increased accordingly to keep the numerical simulations stable. The transformation from velocity space to moment space to relax hydrodynamic quantities directly, and relax them with different relaxation time parameters increased the stability and accuracy of the LBM numerical scheme. But it did not deal with the inaccuracies associated with the violation of the Galilean invariance. What causes this Galilean invariance violation is the constraint of the limited velocity sets used in the lattice models, which somehow hammer the isotropic properties to recover the correct form of N-S equations from the LBE. Therefore, to reduce the degree of violation of the Galilean invariance, Geier et al 26 shifted the moment of the system by its macroscopic velocity u, giving birth to the cascaded central moment LBM scheme. The other important fact to mention here is that the relaxation of the moments in the cascaded scheme is not independent like MRT-LBM. In the cascaded collision scheme, relaxation of high order moments is affected by the relaxation of lower-order moments, creating a cascade effect in relaxation towards high-order moments. The streaming for cascaded LBM is also performed in velocity space, and collisions are performed in central moment space. The definition of the central moment κ reads

20


Page 20 of 105

Page 21 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


κxm yn =

X

(cix − ux )m (ciy − uy )n fi

i

The lattice Boltzmann equation for the cascaded LBM reads

fi (x + ci 4 t, t + 4t) − fi (x, t) = K.~k,

(12)

where K is the collision matrix, ~k is a central moment vector which constructs the cascaded relaxation. This scheme is a semi-implicit scheme, which deviates from the linear schemes such as BGK and MRT. The details regarding cascaded LBM can be found in Ref. 25,26,59,188,189 The hydrodynamic equations from cascaded LBM scheme and its stability, accuracy analysis have been presented in these papers. The cascaded LBM is gradually becoming popular and used recently by various groups. 28–30,59,181 The results show that cascaded collision scheme increases the stability and accuracy beyond the BGK and MRT collision schemes.

Industrial Applications of LBM LBM for Bubble Dynamics Interfacial phenomena are one of the important phenomena which can be solved efficiently by LBM. 190,191 Problems with the multiphase flow are addressed by LBM using three frameworks, (1) color-gradient models, (2) Shan-Chen model, and (3) free energy models. 17,19 Particularly, The Shan-Chen model and the free energy model are the popular approaches to address multiphase-multicomponent flows in the LBM framework. Gunstensen et al. 192 published the pioneering work in which they presented the color-gradient model for immiscible fluids. Recently, Yu et al. 193 adopted the color-gradient model to simulate ternary fluid flows for a full range of interfacial tensions in critical states. The Shan Chen model is named after the authors X. Chan and H. Chen who wrote the most noted work in the 21



multiphase multicomponent application of the LBM 194 . The free energy model, as suggested by the name, contains free energy functional which makes the model thermodynamically consistent. The free energy model was developed by Swift et al. and can be found in the pioneering work by the authors 195 . In one recent important development, Semprebon et al. developed a ternary free-energy LBM where the surface tensions and contact angles are tunable independently 196 . A detailed implementation of this ternary free-energy model using the BGK collision scheme can be found in the paper, which can serve as a base to develop an efficient multiphase-multicomponent thermodynamically consistent model using different collision schemes such as MRT and more sophisticated collision operators. This scheme is suitable for density-match fluids cases and not for high density-viscosity ratios. Taking this work one step further, Wohrwag et al. 197 proposed the entropic implementation of the freeenergy LBM model for high-density ratios. This paper can be considered as an important development in the free-energy LBM model for multiphase flow. Bubbles are gas-liquid contact systems, which are immensely important in heat and mass transfer applications. Hydrodynamics of fluid flow in bubbly flows are highly influenced by bubble behavior and bubble collisions can induce bubble breakage or coalescence. The bubble-bubble interaction for large gas-liquid density ratio can be studied by using pseudopotential based single component multiphase LBM 191,198 . A force term is needed to address the interaction between gas and liquid. It must be noted that the introduction of a force term in LBE changes the equation of state from ideal to non-ideal. The force term to address the interaction force and external force can be given by

F (x) = Fgravity (x) + Finteraction (x),

where Fgravity is the gravitational force which acts as an external force on the bubbly flow.

22


Page 22 of 105

Page 23 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The mass and momentum conservation for such system reads 198

ρ=

X

fi

i

ρu =

X i

1 ci fi + F (x) 2

In LBM, fluid is considered as fictitious particles, and each particle is an ensemble of a large number of molecules, therefore to describe the interaction force, inter-molecular potential such as Lennard-Jones cannot be used. 191 Shan and Chen 194 suggested the nearestneighbor interaction potential to write down Finteraction . The interaction potential can be written such that V (x, x + 4x) = Gwi ψ(x)ψ(x + 4x), where G is the strength of the interparticle interaction, ψ is a particle density function and wi is the weight factor associated with the particular lattice link. The interparticle interaction force can now be written as follows

Finteraction (x) = −Gψ(x)

n X

wi ψ(x + 4x) 4 x,

(13)

i=1

where n is the number of nearest-neighbor nodes. The dimensionless parameters which we need to study bubbles are Morton number, Mo = and Reynolds number, Re =

ρU d . µl

g(ρl −ρg )µ4l , ρ2l σ 3

Eotvos number, Eo =

g(ρl −ρg )d2 σ

The surface tension is denoted by σ, ρl and ρg are liquid and

gas densities, respectively, µl is the viscosity of the liquid, d is the characteristic length, and g is gravitational constant. Liu et al. 191 adopted the BGK collision scheme to study bubble dynamics, while Yu et al. 199 adopted MRT collision scheme to simulate bubble dynamics. In this very important paper 199 , Yu et al. presented an excellent analysis of the work done by some groups on bubble dynamics using LBM 200,201 for a very small liquid-gas density ratio, and also highlighted the benefits of using improved collision schemes such as MRT for better accuracy and numerical stability. As it is a well-established fact that MRT collision scheme 23



provides higher stability and accuracy compared to the BGK collision scheme because the MRT scheme provides a higher degree of freedom, i.e. more parameters in terms of the relaxation times to control the transport coefficients efficiently. The authors stated that the impact of density ratio on bubble behavior was negligible for values larger than 50. 199,202 Therefore, it must be noted that for efficient simulations of bubble dynamics, the density ratio must be larger. Yu et al. 199 investigated bubble rise, bubble drag coefficient, and bubble interaction in the homogeneous regime. In another important development, Yu and Fan 203 increased the grid resolution by using a technique called adaptive mesh refinement (AMR) to provide better a definition of the interfaces. Therefore, for efficient simulations of bubble dynamics with Shan-Chen model with higher accuracy and stability, more sophisticated collision scheme such as MRT, cascaded or entropic LBM schemes can be incorporated with the techniques developed by Yu and Fan to better define the interfaces. We also must remember that the correct implementation of the multiphase effects through the forcing term near the boundary play a great role in the accuracy and stability of the LBM simulations of bubbles flow. The Shan-Chen type pseudopotential methods lack in thermodynamic consistency and local momentum conservation, therefore, it was suggested to use the phase-field or free energy model, which can take into account the thermodynamics efficiently. 19,87 Shu and Yang discussed in their paper 87 that BGK-LBM was the most frequently used scheme to solve bubble dynamics using the phase-field method framework of multiphase LBM. They also highlighted that most groups 204–206 used BGK-LBM and phase field methods for highdensity ratio cases, but for limited Re and Mo numbers. Since bubble dynamics can be considered a low viscosity flow scenario, therefore BGK does not remain the correct choice, and one should use MRT or much-advanced collision schemes such as cascaded or entropic schemes for better accuracy and stability. And also, these studies were limited for lower Re number and high Mo numbers, which fall short in addressing the real-world scenarios present in bubbly flows in the industry. It must be mentioned that Shan-Chen type methods 24


Page 24 of 105

Page 25 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and free energy models struggle in addressing the high viscosity ratio flow due to numerical instability. 190 All these shortcomings were addressed in an important paper by Shu and Yang who implemented the MRT collision scheme along with the free-energy multiphase model to address bubbles flow with large Re numbers and low Mo numbers (more closer to industrial problems) for a high-viscosity ratio. The authors provided simulations for a single bubble, bubble pair, and bubble swarm, for a very large kinematic viscosity ratio of 1 : 1000. Below we present in brief how Ghosh et al., 73 and Shu and Yang 87 combined the LBM with the phase-field method to solve bubbles dynamics and flows. The phase-field LBM models generally employ Cahn-Hilliard equation to capture the gas-liquid interfaces, and the system of equation governing such flow can be presented as a combination of Cahn-Hilliard equation and N-S equations. The Cahn-Hilliard equation reads

∂t φ + ∇ · (φu) = θm ∇2 µφ .

(14)

The mass and momentum conservation equations read ∂t ρ + ∇ · (ρu) = 0 (15) 2

∂t (ρu) + ∇ · (ρuu) = −∇ · p + µ∇ u + F, where θm is the mobility, µφ is chemical potential and φ the scalar quantity. To reduce the numerical oscillations, Lee and Lin, 205 and Zheng et al. 206 constructed a relation for pressure tensor gradient, −∇ · p = −φ∇µφ − ∇p0 , to combine the N-S equations with the Cahn-Hilliard potential. 87 The term p0 is nothing but the equation of state for an ideal gas, i.e. p0 = ρc2s , same as in LBM recovered N-S equations for an ideal gas (without force terms). The addition of the force term in LBM motivates the system to deviate from p0 to p. To solve the above system of equations, the chemical potential must be determined from the

25



Page 26 of 105

relation between free energy and interfacial tension and bulk free energy density 206 such that

µφ = ∂φ ,

where is the free energy density. Now the LBM implementation to solve above system of macroscopic equations was presented by Shu and Yang, where the authors simulated 2D and 3D flows using D2 Q9 and D3 Q15 lattice models. We have already discussed the techniques and recovered N-S equations from the LBM scheme. Here we have two systems of equations to solve, (1) Cahn-Hilliard equation (phase field method) and (2) N-S equations for mass and momentum conservation. The first step to solve this system of equation is to define the corresponding LBEs. The MRT-LBE for Cahn-Hilliard reads

fi (x + ci 4 t, t + 4t) = fi (x, t) − M−1 R(mi (x, t) − meq i (x, t)),

(16)

which is the general and valid form of LBE for any MRT case. Now the most important part is to define how the moments in presence of chemical potential (intefacial and free energy environment) relax towards equilibrium. Therefore, there is a need to define a new expression for f eq , because the previously defined fieq Eq. (5) is not consistent with the Cahn-Hilliard approach, with external and interaction forces between phases. For resting node (i=0), the equilibrium for Cahn-Hilliard reads

f0eq

(1 − w0 )θµφ = φ− . c2s

(17)

For other lattice nodes with non-zero velocities, the equilibrium rule for Cahn-Hilliard can be written as fieq

= wi

θµφ φci · u + c2s c2s

26


(18)

Page 27 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The conservation rule for Cahn-Hilliard equation reads

φ=

X

fi =

X

i

fieq

i

The above system of equations can now be solved on any lattice model. After adopting the C-E expansion or EPDE technique, we can recover the Cahn-Hilliard equation from the MRT-LBM for mobility coefficient θm and stability and accuracy checks. The mobility θm can be determined by comparing the original Cahn-Hilliard equation with the recovered macroscopic equation, which is θm = θ τ −

1 , 2

it should be kept in mind that the relaxation parameter τ is generalized here, and for an MRT scheme and depending on the lattice model we choose, the τ parameter can be written more specifically. After solving the Cahn-Hilliard equation, similar steps must be taken to solve the system of mass and momentum conservation (N-S equations) using MRT-LBM. For simplicity, we can choose a different notation for velocity distribution function to solve N-S equation, say g. It must be noted that the force term is added to the flow field LBE, and can be directly included in the collision step. The MRT-LBE 207 with force term reads gi (x + ci 4 t, t + 4t) = gi (x, t) − M−1 R(mi (x, t) − meq i (x, t)) 1 −1 (ci − u) · F eq gi + I− M R 2 ρc2s

(19)

The conservation rules stand similar to those we already mentioned earlier with forcing term. The N-S equations can now be recovered as usual to determine transport coefficients, and perform stability and accuracy analysis. It was shown that for D2 Q9 and D3 Q15 lattice

27



models, the kinematic viscosity is of the standard form

ν=

c2s

1 τ− . 2

Now the LBM implementation presented above can be used to address various bubble related flows such as static bubble immersed in liquid, the rise of a single, pair, and swarm of bubbles, etc. The effect of discrete phase densities on bubble dynamics for air-water, liquefied LPGwater, and kerosene-water systems are presented below, where the water is the continuous phase. The densities of the discrete phases air 1.3 kg/m3 , liquefied LPG 560 kg/m3 and kerosene 787 kg/m3 . The results of bubble dynamics using BGK-LBM are presented in Fig. (2). Due to higher density difference, the flow field for air bubble in water is faster than the velocity of LPG and kerosene in water. Therefore, for a large density difference system, the lower discrete density moves faster.

Figure 2: Bubble shapes and positions for (a) air-water, (b) liquefied LPG-water, (3) kerosene-water. The water is kept as the continuous phase. Reprinted with permission from Ghosh et al. 73 Copyright (2012) Americal Chemical Society.

The effect of surface tension on bubble dynamics is presented in Fig. (3). The bubble 28


Page 28 of 105

Page 29 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


shape is directly influenced by the surface tension once the bubble of the discrete phase is placed in a column with a continuous phase. Ghosh et al. 73 simulated this effect for kerosene being the discrete phase in water, the continuous phase for three different surface tensions, 0.96, 0.32 and 0.045 N/m. It is evident from Fig. (3) that lower surface tension induces faster movement and stretching of the bubble.

Figure 3: Effect of surface tension on discrete phase kerosene cylindrical bubble in continuous phase water for three different surface tensions, (a) 0.96 N/m, (b) 0.32 N/m and (c) 0.045 N/m after 0.48 s. Reprinted with permission from Ghosh et al. 73 Copyright (2012) Americal Chemical Society. So far we have seen that MRT LBM scheme was incorporated to improve the performance of the numerical schemes used in Shan-Chen and free energy models. The use of more sophisticated collision schemes such as MRT, cascaded, cumulant 208 , and entropic collision schemes open up more possibilities for future research to improve performance. For example, Mazloomi et al. developed an entropic LBM for multiphase flow 209 . Entropic LBM approaches are well known for their unconditional stability and accuracy for complex flows 31 . The multiphase flow solver based on the cascaded LBM was proposed by Lycett-Brown et 29



al. 210 , in which the authors studied the phenomena of binary droplet collision.

LBM for Droplets Applications The interactions between droplets and surfaces are of great importance in microfluidic devices and interfacial phenomena at micro-level. Droplet falling on thin films, solid surfaces, a droplet in the droplet, interaction between droplets and particles, droplet formation in microchannels, droplet coalescence and breakage, mixing inside micro-droplets are some of the events that can occur in various chemical processes. Microfluidic systems have a greater role to play in various streams such as biological, chemical and physical research. 211,212 Droplet formation is the process which is abundantly present in a large number of microfluidic devices because, before droplets coalescence and breakages, the droplets must be formed. Lab-ona-chip devices are the best examples which deal with monodispersed microdroplets. 213 The LBM frameworks of multiphase multicomponent models are also applicable to deal with droplets. Pilch and Edman 214 described six manners based on Bo numbers in which a droplet can break. The droplet dynamics of falling drops or coalescence can also be solved by using similar system of equations, LBM to address mass and momentum conservation equations, and Cahn-Hilliard equation with the phase field method to address the liquid-solid interface efficiently.

30


Page 30 of 105

Page 31 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4: Effect of density ratio and Bo number on falling drop impact on liquid film for (a) ρ∗ = 5 and Bo = 20, and (b) ρ∗ = 100 and Bo = 100 (Wu et al.) 92 Reprinted with permission. Copyright Elsevier 2015.

The falling of droplets and its impact on a liquid film can be solved by using hybrid LBM or pure LBM approach. The N-S equations can be solved by LBM and Cahn-Hilliard equation can be solved by any of the FDM, FVM or FEM methods, or vice versa. Though it is highly recommended to use LBM approach for both the equation systems, because of the benefit of parallelization of the codes and to avoid issues due to coupling between two different numerical approaches. Let us discuss both cases to simulate falling drop phenomenon. Before 31



to solve the problem of a drop falling freely under the influence of gravity using the numerical techniques, we must define some non-dimensional parameters for droplets. The concept of density ratio and viscosity ratio remains the same as for the bubble dynamics. Two new parameters related to the geometrical properties of drops which describe droplets behaviors are Bond number (Bo) and Ohnesorge number (Oh). 92 The expressions for Bo and Oh are as follows

ρ∗ =

ρl , ρs

µ∗ =

µl , µs

4R2 ρl Bo = σ, g

Oh = √

µl , 2Rρl σ

where ρl and ρs are densities of the liquid droplet and ambient fluid, µl and µs are dynamic viscosities of droplet fluid and ambient fluid, σ is the surface tension, and R is the droplet diameter. Wu et al. used a hybrid approach to solve the drop falling and its impact on the liquid film. The BGK-LBM was used to address the flow field and similar mathematical formulation as presented in bubble dynamics case has been used. One new detail comes from the fact that the drop is falling freely under gravity, therefore an extra force exerted by the gravity on the drop, which is Fg = g(ρs −ρl ) is added to the force term F . The Cahn-Hilliard equation was solved by a semi-implicit second-order scheme for temporal discretization and a fifth-order weighted essentially non-oscillation (WENO) scheme for space discretization. 92,215 The effect of density ratio (ρ∗ ) on the shape of a falling drop is presented in Fig. (4). Two density ratio values for two different Bond numbers were chosen, and instantaneous flow

32


Page 32 of 105

Page 33 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


patterns are presented. The results show that for lower density ratio, only coalescence is observed, but for larger density ratio and large Bond number, the shape of the drop remains much intact and after impact, coalescence and splash phenomena are observed. In a pure LBM approach, both the N-S equations and the Cahn-Hilliard equations (phase-field LBM approach) can be solved by LBM collision schemes. Another way to define the interfaces is to adopt the pseudopotential methods discussed earlier, i.e. Shan-Chen model. Now imagine a different scenario when the drop falls on a solid surface, i.e. drop impingement. Due to the presence of a solid surface, the wetting behavior, i.e. contact angle must be defined. The three-phase contact angle Θ is given by 216 −n · ∇ρ π −Θ = , tan 2 |∇ρ − (n · ∇ρ)n| where n is a normal unit vector. The above expression for the contact angle is discretized and the desired wetting angle is achieved. For rough surfaces, Gac and Gradon used LBM to model deposition of small droplets 217 and droplet-particles interaction. 218 Fu et al. used multiphase multicomponent LBM to model Villermaux/Dushman competing parallel reactions 219 and reactive mixing behavior. Yang et al. simulated the droplet formation and cell encapsulation process in a microfluidic device with a T-junction shown in Fig. (5). 220 The authors combined the BGK-LBM for Shan-Chen type model with the immersed boundary method to solve the flow field with encapsulated cells. The droplet break-up regime is known as squeezing regime (see Fig. 6) for the lower capillary region, the dripping regime (see Fig. 7) for the medium capillary region, and the jetting regime (not shown) for the higher capillary regions. 220 In another study, Zhao et al. simulated the mixing of droplets inside micro-channels using pseudopotential LBM method. 221

33



Figure 5: The schematic diagram of the microfluidic T-junction device (Yang et al.) 220 Reprinted with permission. Copyright Elsevier 2012.

Figure 6: Droplet formation in the squeezing regime (Yang et al.) 220 Reprinted with permission. Copyright Elsevier 2012.

34


Page 34 of 105

Page 35 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7: Droplet formation in the dripping regime (Yang et al.) 220 Reprinted with permission. Copyright Elsevier 2012.

Recently, Abadi et al. 222 proposed a pure LBM approach for three-component multiphase flow for high density and viscosity ratios. The interfaces in the three-component system were characterized by phase-field modeling using the Cahn-Hilliard theory. The scheme was free from numerical oscillations and produces high accuracy. The authors adopted the work by Boyer and Lapuerta 223 and Boyer et al. 224 regarding Cahn-Hilliard theory to define the threecomponent free energy and the bulk free energy for a fluid. In an important work, the effects of the incompressibility on parasitic currents elimination in two-phase flow were studied by Lee 225 . Another deserving mention is the multicomponent LBM solver proposed by Liang et al. 226 . This method computes fluid velocity independently of the forcing term, which makes the method suitable for multiscale analysis. Recently, Ridl and Wagner 227 proposed an LBM model for a mixture of van-der-Waals fluids using mean-field forcing originated from the gradient of the chemical potential and not from the pseudo-potential approach like used in the Shan-Chen model. Very interestingly, the authors defined the lattice-free energy of a mixture of such fluids and efficiently recovered the theoretical phase-diagrams 35



for a mixture of two and more van-der Waals fluids. The flow of ternary fluid system with a liquid drop through a regularly textured/patterned surface was developed by Sadullah et al. 228 using the free-energy LBM. The authors studied the partial wetting lubricants ridge to study their behavior such as contact line pinning and viscous dissipation for different angles. The readers interested in modeling of droplets’ contact angle behavior on patterned surfaces are recommended to refer to an important paper by Kusumaatmaja and Yeomans 229 and Varagnolo et al. 230 which discuss contact angle hysteresis on patterned and superhydrophobic surfaces, and sliding of water drops on a chemically heterogeneous surface, respectively.

LBM for Reactors and Reactions Reactive flow is another important research area with huge industrial importance. Bioengineering, fuel cells, combustion, high vacuum heterogeneous catalysis, bioreactors, pharmaceuticals, scalar mixing, fermentation, and catalytic reactors are some examples of reactive flows and reactions. There are four types of numerical models to address chemical reactions, which are (1) macroscopic, (2) mesoscopic, (3) microscopic and (4) multiscale model. 88 The atom-atom or molecule-molecule interaction come under microscopic approach (molecular dynamics), which is immensely computationally costly for hydrodynamics. The macroscopic numerical approaches are techniques such as FDM, FVM, FEM, which are well known to solve fluid flow. These techniques use and solve N-S equations. But these methods are not very efficient when it comes to solving reactive flows, because in reactive flows there are complex boundary conditions which these model cannot take into account to address the complexity of fluid-solid interfaces, and these models are continuum approach which overlooks the microscopic features of the fluid-solid interfaces. 88 We have already discussed that LBM is a kinetic approach and can solve flow at the micro/meso level and can be extended to describe the hydrodynamic description of the flow system. Due to its kinetic nature, it can naturally address the complex boundary conditions and interfacial phenomena. Solving chemical reactions using LBM comes under double distribution function approach 36


Page 36 of 105

Page 37 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(used for thermal and a scalar transport; advection-diffusion equation). One LBE is written to address mass and momentum conservation, and another LBE written to address the scalar transport such as concentration of each component that can react. The chemical reactions can be modeled by defining a source term for products and sink term for reactants and adding them to the collision step. The mathematical formulation for chemical reactions is as follows. The LBE for flow field, mass and momentum conservation, viscosity terms remain the same as all previous cases

fi (x + ci 4 t, t + 4t) − fi (x, t) = −

fi (x, t) − fieq (x, t) , τf

where τf is the relaxation time for flow field. The LBE for concentration 88,231 (which is a scalar quantity obeying advection-diffusion equation) can be written as

gia (x + ci 4 t, t + 4t) − gia (x, t) = −

gia (x, t) − gia,eq (x, t) + wi Ria , τc

where g a is the concentration distribution function of the reactant a, τc is the relaxation time for concentration field, Ra is the reaction rate of the reactant a. The conservation rule for the concentration is C a (x, t) =

X

gia ,

i

where C a is the concentration of the reactant a. For example, for a three-reactants, a, b, and c system, a separate LBE will be dedicated to each component, where concentrations of other reactants can be computed similarly as above. The definition of the reaction rate can be given by a first-order reaction model

Ra = −kC a , −E

where k is the reaction rate constant and given by given by the Arrhenius equation k = Ae RT . The reaction rate defined in the Arrhenius model is a macroscopic parameter, which is not 37



consistent with the mesoscopic description of concentration given by the concentration LBE. Bresolin and Oliveira 232 developed purely mesoscopic collision theory based on collision theory of Eliason and Hirschfelder to address reaction rate. 233 The C-E or EPDE technique is used on concentration LBE to recover macroscopic equation to determine concentration transport parameter; diffusivity. The diffusivity reads

γ=

c2s

1 τc − 2

Abdollahzadeh et al. compared the concentration profile obtained by the Arrhenius model and collision model for a case of reactive flow in a channel with an obstacle in the middle. 88 The schematic diagram of the reaction problem and comparisons of the results are provided in Fig. (8, 9, 10)

Figure 8: Reaction in a channel with cylinder in the middle (Abdollahzadeh et al.) 88 Reprinted with permission. Copyright Elsevier 2018.

38


Page 38 of 105

Page 39 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 9: Concentration profile distribution in channel with cylinder using (a) Arrhenius macro-model (macroscopic approach) and (b) collision scheme by Bresolin and Oliveira (LBM approach) 232 (Abdollahzadeh et al.) 88 Reprinted with permission. Copyright Elsevier 2018.

39



Figure 10: Comparison of the concentration profile above (top) and below (bottom) of the cylinder in the channel using (a) Arrhenius macro-model and (b) collision scheme by Bresolin and Oliveira 232 (Abdollahzadeh et al.) 88 Reprinted with permission. Copyright Elsevier 2018.

LBM for Heterogeneous Reactions Another example of a reactive flow that can be solved by the LBM is flow through the heterogeneous catalyst. The conversion efficiency of catalytic particles can also be characterized using LBM. 234 Another chemical process such as crystallization can also be simulated by LBM. The population balance equations are solved by using LBM to simulate the crystallization process such as size-dependent growth and nucleation 72,235 . The work developed by Majumder et al. provides in-depth discussions on applying entropic LBM to study crystal growth using population balance model 72 . It is interesting to see that the LBM results were compared to the high-resolution (HR) method for crystal growth and an agreement was found. Moreover, the LBM implementation by Majumder et al. 72 possess faster speed compared to the HR method. Salah et al. 75 used multiphase flow LBM 204 to optimize the 40


Page 40 of 105

Page 41 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


design and performance of polymer electrolyte membrane fuel cell (PEMFC). Such fuel cells contain gas diffusion layers and channels and it is extremely important to study the flow of droplets inside those layers or channels for efficient water drainage. Ostadi et al. studied the properties of the structure of the gas diffusion layers using X-ray nanotomography. 126 Grucelski and Pozorski applied LBM to model coal devolatilization. 127 Very recently, a Gray LBM was proposed by Zhao and Wang 84 to study the effect of interlayer heterogeneity on multi-seam coalbed methane production. Production of coalbed methane from multiple producing zones is of great importance from the economic aspect. The case of the gas flow in coal formation is different from the case of the gas flow in conventional reservoirs. The coalbed formations are highly complex porous media with ultra-low permeability. The authors used Gray LBM and non-uniform grid refinement techniques to increase computational efficiency. For large scale applications, especially where complex geometry has a role to play, uniform or regular grids such as square lattices limits the performance of the LBM solver. Interpolation-supplemented LBM models 124 which uses non-uniform mesh grid enhances the scope for field-scale applications 84 .

Figure 11: Schematic diagram of stirred tank of Rushton turbine (Shu and Yang). 128 Reprinted with permission . Copyright Elsevier 2018.

41



LBM for Flow in Stirred Tanks Another LBM application in heavy industry is flow modeling in stirred tanks. Stirred tanks are used on a large scale in various industries to mix fluids, crystallization, stirring by creating turbulence. 236,237 The turbulence increases the efficiency of the mixing process. Various CFD tools have been used to solve stirred tank problems but are computationally very expensive due to the necessity to solve Poisson equation for large flow region in case of algorithms based on SIMPLE or PISO algorithms and their relatives. LBM does not need to solve Poisson equation, and moreover, it can be easily implemented on GP-GPU devices, 238,239 which is not the case for the classical numerical solvers based on finite differences, elements or volumes. For stirred tanks flow, Immersed Boundary Methods (IBM) 240 were implemented which can effectively address the complexity associated with the moving boundaries. Pioneering studies where LBM was used for simulation of stirred tanks are e.g. 241,242 The LBM simulation of non-Newtonian fluid in the impeller stirred tank is performed by Chen et al. 243 Industrialscale LBM simulations can be found in. 244,245 All those mentioned works used multiple CPU cores to speed up the simulations, later GP-GPU version of LBM algorithms were implemented see e.g. 246,247 Shu et al. used Large Eddy Simulation (LES) to address the turbulence and MRT-LBM to address the flow field, the complex moving boundaries were addressed by IBM. The reader is recommended to go through Ref. 128 for details on LES and IBM. The schematic diagram of the stirrer tank is presented in Fig. (11). The instantaneous velocity profiles are presented in Fig. (12) for three different numerical approaches. It can be verified from the profiles that MRT-LBM along with IBM and LES performs better compared to ANSYS Fluent with body-fitted mesh boundary walls. The velocity profile for case (a) is rather unphysical. It is evident that velocity is large in the vicinity of impellers and uniform structure of the flow field is seen (unphysical). Some non-uniformity in velocity field near to impellers blades has been seen in case (b). But In this case, the periodic flows around the baffles are absent. The

42


Page 42 of 105

Page 43 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


case (c), with MRT-LBM and LES, produces better results, the non-uniformity in velocity profile can be seen both close to impellers and baffles, and periodic flows are also present around the baffles. A similar study of flow fields for different impellers based on ANSYS Fluent simulations together with experiments which reveal the performance of the impeller on the mixing efficiency can be found in. 248

Figure 12: Comparison of instantaneous velocity profile computed by various numerical tools. (a) Flow solver- ANSYS Fluent, wall boundary condition- body-fitted mesh, moving boundary treatment- sliding mesh, turbulence model- Smagorinsky LES. (b) Flow solver- ANSYS Fluent, wall boundary condition- body-fitted mesh, moving boundary treatment- sliding mesh, turbulence model- transient SST k-ω. (c) Flow solver- MRT-LBM, wall boundary condition- IBM, moving boundary treatment- IBM, turbulence model- LES. (Shu and Yang 2018). 128 Reprinted with permission. Copyright Elsevier 2018.

LBM for Flow in Porous Media Flow in porous media is a complete research world in itself. A large number of engineering streams deal with the flow in porous media. The flow of bone-marrow through bones, exploration, and production of hydrocarbons from reservoirs, flow through catalysts, ground-water exploration and contamination, filtration, and transport of heat and mass through packed beds are just some of the examples termed as flow in porous media. From the industrial application point of view, exploration and production of hydrocarbons from oil and gas reservoirs can be considered as the most popular field of research in porous media. Modeling of flow through porous media is highly challenged by the complex pore-space geometry present in 43



them. Another challenge associated with such flows is the microscopic interactions between solid-liquid (wettability) and reactive fluid-solid interaction in porous rocks for single-phase flow, and fluid-fluid, fluid-solid interactions in multiphase multicomponent flows. 249 Therefore, to solve flow in porous media effectively, one needs to have (1) efficient techniques to address complex boundary conditions, (2) efficient techniques to define physicochemical properties of fluid-solid system, i.e. wettability in the numerical model, (3) clear definition of static or moving fluid-fluid interface structures, and (4) computation efficiency. After paying close attention to the previous four points, LBM comes up as the most viable numerical tool to address flow in porous media because (1) LBM can incorporate complex boundary conditions naturally in the LBE so no need for mesh refinement (computationally costly), (2) moving boundaries can be taken into account effectively using IBM, (3) due to the kinetic-nature of the LBM, fluid-fluid interfaces can also be defined efficiently capturing maximum of the microscopic details, and (4) the ease of GPU computation implementation of the LBM code saves one from using tens of CPUs and longer simulation time (computationally expensive) because of local nature of data communication in LBM (local nature of collision processes) parallelization is convenient. Before starting modeling flow in porous media using LBM one should note that there exist two approaches to solve such problems, (1) Representative Element Volume (REV) 250,251 approach and (2) Pore-Scale approach 102 . 52,55 The REV approach does take into consideration the petrophysical properties such as porosity and permeability, and Pore-Scale approach are micro-level flows in which microscopic behavior such as flow through capillary channels, fluid-fluid contact systems, and interactions are observed. For pore-scale modeling applications, computational efficiency has been an issue for LBM. Simulating large scale porous samples is highly memory demanding and even massively parallel implementations because the run-time scale is very small for large real flow rate, cannot predict flow through capillary accurately 102 . These limitations reduce the scope of the LBM to correctly study pore-scale multiphase flow in large samples and determining relative permeability is a challenging task. 44


Page 44 of 105

Page 45 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Though by using multi GP-GPU computing this problem can be addressed. Single-phase flow can be predicted efficiently by LBM without making much serious effort regarding computational cost. But real-world applications deal with the multiphase flow in most of the cases, and massively parallel implementation of the LBM on 3D pore-space of porous media generated by imaging techniques such as micro-CT (digital rock analysis) can predict multiphase flow efficiently, but at a high computational cost. Digital rock physics is becoming very popular recently as it takes into produces more realistic flow characteristics in the porous samples. Hao et al. 252 presented one of the first implementations of multiphase flow LBM to study relative permeability in porous media using the pore-scale approach. But the authors did not use real pore-space imaging techniques, instead, the packed-sphere bed and carbon paper gas diffusion layer were used as porous media. This work is of particular importance as it helps construct a good base for relative permeability prediction in porous media with great insights, which can be implemented anytime on images acquired by imaging techniques. For an immiscible two-phase flow system, due to interfaces and surface tension, free-energy multiphase flow LBM approach suits better as it is also thermodynamically consistent. The authors used the free-energy function Cahn-Hilliard theory to model surface tension. The details of the multiphase model used in this study can be found in a separate paper also by Hao and Cheng 253 . The LBM simulation results of the non-wetting phase distribution and wetting phase velocity vectors in sphere-packed bed is shown in Fig. (13). Due to low Capillary number and Re number, Capillary force is the driving force and nonwetting phase (blue color) gets dispersed as trapped droplets inside the pores at small saturation (due to low Capillary number and Re). The authors further state that the pore-network model is incapable of producing such results. The relative permeabilities of the nonwetting and wetting phase are shown in Fig. (14). The LBM approach showed great agreement with the experimental results 254 for the same sphere-packed bed setup. The readers interested in pore-scale applications for more complex pore-space geometry of real minerals (digital rock 45



analysis) can refer to the recent work by Shah et al. 104 and its references.

Figure 13: Body centered sphere bed with velocity distribution and non-wetting phase. Porosity: 0.55, non-wetting phase contact angle: 180◦ and phase saturation: 0.22. The body force at a Capillary number of the order of 10−6 and Re number: 10−4 . Fig. from Hao and Cheng (2010). 252 Reprinted with permission. Copyright Elsevier 2010.

46


Page 46 of 105

Page 47 of 105 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 14: Relative permeabilities for nonwetting and wetting phase and comparison with the experimental results from Bryant and Blunt 254 . kr,w and kr,nw are the relative permeabilities for wetting and nonwetting phase, respectively. SN W is the saturation of the nonwetting phase. Fig. from Hao and Cheng (2010). 252 Reprinted with permission. Copyright Elsevier 2010.

The BGK, MRT, Cascaded or any other collision scheme can be chosen to solve the mass and momentum conservation equations. We already highlighted earlier, that the location of the pore walls depends on the relaxation parameter, therefore for higher stability and accuracy, MRT and cascaded LBM are the better choices compared to the BGK-LBM. Let us throw some light on LBM implementations to address diverse problems. It must be noted that flow in porous media, based on the Re number, can be divided into two categories, nonDarcy flow (large Re) and Darcy flow (small Re). 249 The flow rate through porous media is defined by Darcy law, which reads q=

k ∂P , µ ∂x

where q is the flow rate, P is the pressure, k is the permeability and µ is the dynamic viscosity. Therefore, the regimes Darcy and non-Darcy can be described based on the validity of this

47



expression. For flow regime where Re > 1, the flow field starts to deviate from the linear behavior of Darcy law. 110,255 For a flow with Re

Lattice Boltzmann Methods for Industrial Applications - American

Recommend Documents