# Simulation of bubbly flow in loop reactors using the Eulerian-Eulerian two-fluid method: Analysis and assessment of medium-to-high liquid flow rate effect

M. Badalan, M. Shiea, Y. Haroun, G.E. Vázquez-Angulo, J.E.V. Guzmán

► To cite this version:

M. Badalan, M. Shiea, Y. Haroun, G.E. Vázquez-Angulo, J.E.V. Guzmán. Simulation of bubbly flow in loop reactors using the Eulerian-Eulerian two-fluid method: Analysis and assessment of medium-to-high liquid flow rate effect. Chemical Engineering Research and Design, 2025, 221, pp.497-511. 10.1016/j.cherd.2025.08.001. hal-05265664

HAL Id: hal-05265664

https://ifp.hal.science/hal-05265664v1

Submitted on 17 Sep 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
# Simulation of bubbly flow in loop reactors using the Eulerian-Eulerian two-fluid method: Analysis and assessment of medium-to-high liquid flow rate effect

M. Badalan<sup>a,b</sup>, M. Shiea<sup>b</sup>, Y. Haroun<sup>b</sup>, G.E. Vázquez-Angulo<sup>a</sup>, J.E.V. Guzmán<sup>a,*</sup>

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

## ARTICLE INFO

**Keywords:**

Bubbly flow

High superficial liquid velocity

CFD

Interfacial forces

Reactor

Two-fluid model

## ABSTRACT

Loop reactors in the chemical industry are characterised by gas-liquid flow in a vertical arrangement, including a riser and a downcomer, at relatively high liquid velocities. A standard method in the industry for simulating such flows in large geometries is the Eulerian two-fluid model, which relies on interfacial forces to describe the interactions (momentum transfer) between the two phases. However, most of these models have been validated with data from bubble columns and low-velocity pipe flows. Therefore, we assess the accuracy of these models against experimental data from medium-to-high velocity (upward and downward) bubbly flows in pipes. The OpenFOAM simulation framework is used for this purpose. For upward flows, current models fail to predict the shift of the void fraction peak towards the centre of the pipe at high liquid velocities. This work discusses several development directions, including adjustments to the lift and wall lubrication coefficients, which significantly improve prediction accuracy. For downward flows, the void fraction distribution is highly dependent on inlet conditions, which aligns with previous results from the literature. This study highlights the mispredictions of conventional models and provides valuable insights into optimising interfacial force models for more accurate simulations in chemical loop reactors.

## 1. Introduction

Loop reactors are highly versatile chemical engineering systems that integrate the advantages of continuous stirred-tank reactors (CSTR) and plug-flow reactors (PFR), offering enhanced control over mixing, heat transfer, and residence time distribution. They are characterised by a defined recirculation flow through a pipe loop, which can be driven mechanically or hydrodynamically, and can operate in batch or continuous configurations. Depending on the design, loop reactors can be classified as airlift loop reactors (ALR), propeller loop reactors (PLR), or jet loop reactors (JLR), with further subdivisions into internal-loop and external-loop systems. These reactors are employed across diverse fields such as biotechnology, polymer synthesis, fine chemicals production, and gas-liquid reactions. For instance, they are used in the manufacture of single-cell proteins, azo dyes, polyethylene and polypropylene, and in rapid gas-liquid reactions like hydrogenation and alkane chlorination. Their ability to achieve efficient mixing and high interphase mass-transfer rates, coupled with good heat transfer performance and low fouling tendencies, makes them ideal for handling multiphase systems, including gas-liquid, liquid-liquid, and slurry reactions (Blenke, 1979; Moran et al., 2017).

The propeller loop reactors (PLR) used for polymerisation are generally external loop reactors, characterised by high liquid velocities (ranging from 5 to 20 m/s) and a bubbly flow regime (Fig. 1). The ability to accurately predict flow regimes, bubble size distribution, pressure drop, and mass transfer rates is crucial for the development and sizing of such reactors. In this regard, computational fluid dynamics (CFD) has progressively emerged as a promising tool for predicting the behaviour of bubbly flows in different geometries. In the available literature, the vast majority of the numerical studies on loop reactors focus on ALRs (Becker et al., 1994; Sokolichin and Eigenberger, 1994; Vial et al., 2002; Gao et al., 2019; Roy and Joshi, 2008; Puiman et al., 2022; Talvy et al., 2007; Huang et al., 2010; Moudoud et al., 2018; Burlutskii and Di Felice, 2019; Teli and Mathpati, 2021; Kojic et al., 2022). To the knowledge of the authors, very few CFD studies on PLRs have been published so far. An example of a numerical analysis of a PLR is the work by Orejas (1999, 2001), which placed emphasis on the energy balance and the chemical reactions. However, the characteristics of the bubbly flow were not explicitly resolved.

<footnote><sup>a</sup> Universidad Nacional Autónoma de México, Instituto de Ingeniería, Cicto. Esc. SN, CdMx, 04510, Coyoacán, Mexico</footnote>

<footnote><sup>b</sup> IFP Énergies Nouvelles, Rond-point de l'échangeur de Solaize, Solaize, 69360, France</footnote>

<footnote>* Corresponding author.
E-mail address: jguzmanv@iingen.unam.mx (J.E.V. Guzmán).</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Schematic simplified example of a propeller loop reactor (PLR) for ethylene polymerisation.</figcaption></figure>

Source: Adapted from Norwood (1966).

One potential reason for this gap is the lack of robust CFD models that are developed and validated for this type of flow. Among the available CFD models commonly used to simulate two-phase bubbly flow, the Eulerian two-fluid model has been widely employed for the simulation of industrial scale systems. This computational method relies on a series of interfacial forces for predicting the momentum exchange between the gas and the liquid (e.g. the drag, lift, turbulent dispersion, virtual mass and wall lubrication forces) (Ishii and Hibiki, 2011). The majority of these force models have been developed and validated in bubble columns and gas-liquid pipe flows with liquid velocities below 1 m/s (Khan et al., 2020; Kim et al., 2021). Even large projects that aimed at developing a baseline CFD model for bubbly flows (which would work robustly in different flow setups), as for instance the one developed at HZDR in Germany, have used mainly

low-liquid velocity data for validation (Liao et al., 2019). Furthermore, experimental data for vertical pipe flows with a liquid velocity larger than 4 m/s is scarce, limiting the development and validation of new models.

Thus, despite the advancements in CFD modelling for bubbly flows, there remains a significant gap in the validation of two-fluid models for the specific flow conditions encountered in loop reactors, particularly in PLRs. Current models have been validated for flows with liquid velocities up to 1 m/s, far below the 5–12 m/s observed in industrial loop reactors. This discrepancy highlights the urgent need for models that can accurately predict the behaviour of bubbly flows at industrially relevant velocities. Moreover, in the case of some of the interfacial forces that are highly relevant for upward pipe flows, such as the lift force, the majority of the existing models, primarily developed for laminar or low-turbulence conditions, fail to account for the complex interactions in highly turbulent flows, leading to inaccuracies in predicting void fraction distributions and other critical parameters (Salibindla et al., 2020).

Our work highlights these shortcomings by assessing state-of-the-art two-fluid models for the simulation of flow conditions closer to those encountered in industrial loop reactors, i.e. medium-to-high liquid velocity (4 m/s) vertical upward and downward bubbly flows. Results are compared with experimental data from the literature and discussed. Based on insights from previous experimental and theoretical work by Salibindla et al. (2020), we place a particular emphasis in our discussion on the lift force, and on how it is affected by turbulence. We highlight that, as this effect is not taken into account by state-of-the art models, they fail in predicting the bubble distributions at flows with 4 m/s. Moreover, we propose an optimised model that performs significantly better in medium-to-high liquid velocity bubbly flows, which can serve as a basis for the simulation of loop reactors.

This research not only highlights the need for developing a more robust framework for simulating industrially relevant flow regimes and proposes a theoretical development direction based on the effect of turbulence on the lift force, but it also proposes a numerically optimised model that is shown to perform well for these flow conditions, and that is validated on different independent sets of experimental data.

## 2. Summary of the experimental data used

Several sets of experimental data were used for the evaluation of the numerical models. These sets were selected based on the liquid velocities and the availability of local measurements of flow properties, e.g. void fraction, bubble sizes. Moreover, in the selected experiments, the flow conditions were curated to maintain a consistent flow regime (bubbly flow), thus avoiding changes in flow morphology within the specified superficial velocity range. This ensures that the data represents a consistent flow structure, which is both crucial for the applicability of our model and relevant for the flow in PLRs.

The different experimental setups are summarised in Table 1, and will be described briefly in the following sections. The selection contains two sets of experimental data (one for the riser and one for the downcomer case) at medium-to-high liquid velocities (Ryan and Kim, 2023; Qiao et al., 2017) (up to 4 m/s), which are among the highest velocities available in the literature where detailed local measurements were performed, to the authors' knowledge. While there are other few datasets available in the literature for upward bubbly flow in the same range of liquid velocities (Serizawa et al., 1992; Hibiki et al., 2001; Ishii et al., 2004), the most recent available dataset, i.e. the one by Ryan and Kim (2023), was selected due to constraints on time and resources.

Two other sets of experiments (Worosz et al., 2015; Kashinsky and Randin, 1999) were selected at low-to-medium liquid velocities, to serve as a baseline comparison for the riser and downcomer cases (since most of the numerical models for bubbly flows have been validated with low liquid velocity experiments). The relevance of these experiments for the different parts of the loop reactor is highlighted in Fig. 2. All the experiments in the selected references were conducted using air and water as working fluids.
<table><caption>**Table 1**
Summary of the experimental cases used for validation.</caption><thead><tr><th>Case name</th><th>Liquid superficial velocity (m/s)</th><th>Gas superficial velocity (m/s)</th><th>Bubble diameter range (mm)</th><th>Pipe diameter (mm)</th><th>Pipe length (m)</th><th>Configuration</th><th>Reported measurement uncertainty<sup>a</sup></th><th>Ref.</th></tr></thead><tbody><tr><td>R1</td><td>4.0</td><td>0.111</td><td>1.6–1.9</td><td>25.4</td><td>2.794</td><td>Riser</td><td></td><td></td></tr><tr><td>R2</td><td>4.0</td><td>0.186</td><td>1.5–2.0</td><td>25.4</td><td>2.794</td><td>Riser</td><td>u<sub>G</sub>: ±10%<sup>b</sup></td><td>Ryan and Kim (2023)</td></tr><tr><td>R3</td><td>4.0</td><td>0.333</td><td>1.6–2.0</td><td>25.4</td><td>2.794</td><td>Riser</td><td></td><td></td></tr><tr><td>R4</td><td>1.56</td><td>0.111</td><td>2.2–3.0</td><td>25.4</td><td>2.794</td><td>Riser</td><td>a<sub>G</sub>: ±10%<sup>c</sup></td><td></td></tr><tr><td>R5</td><td>2.5</td><td>0.111</td><td>1.7–2.6</td><td>25.4</td><td>2.794</td><td>Riser</td><td></td><td></td></tr><tr><td>W1</td><td>2.0</td><td>0.18</td><td>3.1–3.4</td><td>50.8</td><td>3.2</td><td>Riser</td><td>a<sub>G</sub>: ±10%<sup>d</sup></td><td>Worosz et al. (2015)</td></tr><tr><td>Q1A</td><td>4.0</td><td>0.108</td><td>1.6–3.2</td><td>50.8</td><td>2.82</td><td>Downcomer inlet type A</td><td></td><td></td></tr><tr><td>Q1B</td><td>4.0</td><td>0.108</td><td>1.8–2.4</td><td>50.8</td><td>2.82</td><td>Downcomer inlet type B</td><td>u<sub>G</sub>: ±10%</td><td>Qiao et al. (2017)</td></tr><tr><td>Q1C</td><td>4.0</td><td>0.108</td><td>1.8–2.3</td><td>50.8</td><td>2.82</td><td>Downcomer inlet type C</td><td>a<sub>G</sub>: ±10%<sup>e</sup></td><td></td></tr><tr><td>K1</td><td>1.0</td><td>0.019</td><td>~1.5<sup>f</sup></td><td>42.3</td><td>4.8</td><td>Downcomer</td><td>a<sub>G</sub>: ±20%</td><td>Kashinsky and Randin (1999)</td></tr><tr><td>K2</td><td>1.0</td><td>0.0403</td><td>~1.5</td><td>42.3</td><td>4.8</td><td>Downcomer</td><td>u<sub>L</sub>: ±4%<sup>g</sup></td><td></td></tr><tr><td>K3</td><td>1.0</td><td>0.0916</td><td>~1.5</td><td>42.3</td><td>4.8</td><td>Downcomer</td><td></td><td></td></tr></tbody></table>

## 2.1. Riser cases - Ryan et al. (2023)

The first experimental riser is a tube with ID of 25.4 mm and with a gas and liquid injection at the bottom, as described in Ryan and Kim (2023) and Ryan et al. (2023). The tube is mounted on a mechanical system that can be inclined from 0° (horizontal) to 90° (vertical). Three measurement ports are located at 30 D, 66 D and 110 D respectively, and are used for the measurement of pressure and profiles of void fraction, gas velocity, bubble Sauter mean diameter and interfacial area. In that work, the authors investigated five different flow regimes (presented in Table 1) in five different pipe orientations. Here, we only simulate the five vertical upward flow (90° pipe orientation) cases, which are particularly relevant for the riser of PLRs.

## 2.2. Riser cases - Worosz et al. (2015)

To establish a comparative baseline, we chose an additional set of experiments (Worosz et al., 2015) conducted at low-to-medium liquid velocities. This selection complements the riser cases (Ryan and Kim, 2023) and lies closer to the velocity ranges typically used for validating numerical models of bubbly flows. In the work published by Worosz et al. (2015), a 50.8 mm vertical pipe is used to investigate the gas-liquid flow in 9 different experimental conditions. Profiles of void fraction and interfacial area are measured at 18 D, 34.5 D and 63 D. Here, we consider only the first run of that work, which is classified in the bubbly flow regime according to the authors.

## 2.3. Downcomer cases - Qiao et al. (2017)

Qiao et al. (2017) investigated the effects of the inlet configuration on downward bubbly and slug flows in vertical pipes of ID 50.8 mm. Four different flow regimes, including the bubbly flow, were investigated using three different inlet types (type A: two-phase flow injector without flow straightener; type B: two-phase flow injector with a flow straightener; type C: horizontal-to-vertical-downward flow through a 90° vertical elbow). Pressure drop, and the profiles of void fraction, interfacial area, gas velocity, bubble Sauter mean diameter and bubble frequency were measured at three locations (7.5 D, 31.5 D and 55.5 D for inlet types A and B; 3 D, 16.5 D and 67.5 D for inlet type C). Here we consider only run 1 from that work, which is clearly in the bubbly flow regime, with the three different inlets.

## 2.4. Downcomer cases - Kashinsky and Randin (1999)

Kashinsky and Randin (1999) investigated the downward gas-liquid flow in a pipe of 42.3 mm inner diameter in 22 different flow regimes. The void fraction profiles and the wall shear stress were measured at 3.4 m and 4.44 m from the two-phase mixer. In this work, we consider only three runs from this reference, corresponding to the largest liquid flow rate (1 m/s) and three different gas flow rates.

## 3. Numerical setup

For the CFD simulation of bubbly flow in this work, we use the Euler-Euler two-fluid model, which has been presented in detail in a number of reference works and will not be repeated here (Drew and Passman, 1999; Yeoh and Tu, 2009; Ishii and Hibiki, 2011). We will use the implementation of the two-fluid model in the open-source software OpenFOAM, specifically the modules MultiphaseEuler and addonMultiphaseEuler (Direct, 2021; Schlegel et al., 2024). The momentum equations of the two-fluid model contain several closure terms representing the interaction (forces) between the two phases, which will be described briefly in the following section.

Both liquid and gas were considered as incompressible fluids, with the liquid density being constant (equal to water density at 20 °C) and the gas density varying according to the ideal gas law (Rzehak et al., 2017). To account for turbulence, we use the k-omega SST turbulence model in the liquid, but we consider the gas as a laminar fluid (assuming negligible turbulence inside the bubbles) (Rzehak et al., 2015, 2017; Liao et al., 2019). Bubble-induced turbulence (BIT) can become important in certain flow conditions (low liquid velocity, large bubble sizes >6 mm, high gas volume fraction) (Liao et al., 2019). However, these conditions are not met in the experiments considered in this work. Therefore, the intrinsic turbulence of the liquid flow is expected to dominate over BIT effects. Moreover, a few preliminary tests including BIT modelling showed that it did not have a notable influence on the results in these cases (see supplementary material SM.1). As a consequence, BIT will not be taken into account in the simulations in this work.

The turbulent boundary layer is taken into account by using standard wall functions, allowing to save computational time.

<footnote>a Estimated by comparing different measurement techniques.</footnote>

<footnote>b Overall, 95% of the measurements agree within ±10%, with an average absolute percent difference (AAPD) of 4.7%.</footnote>

<footnote>c Almost all conditions are shown to agree within ±10%, with an AAPD of 4.6%.</footnote>

<footnote>d Most of the measurements agree within ±10%, with an average disagreement of about –8%.</footnote>

<footnote>e Higher uncertainty was found for high gas flow rate conditions, which were not used in the present work.</footnote>

<footnote>f The authors of Kashinsky and Randin (1999) provide only an average bubble size, but not the entire range of bubble sizes.</footnote>

<footnote>g At a probability level of 95%.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 2. Parts of the loop reactor for which the different experimental sets are relevant.</figcaption></figure>

### 3.1. Interfacial forces

The momentum exchange between the two fluids in the Eulerian two-fluid framework is modelled by interfacial forces. The most relevant ones in the simulation of gas-liquid bubbly flow are the drag force, lift force, virtual mass, turbulent dispersion and wall lubrication. The drag force represents the resistance to the bubbles due to their motion inside the liquid phase:

$$ \vec{F}_{drag} = -\frac{3}{4} C_D \frac{\rho_L \alpha_G}{d_B} |\vec{u}_G - \vec{u}_L| (\vec{u}_G - \vec{u}_L) \quad (1) $$

Here, u refers to velocity, ρ to density, α to volumetric fraction, d_B is the bubble diameter, and the subscripts L and G denote the liquid and the gas respectively. Different correlations exist in the literature for the estimation of the drag coefficient C_D (e.g. Schiller–Naumann, Ishii-Zuber, Tomiyama, etc.) (Khan et al., 2020).

The lift force represents the lateral force experienced by a bubble in a shear flow, and can be expressed as:

$$ \vec{F}_{lift} = -C_L \rho_L \alpha_G (\vec{u}_G - \vec{u}_L) \times \text{rot} (\vec{u}_L) \quad (2) $$

Different models for the lift coefficient C_L from the literature will be tested in this work (Tomiyama, Saffman-Mei, Moraga, Legendre-Magnaudet, etc.) (Tomiyama et al., 2002; Legendre and Magnaudet, 1998; Moraga et al., 1999; Mei, 1992). While most of these models have been developed from experimental or theoretical work on bubbles in laminar flow, Salibindla et al. (2020) have proposed a model taking into account the influence of turbulence on the deformation, and thus on the lift force on bubbles.<fnref n="1" /> As we will see, the application conditions of their model are well met in our experimental data, therefore particular attention will be given to this model in attempting to explain the insufficiency of the classical models in high-liquid-velocity pipe flow.

The virtual mass force accounts for the resistance to the acceleration of the liquid mass when a bubble is accelerated:

$$ \vec{F}_{VM} = -C_{VM} \rho_L \alpha_G \left( \frac{D_G \vec{u}_G}{Dt} - \frac{D_L \vec{u}_L}{Dt} \right) \quad (3) $$

Here, D_G/Dt and D_L/Dt are used to indicate material derivatives of the velocity of the gas and liquid phases, respectively. In the literature, there is a consensus of using a fixed value of 0.5 for the virtual mass coefficient C VM in bubbly flows (Rzehak et al., 2017).

The turbulent dispersion force refers to the influence of turbulent liquid velocity fluctuations on the bubbles. We employ the model proposed by Burns et al. (2004):

$$ \vec{F}_{TD} = -\frac{3}{4} C_D \frac{\alpha_G}{d_B} |\vec{u}_G - \vec{u}_L| \frac{\eta_L^{turb}}{\sigma_{TD}} \left( \frac{1}{\alpha_G} + \frac{1}{\alpha_L} \right) \nabla \alpha_G \quad (4) $$

Here, η_L^{turb} is the turbulent viscosity of the liquid, and σ_{TD} is the turbulent dispersion Schmidt number.

Finally, the wall lubrication force describes the effect of the presence of the wall on a bubble moving close to the wall. According to Antal it is expressed as follows (Antal et al., 1991):

$$ \vec{F}_{WL} = -C_{WL} \alpha_G \rho_L |\vec{u}_L - \vec{u}_G|^2 \vec{n}_W \quad (5) $$

Here, ⃗n_W is the unit normal vector pointing away from the wall. In Antal's model, the wall lubrication coefficient C_{WL} is computed as:

$$ C_{WL} = \max \left( 0, \frac{C_{W1}}{d_B} + \frac{C_{W2}}{y_W} \right) \quad (6) $$

with d_B the bubble diameter, y_W the distance from the wall, and C_{W1} and C_{W2} model constants, with the default values of -0.01 and 0.05 respectively. It is worth noting that this wall force is only active close to the wall, up to a cut-off distance of y_W ≤ (C_{W2}/C_{W1}) d_B. A plot of C_{WL} with distance from the wall is given in the supplementary material SM.12.

A large number of models have been developed for each of the different coefficients used in the expressions of the previously mentioned interfacial forces, and there is currently no consensus with respect to a universal combination of models (Khan et al., 2020). At the HZDR in Germany, a considerable amount of work has been dedicated to finding and validating a baseline model for bubbly flows, which would perform best in a large number of test cases. In this work, we will use this baseline model for bubbly flows (Rzehak et al., 2015; Kriebitzsch and Rzehak, 2016; Rzehak et al., 2017), as a starting set of closure terms. It is summarised in Table 2. We note that more recent versions of the baseline model included the Hosokawa model (instead of the Antal model) for the wall lubrication and the Hessenkemper model (instead of the Tomiyama model) for the lift force. However, a preliminary study showed that these two models provided worse or similar accuracy in the cases studied in this work (see supplementary material SM.2). Moreover, the Hosokawa model does not include a tuning coefficient, and thus cannot be tuned for finding optimal coefficient values for the

<footnote n="1">¹ We note that Salibindla et al. also propose a drag model in their work, but this model brought only negligible difference to the profiles obtained with the baseline model (see SM.13) and is therefore not discussed in this work.</footnote>
<table><caption>**Table 2**
Summary of the interfacial force models in the baseline model for bubbly flows from HZDR.</caption><thead><tr><th>Force</th><th>Model</th><th>References</th></tr></thead><tbody><tr><td><b>F&#x0302;<sub>drag</sub></b></td><td>Ishii and Zuber</td><td>Ishii and Zuber (1979)</td></tr><tr><td><b>F&#x0302;<sub>lift</sub></b></td><td>Tomiyama</td><td>Tomiyama et al. (2002)</td></tr><tr><td><b>F&#x0302;<sub>VM</sub></b></td><td>C<sub>VM</sub> = 0.5</td><td>Crowe et al. (2011)</td></tr><tr><td><b>F&#x0302;<sub>TD</sub></b></td><td>Burns<br>σ<sub>TD</sub> = 0.9</td><td>Burns et al. (2004)</td></tr><tr><td><b>F&#x0302;<sub>WL</sub></b></td><td>Antal<br>C<sub>W1</sub> = -0.01<br>C<sub>W2</sub> = 0.05</td><td>Antal et al. (1991)</td></tr></tbody></table>

types of flow studied here. Therefore, the version of the HZDR baseline model with the Antal and the Tomiyama interfacial forces was retained in this work.

The results presented further in this work will show that, while the gas velocity is generally predicted with satisfactory accuracy by the HZDR baseline model, the greatest challenge in our simulations will be of accurately predicting the radial void fraction distribution in vertical upward bubbly pipe flow. Some preliminary theoretical considerations guided us towards focusing on optimising the models for specific interfacial forces. The lift force is particularly influential as it governs the lateral migration of bubbles, which directly impacts their radial distribution. Additionally, the wall lubrication force plays a critical role by preventing bubbles from attaching to the wall, thereby influencing their spatial arrangement within the pipe. On the other hand, while the drag force primarily affects the vertical velocity of bubbles and the virtual mass force is more relevant in accelerating flows, their influence on the radial distribution in steady vertical flows is expected to be less pronounced. Finally, the turbulent dispersion force does have a significant impact on the radial distribution of bubbles. Nevertheless, its nature is diffusive, acting to spread bubbles from regions of high void fraction to regions of low void fraction, thereby flattening peaks in the void fraction profile. Unlike the lift force, turbulent dispersion does not have a preferred direction and thus does not actively drive bubbles in a specific lateral direction. Given the objective of explaining the different experimentally observed radial distributions, characterised by lateral shifts in void fraction peaks, priority is given to optimising the lift force model rather than the turbulent dispersion model. This approach enables an improved capture of the directional influences that shape the radial void fraction distribution.

Therefore, the interfacial forces primarily considered in this work are the lift force and the wall lubrication force. Various models and various parameter values for these interfacial forces will be tested and discussed in this work.

## 3.2. Geometry and mesh

For most of the cases studied here, the geometry consists of a vertical pipe. For the Ryan and Kim (2023) cases, we first built a 3D O-grid mesh, represented in Fig. 3(b). The gas is injected on a circular ring with the inner and outer radii equal to 0.5 and 0.8 times the pipe radius. The liquid is injected from both the inner circle and the outer ring. The simulations on this mesh were compared with axisymmetric simulations on a 2D wedge mesh (opening angle 2.5°) with similar cell sizing in the radial and axial directions (see supplementary material SM.3). The results were found to be identical, hence the 2D wedge mesh was further used for most of the simulations in this work (Fig. 3(a)). Similar 2D meshes, with similar cell sizes, were used also for the cases from Worosz et al. (2015) and Qiao et al. (2017) (inlet types A and B) and Kashinsky and Randin (1999).

The only exception was the Qiao et al. (2017) case with inlet type C, which includes a horizontal pipe section and an elbow. Since in this case the axisymmetrical assumption cannot be used, we used half of a

3D mesh, with a symmetry assumption on the cutting plane (Figs. 3(c) and 3(d)).

To ensure the grid independence of the numerical results, a grid convergence study was conducted for the case R1, and the results are detailed in the supplementary material SM.4.

Moreover, the y+ value at the walls was found to be between 30 and 80 for all cases. The grid convergence study showed that the results were independent of the resolution at the wall (y+ value) within this range.

## 3.3. Boundary conditions and numerical parameters

At the liquid inlet, we impose the same value $u_{\text{liquid,inlet}}$ for both $u_L$ and $u_G$, whereas the void fraction is set to $\alpha_G = 10^{-6}$. The pressure gradient at both the liquid and gas inlets is imposed in order to maintain the prescribed flux (fixedFluxPressure BC). The turbulent kinetic energy at the inlet is estimated according to the standard method:

$$k_{\text{inlet}} = \frac{3}{2} (I_t \cdot |u_L|)^2 \quad (7)$$

For the turbulent intensity $I_t$, we prescribe a value of 5%, and $|u_L|$ is the superficial liquid velocity in the respective simulation case. The specific turbulent dissipation rate at the inlet is estimated as follows:

$$\omega_{\text{inlet}} = \frac{k_{\text{inlet}}^{1/2}}{C_{\mu} L_{t}} \quad (8)$$

Here, we use the standard value of 0.09 for the shear-induced turbulence coefficient $C_\mu$ and a value of 10% of the pipe diameter for the turbulent length scale $L_t$ (Kriebitzsch and Rzehak, 2016).

At the gas inlet (on the ring of cells between 0.5 and 0.8 of the radius), the void fraction $\alpha_G$ is set to 0.2 (injection of bubbles), $u_L$ is set to 0 and $u_G$ is set according to the gas flux $j_{G,inlet}$. Since the air pressure and density vary with the position in the vertical pipe, the air volumetric flux $j_G$ is not constant along the pipe length. The value of $j_{G,inlet}$ is calculated from the value measured at the measurement port $j_{G,loc}$ located farthest from the inlet, by keeping a constant mass flux and taking into account the pressure difference.

In the case Q1B, in order to simulate the effect of a flow straightener, both inlet regions were merged into a single uniform inlet region. The gas void fraction was set to 0.2, while the liquid velocity and gas velocity were adapted in order to obtain the required flow rates for both phases.

At the wall, we apply a no-slip condition for the liquid velocity, a free-slip condition for the gas velocity, zero-gradient for pressure and void fraction, and use wall functions for the turbulent variables.

Finally, at the outlet we apply an “inletOutlet” condition (which imposes generally a zero-gradient condition and switches to a fixed value in the case of backflow) for $\alpha_G$, $u_G$ and $u_L$, and a fixed value (“prghPressure”) for the pressure.

The temperature of the fluids was set to 20 °C at the inlet, with adiabatic walls.

For the bubble size distribution, previous research (Kriebitzsch and Rzehak, 2016) recommends using a monodisperse assumption when the bubble diameter $d_B$ is smaller than 5 mm, which is the case in all the experimental cases used for the simulations in this work. This assumption is further justified by the observation that the bubble size distributions are relatively narrow in the experiments (see supplementary material SM.10). The constant value for $d_B$ was set as the average bubble diameter measured at the farthest port from the inlet. While it would be desirable to have a model that can accurately predict bubble diameters without relying on experimental data, incorporating such a model would have added significant complexity to this initial phase of the study. Our current focus is on developing a model that provides accurate predictions of the radial distributions of two-phase flow properties, such as void fraction and gas velocity, in developed flow across the different geometries relevant to loop reactors (riser,
<figure><img src="image_5.png" /><figcaption>Fig. 3. Meshes used in this work.</figcaption></figure>

downcomer). Future work will be devoted to validating a predictive model for bubble sizes.

A more detailed overview of the other numerical parameters and of the numerical schemes that were used inside the multiphaseEulerFoam module is given in the supplementary material SM.9.

# 4. Results

## 4.1. Upward flow

### 4.1.1. Structure of the flow

Let us first examine the experimental data for one of the upward flow cases at high liquid velocity (case R3), in order to understand the structure of the flow and the distribution of the phases. In Figs. 4(a) and 4(c), we plot the experimental distributions of void fraction and gas velocity for case R3 from Ryan and Kim (2023), at the position z = 110 D, where the two-phase flow is expected to be fully developed. Experimental profiles were available only for two azimuthal angles, distanced by 22.5°. The data was interpolated linearly between these two radii.

The first observation we make is that the void fraction values are quite low (the maximum value is around 0.1), leading to the interpretation that we have a bubbly flow pattern.<fnref n="2" /> Most of the bubbles are concentrated in the region between 0.4 *R* and 0.8 *R*, with a maximum at around 0.65 *R*. While a significant part of the bubbles also flow in the central region of the pipe (void fraction of around 7%), very few bubbles are found close to the wall, beyond 0.9 *R*. Regarding the gas velocity distribution, we see a velocity peak in the centre of the pipe (roughly 5 m/s), decreasing very gently until 0.8-0.9 *R*, and dropping drastically after 0.9 *R*. Finally, both distributions present only a very small dependence on the azimuthal angle.

In the discussion of the velocity profiles, it is important to consider the impact of the probe size on the accuracy of the measurements. In high superficial velocity flows, the probe may not fully resolve the shear layer near the wall unless special care is taken to miniaturise or align the probe. This limitation can affect the precision of the experimental data, particularly in regions where steep velocity gradients are present. Therefore, while the experimental profiles provide valuable insights

<footnote n="2">² This can be compared e.g. with the flow pattern map by Griffith and Walis, where bubbly flow is found for void fractions below 0.1, or corroborated with the flow pattern map proposed by Ryan et al. (2023).</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Case R3: contour plot of the void fraction and gas velocity distributions in the cross-section of the pipe (only 22.5° of the cross-section are shown) at the position z = 110 D, obtained from the experimental data (Ryan and Kim, 2023) and from the simulation with the HZDR baseline model.</figcaption></figure>

into the flow structure, the results near the wall should be interpreted with caution due to potential measurement uncertainties.

For comparison, in Figs. 4(b) and 4(d), we show the predictions obtained with the baseline HZDR model for the same case R3. In the void fraction distribution we notice a very different structure of the two-phase flow than in the experiment. While the void fraction values remain relatively low (maximum of around 0.15) and correspond to a bubbly flow regime, it appears that the majority of the bubbles are concentrated in the vicinity of the wall, between 0.8 *R* and 0.9 *R*. The maximum concentration of bubbles is about 50% larger than in the experiment, and is located at around 0.9 *R*. Almost no bubbles flow in the central pipe region and up to 0.6 *R*, and the void fraction starts to increase gradually after 0.6 *R*.

In the velocity distribution, we notice a larger degree of similarity to the experiment. We notice a central peak of around 5 m/s and a decrease at a modest rate until 0.95 *R*. The drop in gas velocity is more abrupt than in the experiment and takes place mostly between 0.95 *R* and the wall.

Finally, the simulation results are axisymmetric, which is in relatively good agreement with the behaviour observed in the experiment. Given this axisymmetrical distribution of the flow in both the experiment and the simulation, in the following we will represent the void fraction and gas velocity distributions using radial profiles, rather than contour plots.

In this first analysis we notice that the baseline HZDR model fails to predict the distribution of the phases in the upward pipe flow at high liquid velocities. In the following, we will compare the measured and simulated distributions for all the upward flow cases and further assess the performance of the baseline HZDR model.

### 4.1.2. Results with the baseline HZDR model for all upward cases

In Fig. 5, the numerical results using the baseline model of HZDR for bubbly flows in the upward cases (R1-R3 - high liquid velocities; R4, R5 and W1 - low liquid velocities) are compared with experimental

void fraction and gas velocity profiles from<fnref n="3" /> Ryan and Kim (2023) and Worosz et al. (2015). The normalised mean absolute error (nMAE) for each simulated profile is reported in the figure caption. Other error metrics were calculated for all the simulated profiles and they are reported in the supplementary material SM.11.

The gas velocity is predicted with a reasonably good accuracy, with a nMAE of the order of 10% to 22%, which is in the same order of magnitude as the uncertainty of these measurements. However, the void fraction prediction is rather poor, especially in the cases with higher liquid velocities (nMAE higher than 75% in cases R1, R2 and R3). It is worth noting that the HZDR baseline model has not been validated for such high liquid velocities (Liao et al., 2019).

Let us discuss the different cases in the order of increasing *u*<sub>L</sub>: in case R4, the baseline model fails to predict both the position and the height of the void fraction wall peak, while it largely overpredicts the values of void fraction towards the centre of the pipe. In case W1, the void fraction distribution towards the centre of the pipe is predicted correctly, but the height of the wall peak is largely underpredicted. Case R5 is where the baseline model performs best: the void fraction distribution in the central region of the pipe and the height of the wall peak are predicted correctly, but the predicted peak position is slightly closer to the wall than in the experiments. In cases R1, R2 and R3 (high liquid velocities), the peak in the measured void fraction profiles shifts more towards the centre of the pipe. However, the baseline model continues to predict a peak at around *r*/*R* = 0.9, and the peak height is largely overpredicted. Moreover, the void fraction in the central zone of the pipe is grossly underpredicted.

We notice the following pattern in the experimental data, that is not well reproduced by the simulations: at lower liquid velocities, the wall peak in the radial void fraction profile is higher and closer to the wall.

<footnote n="3">³ The numerical and experimental profiles presented here and throughout the rest of the article, unless otherwise mentioned, are the profiles at the location of the measurement port farthest from the inlet.</footnote>
<figure><img src="image_7.png" /><figcaption>(a) Cases R4, R5 and W1 - void fraction. Normalised mean absolute errors (nMAE): R4 - 48%, R5 - 21%, W1 - 33%.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) Cases R4, R5 and W1 - gas velocity. nMAE: (d) Cases R1, R2 and R3 - gas velocity. nMAE: R1 - 14%, R2 - 17%, R3 - 12%.</figcaption></figure>

<caption>**Fig. 5.** Numerical results with the baseline HZDR model for the upward flow cases - comparison with experimental data.</caption>

At higher liquid velocities, the peak is shifted towards the centre of the pipe and looks flattened.

Besides the theoretical considerations on the nature of interfacial forces presented in Section 3.1, a preliminary numerical study was conducted on case R1, including the drag, lift, turbulent dispersion and wall lubrication models (see supplementary material SM.5) in order to select the interfacial forces that can be modified to reproduce the experimentally observed behaviour. It was found that the most relevant forces that determine the void fraction profile are the lift force and the wall lubrication. While the turbulent dispersion also has a significant effect on the radial void fraction profile, its diffusive nature (without a preferred direction of action) does not make it a good candidate in explaining the lateral shift of void fraction peaks observed experimentally. Therefore, the lift force (including the Salibindla model that takes turbulence into account) and the wall lubrication will be considered in the following model optimisation.

In the following, each case is analysed separately and the adjustments that can be made to the models to obtain better predictions of the void fraction profile are discussed.

### 4.1.3. Detailed discussion of each case and adjustment of the models

**Case R4.** In this case, one idea for improving the prediction of the void fraction distribution would be to reduce the wall force, because this force seems to push the wall peak too far from the wall and reduce its height. However, as it can be seen in Fig. 6(a), this does

not significantly improve the prediction. Reducing the wall force has the effect of moving the wall peak very close to or even immediately at the wall, which is known to be non-physical. The second idea would be to modify the value of the lift coefficient, which has a direct impact on the radial void fraction distribution: a higher lift force tends to push the small bubbles<fnref n="4" /> towards the wall and accentuates the wall peak; a smaller lift force leads to a more uniform distribution and to higher values of the void fraction at the centre of the pipe. In Fig. 6(b), we see that the deviation in the numerical results is indeed strongly linked to the lift force. We noticed that the Tomiyama lift model yields a lift coefficient of approximately 0.3, the Saffman-Mei lift model yields $C_L \approx 0.1$, the Moraga model $C_L \approx 0.25$ and the Legendre-Magnaudet model $C_L \approx 0.5$, in all the cases studied in this work.<fnref n="5" /> For brevity, only the Tomiyama lift model will be labelled as such in the plots (since it is part of the HZDR baseline model), while the others will be labelled by their predicted $C_L$ value.

The experimental data at low liquid velocities (case R4) is much better reproduced by a constant lift coefficient of 0.8 - nMAE drops to about 5%, which is in the range of the experimental uncertainty - than with the value of 0.3 calculated by the baseline model.

<footnote n="4">⁴ In all the cases studied here, bubbles are considered to be small, with less than 3 mm in diameter.</footnote>

<footnote n="5">⁵ See supplementary material SM.6.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Simulation cases R4 and W1 (low to medium liquid velocity). Effect of the wall lubrication and lift models on the void fraction distribution.</figcaption></figure>

Regarding the Salibindla lift model, the condition of applicability of the We-base $C_L$ correlation ($\epsilon \ge 0.5$) is only met in the region very close to the wall, see supplementary material SM.7. This model does not yield a more satisfactory prediction than the baseline model.

Finally, as expected, the lift coefficient proved to have very little or no effect on the gas velocity distribution.

**Case W1.** In case W1, the baseline model has a similar behaviour to case R4 close to the wall, but the void fraction value at the centre of the pipe is predicted correctly. Again, reducing the wall force has only a limited beneficial effect on the prediction, but leads to void fraction profiles with the wall peak immediately at the wall, as shown in Fig. 6(c). As in case R4, increasing the lift coefficient $C_L$ to 0.8 allows to obtain an accurate prediction of the bubbly flow close to the wall (Fig. 6(d)). However, the accurate prediction of the void fraction distribution in the central zone of the pipe is lost when increasing $C_L$. For this case,

we could not find any reasonable combination of interfacial models that would lead to an accurate prediction in both regions.<fnref n="6" />

Again, as it can be seen in the supplementary material SM.7, the Salibindla lift model in case W1 can only be applied very close to the wall. Like in case R4, this model yielded a worse void fraction prediction than the baseline model, as not even the overall shape of the profile is reproduced.

**Case R5.** Among the different configurations investigated here, this is where the baseline model has the best performance (nMAE of 21%). The wall peak height and the void fraction in the central region of the pipe are predicted with a satisfactory accuracy. However, the position of the wall peak is not well predicted. The position of the wall peak in the experimental data seems to shift slightly towards the pipe

<footnote n="6">⁶ We note that the nMAE metric may not be sufficient to evaluate the performance of a model for case W1, because we only have one experimental point in the void fraction wall peak, where actually most of the air flow occurs.</footnote>
<figure><img src="image_10.png" /><figcaption>(a) Case R5: Antal wall lubrication force - $C_{W2}$ coefficient. nMAE: baseline - 21%, $C_{W2}$ 0.075 - 10.5%, $C_{W2}$ 0.1 - 11%, $C_{W2}$ 0.15 - 34%.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) Case R1: Antal wall lubrication force - $C_{W2}$ coefficient. nMAE: baseline - 88%, ($C_{L}$ 0.3 - $C_{W2}$ 0.05) - 91%, ($C_{L}$ 0.3 - $C_{W2}$ 0.2) - 29%, ($C_{L}$ 0.3 - $C_{W2}$ 0.25) - 19%, ($C_{L}$ 0.3 - $C_{W2}$ 0.3) - 25%, ($C_{L}$ 0.1 - $C_{W2}$ 0.1) - 13%.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Case R1: Lift force - $C_{L}$ coefficient. nMAE: baseline - 88%, $C_{L}$ -0.23 - 100%, $C_{L}$ -0.025 - 29%, $C_{L}$ 0.1 - 30%, $C_{L}$ 0.3 - 91%, $C_{L}$ 0.5 - 112%, $C_{L}$ 0.8 - 125%, Salibindla - 79%.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(d) Cases R1, R2 and R3 with optimised $C_{L}$ and $C_{W2}$ coefficients. nMAE: (R1 - $C_{L}$ 0.1 - $C_{W2}$ 0.1) - 13%, (R1 - $C_{L}$ 0.3 - $C_{W2}$ 0.25) - 19%, (R2 - $C_{L}$ 0.1 - $C_{W2}$ 0.1) - 15%, (R2 - $C_{L}$ 0.3 - $C_{W2}$ 0.25) - 18%, (R3 - $C_{L}$ 0.1 - $C_{W2}$ 0.1) - 11%, (R3 - $C_{L}$ 0.3 - $C_{W2}$ 0.25) - 14%.</figcaption></figure>

<caption>**Fig. 7.** Simulation cases R5, R1, R2 and R3 (medium to high liquid velocity). Effect of the wall lubrication and lift models on the void fraction distribution.</caption>

centre, compared to cases R4 and W1 that had a lower liquid velocity. Therefore, we attempt an adjustment of the wall force. In Fig. 7(a) we see that an increase of about 50% in the $C_{W2}$ wall force coefficient allows to obtain an accurate void fraction profile with the numerical model (nMAE drops to 10.5%).

**Cases R1, R2 and R3.** The superficial liquid velocity is equal in these three cases, which seems to be the reason for the similarity between their void fraction profiles. We will focus our analysis on case R1, and then test the improved models on the other two setups. Fig. 7(b) indicates that the lift force can only partially account for the shape of the void fraction distribution in these high-liquid-velocity cases. Indeed, at smaller $C_{L}$, a flattened distribution is obtained, which seems to agree more with the experimental data.<fnref n="7" /> However, the peak of the void fraction in the experimental data is even more shifted towards

the centre of the pipe (at around $r/R = 0.6$) than we have seen in the previous case R5, while the numerical results seem to indicate that the lift force alone is not able to produce such a shift.

A few attempts were made to simulate the flow in case R1 using a negative lift coefficient. However, the void fraction profiles obtained with a negative $C_{L}$ have a fundamental difference compared to the experimental profile: they exhibit a monotonous increase from the wall towards the centre of the pipe, with a peak at the centre. Although with a very small negative $C_{L}$, the profile remains close to the experimental data, this discrepancy in shape and the failure to predict an eccentric peak do not make it a good candidate for an optimised model.

Finally, we note that the turbulence in this case is strong (according to the definition by Salibindla et al. (2020) $\epsilon \gtrsim 0.5$) in the entire pipe (see supplementary material SM.7), so that the Salibindla model is expected to apply to the whole flow. However, this model again yields a very inaccurate void fraction distribution (nMAE 79%).

In Fig. 7(c), we test different values of the wall lubrication coefficient $C_{W2}$, and also a combination of $C_{W2} = 0.1$ with $C_{L} = 0.1$ (the

<footnote n="7">⁷ On the choice of the values of lift coefficient tested here, see also supplementary material SM.6.</footnote>
<figure><img src="image_14.png" /><figcaption>Fig. 8. Contour plot of the void fraction distribution in the cross-section of the pipe (only 22.5° of the cross-section are shown) at the position 110 D, obtained in the simulation with the two optimised models.</figcaption></figure>

value $C_L = 0.1$ is predicted by the Saffman-Mei lift model, as shown in supplementary material SM.6). The wall force can indeed shift the wall peak towards the centre of the pipe, but it also produces an increase of the void fraction at the centre of the pipe. For values of $C_{W2}$ higher than ≈0.25, a secondary void fraction peak is produced at the centre of the pipe, which is not in agreement with the experimental observations. The two combinations that match best with the experimental data are:

* the Tomiyama lift model ($C_L = 0.3$) with a high wall force ($C_{W2} \approx 0.25$, five times larger than in the baseline model) - we will name it OM1 (optimised model 1).

* the Saffman-Mei lift model ($C_L = 0.1$) with an increased wall force ($C_{W2} \approx 0.1$, two times larger than in the baseline model) - we will name it OM2 (optimised model 2).

In Fig. 7(d) we compare the predictions of these two optimised models with the experimental data in cases R1, R2 and R3. Overall, although not yet perfect, we see a net improvement of the predicted void fraction profiles, compared to the baseline model (nMAE in the range of 13 to 15%, in the same order of magnitude with the experimental uncertainty of 10%). The contour plots of void fraction distribution obtained with the two optimised models are shown in Fig. 8. Comparing these with the experimental data from Fig. 4(a), we see that the overall distribution of the phases and the void fraction values are predicted with a significantly higher accuracy than with the baseline HZDR model.

## 4.2. Downward flow

In the downward cases, we assess the performance of the baseline model and of the optimised models found in the previous section, by comparison with the experimental data. At low liquid velocities (cases K1, K2 and K3), the baseline model provides a reasonably good prediction of the void fraction distribution and of the liquid velocity, compared to experimental data from Kashinsky and Randin (1999) (Fig. 9). In the cases with higher gas flux (K2 and K3), the void fraction is underpredicted in the central part of the pipe, but the normalised mean absolute error does not exceed 20%, which can be considered satisfactory considering the simplifying assumptions used. On the other hand, in the cases with lower gas flux (K1 and K2), the simulated liquid velocity profiles deviate from experimental data at $0.4 < r/R < 0.95$, with a local relative error reaching the order of 15%. The “improved” case for low liquid velocities ($C_L = 0.8$) yields globally similar results to the baseline model for these cases, but with a slightly worse prediction in the region $0.7 < r/R < 1$.

Among the cases with higher liquid velocity (Q1 cases), we only discuss the cases with a flow straightener at the inlet (Q1B) and with a 90° elbow at the inlet (Q1C). The case without flow straightener (Q1 A) is more difficult to reproduce, since the geometry of the injector used in the experiments is not known. However, some results for the case

Q1 A are reported in the supplementary material SM.8, which show a good qualitative agreement with the experimental data.

For case Q1B, both the baseline model and the optimised models are able to accurately predict the gas velocity profile (Fig. 10(b), nMAE < 2%). When comparing the void fraction distribution, none of the models represents exactly the experimental pattern (Fig. 10(a)). The baseline model predicts well the height of the central void fraction peak, but underpredicts the void fraction close to the pipe walls, whereas the OM2 underpredicts the void fraction at the pipe centre, but gives a more accurate prediction of the profile in the central region and close to the wall. The prediction by the OM1 model is not accurate in either region of the flow.

In the simulations of Case Q1C, both the void fraction profile (Fig. 10(c)) and the gas velocity profile (Fig. 10(d)) exhibit an asymmetrical distribution with an eccentric peak, even after 67.5 diameters of the downcomer. In contrast, the experimental data indicates that the flow becomes symmetrical again before reaching 67.5 diameters. Among the three models tested, the OM2 model shows the least divergence from the experimental void fraction distribution.

The discrepancy between the simulation results and the experimental observations can be attributed to the flow pattern in the horizontal section, which may differ from a dispersed bubbly flow. Specifically, the simulation revealed an accumulation of gas at the top of the horizontal pipe, indicative of flow stratification. The models employed in this study are not well-suited for such stratified flow conditions, and consequently, they struggle to accurately predict the transition back to a dispersed bubbly flow pattern in the vertical section, which is observed experimentally.

## 5. Discussion

The baseline model is a good starting basis for the simulation of disperse bubbly flows. However, as it has been validated only for vertical pipe flows with liquid velocities up to 1 m/s, it might need adjustments for cases with higher liquid velocities, as highlighted by the simulation results in the present work. In particular, the results in this work showed that the void fraction profiles at high liquid velocities are not well predicted by the baseline model.

At higher liquid flow rates, there seems to occur a shift of the void fraction peak position towards the centre of the pipe, as well as a flattening of this peak. Theoretically, the interfacial forces that can act on the bubbles in the lateral direction and significantly modify the void fraction profile are the lift, turbulent dispersion, and wall lubrication (only close to the wall) forces. Our results pointed out that it was a combination of lift force and wall lubrication that best accounted for these void fraction profiles.

The shift in the void fraction peak could be relatively well captured by two optimised models, OM1 and OM2, which are summarised in Table 3. OM1 implies the use of a very large wall lubrication force (five times larger than in the baseline model), which is active up to
<figure><img src="image_15.png" /></figure>

(a) Void fraction distribution. nMAE: (K1 - baseline) - 23%, (K1 - $C_L$ 0.8) - 35%, (K2 - baseline) - 16%, (K2 - $C_L$ 0.8) - 25%, (K3 - baseline) - 17%, (K3 - $C_L$ 0.8) - 19%.

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

(b) Dimensionless liquid velocity distribution. The liquid velocity is made dimensionless with the value at the centre of the pipe. nMAE: (K1 - baseline) - 17%, (K1 - $C_L$ 0.8) - 17%, (K2 - baseline) - 13%, (K2 - $C_L$ 0.8) - 13%, (K3 - baseline) - 14%, (K3 - $C_L$ 0.8) - 14%.

<caption>**Fig. 9.** Simulation cases K1, K2 and K3: performance of the baseline model and using a constant lift coefficient $C_L = 0.8$.</caption>

<table><caption>**Table 3**
Summary of the interfacial force models in the optimised models OM1 and OM2.</caption><thead><tr><th>Force</th><th>Model OM1</th><th>Model OM2</th></tr></thead><tbody><tr><td><b>F&#x0302;<sub>drag</sub></b></td><td>Ishii and Zuber</td><td>Ishii and Zuber</td></tr><tr><td><b>F&#x0302;<sub>lift</sub></b></td><td>$u_{SL} < 2.5$: const. $C_L$ following Fig. 11<br>$u_{SL} \ge 2.5$: Tomiyama ($C_L \sim 0.3$)</td><td>$u_{SL} < 2.5$: const. $C_L$ following Fig. 11<br>$u_{SL} \sim 2.5$: Tomiyama ($C_L \sim 0.3$)<br>$u_{SL} > 2.5$: Saffman-Mei ($C_L \sim 0.1$)</td></tr><tr><td><b>F&#x0302;<sub>VM</sub></b></td><td>$C_{VM} = 0.5$</td><td>$C_{VM} = 0.5$</td></tr><tr><td><b>F&#x0302;<sub>TD</sub></b></td><td>Burns<br>$\sigma_{TD} = 0.9$</td><td>Burns<br>$\sigma_{TD} = 0.9$</td></tr><tr><td><b>F&#x0302;<sub>WL</sub></b></td><td>Antal with $C_{W1} = -0.01$<br>$u_{SL} \le 2.5$: $C_{W2}$ following Fig. 11<br>$2.5 < u_{SL} < 4$: $C_{W2}$ rising linearly<br>$u_{SL} \sim 4$: $C_{W2} = 0.25$</td><td>Antal with $C_{W1} = -0.01$ $C_{W2}$ following Fig. 11</td></tr></tbody></table>

a large distance from the wall. From a physical point of view, it is hard to justify that the wall lubrication force is active very far from the wall. The OM2 model combines a decrease of the lift force with a slight increase of the wall force (two times larger than in the baseline model) at higher liquid velocities (Fig. 11). Compared to the OM1, the OM2 model implies a much lesser increase in the wall force and in its active range, making it physically more reasonable. Moreover, the combination of results from upflow and especially downflow cases showed that the OM2 has a better and more robust performance than OM1. Therefore, in the following we will focus our discussion on the OM2 model.

While the modifications of lift and wall lubrication forces in the OM2 model have proved to significantly improve the accuracy of the baseline model in high-liquid-velocity upward bubbly flows, more work needs to be done for a thorough understanding and modelling of the underlying phenomena. Most notably, the shift of the void fraction peak towards the centre of the pipe could not be captured by any of the classical lift models. One possible underlying cause is that all of these models were developed based on experiments or theoretical developments in laminar flow, which is more representative for bubble columns or low-liquid-velocity pipe flows. In our case, the liquid flow is highly turbulent, and the mentioned models fail to capture the effect of turbulence on the bubbles.

An interesting recent development in this direction was made by Salibindla et al. (2020). They showed that in highly turbulent flow, even small-sized bubbles (of 1–2 mm in diameter), like the ones in the cases investigated in the present work, are deformed by the turbulent stresses. As a result, the lift force on these bubbles is modified. Therefore, lift models that do not take turbulence into account are expected to be inaccurate in these flow regimes. The authors proposed a new lift model based on the Weber number $We$. Notably, this model predicts a possible change of sign in the lift force for small bubble diameters (1–2 mm) in highly turbulent conditions, as opposed to the classical models, where the change of sign occurs at much larger bubble sizes (6–7 mm for the Tomiyama model).

As generally the turbulence is highest close to the pipe wall, the Salibindla model would qualitatively result in a negative lift coefficient at the wall and possibly a positive lift coefficient in the centre of the pipe (assuming constant bubbly sizes). This would indeed push the void fraction peak away from the wall, towards the region where the lift coefficient changes sign. This provides a strong physical basis for explaining the position of the void fraction peak and its shift towards the centre in highly turbulent bubbly flows, without the need to invoke a less physical wall force that might act far away from the wall.

Therefore, the authors believe that a Salibindla-type lift model, accounting for the effect of turbulent stress on small bubbles and on the lift force, provides currently the most physically sound explanation to the void fraction distributions observed experimentally in high-liquid-velocity flows. Although the lift model proposed by Salibindla et al. (2020) did not offer good predictions in the cases investigated in this work, adjustments to this model would likely significantly increase its performance in high-velocity upward pipe flow. More work is needed in this direction in order to propose adjustments for the Salibindla model in vertical upward pipe flows.

We note that other authors have worked on the development of new lift models that take turbulence effects into account (Sugrue, 2017). However, these models still prove to be insufficiently mature (Shiea et al., 2019).

In the current work, satisfactory engineering predictions for high-velocity cases were obtained by increasing the wall force and adjusting the lift coefficient (optimised model OM2). While these modifications improve the predictive power of the baseline model from an engineering point of view (for the study of PLRs), the authors believe that the
<figure><img src="image_17.png" /><figcaption>(a) Case Q1B: Void fraction distribution. nMAE: baseline - 28%, ($C_L$ 0.1 - $C_{W2}$ 0.1) - 30%, $C_{W2}$ 0.25 - 54%.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(c) Case Q1C: Void fraction distribution. nMAE: baseline - 73%, ($C_L$ 0.1 - $C_{W2}$ 0.1) - 35%, $C_{W2}$ 0.25 - 97%.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) Case Q1B: Gas velocity distribution. nMAE: baseline - 1%, ($C_L$ 0.1 - $C_{W2}$ 0.1) - 1%, $C_{W2}$ 0.25 - 2%.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(d) Case Q1C: Gas velocity distribution. nMAE: baseline - 7%, ($C_L$ 0.1 - $C_{W2}$ 0.1) - 7%, $C_{W2}$ 0.25 - 6%.</figcaption></figure>

<caption>**Fig. 10.** Simulation cases Q1B and Q1C: performance of the baseline model and of the optimised models for high liquid velocity.</caption>

increase of wall force and of its range of action is difficult to justify on a physical basis. The same effects on the flow would probably be obtained with an appropriate model for the lift coefficient (similar to the Salibindla model), but currently no such model is available. Future work will be dedicated to improving lift coefficient models that take into account the effect of turbulence, in the case of high velocity vertical pipe flows.

Another improvement direction would be the addition of a model to predict the evolution of bubble sizes, such as IATE (Kim et al., 2021) or PBM (Liao et al., 2015; Lehnigk et al., 2022). This would be a step towards obtaining models that are predictive also in the flow developing regions.

## 6. Conclusions

We have performed an assessment of two-fluid models in Open-FOAM for the simulation of upward and downward bubbly pipe flows with medium-to-high liquid velocities in vertical pipes, which are relevant in a variety of industrial contexts. Our focus in this paper is on (fast) chemical and biochemical loop reactors. The simulations were

compared with measured data from four different experimental setups from the literature.

In the upward flow cases, the baseline model for bubbly flows correctly calculates the gas velocity distribution. However, it proved inaccurate in predicting the radial distribution of the void fraction. The largest discrepancy originated from an experimentally observed shift in the position of the void fraction peak towards the centre of the pipe at higher liquid velocities, that is not correctly captured by the numerical model. It was found that the combination of modified lift and wall lubrication coefficients was the most efficient way to enhance the radial distribution of the void fraction profile and its peak. A relatively good agreement with the experimental data was found by decreasing the lift coefficient and increasing the $C_{w2}$ coefficient of the Antal wall force as the liquid velocity increases, which implies that the wall force is active on a larger distance from the wall at higher velocities. While this effectively improved the predictions of the baseline model, a plausible physical explanation for this behaviour probably lies in the effect of turbulence on the lift force at high liquid velocities, as pointed out by Salibindla et al. (2020). However, it is important to note that the Salibindla model is specifically developed for contaminated bubbles, and thus its conclusions, while important, are not directly applicable
<figure><img src="image_21.png" /><figcaption>**Fig. 11.** Proposed variation of the baseline model for bubbly flows at medium and high liquid velocities (OM2).</figcaption></figure>

to our current study. Because the model proposed by Salibindla et al. did not produce better predictions for the cases investigated here, more work is needed to encompass the high-velocity vertical pipe flow cases.

For the downward flow cases, a strong dependence of the void fraction profile on the inlet conditions was found. This is in good agreement with experimental observations by Qiao et al. (2017). Additionally, good quantitative agreement was also found in the case with an inlet flow straightener, which evolves more rapidly to a developed flow. The optimised model OM2 performed comparably to our baseline model in the examined downward flows.

The findings of this study offer a significant advancement regarding the simulation of high-velocity vertical pipe bubbly flows, which are critical for the design and optimisation of chemical loop reactors. By proposing an optimised model that adjusts the lift and wall lubrication coefficients, we have demonstrated a marked improvement in predicting void fraction distributions, particularly at higher liquid velocities where traditional models fall short. This optimised approach not only enhances the accuracy of simulations for engineering predictions but also underscores the limitations of existing models in capturing the complex dynamics of bubbly flows under such conditions. Furthermore, our work highlights the pivotal role of turbulence in influencing the lift force at high liquid velocities, aligning with the insights provided by Salibindla et al. (2020). This emphasis on turbulence-induced effects paves the way for future refinements in the modelling of the lift force, ultimately benefiting both academic research and industrial applications in chemical engineering.

Nevertheless, further work is required to extend the available experimental intervals to higher liquid velocities relevant to loop reactors. As well, it would be desirable to implement and validate models that can predict bubble sizes and validate them for PLR-relevant flow regimes, instead of depending on experimental bubble data. Such developments will allow us to finely study the shift in the position of the void's fraction profile (wall) peak for increasing liquid velocity. Moreover, such undertakings must be complemented by theoretical work to enhance the existing bubbly flow models and to extend their applicability to a wider range of liquid velocities.

### CRediT authorship contribution statement

**M. Badalan:** Writing – original draft, Visualization, Validation, Methodology, Investigation, Formal analysis. **M. Shiea:** Writing – review & editing, Validation, Formal analysis. **Y. Haroun:** Writing – review & editing, Supervision, Methodology, Conceptualization. **G.E. Vázquez-Angulo:** Visualization, Investigation. **J.E.V. Guzmán:** Writing – review & editing, Supervision, Conceptualization.

### 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

This work was supported by UNAM Posdoctoral Program (POSDOC). The authors wish to thank Dr. D. Ryan and Dr. S. Kim for kindly sharing their experimental data (Ryan et al., 2023; Ryan and Kim, 2023), which significantly contributed to this study.

### Appendix A. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.ched.2025.08.001.

### References

Antal, S.P., Lahey, R.T., Flaherty, J.E., 1991. Analysis of phase distribution in fully developed laminar bubbly two-phase flow. Int. J. Multiph. Flow 17 (5), 635–652.

Becker, S., Sokolichin, A., Eigenberger, G., 1994. Gas–liquid flow in bubble columns and loop reactors: Part II. Comparison of detailed experiments and flow simulations. Chem. Eng. Sci. 49 (24), 5747–5762.

Blenke, Heinz, 1979. Loop reactors. In: Ghose, T.K., Blakebrough, N., Fiechter, A. (Eds.), Advances in Biochemical Engineering, Volume 13. Springer, Berlin, Heidelberg, pp. 121–214.

Burlutskii, E., Di Felice, R., 2019. Experimental and numerical study of two-phase flow mixing in gas–liquid external-loop airlift reactor. Int. J. Multiph. Flow 119, 1–13.

Burns, Alan D., Frank, Thomas, Hamill, Ian, Shi, Jun-Mei, et al., 2004. The Favre averaged drag model for turbulent dispersion in Eulerian multi-phase flows. In: 5th International Conference on Multiphase Flow. In: 4, number 392, ICMF, pp. 1–17.

Crowe, Clayton T., Schwarzkopf, John D., Sommerfeld, Martin, Tsuji, Yutaka, 2011. Multiphase Flows with Droplets and Particles, second ed. CRC Press, Boca Raton.

Direct, Helmholtz-Zentrum Dresden-Rossendorf (HZDR) and CFD, Disperse Multiphase Flows in OpenFOAM | OpenFOAM.

Drew, Donald A., Passman, Stephen L., 1999. Theory of Multicomponent Fluids. In: Applied Mathematical Sciences, vol. 135, Springer, New York, NY.

Gao, Yongxiang, Hong, Du, Cheng, Youwei, Wang, Lijun, Li, Xi, 2019. CFD simulation for up flow jet-loop reactors by use of bi-dispersed bubble model. Chem. Eng. Res. Des. 141, 66–83.

Hibiki, Takashi, Ishii, Mamoru, Xiao, Zheng, 2001. Axial interfacial area transport of vertical bubbly flows. Int. J. Heat Mass Transfer 44 (10), 1869–1888.

Huang, Qingshan, Yang, Chao, Yu, Gengzhi, Mao, Zai-Sha, 2010. CFD simulation of hydrodynamics and mass transfer in an internal airlift loop reactor using a steady two-fluid model. Chem. Eng. Sci. 65 (20), 5527–5536.

Ishii, Mamoru, Hibiki, Takashi, 2011. Thermo-Fluid Dynamics of Two-Phase Flow. Springer, New York, NY.

Ishii, M., Paranjape, S.S., Kim, S., Sun, X., 2004. Interfacial structures and interfacial area transport in downward two-phase bubbly flow. Int. J. Multiph. Flow 30 (7), 779–801, A Collection of Papers in Honor of Professor G. Yadigaroglu on the Occasion of his 65th Birthday.

Ishii, Mamoru, Zuber, Novak, 1979. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AIChE J. 25 (5), 843–855.

Kashinsky, O.N., Randin, V.V., 1999. Downward bubbly gas–liquid flow in a vertical pipe. Int. J. Multiph. Flow 25 (1), 109–138.

Khan, Irfan, Wang, Mingjun, Zhang, Yapei, Tian, Wenxi, Su, Guanghui, Qiu, Suizheng, 2020. Two-phase bubbly flow simulation using CFD method: A review of models for interfacial forces. Prog. Nucl. Energy 125, 103360.

Kim, Seungjin, Ishii, Mamoru, Kong, Ran, Wang, Guanyi, 2021. Progress in two-phase flow modeling: Interfacial area transport. Nucl. Eng. Des. 373, 111019.

Kojic, Predrag, Kojic, Jovana, Pezo, Milada, Krulj, Jelena, Pezo, Lato, Mirkov, Nikola, 2022. Numerical study of the hydrodynamics and mass transfer in the external loop airlift reactor. Chem. Ind. Chem. Eng. Q. 28 (3), 225–235.

Kriebitzsch, Sebastian, Rzehak, Roland, 2016. Baseline model for bubbly flows: Simulation of monodisperse flow in pipes of different diameters. Fluids 1 (3), 29, Number: 3 Publisher: Multidisciplinary Digital Publishing Institute.

Legendre, Dominique, Magnaudet, Jacques, 1998. The lift force on a spherical bubble in a viscous linear shear flow. J. Fluid Mech. 368, 81–126.

Lehnigk, Ronald, Bainbridge, William, Liao, Yixiang, Lucas, Dirk, Niemi, Timo, Peltola, Juho, Schlegel, Fabian, 2022. An open-source population balance modeling framework for the simulation of polydisperse multiphase flows. AIChE J. 68 (3), e17539.
Liao, Yixiang, Ma, Tian, Krepper, Eckhard, Lucas, Dirk, Fröhlich, Jochen, 2019. Application of a novel model for bubble-induced turbulence to bubbly flows in containers and vertical pipes. Chem. Eng. Sci. 202, 55-69.

Liao, Yixiang, Rzehak, Roland, Lucas, Dirk, Krepper, Eckhard, 2015. Baseline closure model for dispersed bubbly flow: Bubble coalescence and breakup. Chem. Eng. Sci. 122, 336-349.

Mei, R., 1992. An approximate expression for the shear lift force on a spherical particle at finite reynolds number. Int. J. Multiph. Flow 18 (1), 145-147.

Moraga, F.J., Bonetto, F.J., Lahey, R.T., 1999. Lateral forces on spheres in turbulent uniform shear flow. Int. J. Multiph. Flow 25 (6-7), 1321-1372.

Moran, Sean, Middleton, John C., Carpenter, Keith J., 2017. Loop reactors. In: Wiley-V.C.H. Verlag GmbH & Co. KGaA (Ed.), Ullmann's Encyclopedia of Industrial Chemistry. Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, pp. 1-7.

Moudoud, Nacereddine, Rihani, Rachida, Bentahar, Fatiha, Legrand, Jack, 2018. Global hydrodynamic of hybrid external loop airlift reactor: Experiments and CFD modelling. Chem. Eng. Process. - Process. Intensif. 129, 118-130.

Norwood, Donald D., 1966. Method and apparatus for the production of solid polymers of olefins. US3248179A.

Orejas, Joaquin A., 1999. Modelling and simulation of a bubble-column reactor with external loop: Application to the direct chlorination of ethylene. Chem. Eng. Sci. 54 (21), 5299-5309.

Orejas, Joaquin A., 2001. Model evaluation for an industrial process of direct chlorination of ethylene in a bubble-column reactor with external recirculation loop. Chem. Eng. Sci. 56 (2), 513-522.

Puiman, Lars, Abrahamson, Britt, Lans, Rob G.J.M. Van Der Haringa, Cees, Noorman, Henk J., Picioreanu, Cristian, 2022. Alleviating mass transfer limitations in industrial external-loop syngas-to-ethanol fermentation. Chem. Eng. Sci. 259, 117770.

Qiao, Shouxu, Mena, Daniel, Kim, Seungjin, 2017. Inlet effects on vertical-downward air-water two-phase flow. Nucl. Eng. Des. 312, 375-388.

Roy, S., Joshi, J.B., 2008. CFD study of mixing characteristics of bubble column and external loop airlift reactor. Asia-Pac. J. Chem. Eng. 3 (2), 97-105.

Ryan, Drew, Kim, Seungjin, 2023. Characterization of local parameters of inclined upward air-water bubbly two-phase flows in a round pipe. Nucl. Eng. Des. 410, 112378.

Ryan, Drew, Kong, Ran, Kang, David, Dix, Adam, Kim, Seungjin, Bian, Jiawei, 2023. Effects of pipe inclination on global two-phase flow parameters. Nucl. Technol. 209 (10), 1485-1494.

Rzehak, Roland, Krepper, Eckhard, Liao, Yixiang, Ziegenhein, Thomas, Kriebitzsch, Sebastian, Lucas, Dirk, 2015. Baseline model for the simulation of bubbly flows. Chem. Eng. Technol. 38 (11), 1972-1978.

Rzehak, Roland, Ziegenhein, Thomas, Kriebitzsch, Sebastian, Krepper, Eckhard, Lucas, Dirk, 2017. Unified modeling of bubbly flows in pipes, bubble columns, and airlift columns. Chem. Eng. Sci. 157, 147-158.

Salibindla, Ashwanth K.R., Masuk, Ashik Ullah Mohammad, Tan, Shiyong, Ni, Rui, 2020. Lift and drag coefficients of deformable bubbles in intense turbulence determined from bubble rise velocity. J. Fluid Mech. 894, A20.

Schlegel, Fabian, Bilde, Kasper Gram, Draw, Mazen, Evdokimov, Ilya, Hänsch, Susann, Kamble, Vikrant Vinayak, Khan, Haris, Krull, Benjamin, Lehnigk, Ronald, Li, Jiadong, Lyu, Hongmei, Meller, Richard, Petelin, Gasper, Kota, Sesi Preetam, Tekavcic, Matej, 2024. Multiphase code repository by HZDR for OpenFOAM foundation software. Language: eng.

Serizawa, Akimi, Kataoka, Isao, Michiyoshi, Itaru, 1992. Experimental data set no. 24: Phase distribution in bubbly flow. Multiph. Sci. Technol. 6 (1-4), 257-301.

Shiea, Mohsen, Buffo, Antonio, Baglietto, Emilio, Lucas, Dirk, Vanni, Marco, Marchisio, Daniele, 2019. Evaluation of hydrodynamic closures for bubbly regime CFD simulations in developing pipe flow. Chem. Eng. Technol. 42 (8), 1618-1626.

Sokolichin, A., Eigenberger, G., 1994. Gas-liquid flow in bubble columns and loop reactors: Part I. Detailed modelling and numerical simulation. Chem. Eng. Sci. 49 (24), 5735-5746.

Sugrue, Rosemary M., 2017. A Robust Momentum Closure Approach for Multiphase Computational Fluid Dynamics Applications (Thesis). Massachusetts Institute of Technology, Accepted: 2018-02-16T19:26:59Z.

Talvy, Samuel, Cockx, Arnaud, Liné, Alain, 2007. Modeling of oxygen mass transfer in a gas-liquid airlift reactor. AIChE J. 53 (2), 316-326.

Teli, Shivanand M., Mathpati, Channamallikarjun S., 2021. Experimental and numerical study of gas-liquid flow in a sectionalized external-loop airlift reactor. Chin. J. Chem. Eng. 32, 39-60.

Tomiyama, Akio, Tamai, Hidesada, Zun, Iztok, Hosokawa, Shigeo, 2002. Transverse migration of single bubbles in simple shear flows. Chem. Eng. Sci. 57 (11), 1849-1858.

Vial, Ch., Poncin, S., Wild, G., Midoux, N., 2002. Experimental and theoretical analysis of the hydrodynamics in the riser of an external loop airlift reactor. Chem. Eng. Sci. 57 (22-23), 4745-4762.

Worosz, Ted, Kim, Seungjin, Hoxie, Chris, 2015. Experiments in cap-bubbly two-phase flows for two-group IATE development. Nucl. Technol. 190 (3), 264-273.

Yeoh, G.H., Tu, J.Y., 2009. Computational Techniques for Multiphase Flows. Elsevier.