# CFD simulation of isolated bubbles rising in Newtonian or non-Newtonian fluids inside a thin-gap bubble column

Sikandar Almani, Walid Blel, Emilie Gadoin, Caroline Gentric

► To cite this version:

Sikandar Almani, Walid Blel, Emilie Gadoin, Caroline Gentric. CFD simulation of isolated bubbles rising in Newtonian or non-Newtonian fluids inside a thin-gap bubble column. Chemical Engineering Research and Design, 2025, 214, pp.202-218. 10.1016/j.cherd.2024.12.030. hal-04903897

HAL Id: hal-04903897

https://nantes-universite.hal.science/hal-04903897v1

Submitted on 21 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# CFD simulation of isolated bubbles rising in Newtonian or non-Newtonian fluids inside a thin-gap bubble column

Sikandar Almani<sup>a,b</sup> <sup>●</sup>, Walid Blel<sup>a</sup>, Emilie Gadoin<sup>a</sup>, Caroline Gentric<sup>a,*</sup> <sup>●</sup>

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

**Keywords:**

Computational Fluid Dynamics (CFD)

Volume of fluid (VOF)

Isolated bubble

Thin gap bubble column

Shear-thinning liquid

## ABSTRACT

Hydrodynamics of the bubbling process can be complex especially in thin bubble columns, when the gap has the same order of magnitude as the bubble diameter, and with complex fluids. It is then important to understand this phenomenon either by experimental investigation through optical methods such as shadowgraphy and/or Particle Image Velocimetry (PIV) or numerically by Computational Fluid Dynamics (CFD), which, when validated, can allow numerical experimentation in situations which are expensive to implement experimentally or time consuming. In this study, three-dimensional numerical simulations of isolated bubbles rising in Newtonian (water) or non-Newtonian (CarboxyMethyl Cellulose (CMC) and Xanthan Gum (XG) solutions) liquid phases mimicking *Chlorella vulgaris* cultures at 42 g.L⁻¹ concentration inside a 4 mm gap bubble column are performed using the volume of fluid (VOF) model with the ANSYS FLUENT 17.2 code. Results are validated by comparison with shadowgraphy experiments. Bubble terminal velocity, shape, and trajectory are numerically analysed. Wall shear stress (WSS) induced by the bubble, strain rate, viscosity and flow field around the bubble are also discussed. Numerical results show similar trends as experimental ones despite slightly lower terminal velocity and aspect ratio values are observed in comparison to the experimental results. The trajectory of the bubble is non-rectilinear for water and rectilinear for non-Newtonian fluids as observed experimentally. This numerical study highlights the bubble-liquid and bubble-wall interactions that will help to understand the complex phenomena of bubble rise in non-Newtonian media/microalgae suspensions at high concentrations at the local level in thin-gap bubble columns.

## 1. Introduction

Gas-liquid interactions have a critical role when it comes to the efficient performance of bubble columns, airlift reactors, spray towers, falling film columns and packed columns (de Jesus et al., 2017; Fourati et al., 2013; Li et al., 2019; Thobie et al., 2017). In bubble columns, interactions between bubbles and the liquid phase are responsible for the efficient mixing and gas-liquid mass transfer (Khamadieva and Böhm, 2006; Thobie et al., 2017), which strongly depend on the rheological behaviour of the liquid phase. In many chemical, biochemical, pharmaceutical and wastewater treatment plants, the rheological behaviour of the liquid phase changes from Newtonian to

non-Newtonian. For instance, in photobioreactors, a microalgae culture of *Chlorella vulgaris* starts with a Newtonian behaviour (like water) and changes its behaviour to the non-Newtonian (shear thinning) as biomass concentration increases more than 30 g.L⁻¹ (Soulies et al., 2013). Besides that, for achieving high volumetric productivities it is necessary to increase the specific illuminated area that can be achieved by reducing the thickness of the photobioreactors (Cornet and Dussap, 2009). This need for confinement along with rheological modifications makes the gas-liquid interactions more complex. In order to understand this complex phenomena at a global scale, where the rheology and the geometry modifications can alter the gas holdup (Ferrario et al., 2025), mixing time and gas-liquid mass transfer, it is necessary in a first step to understand the rise of a single bubble in such environments and conditions.

**Abbreviations:** CFD, Computational Fluid Dynamics; CLSVOF, Coupled Level Set Volume Of Fluid; CMC, CarboxyMethyl Cellulose; CSF, Continuum Surface Force; DNS, Direct Numerical Simulation; E-E, Euler- Euler; E-L, Euler- Lagrange; HEC, Hydroxyl-Ethyl Cellulose; LFRM, Local Front Reconstruction Method; LS, Level Set; PBM, Population-Based Model; PBR, Photobioreactor; PISO, Pressure-Implicit with the Splitting of Operators; PIV, Particle Image Velocimetry; PLIC, Piecewise-Linear-Interface-Calculation; VOF, Volume Of Fluid; WSS, Wall Shear Stress; XG, Xanthan Gum.

<footnote><sup>a</sup> Nantes Université, Oniris, CNRS, GEPEA, UMR 6144, Saint-Nazaire 44600, France</footnote>

<footnote><sup>b</sup> Process, Simulation & Modelling Research Group, Department of Chemical Engineering, Mehran University of Engineering and Technology, Jamshoro, Pakistan</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: caroline.gentric@univ-nantes.fr (C. Gentric).</footnote>
Nomenclature

d<sub>b</sub> (m), bubble diameter

D<sub>c</sub> (m), column diameter

d<sub>cap</sub> (m), capillary diameter

E (-), aspect ratio

e (m), gap width

f<sub>σ</sub> (N.m<sup>-3</sup>), surface tension source term.

K (Pa.s), power law consistency index

g (m.s<sup>-2</sup>), gravitational acceleration

n (-), power law flow index

p (Pa), pressure

u (m.s<sup>-1</sup>), flow velocity

V<sub>t</sub> (m.s<sup>-1</sup>), terminal velocity

V<sub>ins</sub> (m.s<sup>-1</sup>), instantaneous velocity

Re (-), Reynolds number

S (s<sup>-1</sup>), strain rate tensor

t (s), time

T (°C), temperature

Greek letters

α (-), gas volume fraction

γ (s<sup>-1</sup>), shear rate

κ (m<sup>-1</sup>), curvature of the interface

μ<sub>app</sub> (Pa.s), apparent viscosity

μ<sub>l</sub> (Pa.s), viscosity of liquid phase

μ<sub>g</sub> (Pa.s), viscosity of gas phase

ρ<sub>l</sub> (kg.m<sup>-3</sup>), density of the liquid phase

ρ<sub>g</sub> (kg.m<sup>-3</sup>), density of the gas phase

σ (N.m<sup>-1</sup>), surface tension

Hydrodynamics of bubbling process can be fundamentally understood by experimental investigation through shadowgraphy/Image analysis (Almani et al., 2021; Xu et al., 2024) and particle image velocimetry (PIV) (Thobie, 2018; Thobie et al., 2022a) and numerically by Computational Fluid Dynamics (CFD) (Ahmad et al., 2024; Ghanavati et al., 2023; Mahmoudi et al., 2022; Mei and Cheng, 2024; Rashid et al., 2023; Vaishnavi et al., 2023; Yu et al., 2025), which has proven to be an alternative tool to understand complex phenomena in gas-liquid flows.

Various techniques have been developed over the last few decades to simulate complex two-phase flow. For the analysis of multiphase flows, the Eulerian-Eulerian (E-E), Eulerian-Lagrangian (E-L), and Direct Numerical Simulation (DNS) methods are the three main approaches that are frequently used in CFD. The DNS approach is preferred over others when examining the bubble formation and its ascent in bubble columns since it allows to track the bubble interface. Numerical studies of single and multiple bubbles rising in unconfined domains have used either Volume of Fluid (VOF) method (Hirt and Nichols, 1981) or Level Set (LS) method (Osher and Sethian, 1988) or even Coupled LS and VOF (CLSVOF) method (Sussman and Puckett, 2000). Currently, the VOF method is the most employed one and it is implemented in most commercial CFD codes. It is used along with Continuum Surface Force model (CSF) (Brackbill et al., 1992) to treat surface tension.

Using the VOF model, the rise of single bubbles in water has been investigated to put into evidence the effects of the geometric dimensions or wall effects (Asim et al., 2023; Islam et al., 2015; Krishna et al., 1999; Krishna and Van Baten, 1999; Kumar and Vanka, 2015) or of physico-chemical properties of the liquid phase (Ma et al., 2012; Vaishnavi et al., 2023; Yujie et al., 2012) and bubble formation has also been investigated numerically (Vaishnavi et al., 2023; Zahedi et al., 2014).

The bubble column can be assumed infinite/unconfined, if the ratio between bubble diameter (d<sub>b</sub>) to column diameter (D<sub>c</sub>) becomes smaller than the value of 0.125 (Krishna et al., 1999). Similarly, van Sint Annaland et al. (2005) concluded that no effect will be observed if the domain has lateral dimensions 3–4 times the initial bubble diameter. For the case of **infinite liquid medium**, Krishna et al. (1999) reproduced numerically that the terminal velocity increases with bubble diameter and bubble shape changes from an elliptical to a spherical cap as the diameter increases. Similarly, Bertola et al. (2004) while investigating the terminal velocities, shapes and trajectories of different bubbles (d<sub>b</sub> = 2–10 mm) using the VOF approach observed similar trends as the experimental observations of Clift et al. (1978). Besides water, numerical studies are also reported for glycerol solutions. For instance, Kong et al. (2019) performed an experimental and numerical study on the oscillation dynamics of a 3 mm diameter bubble rising in glycerol solutions. They observed the effect of liquid phase viscosity on bubble

ascent and bubble shape.

Besides infinite medium, very few numerical studies are performed for **confined geometries**. For instance, Keshavarzi et al. (2014) compared VOF and CLSVOF models for the transient rise of 5 and 9 mm bubbles in a confined domain (20 × 200 mm²) having a 5 mm gap. The numerical results are compared to experimental ones obtained by high-speed camera images for the same experimental setup. Another VOF numerical study by Wang et al. (2015) was performed followed by experiments on confined domains (Hele-Shaw cells) having gap thicknesses of 0.5 mm and 1 mm for air-water system representing high Reynolds number flow (Wang et al., 2014). Results suggested that bubble shape has two regimes depending on the bubble diameter: oblate ellipsoid regime is observed for d<sub>b</sub> < 10.5 mm and a spherical shape after d<sub>b</sub> > 12 mm, this observation is similar to the experimental conclusions. Motivated by the goal to clean flat sheet membrane surfaces with bubbling, Böhm et al. (2016b) performed simulations for a single bubble rising in flat rectangular geometry (3 mm gap) and compared them with the experiments conducted in such configurations. From the same research group, Drews et al. (2010) worked on the numerical simulations to highlight the effect of bubble movement and the shear produced by bubbles to clean the membrane walls. The authors pointed out that numerical simulations show shear stresses that are somewhat overestimated compared to the experimental results. Prieske et al. (2010), investigated single bubble ascent using VOF model and pointed out that the highest shear stresses are observed in the smallest gap channels. Gumulya et al. (2021) simulated bubble dynamics (d<sub>b</sub> = 5–7 mm) at high Reynolds numbers (Re ≈1000–2000) in pure water in a rectangular column. The size of the gap is varied based on the diameter of the bubble in order to investigate the impact of confinement ratio (e/d<sub>b</sub> = 1–2) on the velocity field around the bubble. Numerical results show that the bubble attains zigzag rising behaviour in confined conditions in comparison to spiral rise observed in free conditions. The zigzag tendency and the shape deformation decrease as the level of confinement increases.

Numerical studies on non-Newtonian fluids are less reported in the literature in comparison to Newtonian ones. Bubble dynamics in **non-Newtonian fluids** (especially shear-thinning) have been numerically modelled by using power law or Carreau-Yasuda rheological models. Some numerical studies deal with activated sludge (Radaei et al., 2018), viscous sugar solutions (Hassan et al., 2007) or fermentation media (Zhang et al., 2010) which are mimicked using non-Newtonian fluids. Premlata et al. (2017a, 2017b) studied the behaviour of bubbles rising in fluids described by Carreau-Yasuda model using VOF method. Working with water and CMC solutions they concluded that with the increase in shear-thinning tendency (lowering “n” value) in an infinite medium, the rise velocity of the bubble increases and the shape becomes more stable.
In another study, Zhang et al. (2010) modelled the shear-thinning fluids (CMC, XG and Hydroxyl-Ethyl Cellulose (HEC)) with Carreau rheological model using the level set approach and investigated velocity and viscosity distribution around a single bubble. It has been concluded that the local change in viscosity around the bubble depends on the bubble shape and the zero-shear viscosity whereas the shear rate and velocity fields are responsible for the formation of a stagnant and high viscosity region at the rear of the bubble. Concerning the terminal velocity and shape in non-Newtonian fluids, Battistella et al. (2020) investigated large range of viscosities covering both shear-thinning and shear-thickening behaviours (0.5 ≤ n ≤ 1.5) using power law. The simulations show that shear-thickening fluid tries to maintain bubble's spherical shape, whereas, in shear-thinning fluids, the small spherical bubbles lose their sphericity. They also formulated the drag relation with a modified Reynolds number that can predict the terminal velocity with a 20% deviation in non-Newtonian fluids in the range of power-law index (n) studied. Islam et al. (2020) performed simulations on Newtonian and non-Newtonian fluids using VOF coupled with CSF formulations. The influence of dimensionless numbers on bubble rise were characterized and it has been observed that by decreasing Morton number, terminal bubble velocity increases and bubble possesses a non-rectilinear/zigzag trajectory in low Morton number for Newtonian fluids. Whereas, for the non-Newtonian fluids, terminal velocity increases with the increase in the rheological index and shape changes from ellipsoidal to spherical cap with a rectilinear trajectory.

Several numerical investigations have been done for confined configurations with Newtonian liquid phase and some with non-Newtonian fluids in infinite geometries. To the best of the authors' knowledge, very less studies are reported for the confined gap configurations dealing with a non-Newtonian liquid phase. The present work aims to investigate numerically the effect of confinement on the shape and hydrodynamics of isolated bubbles rising in a thin-gap bubble column (gap thickness e = 4 mm) operating with Newtonian or non-Newtonian fluids. Water is used as a low viscosity Newtonian fluid and aqueous solutions of CMC and XG as non-Newtonian ones. Rheological measurements have been performed to select the concentrations of model fluids (CMC and XG) that mimic the rheology of the microalgae, *Chlorella vulgaris*, at high concentration (42 g.L⁻¹) which presents a shear-thinning behaviour (Souliès et al., 2013). Experimentally, the terminal velocity (Vt), the equivalent bubble diameter (db), the bubble shape and the trajectory of bubbles generated from five different capillary diameters (dcap) in three fluids have been characterized by shadowgraphy (Almani et al., 2021). Numerically, simulations are performed using ANSYS FLUENT 17.2 academic version. The VOF method for interface tracking, along with the Continuum Surface Force (CSF) model to treat surface tension is used to simulate the behaviour of bubbles with different diameters in a thin-gap bubble column for each tested fluid. The instantaneous velocity and trajectory of the bubble are determined by tracking its center of gravity. Numerical results are compared to experimental ones to validate the numerical simulation tool, first in an infinite domain and then in this thin-gap bubble column. Results are also analysed for a better understanding of gas/liquid or fluid/wall interactions (i.e. velocity, strain rate and viscosity fields around the bubble are discussed).

Thus, this study is one on the first presenting direct numerical simulation of bubble rising in confined in non-Newtonian solutions. Experimental validation is carried out for terminal velocity, bubble shape, trajectory and shows reasonable agreement. This allows to investigate numerically parameters which have not been measured such as wall shear stress, flow in the bubble wake which can present operational interest in processes such as photobioreactors, or flat membrane modules. After validation, this numerical approach can be used for numerical experiments i.e. investigation of other liquid phases, different gap thicknesses and geometries.

## 2. Computational Fluid Dynamics (CFD) model and fluid properties

### 2.1. Governing equations

#### 2.1.1. Conservation equations

In the VOF method (Hirt and Nichols, 1981), both phases are treated using the same set of momentum and continuity equations:

$$ \text{Continuity} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 \qquad (1) $$

$$ \text{Momentum} \frac{\partial}{\partial t}(\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \mathbf{u}) = - \nabla p + \nabla \cdot \mu (\nabla \mathbf{u} + \nabla \mathbf{u}^T) + \rho \mathbf{g} + \mathbf{f}_{\sigma} \quad (2) $$

The solution of the volume fraction equation is used to track the interface between the phases:

$$ \text{Volume Fraction} \frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = 0 \qquad (3) $$

In the above equations, u, p, t and α are fluid velocity, pressure, time and gas volume fraction in the computational cell. g is the acceleration due to gravity, ρ the average density, μ the average dynamic viscosity and fσ is an additional surface tension source term. The value α = 0 corresponds to the cell being filled completely with liquid whereas α = 1 indicates that the cell is filled with gas. For a cell containing the interface, the value of α is between 0 and 1. ρ and μ are defined in each phase and in cells containing the interface, they can be estimated by:

$$ \rho = \alpha \rho_g + (1 - \alpha) \rho_l \qquad (4) $$

$$ \mu = \alpha \mu_g + (1 - \alpha) \mu_l \qquad (5) $$

The surface tension force along the gas bubble-water interface is calculated using the Continuum Surface Force (CSF) model (Brackbill et al., 1992):

$$ \mathbf{f}_{\sigma} = \sigma \rho \kappa \nabla \alpha \qquad (6) $$

where σ and κ are the surface tension and curvature of the interface respectively. κ is defined in terms of the divergence of the unit normal.

$$ \kappa = \nabla \cdot (\nabla \alpha / |\nabla \alpha|) \qquad (7) $$

where ∇α is the gradient of the volume fraction α and |∇α| its norm.

#### 2.1.2. Viscosity model for non-Newtonian liquids

For incompressible Newtonian fluids, the viscosity μl is constant. For non-Newtonian fluids, the viscosity varies with the shear rate. The power law is adopted here to represent the shear-thinning behaviour of the CMC and XG solutions:

$$ \mu_l = K \dot{\gamma}^{n-1} \qquad (8) $$

where K is the consistency index, n is the power law index and the shear rate γ can be written as a function of the magnitude of the strain rate tensor (S):

$$ \dot{\gamma} = \sqrt{\frac{2}{S(S: S)}} \qquad (9) $$

and the strain rate tensor is defined as:

$$ S = \frac{1}{2} (\nabla \mathbf{u} + \nabla \mathbf{u}^T) \qquad (10) $$

### 2.2. Confined bubble column mesh

The simulation domain corresponding to the confined conditions has dimensions of 4 × 40 × 100 mm³. The mesh size is 0.2 mm at the centre of the geometry to have an accurate interface identification, and is
reduced near the walls to 0.05 mm to study the shear generated by bubbles at walls. This mesh size is chosen to ensure that the bubble should have the sufficient number of grid cells required to fully capture the dynamic motion of the gas-liquid interface. According to Bertola et al. (2004), there is a negligible effect on results after maintaining 300 cells/bubble. For this reason, the mesh size was even reduced to 0.1 mm for the smallest bubbles having a diameter lower than 2 mm.

Five different equivalent bubble diameters, identified experimentally in the three different fluid environments, are simulated in this confined column. The numerical terminal velocities, the bubble shape and the trajectory are then compared with the experimental results of Almani et al. (2021).

## 2.3. Numerical simulation strategies and boundary conditions

The Computational Fluid Dynamics (CFD) software, ANSYS FLUENT 17.2, is used to simulate the ascent of an isolated air bubble in a confined medium with a Newtonian (water) or non-Newtonian (CMC and XG solutions) liquid phase. The physico-chemical and rheological properties of liquids are shown in Table 1. The computational domain is shown in Fig. 1a. The boundary conditions assigned at the bottom and the front and rear walls of the 4 mm gap column are “walls”, the lateral sides of the column are “symmetries” and the bubble column exit is “pressure outlet”. No-slip conditions are applied at the walls with the wettability contact angle of 180° corresponding to hydrophobic behaviour. The operating pressure is equal to the atmospheric pressure, i.e., 101325 Pa and along the Z-direction the gravitational force (g) of -9.81 m.s⁻² is assigned.

The simulation case is run under double precision with a transient, pressure-based solver. The Volume Of Fluid (VOF) model is used to perform interface tracking along with the Continuum Surface Force model (CSF) to simulate the surface tension.

For non-Newtonian solutions, the viscosity is modelled by the non-Newtonian power-law model with minimum and maximum viscosity limits (0.001 Pa.s and 10 Pa.s respectively) and values of consistency index (K) and power-law index (n) obtained from rheological measurements. Bubbles with diameters corresponding to those obtained experimentally (Almani et al., 2021) are initially positioned at a height of 10 mm from the bottom centre of the column by imposing an air volume fraction of 1.0 in a spherical region with the diameter of the simulated bubble. Pressure-velocity coupling is solved using the pressure-implicit with split operator (PISO) method, the pressure discretization method is pressure staggering scheme (PRESTO), and the discretization of momentum equations is performed with the QUICK scheme. An explicit geometric reconstruction scheme that is based on the PLIC (Piecewise-Linear-Interface-Calculation) method is used to reconstruct interface position. The time step is 0.0005 s for large bubble diameters and was reduced to 0.0001 s for the smallest bubble having d<sub>b</sub> less than 2 mm.

<table><caption>**Table 1**
Physico-chemical and rheological properties of Newtonian and non-Newtonian fluids.</caption><thead><tr><th>Fluids</th><th>ρ (kg.m<sup>-3</sup>)</th><th>σ (N.m<sup>-1</sup>)</th><th>K (Pa.s<sup>n</sup>)</th><th>n</th><th>T (°C)</th></tr></thead><tbody><tr><td>Demineralized water</td><td>1000</td><td>0.0728 ± 0.002</td><td>0.0010</td><td>1.00</td><td>22</td></tr><tr><td>CMC solution (2 g.L<sup>-1</sup>) (Weak Shear-thinning power law model) Shear rate range from 1 to 1000 s<sup>-1</sup></td><td>1001.3</td><td>0.0674 ± 0.001</td><td>0.0236</td><td>0.87</td><td>22</td></tr><tr><td>CMC solution (2 g.L<sup>-1</sup>) (Strong Shear-thinning power law model) Shear rate range from 0.01 to 1000 s<sup>-1</sup></td><td>1001.3</td><td>0.0674 ± 0.001</td><td>0.0608</td><td>0.65</td><td>22</td></tr><tr><td>XG solution (1 g.L<sup>-1</sup>)</td><td>998.7</td><td>0.0519 ± 0.002</td><td>0.2000</td><td>0.39</td><td>22</td></tr><tr><td>Chlorella vulgaris suspension (42 g.L<sup>-1</sup>)</td><td>-</td><td>-</td><td>0.0561</td><td>0.31</td><td>22</td></tr></tbody></table>

The real-time for the simulation is 4–5 days for each bubble using the parallel processing mode on a Dell i5-7600 with 16 GB of RAM using four parallel Intel (R) Core (TM) processors each 3.50 GHz.

## 2.4. Physico-chemical and rheological properties of Newtonian and non-Newtonian fluids

The aim of this work is to simulate single bubble ascent in water and in two non Newtonian solutions mimicking a *Chlorella vulgaris* suspension as it was performed experimentally by Almani et al. (2021). For experimental work, the rheological characterization of non Newtonian fluids was performed at 22°C using a rotational stress-controlled rheometer (MCR500, PAAR Physica®) equipped with a cone-plate device (50 mm in diameter, 3 degrees in cone angle). The shear rate was varied in the range of 0.01–1000 s⁻¹. The power-law model (Eq. 8) was used to represent the shear-thinning behaviour of the CMC and XG solutions. Fig. 2a. presents the rheograms of the three fluids along with their power law fit. The values of physico-chemical and rheological properties of Newtonian and non-Newtonian fluids are presented in Table 1. It can be noticed from Fig. 2b that two values of K and n can be used for CMC rheogram description depending on the range of shear rate where they are estimated: whole shear rate domain or 1–1000 s⁻¹. The influence of the choice between these values for numerical simulations will be discussed further in the paper.

## 2.5. Determination of the numerical terminal velocity

The instantaneous velocity of the rising bubble is obtained by calculating the difference in the vertical position of its gravity centre between two consecutive time steps divided by the time interval between them. Fig. 3 shows the instantaneous bubble rising velocity vs flow time, for a 3.36 mm bubble diameter in XG solution. Blue points refer to the velocity estimated at each time step and the red line is the cumulative average during 0.2 seconds. For all cases, the bubble is initially motionless (initial condition), and the bubble rising velocity increases during 0.1–0.2 seconds and then there is no more variation. The asymptotic value thus corresponds to the “numerical” terminal velocity.

## 2.6. Domain size and mesh independency test

For numerical approach validation, the case of a rising bubble in a 3D infinite liquid medium for which experimental and numerical results are available is first considered. The first set of simulations was carried out in 2D to determine the size of the computational domain to ensure that it is sufficiently large to correspond to the condition of an infinite quiescent liquid. For that, the width was varied from 20 mm to 50 mm, whereas the height of the rectangular domain was fixed at 100 mm to leave enough height for the terminal velocity (V<sub>t</sub>) to be reached (V<sub>t</sub> is attained up to 30–40 mm for the bubble diameter of 3 mm according to Islam et al. (2015). It was also reported in the literature that the computed terminal velocities were not affected by the computational domain size, when the lateral dimensions in the horizontal directions (i. e. the x- and y-directions) exceed 3–4 times the bubble diameter (van Sint Annaland et al., 2005). In our numerical simulations, no effect on terminal velocity values was observed from the width of 30 mm for a bubble having a 3 mm diameter (which corresponds to a distance between the bubble and each lateral wall equal to 4 times the bubble diameter).

Then, the objective of the mesh independency test is to be sure that the chosen mesh size does not significantly impact the results, keeping in mind the computational cost for the simulation. Besides that, the bubble must be defined by sufficient computational cells number to represent accurately the bubble shape.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Confined medium computational domain: (a) Isometric view and side view, (b) Mesh refinement structure near the wall.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Rheograms and identified power law models: a. for all fluids, b. for the CMC solution.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Determination of the numerical terminal velocity of the bubble ($d_b$ =3.36 mm in XG solution).</figcaption></figure>

For the mesh independency test, a 3D infinite medium numerical domain was created having 40 mm width, 40 mm length and 100 mm height. The boundary condition for lateral walls was set to be “symmetry”, the bottom wall was considered as “wall” and the top wall was set as “pressure outlet”. Four mesh sizes have been tested in this 3D infinite medium numerical domain, the center zone of this

parallelepiped (15 × 15 × 100 mm³) having a uniform mesh size of 0.2, 0.25, 0.3 or 0.5 mm. These meshes correspond to 1766, 904, 523 and 113 grid cells/bubble respectively. Outside this central zone, a smooth transition takes place until the mesh size reaches 0.5 mm at the domain end. Then, simulations of the ascent of a 3 mm diameter bubble are carried with these meshes. Initial condition consists in “patching” a 3 mm spherical air bubble at a height of 10 mm. The vertical position of the centre of gravity (CG$_z$) of the bubble is tracked and the instantaneous bubble rising velocities are compared for the 4 meshes in Fig. 4. The results show that 0.2 and 0.25 mm mesh sizes have similar instantaneous bubble rising velocities, hence the mesh size of 0.25 mm has been chosen as the largest mesh which satisfies grid independence criterion.

### 3. Results and discussion

#### 3.1. Infinite medium bubble simulation for numerical code validation

For the numerical code validation, the ascent of bubbles of four different diameters ($d_b$ = 3, 4, 7 and 10 mm) in an infinite medium was simulated. The obtained numerical results are compared to those of Bertola et al. (2004), who simulated the rise of a single bubble in a cylindrical column having a diameter equal to 5 times the one of the bubble (5$d_b$) and a height of 8 times the diameter of the bubble (8$d_b$) for different $d_b$ (3, 4, 7 and 10 mm) in water. Comparison is also made with the correlation of Clift et al. (1978) for the terminal velocity of bubbles in pure water. A good agreement is obtained between the results of this study and the literature (Fig. 5).

The validation of the current numerical code with the literature shows satisfying results with only the underestimation of numerical terminal velocity for a 10 mm bubble diameter. It can be commented that regarding the experimental study (Almani et al., 2021), the bubble diameter range is between 1.3 and 3.4 mm (estimated by shadowgraphy technique), a range in which the approach is validated.

Once the numerical parameters adopted for the VOF approach are validated, single bubbles rising in water and shear-thinning fluids are simulated. Numerically, the analysed parameters are bubble terminal velocity, bubble shape and bubble trajectory. Besides this, velocity, strain rate and viscosity fields around the bubble, as well as the wall shear stress (WSS) are also discussed.

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Instantaneous bubble rising velocity vs flow time: Mesh independency test for 3 mm bubble diameter in a 3D infinite bubble column.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 5. Model validation with the results of Bertola et al. (Bertola et al., 2004) and the correlation of Clift et al. (1978).</figcaption></figure>

## 3.2. Terminal velocity of bubbles

### 3.2.1. Terminal velocity of bubbles in water

Terminal velocity of isolated bubbles rising in water inside the confined numerical domain is determined. The bubbles considered in this study have a confinement ratio ($d_b/e$) between 0.34 and 0.84, as the ones studied experimentally. The numerical results are compared to the experimental ones obtained by the shadowgraphy technique (Almani et al., 2021).

Fig. 6 shows the experimental and numerical terminal velocities as a function of the diameter for isolated bubbles rising in water. The experimental results observed in a thin-gap bubble column ($e = 4$ mm) have shown a decrease in bubble terminal velocity as a function of bubble diameter for bubbles between 1.74 mm and 3.12 mm and after that, a sudden increase was observed for 3.28 mm bubbles. The numerical results show that the terminal velocity decreases following the same trend as the experimental one up to 2.78 mm diameter and then increases for the other two bubble diameters. Thus, the wall friction effect seems to dominate for the smallest bubbles with a terminal velocity decrease with increasing diameter and then buoyancy effects prevail with an increasing terminal velocity with the bubble size. A similar decrease followed by an increase in terminal velocity as a

<figure><img src="image_8.png" /><figcaption>Fig. 6. Terminal velocity of an isolated bubble in a thin-gap bubble column vs diameter - demineralized water. Comparison between experimental (Almani et al., 2021) and numerical results.</figcaption></figure>

function of bubble diameter was observed experimentally by Böhm et al. (2014) in 5 and 7 mm gap columns. In their study, this increase was observed when the bubble diameter becomes larger or equal to the gap of the column i.e. when the confinement ratio becomes $d_b/e \ge 1$. This trend was also observed by Prieske et al. (2010) during CFD simulations in confined columns (3–7 mm) using the VOF method.

Table 2 presents the experimental and numerical terminal velocities for five different bubble diameters observed in Almani et al. (2021) and for the three fluids (water, CMC and XG).

The maximum relative error between experimental and numerical terminal velocities for water is 11 %. This error is comparable to the ones observed in similar studies in the literature. Indeed, Böhm et al. (2013b) while comparing the experiment and simulation results in 5 mm gap configuration observed a difference of 4 %, 6 % and 16 % for the 3 mm, 5 mm and 7 mm bubble respectively. In another article, Böhm et al. (2014), compared their experimental results with the CFD simulation results of Prieske et al. (2010) with the same configuration and found an average error of 9.4 %. They explained that this error could be due to the half channel depth in simulations and the existence of a symmetry plane between the x-y walls that can suppress the oscillations of the bubbles especially, those having a lower diameter than the column width.

### 3.2.2. Terminal velocity of bubbles in the non-Newtonian liquid phases

Firstly, a sensitivity study has been performed for CMC to investigate the effect of rheological parameters on the terminal velocity of the bubbles. Then results are discussed for both solutions of CMC and XG.

#### 3.2.2.1. Sensitivity study of the rheological parameters of the CMC solution on the numerical terminal velocities.
The experimental rheogram of the CMC solution (2 g.L⁻¹) can be modelled by the power law model as can be seen in Fig. 2b. The power law can be fitted on the whole shear rate domain (0.01–1000 s⁻¹), which gives a power law index n = 0.65 or can be adjusted for high shear rates only (larger than 1 s⁻¹) which gives n = 0.87. The sensitivity study aims to discuss the effect of this rheological parameter on the velocity of a bubble rising in a shear-thinning fluid.

Simulation of a 3.28 mm diameter bubble in the CMC solution has been performed with n = 0.65 to obtain the local strain rate and viscosity values around the bubble in the liquid phase. Fig. 7 presents their values on a horizontal line cutting the bubble at its centre in the X direction (i.e. along the column width). Fig. 7a shows that the strain rate near the bubble interface is as high as 800 s⁻¹. This strain rate decreases towards the lateral boundaries and reaches a minimum strain rate value of 1 s⁻¹ (Fig. 7a). On the other hand, the viscosity of the liquid near the bubble is observed low and increases towards the lateral walls due to lower shear rate values observed far from the bubble (Fig. 7b). It is therefore concluded that the bubble can be simulated in this CMC solution with values of K and n corresponding to the high shear rate part of the rheogram (from 1 to 1000 s⁻¹), i.e. K= 0.0236 and n = 0.87. In fact, this model, which corresponds to a low shear-thinning fluid, allows a better fitting of the apparent viscosity in the shear rate range encountered in the confined column around the bubble.

Simulations of five different bubble diameters in the CMC solution have been performed with both values of K and n, which correspond to weak shear-thinning (n = 0.87) and strong shear-thinning (n = 0.65) solutions. The numerical terminal velocity results show two domains depending on the bubble diameter similar to experimental results. In the first domain, i.e. for small bubble diameters up to 2.2 mm, the terminal velocity increases with the diameter, and for diameters greater than 2.2 mm, a plateau of terminal velocity appears. The plateau of terminal velocity for n = 0.65 is 15 cm.s⁻¹ and for n = 0.87 is 17 cm.s⁻¹ whereas it was 19 cm.s⁻¹ for the experimental results (Fig. 8). The maximum relative error is acceptable for n = 0.87 as it is 11 %, whereas, it reaches 46 % for a strong shear-thinning solution (Table 2). Hence this
<table><caption>**Table 2**
Experimental (Almani et al., 2021) and numerical results of terminal velocity of a single bubble in water, CMC and XG environments.</caption><thead><tr><th>Fluid</th><th>d<sub>cap</sub><br>(mm)</th><th>d<sub>b</sub><br>(mm)</th><th>Confinement ratio<br>(d<sub>b</sub>/e)</th><th>V<sub>t</sub> (m.s<sup>-1</sup>) Experimental<br>Almani et al. 2021</th><th>V<sub>t</sub> (m.s<sup>-1</sup>) Numerical<br>This Study</th><th>Relative Error % |(V<sub>t(exp)</sub>-V<sub>t</sub>/(num))/V<sub>t(exp)</sub>| × 100</th></tr></thead><tbody><tr><td rowspan="5">Water n = 1</td><td>0.127</td><td>1.74</td><td>0.435</td><td>0.1948</td><td>0.2067</td><td>6.1</td></tr><tr><td>0.254</td><td>1.95</td><td>0.487</td><td>0.2014</td><td>0.1795</td><td>10.8</td></tr><tr><td>0.508</td><td>2.78</td><td>0.695</td><td>0.1853</td><td>0.1684</td><td>9.1</td></tr><tr><td>0.762</td><td>3.12</td><td>0.780</td><td>0.1704</td><td>0.1825</td><td>7.1</td></tr><tr><td>1.016</td><td>3.28</td><td>0.820</td><td>0.2030</td><td>0.1865</td><td>8.1</td></tr><tr><td rowspan="5">CMC (2 g.L<sup>-1</sup>) n = 0.87 Weak<br>Shear-thinning</td><td>0.127</td><td>1.58</td><td>0.395</td><td>0.1301</td><td>0.1180</td><td>9.3</td></tr><tr><td>0.254</td><td>2.22</td><td>0.555</td><td>0.1797</td><td>0.1667</td><td>7.2</td></tr><tr><td>0.508</td><td>2.82</td><td>0.705</td><td>0.1876</td><td>0.1692</td><td>9.8</td></tr><tr><td>0.762</td><td>3.16</td><td>0.790</td><td>0.1883</td><td>0.1677</td><td>10.9</td></tr><tr><td>1.016</td><td>3.28</td><td>0.820</td><td>0.1890</td><td>0.1695</td><td>10.3</td></tr><tr><td rowspan="5">CMC (2 g.L<sup>-1</sup>) n = 0.65 Strong<br>Shear-thinning</td><td>0.127</td><td>1.58</td><td>0.395</td><td>0.1301</td><td>0.0700</td><td>46.1</td></tr><tr><td>0.254</td><td>2.22</td><td>0.555</td><td>0.1797</td><td>0.1366</td><td>23.9</td></tr><tr><td>0.508</td><td>2.82</td><td>0.705</td><td>0.1876</td><td>0.1519</td><td>19.0</td></tr><tr><td>0.762</td><td>3.16</td><td>0.790</td><td>0.1883</td><td>0.1538</td><td>18.3</td></tr><tr><td>1.016</td><td>3.28</td><td>0.820</td><td>0.1890</td><td>0.1513</td><td>19.9</td></tr><tr><td rowspan="5">XG (1 g.L<sup>-1</sup>) n = 0.39</td><td>0.127</td><td>1.34</td><td>0.335</td><td>0.0800</td><td>0.0744</td><td>7.0</td></tr><tr><td>0.254</td><td>2.20</td><td>0.550</td><td>0.1587</td><td>0.1414</td><td>10.9</td></tr><tr><td>0.508</td><td>3.10</td><td>0.775</td><td>0.1621</td><td>0.1397</td><td>13.8</td></tr><tr><td>0.762</td><td>3.14</td><td>0.785</td><td>0.1647</td><td>0.1388</td><td>15.7</td></tr><tr><td>1.016</td><td>3.36</td><td>0.840</td><td>0.1594</td><td>0.1361</td><td>14.6</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Numerical strain rate (a) and viscosity profile (b) on a horizontal line cutting the bubble at its centre in the X direction (i.e. along the column width) at 70 mm height (d<sub>b</sub> = 3.28 mm of 2 g.L<sup>-1</sup> CMC (n = 0.65)).</figcaption></figure>

sensitivity study highlights the importance of selecting/measuring the right parameters depending on the local shear rates in the system. Therefore, the results of weak shear-thinning CMC solution n = 0.87 have been chosen for the hereafter numerical results of this study.

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Effect of rheological parameters of the CMC solution on the numerical terminal velocities.</figcaption></figure>

### 3.2.2.2. Terminal velocity of bubbles in the non-Newtonian liquid phase (CMC and XG solutions).
Fig. 9 presents the evolution of both experimental and numerical terminal velocities of bubbles vs diameter in both non-Newtonian fluids: CMC (n = 0.87) and XG (n = 0.39) solutions. The numerical terminal velocities for both shear-thinning fluids show a similar trend as experimental observations, i.e. for small bubble diameters up to 2.2 mm, the terminal velocity increases with the diameter, and for diameters greater than 2.2 mm, a plateau of terminal velocity appears; but the numerical value of the plateau is slightly lower than the experimental one. For instance, numerical plateau values are 14 cm.s<sup>-1</sup> and 17 cm.s<sup>-1</sup> instead of experimental plateau values of 16 cm.s<sup>-1</sup> and 19 cm.s<sup>-1</sup> for XG and CMC respectively. This effect of the rheological parameters on the terminal velocity of bubbles in non-Newtonian fluids was also observed in the VOF study of Islam et al. (2020), who put into evidence a decrease of the bubble terminal velocity with the decrease of the power law index (n).

In our study, the highest relative error for the non-Newtonian fluids is 10.9 % for CMC and 15.7 % for XG as shown in Table 2, whereas, the average relative error is 9.5 % and 12.4 % for CMC and XG respectively. The error is of the same order of magnitude as what has been observed for water, which can validate models used for simulations in non-
<figure><img src="image_11.png" /><figcaption>Fig. 9. Comparison between experimental (Almani et al., 2021) and numerical terminal velocities of isolated bubbles in CMC and XG solutions in a confined column.</figcaption></figure>

Newtonian fluids.

### 3.3. Bubble shape and trajectory

#### 3.3.1. Bubble Shape

The evolution of a bubble shape as a function of increasing diameter changes from a rigid sphere to an ellipsoidal and finally to a spherical cap (Clift et al., 1978). Here, numerical aspect ratio E is obtained from the image processing of the air volume fraction image with Image J software to determine the major and minor axis of the bubble, and E is defined as:

$$E = \frac{b}{a} \qquad (11)$$

where a is the major axis or bubble's length in horizontal direction and b is the minor axis or bubble's length in vertical direction.

The comparison of the numerical and experimental aspect ratios is shown in Fig. 10 (a, b and c) for water, CMC and XG respectively. The experimental results for water show that experimentally the bubble shape changes from spherical to ellipsoidal after $d_b = 3$ mm, where E becomes lower than 0.9 (Fig. 10a). For the shear-thinning XG and CMC solutions, two shape categories are observed using experimental results (Fig. 10b and c): rigid spheres in the creeping flow regime ($d_b < 3$ mm), and oblate ellipsoid shape for bubble diameters equal or greater than 3 mm, corresponding to the plateau zone of the terminal velocity. For the numerical aspect ratio, these two shape categories are also observed but with lower aspect ratio values. For example, the numerical aspect ratio for a CMC bubble with more than 3 mm diameter is 0.7 instead of 0.8 experimentally, whereas, for XG, it has a value of 0.7 in comparison to the experimental one (0.9) (Fig. 10b and c). Drews et al. (2010) simulated bubble diameters from 3 to 10 mm in different flat plate reactors at gap thicknesses of 3 and 7 mm using the VOF method. The results were compared with experimental ones obtained with high-speed camera, which also showed the same difference in the bubble shapes as observed here. In fact, the numerical bubbles had lower aspect ratio than the experimental ones.

#### 3.3.2. Bubble trajectory

The comparison between the experimental and numerical trajectories of the lowest and the highest bubble diameters in all three fluids are shown in Fig. 11 and Fig. 12.

Fig. 11 shows the numerical bubble contours in a plane located at mid-thickness of the column, the lateral view of bubble rise and the top

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 10. Experimental (Almani et al., 2021) and numerical aspect ratio, E, for different bubble diameters (a) Water (b) CMC ($2 \text{ g.L}^{-1}$) and (c) XG ($1 \text{ g.L}^{-1}$).</figcaption></figure>

view showing the bubble-wall interaction in 4 mm gap in the numerical domain. The numerical bubble trajectory in the mid-plane can be compared to the experimental one obtained by shadowgraphy. The top view shows the bubble position at 70 mm height (Z direction) where terminal velocity is attained. The advantage of numerical simulation in comparison with the experimental results is to be able to display lateral and top views to better understand bubble/wall interactions. It can be seen in the top and lateral views that the wall interactions are less for the smallest bubble diameters compared to the largest ones: the smallest bubbles rise in the liquid without colliding with the walls and still experience slight fluctuations whereas the largest bubble trajectory is constrained between the front and back walls since its major axis has
<figure><img src="image_15.png" /><figcaption>**Fig. 11.** Experimental (Almani et al. 2021) and numerical trajectory evaluation of the lowest and the highest bubble diameters in all three fluids. Numerical: $\Delta t = 0.01$ s (small $d_b$) and $\Delta t = 0.025$ s (large $d_b$); Experimental: $\Delta t = 0.067$ s.</figcaption></figure>

about the same size as the gap.

Fig. 12 compares the experimental and numerical positions of the centre of gravity along the X axis (CGx) during bubble ascent and also the numerical trajectories of five bubble diameters rising in three fluids. A good agreement between experimental and numerical results is observed in Fig. 11. and Fig. 12. Experimental bubble trajectories are well reproduced, especially for the fluctuations and inclinations in the case of the largest diameter in water. These results show the rectilinear trajectory behaviour of bubbles in non-Newtonian fluids, whereas in water, bubbles oscillate more. The oscillation along a zig-zag path observed in air-water system leads to periodic inclinations of the ellipsoidal bubbles, whereas in the case of CMC and XG solutions the rectilinear trajectory imposes an oblate ellipsoid shape with a fixed horizontal major axis (Fig. 11).

These bubble oscillations in water increase with the bubble diameter whereas in non-Newtonian fluids bubbles become more ellipsoidal and stable as the diameter increases. For instance, it can be seen in Fig. 11 that the smallest bubbles rising in water and both shear-thinning fluids have the same rectilinear trajectory but the bubble with the largest diameter in case of water has a zig-zag trajectory with a deformable shape due to more oscillations whereas the largest bubbles rising in CMC and XG solutions have rectilinear trajectories along with a stable shape

and negligible oscillations. While comparing the experimental and CFD results for the rising path of 3 and 5 mm bubble diameters in water inside a 5 mm gap column, Böhm et al. (2013b) concluded that there are more oscillations for a bubble of 5 mm diameter in comparison to 3 mm one. In Fig. 12, it can be observed that for the air-water system, the passage from the rectilinear trajectory towards an oscillating one is observed when the bubble shape becomes more ellipsoidal just after the first bubble diameter $d_b = 1.74$ mm, when the numerical aspect ratio E becomes lower than 0.76. For the CMC and the XG solutions, the viscous forces are very important, therefore the trajectory remains rectilinear whatever the diameter of the bubbles.

### 3.4. Velocity vector field around a bubble

Fig. 13 shows the velocity vector field along with streamlines for the smallest and the largest bubble diameters in a 1.6 × 1 cm² image, for the three fluids, when terminal velocity is reached. High velocities vectors are observed around the bubble and in the immediate wake. Wake structure is different between Newtonian (water) and non-Newtonian (CMC and XG) fluids.

The velocity field around the bubble rising in the water (Fig. 13a and b) presents two zones: the first, in which the bubble moves the liquid
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Experimental (Almani et al. 2021) and numerical bubble trajectories for the 5 bubble diameters and the 3 tested fluids. (a) water, (b) CMC (2 g.L⁻¹) and (c) XG (1 g.L⁻¹).</figcaption></figure>

ahead of it in the upstream flow and the liquid from upstream flows back to the wake from both sides of the bubble and the second zone, where bubble oscillations produce regular spiral vortices in the near wake of the bubble. It can be seen that the magnitude of regular spiral vortices depends on the bubble diameter. It is important to mention here that these spiral vortices have high velocities similar to the velocities observed around the bubble which can explain the contribution of wakes in mixing performances as stated by Alméras et al. (2018) for a 1 mm gap bubble column. Due to these instabilities in the central wake, the bubble goes through cycles of acceleration and deceleration hence as a result bubbles in water have zig-zag/non-rectilinear trajectories compared to the ones in non-Newtonian medium. In their PIV study in a 5 mm gap channel, Böhm et al. (2016a) have shown that bubbles with a diameter of 5 mm rising in water follow an oscillating motion with bubble deformations, and Von Karman vortex street in their wake.

For non-Newtonian fluids (Fig. 13 c-f), the flow field has a typical pattern, in which the bubble moves the liquid ahead of it in the upstream flow (zone-I) and liquid from upstream flows back to the wake from both sides of the bubble creating a pseudo-plug flow (zone-II), as stated by Li et al. (2012) for CMC solutions (2–10 g.L⁻¹). According to Li et al. (2012), at low concentrations of CMC, a spiral vortex rope appears in wake of the bubble, whereas as the concentration of CMC increases the velocity of the bubble reduces and the trajectory of bubble becomes rectilinear. Thus, at high CMC concentrations, the spiral vortex rope structure is replaced by a vortex ring, which follows the bubble due to upward pressure induced by the velocity gradient. In the study of Li et al. (2012), a CMC concentration of 2 g.L⁻¹ showed a spiral vortex rope but here the bubble in the same concentration is showing a vortex ring (Fig. 13d), one possible reason of this difference can be the geometry of column.

For bubbles rising in XG solution, a similar flow structure as the bubbles in CMC solution is observed (zones I and II). In general, for both shear-thinning solutions, the simulations have shown that the increase in bubble diameter enlarges the zones I and II. The shape of the vortex

ring (zone II) changes from circular to elliptical with vertical elongation as the bubble diameter increases. The flow in the wake remains always upward throughout the length of the numerical domain.

### 3.5. Strain rate and viscosity around the bubble

Fig. 14 and Fig. 15 present the numerical strain rate and viscosity fields respectively for the smallest and largest bubble diameters in a 2 × 2 cm² image for the three fluids, at a location where terminal velocity is attained.

For water, the viscosity is obviously constant. The strain rate field depends on the bubble size: for $d_b = 1.74$ mm (Fig. 14d), it is only slightly meandering due to minor bubble fluctuations (Fig. 11). The contour scale for the strain rate is between 0 and 300 s⁻¹. It can be seen that the high strain rate region observed around the bubble have values of more than 300 s⁻¹ and lower values of this parameter up to 180 s⁻¹ are noticed in the wake of the bubble. On the other hand, the 3.28 mm bubble case shown in Fig. 14d has high strain rate values more than 300 s⁻¹ observed around the bubble but also in the wake due to non-uniform wake structure (presence of vortices), non-rectilinear trajectories and more frequent wall interactions (as shown in Fig. 11).

The bubble rise in a non-Newtonian fluid contributes to the change of the liquid viscosity around the bubble. A direct link between the magnitude of the strain rate and the apparent viscosity can be seen in the numerical contours obtained from simulations, where the viscosity of the shear-thinning liquid is small in the high strain rate region and vice versa, following the power-law. As discussed earlier, the bubble trajectory in shear-thinning fluids is rectilinear with a stable shape and very few oscillations resulting in no meandering motion in the wake. Fig. 14 (b,e) and Fig. 15(b,e) show the strain rate and viscosity fields respectively for CMC (2 g.L⁻¹). It can be seen that in comparison to water, the bubble rising in CMC solution induces a lower length of the wake for the largest bubble diameter, hence the highest values of strain rate more than 300 s⁻¹ are observed just at the interface of the bubble and the
Streamtraces with velocity vectors showing fluid flow

Water ($d_b$ = 1.74 mm)

<figure><img src="image_17.png" /></figure>

(a)

Water ($d_b$ = 3.28 mm)

<figure><img src="image_18.png" /></figure>

(b)

CMC ($d_b$ = 1.58 mm)

<figure><img src="image_19.png" /></figure>

(c)

CMC ($d_b$ = 3.28 mm)

<figure><img src="image_20.png" /></figure>

(d)

XG ($d_b$ = 1.34 mm)

<figure><img src="image_21.png" /></figure>

(e)

XG ($d_b$ = 3.36 mm)

<figure><img src="image_22.png" /></figure>

(f)

<caption>Fig. 13. Streamlines superimposed with velocity vectors of smallest and largest bubble diameters in the middle vertical plane.</caption>
<figure><img src="image_23.png" /><figcaption>Fig. 14. Strain rate fields for the smallest and largest bubble rising in three fluids: front and top views.</figcaption></figure>

immediate wake (forming a lower viscosity tail) after the bubble whereas the lowest values up to 120 s⁻¹ are observed in the 1 cm periphery of the bubble. For the lowest bubble diameter, d_b = 1.58 mm, high strain rates are around the bubble.

For XG, which has more shear-thinning behaviour with n = 0.39, the results of d_b = 1.34 mm, corresponding to the smallest bubble in this study, shows that the strain rate is always high around the bubble having a value of 300 s⁻¹ with a small wake size. The viscosity behaviour for this bubble diameter shows a lower viscosity region with 0.005 Pa.s value near the bubble interface and 0.011 Pa.s in the immediate wake

(Fig. 15c). For the largest bubble of this study, which is an ellipsoidal bubble (d_b = 3.36 mm), Fig. 14f shows a strain rate of 300 s⁻¹ around the bubble and a high strain rate wake behind the bubble. It can also be seen for this bubble size that the minimum viscosity (0.005 Pa.s) is achieved near the bubble interface and the maximum viscosity (0.019 Pa.s) region around the bubble and in the wake. Similar viscosity distribution was observed by Premlata et al. (2017b) for shear-thinning fluids with power law index n = 0.8 and 0.6 using VOF method with the Carreau-Yasuda model to treat the viscosity of shear-thinning fluids in a confined cylindrical domain.
<figure><img src="image_24.png" /><figcaption>Fig. 15. Viscosity fields for the smallest and largest bubble rising in the three fluids: front and top views.</figcaption></figure>

It is worth mentioning that the top view of the 4 mm gap of the column in Fig. 14 and Fig. 15, shows that the highest strain rate and lowest viscosity region are located near the bubble interface and in the liquid film between the wall and bubble for the largest diameter. This will be of importance for the wall shear stress induced by the bubble for this configuration.

### 3.6. Wall shear stress (WSS) induced by bubble

The rise of an isolated bubble in a confined column generates shear forces at the wall. Particularly for PBR application, these shear stress

forces can contribute to mitigating the problem of biofilm establishment at the wall. For thin-gap configuration/hele-shaw cells, it was also observed by Roudet et al. (2017), that in these confined domains the mass transfer is more important due to the contribution of thin liquid films between the bubbles and the walls. In the literature, few numerical parametric studies on wall shear stress generated in a confined domain are available. Prieske et al. (2010) performed VOF simulation of single bubbles rising in water inside flat channels with gap thickness between 3 and 7 mm. Their study aimed at quantifying the maximum wall shear stress exerted by single bubble rising. They observed that the highest WSS values have been obtained for the smallest gap channel. For
instance, WSS of 3 mm and 5 mm bubbles in 3 mm gap configuration were 1.5 Pa and 3.3 Pa respectively in water as an operating fluid. In addition to that, the WSS increases with the increase of the bubble diameter up to a critical diameter and above which, the bubble did not generate higher shear stresses. Another similar conclusion has been also reported by Ndinisa et al. (2006), who observed that the cleaning efficiency increased with the increase of bubble diameter. However, when the bubble size becomes larger than the gap/width size of the column, the cleaning efficiency has a minor effect with the further increase in bubble size. The maximum shear stress reported by Ndinisa et al. (2006) is 0.7 Pa. Böhm et al. (2013a), motivated by the goal of removing the deposition layers from the membrane surfaces performed shear stress statistical analysis for single bubbles rising in a flat sheet confined configuration having gap thickness of 3–7 mm and using the electro-diffusion method. The highest maximum shear stress values for the 9 mm bubble diameter rising in a 5 mm gap were reported as 1.3 Pa (for

steady analysis) and 4.7 Pa (for transient corrected analysis). Recently, Thobie et al. (2022b) estimated wall shear stress via the electro-diffusion method in a channel with cross section of 4 × 4 mm² and reported a maximum shear stress of around 0.8–1 Pa.

Wall shear stress induced by bubble passage, for the lowest and the largest bubble diameters in three fluids is shown in Fig. 16 and the maximum wall shear stress for each case is presented in Table 3. Alike strain rate, the wall shear stress also shows that for water, due to meandering motion, the WSS is strongly present in wake of the bubble showing values maximum up to 23.5 Pa observed around the largest bubble (Table 3). It can also be seen in Fig. 13b that for water the velocity field for the largest bubble diameter ($d_b$ = 3.28 mm) shows a meandering motion due to bubble oscillations, and a higher number of interactions with the walls (Fig. 11), the impact of these interactions is evident on WSS (Fig. 16d). Besides that, the velocity profiles for sheath-thinning fluids have a symmetrical distribution around the bubble and

<figure><img src="image_25.png" /><figcaption>Fig. 16. Wall shear stress field for the smallest and largest bubble rising in three fluids.</figcaption></figure>
<table><caption>**Table 3**
Maximum wall shear stress induced by the bubble in three fluid environments on front wall.</caption><thead><tr><th>Fluid</th><th>Bubble Diameter (mm)</th><th>Maximum shear stress range (Pa)</th></tr></thead><tbody><tr><td rowspan="2">Water (n = 1)</td><td>1.74</td><td>0.80</td></tr><tr><td>3.28</td><td>23.5</td></tr><tr><td rowspan="2">CMC (2 g.L<sup>-1</sup>)<br>(n = 0.87)</td><td>1.58</td><td>0.70</td></tr><tr><td>3.28</td><td>6.50</td></tr><tr><td rowspan="2">XG (1 g.L<sup>-1</sup>)<br>(n = 0.39)</td><td>1.34</td><td>0.60</td></tr><tr><td>3.36</td><td>4.20</td></tr></tbody></table>

the wake is forming a tail. This information can be linked to the rectilinear trajectory of bubbles in shear-thinning fluids. The values of the WSS are maximum (6.5 Pa for CMC and 4.2 Pa for XG) around the bubble and in the tail behind it (Table 3). This observation is for both the front and the back walls. Whereas for the lowest bubble diameters, as the confinement ratio is very low, the bubble usually has fewer interactions or interaction with only one wall, in that manner, the WSS has a maximum value of 0.7 Pa.

This must be noted that these analyses of strain rate, viscosity and wall shear stress are based on a relative position of the bubble in the gap, so the results can slightly differ if another instant is considered but the general physical phenomena are conserved.

## 4. Conclusion

In this study, the Volume of Fluid (VOF) method was used to simulate single bubbles rising in Newtonian and non-Newtonian fluids inside a thin-gap bubble column with a 4 mm gap. Different bubble sizes rising in water, CMC and XG aqueous solutions have been simulated. The numerical results were compared to the experimental bubble terminal velocity, shape and trajectory measured experimentally through the shadowgraphy technique by Almani et al. (2021) and show similar trends. The terminal velocity in water decreases with the increase of bubble diameter and then increases as in experimental results, whereas for non-Newtonian fluids a plateau is achieved after $d_b = 2.2$ mm. However, slightly lower terminal velocities have been observed in comparison to the experimental results. The numerical values of the aspect ratio of the bubbles are lower than that of the experimental ones for all fluids considered in this study. The trajectory of the bubble shows the same trend as seen experimentally: non-rectilinear for water and rectilinear for non-Newtonian fluids.

Besides that, the flow field around the bubbles, strain rate and viscosity profiles around the bubble as well as wall shear stress induced by the bubble were also studied numerically. Bubbles in water have more oscillations and interaction with the wall causing meandering motion that influences the strain rate, wall shear stress and the velocity profiles near the bubble and also in the wake. Besides that, for water, the flow field also shows the Von Karman vortex street structure behind the bubble. Whereas bubbles in the non-Newtonian fluids have relatively lower terminal velocity, stable bubble shape and rectilinear trajectory and show stable profiles of strain rate, wall shear stress and velocity, where the highest values of these parameters are achieved near the bubble or in the immediate wake. It is also observed that the length of the wake decreases as the viscosity of the liquid phase increases.

This numerical study has enabled to understand the bubble/liquid and bubble/wall interactions more deeply, which was not possible experimentally, such as the visualization of the strain rate, viscosity and lateral/ top views of the column during the bubble rise. These informations are of operational interest in processes such as intensified photobioreactors, or flat membrane modules and can be used for further numerical investigation of such processes.

## CRediT authorship contribution statement

**Caroline Gentric:** Writing – review & editing, Validation, Supervision, Methodology, Funding acquisition, Conceptualization. **Sikandar Almani:** Writing – original draft, Methodology, Investigation. **Walid Blel:** Writing – review & editing, Supervision, Methodology. **Emilie Gadoin:** Writing – review & editing, Supervision, Methodology.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

The authors would like to gratefully acknowledge the financial support from the Higher Education Commission (HEC), Pakistan Grant Numbers: Ref-50027885/PM/UESTPs/HRD/HEC and also to Campus France for the kind cooperation in Sikandar ALMANI's Ph.D.

## References

Ahmad, I., Husain, S., Hassan, L.H.S., Khan, N.A., Zahmatkesh, S., Ni, B.J., 2024. Computational fluid dynamics (CFD) simulation modeling for the cultivation of microalgal monoculture in axenic enclosed bubble column photobioreactor (BCPBR). Algal Res. 83, 103725. https://doi.org/10.1016/j.algal.2024.103725.

Almani, S., Blel, W., Gadoin, E., Gentric, C., 2021. Investigation of single bubbles rising in Newtonian and non-Newtonian fluids inside a thin-gap bubble column intended for microalgae cultivation. Chem. Eng. Res. Des. 167, 218–230. https://doi.org/10.1016/j.chedr.2021.01.010.

Alméras, E., Risso, F., Roig, V., Plais, C., Augier, F., 2018. Mixing mechanism in a two-dimensional bubble column. Phys. Rev. Fluids 7, 1–19. https://doi.org/10.1103/PhysRevFluids.3.074307.

Asim, M., Almani, S., Unar, I.N., Abro, M., Siddiqui, S., Qureshi, K., Solangi, Z.A., Soomro, S.A., 2023. Isolated bubble ascent in a non-Newtonian media inside an infinite bubble column: a CFD study. Braz. J. Chem. Eng. https://doi.org/10.1007/s43153-023-00393-0.

Battistella, A., van Schijndel, S.J.G., Baltussen, M.W., Roghair, I., van Sint Annaland, M., 2020. On the terminal velocity of single bubbles rising in non-Newtonian power-law liquids. J. Nonnewton. Fluid Mech. 278, 104249. https://doi.org/10.1016/j.jnmfm.2020.104249.

Bertola, F., Baldi, G., Marchisio, D., Vanni, M., 2004. Momentum transfer in a swarm of bubbles: Estimates from fluid-dynamic simulations. Chem. Eng. Sci. 59, 5209–5215. https://doi.org/10.1016/j.ces.2004.07.134.

Böhm, L., Brehmer, M., Kraume, M., 2016a. Comparison of the single bubble ascent in a Newtonian and a Non-Newtonian liquid: A phenomenological PIV study. Chem.-Ing.-Tech. 88, 93–106. https://doi.org/10.1002/cite.201500105.

Böhm, L., Drews, A., Kraume, M., 2013a. Bubble induced shear stress in flat sheet membrane systems—Serial examination of single bubble experiments with the electrodiffusion method. J. Memb. Sci. 437, 131–140. https://doi.org/10.1016/j.memsci.2013.02.036.

Böhm, L., Kolano, M., Kraume, M., 2016b. Simulation of the single bubble ascent with openfoam. Tech. Trans. 25–32.

Böhm, L., Kurita, T., Kimura, K., Kraume, M., 2014. Rising behaviour of single bubbles in narrow rectangular channels in Newtonian and non-Newtonian liquids. Int. J. Multiph. Flow. 65, 11–23. https://doi.org/10.1016/j.ijmultiphaseflow.2014.05.001.

Böhm, L., Prieske, H., Kraume, M., 2013b. Fluid dynamic optimization of flat sheet membrane modules - movement of bubbles in vertical channels. Chem. Eng. Trans. 32, 1501–1506. https://doi.org/10.33032/CET1332251.

Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surface tension. J. Comput. Phys. 100, 335–354. https://doi.org/10.1016/0021-9991(92)90240-Y.

Clift, R., Grace, J.R., Weber, M.E., 1978. Bubbles, drops and particles. Academic Press. Academic Press, NY.

Cornet, J.F., Dussap, C.G., 2009. A simple and reliable formula for assessment of maximum volumetric productivities in photobioreactors. Biotechnol. Prog. 25, 424–435. https://doi.org/10.1002/btpr.138.

de Jesus, S.S., Moreira Neto, J., Maciel Filho, R., 2017. Hydrodynamics and mass transfer in bubble column, conventional airlift, stirred airlift and stirred tank bioreactors, using viscous fluid: A comparative study. Biochem. Eng. J. 118, 70–81. https://doi.org/10.1016/j.bej.2016.11.019.

Drews, A., Prieske, H., Meyer, E.L., Senger, G., Kraume, M., 2010. Advantageous and detrimental effects of air sparging in membrane filtration: Bubble movement, exerted shear and particle classification. Desalination 250, 1083–1086. https://doi.org/10.1016/j.desal.2009.09.113.

Ferrario, A., Varallo, N., Besagni, G., Mereu, R., 2025. Influence of the gas phase on a large-scale bubble column fluid dynamics: Gas holdup, flow regime transitions, and
bubble size distributions. Chem. Eng. Sci. 302, 120792. https://doi.org/10.1016/j.ces.2024.120792.

Fourati, M., Roig, V., Raynal, L., 2013. Liquid dispersion in packed columns: Experiments and numerical modeling. Chem. Eng. Sci. 100, 266-278. https://doi.org/10.1016/j.ces.2013.02.041.

Ghanavati, A., Khodadadi, S., Taleghani, M.H., Gorji-Bandpy, M., Domiri Ganji, D., 2023. Numerical simulation of the motion and interaction of bubble pair rising in a quiescent liquid. Appl. Ocean Res. 141, 103769. https://doi.org/10.1016/j.apor.2023.103769.

Gumulya, M., Utikar, R.P., Pareek, V.K., Evans, G.M., Joshi, J.B., 2021. Dynamics of bubbles rising in pseudo-2D bubble column: Effect of confinement and inertia. Chem. Eng. J. 405, 126615. https://doi.org/10.1016/j.cej.2020.126615.

Hassan, N.M.S., Khan, M.M.K., Rasul, M.G., Rackemann, D.W., 2007. An experimental study of bubble rise characteristics in non - Newtonian (Power-Law) fluids. th Australas. Fluid Mech. Conf. 16, 1315-1320.

Hirt, C.W., Nichols, B.D., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 39, 201-225. https://doi.org/10.1016/0021-9991(81) 90145-5.

Islam, M.T., Ganesan, P., Sahu, J.N., Uddin, M.N., Mannan, A., 2015. A Single Air Bubble Rise in Water: a CFD Study. Mech. Eng. Res. J. 9, 1-6.

Islam, M.T., Ganesan, P.B., Cheng, J., SalahUddin, M., 2020. Single bubble rising behaviors in Newtonian and non-Newtonian fluids with validation of empirical correlations: A computational fluid dynamics study. Eng. Rep. 2 (1), e12100. https://doi.org/10.1002/eng2.12100.

Keshavarzi, G., Pawell, R.S., Barber, T.J., Yeoh, G.H., 2014. Transient analysis of a single rising bubble used for numerical validation for multiphase flow. Chem. Eng. Sci. 112, 25-34. https://doi.org/10.1016/j.ces.2014.02.027.

Khamadieva, R., Böhm, U., 2006. Mass transfer to the wall of a packed and unpacked bubble column operating with Newtonian and non-Newtonian liquids. Chem. Eng. J. 116, 105-113. https://doi.org/10.1016/j.cej.2005.11.003.

Kong, G., Mirsandi, H., Buist, K.A., Peters, E.A.J.F., Baltussen, M.W., Kuipers, J.A.M., 2019. Oscillation dynamics of a bubble rising in viscous liquid. Exp. Fluids 60, 1-13. https://doi.org/10.1007/s00348-019-2779-1.

Krishna, R., Urseanu, M., van Baten, J.M., Ellenberger, J., 1999. Wall effects on the rise of single gas bubbles in liquids. Int. Commun. Heat. Mass Transf. 26, 781-790. https://doi.org/10.1017/CBO9781107415324.004.

Krishna, R., Van Baten, J.M., 1999. Rise characteristics of gas bubbles in a 2D rectangular column: VOF simulations vs experiments. Int. Commun. Heat. Mass Transf. 26, 965-974. https://doi.org/10.1016/S0735-1933(99)00086-X.

Kumar, P., Vanka, S.P., 2015. Effects of confinement on bubble dynamics in a square duct. Int. J. Multiph. Flow. 77, 32-47. https://doi.org/10.1016/j.ijmultiphaseflow.2015.06.014.

Li, S., Ma, Y., Fu, T., Zhu, C., Li, H., 2012. The viscosity distribution around a rising bubble in shear-thinning non-newtonian fluids. Braz. J. Chem. Eng. 29, 265-274. https://doi.org/10.1590/S0104-66322012000200007.

Li, X., Wang, W., Zhang, P., Li, J., Chen, G., 2019. Interactions between gas-liquid mass transfer and bubble behaviours. R. Soc. Open Sci. 6. https://doi.org/10.1098/rros.190136.

Ma, D., Liu, M., Zu, Y., Tang, C., 2012. Two-dimensional volume of fluid simulation studies on single bubble formation and dynamics in bubble columns. Chem. Eng. Sci. 72, 61-77. https://doi.org/10.1016/j.ces.2012.01.013.

Mahmoudi, S., Hemmatian, F., Dahkaee, K.P., Hlawitschka, M.W., Kantzas, A., 2022. Detailed study of single bubble behavior and drag correlations in Newtonian and non-Newtonian liquids for the design of bubble columns. Chem. Eng. Res. Des. 179, 119-129. https://doi.org/10.1016/j.cherd.2022.01.013.

Mei, Z., Cheng, X., 2024. Analysis of single bubble dynamics based on VOF method. Ann. Nucl. Energy 208, 110757. https://doi.org/10.1016/j.anucene.2024.110757.

Ndinisa, N.V., Fane, A.G., Wiley, D.E., 2006. Fouling control in a submerged flat sheet membrane system: Part I - Bubbling and hydrodynamic effects. Sep. Sci. Technol. 41, 1383-1409. https://doi.org/10.1080/01496390600633873.

Osher, S., Sethian, J.A., 1988. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. J. Comput. Phys. 79, 12-49. https://doi.org/10.1016/0021-9991(88)90002-2.

Premlata, A.R., Tripathi, M.K., Karri, B., Sahu, K.C., 2017a. Numerical and experimental investigations of an air bubble rising in a Carreau-Yasuda shear-thinning liquid. Phys. Fluids 29. https://doi.org/10.1063/1.4979136.

Premlata, A.R., Tripathi, M.K., Karri, B., Sahu, K.C., 2017b. Dynamics of an air bubble rising in a non-Newtonian liquid in the axisymmetric regime. J. Nonnewton. Fluid Mech. 239, 53-61. https://doi.org/10.1016/j.jnnfm.2016.12.003.

Prieske, H., Böhm, L., Drews, A., Kraume, M., 2010. Optimised hydrodynamics for membrane bioreactors with immersed flat sheet membrane modules. Desalin. Water Treat. 18, 270-276. https://doi.org/10.5004/dwt.2010.1784.

Radaei, E., Liu, X., Tng, K.H., Wang, Y., Trujillo, F.J., Leslie, G., 2018. Insights on pulsed bubble control of membrane fouling: Effect of bubble size and frequency. J. Memb. Sci. 554, 59-70. https://doi.org/10.1016/j.memsci.2018.02.058.

Rashid, F.L., Hussein, E., Hussein, A.K., Younis, O., 2023. Parametric Study of Single Air Bubble Rising Through Different Salinity Water Column Using Volume of Fluid (Vof) Technique. J. Eng. Sci. Technol. 18, 671-684.

Roudet, M., Billet, A.M., Cazin, S., Risso, F., Roig, V., 2017. Experimental investigation of interfacial mass transfer mechanisms for a confined high-reynolds-number bubble rising in a thin gap. AIChE J. 63, 2394-2408. https://doi.org/10.1002/aic.15562.

Souliès, A., Pruvost, J., Legrand, J., Castelain, C., Burghelea, T.I., 2013. Rheological properties of suspensions of the green microalga Chlorella vulgaris at various volume fractions. Rheol. Acta 52, 589-605. https://doi.org/10.1007/s00397-013-0700-z.

Sussman, M., Puckett, E.G., 2000. A Coupled Level Set and Volume-of-Fluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows. J. Comput. Phys. 162, 301-337. https://doi.org/10.1006/jcp.2000.6537.

Thobie, C., 2018. Caractérisation de l'hydrodynamique et des transferts gaz-liquide dans un photobioreacteur intensifié. Etude de l'effet du bullage sur le développement de biofilm micro-algal. Thèse de Doctorat,Université de Nantes.

Thobie, C., Blel, W., Hauser, J.L., Pruvost, J., Gentric, C., 2022a. Different types of bubbly flows in a confined channel with the aim of limiting microalgae biofilm development - Part I: Hydrodynamic study. Chem. Eng. Process. - Process. Intensif. 173. https://doi.org/10.1016/j.cep.2022.108844.

Thobie, C., Blel, W., Hauser, J.L., Pruvost, J., Gentric, C., 2022b. Different types of bubbly flows in a confined channel with the aim of limiting microalgae biofilm development - Part I: Hydrodynamic study. Chem. Eng. Process. - Process. Intensif. 173. https://doi.org/10.1016/j.cep.2022.108844.

Thobie, C., Gadoin, E., Blel, W., Pruvost, J., Gentric, C., 2017. Chemical Engineering & Processing: Process Intensification Global characterization of hydrodynamics and gas-liquid mass transfer in a thin-gap bubble column intended for microalgae cultivation. Chem. Eng. Process. Process. Intensif. 122, 76-89. https://doi.org/10.1016/j.cep.2017.10.009.

Vaishnavi, G.N.V.S., Ramarajan, J., Jayavel, S., 2023. Numerical studies of bubble formation dynamics in gas-liquid interaction using Volume of Fluid (VOF) method. Therm. Sci. Eng. Prog. 39, 101718. https://doi.org/10.1016/j.tsep.2023.101718.

van Sint Annaland, M., Deen, N.G., Kuipers, J.A.M., 2005. Numerical simulation of gas bubbles behaviour using a three-dimensional volume of fluid method. Chem. Eng. Sci. 60, 2999-3011. https://doi.org/10.1016/j.ces.2005.01.031.

Wang, X., Blanpain, B., Degrève, J., 2015. Numerical Simulation of Two-Dimensional Bubble Dynamics and Evaporation.

Wang, X., Klaasen, B., Degrève, J., Blanpain, B., Verhaeghe, F., 2014. Experimental and numerical study of buoyancy-driven single bubble dynamics in a vertical Hele-Shaw cell. Phys. Fluids 26, 0-18. https://doi.org/10.1063/1.4903488.

Xu, X., Li, X., Li, Z., Zhang, D., 2024. Rising behavior and mass transfer of CO2 bubbles in shear-thinning non-Newtonian fluids. Chem. Eng. Res. Des. 207, 232-242. https://doi.org/10.1016/j.cherd.2024.06.003.

Yu, H., Yang, S., Wang, H., 2025. Investigation of multi-inlet arrangement effect on mesoscale bubble dynamics in the gas-liquid-solid reactor via VOF-DEM. Chem. Eng. Sci. 301, 120710. https://doi.org/10.1016/j.ces.2024.120710.

Yujie, Z., Mingyan, L., Yonggui, X., Can, T., 2012. Three-dimensional volume of fluid simulations on bubble formation and dynamics in bubble columns. Chem. Eng. Sci. 73, 55-78. https://doi.org/10.1016/j.ces.2012.01.012.

Zahedi, P., Saleh, R., Moreno-Atanasio, R., Yousefi, K., 2014. Influence of fluid properties on bubble formation, detachment, rising and collapse; Investigation using volume of fluid method. Korean J. Chem. Eng. 31, 1349-1361. https://doi.org/10.1007/s11814-014-0063-x.

Zhang, L., Yang, C., Mao, Z.S., 2010. Numerical simulation of a bubble rising in shear- thinning fluids. J. Nonnewton. Fluid Mech. 165, 555-567. https://doi.org/10.1016/j.jnnfm.2010.02.012.