Multiphase flows are ubiquitous, and the dynamics of droplets and bubbles is one of the core areas of research in the field of fluid mechanics. Understanding the dynamics of these fluids is vital to the design and optimization of numerous engineering applications, such as liquid-fueled rocket engines, blood pumps, diesel engines, gas turbines, cooling systems, and inkjet printers. In addition, the study of bubble detachment and droplet formation is of great significance for applications such as electronic spray coating and pharmaceutical applications.

Comprehensive studies on the dynamics of bubbles and droplet deformation were compiled by Clift et al. [1]. The authors primarily presented experimental observations, along with some computational results. Over the last few decades, the Lattice Boltzmann Method (LBM) has been recognized as a promising approach to the simulation of multiphase flows, offering advantages for the complex interfacial dynamics between different phases. Several LBM models have been proposed to simulate multiphase flows. The chromodynamic model, originally proposed by Gunstensen et al. [2] for the simulation of immiscible fluids using the lattice gas method developed by Rothman [3], and later extended by Grunau et al. [4], is one of the earliest models. Later, Swift et al. [5] proposed the free-energy model, where the free-energy concept is incorporated into the LBM. Among the multiphase Lattice Boltzmann models available, the pseudo-potential model of Shan and Chen [6] has become extremely popular due to its simplicity, flexibility, and ease of implementing force interactions. In the pseudo-potential model, the non-local interactions of the particles are computed based on the additional force applied to the nodes of the lattice. Thus, the pseudo-potential model enables the simulation of phase separation phenomena. An important advantage of the Shan-Chen model is that it allows simulation of different phase viscosity ratios using multiple distribution functions.

Takada et al. [7] used a free-energy model to simulate the motion of buoyancy-driven bubbles and validated the results using the Volume of Fluid (VOF) method. Gupta and Kumar [8] used the Shan-Chen model for simulating the dynamics of single and multiple bubbles under the influence of gravity and buoyancy forces. Another form of LBM for interface dynamics has been proposed by He, Shan, and Doolen [9]; this model is known as the HSD model. Later, He et al. [10] extended their model for incompressible flows by using two distinct distribution functions to solve for the pressure and velocity fields and to track interface dynamics. Qin et al. [11] proposed a multi-component multiphase LBM with a Multi-Relaxation Time (MRT) collision operator to investigate bubble breakup in T-shaped microchannels for high density ratios. The findings of this numerical simulation indicate that there are two types of bubble breakup regimes: one due to bubble– wall contact and highlight the significant role of liquid viscosity in accelerating bubble breakup beyond a particular threshold. Chen et al. [12] used a three-dimensional Shan-Chen multi-component LBM to simulate the coalescence of compound droplets with unequal-sized inner cores. Numerical findings of this simulation indicate that there exist discrete processes of liquid bridge growth, relaxation, and breakup/re-coalescence, with inertial scaling of Rb ∝ t1/2 for the outer liquid bridge. Jing et al. [13] used a modified phase-field LBM coupled with the smoothed boundary method to investigate bubble dynamics in a suspension of obstacles and dendrites in a viscous electrolyte solution. Numerical results from this simulation indicate that wetting of bubble surfaces reduces bubble–obstacle contact, enabling low-deformation sliding. In contrast, wetting of bubble surfaces increases adhesion, resulting in a shear-induced necking instability.

Wan et al. [14] used a phase-field LBM based on the Allen-Cahn equation to simulate bubble dynamics in microchannels with multi-hole orifice plates. Numerical findings of this simulation indicate that there exist two critical Weber numbers, beyond which three bubble dynamics regimes exist, and that increasing contact angles substantially increases bubble passage time, particularly for high Weber number cases. Numerical findings of a numerical simulation by Li et al [15] of bubble dynamics under forced flow in microchannels using a phase-field LBM indicate that bubbles can be detached from a channel wall for high Reynolds numbers, small initial contact angles, and hydrophilic channel walls. Furthermore, bubbles can be split from a bubble resting on a superhydrophilic surface at high Reynolds numbers. Research on multiscale bubble breakup dynamics adjacent to a blade in unsteady turbulence was conducted by Ji et al. [16] using large eddy simulation coupled with the VOF method. Their findings indicated that the dispersion behavior of bubbles in high-velocity-gradient regions surrounding the blade is more pronounced than their ascent behavior, with the bubble fragmentation process governed by the Weber number.

Ding et al. [17] used the LBM coupled with a free-surface approach to study the mechanism of oxygen bubble transport in the porous transport layers and the flow channels of proton exchange membrane water electrolyzers. The results show that hydrophilic porous transport layers favor upward bubble motion over an optimal range of contact angles. At the same time, fibrous structures are more favorable at high current densities than spherical structures.

A mixed Shan-Chen thermal approach for two-phase Rayleigh-Benard convective heat transfer problems is proposed by Al-Zahiwat et al. [18]. Numerical calculations at relatively high Rayleigh numbers (105) show that the convective heat transfer between the hot and cold walls is enhanced by increasing the Rayleigh number. At the same time, the variation of the interface is dependent on the thermal expansion coefficient and the temperature difference.

In the study by Ho and Vu [19], the authors performed numerical calculations of the head-on collision of equal-sized compound droplets, while Zhao et al. [20] extended the calculations to unequal-sized compound droplets. These studies have given detailed insights into the coalescence characteristics and the deformation of compound droplets under different impact conditions.

Despite the achievements, the influence of the boundary on bubble dynamics within the channels remains poorly understood. For example, the interaction between the rising bubble in the channel and the boundary is not well understood, especially in microfluidic systems, enhanced oil recovery, and biomedical applications. Although some studies on the bubble dynamics inside the microchannels with different geometries have been carried out recently [11], [13], [15], the specific influence of the channel walls on the bubble dynamics, such as the trajectory, shape, and breakup inside quiescent fluids, is not yet investigated.

The present study addresses this gap by numerically investigating the buoyancy-driven ascent of a single bubble within a quiescent fluid confined by channel walls. Simulations are performed using the Shan-Chen multiphase LBM framework. A systematic analysis is conducted to evaluate the effects of wall proximity and Eötvös number on bubble deformation, trajectory, and breakup. To the authors' best knowledge, this is the first study to systematically examine the influence of channel walls on bubble dynamics using the LBM approach.

The multicomponent-multiphase interparticle potential model was proposed by Shan and Chen in 1993 [6] based on the BGK collision model [21]. In this model, two components (each fluid type corresponds to one component) interact via a force law based on the interparticle potential. The discrete lattice Boltzmann equation for each component is as follows [6]:

σ σ 1 σ σ(eq) σ σ

fi (x + eiδt, t + δt) − fi (x, t) = − σ (fi (x, t) − fi (n , ueq)). (1)

The term τσ is the relaxation time parameter related to component σ, which is related to the kinematic viscosity as follows:

. (2)

νσ is the kinematic viscosity of component σ and Cs = R ̄ T = 3 The equilibrium distribution function of each component is obtained as follows:

σ σ σ 2

σ(eq) σ ei. ueq ei. ueq 2 ueq

fi = win [1 + 2 + ( 4 ) − 2 ]. (3)

cs 2cs 2cs

There are weighting coefficients, which are defined as follows:

In this model, ci represents the lattice velocity in direction i and is expressed as follows:

The number density (nσ) and the velocity value of component σ (uσ) at location x and time t are obtained, respectively, as follows:

nσ(x, t) = ∑ fiσ(x, t) ,

i (6)

nσuσ(x, t) = ∑ eifiσ(x, t).

The total interparticle force (forces between two fluids) that is applied to particles of component σ at location x is obtained as follows [22]:

Ffσ(x) = Fσ(x) = −ψσ(x) ∑ Gσσ̄ ∑ ψσ̄ (x + eiδt)ei. (7)

σ,σ̄ i

Ffσ represents the total force exerted on component σ at location x from fluid particles at neighboring points (interparticle force between two fluids). After the collision step, the new momentum of component σ at each lattice node is modified as follows [22]:

ρσuσ eq = ρσu′ + τσFσ. (8)

uσeq is the equilibrium velocity of component σ, and Fσ is the total force exerted on component σ. u' is the composite macroscopic velocity, which is very important for analyzing the overall flow. To conserve momentum at each collision step, u' is obtained from the following relation:

′ ∑σ ρσuσ/τσ

u = σ σ . (9)

∑σ ρ /τ The total density, macroscopic momentum, and total pressure at each point of the flow are obtained using the following relation:

ρ(x, t) = ∑ mσ ∑ fiσ(x, t) ,

σ σ 1 σ

ρ(x, t)U(x, t) = ∑ m ∑ eifi (x, t) + 2 ∑ F (x, t) , (10)

σ i σ

R ̄ T σ σ̄ P(x, t) = R ̄ T ∑ mσnσ + ∑ Gσσ̄ ψ ψ .

σ σσ̄

For the total effective force on component σ, we can write:

(11)

In this equation, Fσ is the total force applied to component σ, is the total force of the interactions of neighboring points between fluids that affects component σ, is the gravitational force, and is the force that is exerted from the solid surface on component σ, which, in the present paper, due to the negligible effects of the walls on the droplet behavior, is neglected.

In the present simulation, to present the results in a standardized manner, the following dimensionless numbers have been used:

The Eötvös number (Eo), which represents the ratio of gravitational force to surface tension force:

(13)

The Morton number (Mo), which actually describes the properties of the fluid surrounding the bubble:

(14)

Dimensionless time:

t = 1,

D 2 (15)

where g is the gravitational acceleration, and are the density of the bubble and the density of the fluid surrounding the bubble, respectively, D is the initial diameter of the bubble, is the surface tension between the two fluids, and are the kinematic viscosity of the bubble and the fluid surrounding the bubble, respectively, and is the lattice time step. In all simulations presented in this paper, the Mo is very small compared to the Eötvös number; therefore, the main determinant of the flow regime is the Eötvös number, and variations in the Mo do not significantly affect the bubble's dynamic behavior.

The dynamic behavior of a buoyancy-driven bubble rising in a confined channel has been numerically simulated using the Shan-Chen LBM. The important findings of the simulation are as follows:

I. The bubble deformation is found to increase with the Eo number: at low Eo numbers, the bubble is slightly flattened into an ellipse; at higher Eo numbers, it becomes a disk; and at very high Eo numbers, it breaks up .

II. The flow pattern also changes with the Eo number: at low Eo numbers, the flow is laminar with trapped vortices; at higher Eo numbers, wake vortices and turbulent flow outside the wake occur due to the lack of viscous damping .

III. The proximity of the bubble to the channel wall also leads to an oscillatory path, with bubbles released near the wall exhibiting oscillations towards and away from the wall, and the oscillation amplitude increases with the Eo number .

IV. There is also a critical Eo beyond which the bubble breaks up, irrespective of proximity to the wall, with a range of 30-40 .

V. The influence of the Mo is found to be negligible, as it is much smaller than the Eo number in all cases . The Shan-Chen LBM has shown promise for simulating the complex behavior of buoyancy-driven bubbles in confined channels, and the results will be useful for developing microfluidic devices.

The authors declare no conflict of interest.

All data are included in the text.

This research received no specific grant from funding agencies in the public, commercial, or not-for-profit sectors.