Numerical Study on Bubble Rising in Complex Channels Saturated with Liquid Using a Phase-Field Lattice-Boltzmann Method

Abstract

Packed bed reactors have been widely applied in industrial production, such as for catalytic hydrogenation. Numerical simulations are essential for the design and scale-up of packed beds, especially direct numerical simulation (DNS) methods, such as the lattice-Boltzmann method (LBM), which are the focus of future researches. However, the large density difference between gas and liquid in packed beds often leads to numerical instability near phase interface when using LBM. In this paper, a lattice-Boltzmann (LB) model based on diffuse-interface phase-field is employed to simulate bubble rising in complex channels saturated with liquid, while the numerical problems caused by large liquid-to-gas density ratio are solved. Among them, the channel boundaries are constructed with regularly arranged circles and semicircles, and the bubbles pass through the channels accompanied by deformation, breakup, and coalescence behaviors. The phase-field LB model is found to exhibit good numerical stability and accuracy in handing the problem of the bubbles rising through the high-density liquid. The effects of channel structures, gas-liquid physical properties, and operating conditions on bubble deformation, motion velocity, and drag coefficient are simulated in detail. Moreover, different flow patterns are distinguished according to bubble behavior and are found to be associated with channel structure parameters, gravity Reynolds number (Re_Gr), and Eötvös number (Eo).

1. Introduction

Packed bed reactors have the advantages of simple structure, convenient operation, low operating cost, good heat, and mass transfer performance, and have thus been widely used in industrial processes such as catalytic hydrogenation [1], oxidation reaction [2], nitrification reaction [3], and wastewater treatment [4], among others. For gas-liquid-solid three-phase packed bed reactors, the dynamic behaviors of the dispersed phase directly affect the mixing, mass transfer, and reaction efficiency between gas and liquid, which are extremely important for the design, optimization, and scale-up of the packed beds. However, the detailed dynamics of the dispersed phase and the complicated topological evolution of the gas-liquid interface are difficult to fully understand through experiments or numerical simulations on the reactor scale. The direct numerical simulation (DNS) provides us with a promising solution, in which the lattice-Boltzmann method (LBM) has attracted more and more attention due to its simple computation and clear physical background. At present, LBM has realized the precise tracking of the gas-liquid interface and the accurate prediction of the bubble (dispersed phase) dynamics from the mesoscopic scale [5]. However, the density ratios of liquid to gas in the actual packed beds are often as high as 1000, which tends to cause numerical instability at the gas-liquid interface in the simulations using LBM.

In recent decades, various lattice-Boltzmann (LB) models have been proposed for multiphase flows. Based on the approaches of describing the interactions between fluids, these models can be divided into four categories, including the color model [6], the pseudo-potential model [7], the free-energy model [8], and the kinetic model [9]. Undesirably, the original forms of these models are all accompanied by numerical instability defects when applied to multiphase systems with large density ratios in the range of 500–1000. In order to remedy the limitation, Inamuro et al. [10] developed an LB model that satisfies Galileo invariance for incompressible two-phase flows with density ratios of 50–1000 based on the free-energy model, but the implementation of this method requires high computational costs for solving the introduced pressure-Poisson equation, which destroys the simplicity of LBM. Lee et al. [11,12] successively proposed innovative LB models with three-step algorithms and second-order mixed difference schemes based on the Cahn-Hilliard model [13], achieving numerical stability at large density ratios up to 1000. In addition to the complicated and inefficient algorithms, some researchers [14,15] pointed out that the models of the Lee’s group [11,12] lacked accurate mass conservation due to the inconsistent invocations of central and biased finite differences. Zheng et al. [16] claimed to develop an LB model for multiphase systems with large density ratios by improving the free-energy model, but Fakhari and Rahimian [17] thought that their model was restricted to density-matched binary fluids, and incapable of modeling actual multiphase flows with large density ratios.

Over the past ten years, Fakhari and his group [17,18,19,20,21,22,23,24] have been dedicated to the LBM studies of multiphase systems with large density ratios, and have proposed various improved LB models and made a series of research progress. In particular, a phase-field LB model suitable for large density ratios has been constantly modified and optimized, which already has the characteristics of exact mass conservation, superior accuracy, and stability [19,20]. Until 2019, this LBM was used to solve various numerical problems such as bubble rising [17,19,21,24], droplet falling [17,19,20,22], von Kármán vortex street flows [18,20], and gas-liquid flows in porous media [23]. However, due to the complicated topological evolutions of the gas-liquid interface, including the bubble deformation, splitting, coalescence, and interplay with obstacles, no researchers have studied the bubble rising in the packed bed by this method so far. In order to save computing resources, the packed bed is simplified and replaced by two complex channels in our work, the bubble rising problem in the channels saturated with liquid will be investigated using the phase-field LB model.

In this paper, the phase-field LB model is tested by several benchmark problems and proved to indeed be a numerical method with high numerical stability and accuracy. Subsequently, this model is employed to simulate bubble rising in complex channels saturated with liquid at large density ratios, in which the problem of specifying the contact angle is also taken into account. The numerical results reveal the effects of channel structures, gas-liquid physical properties, and operating conditions on bubble deformation, motion velocity, and drag coefficient. Additionally, different flow patterns are recognized and found being dependent of channel structure parameters, gravity Reynolds number (Re_Gr), and Eötvös number (Eo). In the following, the specific numerical algorithm of phase-field LB model and its implementation details are presented in Section 2. Several verification examples used to test the phase-field LB model are listed in Section 3. Then, detailed numerical results and the discussions are given in Section 4. Finally, a summary and conclusions for the whole article are given in Section 5.

2. Numerical Method

In this paper, a phase-field LB model [20] is employed to solve the problem of bubble rising in complex channels saturated with liquid. This model is characterized by good accuracy and strict mass conservation, and has smaller spurious velocities and performs better in handling multiphase flows with large density and viscosity ratios than color and pseudo-potential models. In addition, the method also takes into account the three-phase contact line dynamics.

2.1. Phase-Field LB Model

The phase-field LB model we currently use is proposed by Fakhari et al. [20], in which the lattice-Boltzmann equations (LBE) combine a modified interface tracking equation and improved hydrodynamic evolution equations to achieve numerical computations of isothermal incompressible multiphase fluid systems.

2.1.1. Macroscopic Governing Equations

In the phase-field LB model, a phase-field variable ϕ is defined and its value is zero in the light fluid (gas bubbles) while one in the heavy fluid (liquid), which changes smoothly across the phase interface between two different fluids. Unlike those phase-field models based on the Cahn-Hilliard equation [13,17,25], a so-called conservative phase-field equation is used in the current model for interface tracking [26]:

$$\frac{\partial \varphi }{\partial t}+\nabla \cdot \left(\varphi \mathit{u}\right)=\nabla \cdot \left[M\left(\nabla \varphi -\frac{4}{\xi }\varphi \left(1-\varphi \right)\widehat{\mathit{n}}\right)\right]$$

(1)

where t is the time, u is the macroscopic velocity vector, M is the mobility, ξ is the interface thickness, and $\widehat{\mathit{n}}$ represents the unit vector normal to the gas-liquid interface, whose direction is away from the liquid side and points to the gas bubble side, that is

$$\widehat{\mathit{n}}=\frac{\nabla \varphi }{\left|\nabla \varphi \right|}$$

(2)

The isothermal incompressible multiphase flows are governed by Navier-Stokes equations, including the continuity equation:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \rho \mathit{u}=0$$

(3)

and the momentum conservation equation,

$$\rho \left(\frac{\partial \mathit{u}}{\partial t}+\mathit{u}\cdot \nabla \mathit{u}\right)=-\nabla p+\nabla \cdot \left(\mu \left[\nabla \mathit{u}+{\left(\nabla \mathit{u}\right)}^T\right]\right)+{\mathit{F}}_{\mathrm{s}}+{\mathit{F}}_{\mathrm{b}}$$

(4)

where ρ and μ represent the fluid density and viscosity respectively, p is the macroscopic pressure, F_b is the body force, and F_s is the surface tension force due to the presence of interface, which can be calculated by [27]:

$${\mathit{F}}_{\mathrm{s}}={\mu }_{\varphi }\nabla \varphi$$

(5)

in which μ_ϕ is the chemical potential of the binary-fluid system, which is defined by the derivative of the volumetric free energy E_f with respect to the phase field ϕ:

$${\mu }_{\varphi }=\frac{\delta E_{\mathrm{f}}}{\delta \varphi }=4\beta \varphi \left(\varphi -1\right)\left(\varphi -1/2\right)-\kappa {\nabla }^2\varphi$$

(6)

where the volumetric free energy is given by [27,28]:

$$E_{\mathrm{f}}={\displaystyle \oint \left(\beta \Psi \left(\varphi \right)+\frac{1}{2}\kappa {\left|\nabla \varphi \right|}^2\right)}dV$$

(7)

in which the coefficients $\beta =12\sigma /\xi$ and $\kappa =3\sigma \xi /2$ are related to the surface tension σ and the interface thickness ξ, and Ψ(ϕ) is the bulk free energy and expressed as:

$$\Psi \left(\varphi \right)={\varphi }^2{\left(1-\varphi \right)}^2$$

(8)

It should be noted that in the above equations, the gradient and Laplacian of the phase field ϕ are computed using a second-order isotropic central difference scheme. This scheme avoids the occurrence of a fourth-order derivative when solving the Laplacian in Equation (6) [13,25], and the combination of central and biased differences that was employed in previous studies [11,12], which weakens the numerical dispersion and maintains the conservation of mass and momentum. Therefore, the phase-field LB model is prone to numerical stability at large density ratios.

2.1.2. LBE for Interface Tracking

The following LBE has been recovered into Equation (1) by Geier et al. [29], which can be used to track the gas-liquid interface:

$$h_{\alpha }\left(\mathit{x}+{\mathit{e}}_{\alpha }\delta t,t+\delta t\right)=h_{\alpha }\left(\mathit{x},t\right)-\frac{h_{\alpha }\left(\mathit{x},t\right)-h_{\alpha }^{\mathrm{eq}}\left(\mathit{x},t\right)}{{\tau }_{\varphi }+1/2}$$

(9)

where h_α is the phase-field distribution function, τ_ϕ is the phase-field relaxation time, e_α represents the lattice-related mesoscopic velocity set. For the lattice of D2Q9 model, e_α is given by [30]:

$${\mathit{e}}_{\alpha }=c\left\{\begin{array}{lll}\left(0,0\right),& & \alpha =0\\ \left(\cos {\theta }_{\alpha },\sin {\theta }_{\alpha }\right),& {\theta }_{\alpha }=\left(\alpha -1\right)\pi /2,& \alpha =1~4\\ \left(\cos {\theta }_{\alpha },\sin {\theta }_{\alpha }\right)\sqrt{2},& {\theta }_{\alpha }=\left(2\alpha -9\right)\pi /4,& \alpha =5~8\end{array}\right.$$

(10)

where $c=\delta \mathit{x}/\delta t=1$, δx and δt represent the unit lattice length and unit time, respectively, and both δx and δt are set to one in uniform grids. $h_{\alpha }^{\mathrm{eq}}$ is the equilibrium phase-field distribution function, which is expressed as:

$$h_{\alpha }^{\mathrm{eq}}=\varphi {\Gamma }_{\alpha }+w_{\alpha }\frac{M}{c_s^2}\left[\frac{4}{\xi }\varphi \left(1-\varphi \right)\right]\left({\mathit{e}}_{\alpha }\cdot \widehat{\mathit{n}}\right)$$

(11)

in which

$${\Gamma }_{\alpha }=w_{\alpha }\left[1+\frac{{\mathit{e}}_{\alpha }\cdot \mathit{u}}{c_s^2}+\frac{{\left({\mathit{e}}_{\alpha }\cdot \mathit{u}\right)}^2}{2c_s^4}-\frac{\mathit{u}\cdot \mathit{u}}{2c_s^2}\right]$$

(12)

In Equations (11) and (12), c_s is the lattice sound speed and $c_s=c/\sqrt{3}$, and w_α represents the lattice-related weight coefficient set [30], where w₀ = 4/9, w_1–4 = 1/9, w_5–8 = 1/36. Mobility M is positively correlated with the phase-field relaxation time τ_ϕ as:

$$M={\tau }_{\varphi }{c}_s^2\delta t$$

(13)

After the calculation of a two-step collision-streaming sequence, the phase field is updated by taking the zeroth moment of the phase-field distribution function:

$$\varphi ={\displaystyle \sum _{\alpha }{h}_{\alpha }}$$

(14)

and subsequently the density ρ can be obtained by linear interpolation:

$$\rho ={\rho }_{\mathrm{g}}+\varphi \left({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}}\right)$$

(15)

where ρ_g and ρ_l represent the densities of the gas and liquid, respectively.

2.1.3. LBE for Hydrodynamics

An improved hydrodynamic LBE based on HCZ model [9] was developed by Fakhari et al. [17] to update pressure and velocity fields in nearly incompressible multiphase flows:

$${\overline{g}}_{\alpha }\left(\mathit{x}+{\mathit{e}}_{\alpha }\delta t,t+\delta t\right)={\overline{g}}_{\alpha }\left(\mathit{x},t\right)+{\Omega }_{\alpha }\left(\mathit{x},t\right)+F_{\alpha }\left(\mathit{x},t\right)$$

(16)

where ${\overline{g}}_{\alpha }$ is called the modified hydrodynamic distribution function, and F_α is the forcing term, which is calculated by [20]:

$$F_{\alpha }\left(\mathit{x},t\right)=\delta t\left[\left({\Gamma }_{\alpha }-w_{\alpha }\right)\left({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}}\right)c_s^2+{\Gamma }_{\alpha }{\mu }_{\varphi }\right]\left({\mathit{e}}_{\alpha }-\mathit{u}\right)\cdot \nabla \varphi +\delta t{\Gamma }_{\alpha }\left({\mathit{e}}_{\alpha }-\mathit{u}\right)\cdot {\mathit{F}}_{\mathrm{b}}$$

(17)

Ω_α is the collision operator, here the collision operator with a multiple-relaxation-time (MRT) model [31] is employed because its performance is more stable in the implementation than the Bhatnagar-Gross-Krook (BGK) model:

$${\Omega }_{\alpha }\left(\mathit{x},t\right)=-M^{-1}\widehat{S}M\left({\overline{g}}_{\alpha }-{\overline{g}}_{\alpha }^{\mathrm{eq}}\right)$$

(18)

where ${\overline{g}}_{\alpha }^{\mathrm{eq}}$ is called the modified equilibrium distribution function and is given by:

$${\overline{g}}_{\alpha }^{\mathrm{eq}}=g_{\alpha }^{\mathrm{eq}}-\frac{1}{2}{F}_{\alpha }$$

(19)

in which the equilibrium distribution function $g_{\alpha }^{\mathrm{eq}}$ is defined as:

$$g_{\alpha }^{\mathrm{eq}}=p{w}_{\alpha }+\rho c_s^2\left({\Gamma }_{\alpha }-w_{\alpha }\right)$$

(20)

M is an orthogonal transformation matrix to transform the distribution functions from physical space into moment space [31] and M⁻¹ is its inverse matrix, and $\widehat{S}$ is a diagonal relaxation matrix. For the D2Q9 lattice, $\widehat{S}$ can be selected as:

$$\widehat{S}=\mathrm{diag}\left(1, 1, 1, 1, 1, 1, 1, s_v, s_v\right)$$

(21)

in which

$$s_v=\frac{1}{\tau +1/2}$$

(22)

where τ is the hydrodynamic relaxation time, which is calculated according to the dynamic viscosity derived by linear interpolation:

$$\mu ={\mu }_{\mathrm{g}}+\varphi \left({\mu }_{\mathrm{l}}-{\mu }_{\mathrm{g}}\right)$$

(23)

and

$$\tau =\frac{\mu }{\rho c_s^2}$$

(24)

Among them, μ_g and μ_l represent the dynamic viscosities of the gas and liquid, respectively.

After solving the modified hydrodynamic distribution function according to the two-step calculation sequence of collision-streaming, the hydrodynamic properties are obtained by:

$$\mathit{u}=\frac{1}{\rho c_s^2}{\displaystyle \sum _{\alpha }{\overline{g}}_{\alpha }}{\mathit{e}}_{\alpha }+\frac{\delta t}{2\rho }\left({\mathit{F}}_{\mathrm{s}}+{\mathit{F}}_{\mathrm{b}}\right)$$

(25)

$$p={\displaystyle \sum _{\alpha }{\overline{g}}_{\alpha }}+\frac{\delta t}{2}\left({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}}\right)c_s^2\mathit{u}\cdot \nabla \varphi$$

(26)

Note that the calculation of pressure p requires the updated velocity u, so the update sequence of velocity must precede the pressure.

2.2. Numerical Implementation

2.2.1. Discretization

As described in Section 2.1, a second-order isotropic central difference scheme for the phase field is adopted in the phase-field LB model, which is conducive to stable implementation while retaining the second-order accuracy. In detail, the gradient terms used in Equations (5), (17), and (26) and the Laplacian term in Equation (6) are computed according to the following discrete schemes:

$$\nabla \varphi =\frac{1}{2c_s^2\delta t}{\displaystyle \sum _{\alpha =1}^8{\mathit{e}}_{\alpha }{w}_{\alpha }\left[\varphi \left(\mathit{x}+{\mathit{e}}_{\alpha }\delta t\right)-\varphi \left(\mathit{x}-{\mathit{e}}_{\alpha }\delta t\right)\right]}$$

(27)

$${\nabla }^2\varphi =\frac{1}{c_s^2\delta t^2}{\displaystyle \sum _{\alpha =1}^8{w}_{\alpha }\left[\varphi \left(\mathit{x}+{\mathit{e}}_{\alpha }\delta t\right)-2\varphi \left(\mathit{x}\right)+\varphi \left(\mathit{x}-{\mathit{e}}_{\alpha }\delta t\right)\right]}$$

(28)

2.2.2. Curved Boundary Treatment

If the solid walls are set in the computational domain, the interaction between the fluids and the solid walls can be described by imposing a specified contact angle at the solid boundary [28]:

$${\left.{\widehat{\mathit{n}}}_{\mathrm{w}}\cdot \nabla \varphi \right|}_{{\mathit{x}}_{\mathrm{w}}}=\Theta {\varphi }_{\mathrm{w}}\left(1-{\varphi }_{\mathrm{w}}\right)$$

(29)

where x_w is the position of a point on the solid boundary, ${\widehat{\mathit{n}}}_{\mathrm{w}}$ is a unit vector normal to the solid boundary, with its direction pointing away from the solid wall, and ϕ_w is the phase field value of the point on the solid boundary. Θ is related to the equilibrium contact angle θ and is written as:

$$\Theta =-\sqrt{\frac{2\beta }{\kappa }}\cos \theta$$

(30)

As shown in Figure 1, the solid (black dots) is presented with a curved solid boundary (black solid line); then Equation (29) is modified to [20]:

$${\left.{\widehat{\mathit{n}}}_w\cdot \nabla \varphi \right|}_{{\mathit{x}}_{\mathrm{w}}}={\left.\frac{\partial \varphi }{\partial {\mathit{n}}_{\mathrm{w}}}\right|}_{{\mathit{x}}_{\mathrm{w}}}=\frac{{\varphi }_{\mathrm{m}}-{\varphi }_{i,j}}{2h}=\Theta {\varphi }_{\mathrm{w}}\left(1-{\varphi }_{\mathrm{w}}\right)$$

(31)

where ϕ_m is the phase field value of the interpolated point in the fluid (blue dots), and $h=\left|{\mathit{x}}_{\mathrm{m}}-{\mathit{x}}_{\mathrm{w}}\right|$ is the distance from the interpolated point to the solid boundary. ϕ_w is estimated by interpolation:

$${\varphi }_{\mathrm{w}}=\frac{{\varphi }_{\mathrm{m}}+{\varphi }_{i,j}}{2}$$

(32)

Equation (31) can be solved by replacing ϕ_w with Equation (32):

$${\varphi }_{i,j}=\frac{1}{a}\left(1+a-\sqrt{{\left(1+a\right)}^2-4a{\varphi }_{\mathrm{m}}}\right)-{\varphi }_{\mathrm{m}}$$

(33)

where

$$a=h\Theta =-h\sqrt{\frac{2\beta }{\kappa }}\cos \theta \ne 0\left(\theta \ne {90}^{\circ }\right)$$

(34)

As for the neutral wetting conditions (θ = 90°), Equation (31) has a trivial solution ϕ_i,j = ϕ_m. To find ϕ_i,j and impose the contact angle conditions on the solid boundary, ϕ_m is the only unknown parameter, which can be obtained by the bidirectional interpolation scheme by utilizing the phase field values of four nearby nodes, and the coordinate positions of nearby nodes are depicted in Figure 1. If the unknown ϕ_i,j is required for interpolation, it is replaced by the ϕ_i,j from the previous time step.

After specifying the wetting boundary conditions, the unknown incoming distribution functions from the solid boundary nodes towards the fluid nodes can be determined by referring to the model of Yu et al. [32]:

$$h_{{\alpha }^{-}}^{\ast }\left({\mathit{x}}_{\mathrm{s}}\right)=h_{\alpha }^{\ast }\left({\mathit{x}}_{\mathrm{f}}\right)$$

(35)

$${\overline{g}}_{{\alpha }^{-}}^{\ast }\left({\mathit{x}}_{\mathrm{s}}\right)=\frac{\Delta }{1+\Delta }\left[{\overline{g}}_{{\alpha }^{-}}^{\ast }\left({\mathit{x}}_{\mathrm{f}}\right)+{\overline{g}}_{\alpha }^{\ast }\left({\mathit{x}}_{\mathrm{f}}\right)\right]+\frac{1-\Delta }{1+\Delta }{\overline{g}}_{\alpha }^{\ast }\left({\mathit{x}}_{\mathrm{ff}}\right)$$

(36)

where

$$\Delta =\frac{\left|{\mathit{x}}_{\mathrm{f}}-{\mathit{x}}_{\mathrm{w}}\right|}{\left|{\mathit{x}}_{\mathrm{f}}-{\mathit{x}}_{\mathrm{s}}\right|}$$

(37)

The subscript α⁻ denotes the incoming distribution functions, such that ${\mathit{e}}_{{\alpha }^{-}}=-{\mathit{e}}_{\alpha }$, and the superscript asterisk denotes the pre-streaming or post-collision state of the distribution functions.

3. Numerical Validation

3.1. Laplace Law

Laplace law states that the pressure difference between the inside and outside of a stationary bubble is proportional to the surface tension and inversely proportional to the bubble radius, which is often used to verify the numerical accuracy of the multiphase LB model [8,11,16,17]. The Laplace law for a two-dimensional bubble is written as:

$$\Delta p=p_{i\mathit{n}}-p_{o\mathit{u}t}=\frac{\sigma }{R_{\mathrm{b}}}$$

(38)

where R_b is the radius of the bubble. The computational domain of this verification is a square domain with 201 × 201 grids, and periodic boundary conditions are applied at its four boundaries. At the initial moment of the simulation, the center of the single bubble is located at (101,101), and the bubble is immersed in the liquid phase without any external force. The physical properties of the liquid and bubble are fixed as ρ_l/ρ_g = 1000 and μ_l/μ_g = 100. Zheng et al. [16] found that the simulation results were in better agreement with the analytical solutions and the spurious current was smaller when the interface thickness was greater than 4.5 lattice units (lu). Thus, the interface thickness is fixed to 6 lu in this Laplace-law test, and the bubble radius and surface tension are adjusted in the range of 10–50 lu and 0.01–0.1, respectively.

Figure 2a,b show the phase field contour of a bubble that has reached a stable state and the distribution curve of the phase field in the radial direction through the bubble, respectively. The analytical curve of the phase field distribution in Figure 2b is derived from the equilibrium phase field profile of diffuse-interface models:

$$\varphi \left(\mathit{x}\right)=\frac{1}{2}\left[1-\tanh \left(\frac{\left|\mathit{x}-{\mathit{x}}_0\right|}{\xi /2}\right)\right]$$

(39)

where x₀ is the location of the bubble interface. By comparison, it can be seen that the numerical results are almost the same as the analytical solutions.

The relationship between pressure difference, surface tension, and radius of the bubble is plotted in Figure 3, which is found to be highly consistent with Laplace law with deviations less than 3% and proves the accuracy of the current phase-field LB model quantitatively.

3.2. Bubble Deformation

The rising process of a single bubble under buoyancy action in a rectangular liquid-filled channel is often studied and used to test numerical models [33,34,35,36]. The single bubble will reach a relatively stable state with a constant rising velocity and terminal shape after rising for some time. In our simulations, a static circular bubble with a diameter D_b = N_x/5 is initially placed at (N_x/2, N_y/4) in a rectangular computational domain discretized with N_x × N_y = 512 × 1536 grid cells. A volumetric buoyancy force ${\mathit{F}}_{\mathrm{b}}=\left(\rho -{\rho }_{\mathrm{l}}\right)G_y\widehat{\mathit{y}}$, where G_y is the magnitude of the gravitational acceleration and $\widehat{\mathit{y}}$ is a unit vector with a vertical downward direction, is continuously imposed on the bubble. The periodic boundary conditions are used at the left and right boundaries of the computational domain, while the bounce-back boundary schemes are applied at the top and bottom boundaries. The density and viscosity ratios of the liquid and bubble are fixed as ρ_l/ρ_g = 1000 and μ_l/μ_g = 100. Besides, there are several dimensionless parameters for characterizing bubble dynamics:

(1)

The gravity Reynolds number (Re_Gr),

$$R{e}_{\mathrm{Gr}}=\frac{\sqrt{G_y{\rho }_{\mathrm{l}}\left({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}}\right)D_{\mathrm{b}}^3}}{{\mu }_{\mathrm{l}}}$$

(40)

(2)

The Eötvös number (Eo),

$$Eo=\frac{G_y\left({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}}\right)D_{\mathrm{b}}^2}{\sigma }$$

(41)

(3)

The Morton number (Mo),

$$Mo=\frac{G_y\left({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}}\right){\mu }_{\mathrm{l}}^4}{{\sigma }^3{\rho }_{\mathrm{l}}^2}$$

(42)

(4)

The above three are not independent, since $Mo=E{o}^3/R{e}_{\mathrm{Gr}}^4$.

Table 1. Comparison of terminal bubble shapes observed in experiments and predicted by front tracking method, previous and present LBMs.

Case	ReGr	Eo	Experiments (Bhaga and Weber, 1981) [37]	Front Tracking Method (Hua and Lou, 2007) [33]	LBM (Liang et al., 2019) [36]	Present LBM
A1	1.67	17.7
A2	79.88	32.2
A3	134.63	115
A4	30.83	339
A5	49.72	641

lists the terminal bubble shapes under different Re_Gr and Eo obtained by present LBM and other experimental or numerical methods [33,34,35,36,37]. By comparison, the terminal bubble shapes in our simulations can be observed to be highly similar to other research results, which achieves numerical stability under large liquid-to-gas density ratio conditions and demonstrates the qualitative accuracy of the current LB model.

During the bubble rising, the Reynolds number of the bubble $Re=u_{\mathrm{g}}{D}_{\mathrm{b}}{\rho }_{\mathrm{l}}/{\mu }_{\mathrm{l}}$ is continuously recorded. The comparison between the Re obtained from the present LBM and the experimental results is listed in

Table 2. Comparison of Re obtained from the present LBM and the experimental results.

Case	Re of Experiments [37]	Re of Present LBM	Relative Error (%)
A1	0.232	0.211	9.05
A2	55.3	47.8	13.56
A3	94.0	87.5	6.91
A4	18.3	16.4	10.38
A5	30.3	27.5	9.24

. The acceptable relative errors further confirm the numerical accuracy of the phase-field LB model.

4. Numerical Results and Discussion

4.1. Channel Construction and Numerical Initialization

Figure 4 shows the schematics of two computational domains in 2D Cartesian systems for the bubble rising problem in complex channels, in which the lateral boundaries are set as periodic boundaries, while the bounce-back boundary schemes are adopted at the top and bottom. Among them, the blue circular areas represent the gas bubbles, the red areas are the liquid that fills the channels and surrounds the bubbles, and the white circular and semicircular areas are the artificially constructed bounce-back domains, which are referred to as solid particles and not involved in the iterative computations of phase field and hydrodynamic properties.

As depicted in Figure 4, two types of symmetric channels are formed by arranging the solid particles in rectangular and triangular matrices, and they are named wavy vertical channel and S-shaped curved channel, respectively. D_b and D_p marked in the figures are the bubble and particle diameters, the horizontal spacing L of the particles refers to the shortest distance between two horizontally adjacent particles and similarly the vertical spacing H is the shortest distance between two vertically adjacent particles, and S represents the shortest distance between two diagonally adjacent particles. Hereinafter we use the particle spacing normalized by the bubble diameter to characterize the channel width.

Initially, a static and stable circular bubble is placed in the middle of the channel and submerged in the liquid with a diameter D_b and a volumetric buoyancy force ${\mathit{F}}_{\mathrm{b}}=\left(\rho -{\rho }_{\mathrm{l}}\right)G_y\widehat{\mathit{y}}$. The phase field in the solid particles is fixed to 0, and in other computational domains the phase field is initialized by the following hyperbolic tangent profile:

$$\varphi \left(\mathit{x}\right)=\frac{1}{2}\left[1+\tanh \left(\frac{R\left(\mathit{x}\right)-R_0}{\xi /2}\right)\right]$$

(43)

where R₀ is the initial bubble radius and represents the position of the gas-liquid interface in the radial direction of the bubble, and R(x) is the distance from any position x in the flow domains to the center of the bubble. Unless otherwise specified, in the following simulations, the initial bubble diameter and particle diameter are set as D_b = D_p = 51 lu. The density and viscosity ratios of the liquid and gas are fixed to ρ_l/ρ_g = 1000 and μ_l/μ_g = 100, respectively. Taking into account the balance of numerical stability and accuracy at large density ratios [22], the mobility and interface thickness are taken as M = 0.03 and ξ = 5 lu in the present simulations, respectively. As for the wetting properties, the contact angle of gas-liquid interface on particle surfaces is specified as 40°.

After initialization, the bubble rises through the channel under the action of buoyancy, and simultaneously undergoes a series of evolutions such as deformation, fragmentation, and coalescence. During these processes, as an important dimensionless parameter related to the bubble dynamics, the drag coefficient (C_D) is constantly monitored:

$$C_D=\frac{4G_y\left({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}}\right)D_{\mathrm{b}}}{3{\rho }_{\mathrm{l}}{u}_{\mathrm{g}}^2}$$

(44)

where u_g is the average velocity of bubbles rising, which is obtained by averaging the instantaneous y-direction velocities over all gas-phase nodes ($\rho <\left({\rho }_{\mathrm{l}}+{\rho }_{\mathrm{g}}\right)/2$).

4.2. Grid Independence

The channel structure and the effect of the channel on the bubble are greatly affected by the grid, so the case of the bubble rising in a rectangular channel without particles at Re_Gr = 100 and Eo = 20 is selected for testing the grid independence. Five different grids with N_x × N_y = 96 × 1536, 128 × 1536, 176 × 1536, 256 × 1536 and 512 × 1536 grid cells are applied in the bubble rising problem, respectively. The number of grid cells in the y-direction is sufficiently large and constant because its influence on the bubble steady velocity is negligible. The rising velocity variation of the bubble with D_b = N_x/5 over the dimensionless time is present in Figure 5, and the gravity-based dimensionless time is defined as $t*=t\sqrt{G_y/D_{\mathrm{b}}}$. It can be observed that when the grid cells are less than 176 × 1536, the bubble terminal velocity is unstable, and the calculation deviation is unacceptable. Thus, the minimum number of grid cells used in the following simulations is set to 176 × 1536.

According to the grid independence test results, considering that the channel width is one of the investigated variables affecting the bubble movement, the two computational domains of wavy and S-shaped channels are discretized using 176 × 1536 to 512 × 1536 and 192 × 1536 to 512 × 1536 grid cells, respectively.

4.3. Mass Conservation

In order to verify the mass conservation of the current model, the evolution of the total system mass versus the dimensionless time is recorded in Figure 6 during a numerical simulation of bubble rising in a wavy vertical channel with L/D_b = 0.719 and H/D_b = 0.875 at Re_Gr = 100 and Eo = 20. It is observed that the variation of the normalized total system mass M/M₀ is much less than 10⁻⁶ over a long period of time, indicating that the mass of the gas-liquid two-phase system is conserved well using current phase-field LB model.

4.4. Channel Width Effect

Undoubtedly, one of the variables that has the most obvious impact on bubble movement in the aforementioned channels is the channel width.

Figure 7 shows the evolution curves of bubbles rising velocity and drag coefficient versus the dimensionless time at Re_Gr = 100 and Eo = 20. In the wavy vertical channel, the bubble is prone to a stabilized state with a periodically fluctuating terminal velocity, along with C_D fluctuating within the certain ranges. In Figure 7, compared with the channel without particles, the rising velocity of the bubble is lower, and the drag coefficient is larger in the channel when the particles are arranged, which indicates that the presence of particles has a significant hindrance to the movement of the bubble.

u_g and C_D under different channel widths are given in Figure 8, which are calculated by averaging the regularly fluctuating u_g and C_D over a certain period of time in Figure 7. C_D is found to be lower as the channel width becomes wider, and the bubble rises faster. The same results are observed by Patel et al. [38] in the studies on the effects of the amplitude of the sinusoidal channel walls on bubble dynamics. In addition, the horizontal spacing between particles has a more severe impact on the bubble movement than the vertical spacing of the particles. Because the horizontal spacing of the particles directly affects the interaction between the bubble and the particles, and determines the difficulty for the bubble to pass through each gap, while the vertical spacing affects the fluctuation frequency of the bubble drag force, which can only indirectly act on the drag force. In view of this, unless otherwise stated, the following numerical results are obtained from the simulations in the channels with H/D_b = 0.875.

4.5. Surface Tension Effect

Surface tension is an important parameter to characterize the deformation properties of the gas-liquid interface. As for the bubble, a large surface tension means that the bubble is not easy to deform and break. In our study, since the total bubble mass is well conserved, we count the ratio of the total length of the gas-liquid interface to the area occupied by the bubble to measure the deformation rate of the bubble. As shown in Figure 9, the bubble deformation rate is proved to reduce with the increase of surface tension, which confirms that the bubble with large surface tension tends to maintain its original shape and resist deformation.

Furthermore, the declining u_g and ascending C_D are observed in Figure 10 when the surface tension is gradually increasing. This is because the bubble resists deformation more strongly when colliding with the particles, which restricts the bubble rising velocity. In contrast, the bubble with smaller surface tension is more likely to pass through the narrow gaps in the channel by deforming or even breaking. This phenomenon is consistent with the numerical results of Patel et al. [38], and similar to the studies on the effect of surface tension on liquid penetration by Shi et al. [39].

4.6. Bubble Diameter Effect

The effects of bubble diameter on bubble dynamics are depicted in Figure 11. The bubble diameter in this part is not equivalent to the particle diameter; it varies in the range of 25–204 lu, but the particle diameter is still fixed at 51 lu. Figure 11a indicates that it is more difficult for a bubble with larger diameter to pass through the channel due to the greater blocking effects of the particles, so the bubble rising velocity is slowed down, resulting in the ascent of drag coefficient, as demonstrated in Figure 11b.

4.7. Driving Force Effect

An additional driving force source term F_d can be added to the right side of the governing Equation (4) to simulate the changes in the pressure difference of the packed bed in the actual industry. Here we use the dimensionless ratio F_d/(G_y·ρ_l) to measure the magnitude of the driving force.

As illustrated in Figure 12, when the incremental additional driving forces are imposed on the bubble based on the presence of buoyancy force, u_g increases linearly while C_D decreases rapidly. Because the additional driving force pushes the bubble to pass through the channel faster and more smoothly, the rising velocity is significantly increased, so that u_g increases and C_D decreases.

4.8. Bubble Flow Pattern

During the numerical simulations, we found that the fluid properties and operating conditions set in the numerical initializations have significant impacts on the evolutions of the bubble. In order to classify the different bubble evolution processes, the gravity Reynolds number and Eötvös number that include fluid property and operating condition parameters are used as the classification criteria.

For the wavy vertical channel, mainly four types of bubble flow patterns are identified, and their detailed bubble evolution diagrams are listed in Figure 13. The four flow patterns are named Aw, Bw, Cw, and Dw, respectively; their subscript “w” represents the wavy vertical channel, and the capital letters are used to distinguish different flow patterns.

By observing Figure 13, the flow pattern Aw describes that the bubble is completely blocked by the two particles at the entrance of the channel, and even cannot enter the channel. It is inferred from the small Eo that the bubble surface tension in this flow pattern is relatively large, which makes it difficult for the bubble to enter the channel through deformation when encountering obstacles. In contrast, the bubble in the flow pattern Bw will become elongated when it flows into the narrow gaps thanks to the reduced surface tension, thus passing through the channel integrally. In flow pattern Cw, whenever the bubble crosses the gaps and collides with the particles, some secondary bubbles are split from the tail of the parent bubble. Then these secondary bubbles will gradually merge and expand following the parent bubble, and the coalesced secondary bubble with smaller rising resistance will catch up with and merge into the parent bubble finally. As for the flow pattern Dw, compared to the flow pattern Cw, more secondary bubbles are generated after the parent bubble interacts with the particles due to the extremely small surface tension, and the parent bubble will eventually be split into multiple secondary bubbles and dispersed in the liquid. As a summary, a lower Eo leads to a steady rise of the bubble, while a higher Eo instigates the dynamics instability and causes a multiple breakup in the form of secondary bubbles, which is in good agreement with the numerical results of Patel et al. [38] on bubble dynamics in sinusoidal channels using the level set method.

For the S-shaped curved channel, according to different bubble evolutions, five flow patterns As, Bs, Cs, Ds, and Es are distinguished, as illustrated in Figure 14. Among them, flow pattern As is described as the bubble cannot enter and pass through the channel due to the obstruction of the narrow channel. In flow pattern Bs, the bubble enters the channel by being squeezed and elongated, and it shuttles left and right to rise through the channel. For flow patterns As and Bs, the bubbles retain good integrity because of the higher surface tension. As for the flow pattern Cs, the bubble first splits into two after colliding with the particle, and then the two secondary bubbles will coalesce under their interaction. In another case, the two separated secondary bubbles will rise simultaneously with negligible mutual interference due to their far distance, which is named flow pattern Es. Similar to the flow pattern Dw, the bubble in flow pattern Ds is gradually broken into multiple small secondary bubbles and they are dispersed in the liquid phase, continuously breaking and coalescing. It is worth mentioning that the bubbles in the flow patterns Dw and Ds have better dispersibility and higher specific surface area, which are more conducive to gas-liquid mixing and mass transfer in actual industrial production.

In our studies, multiple sets of simulations have been executed within the Re_Gr range of 0–350 and Eo range of 0–250. Figure 15 depicts the divisions of flow patterns while changing Re_Gr and Eo with fixed channel widths L/D_b = 0.719 and H/D_b = 0.875 in the wavy vertical channel and L/D_b = 1.344 and H/D_b = 0.875 in the S-shaped curved channel. Eo is found to have more significant effects on the flow patterns compared to Re_Gr, and lower Eo numbers often correspond to the bubbles with higher integrity because the surface tension that is related to Eo number plays a leading role in bubble deformation. In the S-shaped channel, the flow pattern Ds occupies a wider range than the flow pattern Dw in the wavy channel, indicating that the bubbles are more likely to reach the state of fragmentation and dispersion because of more frequent collisions with particles.

It is worth reminding that the trend of the boundaries between different flow patterns in an S-shaped curved channel is similar to the shape regime curve of the single bubble under gravitational motion in the work of Clift et al. [40].

It can be qualitatively summarized from Figure 15 that the bubble becomes more broken and dispersed as Re_Gr and Eo increase, and this is attributed to the changes in liquid viscosity and bubble surface tension that have significant impacts on the gas-liquid interface. The increase of Re_Gr and Eo leads to the reduction of liquid viscosity and surface tension, respectively. The reduced liquid viscosity causes the weakening of the viscous effects, which results in the increase of bubble wobbling [38], and the lowered surface tension further promotes bubble deformation and breakup.

In addition to the Re_Gr and Eo, channel width also has an important influence on the bubble flow patterns. The flow pattern divisions according to the channel widths of two types of channels are demonstrated in Figure 16. The horizontal spacing between particles is found to have more obvious impacts on the flow of bubbles compared to the vertical spacing. The horizontal spacing of the particles directly affects the interaction between the bubble and the particles, thereby affecting the bubble shapes and flow patterns.

On the other hand, in the wavy channel, the changes in channel width do not contribute much to the transitions of bubble flow patterns while maintaining the same Re_Gr and Eo. Conversely, in the relatively complex S-shaped channel, the transitions of the flow patterns are more sensitive to the variations of channel width; four different flow patterns are revealed by changing the channel width. Under this Re_Gr and Eo condition, the flow pattern Cs is more conducive to the even distribution of the dispersed phase in practice because too-small channel width is not prone to the disintegration and flow of the dispersed phase, and too large channel width reduces the contact between the dispersed phase and the particles.

5. Conclusions

In this paper, the numerical studies on two-dimensional bubble rising in complex channels saturated with liquid at large density ratios within a wide range of gravity Reynolds numbers and Eötvös numbers have been implemented using phase-field LB model. The main research conclusions of this work are as follows:

(1)

The present LB model is tested through three aspects of Laplace law, bubble deformation, and mass conservation, and it has been proven to have good stability, accuracy, and conservation from both qualitative and quantitative perspectives.

(2)

In the simulations of bubble rising in complex channels, the effects of channel width, surface tension, bubble diameter and additional driving force on bubble motion are investigated in detail. The larger channel width and additional driving force as well as smaller bubble diameter and surface tension lead to lower drag coefficients, which are conducive to smooth passage through the channels for the bubble.

(3)

Four and five types of bubble flow patterns are divided according to different bubble evolution processes under different Re_Gr, Eo and channel structures conditions in the wavy vertical channel and S-shaped curved channel, respectively. The detailed flow pattern diagrams are drawn for flow pattern recognition. To some extent, this study has some guiding significance for the regulation of bubble flow patterns in the industrial packed beds.