Comparison of multiphase SPH and LBM approaches for the simulation of intermittent flows
Thomas Douillet-Grellier, Sébastien Leclaire, David Vidal, François Bertrand, Florian De Vuyst
CComparison of multiphase SPH and LBM approaches forthe simulation of intermittent flows
Thomas Douillet-Grellier a,b, ∗ , Sébastien Leclaire c,d , David Vidal d , FrançoisBertrand d , Florian De Vuyst e a CMLA, CNRS, ENS Paris-Saclay, Université Paris-Saclay, Cachan, France. b Total S.A., Tour Coupole, Paris La Défense, France c Department of Mechanical Engineering, Polytechnique Montréal, Québec, Canada d URPEI, Department of Chemical Engineering, Polytechnique Montréal, Québec, Canada e LMAC, Université de Technologie de Compiègne, Compiègne, France
Abstract
Smoothed Particle Hydrodynamics (SPH) and Lattice Boltzmann Method(LBM) are increasingly popular and attractive methods that propose efficientmultiphase formulations, each one with its own strengths and weaknesses. Inthis context, when it comes to study a given multi-fluid problem, it is helpfulto rely on a quantitative comparison to decide which approach should be usedand in which context. In particular, the simulation of intermittent two-phaseflows in pipes such as slug flows is a complex problem involving moving andintersecting interfaces for which both SPH and LBM could be considered. It isa problem of interest in petroleum applications since the formation of slug flowsthat can occur in submarine pipelines connecting the wells to the production fa-cility can cause undesired behaviors with hazardous consequences. In this work,we compare SPH and LBM multiphase formulations where surface tension ef-fects are modeled respectively using the continuum surface force and the colorgradient approaches on a collection of standard test cases, and on the simula-tion of intermittent flows in 2D. This paper aims to highlight the contributionsand limitations of SPH and LBM when applied to these problems. First, wecompare our implementations on static bubble problems with different densityand viscosity ratios. Then, we focus on gravity driven simulations of slug flowsin pipes for several Reynolds numbers. Finally, we conclude with simulationsof slug flows with inlet/outlet boundary conditions. According to the resultspresented in this study, we confirm that the SPH approach is more robust andversatile whereas the LBM formulation is more accurate and faster.
Keywords:
SPH, LBM, multiphase, boundary conditions, slug
1. Introduction
Two-phase flows in pipes are encountered in various industrial applicationsand research areas. For instance, gas-liquid flows are present in chemical andheat transfer systems such as evaporators or refrigerators, but they can also befound in the transportation of oil and gas through pipelines in the petroleum ∗ [email protected] Preprint submitted to arXiv March 5, 2019 a r X i v : . [ phy s i c s . c o m p - ph ] M a r ndustry. In particular, intermittent flow patterns (also known as slug or plugflow [1]) are undesirable in pipeline networks. Those slug patterns, that canmeasure up to tens of meters, are known to damage facilities (separator flooding,compressor starving, water hammer phenomenon, vibrations and fatigue) andto reduce flow efficiency. Therefore, it is crucial to be able to predict whetheror not a slug flow regime will occur. If it does, the quantities of interest are thesize of slugs, their transit time and their frequency.In this context, numerical simulations have emerged as a tool of choice [2, 3]to understand the formation of intermittent flow patterns and to predict theappearance of a slug flow regime. At the industrial level, the first simula-tions [4, 5, 6] were done in the mid 80’s. Nowadays, in the oil and gas indus-try, two commercial software products are competing for slugging simulation :OLGA developed by SPT group [7] and LedaFlow, proposed by Kongsberg [8].A detailed comparison of both codes [9] concludes that although performingequally well on simple cases, they have trouble to simulate complex cases with adominant gas phase. From an academic perspective, different models and meth-ods have been used for slug flow modeling, for example volume-of-fluid [10, 11],level-set [12, 13], lattice Boltzmann method (LBM) [14], smoothed particle hy-drodynamics (SPH) [15, 16] or phase field [17, 18], but these are applied onmicrofluidic problems for the most part. In this work, we will focus on two ofthem : SPH and LBM, because while they are very different numerical methods,both in their origin and in their nature, they have shown a strong potential tomodel multiphase flows [19, 20, 21].SPH is a Lagrangian meshfree method introduced in the late 70’s for astro-physics applications [22, 23] and later expanded to standard fluid mechanics.SPH can be seen from two different perspectives. On one hand, it is a way todiscretize partial differential equations such as the Navier-Stokes equations. Onthe other hand, it is a discrete Hamiltonian system composed of material pointsof constant mass that are tracked in time. Its pure meshless nature allows to getrid of many issues associated with meshing. However, it comes at some expensestoo. Since there is no connectivity between particles, a neighbor searching pro-cedure has to be carried for every particle at every time step, which is a seriousbottleneck for code efficiency. Among the numerous applications of SPH, we canmention astrophysics [24], hydrodynamics [25], geophysics [26, 27, 28, 29] andcomputer graphics [30]. Some extensive reviews have been published [31, 32, 33].LBM originates from two distinct approaches: the kinetic gas theory withdiscrete velocities and the lattice gas cellular automata method (LGCA orLGA [34]). In the late 80’s, probability density functions of particles were intro-duced within LGCA, giving birth to LBM. It consists in solving the Boltzmannequation in a discrete velocity space, which has been proven to be equivalentto solving the Navier-Stokes equations in the limit of low Mach and Knudsennumbers (as can be shown by a multiscale Chapman-Enskog expansion [35]). Inpractice, LBM distinguishes itself from other methods for several reasons. First,the LBM operates on an uniform lattice (mostly square or hexagonal lattices).Then, LBM offers a local algorithm involving simple arithmetic operations withno differential terms, which makes it short, fast and particularly well suitedfor parallel execution [36]. Finally, traditional CFD methods are based on thecalculation of macroscopic variables (velocity, pressure, density) whereas LBMtracks the evolution of probability distribution functions of particles [37]. LBMhas been used for decades for flow simulations in complex geometries, especially2n porous media [38, 39, 40, 41].Our goal here is to propose a comparison of the multiphase SPH formulationpresented in [42] and the color-gradient multiphase LBM formulation introducedin [43], on a collection of standard 2D test cases and on the simulation of slugflow regimes with periodic and inlet/outlet boundary conditions. To the best ofour knowledge, this the first time such a comparison is presented.We first detail the multiphase LBM and SPH formulations used in this workincluding surface tension models and boundary conditions. Then, we compareboth approaches on static bubble tests with different density and viscosity ra-tios. Finally, we extend the comparison to two cases of slug flows, one inducedby gravity with periodic boundary conditions and the other one based on in-let/outlet boundary conditions.In addition, we provide in Appendix Appendix A a comparison of SPHand LBM on two test cases where one is a single-phase flow and the other doesnot involve surface tension : the lid-driven cavity flow and the Rayleigh-Taylorinstability problems.
2. LBM immiscible multiphase model
Four main multiphase formulations are available for LBM : the pseudo-potential model [44], the free energy model [45], the mean field model [46]and the color gradient model [47]. We recommend the reading of [48, 49] forthose looking for a detailed comparison of these techniques. In this work, wechoose to work with the color gradient model because among the diffuse in-terface approaches, it is the one with a thin interface compared to the pseudopotential approach for example. In addition, the pseudo-potential model iscumbersome to use and to parametrize because there is coupling between thebasic parameters [49]. The free energy model requires solving a Poisson equa-tion at every time step which is time consuming and the mean field approachis limited to small density ratios [50]. Moreover, the color gradient model ben-efits from the large body of verification and validation cases available in theliterature [51, 52, 53, 54, 55, 56, 49].The present LBM approach is the two-phase model introduced in [43] andcompleted with the improvements proposed in [52, 51] for the recoloring operatorand the color gradient. In addition, the contact angle ajustment is based on [56,57] and the corrective procedure to properly recover Navier-Stokes equations isborrowed from [58]. We work with sets of distribution functions, one for eachfluid, moving on a D2Q9 lattice. The associated velocity vectors are c i with i ∈ [0 , . As traditionally done in LBM, we set the lattice time step ∆ t andthe lattice space step ∆ x to . We can then define the 2D velocity vectors asfollows : c i = (0 , , i = 0[sin( θ i ) , cos( θ i )] , i = 1 , , , θ i ) , cos( θ i )] √ , i = 2 , , , , with θ i = π − i ) . (2.1)The distribution functions for a fluid of color k (e.g. k = r for red and k = b for blue) are denoted N ki ( x , t ) , while N i ( x , t ) = N ri ( x , t ) + N bi ( x , t ) is used forthe color-blind distribution function. In the rest of this section, the integersubscript i is varying between [0 , while k is referring to the color blue b or3ed r of the fluid. The lattice Boltzmann equation that describes the evolutionof the system is then : N ki ( x + c i , t + 1) = N ki ( x , t ) + Ω ki (cid:0) N ki ( x , t ) (cid:1) , (2.2)where the collision operator Ω ki is the result of the combination of three suboperators (as previously done in [59]) : Ω ki = (Ω ki ) (3) (cid:104) (Ω ki ) (1) + (Ω ki ) (2) (cid:105) . (2.3)The Eq. (2.2) is solved using four consecutive steps which make use of thefollowing operators :1. Single phase collision step (see Sec. 2.1) : N ki ( x , t ∗ ) = (Ω ki ) (1) (cid:0) N ki ( x , t ) (cid:1) . (2.4)2. Perturbation step (multiphase collision 1/2) (see Sec. 2.2) : N ki ( x , t ∗∗ ) = (Ω ki ) (2) (cid:0) N ki ( x , t ∗ ) (cid:1) . (2.5)3. Recoloring step (multiphase collision 2/2) (see Sec. 2.3) : N ki ( x , t ∗∗∗ ) = (Ω ki ) (3) (cid:0) N ki ( x , t ∗∗ ) (cid:1) . (2.6)4. Streaming step : N ki ( x + c i , t + 1) = N ki ( x , t ∗∗∗ ) . (2.7)We will now examine in detail these steps as well as the specific treatmentsfor the imposition of static contact angles and boundary conditions. The first sub operator, (Ω ki ) (1) of Eq. (2.4), is the usual BGK operatorfor single phase LBM. The distribution functions are relaxed towards a localequilibrium as follows : (Ω ki ) (1) ( N ki ) = N ki − ω eff (cid:16) N ki − N k ( e ) i (cid:17) , (2.8)where ω eff is the effective relaxation rate. In practice, we first calculate the fluiddensities based on the th moment of the distribution functions : ρ k = (cid:88) i N ki = (cid:88) i N k ( e ) i , (2.9)where N k ( e ) i are the equilibrium distribution functions. The total fluid density isgiven by ρ = ρ r + ρ b , while the total momentum is computed as the st momentof the distribution functions : ρ u = (cid:88) i (cid:88) k N ki c i = (cid:88) i (cid:88) k N k ( e ) i c i , (2.10)4here u is the density weighted arithmetic average velocity of the fluid. Theequilibrium distribution functions N k ( e ) i are defined in [43] by : N k ( e ) i ( ρ k , u , α k ) = ρ k (cid:18) φ ki + W i (cid:20) c i · u + 92 ( c i · u ) −
32 ( u ) (cid:21)(cid:19) . (2.11)These equilibrium distribution functions N k ( e ) i are chosen to satisfy massand momentum conservation and are based on the Maxwell-Boltzmann distri-bution. The weights W i are those of a standard D2Q9 lattice while the values φ ki depend on the density ratio. They are expressed as follows : W i = / , i = 01 / , i = 1 , , , / , i = 2 , , , φ ki = α k , i = 0(1 − α k ) / , i = 1 , , , − α k ) / , i = 2 , , , . As introduced in [60] for two-phase flows, the density ratio between red andblue fluids is γ , and must be computed as follows to obtain a stable interface : γ = ρ r ρ b = 1 − α r − α b , (2.12)where the superscript refers to the initial value of the density. Besides, thepressure of the fluid of color k is : p k = 3 ρ k (1 − α k )5 = ρ k ( c ks ) . (2.13)In Eq. (2.12), one of the α k is actually a free parameter. In practice, if weassume that the least dense fluid is the blue one, we set a positive value for α b > ( . in this paper) and so we are certain that < α b ≤ α r < usingEq. (2.12). These parameters define the value of the sound speed c ks in eachfluid of color k [43].The effective relaxation parameter ω eff is chosen so that the evolution Eq. (2.2)allows to recover the macroscopic Navier-Stokes equations for a single-phase flowin the single-phase areas. This parameter depends on the fluid kinematic vis-cosity ν k through the following formula: ω k = 1 (cid:14)(cid:0) ν k + (cid:1) . However, when theviscosities of the fluids are different, the relaxation parameters are also differentand an interpolation procedure is needed to define an effective relaxation param-eter ω eff at the interface. It is common to use a quadratic interpolation [43, 60].In order to detect in which area we are located (pure red fluid, pure blue fluidor interface), it is necessary to introduce a color field : ψ = (cid:18) ρ r ρ r − ρ b ρ b (cid:19)(cid:30)(cid:18) ρ r ρ r + ρ b ρ b (cid:19) . (2.14)The color field ψ lies between − and +1 . Within an interface, the colorfield is strictly between − and +1 whereas it is equal to − or +1 when locatedin an area that contains respectively only red fluid or blue fluid. The relaxationfactor ω eff in Eq. (2.2) is then evaluated as follows : ω eff = ω r , ψ > δf r ( ψ ) , δ ≥ ψ > f b ( ψ ) , ≥ ψ ≥ − δω b , ψ < − δ , (2.15)5n which δ is a free parameter and f r ( ψ ) = χ + ηψ + κψ ,f b ( ψ ) = χ + λψ + νψ , (2.16)with : χ = (2 ω r ω b )/( ω r + ω b ) ,η = (2( ω r − χ ))/ δ,κ = − η /(2 δ ) ,λ = (2( χ − ω b ))/ δ,ν = λ /(2 δ ) . (2.17)The parameter δ influences the thickness of the interface when the fluid viscosi-ties are different. The larger δ , the thicker the fluid interface. There is a tradeoff to find between robustness and interface sharpness. It is set to . for allsimulations in this paper. If the fluid viscosities are identical, δ is ineffective,because ω eff = ω r = ω b . Alternatively, for the first sub operator, (Ω ki ) (1) , one can use the MultipleRelaxation Time (MRT) operator instead of the BGK operator. The MRTapproach is more stable than its BGK counterpart [61, 58]. It reads as follows: (Ω ki ) (1) ( N ki ) = N ki − M − SM (cid:16) N ki − N k ( e ) i (cid:17) , (2.18)where S is a diagonal matrix given by : S = diag ( s , s , s , s , s , s , s , s , s ) . (2.19)The elements s i , i ∈ [0 . . . are the relaxation parameters. Following [54], ∀ i ∈ [0 . . . , s i = λω eff while s = s = ω eff . As in [54], we choose λ = 4 / .Moreover, M is a linear orthogonal transformation matrix that projects thedistribution functions into the moment space. The discrete moment matrices M and M − are explicitly given in Appendix Appendix B for a D2Q9 lattice. It has been emphasized in several papers [62, 63, 53] that the color gradientmodel does not fully recover the Navier-Stokes equations. An unwanted errorterm arises in the momentum equations when the density ratio is not one. Dif-ferent techniques have been proposed to attenuate this issue [53, 64, 63, 58]. Inthe present work, we use the correction introduced in [58] for the MRT approach.It consists in two additions. First, a modified equilibrium distribution functionsbased on the rd order Hermite expansion of the Maxwellian distribution [65, 66]is used instead of Eq. (2.11). It is expressed as follows : N k ( e ) i ( ρ k , u , α k ) = ρ k φ ki + ρ k W i (cid:18)(cid:20) c i · u (cid:20) (cid:0) c ks ) − (cid:1) (3 c i − (cid:21) + 92 ( c i · u ) −
32 ( u ) (cid:21)(cid:19) . (2.20)6econd, a source term U k is added in Eq. (2.18). It reads : U k = M − C k , (2.21)where C k = (cid:2) , C k , , , , , , C k , (cid:3) T . The components C k and C k are com-puted as follows : C k = 3(1 − s / ∂ x Q kx + ∂ y Q ky ) ,C k = 3(1 − ω eff / ∂ x Q kx + ∂ y Q ky ) , with Q kx = (1 . α k − . ρ k u x ,Q ky = (1 . α k − . ρ k u y . (2.22)In particular, the derivatives are evaluated using a -point isotropic finite dif-ference approximation described shortly afterwards in Eq. (2.24). In the color gradient model, surface tension forces are introduced by meansof a perturbation operator shown in Eq. (2.5) [67, 43, 47]. To begin, it is neededto introduce the color gradient F of the color field ψ (see Eq. (2.14)) thatapproximates the fluid-fluid interface normal. It is written as : F = ∇ ψ. (2.23)In this work, a bi-dimensional S I (spatial order S = 2 , isotropic order I = 4 )discrete gradient operator is used [68]. As shown in [51], this kind of gradientoperator enhances the accuracy of the color-gradient model significantly whilehaving the advantage to only rest on the nearest neighboring nodes. It takesthe following form : ∇ f ( x ) ≈ (cid:88) i ξ i c i f ( x + c i ) , with ξ i = , i = 01 / , i = 1 , , , / , i = 2 , , , . (2.24)The perturbation operator, for the fluid k , is defined by : (Ω ki ) (2) ( N ki ) = N ki + (cid:88) ll (cid:54) = k A | F kl | (cid:20) W i ( F · c i ) | F | − B i (cid:21) , (2.25)where : B i = − / , i = 02 / , i = 1 , , , / , i = 2 , , , . (2.26)The parameter A which can vary in space and time handles the couplingbetween both fluids and is therefore linked with the surface tension coefficient σ . It is possible to predict the surface tension σ between the two fluids using onlythe basic parameters of the model. As explained in [52], knowing the form of theexpression describing the surface tension and performing simulations on planarinterfaces, one can derive an expression that links the surface tension across aninterface to the model parameters. For isotropic color gradients defined as inEq. (2.24), the surface tension is set as follows : σ = 49 Aω eff . (2.27)7f σ and ω eff are fixed, one can obtain the value of A . Note that Eq. (2.27) isnot universal and is susceptible to change if one uses a different color gradientor a different gradient operator. It has been shown in [43] that the perturbationoperator allows to recover, within the macroscopic limit, the capillary stresstensor responsible for the surface tension forces in the macroscopic two-phaseflow equations. The recoloring operator (Ω ki ) (3) of Eq. (2.6) is used to maximize the amountof fluid k at the interface that is sent to the k fluid region. It guaranteesthe fluid’s immiscibility. It respects the principles of mass and momentumconservation. The form of recoloring used in this paper is a combination ofideas taken from [69] and [70] and fully developed in [52]. It reads : (Ω ri ) (3) ( N ri ) = ρ r ρ N i + β ρ r ρ b ρ cos( ϕ i ) (cid:80) k N k ( e ) i ( ρ k , , α k ) , (Ω bi ) (3) ( N bi ) = ρ b ρ N i − β ρ r ρ b ρ cos( ϕ i ) (cid:80) k N k ( e ) i ( ρ k , , α k ) , (2.28)where β ∈ [0 . . . is a parameter controlling the interface thickness [69] thatwill be set to . unless otherwise stated. cos( ϕ i ) = ( c i · F )/( | c i || F | ) is thecosine of the angle between the color gradient F and the lattice direction vector c i . Note that for i = 0 or | F | = 0 , there is a division by . In such a case, weset the whole term equal to 0 to respect mass conservation. Based on [56], the imposition of a contact angle is performed using themethod described in [57]. In order to properly describe the wetting boundaryconditions, we divide the lattice nodes in two categories C f and C s , each categorybeing also subdivided in two subcategories C fs , C ff , C sf and C ss .- C f : the set of fluid lattice nodes – C fs : fluid lattice nodes in contact with at least one solid lattice node – C ff : fluid lattice nodes not in contact with any solid lattice node- C s : the set of solid lattice nodes – C sf : solid lattice nodes in contact with at least one fluid lattice node – C ss : solid lattice nodes not in contact with any fluid lattice nodeWhen computing the color gradient in the fluid (i.e. for lattice nodes ∈ C f ),the knowledge of the color field at the boundary is required (i.e. for lattice nodes ∈ C sf ) because Eq. (2.24) involves the neighboring lattice nodes. Therefore, weneed to extrapolate the color field to the boundary nodes, we do so using thefollowing expression : ∀ x ∈ C sf , φ ( x ) = (cid:88) i x + c i ∈ C fs W i φ ( x + c i ) c iα (cid:44) (cid:88) i x + c i ∈ C fs W i . (2.29)8t is now possible to evaluate the orientation of the color gradient in the fluidusing the expression hereafter : n ∗ = ∇ φ ( x ) |∇ φ ( x ) | . (2.30)In [57], they use an th order isotropic stencil to compute the surface normal n s . In the subsequent simulations, boundary surfaces are flat and normalsare known so we directly input the exact value. The correct color gradientorientation n depends on the prescribed contact angle θ and is evaluated asfollows : n = n , if D < D n , if D > D n s , if D = D , , (2.31)where D and D are the Euclidean distances between the current unit normalvector n ∗ and the two possible theoretical unit normal vectors n and n of theinterface at the contact line and are given by : D = | n ∗ − n | ,D = | n ∗ − n | , (2.32)with : n = (cid:0) n sx cos θ − n sy sin θ, n sy cos θ + n sx sin θ (cid:1) , n = (cid:0) n sx cos θ + n sy sin θ, n sy cos θ − n sx sin θ (cid:1) . (2.33)Finally, once n is known, the corrected color gradient value is computed bytaking : ∇ φ ( x ) = |∇ φ ( x ) | n . (2.34) An overview of the available techniques to impose boundary conditions forsingle phase LBM can be found in [71]. For multiphase LBM, the literatureis less abundant. In particular, the case of inlet/outlet boundary conditionsis particularly difficult because specific treatments are needed to handle theinterface when the fluids are entering and/or leaving the domain. Being able tohave efficient inlet/outlet boundary conditions is attractive because it extendsthe range of two-phase flow simulations that could be explored [72, 73, 74, 75, 76]and is mandatory for open channels. Three kinds of boundary conditions areused in this work :1. No slip boundary conditions are imposed using full bounceback [77]. Then,free slip boundary conditions are also imposed using full bounceback except thatthe diagonally traveling distribution functions are sent forward along the wallrather than reflected back the way they came. Finally, velocity and/or pressureboundary conditions are imposed following [78].9. Periodic boundary conditions is a very useful tool in computational simula-tions because it allows to reduce the size of the simulated domain. The im-plementation of these boundary conditions consists in sending the distributionfunctions that are leaving the domain on one end to the other end of the domainas if the two sides of the domain were connected.3. Inlet/outlet boundary conditions are achieved using a modified version ofZou-He boundary conditions [78]. The approach detailed in this paper is verysimilar to what is described for two-phase pressure boundary conditions in [74].First, we will describe how we inject two phase flows with two different velocities u b inlet and u r inlet from the north wall. The prescribed velocity fields v b, in and v r, in are designed so that v b, in = u b inlet on blue boundary lattice nodes and v b, in = 0 on red boundary lattice nodes. Similarly, v r, in = u r inlet on red boundary latticenodes and v r, in = 0 on blue boundary lattice nodes. It is then possible togenerate a color-blind prescribed velocity field v in = v b, in + v r, in . Following theclassic Zou-He procedure described in [78], we can then compute the modifieddensity and distribution functions. It reads : ρ in = v in y ( N + N + N + 2 ( N + N + N )) ,N = N − ρ in v in y + H in ,N = N + ( N − N ) − ρ in v in x − ρ in v in y − H in ,N = N + ( N − N ) + ρ in v in x − ρ in v in y − H in , (2.35)where H in is a corrective term that depends if we use Eq. (2.11) or Eq. (2.20)for the equilibrium. In practice, to derive Zou-He boundary conditions, one hasto solve a linear system where one term comes from the equilibrium distributionfunction. If we use Eq. (2.20), we obtain extra terms compared to the classicalZou-He approach due to the Hermite expansion. It is computed as follows : H in = if we use Eq. (2.11) , (cid:104) φ +12 ρ in ( c bs ) + (cid:16) − φ +12 (cid:17) ρ in ( c rs ) − ρ in (cid:105) v in y if we use Eq. (2.20) . (2.36)It is then needed to redistribute these quantities in function of the color fieldvalue : ρ b = φ +12 ρ in , ρ r = (cid:16) − φ +12 (cid:17) ρ in ,N b = φ +12 N , N r = (cid:16) − φ +12 (cid:17) N ,N b = φ +12 N , N r = (cid:16) − φ +12 (cid:17) N ,N b = φ +12 N , N r = (cid:16) − φ +12 (cid:17) N . (2.37)Second, we will describe how we impose a constant pressure p out at the outletlocated on the south wall. The corresponding prescribed densities are ρ b,out = p out (cid:14) ( c bs ) and ρ r,out = p out (cid:14) ( c rs ) . The color-blind prescribed density is then10 out = ρ b,out + ρ r,out . In addition, we also have : v out x = 0 v out y = ρ out ( N + N + N + 2 ( N + N + N )) ,N = N + ρ out v out y + H out ,N = N − ( N − N ) − ρ out v out x + ρ out v out y − H out ,N = N − ( N − N ) + ρ out v out x + ρ out v out y − H out , (2.38)where H out is evaluated as follows : H out = if we use Eq. (2.11) , − (cid:104) φ +12 ρ out ( c bs ) + (cid:16) − φ +12 (cid:17) ρ out ( c rs ) − ρ out (cid:105) v out y if we use Eq. (2.20) . (2.39)We can then redistribute these quantities similarly with what was done inEq. (2.37) : ρ b = φ +12 ρ out , ρ r = (cid:16) − φ +12 (cid:17) ρ out ,N b = φ +12 N , N r = (cid:16) − φ +12 (cid:17) N ,N b = φ +12 N , N r = (cid:16) − φ +12 (cid:17) N ,N b = φ +12 N , N r = (cid:16) − φ +12 (cid:17) N . (2.40)A test case has been setup to check the performance of these boundaryconditions. Initial and final configurations can be found in Fig. 1. In Fig. 2 isshown the evolution of the inlet velocities and outlet pressure with the numberof iterations. Note that quantities have been averaged along the height of thepipe. We can see that that we are recovering the prescribed values at theinlet and at the outlet after a transient period with a maximum error ≤ .In Fig. 3, we see the distribution of the color field, the inlet velocity and theoutlet pressure along the height of the pipe. It is possible to observe a velocitypeak and a pressure peak located at the interface. This is likely due to thefact that fluids are mixed at the interface resulting in governing equations notbeing properly solved at this location. Moreover, slight discrepancies can beseen at the walls due to boundary conditions. Overall, the previously describedboundary conditions are giving satisfactory results and will be used later in thepaper.
3. SPH immiscible multiphase model
In this section, we will describe in details the SPH immisible multiphasemodel used in this work. This formulation and the associated open boundaryconditions are taken from [42, 16].
For an incompressible fluid with a constant viscosity, the mass and momen-tum equation (completed with the equation of state) in a Lagrangian system11 H v inlight v inheavy = 6 v inlight p out Light FluidHeavy Fluid (a) (b)Figure 1: (a) Initial configuration sketch. (b) Phases distribution after iterations.
0k 10k 20k 30k 40k . . N iter v li g h t / v i n li g h t (a) Inlet light phase velocity
0k 10k 20k 30k 40k . . N iter v h e a vy / v i nh e a vy (b) Inlet heavy phase velocity
0k 10k 20k 30k 40k . . N iter p / p o u t (c) Outlet pressureFigure 2: Evolution of selected quantities with the number of iterations − − . . Y/H φ (a) Color field . . . . Y/H v / v i n li g h t (b) Inlet velocity . . Y/H p / p o u t (c) Outlet pressureFigure 3: Variation of selected quantities along the height of the pipe at steady state DρDt = − ρ ∇ · u , (3.1) D u Dt = − ∇ pρ + ν ∇ u + F st ρ + g , (3.2) p = c ρ γ (cid:20)(cid:18) ρρ (cid:19) γ − (cid:21) + p , (3.3)with u fluid velocity, ρ fluid density, µ the fluid viscosity, ν = µρ the kinematicviscosity, p fluid pressure, µ fluid dynamic viscosity, g gravity, c fluid speedof sound (here constant), γ fluid adiabatic index, ρ fluid initial density, p background pressure, F st is the surface tension force and D/Dt denotes thematerial derivative following the motion.The Tait’s equation of state (3.3) is added to Eqs. (3.1)-(3.2) to close thesystem. In this work, γ will be set to . for all fluids considered in all subse-quent simulations. This is the so-called Weakly Compressible SPH formulation(WCSPH). Just like LBM, it is not a truly incompressible approach since thedensity is allowed to vary. This artificial compressibility has to be as weak aspossible and is controlled by the speed of sound c . In this paper, given a refer-ence length L ref and a reference speed U ref , we used the following formulas [79]to set a value for c and p : c α = max (cid:18) U ref √ ∆ ρ , (cid:113) | g | L ref ∆ ρ , (cid:113) σ αβ ρ α L ref , (cid:113) µ α U ref ρ α L ref ∆ ρ (cid:19) , ∀ α ∈ { , . . . , N phases } ,p = max α ∈{ ,...,N phases } c α ρ α γ α , (3.4)with ∆ ρ = 0 . to enforce (not strictly) a maximum variation of of thedensity field and σ αβ the surface tension coefficient between phases α and β .The integer N phases is the number of different phases.In order to model surface tension between fluids, an interaction force F st isadded to the momentum Eq. (3.2). Following [80], the continuum surface stressmethod introduces a surface tension force per unit volume that is expressed asthe divergence of the capillary pressure tensor : F st = −∇ · Π , (3.5)with Π the capillary pressure tensor defined by : Π = (cid:88) α,β | α<β Π αβ , (3.6)where α, β ∈ { , . . . , N phases } and Π αβ is expressed as : Π αβ = − σ αβ (cid:0) I − ˜ n αβ ⊗ ˜ n αβ (cid:1) δ αβ , (3.7)with ˜ n αβ the unit normal vector from phase α to phase β , σ αβ the surfacetension coefficient between phase α and phase β , δ αβ a well-chosen surface deltafunction and I the identity matrix. For example, in the case of a three-phasesystem with a wetting phase s , a non-wetting phase n and a solid phase s , thestress tensor reads Π = Π ns + Π ws + Π nw .13 .2. SPH formulation The SPH formulation adopted for this paper is identical to the one usedin [16]. We will not here recall in detail the SPH discretization process andwill assume that the reader already has experience with SPH. For those whoare looking for an exhaustive description of this method, we recommend forexample the reading of [81].An SPH discretization consists of a set of points with fixed mass, whichpossess material properties and interact with its neighboring particles througha weighting function (or smoothing kernel) [23]. A particle’s support domain, Λ , is defined by its smoothing length, h , which is the radius of the smoothingkernel W . In all simulations presented in this paper, h = 2∆ r where ∆ r isthe initial particle spacing. To obtain the value of a function at a given par-ticle location, values of that function are found by taking a weighted (by thesmoothing function) interpolation from all particles within the given particle’ssupport domain. An analytical differentiation of the smoothing kernel is usedto find gradients of this function.The interpolated value of a function A at the position x a of particle a canbe expressed using SPH smoothing as : A ( x a ) = (cid:88) b ∈ Λ a A b m b ρ b W ( x a − x b , h ) = (cid:88) b ∈ Λ a A b m b ρ b W ab , (3.8)in which A b = A ( x b ) , m b and ρ b are the mass and the density of neighboringparticle b . The set of particles Λ a = { b ∈ N | | x a − x b | ≤ κh } contains neighborsof particle a that lie within its defined support domain. The coefficient κ dependson the choice of the kernel, it is equal to for the th order C Wendland kernelfunction [82, 83]) used in this paper. For the sake of clarity, W ( x a − x b , h ) hasbeen denoted W ab . In 2D, this kernel is expressed as follows : W ( q, h ) = 74 πh (cid:16) − q (cid:17) (1 + 2 q ) , with q = | x a − x b | h . (3.9)Several multiphase formulations [84, 85, 86] have been proposed for SPHthroughout the years. In this work, the formalism presented in [42] has beenadopted. The density is directly evaluated through a kernel summation whichgives an exact solution to the continuity equation. It reads : ρ a = m a (cid:88) b ∈ Λ a W ab . (3.10)Discrete gradient and divergence operators in this formalism are given by : ∇ A ( x a ) = (cid:80) b ∈ Λ a (cid:16) A a Θ a + A b Θ b (cid:17) Θ a ∇ a W ab , ∇ · A ( x a ) = (cid:80) b ∈ Λ a (cid:16) A a Θ a + A b Θ b (cid:17) Θ a ∇ a · W ab , (3.11)where Θ a = ρ a / m a . It follows that the full multiphase SPH formulation for aparticle a is : 14 a = m a (cid:80) b ∈ Λ a W ab , D u Dt = − m a (cid:80) b ∈ Λ a (cid:16) p a I +Π a Θ a + p b I +Π b Θ b (cid:17) ∇ a W ab + m a (cid:80) b ∈ Λ a µ a µ b µ a + µ b (cid:16) a + b (cid:17) x ab · ∇ a W ab | x ab | + η u ab + g ,p a = c a ρ a γ a (cid:104)(cid:16) ρ a ρ a (cid:17) γ a − (cid:105) + p , D x a Dt = u a , (3.12)where the viscous term ν ∇ u has been discretized using the inter-particle av-eraged shear stress [87] and η = 0 . h is a safety factor to avoid a divisionby zero. Moreover, the evaluation of normal vectors is performed through thecomputation of the gradient of a color function χ defined for a given particle a and a given phase α as : χ αa = (cid:40) if a ∈ phase α, else . (3.13)The normal vector n αβa of particle a belonging to phase α to the interface αβ is then evaluated by : n αβa = ∇ χ αβa = (cid:88) b ∈ Λ a (cid:32) χ βa Θ a + χ βb Θ b (cid:33) Θ a ∇ a W ab . (3.14)The surface delta function δ αβa is chosen to be equal to | n αβa | and ˜ n αβ = n αβa (cid:14) | n αβa | for use in Eq. (3.7). Corrective terms are commonly used in SPH to remediate intrinsic issuesof this formulation such as particles disorder or micro-mixing at the interface.Three different SPH corrections have been used in this work :1. As suggested in [88], the kernel gradient is enhanced in order to restoreconsistency. For a given particle a , it reads : (cid:103) ∇ W ab = L a ∇ W ab , (3.15)where L a = (cid:16)(cid:80) b ∈ Λ a m a ρ a ∇ W ab ⊗ ( x b − x a ) (cid:17) − . Note that the tilde notationwill be dropped in the rest of paper although the kernel gradient correction willbe always used.2. In order to maintain a good spatial distribution of particles and ensure abetter accuracy, a shifting technique for multiphase flows has been used [89].At the end of every timestep, all particles are shifted by a distance δ r s fromtheir original position. This shifting distance of a particle a is implementedthrough : δ r sa = (cid:26) − D a ∇ C a , if a ∈ light phase − D a (cid:0) ∂C a ∂s s + α n (cid:0) ∂C a ∂n n − β n (cid:1)(cid:1) , if a ∈ heavy phase (3.16)15here C a = (cid:80) b ∈ Λ a m a ρ a W ab is the particle concentration, ∇ C a = (cid:80) b ∈ Λ a m a ρ a ( C b − C a ) ∇ W ab is the particle concentration gradient, D a is the diffusion coefficient, s and n are respectively the tangent and normal vectors to the interface light/heavyphase (with n oriented towards the light phase), β n is a reference concentra-tion gradient (taken equal to its initial value) and α n is the normal diffusionparameter set equal to . . The diffusion coefficient D a is computed as follows: D a = A s | u a | ∆ t, (3.17)where A s is a parameter set to , u a is the velocity of particle a , and ∆ t is thetime step.3. Multiphase SPH can suffer from sub-kernel micro-mixing phenomena as high-lighted in numerous publications [84, 85, 90, 91]. Around the interface, withina distance corresponding to the range of the smoothing length, particles tendto mix. It is because there is no mechanism that guarantees phases immisci-bility in the surface tension’s continuum surface stress model. As suggested bythe previously mentioned authors, we introduce a small repulsive force betweenphases as follows : F corra = ε (cid:88) b ∈ Λ a ,b/ ∈ Ω a (cid:18) a + 1Θ b (cid:19) ∇ a W ab , (3.18)where ε = L ref / h for all simulations as suggested in [92] and where L ref is areference length, typically the diameter of the pipe. The impact of this correctiveforce on the simulation of intermittent flows is evaluated in [93]. Concerning time integration in SPH, different schemes are eligible. Amongthem, one can mention the Runge-Kutta, the Verlet or the Leapfrog schemes.In this work, the Predictor-Corrector Leapfrog scheme was adopted. It is asymplectic integrator which is recommended for SPH because of its conservativenature [81]. Indeed, SPH generally requires very small time steps resulting in ahigh number of iterations. The algorithm proceeds through the following steps.For every particle a ,1. Predictor step : u n = (cid:26) u , if t = 0 , u n − + ∆ t D u Dt n − , if t > .
2. Compute ρ na and p na using the corresponding expressions in Eq. (3.12).3. Evaluate D u Dt n using the momentum equation in Eq. (3.12).4. Corrector step : u n + = (cid:26) u n + ∆ t D u Dt n , if t = 0 , u n − + ∆ t D u Dt n , if t > , x n +1 = x n + ∆ t u n + . ∆ t has to respect the Courant-Friedrichs-Lewy (CFL)criteria to ensure a stable evolution of the system, e.g. ∆ t = min (∆ t visc , ∆ t grav , ∆ t speed , ∆ t st ) , (3.19)where, following [79], we have : ∆ t visc = 0 .
125 min α ∈{ ,...,N phases } h ρ α µ α , ∆ t grav = 0 . (cid:113) h | g | , ∆ t speed = 0 .
25 min α ∈{ ,...,N phases } hc α , ∆ t st = 0 .
25 min α,β ∈{ ,...,N phases } (cid:113) h ρ α πσ αβ . (3.20)A recent article [94] proposes a detailed investigation on maximum admissi-ble time steps for WCSPH. Ω f Ω g Λ g Γ fg Figure 4: Schematic of a ghost particle g (in black) and its associated support domain Λ g (hatched area) intersecting with the fluid domain Ω f (gray particles) and the ghost domain Ω g (white particles) separated by Γ fg . In the subsequent simulations, we used three types of boundary conditions(BC): no-slip wall (superscript w ), periodic and inlet/outlet (superscripts in and out ).For no-slip wall BC, the ghost particle method has been used along with thefollowing prescribed values for the pressure p w , density ρ w and velocity v w fora given ghost particle g : p wg = 1 V ga (cid:88) a ∈ Ω f ∩ Λ g p a m a ρ a W ga , (3.21) ρ wg = 1 V ga (cid:88) a ∈ Ω f ∩ Λ g ρ a m a ρ a W ga , (3.22) v wg = − V ga (cid:88) a ∈ Ω f ∩ Λ g v a m a ρ a W ga , (3.23)with V ga = (cid:80) a ∈ Ω f ∩ Λ g m a ρ a W ga , Ω f the set of fluid particles and Λ g the set ofneighboring particles of ghost particle g . A schematic drawn in Fig. 4 helps tovisualize what is the intersection Ω f ∩ Λ g . Note that, for free-slip wall BC, one17an use the same equations and change the sign of Eq. (3.23) for the directionwhere the free slip is allowed. Additionally, for velocity wall BC, the term v (where v is the prescribed velocity) can be added to Eq. (3.23).For periodic BC, different variants are available in the literature. In essence,particles close to a periodic boundary are allowed to interact with particles nearan associated boundary. Quantities are exchanged both ways and if a particleleaves the domain through one side, it reenters the domain from the other side.For inlet/outlet BC, the method extensively described in [16] (combiningideas from [95, 96]) has been used. The inlet and outlet boundaries are extendedwith a buffer layer of size κh to ensure a full kernel support. At the inlet, the goalis to inject particles with a prescribed velocity profile. At the outlet, particlesneed to leave the domain smoothly while imposing a prescribed pressure profile(or density since they are connected through Eq. (3.3)).On one hand, a particle i in the inlet buffer is moving with a prescribedvelocity profile v p and it carries the following values of pressure p in , density ρ in and velocity v in : p in i = 1 V ia (cid:88) a ∈ Ω f ∩ Λ i p a m a ρ a W ia , (3.24) ρ in i = 1 V ia (cid:88) a ∈ Ω f ∩ Λ i ρ a m a ρ a W ia , (3.25) v in i = 2 v p − V ia (cid:88) a ∈ Ω f ∩ Λ i v a m a ρ a W ia . (3.26)with V ia = (cid:80) a ∈ Ω f ∩ Λ i m a ρ a W ia and Λ i the set of neighboring particles of inletparticle i .On the other hand, at the outlet, a particle o in the buffer is moved accordingto a smoothed convective velocity v conv . For example, if the outlet boundary isvertical and the flow leaves along the x direction, it reads : v out,convo = 1 V (cid:48) oa (cid:88) a ∈ Λ o v a m a ρ a W oa , (3.27)with V (cid:48) oa = (cid:80) a ∈ Λ o m a ρ a W oa the set of neighboring particles of outlet particle o . Note that in Eq. (3.27), the summation is over the full kernel support Λ o including fluid and outlet particles and not only over the intersection Ω f ∩ Λ o .Besides, particle o also carries the following values of pressure p out , density ρ out and velocity v out : 18 out o = 2 p p − V oa (cid:88) a ∈ Ω f ∩ Λ o p a m a ρ a W oa , (3.28) ρ out o = 2 ρ p − V oa (cid:88) a ∈ Ω f ∩ Λ o ρ a m a ρ a W oa , (3.29) v out o,x = 1 V oa (cid:88) a ∈ Ω f ∩ Λ o v a,x m a ρ a W oa , (3.30) v out o,y = − V oa (cid:88) a ∈ Ω f ∩ Λ o v a,y m a ρ a W oa , (3.31)with V oa = (cid:80) a ∈ Ω f ∩ Λ o m a ρ a W oa , p p and ρ p the prescribed pressure and density.Concerning the velocity, null cross velocities (here v y ) are enforced to ensurea divergence free velocity field at the outlet. An evaluation of these boundaryconditions on a collection of test cases can be found in [16].
4. Static bubble tests
In this section, the goal is to validate and compare the implementation ofLBM and SPH surface tension models respectively described in Secs. 2.2 and 3.1.
Square-to-droplet case..
The standard square-to-droplet test case is simulatedand when a steady state is reached, the pressure difference between the exteriorand the interior of the bubble is measured and compared to Laplace’s formula: ∆ P = σR = σ √ πa , (4.1)with ∆ P the pressure difference, σ the surface tension coefficient, R the bubble’sradius and a the lateral dimension of the initial square bubble. Simulations areperformed for three different resolutions : × , × and × nodes/particles. Four different combinations of density and viscosity ratios weretested : ( ρ heavy ρ light , µ heavy µ light ) = (1 , , (5 , , (1000 , and (1 , . The surfacetension coefficient is σ = 1 .
88 N / m for SPH and σ = 0 . l.u. for LBM (where l.u. stands for lattice units). The whole domain is × and the lateralsize of the initial square droplet is a = 0 .
33 m . The time is normalized by t σ = (cid:112) ρa /σ . Note that, following [55], parameter β in Eq.2.28 is adjustedwhen the resolution is increased taking the lowest resolution as a reference i.e.: β = β × (cid:18) ∆ x ∆ x × (cid:19) / , (4.2)with β × = 0 . .Initially, when the density and viscosity ratios are set to one, one can observein Fig. 5 the deformation of an initial square bubble towards a circular bubbleunder the influence of the surface tension. The Lagrangian/Eulerian differencebetween SPH and LBM is magnified in Fig. 5. We clearly see that, in SPH,particles of each phase move to form a circular bubble over time whereas inLBM, nodes are fixed and they switch phase to form the expect circular bubble.19esides, when the circular bubble is stabilized, both methods present residualvelocities around the interface as shown in Fig. 6. However, those spurious cur-rents are much more spread into the domain in SPH compared to LBM wherethey are localized around the interface. Note that, in the LBM color gradi-ent framework, it is possible to significantly reduce the amplitude of spuriouscurrents by choosing a more isotropic gradient operator [51] but, as it involvessecond range neighbors, it is more computationally expensive. it leads to. ForSPH, Hu and Adams’ formulation [42] used in this paper has been reported togenerate stronger spurious currents than other formulations [97]. (a) t/t σ = 0 . (b) t/t σ = 0 . (c) t/t σ = 1 . (d) t/t σ = 0 . (e) t/t σ = 0 . (f) t/t σ = 1 . Figure 5: Evolution of the initial square bubble at selected timesteps. (a,b,c) SPH. (d,e,f)LBM.Figure 6: Normalized velocity field | v || v max | for LBM (left) and SPH (right) at t/t σ = 1 . . In Fig. 7, the pressure profiles at steady state for the different resolutions andthe different density and viscosity ratios considered are shown along with thecorresponding L error plots. First, one can clearly see that LBM is returningincorrect pressure values at the interface. Indeed, LBM presents non-physical20ressure peaks at the bubble’s interface that tend to grow with the number ofnodes. In fact, in the LBM color gradient method, the pressure is not well-defined at the interface. The pressure formula of Eq. (2.13) does not makesense at the interface where fluids are mixed and there is no mixture pressuredefined in the considered framework. Hence, summing the fluid pressures isjust an analytical construction that depends on the density profiles and theinterface width. Thus, if we look at the LBM L error along the whole horizontalcenterline, we do not have mesh convergence since the error is growing at theinterface. However, when we restrict the calculation of the L error inside thebubble (i.e. when . < X < . ), we do obtain a negative slope indicatingmesh convergence.Next, analyzing the impact of the density and viscosity ratios, for the casewhere ( ρ heavy ρ light , µ heavy µ light ) = (1 , in Figs. 7a and 7b, we get approximately thesame order of convergence for both methods ( . and . for LBM andSPH respectively) and the same error levels ( ≤ at the bubble’s center)even though SPH always has a slightly higher error level. Next, when thedensity and viscosity ratios remains moderate (i.e. respectively up to and inFigs. 7c and 7d), both methods are under . error compared to the referencesolution. Additionally, we see that LBM offers a better order of convergencethan SPH. In fact, LBM sees its order of convergence increased (from . to . ) compared to the previous case unlike SPH where it decreases (from . to . ). Moreover, LBM is less accurate than SPH for the two lowestresolutions × and × but performs better for the × casethanks to its higher order of convergence. Overall, although both methodsare returning satisfactory results for this case, we begin to observe a fall inperformance whether it is for the order of convergence (for SPH) or for the errorlevels (for LBM and SPH) because gradients at stake are steeper. Then, for thehigh density ratio case where ( ρ heavy ρ light , µ heavy µ light ) = (1000 , shown in Figs. 7eand 7f, we can see that the order of convergence of SPH remains roughly thesame compared to the first case ( . vs . ). The maximum error levelis ≤ i.e. higher than the first case. This tends to indicate that when thedensity ratio increases, SPH is a quite robust and offers a reasonable accuracyfor the same order of convergence. On the other hand, LBM sees its order ofconvergence strongly affected by the presence of this density ratio ( . vs . ) while maintaining approximately the same error level. Finally, we lookedat one last case, shown in Figs. 7g and 7h, where the density ratio is equal to and the viscosity ratio increased up to . One can immediately note that,for both methods, the pressure profiles are heavily impacted at the interface(oscillations) in particular for the × case. However, when looking at thepressure jump at the center of the bubble, LBM appears very robust to thepresence of such a strong viscosity ratio. Indeed, we see that the error levels areof the same order than those of the first case and that the order of convergenceis even higher ( . vs . ). It shows that refining the lattice strongly helpsto stabilize the pressure field. On the contrary, SPH appears more affected.The error levels are the highest of all four cases considered and the order ofconvergence is inferior to the first case ( . vs . ). Moreover, the errordoes not seem to decrease anymore exponentially with the number of particlesalthough more simulations would be needed to further check that statement.To sum up, for limited density ratios and viscosity ratios, both methods are21ble to reproduce the pressure jump predicted by Eq. (4.1) with a good accuracyand with steep and clean pressure/density profiles. When the the density ratioincreases up to , SPH seems to be more resilient than LBM in the sense thatits order of convergence is not impacted by the presence of such a strong densityratio. LBM seems less robust in the same situation. On the contrary, when theviscosity ratio goes up to , both methods render perturbed pressure profiles.However, it is SPH that appears to have more problem to handle a strongviscosity ratio whereas LBM maintains its performance level. Contact angle case..
In addition, we compared the ability of the previously de-scribed implementations of SPH and LBM to prescribe a contact angle betweena wetting phase, a non-wetting phase and a solid phase. Simulations are donewith × nodes/particles. The density and viscosity ratios are both equalto one. The whole domain is × and the lateral dimensions of the initialrectangular droplet is .
33 m × .
165 m . For SPH, the surface tension coefficientbetween the wetting and non-wetting phase is σ nw = 1 .
88 N / m whereas the onebetween the wetting and the solid phase is set to σ sw = 0 N / m and the onebetween the non-wetting and the solid phase σ sn is adjusted to prescribed thedesired contact angle θ prescribed c using the Young-Laplace equation θ c = σ sw − σ sn σ nw .For LBM, we follow the procedure described in Sec. 2.4. The surface tensioncoefficient is set to σ = 0 . l.u. . At steady state, the observed contact angle θ observed c is measured and reported in Fig. 8. The coefficient of determinationis ≥ . for both methods confirming that they can accurately reproduce aprescribed contact angle. Finally, one can observe the normalized velocity fieldat steady state for θ prescribed c = 150 ◦ . The same comments made before are stillvalid, LBM spurious currents are less spread throughout the domain than inSPH. This is likely due to the Lagrangian nature of SPH where particles haveto rearrange to match the simulated physics.
5. Intermittent two-phase flows in pipes
In the following section, we extend our comparative study to two cases ofintermittent two-phase flows in pipes for different Reynolds numbers. The firstone is periodic and gravity-driven while the second one is generated by a velocityinlet and a pressure boundary condition respectively at the inlet and outlet ofthe pipe.
In this section, we study the establishment of different periodic two-phaseflow patterns under the influence of gravity g starting from a given bubbly flow.Following [15], the initial configuration is composed of . of light phase and . of heavy phase and is described in Fig. 10. All physical properties andsimulation properties are in Tab. 1. The heavy phase viscosity µ l is adjustedas function of the Reynolds number Re = gH ν l . The initial velocity field V is V ( x, y ) = g ν l y ( H − y ) . Viscosity values and dimensionless numbers for eachcase are reported in Tabs. 2 and 3. No-slip boundary conditions are appliedto the walls. The simulations is done with 50000 nodes/particles for t = 30 s .Four different Reynolds numbers were tested : Re = 10 , , and . The22 . . . . . . . X ∆ P / ∆ P a n a l y t i c a l × -SPH × -SPH × -SPH × -LBM × -LBM × -LBMAnalytical (a) ( ρ heavy ρ light , µ heavy µ light ) = (1 , − − − − Number of nodes/particles L E rr o r SPH fullcase ( S = − . )SPH bubble only ( S = − . )LBM fullcase ( S = 0 . )LBM bubble only ( S = − . ) (b) ( ρ heavy ρ light , µ heavy µ light ) = (1 , . . . . . . . X ∆ P / ∆ P a n a l y t i c a l × -SPH × -SPH × -SPH × -LBM × -LBM × -LBMAnalytical (c) ( ρ heavy ρ light , µ heavy µ light ) = (5 , − − − − Number of nodes/particles L E rr o r SPH fullcase ( S = − . )SPH bubble only ( S = − . )LBM fullcase ( S = 0 . )LBM bubble only ( S = − . ) (d) ( ρ heavy ρ light , µ heavy µ light ) = (5 , . . . . . . . X ∆ P / ∆ P a n a l y t i c a l × -SPH × -SPH × -SPH × -LBM × -LBM × -LBMAnalytical (e) ( ρ heavy ρ light , µ heavy µ light ) = (1000 , − − − − Number of nodes/particles L E rr o r SPH fullcase ( S = − . )SPH bubble only ( S = − . )LBM fullcase ( S = 0 . )LBM bubble only ( S = − . ) (f) ( ρ heavy ρ light , µ heavy µ light ) = (1000 , . . . . . . . X ∆ P / ∆ P a n a l y t i c a l × -SPH × -SPH × -SPH × -LBM × -LBM × -LBMAnalytical (g) ( ρ heavy ρ light , µ heavy µ light ) = (1 , − − − − Number of nodes/particles L E rr o r SPH fullcase ( S = − . )SPH bubble only ( S = − . )LBM fullcase ( S = 0 . )LBM bubble only ( S = − . ) (h) ( ρ heavy ρ light , µ heavy µ light ) = (1 , Figure 7: (a,c,e,g) Pressure profiles at steady state for different resolutions and different den-sity and viscosity ratios. (b,d,f,h) Log-log L error plots as function of resolution (superposedwith linear regressions of slope S ).
30 60 90 120 150 1800306090120150180 θ prescribed c ( ◦ ) θ o b s e r v e d c ( ◦ ) SPHLBM
Figure 8: Comparison between prescribed and observed contact angles for both SPH and LBM(density and viscosity ratios are both equal to one).Figure 9: Normalized velocity field for LBM (left) and SPH (right) for θ prescribed c = 150 ◦ . L = 6 HH L/ H/ H/ L/ (periodic) g Figure 10: Initial configuration sketch. The large bubbles radius is R = 0 . H and the smallbubbles’ radius is R = R/ . The pipe height is H = 0 . . Property Light Phase Heavy Phase UnitsDensity ( ρ ) . . / m Viscosity ( µ ) µ g = µ l / µ l Pa . s Contact Angle ( θ c ) ◦ Surface Tension ( σ nw ) . × − N / m Gravity ( g x ) . / s Space step ( ∆ x ) . × − m Domain size ( L x × L y ) . × . (a) SPH Light Phase Heavy Phase Units . . l.u.µ g = µ l / µ l l.u. ◦ . × − l.u. . × − l.u. . × − l.u. × l.u. (b) LBMTable 1: Simulation parameters. ase µ g ( Pa . s ) µ l ( Pa . s ) Re = 10 ) . × − . × − Re = 50 ) . × − . × − Re = 100 ) . × − . × − Re = 500 ) . × − . × − (a) SPH µ g ( l.u. ) µ l ( l.u. ) τ g τ l . × − . × − .
014 0 . . × − . × − .
730 0 . . × − . × − .
663 0 . . × − . × − .
572 0 . (b) LBMTable 2: Viscosity values for each case. Case Re = gL y ν l La = σρ l µ l Bo = ∆ ρgL y σ
10 82 0 .
50 408 0 .
100 816 0 .
500 4077 0 . Table 3: Reynolds, Laplace and Bond numbers for each case. Re = 10 Re = 50 Re = 100 Re = 500 Figure 11: Steady state for different Re . SPH and LBM results are on the left and rightcolumns respectively. Re numbers considered. For ≤ Re ≤ , we obtain abubbly flow composed of three different Taylor bubbles whereas for Re = 500 ,we have an annular flow where the heavy phase is in contact with the pipeand the light phase travels in the middle. For Re = 10 , we observe that SPHpresent a bubbly flow where one bubble is clearly smaller than the two others.It is not the case in LBM where all bubbles are identical within each case.Besides, for this case, we have heavy phase droplets than are captured insidelight phase bubbles. Note that these small bubbles are to be absorbed by themain flow if the simulation lasted longer because they are moving slower thantheir environment. For Re = 50 and Re = 100 , we obtain in all cases thesame pattern made of three identical Taylor bubbles. Moreover, as shown inFig. 12, the bubbles’ shapes between SPH and LBM for ≤ Re ≤ are verysimilar. Finally, as Re grows, the Taylor bubbles are getting slightly shorterand higher in size. For Re = 500 , we again see that LBM offers a perfectlysymmetric annular pattern. On the contrary, the bottom heavy phase layer inSPH is thicker than the top one. In general, LBM provides more symmetricresults than SPH because of its Eulerian nature. Re = 10 Re = 50 Re = 100 Figure 12: Superposition of bubbles’ shapes obtained with SPH (blue) and LBM (red) fordifferent Reynolds numbers.
In Fig. 13, we can see that for Re = 10 , the pressure fields is dominatedby the captured droplets of heavy phase inside the bubbles. For Re = 50 and , the pressure field reaches a maximum for SPH inside the bubbles whereasfor LBM it is at the interface. Nevertheless, as predicted by Laplace’s law, thepressure is higher inside the light phase’s bubbles than in the heavy phase bulk.When looking at the velocity fields in Fig. 14, we see that they are also verysimilar. The same patterns surrounding the bubbles can be observed. For theannular case where Re = 500 , it is possible to see that the no-slip condition onthe walls affects the flow more strongly in LBM than in SPH which results in aflatter velocity profile for the latter.For Re = 50 and Re = 100 , where the SPH and LBM patterns are theclosest, we compared the density, velocity and pressure fields along the centerlineon Figs. 15 and 16. Note that because the bubbles do not have the exact sameposition, we have shifted the LBM profiles from a fixed distance to be able tosuperpose the profiles. On the density plots of Figs. 15a and 16a, the differentdensity treatment in both methods clearly appears. In LBM, the density issmoothed at the interface whereas in SPH, thanks to its Lagrangian nature,26 e = 10 Re = 50 Re = 100 Re = 500 Figure 13: Normalized pressure fields at steady state for different Re . SPH and LBM resultsare on the left and right columns respectively. Pressure is normalized by the maximum pressureinside the bubbles. Re = 10 Re = 50 Re = 100 Re = 500 Figure 14: Normalized velocity fields at steady state for different Re . SPH and LBM resultsare on the left and right columns respectively. faster).At each interface, the velocity reaches a local minimum. The only differencesbetween both profiles is that in certain areas, the SPH velocity peaks have asmaller amplitude than in LBM. For example, the bubble velocity is the samefor all three bubbles in LBM for Re = 50 whereas for SPH the last bubbletravels about faster than the other ones. One last comment is that in LBMat the interface, the velocity field present non-physical oscillations due to thefluids mixing at the interface. It is not the case in SPH because the pressurefield does not suffer from pressure overshoots and fluids are not mixed at theinterface. . . . . X − X min ) /X max ρ LBM SPH (a) Density . . . . . . . X − X min ) /X max ( p − p m i n ) / ( p | bubb l e m a x − p m i n ) LBM SPH (b) Normalized Pressure . . . . . . . .
81 ( X − X min ) /X max ( V x − V x m i n ) / ( V x | bubb l e s m a x − V x m i n ) LBM SPH (c) Normalized VelocityFigure 15: Case Re = 50 . (a) Superposed densities. (b) Superposed normalized pressures. (c)Superposed normalized velocities. . . . . X − X min ) /X max ρ LBM SPH (a) Density . . . . . . . X − X min ) /X max ( p − p m i n ) / ( p | bubb l e m a x − p m i n ) LBM SPH (b) Normalized Pressure . . . . . . . .
81 ( X − X min ) /X max ( V x − V x m i n ) / ( V x | bubb l e s m a x − V x m i n ) LBM SPH (c) Normalized VelocityFigure 16: Case Re = 100 . (a) Superposed densities. (b) Superposed normalized pressures.(c) Superposed normalized velocities.
28o conclude, we can add that SPH and LBM are both well capable of simu-lating the transition from a given bubbly flow to a slug flow composed of Taylorbubbles for Re ≤ . To further assess their relative performance, an extendedcomparison with other methods or with experimental data would be of greatinterest. Note that we have limited our study to Re ≤ because, for highervelocities and/or smaller viscosities, we lie outside LBM stability region whetherbecause the low Mach rule is violated or because the relaxation time is too closefrom . . In this section, we study the ability of both methods to simulate a predictedintermittent flow regime. We consider an horizontal pipe of diameter D = 1 m and length L = 10 D . The light phase and heavy phase are denoted with a g and l subscript respectively. The flow enters from the inlet (left) and is assumed tobe stratified with given volume fractions for each phase α g = 0 . and α l = 0 . .All the physical properties are summarized in Tab. 4. Using these properties, itis possible to plot the flow regime map, see Fig. 18 , and to pick an area to beinvestigated in the intermittent region. In this area, we adjust the viscosity tochoose two cases that correspond to Re = 125 and Re = 312 . . Viscosity valuesand dimensionless numbers for each case are reported in Tabs. 5 and 6. Bothcases were simulated in 2D with nodes/particles. The simulation time was
30 s . At the inlet, each phase is injected with a constant velocity correspondingto its superficial velocity u sg,l = α g,l u g,l . At the outlet, a constant pressureequal to the initial pressure p is prescribed. Free slip boundary conditions areapplied on the top and bottom walls (we were not able to generate a slug flow inLBM with no-slip boundary conditions under the same conditions). The initialsetup is presented in Fig. 17. The phases distributions for each case are shownin Fig. 19 along with the associated plots showing the volume fraction evolutionover time, the average pressure drop evolution over time, the heavy/light phasevelocity evolution over time respectively in Fig. 20, Fig. 22 and Fig. 23. L = 10m D = 1m u g u l p Light Fluid ( α g )Heavy Fluid ( α l ) Figure 17: Initial configuration sketch.
In Fig. 19, one can see snapshots of the phases’ distribution for both methodsat selected timesteps. It is clear that both methods are able to generate a slugflow as predicted by the flow map of Fig. 18. However, LBM produces a muchmore regular intermittent flow pattern with a higher slug frequency than SPH.This can also be seen in Figs. 20 and 21. For example, for Re = 312 . , theLBM slug frequency at the outlet is approximately .
85 Hz whereas for SPH itis close to .
52 Hz . The slug frequency seems to remain roughly stable or toslightly increase (about +1 Hz for LBM and about +0 .
16 Hz for SPH for the In order to plot the map, one has to compute the Lockhart-Martelli [98] parameter whichdepends on n , m , C g and C l . In this study, we used n = m = 2 and C g = C l = 0 . roperty Light Phase Heavy Phase UnitsDensity ( ρ ) . . / m Viscosity ( µ ) µ g = µ l / µ l Pa . s Sound speed ( c s ) .
73 68 .
75 m / s Surface Tension ( σ nw ) .
01 N / m Contact Angle ( θ c ) ◦ Gravity ( g z ) .
556 m / s Space step ( ∆ x ) .
02 m
Time step ( ∆ t ) . × − s Domain size ( L x × L y ) × Inlet velocity ( u ) . .
25 m / s (a) SPH Light Phase Heavy Phase Units . . l.u.µ g = µ l / µ l l.u. . l.u. . × − l.u. ◦ . × − l.u. . × − l.u. . × − l.u. × l.u. . × − . × − l.u. (b) LBMTable 4: Simulation parameters. − − − − Superficial gas velocity ( m / s ) Sup e r fi c i a lli q u i d v e l o c i t y ( m / s ) IntermittentMistDispersedStratifiedConsidered area
Figure 18: Flow regime maps
Case µ g ( Pa . s ) µ l ( Pa . s ) Re = 125 ) . × − . × − Re = 312 . ) × − × − (a) SPH µ g ( l.u. ) µ l ( l.u. ) τ g τ l . × − . × − .
725 0 . . × − . × − .
59 0 . (b) LBMTable 5: Viscosity values for each case. Case Re = L y u l ν l La = σρ l µ l Bo = ∆ ρgL y σ
125 0 . . Table 6: Reynolds, Laplace and Bond numbers for each case. a) SPH - Re = 125 (b) LBM - Re = 125 (c) SPH - Re = 312 . (d) LBM - Re = 312 . Figure 19: (a, b) From top to bottom : snapshots at t = 4 . , . ,
25 s . (c, d) From top tobottom : snapshots at t = 4 . , . , . . highest peak) when Re changes from to . although it is less obviousin SPH periodograms due to the noise and composition of the signal. It isexpected that the slug frequency increases when Re increases but we could notraise Re higher without making LBM simulations unstable ( τ → or M a → ). In addition, SPH periodograms are noisier with to major frequencycomponents unlike LBM where one frequency clearly emerges. It indicates thatSPH volume fraction signals are less regular and are composed of signals withdifferent frequencies. Moreover, we can observe in both methods that the pointwhere the first slug appears is in general closer from the pipe entry when Re issmaller. This is expected since when velocities are smaller at the entry, the firstslug tends to form earlier in the pipe.In Fig. 22, the raw and average pressure drops evolution over time for all thedifferent cases studied are shown. One can notice that in all cases, the pressuredrops are ≤ on average (even though strong oscillations are observed). Itmeans that the average pressure along the pipe’s height is higher at the outletthan at the inlet. This phenomenon was already seen in [16] and may be due toseveral factors : the quality of the boundary conditions, the recirculation areasat the entry that tend to lower the pressure and the averaging area chosen (i.e. . (from the inlet) × D ). Nevertheless, we observed that both methods arereturning the same average pressure drop level around ≈ − . p . SPH plotspresent much stronger oscillations than LBM ones (up to times the meanvalue for SPH against times the mean value for LBM). This is a known issuein weakly compressible SPH where the particles spatial distribution is directlylinked to the density/pressure.In Fig. 23, the evolution over time of the heavy phase and light phase veloc-ities within the whole pipe are shown for all cases considered. After a transientperiod of ≈ , the light and heavy phase velocities stabilize around a fixedvalue when a periodic state is reached. The main difference between SPH andLBM on this aspect is that the final velocity values are higher in LBM than31 . . . Time ( s ) V o l u m e F r a c t i o n Heavy Phase Light Phase (a) SPH - Re = 125 . . . Time ( s ) V o l u m e F r a c t i o n Heavy Phase Light Phase (b) LBM - Re = 125 . . . Time ( s ) V o l u m e F r a c t i o n Heavy Phase Light Phase (c) SPH - Re = 312 . . . . Time ( s ) V o l u m e F r a c t i o n Heavy Phase Light Phase (d) LBM - Re = 312 . Figure 20: Evolution of the volume fractions at the outlet over time.
Frequency f ( Hz ) | F ( f ) | (a) SPH - Re = 125 Frequency f ( Hz ) | F ( f ) | (b) LBM - Re = 125 Frequency f ( Hz ) | F ( f ) | (c) SPH - Re = 312 . Frequency f ( Hz ) | F ( f ) | (d) LBM - Re = 312 . Figure 21: Periodograms of the heavy phase volume fraction time series (in blue in Fig. 20).The mean has been removed from the signal before performing the Fourier transform to removethe frequency component. − − · − Time (s) A v e r ag e P r e ss u r e D r o p ( P a ) Raw Average (a) SPH - Re = 125 − − · − Time (s) N o r m a li ze d P r e ss u r e D r o p Raw Average (b) LBM - Re = 125 − − · − Time (s) N o r m a li ze d P r e ss u r e D r o p Raw Average (c) SPH - Re = 312 . − − · − Time (s) N o r m a li ze d P r e ss u r e D r o p Raw Average (d) LBM - Re = 312 . Figure 22: Evolution of the normalized pressure drops over time.
Time (s) N o r m a li ze d V e l o c i t y Heavy phase (Raw) Heavy Phase (Average)Light Phase (Raw) Light Phase (Average) (a) SPH - Re = 125 Time (s) N o r m a li ze d V e l o c i t y Heavy Phase (Raw) Heavy phase (Average)Light Phase (Raw) Light Phase (Average) (b) LBM - Re = 125 Time (s) N o r m a li ze d V e l o c i t y Heavy Phase (Raw) Heavy Phase (Average)Light Phase (Raw) Light Phase (Average) (c) SPH - Re = 312 . Time (s) N o r m a li ze d V e l o c i t y Heavy Phase (Raw) Heavy Phase (Average)Light Phase (Raw) Light Phase (Average) (d) LBM - Re = 312 . Figure 23: Evolution of the normalized velocities over time. v in higherthan in SPH. In our opinion, this is due to the wall boundary conditions thatare handled in a very different way in both methods (full bounce back approachin LBM and interpolation-based approach in SPH) and to the wetting bound-ary conditions that is also handled differently. Those can strongly affect theflow, especially in dynamic cases like the ones considered in this section. Froma general point of view, one can add that the oscillations observed in Figs. 22and 23 are gaining amplitude when Re increases in SPH whereas it is not thecase in LBM for which they tend to reduce.In a nutshell, as predicted by Taitel and Dukler’s flow map, both methods arecapable to generate a slug flow pattern starting from the exact same simulationsetup. Nevertheless, unlike all previous test cases for which the results wereglobally similar, we obtain significantly different flow patterns. Indeed, SPHproduces bigger bubbles and in a more irregular way compared to LBM. Wetend to believe that it is due to boundary conditions. In addition, we could notextend our study to higher Re number because LBM was not stable anymore.Although, we tend to think that LBM is entitled to propose the best solution,in particular because of its superior accuracy, its narrow stability range may bea serious drawback to simulate two-phase flows in pipes at realistic Re numbers.One last point that has not been addressed yet is the computational efficiencyof both methods. On that aspect and despite very important progresses madein the SPH community, LBM remains superior in terms of speed, thanks to itslocal lattice-based algorithm that does not require a nearest neighbor searchat every time step. For this paper, we have implemented both methods inthe same framework in Fortran combined with an OpenMP library to handlemulti-threading. For a simulation involving nodes/particles, LBM wasabout times faster than SPH on a laptop equipped with a . IntelCore i7 processor with cores and Gb of RAM. This number is given as anindicative number and should be taken with caution since the codes were notfully optimized.
6. Conclusions
In this paper, we propose a 2D multiphase comparison of two particle meth-ods, SPH and LBM, that are very different by nature. Despite these differences,both methods have been extensively used to model multiphase flows with suc-cess. We have chosen a multiphase formulation for each method among thoseavailable in the literature : the continuum surface force technique for SPH andthe color gradient method for LBM. Then, we have simulated a collection ofstatic bubble tests with different density and viscosity ratios with both meth-ods. Finally, we have prolonged the comparison to more realistic cases, i.e. tothe generation of slug flows in pipes. To this end, we have extended LBM Zou-He boundary conditions to be able to handle velocity inlet and pressure outletconditions in a multiphase context.From a general point of view, we have confirmed that LBM offers a betterorder of convergence and a better accuracy than SPH although it suffers from amore narrow stability range than SPH. In many situations for which the Machnumber is too high or the viscosity is too low, LBM will be unstable contraryto SPH which is only controlled by the CFL condition. Note that research onthe extension of LBM to high Mach numbers is very active. In addition, SPH34ends to generate pressure fields that are noisier than with LBM because of theLagrangian behavior of particles whose position is directly linked to pressurethrough the density evaluation. This problem has been the subject of manyinvestigations and several treatments are available. Moreover, LBM is morecomputationally efficient than SPH by construction.On the multiphase aspects, both methods are very capable to simulate avariety of dynamic incompressible multiphase flow problems with good precisionin the case of moderate density and viscosity ratios. However, we have notedthat, on one hand, SPH appears more robust for high density ratios than LBMand, on the other hand, SPH has more trouble to handle high viscosity ratiosthan LBM. Another difference is that, at the interface, the fluids are mixedresulting in a diffuse interface whereas with SPH, particles are affected to onephase or the other without any mixing. More specifically, both methods havebeen able to simulate slug flows where expected but, due to boundary conditions,there might be differences in slug frequency and/or slug sizes.To conclude, according to the results presented in this paper, our recommen-dation would be to use LBM when stability is not an issue and SPH otherwise.In future works, we would like to extend our comparison to other methods,experimental measurements and commercial software. This would help to char-acterize more precisely which method is best adapted to a given situation.
Acknowledgements
The authors would like to thank Total for its scientific and financial support.
ReferencesReferences [1] J Fabre and A Line. Modeling of two-phase slug flow.
Annual Review ofFluid Mechanics , 24(1):21–46, jan 1992.[2] Min Lu.
Experimental and computational study of two-phase slug flow .PhD thesis, Imperial College London, 2015.[3] Simon Pedersen, Petar Durdevic, and Zhenyu Yang. Challenges in slugmodeling and control for offshore oil and gas productions: A review study.
International Journal of Multiphase Flow , 88:270 – 284, 2017.[4] Yehuda Taitel, Dvora Bornea, and A. E. Dukler. Modelling flow patterntransitions for steady upward gas-liquid flow in vertical tubes.
AIChEJournal , 26(3):345–354, 1980.[5] M. Viggiani, O. Mariani, V. Battarra, A. Annunziato, and U. Bollettini.A model to verify the onset of severe slugging. In
PSIG Annual Meeting ,Toronto, Ontario, 1988. Pipeline Simulation Interest Group.[6] C. Sarica and O. Shoham. A simplified transient model for pipeline-risersystems.
Chemical Engineering Science , 46(9):2167 – 2179, 1991.[7] Kjell H. Bendiksen, Dag Maines, Randi Moe, and Sven Nuland. Thedynamic two-fluid model olga: Theory and application.
SPE , 1991.358] Kongsberg.
LedaFlow - The new multiphase simulator: User guide . Kongs-berg, 2014.[9] R. Belt, E. Duret, D. Larrey, B. Djoric, and S. Kalali. Comparison of com-mercial multiphase flow simulators with experimental and field databases.In ,Cannes, France, 2011. BHR Group.[10] T. Taha and Z.F. Cui. Hydrodynamics of slug flow inside capillaries.
Chemical Engineering Science , 59(6):1181–1190, mar 2004.[11] Zahid I Al-Hashimy, Hussain H Al-Kayiem, Rune W Time, and Zena KKadhim. Numerical characterisation of slug flow in horizontal air/waterpipe flow.
International Journal of Computational Methods and Experi-mental Measurements , 4(2):114–130, 2016.[12] Koji Fukagata, Nobuhide Kasagi, Poychat Ua-arayaporn, and TakehiroHimeno. Numerical simulation of gas–liquid two-phase flow and convectiveheat transfer in a micro tube.
International Journal of Heat and FluidFlow , 28(1):72–82, feb 2007.[13] Enrique Lizarraga-García.
A study of Taylor bubbles in vertical and in-clined slug flow using multiphase CFD with level set . PhD thesis, Mas-sachusetts Institute of Technology, 2016.[14] Zhao Yu, Orin Hemminger, and Liang-Shih Fan. Experiment and lat-tice boltzmann simulation of two-phase gas–liquid flows in microchannels.
Chemical Engineering Science , 62(24):7172–7183, dec 2007.[15] Jean-Pierre Minier. Simulation of two-phase flow patterns with a newapproach based on smoothed particle hydrodynamics. NUGENIA ProjectSPH-2PHASEFLOW Presentation Slides, 2016.[16] Thomas Douillet-Grellier, Florian De Vuyst, Henri Calandra, and PhilippeRicoux. Simulations of intermittent two-phase flows in pipes usingsmoothed particle hydrodynamics.
Computers & Fluids , 177:101–122, nov2018.[17] Qunwu He and Nobuhide Kasagi. Phase-field simulation of small capillary-number two-phase flow in a microtube.
Fluid Dynamics Research , 40(7-8):497–509, jul 2008.[18] Fangfang Xie, Xiaoning Zheng, Michael S. Triantafyllou, Yiannis Con-stantinides, Yao Zheng, and George Em Karniadakis. Direct numericalsimulations of two-phase flow in an inclined pipe.
Journal of Fluid Me-chanics , 825:189?207, 2017.[19] Haibo Huang, Michael C. Sukop, and Xi-Yun Lu.
Multiphase LatticeBoltzmann Methods: Theory and Application . John Wiley & Sons, Ltd,jul 2015.[20] Zhi-Bin Wang, Rong Chen, Hong Wang, Qiang Liao, Xun Zhu, and Shu-Zhe Li. An overview of smoothed particle hydrodynamics for simulatingmultiphase flow.
Applied Mathematical Modelling , 40(23-24):9625–9655,dec 2016. 3621] Q. Li, K.H. Luo, Q.J. Kang, Y.L. He, Q. Chen, and Q. Liu. Latticeboltzmann methods for multiphase flow and phase-change heat transfer.
Progress in Energy and Combustion Science , 52:62–105, feb 2016.[22] Leon B Lucy. A numerical approach to the testing of the fission hypothesis.
The astronomical journal , 82:1013–1024, 1977.[23] R.A. Gingold and J.J. Monaghan. Smoothed particle hydrodynamics-theory and application to non-spherical stars.
Monthly Notices of theRoyal Astronomical Society , 181:375–389, 1977.[24] Volker Springel. Smoothed particle hydrodynamics in astrophysics.
An-nual Review of Astronomy and Astrophysics , 48(1):391–430, aug 2010.[25] Damien Violeau and Benedict D Rogers. Smoothed particle hydrody-namics (sph) for free-surface flows: past, present and future.
Journal ofHydraulic Research , 54(1):1–26, 2016.[26] Larry D. Libersky and A. G. Petschek. Advances in the free-lagrangemethod including contributions on adaptive gridding and the smooth par-ticle hydrodynamics method: Proceedings of the next free-lagrange con-ference held at jackson lake lodge, moran, wy, usa 3–7 june 1990. pages248–257, Berlin, Heidelberg, 1991. Springer Berlin Heidelberg.[27] H.H. Bui, R. Fukagawa, K. Sako, and S. Ohno. Lagrangian meshfree parti-cles method (SPH) for large deformation and failure flows of geomaterialusing elastic–plastic soil constitutive model.
International Journal forNumerical and Analytical Methods in Geomechanics , 32(12):1537–1570,2008.[28] Thomas Douillet-Grellier, Ranjan Pramanik, Kai Pan, Aziz Albaiz,Bruce David Jones, Hamid Pourpak, and John Richard Williams. Mesh-free numerical simulation of pressure-driven fractures in brittle rocks.
SPEHydraulic Fracturing Technology Conference , 2016.[29] Thomas Douillet-Grellier, Bruce D. Jones, Ranjan Pramanik, Kai Pan,Abdulaziz Albaiz, and John R. Williams. Mixed-mode fracture model-ing with smoothed particle hydrodynamics.
Computers and Geotechnics ,79:73 – 85, 2016.[30] Markus Ihmsen, Jens Orthmann, Barbara Solenthaler, Andreas Kolb, andMatthias Teschner. Sph fluids in computer graphics. The EurographicsAssociation, 2014.[31] J.J. Monaghan. Smoothed particle hydrodynamics and its diverse appli-cations.
Annual Review of Fluid Mechanics , 44(1):323–346, 2012.[32] Daniel J. Price. Smoothed particle hydrodynamics and magnetohydrody-namics.
Journal of Computational Physics , 231(3):759 – 794, 2012. Spe-cial Issue: Computational Plasma PhysicsSpecial Issue: ComputationalPlasma Physics.[33] MS Shadloo, G Oger, and D Le Touzé. Smoothed particle hydrodynam-ics method for fluid flows, towards industrial applications: Motivations,current state, and challenges.
Computers & Fluids , 136:11–34, 2016.3734] U. Frisch, B. Hasslacher, and Y. Pomeau. Lattice-gas automata for thenavier-stokes equation.
Phys. Rev. Lett. , 56:1505–1508, Apr 1986.[35] Erlend Magnus Viggen.
The Lattice Boltzmann Method with Applicationsin Acoustics . PhD thesis, Norwegian University of Science and Technology,2009.[36] Jens Harting, Jonathan Chin, Maddalena Venturoli, and Peter V Coveney.Large-scale lattice boltzmann simulations of complex fluids: advancesthrough the advent of computational grids.
Philosophical Transactionsof the Royal Society of London A: Mathematical, Physical and Engineer-ing Sciences , 363(1833):1895–1915, 2005.[37] R. Benzi, S. Succi, and M. Vergassola. The lattice boltzmann equation:theory and applications.
Physics Reports , 222(3):145 – 197, 1992.[38] Antonio Cancelliere, Celeste Chang, Enrico Foti, Daniel H. Rothman,and Sauro Succi. The permeability of a random medium: Comparison ofsimulation with theory.
Physics of Fluids A , 2(12):2085–2088, 1990.[39] Bruno Ferreol and Daniel H. Rothman. Lattice-boltzmann simulationsof flow through fontainebleau sandstone.
Transport in Porous Media ,20(1):3–20, 1995.[40] Benjamin Ahrenholz, Jonas Tilke, and Manfred Krafczyk. Lattice-boltzmann simulations in reconstructed parametrized porous media.
In-ternational Journal of Computational Fluid Dynamics , 20(6):369–377,2006.[41] Zhaoli Guo and T. S. Zhao. Lattice boltzmann model for incompressibleflows through porous media.
Phys. Rev. E , 66:036304, Sep 2002.[42] X.Y. Hu and N.A. Adams. A multi-phase SPH method for macroscopicand mesoscopic flows.
Journal of Computational Physics , 213(2):844 –861, 2006.[43] T Reis and T N Phillips. Lattice boltzmann model for simulating immisci-ble two-phase flows.
Journal of Physics A: Mathematical and Theoretical ,40(14):4033, 2007.[44] Xiaowen Shan and Hudong Chen. Lattice boltzmann model for simu-lating flows with multiple phases and components.
Physical Review E ,47(3):1815–1819, mar 1993.[45] Michael R. Swift, W. R. Osborn, and J. M. Yeomans. Lattice boltzmannsimulation of nonideal fluids.
Physical Review Letters , 75(5):830–833, jul1995.[46] Xiaoyi He, Shiyi Chen, and Raoyang Zhang. A lattice boltzmann schemefor incompressible multiphase flow and its application in simulation ofrayleigh–taylor instability.
Journal of Computational Physics , 152(2):642–663, jul 1999. 3847] Andrew K. Gunstensen, Daniel H. Rothman, Stéphane Zaleski, and Gi-anluigi Zanetti. Lattice boltzmann model of immiscible fluids.
Phys. Rev.A , 43:4320–4327, Apr 1991.[48] Haihu Liu, Qinjun Kang, Christopher R. Leonardi, Sebastian Schmi-eschek, Ariel Narváez, Bruce D. Jones, John R. Williams, Albert J. Val-occhi, and Jens Harting. Multiphase lattice boltzmann simulations forporous media applications.
Computational Geosciences , 20(4):777–805,dec 2015.[49] Sébastien Leclaire, Andrea Parmigiani, Bastien Chopard, and Jonas Latt.Three-dimensional lattice boltzmann method benchmarks between color-gradient and pseudo-potential immiscible multi-component models.
In-ternational Journal of Modern Physics C , 28(07):1750085, jul 2017.[50]
Computational Methods for Multiphase Flow . Cambridge University Press,2009.[51] Sébastien Leclaire, Marcelo Reggio, and Jean-Yves Trépanier. Isotropiccolor gradient for simulating very high-density ratios with a two-phaseflow lattice boltzmann model.
Computers & Fluids , 48(1):98–112, sep2011.[52] Sébastien Leclaire, Marcelo Reggio, and Jean-Yves Trépanier. Numericalevaluation of two recoloring operators for an immiscible two-phase flowlattice boltzmann model.
Applied Mathematical Modelling , 36(5):2237–2252, may 2012.[53] Sébastien Leclaire, Nicolas Pellerin, Marcelo Reggio, and Jean-Yves Tré-panier. Enhanced equilibrium distribution functions for simulating im-miscible multiphase flows with variable density ratios in a class of latticeboltzmann models.
International Journal of Multiphase Flow , 57:159–168,dec 2013.[54] S Leclaire, N Pellerin, M Reggio, and J-Y Trépanier. Unsteady im-miscible multiphase flow validation of a multiple-relaxation-time latticeboltzmann method.
Journal of Physics A: Mathematical and Theoretical ,47(10):105501, feb 2014.[55] Sébastien Leclaire, Nicolas Pellerin, Marcelo Reggio, and Jean-Yves Tré-panier. An approach to control the spurious currents in a multiphaselattice boltzmann method and to improve the implementation of ini-tial condition.
International Journal for Numerical Methods in Fluids ,77(12):732–746, jan 2015.[56] Sébastien Leclaire, Kamilia Abahri, Rafik Belarbi, and Rachid Bennacer.Modeling of static contact angles with curved boundaries using a multi-phase lattice boltzmann method with variable density and viscosity ratios.
International Journal for Numerical Methods in Fluids , 82(8):451–470, oct2016.[57] Zhiyuan Xu, Haihu Liu, and Albert J. Valocchi. Lattice boltzmann sim-ulation of immiscible two-phase flow with capillary valve effect in porousmedia.
Water Resources Research , 53(5):3770–3790, 2017.3958] Yan Ba, Haihu Liu, Qing Li, Qinjun Kang, and Jinju Sun. Multiple-relaxation-time color-gradient lattice boltzmann model for simulating two-phase flows with high density ratio.
Phys. Rev. E , 94:023310, Aug 2016.[59] J. Tolke, M. Krafczyk, M. Schulz, and E. Rank. Lattice boltzmann sim-ulations of binary fluid flow through porous media.
Philosophical Trans-actions of the Royal Society A: Mathematical, Physical and EngineeringSciences , 360(1792):535–545, mar 2002.[60] Daryl Grunau, Shiyi Chen, and Kenneth Eggert. A lattice boltzmannmodel for multiphase fluid flows.
Physics of Fluids A: Fluid Dynamics ,5(10):2557–2562, oct 1993.[61] Pierre Lallemand and Li-Shi Luo. Theory of the lattice boltzmann method:Dispersion, dissipation, isotropy, galilean invariance, and stability.
Phys.Rev. E , 61:6546–6562, Jun 2000.[62] Haihu Liu, Albert J. Valocchi, and Qinjun Kang. Three-dimensional lat-tice boltzmann model for immiscible two-phase flow simulations.
Phys.Rev. E , 85:046309, Apr 2012.[63] Hatro Huang, Jun-Jie Huang, Xi-Yun LU, and Michael C. Sukop. Onsimulations of high-density ratio flows using color-gradient multiphaselattice Boltzmann models.
International Journal of Modern Physics C ,24(04):1350021, apr 2013.[64] D. J. Holdych, D. Rovas, J. G. Georgiadis, and R. O. Buckius. An im-proved hydrodynamics formulation for multiphase flow lattice-boltzmannmodels.
International Journal of Modern Physics C , 9(8):1393–1404, 121998.[65] Xiaowen Shan, Xue Feng Yuan, and Hudong Chen. Kinetic theory rep-resentation of hydrodynamics: A way beyond the navier-stokes equation.
Journal of Fluid Mechanics , 550:413–441, 3 2006.[66] Q. Li, K. H. Luo, Y. L. He, Y. J. Gao, and W. Q. Tao. Coupling latticeboltzmann model for simulation of thermal flows on standard lattices.
Phys. Rev. E , 85:016710, Jan 2012.[67] I. Halliday, S. P. Thompson, and C. M. Care. Macroscopic surface tensionin a lattice bhatnagar-gross-krook model of two immiscible fluids.
PhysicalReview E , 57(1):514–523, jan 1998.[68] Sébastien Leclaire, Maud El-Hachem, Jean-Yves Trépanier, and MarceloReggio. High order spatial generalization of 2d and 3d isotropic discretegradient operators with fast evaluation on GPUs.
Journal of ScientificComputing , 59(3):545–573, sep 2013.[69] M. Latva-Kokko and Daniel H. Rothman. Diffusion properties of gradient-based lattice boltzmann models of immiscible fluids.
Physical Review E ,71(5), may 2005.[70] I. Halliday, A. P. Hollis, and C. M. Care. Lattice boltzmann algorithm forcontinuum multicomponent flow.
Physical Review E , 76(2), aug 2007.4071] Zhaoli Guo and Chang Shu.
Lattice Boltzmann method and its applicationsin engineering , volume 3. World Scientific, 2013.[72] Qin Lou, Zhaoli Guo, and Baochang Shi. Evaluation of outflow boundaryconditions for two-phase lattice boltzmann equation.
Physical Review E ,87(6), jun 2013.[73] Long Li, Xiaodong Jia, and Yongwen Liu. Modified outlet boundary con-dition schemes for large density ratio lattice boltzmann models.
Journalof Heat Transfer , 139(5):052003, mar 2017.[74] Jingwei Huang, Feng Xiao, and Xiaolong Yin. Lattice boltzmann sim-ulation of pressure-driven two-phase flows in capillary tube and porousmedium.
Computers & Fluids , 155:134–145, sep 2017.[75] Yuze Hou, Hao Deng, Qing Du, and Kui Jiao. Multi-component multi-phase lattice boltzmann modeling of droplet coalescence in flow channelof fuel cell.
Journal of Power Sources , 393:83–91, jul 2018.[76] Victor W. Azizi Tarksalooyeh, Gábor Závodszky, Britt J. M. van Rooij,and Alfons G. Hoekstra. Inflow and outflow boundary conditions for 2dsuspension simulations with the immersed boundary lattice boltzmannmethod.
Computers & Fluids , 172:312–317, aug 2018.[77] Xiaoyi He, Qisu Zou, Li-Shi Luo, and Micah Dembo. Analytic solutionsof simple flows and analysis of nonslip boundary conditions for the latticeboltzmann BGK model.
Journal of Statistical Physics , 87(1-2):115–136,apr 1997.[78] Qisu Zou and Xiaoyi He. On pressure and velocity boundary conditionsfor the lattice boltzmann BGK model.
Physics of Fluids , 9(6):1591–1598,1997.[79] Joseph P. Morris. Simulating surface tension with smoothed particle hy-drodynamics.
International Journal for Numerical Methods in Fluids ,33(3):333–353, 2000.[80] Bruno Lafaurie, Carlo Nardone, Ruben Scardovelli, Stéphane Zaleski, andGianluigi Zanetti. Modelling merging and fragmentation in multiphaseflows with surfer.
Journal of Computational Physics , 113(1):134 – 147,1994.[81] Damien Violeau.
Fluid Mechanics and the SPH method theory and appli-cations . Oxford University Press, 2012.[82] Holger Wendland. Piecewise polynomial, positive definite and compactlysupported radial functions of minimal degree.
Advances in ComputationalMathematics , 4(1):389–396, Dec 1995.[83] Walter Dehnen and Hossam Aly. Improving convergence in smoothedparticle hydrodynamics simulations without pairing instability.
MonthlyNotices of the Royal Astronomical Society , 425(2):1068–1082, 2012.4184] Andrea Colagrossi and Maurizio Landrini. Numerical simulation of interfa-cial flows by smoothed particle hydrodynamics.
Journal of ComputationalPhysics , 191(2):448 – 475, 2003.[85] N. Grenier, M. Antuono, A. Colagrossi, D. Le Touzé, and B. Alessandrini.An hamiltonian interface sph formulation for multi-fluid and free surfaceflows.
Journal of Computational Physics , 228(22):8380 – 8393, 2009.[86] Nima Tofighi and Mehmet Yildiz. Numerical simulation of single dropletdynamics in three-phase flows using isph.
Computers & Mathematics withApplications , 66(4):525 – 536, 2013.[87] Eirik G. Flekkoy, Peter V. Coveney, and Gianni De Fabritiis. Foundationsof dissipative particle dynamics.
Phys. Rev. E , 62:2140–2157, Aug 2000.[88] J. Bonet and T.-S.L. Lok. Variational and momentum preservation aspectsof smooth particle hydrodynamic formulations.
Computer Methods inApplied Mechanics and Engineering , 180(1):97 – 115, 1999.[89] Athanasios Mokos, Benedict D. Rogers, and Peter K. Stansby. A multi-phase particle shifting algorithm for SPH simulations of violent hydrody-namics with a large number of particles.
Journal of Hydraulic Research ,55(2):143–162, sep 2016.[90] Kamil Szewc.
Développement d’une approche particulaire de type SPHpour la modélisation des écoulements multiphasiques avec interfaces vari-ables . PhD thesis, Université de Lorraine, June 2013.[91] Alex Ghaitanellis.
Modélisation du charriage sédimentaire par une ap-proche granulaire avec SPH . PhD thesis, Université Paris-Est, 2017. Thèsede doctorat dirigée par Violeau, Damien Mécanique des fluides Paris Est2017.[92] Kamil Szewc and Micha? Tadeusz Lewandowski. Further investigationof the spurious interface fragmentation in multiphase smoothed particlehydrodynamics, 2016.[93] T. Douillet-Grellier, F. De Vuyst, H. Calandra, and P. Ricoux. Influence ofthe spurious interface fragmentation correction on the simulation of flowregimes. In
Proceedings of the International 13th SPHERIC Workshop,June 26-28, Galway, Ireland , 2018.[94] Damien Violeau and Agnes Leroy. On the maximum time step in weaklycompressible SPH.
Journal of Computational Physics , 256:388 – 415, 2014.[95] A. Tafuni, J.M. Domínguez, R. Vacondio, and A.J.C. Crespo. A versa-tile algorithm for the treatment of open boundary conditions in smoothedparticle hydrodynamics gpu models.
Computer Methods in Applied Me-chanics and Engineering , 2018.[96] Carlos E. Alvarado-Rodríguez, Jaime Klapp, Leonardo Di G. Sigalotti,José M. Domínguez, and Eduardo de la Cruz Sánchez. Nonreflectingoutlet boundary conditions for incompressible flows using sph.
Computers& Fluids , 159:177 – 188, 2017. 4297] P. Kunz, I. M. Zarikos, N. K. Karadimitriou, M. Huber, U. Nieken, andS. M. Hassanizadeh. Study of multi-phase flow in porous media: Com-parison of SPH simulations with micro-model experiments.
Transport inPorous Media , 114(2):581–600, nov 2015.[98] R. W. Lockhart and R. C. Martinelli. Proposed correlation of data forisothermal two-phase, two-component flow in pipes.
Chemical EngineeringProgress , 45(1):39–48, 1949.[99] U Ghia, K.N Ghia, and C.T Shin. High-re solutions for incompressibleflow using the navier-stokes equations and a multigrid method.
Journalof Computational Physics , 48(3):387 – 411, 1982.[100] Jonas Latt, Bastien Chopard, Orestis Malaspinas, Michel Deville, andAndreas Michler. Straight velocity boundaries in the lattice boltzmannmethod.
Physical Review E , 77(5), may 2008.[101] Amir Banari, Christian F. Janssen, and Stephan T. Grilli. An improvedtwo-phase lattice boltzmann model for high density ratios : applicationto wave breaking. 2012.
Appendix A. Additional test cases
For all the following additional test cases, no kernel gradient correction orshifting or interface correction were used.
Appendix A.1. Lid-driven Cavity Flow
The goal of this section is to validate the implementation of SPH and LBMfor the single phase Navier-Stokes case. The test case chosen for this purpose isthe well-known 2D lid-driven cavity flow problem shown in Fig. A.24. This isa common problem in the fluid mechanics community and numerous referencesolutions performed with different numerical methods are available in the liter-ature. In this case, we use Ghia et al. solution as a reference [99]. Note thatGhia’s solution is also numerical. (0 , ,
1) (1 , , xy v x , v y = x , v y = x = U lid , v y = x , v y = Figure A.24: The 2D Lid-driven Cavity Flow problem Re = U lid Lν where U lid is the velocity of the imposed at the top boundary, ν is the kinematicviscosity and L is the characteristic length of the problem. The simulationswere performed for Re = 100 , , and and for ×
50 = 2500 , ×
100 = 10000 and ×
200 = 40000 particles or nodes (respectively forSPH and LBM). The density is set to / m , the velocity of the lid is U lid = 0 . l.u. for LBM and / s for SPH, the domain is L x × L y = 1 m × and the viscosity ν is adjusted to reach the desired Reynolds number.For LBM, due to stability issues, the MRT collision operator was used. Thestandard set of relaxation times S defined in Eq. (2.19). In order to have stableresults for at least one lattice size for every Reynolds number, a specific setupwas used where indicated (referred as LBM ∗ ). The relaxation times are thefollowing : S ∗ = diag (1 . , . , . , . , . , . , . , ω eff , ω eff ) (A.1)and the lid velocity is increased U ∗ lid = 0 . l.u. .The velocity boundary condition at the top boundary has been applied usingthe procedure described in Sec. 3.5 for SPH and in Sec. 2.5 for LBM. For theother boundaries, a no-slip boundary conditions has been applied. The simula-tions are terminated when a steady state is reached (i.e. (cid:114)(cid:80) a | ρ n +1 a − ρ na | ρ na < e − or after
60 s of real simulated time).
Appendix A.1.1. Re = 100 × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) Figure A.25: SPH results for Re = 100 × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) Figure A.26: LBM results for Re = 100 When Re = 100 , the MRT operator for LBM with the standard relaxationtimes S is able to simulate the test case for all grid resolutions that were con-44idered. As shown in Fig. A.29, both LBM and SPH are able to reproduce thevelocity field more and more accurately as the lattice/particles resolution is in-creased. However, LBM always present a higher order of convergence ( ≈ timesfaster). Moreover, LBM is the method that offers the best accuracy comparedwith Ghia et al.’s solution with an L discrepancy of ≤ . for the × lattice resolution. On the other hand, the SPH method shows a higher L discrepancy (at the boundaries in particular) with a maximum discrepancy of (cid:46) . for the × particles resolution.Concerning the spatial distribution of the flow shown in Figs. A.26 and A.25,LBM shows the appearance of two vortexes at the two bottom corners of thedomain which is in accordance with the theory. On the contrary, SPH is notable to reproduce those two vertexes but instead has flow perturbations in theconcerned areas.In fact the two expected vertexes at the corners are appearing during theSPH simulations but they are highly unstable. They keep forming (togetheror independently) and vanishing as the simulation progresses. It indicates thatSPH captures an instability in the correct areas but fails to reach a steady statethus the formation of spurious perturbations. Those vertexes being of smallintensity, their formation is probably affected by the boundary conditions. × - t = 2 .
48s 100 × - t =46 .
66s 200 × - t =54 . Figure A.27: SPH streamlines for Re = 100 at selected timesteps In Fig. A.27, one can note that the two expected vertexes at the cornersare in fact appearing during the SPH simulations but they are highly unstable.They keep forming (together or independently) and vanishing as the simulationprogresses. It indicates that SPH captures an instability in the correct areasbut fails to reach a steady state thus the formation of spurious perturbations.Those vertexes being of small intensity, their formation is probably affected bythe boundary conditions.
LBM - × SPH - × Figure A.28: Density fields for Re = 100 The density fields of the two methods for Re = 100 in Fig. A.28 show thatLBM has a smoother density field compared with SPH. As expected, due to the45hoice to use the weakly compressible approach, SPH presents a noisy densityfield. It is expected that an incompressible approach (ISPH) with a Poissonsolver would improve the quality of the density (and thus pressure) field. Theseobservations are valid for all four Reynolds numbers studied in this section. − . − . . . . . . . . . V x /U lid Y / L × -LBM × -LBM × -LBM × -SPH × -SPH × -SPHGhia et al. − − Number of nodes/particles L D i s c r e p a n c y SPH ( S = − . )LBM ( S = − . ) − . − . . . . . . . . . V y /U lid X / L × -LBM × -LBM × -LBM × -SPH × -SPH × -SPHGhia et al. − − Number of nodes/particles L D i s c r e p a n c y SPH ( S = − . )LBM ( S = − . ) Figure A.29: Re = 100 Appendix A.1.2. Re = 400 For Re = 400 , the MRT operator for LBM with the standard relaxation times S is able to simulate the test case for all grid resolutions that were considered.The superiority of LBM in terms of accuracy is magnified in that case.The LBM method shows more accurate results that SPH for all resolutionsconsidered as shown in Fig. A.33. The maximum discrepancy is always locatedat the domain’s boundaries. As an example, for the × resolution, bothLBM and SPH have a maximum difference with Ghia’s reference ≥ at theright boundary whereas it is ≤ and ≤ inside the domain for LBM andSPH respectively.Once again, LBM shows a better global accuracy for the same resolution anda higher order of convergence than SPH as shown in Fig. A.33. In particular,for the × resolution, the LBM L discrepancy on V x along the verticalcenterline is more than times lower than the SPH one. For V y along thehorizontal centerline, both methods have a comparable discrepancy.The streamlines plots of Figs. A.31 and A.30 are showing that LBM correctlyreproduces the existence of two vertexes at the bottom corners of the domain.For the × case, a spurious vertex appears at the top left corner of thedomain and is likely due to the combination of boundary conditions (bounceback+ Zou-He) at this location as it is smoothed out when the resolution increases.46 × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) Figure A.30: SPH results for Re = 400 × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) Figure A.31: LBM results for Re = 400 On the other hand, the SPH results are not able to simulate an establishedvertex pattern at the bottom corners. In the bottom right corner where thevertex is the strongest, for the × and × cases, a vertex appears tobe growing but is either not at the correct location or not with the correct am-plitude. In fact, when looking at earlier streamlines plots as shown in Fig. A.32,one can see that SPH does generate vortexes in the correct areas at selectedinstants during the simulation but they fail to stabilize and are continuouslyappearing and disappearing. × - t = 31 . s × - t =58 . s × - t = 52 . s Figure A.32: SPH streamlines for Re = 400 at selected timesteps Appendix A.1.3. Re = 1000 For Re = 1000 , the MRT operator with the standard relaxation times S failsto give stable results for the × resolution. However, when using anotherset of relaxation times S ∗ , one can obtain a stable solution. The impact ofempirically adjusting the relaxation times to "make it work" remains to beinvestigated. 47 . − . . . . . . . . . V x /U lid Y / L × -LBM × -LBM × -LBM × -SPH × -SPH × -SPHGhia et al. − − Number of nodes/particles L D i s c r e p a n c y SPH ( S = − . )LBM ( S = − . ) − . − . . . . . . . . . V y /U lid X / L × -LBM × -LBM × -LBM × -SPH × -SPH × -SPHGhia et al. − − Number of nodes/particles L D i s c r e p a n c y SPH ( S = − . )LBM ( S = − . ) Figure A.33: Re = 400 × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) Figure A.34: SPH results for Re = 1000 × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) Figure A.35: LBM results for Re = 1000
48s in the previous cases, one can observe in Fig. A.38 that LBM exhibits abetter accuracy than SPH for almost all resolutions. For the highest resolution,LBM has a maximum L discrepancy of ≈ . . For the same resolution, SPHgives a L discrepancy of ≈ . . Besides, the LBM order of convergence is stillup to − times higher than the SPH one.For this Reynolds number, it can be seen in Figs. A.35 and A.34 that SPHis capable of generating a vertex pattern at the bottom right corner for thetwo highest resolutions but it is unstable for the smallest resolution. Moderatedeviations of the flow indicating a potential growing vortex can be seen at thebottom left corner. When looking at the streamlines of the SPH simulations,we observe that all three resolutions are generating vertexes in the correct spotsat selected instants but only the × case manage to stabilize one at thebottom right corner.As previously said for the smaller Reynolds numbers, LBM is again showingthe appearance of the two vertexes at the correct locations. An instability isgrowing at the top left corner but disappears at the highest resolution. × - t =52 . s × - t =50 . s × - t =52 . s Figure A.36: SPH streamlines for Re = 1000 at selected timesteps For this Reynolds number, it can be seen in Figs. A.35 and A.34 that SPHis capable of generating a vertex pattern at the bottom right corner for thetwo highest resolutions but it is unstable for the smallest resolution. Moderatedeviations of the flow indicating a potential growing vortex can be seen at thebottom left corner. When computing the streamlines for selected timesteps ofthe SPH simulations as shown in Fig. A.37, it is seen that all three resolutionsare generating vertexes in the correct spots but only the × case manageto stabilize one at the bottom right corner.As previously said for the smaller Reynolds numbers, LBM is again showingthe appearance of the two vertexes at the correct locations. An instability isgrowing at the top left corner but disappears at the highest resolution. × - t =52 . s × - t =50 . s × - t =52 . s Figure A.37: SPH streamlines for Re = 1000 at selected timesteps . − . . . . . . . . . V x /U lid Y / L × -LBM* × -LBM × -LBM × -SPH × -SPH × -SPHGhia et al. − − Number of nodes/particles L D i s c r e p a n c y SPH ( S = − . )LBM ( S = − . ) − . − . . . . . . . . . V y /U lid X / L × -LBM* × -LBM × -LBM × -SPH × -SPH × -SPHGhia et al. − − Number of nodes/particles L D i s c r e p a n c y SPH ( S = − . )LBM ( S = − . ) Figure A.38: Re = 1000 Appendix A.1.4. Re = 10000 For Re = 10000 , the MRT operator despite its superior stability propertiescompared to BGK is unable to give stable results for none of the consideredlattice resolutions. Even using the set of relaxation times S ∗ , only the highestlattice resolution × prevents the simulation to blow up. In fact, Zou-heboundary conditions are known to be unstable at high Re and this is likely to beone the reasons the LBM simulations fail for the lowest resolutions considered. Itis possible to enhance the stability of the velocity boundary conditions, see [100]for example.For high Reynolds numbers, the flow is typically considered turbulent. Sinceno LBM nor SPH models considered in this study include the effect of turbu-lence, results are to be taken with caution. In consequence, both methods areshowing larger errors than in the previous cases where Re was much smaller.Nevertheless, the pattern is the same. As observed in Fig. A.43, LBM alwaysoffers a much better accuracy than SPH for the × resolution.In Fig. A.40, the LBM results are showing a high number of vertexes at thebottom right corner (5 vertexes), the bottom left corner (3 vertexes) and the topleft corner (2 vertexes). This is not agreeing with the theory where only 1 vertexis reported at the top left and bottom left corners and 2 vertexes at bottom rightcorner. These spurious vertexes could be due to the use of the MRT operatorwith relaxation times tuned based on a trial-and-error approach. The numberof vortexes is variable during the simulation as shown in Fig. A.42 where thenumber of vortexes is correct. Extra vortexes keep appearing and disappearingthroughout the simulation. No steady state is reached by the LBM in this case.The SPH streamlines plots of Fig. A.39 are not showing any vertex pattern50 × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) × - Velocity (left) - Streamlines (right) Figure A.39: SPH results for Re = 10000 × - Velocity (left) - Streamlines (right) Figure A.40: LBM results for Re = 10000 until the highest resolution is reached. For this × case, one can notethe appearance of a vertex at the top left corner, a small growing vertex at thebottom left corner and a growing vertex next to two very small vertexes at thebottom right corner. Those vertexes are stable through the simulation unlikethe one at the bottom left corner as suggested by Figs. A.41. Those figures alsoshow that at a smaller resolution, none of the expected vertexes are stable. Appendix A.2. Rayleigh-Taylor Instability
The Rayleigh-Taylor instability is a well-known two-phase problem in whicha heavy fluid is placed on top of a light fluid with a given interface shape andsubmitted to gravity. Several previous works have reproduced this case withSPH or LBM, for example [85, 90, 101, 58]. The test case and its parametersare borrowed from [85]. The computational domain is twice as high as long, H × L with H = 2 L and populated with nodes/particles. The densityratio is . while the viscosity ratio is . Gravity is set g = 9 .
81 m / s − forSPH and g = 1 × − l.u. and oriented downwards. Therefore, the viscosity ν is adjusted to match the desired Reynolds number Re = (cid:113) ( H/ gν = 420 .No surface tension is used. No slip boundary conditions are applied to thewalls. The interface is initialized as follows : y = 1 − sin(2 πx ) . Time t isnon-dimensionalized by t g = 1 / (cid:112) g/H . The distribution of the two phases isshown at selected timesteps in Fig. A.44, superposed with results from [85].Both methods are able to simulate the instability patterns as expected. Somedifferences are observable when t/t g ≥ in particular when the interface isstrongly distorted. LBM grows instabilities slightly faster than SPH and is closerto the behavior of the superposed Level-Set interface. On the other hand, our51 × - t =59 . s × - t =46 . s × - t =51 . s Figure A.41: SPH streamlines for Re =10000 at selected timesteps × Figure A.42: LBM streamlines for Re =10000 at selected timestep − . − . . . . . . . . . V x /U lid Y / L × -LBM* × -SPH × -SPH × -SPHGhia et al. − Number of nodes/particles L D i s c r e p a n c y LBMSPH ( S = − . ) − . − . − . . . . . . . . . V y /U lid X / L × -LBM* × -SPH × -SPH × -SPHGhia et al. − Number of nodes/particles L D i s c r e p a n c y LBMSPH ( S = − . ) Figure A.43: Re = 10000 t/t g = 5 located on both ends of the mushroom-likeshapes, but at an higher computational cost. (a) SPH (b) LBMFigure A.44: Phase distribution of Rayleigh-Taylor instability at selected timesteps : t/t g = 1 , and . Superposed with SPH interface (in black) and with Level-Set interface (in red) bothextracted from [85]. Appendix B. Transformation matrices M = +1 +1 +1 +1 +1 +1 +1 +1 +1 − − − − − − − − − − − − − − − − − − − − − − − − − , M − = +4 − − − − − − − − − − − − − − − − − − − − − − − − −