Comparison of multiphase SPH and LBM approaches for the simulation of intermittent flows

TL;DR

This paper compares multiphase SPH and LBM methods for simulating complex intermittent two-phase flows in pipes, highlighting their respective strengths and limitations through various test cases.

Contribution

It provides a quantitative comparison of SPH and LBM for multiphase flow simulation, focusing on their accuracy, robustness, and computational efficiency in pipe flow scenarios.

Findings

SPH is more robust and versatile.

LBM is more accurate and faster.

Both methods have specific strengths depending on the application.

Abstract

Smoothed Particle Hydrodynamics (SPH) and Lattice Boltzmann Method (LBM) are increasingly popular and attractive methods that propose efficient multiphase formulations, each one with its own strengths and weaknesses. In this context, when it comes to study a given multi-fluid problem, it is helpful to rely on a quantitative comparison to decide which approach should be used and in which context. In particular, the simulation of intermittent two-phase flows in pipes such as slug flows is a complex problem involving moving and intersecting interfaces for which both SPH and LBM could be considered. It is a problem of interest in petroleum applications since the formation of slug flows that can occur in submarine pipelines connecting the wells to the production facility can cause undesired behaviors with hazardous consequences. In this work, we compare SPH and LBM multiphase formulations where surface tension effects are modeled respectively using the continuum surface force and the color gradient approaches on a collection of standard test cases, and on the simulation of intermittent flows in 2D. This paper aims to highlight the contributions and limitations of SPH and LBM when applied to these problems. First, we compare our implementations on static bubble problems with different density and viscosity ratios. Then, we focus on gravity driven simulations of slug flows in pipes for several Reynolds numbers. Finally, we conclude with simulations of slug flows with inlet/outlet boundary conditions. According to the results presented in this study, we confirm that the SPH approach is more robust and versatile whereas the LBM formulation is more accurate and faster.

Equations

$$\displaystyle\bm{c}_{i}=\left\{\begin{array}[]{ll}(0,0),&i=0\\ \left[\sin(\theta_{i}),\cos(\theta_{i})\right],&i=1,3,5,7\\ \left[\sin(\theta_{i}),\cos(\theta_{i})\right]\sqrt{2},&i=2,4,6,8\end{array},\text{ with }\theta_{i}=\frac{\pi}{4}(4-i).\right.$$

$$N^{k}_{i}(\bm{x}+\bm{c}_{i},t+1)=N^{k}_{i}(\bm{x},t)+\Omega^{k}_{i}\big{(}N^{k}_{i}(\bm{x},t)\big{)},$$

$${\Omega }_i^k=({\Omega }_i^k{)}^{(3)}\left[({\Omega }_i^k{)}^{(1)}+({\Omega }_i^k{)}^{(2)}\right].$$

$$N^{k}_{i}(\bm{x},t_{*})=(\Omega^{k}_{i})^{(1)}\big{(}N^{k}_{i}(\bm{x},t)\big{)}.$$

$$N^{k}_{i}(\bm{x},t_{**})=(\Omega^{k}_{i})^{(2)}\big{(}N^{k}_{i}(\bm{x},t_{*})\big{)}.$$

$$N^{k}_{i}(\bm{x},t_{***})=(\Omega^{k}_{i})^{(3)}\big{(}N^{k}_{i}(\bm{x},t_{**})\big{)}.$$

$$N_i^k(\mathbf{x}+{\mathbf{c}}_i,t+1)=N_i^k(\mathbf{x},t_{\ast \ast \ast }).$$

$$({\Omega }_i^k{)}^{(1)}(N_i^k)=N_i^k-{\omega }_{\text{eff}}\left(N_i^k-N_i^{k(e)}\right),$$

$${\rho }_k=\sum _i{N}_i^k=\sum _i{N}_i^{k(e)},$$

$$\rho \mathbf{u}=\sum _i\sum _k{N}_i^k{\mathbf{c}}_i=\sum _i\sum _k{N}_i^{k(e)}{\mathbf{c}}_i,$$

$$N_i^{k(e)}({\rho }_k,\mathbf{u},{\alpha }_k)={\rho }_k\left({\varphi }_i^k+W_i\left[3{\mathbf{c}}_i\cdot \mathbf{u}+\frac{9}{2}({\mathbf{c}}_i\cdot \mathbf{u}{)}^2-\frac{3}{2}(\mathbf{u}{)}^2\right]\right).$$

$$\displaystyle W_{i}=\left\{\begin{array}[]{ll}4/9,&i=0\\ 1/9,&i=1,3,5,7\\ 1/36,&i=2,4,6,8\end{array}\right.$$

$$\displaystyle\phi^{k}_{i}=\left\{\begin{array}[]{ll}\alpha_{k},&i=0\\ (1-\alpha_{k})/5,&i=1,3,5,7\\ (1-\alpha_{k})/20,&i=2,4,6,8\end{array}.\right.$$

$$\gamma =\frac{{\rho }_r^0}{{\rho }_b^0}=\frac{1-{\alpha }_r}{1-{\alpha }_b},$$

$$p^k=\frac{3{\rho }_k(1-{\alpha }_k)}{5}={\rho }_k(c_s^k{)}^2.$$

$$\psi =\left(\frac{{\rho }_r}{{\rho }_r^0}-\frac{{\rho }_b}{{\rho }_b^0}\right)/\left(\frac{{\rho }_r}{{\rho }_r^0}+\frac{{\rho }_b}{{\rho }_b^0}\right).$$

$$\displaystyle\omega_{\text{eff}}=\left\{\begin{array}[]{ll}\omega_{r},&\psi>\delta\\ f_{r}(\psi),&\delta\geq\psi>0\\ f_{b}(\psi),&0\geq\psi\geq-\delta\\ \omega_{b},&\psi<-\delta\end{array},\right.$$

$$\displaystyle\begin{array}[]{ll}f_{r}(\psi)=\chi+\eta\psi+\kappa\psi^{2},\\ f_{b}(\psi)=\chi+\lambda\psi+\nu\psi^{2},\\ \end{array}$$

$$\displaystyle\begin{array}[]{ll}\chi=\left.\left(2\omega_{r}\omega_{b}\right)\middle/\left(\omega_{r}+\omega_{b}\right),\right.\\ \eta=\left.\left(2(\omega_{r}-\chi)\right)\middle/\delta,\right.\\ \kappa=\left.-\eta\middle/\left(2\delta\right),\right.\\ \lambda=\left.\left(2(\chi-\omega_{b})\right)\middle/\delta,\right.\\ \nu=\left.\lambda\middle/\left(2\delta\right).\right.\end{array}$$

$$({\Omega }_i^k{)}^{(1)}(N_i^k)=N_i^k-{\mathbf{M}}^{-1}\mathbf{S}\mathbf{M}\left(N_i^k-N_i^{k(e)}\right),$$

$$\mathbf{S}=\text{diag}\left(s_0,s_1,s_2,s_3,s_4,s_5,s_6,s_7,s_8\right).$$

$$N_i^{k(e)}({\rho }_k,\mathbf{u},{\alpha }_k)={\rho }_k{\varphi }_i^k+$$

$${\rho }_k{W}_i\left(\left[3{\mathbf{c}}_i\cdot \mathbf{u}\left[1+\frac{1}{2}\left(3(c_s^k{)}^2-1\right)\left(3{\mathbf{c}}_i-4\right)\right]+\frac{9}{2}({\mathbf{c}}_i\cdot \mathbf{u}{)}^2-\frac{3}{2}(\mathbf{u}{)}^2\right]\right).$$

$${\mathbf{U}}^k={\mathbf{M}}^{-1}{\mathbf{C}}^k,$$

$$\displaystyle\begin{array}[]{ll}C^{k}_{1}=3(1-s_{1}/2)(\partial_{x}Q^{k}_{x}+\partial_{y}Q^{k}_{y}),\\ C^{k}_{7}=3(1-\omega_{\text{eff}}/2)(\partial_{x}Q^{k}_{x}+\partial_{y}Q^{k}_{y}),\end{array}\text{ with}$$

$$\displaystyle\begin{array}[]{ll}Q^{k}_{x}=(1.8\alpha^{k}-0.8)\rho^{k}u_{x},\\ Q^{k}_{y}=(1.8\alpha^{k}-0.8)\rho^{k}u_{y}.\end{array}$$

$$\mathbf{F}=\nabla \psi .$$

$$\displaystyle\bm{\nabla}f(\bm{x})\approx\sum_{i}\xi_{i}\bm{c}_{i}f(\bm{x}+\bm{c}_{i}),\text{ with }\xi_{i}=\left\{\begin{array}[]{ll}0,&i=0\\ 1/3,&i=1,3,5,7\\ 1/12,&i=2,4,6,8\end{array}.\right.$$

$$({\Omega }_i^k{)}^{(2)}(N_i^k)=N_i^k+\sum _{\begin{array}{c}l\\ l\ne k\end{array}}\frac{A}{2}|{\mathbf{F}}_{kl}|\left[W_i\frac{(\mathbf{F}\cdot {\mathbf{c}}_i{)}^2}{|\mathbf{F}{|}^2}-B_i\right],$$

$$\displaystyle B_{i}=\left\{\begin{array}[]{ll}-4/27,&i=0\\ 2/27,&i=1,3,5,7\\ 5/108,&i=2,4,6,8\end{array}.\right.$$

$$\sigma =\frac{4}{9}\frac{A}{{\omega }_{\text{eff}}}.$$

Full text

††journal: arXiv

1 Introduction

Two-phase flows in pipes are encountered in various industrial applications and research areas. For instance, gas-liquid flows are present in chemical and heat transfer systems such as evaporators or refrigerators, but they can also be found in the transportation of oil and gas through pipelines in the petroleum industry. In particular, intermittent flow patterns (also known as slug or plug flow [1]) are undesirable in pipeline networks. Those slug patterns, that can measure up to tens of meters, are known to damage facilities (separator flooding, compressor starving, water hammer phenomenon, vibrations and fatigue) and to reduce flow efficiency. Therefore, it is crucial to be able to predict whether or not a slug flow regime will occur. If it does, the quantities of interest are the size 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 the appearance of a slug flow regime. At the industrial level, the first simulations [4, 5, 6] were done in the mid 80’s. Nowadays, in the oil and gas industry, 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 performing equally well on simple cases, they have trouble to simulate complex cases with a dominant gas phase. From an academic perspective, different models and methods 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 hydrodynamics (SPH) [15, 16] or phase field [17, 18], but these are applied on microfluidic problems for the most part. In this work, we will focus on two of them : 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 to model multiphase flows [19, 20, 21].

SPH is a Lagrangian meshfree method introduced in the late 70’s for astrophysics 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 to discretize partial differential equations such as the Navier-Stokes equations. On the other hand, it is a discrete Hamiltonian system composed of material points of constant mass that are tracked in time. Its pure meshless nature allows to get rid of many issues associated with meshing. However, it comes at some expenses too. Since there is no connectivity between particles, a neighbor searching procedure has to be carried for every particle at every time step, which is a serious bottleneck for code efficiency. Among the numerous applications of SPH, we can mention astrophysics [24], hydrodynamics [25], geophysics [26, 27, 28, 29] and computer graphics [30]. Some extensive reviews have been published [31, 32, 33].

LBM originates from two distinct approaches: the kinetic gas theory with discrete velocities and the lattice gas cellular automata method (LGCA or LGA [34]). In the late 80’s, probability density functions of particles were introduced within LGCA, giving birth to LBM. It consists in solving the Boltzmann equation in a discrete velocity space, which has been proven to be equivalent to solving the Navier-Stokes equations in the limit of low Mach and Knudsen numbers (as can be shown by a multiscale Chapman-Enskog expansion [35]). In practice, 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 with no differential terms, which makes it short, fast and particularly well suited for parallel execution [36]. Finally, traditional CFD methods are based on the calculation of macroscopic variables (velocity, pressure, density) whereas LBM tracks the evolution of probability distribution functions of particles [37]. LBM has been used for decades for flow simulations in complex geometries, especially in porous media [38, 39, 40, 41].

Our goal here is to propose a comparison of the multiphase SPH formulation presented in [42] and the color-gradient multiphase LBM formulation introduced in [43], on a collection of standard 2D test cases and on the simulation of slug flow regimes with periodic and inlet/outlet boundary conditions. To the best of our knowledge, this the first time such a comparison is presented.

We first detail the multiphase LBM and SPH formulations used in this work including surface tension models and boundary conditions. Then, we compare both approaches on static bubble tests with different density and viscosity ratios. Finally, we extend the comparison to two cases of slug flows, one induced by gravity with periodic boundary conditions and the other one based on inlet/outlet boundary conditions.

In addition, we provide in Appendix A a comparison of SPH and LBM on two test cases where one is a single-phase flow and the other does not involve surface tension : the lid-driven cavity flow and the Rayleigh-Taylor instability 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] for those looking for a detailed comparison of these techniques. In this work, we choose to work with the color gradient model because among the diffuse interface approaches, it is the one with a thin interface compared to the pseudo potential approach for example. In addition, the pseudo-potential model is cumbersome to use and to parametrize because there is coupling between the basic parameters [49]. The free energy model requires solving a Poisson equation at every time step which is time consuming and the mean field approach is limited to small density ratios [50]. Moreover, the color gradient model benefits from the large body of verification and validation cases available in the literature [51, 52, 53, 54, 55, 56, 49].

The present LBM approach is the two-phase model introduced in [43] and completed with the improvements proposed in [52, 51] for the recoloring operator and the color gradient. In addition, the contact angle ajustment is based on [56, 57] and the corrective procedure to properly recover Navier-Stokes equations is borrowed from [58]. We work with $2$ sets of distribution functions, one for each fluid, moving on a D2Q9 lattice. The associated velocity vectors are $\bm{c}_{i}$ with $i\in\left[0,8\right]$. As traditionally done in LBM, we set the lattice time step $\Delta t$ and the lattice space step $\Delta x$ to $1$. We can then define the 2D velocity vectors as follows :

[TABLE]

The distribution functions for a fluid of color $k$ (e.g. $k=r$ for red and $k=b$ for blue) are denoted $N^{k}_{i}(\bm{x},t)$, while $N_{i}(\bm{x},t)=N^{r}_{i}(\bm{x},t)+N^{b}_{i}(\bm{x},t)$ is used for the color-blind distribution function. In the rest of this section, the integer subscript $i$ is varying between $\left[0,8\right]$ while $k$ is referring to the color blue $b$ or red $r$ of the fluid. The lattice Boltzmann equation that describes the evolution of the system is then :

[TABLE]

where the collision operator $\Omega^{k}_{i}$ is the result of the combination of three sub operators (as previously done in [59]) :

[TABLE]

The Eq. (2.5) is solved using four consecutive steps which make use of the following operators :

Single phase collision step (see Sec. 2.1) :

[TABLE] 2. 2.

Perturbation step (multiphase collision 1/2) (see Sec. 2.2) :

[TABLE] 3. 3.

Recoloring step (multiphase collision 2/2) (see Sec. 2.3) :

[TABLE] 4. 4.

Streaming step :

[TABLE]

We will now examine in detail these steps as well as the specific treatments for the imposition of static contact angles and boundary conditions.

2.1 Single phase collision operators

2.1.1 BGK operator

The first sub operator, $(\Omega^{k}_{i})^{(1)}$ of Eq. (2.7), is the usual BGK operator for single phase LBM. The distribution functions are relaxed towards a local equilibrium as follows :

[TABLE]

where $\omega_{\text{eff}}$ is the effective relaxation rate. In practice, we first calculate the fluid densities based on the $0^{th}$ moment of the distribution functions :

[TABLE]

where $N^{k(e)}_{i}$ are the equilibrium distribution functions. The total fluid density is given by $\rho=\rho_{r}+\rho_{b}$, while the total momentum is computed as the $1^{st}$ moment of the distribution functions :

[TABLE]

where $\bm{u}$ is the density weighted arithmetic average velocity of the fluid. The equilibrium distribution functions $N^{k(e)}_{i}$ are defined in [43] by :

[TABLE]

These equilibrium distribution functions $N^{k(e)}_{i}$ are chosen to satisfy mass and momentum conservation and are based on the Maxwell-Boltzmann distribution. The weights $W_{i}$ are those of a standard D2Q9 lattice while the values $\phi^{k}_{i}$ depend on the density ratio. They are expressed as follows :

\displaystyle W_{i}=\left\{\begin{array}[]{ll}4/9,&i=0\\ 1/9,&i=1,3,5,7\\ 1/36,&i=2,4,6,8\end{array}\right.

\displaystyle\phi^{k}_{i}=\left\{\begin{array}[]{ll}\alpha_{k},&i=0\\ (1-\alpha_{k})/5,&i=1,3,5,7\\ (1-\alpha_{k})/20,&i=2,4,6,8\end{array}.\right.

As introduced in [60] for two-phase flows, the density ratio between red and blue fluids is $\gamma$, and must be computed as follows to obtain a stable interface :

[TABLE]

where the superscript [math] refers to the initial value of the density. Besides, the pressure of the fluid of color $k$ is :

[TABLE]

In Eq. (2.15), one of the $\alpha_{k}$ is actually a free parameter. In practice, if we assume that the least dense fluid is the blue one, we set a positive value for $\alpha_{b}>0$ ($0.2$ in this paper) and so we are certain that $0<\alpha_{b}\leq\alpha_{r}<1$ using Eq. (2.15). These parameters define the value of the sound speed $c_{s}^{k}$ in each fluid of color $k$ [43].

The effective relaxation parameter $\omega_{\text{eff}}$ is chosen so that the evolution Eq. (2.5) allows to recover the macroscopic Navier-Stokes equations for a single-phase flow in the single-phase areas. This parameter depends on the fluid kinematic viscosity $\nu_{k}$ through the following formula: $\omega_{k}=\left.1\middle/\left(3\nu_{k}+\frac{1}{2}\right)\right.$. However, when the viscosities of the fluids are different, the relaxation parameters are also different and an interpolation procedure is needed to define an effective relaxation parameter $\omega_{\text{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 fluid or interface), it is necessary to introduce a color field :

[TABLE]

The color field $\psi$ lies between $-1$ and $+1$. Within an interface, the color field is strictly between $-1$ and $+1$ whereas it is equal to $-1$ or $+1$ when located in an area that contains respectively only red fluid or blue fluid. The relaxation factor $\omega_{\text{eff}}$ in Eq. (2.5) is then evaluated as follows :

[TABLE]

in which $\delta$ is a free parameter and

[TABLE]

with :

[TABLE]

The parameter $\delta$ influences the thickness of the interface when the fluid viscosities are different. The larger $\delta$, the thicker the fluid interface. There is a trade off to find between robustness and interface sharpness. It is set to $0.1$ for all simulations in this paper. If the fluid viscosities are identical, $\delta$ is ineffective, because $\omega_{\text{eff}}=\omega_{r}=\omega_{b}$.

2.1.2 MRT operator

Alternatively, for the first sub operator, $(\Omega^{k}_{i})^{(1)}$, one can use the Multiple Relaxation Time (MRT) operator instead of the BGK operator. The MRT approach is more stable than its BGK counterpart [61, 58]. It reads as follows :

[TABLE]

where $\bm{S}$ is a diagonal matrix given by :

[TABLE]

The elements $s_{i},i\in\left[0\ldots 8\right]$ are the relaxation parameters. Following [54], $\forall i\in\left[0\ldots 6\right],s_{i}=\lambda\omega_{\text{eff}}$ while $s_{7}=s_{8}=\omega_{\text{eff}}$. As in [54], we choose $\lambda=4/5$. Moreover, $\bm{M}$ is a linear orthogonal transformation matrix that projects the distribution functions into the moment space. The discrete moment matrices $\bm{M}$ and $\bm{M}^{-1}$ are explicitly given in Appendix B for a D2Q9 lattice.

2.1.3 Proper recovery of Navier-Stokes equations

It has been emphasized in several papers [62, 63, 53] that the color gradient model does not fully recover the Navier-Stokes equations. An unwanted error term arises in the momentum equations when the density ratio is not one. Different techniques have been proposed to attenuate this issue [53, 64, 63, 58]. In the present work, we use the correction introduced in [58] for the MRT approach. It consists in two additions. First, a modified equilibrium distribution functions based on the $3^{rd}$ order Hermite expansion of the Maxwellian distribution [65, 66] is used instead of Eq. (2.14). It is expressed as follows :

[TABLE]

Second, a source term $\bm{U}^{k}$ is added in Eq. (2.32). It reads :

[TABLE]

where $\bm{C}^{k}=\left[0,C^{k}_{1},0,0,0,0,0,C^{k}_{7},0\right]^{T}$. The components $C^{k}_{1}$ and $C^{k}_{7}$ are computed as follows :

\displaystyle\begin{array}[]{ll}C^{k}_{1}=3(1-s_{1}/2)(\partial_{x}Q^{k}_{x}+\partial_{y}Q^{k}_{y}),\\ C^{k}_{7}=3(1-\omega_{\text{eff}}/2)(\partial_{x}Q^{k}_{x}+\partial_{y}Q^{k}_{y}),\end{array}\text{ with}

(2.38)

\displaystyle\begin{array}[]{ll}Q^{k}_{x}=(1.8\alpha^{k}-0.8)\rho^{k}u_{x},\\ Q^{k}_{y}=(1.8\alpha^{k}-0.8)\rho^{k}u_{y}.\end{array}

(2.41)

In particular, the derivatives are evaluated using a $9$-point isotropic finite difference approximation described shortly afterwards in Eq. (2.46).

2.2 Perturbation operator

In the color gradient model, surface tension forces are introduced by means of a perturbation operator shown in Eq. (2.8) [67, 43, 47]. To begin, it is needed to introduce the color gradient $\bm{F}$ of the color field $\psi$ (see Eq. (2.17)) that approximates the fluid-fluid interface normal. It is written as :

[TABLE]

In this work, a bi-dimensional S$2$I$4$ (spatial order $S=2$, isotropic order $I=4$) discrete gradient operator is used [68]. As shown in [51], this kind of gradient operator enhances the accuracy of the color-gradient model significantly while having the advantage to only rest on the nearest neighboring nodes. It takes the following form :

[TABLE]

The perturbation operator, for the fluid $k$, is defined by :

[TABLE]

where :

[TABLE]

The parameter $A$ which can vary in space and time handles the coupling between both fluids and is therefore linked with the surface tension coefficient $\sigma$. It is possible to predict the surface tension $\sigma$ between the two fluids using only the basic parameters of the model. As explained in [52], knowing the form of the expression describing the surface tension and performing simulations on planar interfaces, one can derive an expression that links the surface tension across an interface to the model parameters. For isotropic color gradients defined as in Eq. (2.46), the surface tension is set as follows :

[TABLE]

If $\sigma$ and $\omega_{\text{eff}}$ are fixed, one can obtain the value of $A$. Note that Eq. (2.52) is not universal and is susceptible to change if one uses a different color gradient or a different gradient operator. It has been shown in [43] that the perturbation operator allows to recover, within the macroscopic limit, the capillary stress tensor responsible for the surface tension forces in the macroscopic two-phase flow equations.

2.3 Recoloring operator

The recoloring operator $(\Omega_{i}^{k})^{(3)}$ of Eq. (2.9) is used to maximize the amount of fluid $k$ at the interface that is sent to the $k$ fluid region. It guarantees the fluid’s immiscibility. It respects the principles of mass and momentum conservation. The form of recoloring used in this paper is a combination of ideas taken from [69] and [70] and fully developed in [52]. It reads :

[TABLE]

where $\beta\in[0\ldots 1]$ is a parameter controlling the interface thickness [69] that will be set to $0.99$ unless otherwise stated. $\cos(\varphi_{i})=\left.\left(\bm{c}_{i}\cdot\bm{F}\right)\middle/\left(\lvert\bm{c}_{i}\rvert\lvert\bm{F}\rvert\right)\right.$ is the cosine of the angle between the color gradient $\bm{F}$ and the lattice direction vector $\bm{c}_{i}$. Note that for $i=0$ or $\lvert\bm{F}\rvert=0$, there is a division by [math]. In such a case, we set the whole term equal to 0 to respect mass conservation.

2.4 Adjustment of the color gradient for static contact angles

Based on [56], the imposition of a contact angle is performed using the method described in [57]. In order to properly describe the wetting boundary conditions, we divide the lattice nodes in two categories $C_{\text{f}}$ and $C_{\text{s}}$, each category being also subdivided in two subcategories $C_{\text{fs}}$, $C_{\text{ff}}$, $C_{\text{sf}}$ and $C_{\text{ss}}$.

$C_{\text{f}}$ : the set of fluid lattice nodes

• –

$C_{\text{fs}}$ : fluid lattice nodes in contact with at least one solid lattice node

• –

$C_{\text{ff}}$ : fluid lattice nodes not in contact with any solid lattice node

• -

$C_{\text{s}}$ : the set of solid lattice nodes

• –

$C_{\text{sf}}$ : solid lattice nodes in contact with at least one fluid lattice node

• –

$C_{\text{ss}}$ : solid lattice nodes not in contact with any fluid lattice node

When computing the color gradient in the fluid (i.e. for lattice nodes $\in C_{\text{f}}$), the knowledge of the color field at the boundary is required (i.e. for lattice nodes $\in C_{\text{sf}}$) because Eq. (2.46) involves the neighboring lattice nodes. Therefore, we need to extrapolate the color field to the boundary nodes, we do so using the following expression :

[TABLE]

It is now possible to evaluate the orientation of the color gradient in the fluid using the expression hereafter :

[TABLE]

In [57], they use an $8^{th}$ order isotropic stencil to compute the surface normal $\bm{n}^{s}$. In the subsequent simulations, boundary surfaces are flat and normals are known so we directly input the exact value. The correct color gradient orientation $\bm{n}$ depends on the prescribed contact angle $\theta$ and is evaluated as follows :

[TABLE]

where $D_{1}$ and $D_{2}$ are the Euclidean distances between the current unit normal vector $\bm{n}^{*}$ and the two possible theoretical unit normal vectors $\bm{n}^{1}$ and $\bm{n}^{2}$ of the interface at the contact line and are given by :

[TABLE]

with :

[TABLE]

Finally, once $\bm{n}$ is known, the corrected color gradient value is computed by taking :

[TABLE]

2.5 Boundary conditions

An overview of the available techniques to impose boundary conditions for single phase LBM can be found in [71]. For multiphase LBM, the literature is less abundant. In particular, the case of inlet/outlet boundary conditions is particularly difficult because specific treatments are needed to handle the interface when the fluids are entering and/or leaving the domain. Being able to have efficient inlet/outlet boundary conditions is attractive because it extends the 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 are used in this work :

No slip boundary conditions are imposed using full bounceback [77]. Then, free slip boundary conditions are also imposed using full bounceback except that the diagonally traveling distribution functions are sent forward along the wall rather than reflected back the way they came. Finally, velocity and/or pressure boundary conditions are imposed following [78]. 2. 2.

Periodic boundary conditions is a very useful tool in computational simulations because it allows to reduce the size of the simulated domain. The implementation of these boundary conditions consists in sending the distribution functions that are leaving the domain on one end to the other end of the domain as if the two sides of the domain were connected. 3. 3.

Inlet/outlet boundary conditions are achieved using a modified version of Zou-He boundary conditions [78]. The approach detailed in this paper is very similar 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 $\bm{u}^{b}_{\text{inlet}}$ and $\bm{u}^{r}_{\text{inlet}}$ from the north wall. The prescribed velocity fields $\bm{v}^{b,\text{in}}$ and $\bm{v}^{r,\text{in}}$ are designed so that $\bm{v}^{b,\text{in}}=\bm{u}^{b}_{\text{inlet}}$ on blue boundary lattice nodes and $\bm{v}^{b,\text{in}}=0$ on red boundary lattice nodes. Similarly, $\bm{v}^{r,\text{in}}=\bm{u}^{r}_{\text{inlet}}$ on red boundary lattice nodes and $\bm{v}^{r,\text{in}}=0$ on blue boundary lattice nodes. It is then possible to generate a color-blind prescribed velocity field $\bm{v}^{\text{in}}=\bm{v}^{b,\text{in}}+\bm{v}^{r,\text{in}}$. Following the classic Zou-He procedure described in [78], we can then compute the modified density and distribution functions. It reads :

[TABLE]

where $H^{\text{in}}$ is a corrective term that depends if we use Eq. (2.14) or Eq. (2.34) for the equilibrium. In practice, to derive Zou-He boundary conditions, one has to solve a linear system where one term comes from the equilibrium distribution function. If we use Eq. (2.34), we obtain extra terms compared to the classical Zou-He approach due to the Hermite expansion. It is computed as follows :

[TABLE]

It is then needed to redistribute these quantities in function of the color field value :

[TABLE]

Second, we will describe how we impose a constant pressure $p^{\text{out}}$ at the outlet located on the south wall. The corresponding prescribed densities are $\rho^{b,out}=\left.p^{\text{out}}\middle/(c_{s}^{b})^{2}\right.$ and $\rho^{r,out}=\left.p^{\text{out}}\middle/(c_{s}^{r})^{2}\right.$. The color-blind prescribed density is then $\rho^{\text{out}}=\rho^{b,out}+\rho^{r,out}$. In addition, we also have :

[TABLE]

where $H^{\text{out}}$ is evaluated as follows :

[TABLE]

We can then redistribute these quantities similarly with what was done in Eq. (2.79) :

[TABLE]

A test case has been setup to check the performance of these boundary conditions. Initial and final configurations can be found in Fig. 1. In Fig. 2 is shown the evolution of the inlet velocities and outlet pressure with the number of iterations. Note that quantities have been averaged along the height of the pipe. We can see that that we are recovering the prescribed values at the inlet and at the outlet after a transient period with a maximum error $\leq 2\%$. In Fig. 3, we see the distribution of the color field, the inlet velocity and the outlet pressure along the height of the pipe. It is possible to observe a velocity peak and a pressure peak located at the interface. This is likely due to the fact that fluids are mixed at the interface resulting in governing equations not being properly solved at this location. Moreover, slight discrepancies can be seen at the walls due to boundary conditions. Overall, the previously described boundary conditions are giving satisfactory results and will be used later in the paper.

3 SPH immiscible multiphase model

In this section, we will describe in details the SPH immisible multiphase model used in this work. This formulation and the associated open boundary conditions are taken from [42, 16].

3.1 Governing equations

For an incompressible fluid with a constant viscosity, the mass and momentum equation (completed with the equation of state) in a Lagrangian system are given as :

[TABLE]

with $\bm{u}$ fluid velocity, $\rho$ fluid density, $\mu$ the fluid viscosity, $\nu=\frac{\mu}{\rho}$ the kinematic viscosity, $p$ fluid pressure, $\mu$ fluid dynamic viscosity, $\bm{g}$ gravity, $c$ fluid speed of sound (here constant), $\gamma$ fluid adiabatic index, $\rho_{0}$ fluid initial density, $p_{0}$ background pressure, $\bm{F}^{st}$ is the surface tension force and $D/Dt$ denotes the material derivative following the motion.

The Tait’s equation of state (3.3) is added to Eqs. (3.1)-(3.2) to close the system. In this work, $\gamma$ will be set to $7.0$ for all fluids considered in all subsequent simulations. This is the so-called Weakly Compressible SPH formulation (WCSPH). Just like LBM, it is not a truly incompressible approach since the density is allowed to vary. This artificial compressibility has to be as weak as possible and is controlled by the speed of sound $c$. In this paper, given a reference length $L_{\text{ref}}$ and a reference speed $U_{\text{ref}}$, we used the following formulas [79] to set a value for $c$ and $p_{0}$ :

[TABLE]

with $\Delta\rho=0.01$ to enforce (not strictly) a maximum variation of $1\%$ of the density field and $\sigma^{\alpha\beta}$ the surface tension coefficient between phases $\alpha$ and $\beta$. The integer $N_{\text{phases}}$ is the number of different phases.

In order to model surface tension between fluids, an interaction force $\bm{F}^{st}$ is added to the momentum Eq. (3.2). Following [80], the continuum surface stress method introduces a surface tension force per unit volume that is expressed as the divergence of the capillary pressure tensor :

[TABLE]

with $\Pi$ the capillary pressure tensor defined by :

[TABLE]

where $\alpha,\beta\in\{1,\ldots,N_{\text{phases}}\}$ and $\Pi^{\alpha\beta}$ is expressed as :

[TABLE]

with $\tilde{\bm{n}}^{\alpha\beta}$ the unit normal vector from phase $\alpha$ to phase $\beta$, $\sigma^{\alpha\beta}$ the surface tension coefficient between phase $\alpha$ and phase $\beta$, $\delta^{\alpha\beta}$ a well-chosen surface delta function and $\bm{I}$ the identity matrix. For example, in the case of a three-phase system with a wetting phase $s$, a non-wetting phase $n$ and a solid phase $s$, the stress tensor reads $\Pi=\Pi^{ns}+\Pi^{ws}+\Pi^{nw}$.

3.2 SPH formulation

The SPH formulation adopted for this paper is identical to the one used in [16]. We will not here recall in detail the SPH discretization process and will assume that the reader already has experience with SPH. For those who are looking for an exhaustive description of this method, we recommend for example the reading of [81].

An SPH discretization consists of a set of points with fixed mass, which possess material properties and interact with its neighboring particles through a weighting function (or smoothing kernel) [23]. A particle’s support domain, $\Lambda$, is defined by its smoothing length, $h$, which is the radius of the smoothing kernel $W$. In all simulations presented in this paper, $h=2\Delta r$ where $\Delta r$ is the initial particle spacing. To obtain the value of a function at a given particle location, values of that function are found by taking a weighted (by the smoothing function) interpolation from all particles within the given particle’s support domain. An analytical differentiation of the smoothing kernel is used to find gradients of this function.

The interpolated value of a function $\bm{A}$ at the position $\bm{x}_{a}$ of particle $a$ can be expressed using SPH smoothing as :

[TABLE]

in which $\bm{A}_{b}=\bm{A}(\bm{x}_{b})$, $m_{b}$ and $\rho_{b}$ are the mass and the density of neighboring particle $b$. The set of particles $\Lambda_{a}=\{b\in\mathbb{N}\mid\lvert\bm{x}_{a}-\bm{x}_{b}\rvert\leq\kappa h\}$ contains neighbors of particle $a$ that lie within its defined support domain. The coefficient $\kappa$ depends on the choice of the kernel, it is equal to $2$ for the $5^{th}$ order $\mathcal{C}^{2}$ Wendland kernel function [82, 83]) used in this paper. For the sake of clarity, $W(\bm{x}_{a}-\bm{x}_{b},h)$ has been denoted $W_{ab}$. In 2D, this kernel is expressed as follows :

[TABLE]

Several multiphase formulations [84, 85, 86] have been proposed for SPH throughout the years. In this work, the formalism presented in [42] has been adopted. The density is directly evaluated through a kernel summation which gives an exact solution to the continuity equation. It reads :

[TABLE]

Discrete gradient and divergence operators in this formalism are given by :

[TABLE]

where $\Theta_{a}=\left.\rho_{a}\middle/m_{a}\right.$. It follows that the full multiphase SPH formulation for a particle $a$ is :

[TABLE]

where the viscous term $\nu\nabla^{2}\bm{u}$ has been discretized using the inter-particle averaged shear stress [87] and $\eta=0.01h$ is a safety factor to avoid a division by zero. Moreover, the evaluation of normal vectors is performed through the computation of the gradient of a color function $\chi$ defined for a given particle $a$ and a given phase $\alpha$ as :

[TABLE]

The normal vector $\bm{n}^{\alpha\beta}_{a}$ of particle $a$ belonging to phase $\alpha$ to the interface $\alpha\beta$ is then evaluated by :

[TABLE]

The surface delta function $\delta^{\alpha\beta}_{a}$ is chosen to be equal to $\lvert\bm{n}^{\alpha\beta}_{a}\rvert$ and $\tilde{\bm{n}}^{\alpha\beta}=\left.\bm{n}^{\alpha\beta}_{a}\middle/\lvert\bm{n}^{\alpha\beta}_{a}\rvert\right.$ for use in Eq. (3.7).

3.3 Corrective terms

Corrective terms are commonly used in SPH to remediate intrinsic issues of this formulation such as particles disorder or micro-mixing at the interface. Three different SPH corrections have been used in this work :

As suggested in [88], the kernel gradient is enhanced in order to restore consistency. For a given particle $a$, it reads :

[TABLE]

where $\bm{L}_{a}=\left(\sum_{b\in\Lambda_{a}}\frac{m_{a}}{\rho_{a}}\nabla W_{ab}\otimes(\bm{x}_{b}-\bm{x}_{a})\right)^{-1}$. Note that the tilde notation will be dropped in the rest of paper although the kernel gradient correction will be always used. 2. 2.

In order to maintain a good spatial distribution of particles and ensure a better accuracy, a shifting technique for multiphase flows has been used [89]. At the end of every timestep, all particles are shifted by a distance $\delta\bm{r}^{s}$ from their original position. This shifting distance of a particle $a$ is implemented through :

[TABLE]

where $C_{a}=\sum_{b\in\Lambda_{a}}\frac{m_{a}}{\rho_{a}}W_{ab}$ is the particle concentration, $\nabla C_{a}=\sum_{b\in\Lambda_{a}}\frac{m_{a}}{\rho_{a}}(C_{b}-C_{a})\nabla W_{ab}$ is the particle concentration gradient, $D_{a}$ is the diffusion coefficient, $\bm{s}$ and $\bm{n}$ are respectively the tangent and normal vectors to the interface light/heavy phase (with $\bm{n}$ oriented towards the light phase), $\beta_{n}$ is a reference concentration gradient (taken equal to its initial value) and $\alpha_{n}$ is the normal diffusion parameter set equal to $0.1$. The diffusion coefficient $D_{a}$ is computed as follows :

[TABLE]

where $A_{s}$ is a parameter set to $2$, $\bm{u}_{a}$ is the velocity of particle $a$, and $\Delta t$ is the time step. 3. 3.

Multiphase SPH can suffer from sub-kernel micro-mixing phenomena as highlighted in numerous publications [84, 85, 90, 91]. Around the interface, within a distance corresponding to the range of the smoothing length, particles tend to mix. It is because there is no mechanism that guarantees phases immiscibility in the surface tension’s continuum surface stress model. As suggested by the previously mentioned authors, we introduce a small repulsive force between phases as follows :

[TABLE]

where $\varepsilon=\left.L_{\text{ref}}\middle/h\right.$ for all simulations as suggested in [92] and where $L_{\text{ref}}$ is a reference length, typically the diameter of the pipe. The impact of this corrective force on the simulation of intermittent flows is evaluated in [93].

3.4 Time discretization and integration

Concerning time integration in SPH, different schemes are eligible. Among them, one can mention the Runge-Kutta, the Verlet or the Leapfrog schemes. In this work, the Predictor-Corrector Leapfrog scheme was adopted. It is a symplectic integrator which is recommended for SPH because of its conservative nature [81]. Indeed, SPH generally requires very small time steps resulting in a high number of iterations. The algorithm proceeds through the following steps. For every particle $a$,

Predictor step :

[TABLE] 2. 2.

Compute $\rho^{n}_{a}$ and $p^{n}_{a}$ using the corresponding expressions in Eq. (3.18). 3. 3.

Evaluate $\frac{D\bm{u}}{Dt}^{n}$ using the momentum equation in Eq. (3.18). 4. 4.

Corrector step :

[TABLE]

[TABLE]

The constant time step $\Delta t$ has to respect the Courant-Friedrichs-Lewy (CFL) criteria to ensure a stable evolution of the system, e.g.

[TABLE]

where, following [79], we have :

[TABLE]

A recent article [94] proposes a detailed investigation on maximum admissible time steps for WCSPH.

3.5 Boundary conditions

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 the following prescribed values for the pressure $p^{w}$, density $\rho^{w}$ and velocity $\bm{v}^{w}$ for a given ghost particle $g$ :

[TABLE]

with $V_{ga}=\sum_{a\in\Omega_{f}\cap\Lambda_{g}}\frac{m_{a}}{\rho_{a}}W_{ga}$, $\Omega_{f}$ the set of fluid particles and $\Lambda_{g}$ the set of neighboring particles of ghost particle $g$. A schematic drawn in Fig. 4 helps to visualize what is the intersection $\Omega_{f}\cap\Lambda_{g}$. Note that, for free-slip wall BC, one can use the same equations and change the sign of Eq. (3.29) for the direction where the free slip is allowed. Additionally, for velocity wall BC, the term $2\bm{v}_{0}$ (where $\bm{v}_{0}$ is the prescribed velocity) can be added to Eq. (3.29).

For periodic BC, different variants are available in the literature. In essence, particles close to a periodic boundary are allowed to interact with particles near an associated boundary. Quantities are exchanged both ways and if a particle leaves the domain through one side, it reenters the domain from the other side.

For inlet/outlet BC, the method extensively described in [16] (combining ideas from [95, 96]) has been used. The inlet and outlet boundaries are extended with a buffer layer of size $\kappa h$ to ensure a full kernel support. At the inlet, the goal is to inject particles with a prescribed velocity profile. At the outlet, particles need 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 prescribed velocity profile $\bm{v}^{p}$ and it carries the following values of pressure $p^{\text{in}}$, density $\rho^{\text{in}}$ and velocity $\bm{v}^{\text{in}}$ :

[TABLE]

with $V_{ia}=\sum_{a\in\Omega_{f}\cap\Lambda_{i}}\frac{m_{a}}{\rho_{a}}W_{ia}$ and $\Lambda_{i}$ the set of neighboring particles of inlet particle $i$.

On the other hand, at the outlet, a particle $o$ in the buffer is moved according to a smoothed convective velocity $\bm{v}^{conv}$. For example, if the outlet boundary is vertical and the flow leaves along the $x$ direction, it reads :

[TABLE]

with $V^{\prime}_{oa}=\sum_{a\in\Lambda_{o}}\frac{m_{a}}{\rho_{a}}W_{oa}$ the set of neighboring particles of outlet particle $o$. Note that in Eq. (3.33), the summation is over the full kernel support $\Lambda_{o}$ including fluid and outlet particles and not only over the intersection $\Omega_{f}\cap\Lambda_{o}$. Besides, particle $o$ also carries the following values of pressure $p^{\text{out}}$, density $\rho^{\text{out}}$ and velocity $\bm{v}^{\text{out}}$ :

[TABLE]

with $V_{oa}=\sum_{a\in\Omega_{f}\cap\Lambda_{o}}\frac{m_{a}}{\rho_{a}}W_{oa}$, $p^{p}$ and $\rho^{p}$ the prescribed pressure and density. Concerning the velocity, null cross velocities (here $v_{y}$) are enforced to ensure a divergence free velocity field at the outlet. An evaluation of these boundary conditions 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 of LBM 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 simulated and when a steady state is reached, the pressure difference between the exterior and the interior of the bubble is measured and compared to Laplace’s formula :

[TABLE]

with $\Delta P$ the pressure difference, $\sigma$ the surface tension coefficient, $R$ the bubble’s radius and $a$ the lateral dimension of the initial square bubble. Simulations are performed for three different resolutions : $60\times 60$, $100\times 100$ and $200\times 200$ nodes/particles. Four different combinations of density and viscosity ratios were tested : $(\frac{\rho_{\text{heavy}}}{\rho_{\text{light}}},\frac{\mu_{\text{heavy}}}{\mu_{\text{light}}})=(1,1)$, $(5,2)$, $(1000,1000)$ and $(1,100)$. The surface tension coefficient is $\sigma=1.88\leavevmode\nobreak\ \mathrm{N}\mathrm{/}\mathrm{m}$ for SPH and $\sigma=0.01\leavevmode\nobreak\ l.u.$ for LBM (where $l.u.$ stands for lattice units). The whole domain is $1\leavevmode\nobreak\ \mathrm{m}\times 1\leavevmode\nobreak\ \mathrm{m}$ and the lateral size of the initial square droplet is $a=0.33\leavevmode\nobreak\ \mathrm{m}$. The time is normalized by $t_{\sigma}=\sqrt{\rho a^{3}/\sigma}$. Note that, following [55], parameter $\beta$ in Eq.2.55 is adjusted when the resolution is increased taking the lowest resolution as a reference i.e. :

[TABLE]

with $\beta_{60\times 60}=0.99$.

Initially, when the density and viscosity ratios are set to one, one can observe in Fig. 5 the deformation of an initial square bubble towards a circular bubble under the influence of the surface tension. The Lagrangian/Eulerian difference between 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 in LBM, nodes are fixed and they switch phase to form the expect circular bubble. Besides, when the circular bubble is stabilized, both methods present residual velocities around the interface as shown in Fig. 6. However, those spurious currents are much more spread into the domain in SPH compared to LBM where they are localized around the interface. Note that, in the LBM color gradient framework, it is possible to significantly reduce the amplitude of spurious currents by choosing a more isotropic gradient operator [51] but, as it involves second range neighbors, it is more computationally expensive. it leads to. For SPH, Hu and Adams’ formulation [42] used in this paper has been reported to generate stronger spurious currents than other formulations [97].

In Fig. 7, the pressure profiles at steady state for the different resolutions and the different density and viscosity ratios considered are shown along with the corresponding $L_{2}$ error plots. First, one can clearly see that LBM is returning incorrect pressure values at the interface. Indeed, LBM presents non-physical pressure peaks at the bubble’s interface that tend to grow with the number of nodes. In fact, in the LBM color gradient method, the pressure is not well-defined at the interface. The pressure formula of Eq. (2.16) does not make sense at the interface where fluids are mixed and there is no mixture pressure defined in the considered framework. Hence, summing the fluid pressures is just an analytical construction that depends on the density profiles and the interface width. Thus, if we look at the LBM $L_{2}$ error along the whole horizontal centerline, we do not have mesh convergence since the error is growing at the interface. However, when we restrict the calculation of the $L_{2}$ error inside the bubble (i.e. when $0.4\leavevmode\nobreak\ \mathrm{m}<X<0.6\leavevmode\nobreak\ \mathrm{m}$), we do obtain a negative slope indicating mesh convergence.

Next, analyzing the impact of the density and viscosity ratios, for the case where $(\frac{\rho_{\text{heavy}}}{\rho_{\text{light}}},\frac{\mu_{\text{heavy}}}{\mu_{\text{light}}})=(1,1)$ in Figs. 7a and 7b, we get approximately the same order of convergence for both methods ($0.718$ and $0.838$ for LBM and SPH respectively) and the same error levels ($\leq 1\%$ at the bubble’s center) even though SPH always has a slightly higher error level. Next, when the density and viscosity ratios remains moderate (i.e. respectively up to $5$ and $2$ in Figs. 7c and 7d), both methods are under $2.5\%$ error compared to the reference solution. Additionally, we see that LBM offers a better order of convergence than SPH. In fact, LBM sees its order of convergence increased (from $0.718$ to $1.577$) compared to the previous case unlike SPH where it decreases (from $0.838$ to $0.392$). Moreover, LBM is less accurate than SPH for the two lowest resolutions $60\times 60$ and $100\times 100$ but performs better for the $200\times 200$ case thanks to its higher order of convergence. Overall, although both methods are returning satisfactory results for this case, we begin to observe a fall in performance whether it is for the order of convergence (for SPH) or for the error levels (for LBM and SPH) because gradients at stake are steeper. Then, for the high density ratio case where $(\frac{\rho_{\text{heavy}}}{\rho_{\text{light}}},\frac{\mu_{\text{heavy}}}{\mu_{\text{light}}})=(1000,1000)$ shown in Figs. 7e and 7f, we can see that the order of convergence of SPH remains roughly the same compared to the first case ($0.838$ vs $0.901$). The maximum error level is $\leq 6\%$ i.e. higher than the first case. This tends to indicate that when the density ratio increases, SPH is a quite robust and offers a reasonable accuracy for the same order of convergence. On the other hand, LBM sees its order of convergence strongly affected by the presence of this density ratio ($0.438$ vs $0.718$) while maintaining approximately the same error level. Finally, we looked at one last case, shown in Figs. 7g and 7h, where the density ratio is equal to $1$ and the viscosity ratio increased up to $100$. One can immediately note that, for both methods, the pressure profiles are heavily impacted at the interface (oscillations) in particular for the $60\times 60$ case. However, when looking at the pressure jump at the center of the bubble, LBM appears very robust to the presence of such a strong viscosity ratio. Indeed, we see that the error levels are of the same order than those of the first case and that the order of convergence is even higher ($0.718$ vs $1.604$). It shows that refining the lattice strongly helps to 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 of convergence is inferior to the first case ($0.838$ vs $0.658$). Moreover, the error does not seem to decrease anymore exponentially with the number of particles although more simulations would be needed to further check that statement.

To sum up, for limited density ratios and viscosity ratios, both methods are able to reproduce the pressure jump predicted by Eq. (4.1) with a good accuracy and with steep and clean pressure/density profiles. When the the density ratio increases up to $1000$, SPH seems to be more resilient than LBM in the sense that its order of convergence is not impacted by the presence of such a strong density ratio. LBM seems less robust in the same situation. On the contrary, when the viscosity ratio goes up to $100$, both methods render perturbed pressure profiles. However, it is SPH that appears to have more problem to handle a strong viscosity ratio whereas LBM maintains its performance level.

Contact angle case.

In addition, we compared the ability of the previously described implementations of SPH and LBM to prescribe a contact angle between a wetting phase, a non-wetting phase and a solid phase. Simulations are done with $100\times 100$ nodes/particles. The density and viscosity ratios are both equal to one. The whole domain is $1\leavevmode\nobreak\ \mathrm{m}\times 1\leavevmode\nobreak\ \mathrm{m}$ and the lateral dimensions of the initial rectangular droplet is $0.33\leavevmode\nobreak\ \mathrm{m}\times 0.165\leavevmode\nobreak\ \mathrm{m}$. For SPH, the surface tension coefficient between the wetting and non-wetting phase is $\sigma^{nw}=1.88\leavevmode\nobreak\ \mathrm{N}\mathrm{/}\mathrm{m}$ whereas the one between the wetting and the solid phase is set to $\sigma^{sw}=0\leavevmode\nobreak\ \mathrm{N}\mathrm{/}\mathrm{m}$ and the one between the non-wetting and the solid phase $\sigma^{sn}$ is adjusted to prescribed the desired contact angle $\theta^{\text{prescribed}}_{c}$ using the Young-Laplace equation $\theta_{c}=\frac{\sigma^{sw}-\sigma^{sn}}{\sigma^{nw}}$. For LBM, we follow the procedure described in Sec. 2.4. The surface tension coefficient is set to $\sigma=0.01\leavevmode\nobreak\ l.u.$. At steady state, the observed contact angle $\theta^{\text{observed}}_{c}$ is measured and reported in Fig. 8. The coefficient of determination is $\geq 0.99$ for both methods confirming that they can accurately reproduce a prescribed contact angle. Finally, one can observe the normalized velocity field at steady state for $\theta^{\text{prescribed}}_{c}=150^{\circ}$. The same comments made before are still valid, LBM spurious currents are less spread throughout the domain than in SPH. This is likely due to the Lagrangian nature of SPH where particles have to 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 of intermittent two-phase flows in pipes for different Reynolds numbers. The first one is periodic and gravity-driven while the second one is generated by a velocity inlet and a pressure boundary condition respectively at the inlet and outlet of the pipe.

5.1 Periodic case

In this section, we study the establishment of different periodic two-phase flow patterns under the influence of gravity $\bm{g}$ starting from a given bubbly flow. Following [15], the initial configuration is composed of $27.25\%$ of light phase and $72.75\%$ of heavy phase and is described in Fig. 10. All physical properties and simulation properties are in Tab. 1. The heavy phase viscosity $\mu_{l}$ is adjusted as function of the Reynolds number $Re=\frac{gH^{3}}{8\nu_{l}^{2}}$. The initial velocity field $\bm{V}$ is $\bm{V}(x,y)=\frac{g}{2\nu_{l}}y(H-y)$. Viscosity values and dimensionless numbers for each case are reported in Tabs. 2 and 3. No-slip boundary conditions are applied to the walls. The simulations is done with 50000 nodes/particles for $t=30\leavevmode\nobreak\ \mathrm{s}$. Four different Reynolds numbers were tested : $Re=10$, $50$, $100$ and $500$. The phase distributions, pressure fields and velocity fields at final state are shown in Figs. 11, 13 and 14 respectively.

In Fig. 11, it is possible to see that both methods reproduce the same flow pattern for all four $Re$ numbers considered. For $10\leq Re\leq 100$, we obtain a bubbly 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 pipe and the light phase travels in the middle. For $Re=10$, we observe that SPH present 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 inside light phase bubbles. Note that these small bubbles are to be absorbed by the main flow if the simulation lasted longer because they are moving slower than their environment. For $Re=50$ and $Re=100$, we obtain in all cases the same pattern made of three identical Taylor bubbles. Moreover, as shown in Fig. 12, the bubbles’ shapes between SPH and LBM for $10\leq Re\leq 100$ are very similar. Finally, as $Re$ grows, the Taylor bubbles are getting slightly shorter and higher in size. For $Re=500$, we again see that LBM offers a perfectly symmetric annular pattern. On the contrary, the bottom heavy phase layer in SPH is thicker than the top one. In general, LBM provides more symmetric results than SPH because of its Eulerian nature.

In Fig. 13, we can see that for $Re=10$, the pressure fields is dominated by the captured droplets of heavy phase inside the bubbles. For $Re=50$ and $100$, the pressure field reaches a maximum for SPH inside the bubbles whereas for LBM it is at the interface. Nevertheless, as predicted by Laplace’s law, the pressure 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 very similar. The same patterns surrounding the bubbles can be observed. For the annular case where $Re=500$, it is possible to see that the no-slip condition on the walls affects the flow more strongly in LBM than in SPH which results in a flatter velocity profile for the latter.

For $Re=50$ and $Re=100$, where the SPH and LBM patterns are the closest, we compared the density, velocity and pressure fields along the centerline on Figs. 15 and 16. Note that because the bubbles do not have the exact same position, we have shifted the LBM profiles from a fixed distance to be able to superpose the profiles. On the density plots of Figs. 15a and 16a, the different density treatment in both methods clearly appears. In LBM, the density is smoothed at the interface whereas in SPH, thanks to its Lagrangian nature, there is no interface smoothing in the density field because a given particle belongs to one phase or not, there is no intermediate state. Concerning the pressure fields shown in Figs. 15b and 16b, we observe that LBM suffers from the same overshoots at the interface that were described and explained in Sec. 4. On the other hand, the SPH pressure field is polluted with noise. Despite these discrepancies, both profiles are very close. Finally, in Figs. 15c and 16c, we see the velocity profiles in both methods have the same shape. The bubbles are moving at a much higher speed than the surrounding fluid (about $30\%$ faster). At each interface, the velocity reaches a local minimum. The only differences between both profiles is that in certain areas, the SPH velocity peaks have a smaller amplitude than in LBM. For example, the bubble velocity is the same for all three bubbles in LBM for $Re=50$ whereas for SPH the last bubble travels about $10\%$ faster than the other ones. One last comment is that in LBM at the interface, the velocity field present non-physical oscillations due to the fluids mixing at the interface. It is not the case in SPH because the pressure field does not suffer from pressure overshoots and fluids are not mixed at the interface.

To conclude, we can add that SPH and LBM are both well capable of simulating the transition from a given bubbly flow to a slug flow composed of Taylor bubbles for $Re\leq 500$. To further assess their relative performance, an extended comparison with other methods or with experimental data would be of great interest. Note that we have limited our study to $Re\leq 500$ because, for higher velocities and/or smaller viscosities, we lie outside LBM stability region whether because the low Mach rule is violated or because the relaxation time is too close from $0.5$.

5.2 Open channel case (inlet/outlet)

In this section, we study the ability of both methods to simulate a predicted intermittent flow regime. We consider an horizontal pipe of diameter $D=1\leavevmode\nobreak\ \mathrm{m}$ and length $L=10D$. 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 to be stratified with given volume fractions for each phase $\alpha_{g}=0.2$ and $\alpha_{l}=0.8$. All the physical properties are summarized in Tab. 4. Using these properties, it is possible to plot the flow regime map, see Fig. 18111In order to plot the map, one has to compute the Lockhart-Martelli [98] parameter which depends on $n$, $m$, $C_{g}$ and $C_{l}$. In this study, we used $n=m=2$ and $C_{g}=C_{l}=0.042$, and to pick an area to be investigated in the intermittent region. In this area, we adjust the viscosity to choose two cases that correspond to $Re=125$ and $Re=312.5$. Viscosity values and dimensionless numbers for each case are reported in Tabs. 5 and 6. Both cases were simulated in 2D with $25000$ nodes/particles. The simulation time was $30\leavevmode\nobreak\ \mathrm{s}$. At the inlet, each phase is injected with a constant velocity corresponding to its superficial velocity $u_{g,l}^{s}=\alpha_{g,l}u_{g,l}$. At the outlet, a constant pressure equal to the initial pressure $p_{0}$ is prescribed. Free slip boundary conditions are applied on the top and bottom walls (we were not able to generate a slug flow in LBM with no-slip boundary conditions under the same conditions). The initial setup is presented in Fig. 17. The phases distributions for each case are shown in Fig. 19 along with the associated plots showing the volume fraction evolution over time, the average pressure drop evolution over time, the heavy/light phase velocity evolution over time respectively in Fig. 20b, Fig. 22b and Fig. 23b.

In Fig. 19, one can see snapshots of the phases’ distribution for both methods at selected timesteps. It is clear that both methods are able to generate a slug flow as predicted by the flow map of Fig. 18. However, LBM produces a much more regular intermittent flow pattern with a higher slug frequency than SPH. This can also be seen in Figs. 20b and 21b. For example, for $Re=312.5$, the LBM slug frequency at the outlet is approximately $4.85\leavevmode\nobreak\ \mathrm{Hz}$ whereas for SPH it is close to $1.52\leavevmode\nobreak\ \mathrm{Hz}$. The slug frequency seems to remain roughly stable or to slightly increase (about $+1\leavevmode\nobreak\ \mathrm{Hz}$ for LBM and about $+0.16\leavevmode\nobreak\ \mathrm{Hz}$ for SPH for the highest peak) when $Re$ changes from $125$ to $312.5$ although it is less obvious in SPH periodograms due to the noise and composition of the signal. It is expected that the slug frequency increases when $Re$ increases but we could not raise $Re$ higher without making LBM simulations unstable ($\tau\to 0$ or $Ma\to 1$). In addition, SPH periodograms are noisier with $3$ to $4$ major frequency components unlike LBM where one frequency clearly emerges. It indicates that SPH volume fraction signals are less regular and are composed of signals with different frequencies. Moreover, we can observe in both methods that the point where the first slug appears is in general closer from the pipe entry when $Re$ is smaller. This is expected since when velocities are smaller at the entry, the first slug tends to form earlier in the pipe.

In Fig. 22b, the raw and average pressure drops evolution over time for all the different cases studied are shown. One can notice that in all cases, the pressure drops are $\leq 0$ on average (even though strong oscillations are observed). It means that the average pressure along the pipe’s height is higher at the outlet than at the inlet. This phenomenon was already seen in [16] and may be due to several factors : the quality of the boundary conditions, the recirculation areas at the entry that tend to lower the pressure and the averaging area chosen (i.e. $0.1\leavevmode\nobreak\ \mathrm{m}$ (from the inlet) $\times D$). Nevertheless, we observed that both methods are returning the same average pressure drop level around $\approx-0.05p_{0}$. SPH plots present much stronger oscillations than LBM ones (up to $7$ times the mean value for SPH against $3$ times the mean value for LBM). This is a known issue in weakly compressible SPH where the particles spatial distribution is directly linked to the density/pressure.

In Fig. 23b, the evolution over time of the heavy phase and light phase velocities within the whole pipe are shown for all cases considered. After a transient period of $\approx 4\leavevmode\nobreak\ \mathrm{s}$, the light and heavy phase velocities stabilize around a fixed value when a periodic state is reached. The main difference between SPH and LBM on this aspect is that the final velocity values are higher in LBM than they are in SPH. Overall, the average velocity field in LBM is $5\%v^{\text{in}}$ higher than in SPH. In our opinion, this is due to the wall boundary conditions that are handled in a very different way in both methods (full bounce back approach in LBM and interpolation-based approach in SPH) and to the wetting boundary conditions that is also handled differently. Those can strongly affect the flow, especially in dynamic cases like the ones considered in this section. From a general point of view, one can add that the oscillations observed in Figs. 22b and 23b are gaining amplitude when $Re$ increases in SPH whereas it is not the case in LBM for which they tend to reduce.

In a nutshell, as predicted by Taitel and Dukler’s flow map, both methods are capable to generate a slug flow pattern starting from the exact same simulation setup. Nevertheless, unlike all previous test cases for which the results were globally similar, we obtain significantly different flow patterns. Indeed, SPH produces bigger bubbles and in a more irregular way compared to LBM. We tend to believe that it is due to boundary conditions. In addition, we could not extend 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 be a 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 efficiency of both methods. On that aspect and despite very important progresses made in the SPH community, LBM remains superior in terms of speed, thanks to its local lattice-based algorithm that does not require a nearest neighbor search at every time step. For this paper, we have implemented both methods in the same framework in Fortran combined with an OpenMP library to handle multi-threading. For a simulation involving $40000$ nodes/particles, LBM was about $4$ times faster than SPH on a laptop equipped with a $2.9\leavevmode\nobreak\ \mathrm{G}\mathrm{Hz}$ Intel Core i7 processor with $4$ cores and $16$ Gb of RAM. This number is given as an indicative number and should be taken with caution since the codes were not fully optimized.

6 Conclusions

In this paper, we propose a 2D multiphase comparison of two particle methods, SPH and LBM, that are very different by nature. Despite these differences, both methods have been extensively used to model multiphase flows with success. We have chosen a multiphase formulation for each method among those available in the literature : the continuum surface force technique for SPH and the color gradient method for LBM. Then, we have simulated a collection of static bubble tests with different density and viscosity ratios with both methods. Finally, we have prolonged the comparison to more realistic cases, i.e. to the 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 outlet conditions in a multiphase context.

From a general point of view, we have confirmed that LBM offers a better order of convergence and a better accuracy than SPH although it suffers from a more narrow stability range than SPH. In many situations for which the Mach number is too high or the viscosity is too low, LBM will be unstable contrary to SPH which is only controlled by the CFL condition. Note that research on the extension of LBM to high Mach numbers is very active. In addition, SPH tends to generate pressure fields that are noisier than with LBM because of the Lagrangian behavior of particles whose position is directly linked to pressure through the density evaluation. This problem has been the subject of many investigations and several treatments are available. Moreover, LBM is more computationally efficient than SPH by construction.

On the multiphase aspects, both methods are very capable to simulate a variety of dynamic incompressible multiphase flow problems with good precision in the case of moderate density and viscosity ratios. However, we have noted that, on one hand, SPH appears more robust for high density ratios than LBM and, on the other hand, SPH has more trouble to handle high viscosity ratios than LBM. Another difference is that, at the interface, the fluids are mixed resulting in a diffuse interface whereas with SPH, particles are affected to one phase or the other without any mixing. More specifically, both methods have been 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 recommendation 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 characterize 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.