Euler-Euler numerical simulations of upward turbulent bubbly flows in vertical pipes with low-Reynolds-number model

Abstract

In this work, numerical simulations of upward turbulent bubbly flows in vertical pipes are conducted in the Euler-Euler framework with a low-Reynolds-number (LRN) model for liquid. It is found that the existing correlation for the drag coefficient of a single bubble in a shear flow, which has been successfully used along with the high-Reynolds-number (HRN) models, cannot be used with the LRN model. The reason is that drag forces of bubbles calculated using such correlation are enormous in near-wall cells (order of ${10}^{12}\mathrm{N}/{\mathrm{m}}^3$ ), which causes a divergence of numerical simulations with LRN. Therefore, a modified correlation for the drag coefficient of a single bubble in shear flow, that can be used successfully with the LRN model, has been proposed. The results of numerical simulations, performed with a new correlation for the drag coefficient, are analysed and compared to selected experimental measurements for different pipe diameters and different flow conditions of gas and liquid. It is shown that the largest effect of the application of the new correlation for the drag coefficient of a single bubble in shear flow in numerical simulations can be achieved on the reduction of the relative velocity between gas and liquid. The degree of this reduction depends on the pipe diameter and liquid volumetric flux.

Introduction

Turbulent bubbly wall-confined vertical flows have many important industrial applications, such as in gas-liquid chemical reactors, boilers, airlift pumps and oil production and transport. The presence of bubbles in these flows can significantly influence skin friction compared to single-phase liquid flow, increasing mixing and heat transfer between phases. Owing to these reasons these flows have been extensively studied in the last 40 years. However, due to the great complexity of bubble-liquid and bubble-bubble interactions and the presence of wall confinement, the study of these flows still represents a challenge for researchers.

Direct numerical simulation studies with interface tracking of every bubble provide the highest level of numerical description of turbulent bubbly flows (Lu and Tryggvason,^1,2 Adoua et al., ³ Bois and du Cluzeau ⁴ ). However, due to high computational expense, at the present, these studies are possible for relatively low Reynolds number flows with limited number of bubbles (around few hundreds at most).

Due to the reasonable computational price of large-scale computations, most numerical studies in the Euler-Euler framework for upward turbulent bubbly flows in vertical pipes (Hosokawa and Tomiyama, ⁵ Lubchenko et al., ⁶ Magolan et al., ⁷ Parekh and Rzehak, ⁸ Colombo et al., ⁹ Lopez de Bertodano et al., ¹⁰ Rzehak and Krepper,^11–14 Rzehak and Kriebitzsch, ¹⁵ Rzehak et al., ¹⁶ Troshko and Hassan ¹⁷ ) are conducted using high-Reynolds-number (HRN) models for liquid phase. In these models, the turbulence in the near-wall region is modelled using various empirical equations known as wall functions, and boundary conditions for the momentum and turbulence transport equations are specified somewhere in the near-wall region, and not at the wall itself.

There is no clear conclusion if wall functions in the context of two-phase flows should always be the same as in single-phase flows (Troshko and Hassan, ¹⁷ Colin et al. ¹⁸ ). Still, if a HRN model is used in Euler-Euler numerical studies of wall confined turbulent bubbly flows, wall functions used for single-phase flow are usually applied without any change.

Alternatives to HRN models are low-Reynolds-number (LRN) models, which fully resolve flow fields up to the wall. Since many physical quantities have strong gradients in near-wall region, the mesh resolution in LRN approach has to be high enough to capture these changes ( $y^{+}<1$ ). This leads to a significantly larger number of cells and higher computational costs compared to HRN approach, but with the increased accuracy in the prediction of flow physics, especially in near-wall region. Colombo et al. ⁹ pointed out that modelling of the near-wall region and bubble-induced turbulence (BIT) are identified as the areas where further development is needed .

In the upward turbulent bubbly flows, due to the buoyancy that affects the bubbles, air velocity is higher than water velocity. This difference in phase velocities modifies the turbulence intensity in the liquid phase. This can be taken into account by including source terms (BIT models) in the single-phase two-equation turbulence models. These source terms are modelled but, no consensus for their form has been reached yet (Magolan et al. ⁷ ). Rzehak and Krepper ¹³ have reported a comparison of different BIT models. They found that it is difficult to report general conclusions and that further research on this topic is needed. To reach a general solution, they proposed two directions: the first, the formation of a more complex BIT model, and the second, the introduction of a separate equation for the bubble-induced turbulent kinetic energy (Chahed et al. ¹⁹ ).

A problem that occupies the most attention, when simulating a two-phase flow, is the modelling of interfacial forces. They reproduce the different physical effects that are part of the momentum interfacial exchange between phases.

There are many models of interfacial forces. Different models were developed based on the analysis of sets of experimental results that belong to a specific experimental case. For such or similar cases, these models will predict the results well numerically, but when trying to apply these models to a case that differs significantly from that used to optimise the model, unsatisfactory results are frequently obtained. This means that generalised interfacial force models, that will give good results in various cases of two-phase flow, still do not exist. For each force model, it is necessary to know what its area of application is. The modelling of these interfacial forces is still an open question.

In upward turbulent bubbly flows the radial distribution of void fraction depends on the competition between lateral forces – lift force, turbulent dispersion force and wall lubrication force. The wall lubrication force affects the flow in the distance of less than one half bubble diameter from the wall. Outside this area, the void fraction profile is determined by the ratio of the turbulent dispersion force and the lift force. Shawkat et al. ²⁰ shows that it is very difficult to distinguish the conditions under which a wall peak or core peak void fraction profile occurs in bubbly flows and that it depends on numerous parameters such as the ratio of volumetric fluxes of phases, the ratio of pipe diameter and bubble diameter, slip ratio, Eötvös number…

Accurate determination of the bubble drag coefficient is of great importance for the successful simulation of two-phase bubbly flow, since the drag force is the dominant interfacial force in the streamwise direction. The unsteady, three-dimensional flow around a spherical bubble moving steadily in viscous linear shear is analysed numerically by Legendre and Magnaudet ²¹ by solving full Navier-Stokes equations. They found that the drag coefficient of a single bubble in a shear flow increases with the increase in the ratio between the liquid velocity difference across the bubble and the relative velocity of gas and liquid. This increase in the drag coefficient is caused by pressure modification around the bubble.

Hosokawa and Tomiyama ⁵ studied experimentally and numerically (with HRN models) upward turbulent bubbly flow in a vertical pipe. They found that the application of a correlation proposed by Legendre and Magnaudet ²¹ for a drag coefficient of a single bubble in shear flow in numerical simulations reduces the relative velocity of gas and liquid. This led to a better agreement of experimental and numerical results.

The present study examines upward turbulent bubbly flows in the Euler-Euler framework using the LRN model for liquid. To the authors’ knowledge, up to now, there are a limited number of such studies (Colombo and Fairweather, ²² Colombo et al. ⁹ ). Also, there is no paper in the available literature to verify the effect of the drag coefficient of a single bubble in the shear flow in combination with LRN models in upward turbulent bubbly flows in vertical pipes. This paper aims to achieve this.

Description of a numerical method

Governing equations of two-phase flows

Using the index $k$  =  $\mathrm{L}$ , $\mathrm{G}$ to denote the liquid and gas phase, respectively, assuming adiabatic flow, two-fluid model continuity, and momentum equations are given by Ishii and Hibiki ²³ :

$$\frac{\partial }{\partial t}({\alpha }_k{\rho }_k)+\nabla ·({\alpha }_k{\rho }_k{\mathit{U}}_k)=0$$

(1)

and

$$\begin{array}{c}\frac{\partial }{\partial t}({\alpha }_k{\rho }_k{\mathit{U}}_k)+\nabla ·({\alpha }_k{\rho }_k{\mathit{U}}_k\otimes {\mathit{U}}_k)\hfill \\ =-{\alpha }_k\nabla p_k+\nabla ·({\alpha }_k{\mathbf{T}}_k)+{\mathit{F}}_k^{\mathrm{inter}}+{\alpha }_k{\rho }_k\mathit{g},\hfill \end{array}$$

(2)

where ${\alpha }_k$ is void fraction, ${\rho }_k$ the phase density, ${\mathit{U}}_k$ the phase velocity and $p_k$ the phase pressure. ${\mathbf{T}}_k$ is a phase stress tensor that comprises viscous and turbulent stresses, while $\mathit{g}$ denotes the gravitational acceleration. The term ${\mathit{F}}_k^{\mathrm{inter}}$ represents the interfacial force that accounts for the momentum transfer between phases. This force requires modelling.

Modelling of interfacial force

In this paper, the interfacial force ${\mathit{F}}_k^{\mathrm{inter}}$ is represented as the sum of independent forces:

$${\mathit{F}}_k^{\mathrm{inter}}={\mathit{F}}_{\mathrm{D}}+{\mathit{F}}_{\mathrm{L}}+{\mathit{F}}_{\mathrm{TD}}+{\mathit{F}}_{\mathrm{WL}}+{\mathit{F}}_{\mathrm{VM}},$$

(3)

where ${\mathit{F}}_{\mathrm{D}}$ denotes the drag force, ${\mathit{F}}_{\mathrm{L}}$ the lift force, ${\mathit{F}}_{\mathrm{TD}}$ the turbulent dispersion force, ${\mathit{F}}_{\mathrm{WL}}$ the wall lubrication force and ${\mathit{F}}_{\mathrm{VM}}$ the virtual mass force. These forces with closures are briefly presented below.

Drag force

The drag force is given by:

$${\mathit{F}}_{\mathrm{D}}=\frac{3}{4}\frac{C_{\mathrm{D}}}{D_{\mathrm{B}}}{\alpha }_{\mathrm{G}}{\rho }_{\mathrm{L}}|{\mathit{U}}_{\mathrm{R}}|{\mathit{U}}_{\mathrm{R}}.$$

(4)

The drag coefficient $C_{\mathrm{DU}}$ for a single bubble in a uniform stream is defined according to Ishii and Zuber, ²⁴

$$C_{\mathrm{DU}}=\max [C_{\mathrm{D},\mathrm{sphere}},\min (C_{\mathrm{D},\mathrm{ellipse}},C_{\mathrm{D},\mathrm{cap}})],$$

(5)

where:

$$\begin{array}{c}{C}_{\mathrm{D},\mathrm{sphere}}=\frac{24}{\mathrm{R}{\mathrm{e}}_{\mathrm{B}}}(1+0.1\mathrm{R}{\mathrm{e}}_{\mathrm{B}}),\hfill \\ C_{\mathrm{D},\mathrm{ellipse}}=\frac{2}{3}\sqrt{\mathrm{Eo}},\hfill \\ C_{\mathrm{D},\mathrm{cap}}=\frac{8}{3},\hfill \end{array}$$

(6)

where the bubble Reynolds number $\mathrm{R}{\mathrm{e}}_{\mathrm{B}}$ , relative velocity ${\mathit{U}}_{\mathrm{R}}$ between gas and liquid and Eötvös number $\mathrm{Eo}$ are defined as:

$$\mathrm{R}{\mathrm{e}}_{\mathrm{B}}=\frac{{\rho }_{\mathrm{L}}|{\mathit{U}}_{\mathrm{R}}|D_{\mathrm{B}}}{{\mu }_{\mathrm{L}}},$$

(7)

$${\mathit{U}}_{\mathrm{R}}={\mathit{U}}_{\mathrm{G}}-{\mathit{U}}_{\mathrm{L}},$$

(8)

$$\mathrm{Eo}=\frac{g({\rho }_{\mathrm{L}}-{\rho }_{\mathrm{G}})D_{\mathrm{B}}^2}{\sigma }.$$

(9)

To account for the effect of shear rate on the increase in drag coefficient of a single spherical bubble, Legendre and Magnaudet ²¹ proposed the following correlation:

$$C_{\mathrm{D}}(\mathrm{R}{\mathrm{e}}_{\mathrm{B}},\mathrm{Sr})=C_{\mathrm{DU}}(\mathrm{R}{\mathrm{e}}_{\mathrm{B}})(1+0.55\mathrm{S}{\mathrm{r}}^2),$$

(10)

where $\mathrm{Sr}$ is the non-dimensional shear rate given by:

$$\mathrm{Sr}=\frac{D_{\mathrm{B}}\omega }{|{\mathit{U}}_{\mathrm{R}}|},$$

(11)

and $\omega$ is the derivative of the liquid streamwise velocity $U_{\mathrm{L}}$ with respect to wall-normal coordinate $y$ . This quantity compares liquid velocity difference across the bubble to the intensity of a relative velocity between gas and water.

Equation (10) is plotted in Figure 1, together with corresponding experimental data from Hosokawa and Tomiyama. ⁵

Figure 1.

Dependence of the ratio of bubble drag coefficient in shear flow $C_{\mathrm{D}}(\mathrm{R}{\mathrm{e}}_{\mathrm{B}},\mathrm{Sr})$ and bubble drag coefficient in uniform flow $C_{\mathrm{DU}}(\mathrm{R}{\mathrm{e}}_{\mathrm{B}})$ on the non-dimensional shear rate $\mathrm{Sr}$ : () experimental data from Hosokawa and Tomiyama ⁵ ; () equation (10) of Legendre and Magnadeut ²¹ ; () proposed equation (12). The dependence is shown for different ranges: (a) up to the 2 and (b) up to the 20.

The agreement of equation (10) against experimental data is relatively good for non-dimensional shear rate $\mathrm{Sr}$ less than $2$ , although there is some scatter in the experimental data (see Figure 1(a)). For values of $\mathrm{Sr}$ larger than $2$ (see Figure 1(b)) there are no experimental data. Those values of $\mathrm{Sr}$ correspond to the near-wall region, in which it is practically impossible to get reliable measurement data.

With respect to equation (11), the non-dimensional shear rate $\mathrm{Sr}$ is linearly proportional to velocity gradient. Considering that velocity gradient strongly increases toward the wall, and relative velocity ${\mathit{U}}_{\mathrm{R}}$ goes to zero (if no-slip boundary conditions are applied for air and water), large values of $\mathrm{Sr}$ are expected in the near-wall region. According to equation (10), the drag coefficient of a single bubble in shear flow increases with the square of $\mathrm{Sr}$ . Therefore, equation (10) predicts very large value of drag coefficient near the wall, as shown in Figure 1(b). Since the drag force ${\mathit{F}}_{\mathrm{D}}$ defined by the equation (4) is directly proportional to drag coefficient, it follows that ${\mathit{F}}_{\mathrm{D}}$ also tends to very large values in the near-wall region.

If the flow is analysed in the context of LRN models, the centroids of the adjacent wall cells are at very small distance from the wall. In typical simulations in the Euler-Euler framework of upward turbulent bubbly flows in vertical pipes with a LRN model for liquid and bubble drag coefficient defined by the equation (10), it is observed that drag force on a single bubble can achieve values of order ${10}^{12}\mathrm{N}/{\mathrm{m}}^3$ . These large values of drag force, as a rule, produce instabilities and divergence of numerical simulations.

This analysis leads to the conclusion that the drag coefficient of a single bubble in shear flow given by equation (10) of Legendre and Magnaudet ²¹ cannot be used along with the LRN models. On the other hand, equation (10) can be used with HRN models, since the centroids of the adjacent wall cells are at much greater distance from the wall ( $y^{+}>30$ ), which leads to physically acceptable values of drag force. Such numerical study is conducted in Hosokawa and Tomiyama. ⁵

To avoid extremely large values of bubble drag coefficient predicted by equation (10) near the wall in the context of LRN models, we propose the following correlation for the drag coefficient of a single bubble in a shear flow, based on linear regression with experimental data of Hosokawa and Tomiyama ⁵ :

$$C_{\mathrm{D}}(\mathrm{R}{\mathrm{e}}_{\mathrm{B}},\mathrm{Sr})=C_{\mathrm{DU}}(\mathrm{R}{\mathrm{e}}_{\mathrm{B}})(1+0.45\mathrm{Sr}),$$

(12)

where it is assumed $\mathrm{Sr}>0$ . As can be seen from Figure 1, this correlation agrees very well with equation (10) up to $\mathrm{Sr}=2$ but gives significantly lower values of drag coefficients than equation (10) near the wall.

It can be seen in equation (12) that the dependence of the new correlation on $\mathrm{Sr}$ is not quadratic as in equation (10), but linear. Also, the value of the coefficient multiplying $\mathrm{Sr}$ is smaller in equation (12). This allows the new correlation for the drag coefficient to be applied in the near-wall region where $\mathrm{Sr}$ has high values. The proposed correlation for the drag coefficient equation (12), unlike the existing expression given by equation (10), can be used with LRN models.

Lift force

The lift force is given by:

$${\mathit{F}}_{\mathrm{L}}=C_{\mathrm{L}}{\alpha }_{\mathrm{G}}{\rho }_{\mathrm{L}}({\mathit{U}}_{\mathrm{G}}-{\mathit{U}}_{\mathrm{L}})\times (\nabla \times {\mathit{U}}_{\mathrm{L}})$$

(13)

where $C_{\mathrm{L}}$ is a lift coefficient. To calculate this coefficient, in this paper, relation given by Shaver and Podowski ²⁵ is used:

$$C_{\mathrm{L}}=\{\begin{array}{c}0,y/D_{\mathrm{B}}<0.5\hfill \\ \multicolumn{1}{c}{C_{\mathrm{L}0}{f}_{\mathrm{L}}(y)0.5\le y/D_{\mathrm{B}}\le 1}\\ \multicolumn{1}{c}{C_{\mathrm{L}0},1<y/D_{\mathrm{B}}}\end{array}$$

(14)

where $f_{\mathrm{L}}(y)=[3{(\frac{2y}{D_{\mathrm{B}}}-1)}^2-2{(\frac{2y}{D_{\mathrm{B}}}-1)}^3]$ . In numerical simulations, constant lift coefficient $C_{\mathrm{L}0}$ has the value $0.03$ .

Turbulent dispersion force

The turbulent dispersion force is modelled by a correlation given by Burns et al. ²⁶ :

$${\mathit{F}}_{\mathrm{TD}}=\frac{3}{4}\frac{C_{\mathrm{D}}}{D_{\mathrm{B}}}{\alpha }_{\mathrm{G}}|{\mathit{U}}_{\mathrm{G}}-{\mathit{U}}_{\mathrm{L}}|\frac{{\mu }_{\mathrm{tL}}}{{\sigma }_{\mathrm{TD}}}(\frac{1}{{\alpha }_{\mathrm{G}}}+\frac{1}{1-{\alpha }_{\mathrm{G}}})\nabla {\alpha }_{\mathrm{G}}$$

(15)

where ${\sigma }_{\mathrm{TD}}$ represents Schmidt number and ${\mu }_{\mathrm{tL}}$ is turbulent dynamic viscosity coefficient determined by turbulent model. In numerical simulations ${\sigma }_{\mathrm{TD}}$ is equal to $1$ .

Wall lubrication force

To reduce gas volume fraction in near-wall region (distance of one bubble radius from the wall), a wall lubrication force is calculated according to Lubchenko et al., ⁶

$$\mathrm{if}y/D_{\mathrm{B}}<0.5:$$

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{F}}_{\mathrm{WL}}=-\frac{3}{4}\frac{C_{\mathrm{D}}}{D_{\mathrm{B}}}{\alpha }_{\mathrm{G}}|{\mathit{U}}_{\mathrm{G}}-{\mathit{U}}_{\mathrm{L}}|\frac{{\mu }_{\mathrm{t},\mathrm{L}}}{{\sigma }_{\mathrm{TD}}}(\frac{1}{{\alpha }_{\mathrm{G}}}+\frac{1}{1-{\alpha }_{\mathrm{G}}})}\\ \multicolumn{1}{c}{{\alpha }_{\mathrm{G}}\frac{1}{y}\frac{D_{\mathrm{B}}-2y}{D_{\mathrm{B}}-y}\mathbf{n},}\end{array}$$

(16)

$$\mathrm{if}y/D_{\mathrm{B}}>0.5:$$

$${\mathit{F}}_{\mathrm{WL}}=0$$

(17)

where $y$ is a normal-wall distance.

Virtual mass force

The virtual mass force is calculated according to:

$${\mathit{F}}_{\mathrm{VM}}={\alpha }_{\mathrm{G}}{\rho }_{\mathrm{L}}{C}_{\mathrm{VM}}(\frac{\mathrm{D}{\mathit{U}}_{\mathrm{G}}}{\mathrm{D}t}-\frac{\mathrm{D}{\mathit{U}}_{\mathrm{L}}}{\mathrm{D}t})$$

(18)

where virtual mass coefficient $C_{\mathrm{VM}}$ is equal to $0.5$ for a spherical particle.

Turbulence modelling

Turbulence is modelled only in the liquid phase, assuming that the flow of gas phase is laminar.

To model turbulent stresses in the liquid phase, Boussinesq eddy viscosity assumption is applied. It assumes that Reynolds stress tensor ${\mathbf{T}}_{\mathrm{L}}^{\mathrm{turb}}$ is linearly proportional to the mean strain-rate tensor ${\mathbf{S}}_{\mathrm{L}}$ :

$${\mathbf{T}}_{\mathrm{L}}^{\mathrm{turb}}=2{\mu }_{\mathrm{tL}}{\mathbf{S}}_{\mathrm{L}}-\frac{2}{3}\mathbf{I}{\rho }_{\mathrm{L}}{k}_{\mathrm{L}},$$

(19)

where ${\mu }_{\mathrm{tL}}$ is turbulent dynamic viscosity coefficient, $\mathbf{I}$ is unit tensor, and ${\rho }_{\mathrm{L}}$ and $k_{\mathrm{L}}$ are density and turbulent kinetic energy of liquid phase, respectively.

Based on the comparison of single-phase liquid turbulent vertical pipe flows (not shown in this paper) using different LRN models already implemented in OpenFOAM, Launder-Sharma $k$ - $\epsilon$ provided the best agreement with experimental results. Therefore, the low-Reynolds-number model $k$ - $\epsilon$ Launder-Sharma^27,28 is used as a turbulent model in the present study.

In this model, turbulent viscosity is given as:

$${\mu }_{\mathrm{tL}}=C_{\mu }{f}_{\mu }{\rho }_{\mathrm{L}}\frac{k_{\mathrm{L}}^2}{{\epsilon }_{\mathrm{L}}},$$

(20)

where $C_{\mu }=0.09$ , and $f_{\mu }$ represents damping function given as:

$$f_{\mu }=\exp \left[\frac{-3.4}{{(1+\mathrm{R}{\mathrm{e}}_{\mathrm{t}}/50)}^2}\right],\mathrm{R}{\mathrm{e}}_{\mathrm{t}}=\frac{k^2}{{\nu }_{\mathrm{L}}\epsilon }$$

(21)

where $\mathrm{R}{\mathrm{e}}_{\mathrm{t}}$ is turbulent Reynolds number and ${\nu }_{\mathrm{L}}={\mu }_{\mathrm{L}}/{\rho }_{\mathrm{L}}$ is the kinematic viscosity of the liquid.

Turbulent kinetic energy $k_{\mathrm{L}}$ and dissipation rate ${\epsilon }_{\mathrm{L}}$ are obtained by solving equations (22) and (23):

$$\begin{array}{c}\frac{\partial }{\partial t}({\rho }_{\mathrm{L}}{\alpha }_{\mathrm{L}}{k}_{\mathrm{L}})+\nabla ·[{\rho }_{\mathrm{L}}{\alpha }_{\mathrm{L}}{k}_{\mathrm{L}}{\mathit{U}}_{\mathrm{L}}-({\mu }_{\mathrm{L}}+\frac{{\mu }_{\mathrm{tL}}}{{\sigma }_{\mathrm{k}}})\nabla k_{\mathrm{L}}]=\hfill \\ P-{\rho }_{\mathrm{L}}({\epsilon }_{\mathrm{L}}+D_{\epsilon })+S_{\mathrm{k}},\hfill \end{array}$$

(22)

$$\begin{array}{c}\frac{\partial }{\partial t}({\rho }_{\mathrm{L}}{\alpha }_{\mathrm{L}}{\epsilon }_{\mathrm{L}})+\nabla ·[{\rho }_{\mathrm{L}}{\epsilon }_{\mathrm{L}}{\mathit{U}}_{\mathrm{L}}-({\mu }_{\mathrm{L}}+\frac{{\mu }_{\mathrm{tL}}}{{\sigma }_{\mathrm{k}}})\nabla {\epsilon }_{\mathrm{L}}]=\hfill \\ (C_{{\epsilon }_1}{f}_{1}P-C_{{\epsilon }_2}{f}_2{\rho }_{\mathrm{L}}{\epsilon }_{\mathrm{L}})\frac{{\epsilon }_{\mathrm{L}}}{k_{\mathrm{L}}}+{\rho }_{\mathrm{L}}{E}_{\epsilon }+S_{\epsilon },\hfill \end{array}$$

(23)

where $P$ represents the rate of production of turbulent kinetic energy by mean flow gradients and it is defined as:

$$P={\mathbf{T}}_{\mathrm{L}}^{\mathrm{turb}}:\nabla {\mathit{U}}_{\mathrm{L}}.$$

(24)

Model constants have the values:

$C_{{\epsilon }_1}=1.44$ , $C_{{\epsilon }_2}=1.92$ , ${\sigma }_k=1$ and ${\sigma }_{\epsilon }=1.3$ . Damping functions $f_1$ , $f_2$ and extra source terms $D_{\epsilon }$ and $E_{\epsilon }$ are given as:

$$\begin{array}{l}{f}_1=1,f_2=1-0.3\exp \left[-\min \left({\mathrm{Re}}_{\mathrm{t}}^2,50\right)\right],\\ D_{\epsilon }=2{\nu }_{\mathrm{L}}|\nabla (\sqrt{k_{\mathrm{L}}}{|}^2,E_{\epsilon }=2{\nu }_{\mathrm{L}}{\nu }_{\mathrm{tL}}{\left|\nabla (\nabla {\mathit{U}}_{\mathrm{L}})\right|}^2,\end{array}$$

(25)

where ${\nu }_{\mathrm{tL}}={\mu }_{\mathrm{tL}}/{\rho }_{\mathrm{L}}$ is turbulent kinematic viscosity.

It should be noted that in Launder-Sharma model ${\epsilon }_{\mathrm{L}}$ represents modified dissipation rate of turbulent kinetic energy, which is related to originally-defined or physical dissipation via sum ${\epsilon }_{\mathrm{L}}+D_{\epsilon }$ . This sum appears on the right-hand side of the transport equation for turbulent kinetic energy (22). Reason for solving the equation for modified dissipation is purely computational, since there are disadvantages in implementation of the non-zero boundary condition for originally-defined dissipation. The transport equation for modified ${\epsilon }_{\mathrm{L}}$ (23) is solved with the value of zero (or very close to zero) at the wall, and then the correct form of physical ${\epsilon }_{\mathrm{L}}$ is recovered through the sum ${\epsilon }_{\mathrm{L}}+D_{\epsilon }$ , which is then used in $k_{\mathrm{L}}$ transport equation. The value of $D_{\epsilon }$ is significant near the wall, while it is negligible away from the wall.

Term $S_{\mathrm{k}}$ in equation (22) represents the interfacial transfer of turbulent kinetic energy, while the term $S_{\epsilon }$ in equation (23) represents the interfacial transfer of dissipation rate. These terms are known as bubble induced turbulence, and need to be modelled.

Bubble induced turbulence (BIT) closure modelling

Due to the presence of the bubbles and their interaction with the liquid phase, the profiles of liquid turbulent kinetic energy and liquid turbulent dissipation rate are modified compared to the case of pure liquid flow. To take into account this modification, source terms $S_{\mathrm{k}}$ and $S_{\epsilon }$ are introduced in equations (22) and (23), respectively. These source terms are expressed as:

$$S_{\mathrm{k}}=K_{\mathrm{BI}}{\mathit{F}}_{\mathrm{D}}·{\mathit{U}}_{\mathrm{R}},$$

(26)

$$S_{\epsilon }=C_{\epsilon \mathrm{B}}\frac{S_{\mathrm{k}}}{\tau },$$

(27)

where ${\mathit{F}}_{\mathrm{D}}$ represents the drag force on bubbles, ${\mathit{U}}_{\mathrm{R}}$ represents the relative velocity between phases, constant $C_{\epsilon \mathrm{B}}=1$ and $\tau$ denotes the characteristic time scale of bubble-induced turbulence. $\tau$ is calculated as a ratio of the bubble diameter and the square root of turbulent kinetic energy, $\tau =D_{\mathrm{B}}/\sqrt{k_{\mathrm{L}}}$ (Rzehak and Krepper, ¹³ Colombo and Fairweather ²⁹ ).

Parameter $K_{\mathrm{BI}}$ modulates turbulent kinetic energy generation. Colombo and Fairweather ²⁹ proposed $K_{\mathrm{BI}}=0.25$ and this value is used in the present paper.

Summary of experimental data

To examine described numerical models for upward turbulent bubbly flow in a vertical pipe, comparisons of numerical simulation results against experimental databases for such flows are necessary. The selected databases are: Hosokawa and Tomiyama, ⁵ Liu ³⁰ and Shawcat et al. ²⁰ Main parameters of these experiments are listed in Table 1. The experimental test facilities used in these three cases are similar but still differ from each other. Therefore, not all of them are presented separately in this paper. They are graphically presented and described in detail in the mentioned papers. Also, more specific details about the measurement procedure and the measuring equipment used can be found there. A brief description of each experiment follows.

Table 1.

Summary of experimental sets: pipe diameter $D$ , pipe length $L$ , distance from the pipe inlet to measuring section $L_{\mathrm{m}}$ , gas and liquid volumetric fluxes, $J_{\mathrm{G}}$ and $J_{\mathrm{L}}$ , respectively, mean void fraction $〈{\alpha }_{\mathrm{G}}〉$ and average bubble diameter $〈D_{\mathrm{B}}〉$ .

Case	$D[\mathrm{mm}]$	$L[\mathrm{mm}]$	$L_{\mathrm{m}}[\mathrm{mm}]$	$J_{\mathrm{G}}[\mathrm{m}/\mathrm{s}]$	$J_{\mathrm{L}}[\mathrm{m}/\mathrm{s}]$	$〈{\alpha }_{\mathrm{G}}〉[\%]$	$〈D_{\mathrm{B}}〉[\mathrm{mm}]$
H1	25	2000	1700	0.018	0.500	2.31	3.21
H2	25	2000	1700	0.025	0.500	3.99	4.25
H3	25	2000	1700	0.020	1.000	1.46	3.52
H4	25	2000	1700	0.036	1.000	3.30	3.66
L1	38	2800	1368	0.027	0.376	4.11	2.68
L15	38	2800	1368	0.027	0.753	2.63	2.43
L29	38	2800	1368	0.027	1.087	1.91	2.07
S1	200	9560	8400	0.015	0.450	2.30	3.85
S2	200	9560	8400	0.030	0.450	3.55	4.03

‘H*’, ‘L*’ and ‘S*’ denote experiments from Hosokawa and Tomiyama, ⁵ Liu ³⁰ and Shawcat et al., ²⁰ respectively, and * represents the number of experiment.

Hosokawa and Tomiyama experimental database

In experiment of Hosokawa and Tomiyama, ⁵ the inner diameter of the vertical pipe was $D=25\mathrm{mm}$ , and the pipe length was $L=2\mathrm{m}$ . Measurements were performed at a section placed $L_{\mathrm{m}}=1.7\mathrm{m}$ ( $L_{\mathrm{m}}/D=68$ ) above the pipe inlet.

The mixing section of air and water was at the bottom of the vertical pipe. Gas and liquid volumetric flows were measured by flowmeters before the mixing section. The liquid and gas volumetric fluxes $J_{\mathrm{L}}$ and $J_{\mathrm{G}}$ were calculated as the ratio of the measured flow rates and cross-sectional area of the pipe.

A laser Doppler velocimetry system (LDV, DANTEC, optics: 60 × 83, processor: 68N10) was used for measuring liquid velocities. Around 50,000 velocity data were taken for each measurement point. The uncertainty estimated at $95\mathrm{\%}$ confidence for $J_{\mathrm{L}}$ was $1\mathrm{\%}$ . Simultaneously, two high-speed cameras (Kodak Motion Coder Analyzer SR-1000, shutter speed: $1/5000\mathrm{s}$ , frame rate: 250 fps) were used to obtain stereoscopic images of bubbles. The spatial resolution of the image was $0.07\mathrm{mm}/\mathrm{pixel}$ . From these images, the three-dimensional bubble shapes were reconstructed. This reconstruction allows the determination of several quantities: bubble volumes $V_{\mathrm{B}}$ , sphere-volume equivalent diameter $D_{\mathrm{B}}$ , aspect ratio $E$ , the bubble size distribution and void fraction ${\alpha }_{\mathrm{G}}$ . The bubble velocity was obtained by calculating the displacement of the bubble centre on two consecutive images. The time between taking two consecutive images is known. For each flow condition, about 500 bubbles were processed. The uncertainties estimated at $95\mathrm{\%}$ confidence for measured ${\alpha }_{\mathrm{G}}$ , $U_{\mathrm{G}}$ , $D_{\mathrm{B}}$ and $E$ were $10\mathrm{\%}$ , $3\mathrm{\%}$ , $4\mathrm{\%}$ and $10\mathrm{\%}$ respectively.

Experiments were performed for four cases, in the present study referred as H1, H2, H3 and H4. This experimental database provides results of measurements of the relative velocity between gas and water ${\mathit{U}}_{\mathrm{R}}={\mathit{U}}_{\mathrm{G}}-{\mathit{U}}_{\mathrm{L}}$ , turbulent kinetic energy $k_{\mathrm{L}}$ of liquid phase and void fraction ${\alpha }_{\mathrm{G}}$ . Individual values of phase velocities are not available in this database.

Liu experimental database

The results of an extensive set of experimental investigations of an upward two-phase bubble flow through the vertical pipe are presented in Liu. ³⁰

The pipe was $L=2800\mathrm{mm}$ long, with the inner diameter of $D=38\mathrm{mm}$ . The measuring section was at $L_{\mathrm{m}}=1368\mathrm{mm}$ , $(L_{\mathrm{m}}/D=36)$ above the pipe inlet.

For injecting air into the water flow, a bundle of $64$ equally-spaced hypodermic needles was used. The diameter of the needle was $0.1\mathrm{mm}$ . This method of mixing was chosen to generate small bubbles and minimise bubble coalescence.

At the measuring station, both one-dimensional hot film sensor (TSI 1218-20W, $0.05\mathrm{mm}$ o.d. ×  $1\mathrm{mm}$ ) and two-dimensional hot film anemometer probe (TSI 1246-60W, X-type, $0.15\mathrm{mm}$ o.d. ×  $2\mathrm{mm}$ long) were used to measure liquid velocity $U_{\mathrm{L}}$ , and Reynolds stress ( $\overline{\mathrm{uu}}$ , $\overline{\mathrm{vv}}$ , $\overline{\mathrm{uv}}$ ). The uncertainty for every flow condition is within $\pm 5\mathrm{\%}$ . Most of the data for liquid-phase parameters were checked more than twice. The consistency of liquid velocity was within $\pm 1.5\mathrm{\%}$ , and the turbulent stresses were within $\pm 3\mathrm{\%}$ . By analysing a hot-film anemometer response, the presence of a liquid-bubble interface was recognised.

A dual-sensor resistivity probe was used for measuring the bubble velocity $U_{\mathrm{G}}$ , bubble size $D_{\mathrm{B}}$ and radial profiles of the void fraction ${\alpha }_{\mathrm{G}}$ . The uncertainty for bubble velocity measurements was within $\pm 5\mathrm{\%}$ . The reproducibilities of the local bubble velocity and the void fraction measured by the resistivity probe were within $\pm 2\mathrm{\%}$ .

Because of using hot-film anemometry, it was important to precisely control the water temperature. Water was cooled in the storage tank to achieve water temperature $10\pm 0.01°\mathrm{C}$ at the measuring station.

From $42$ sets of experimental setups, three cases L1, L15 and L29 are chosen for analysis in the present paper. Experimental results of the radial distribution of liquid and gas velocities, $U_{\mathrm{L}}$ and $U_{\mathrm{G}}$ , respectively, kinetic turbulence energy of liquid phase $k_{\mathrm{L}}$ and void fraction ${\alpha }_{\mathrm{G}}$ are available in this database.

Shawkat et al. experimental database

The ratio of pipe diameter and mean bubble diameter has a significant influence on the flow patterns and phase distribution. ²⁰ Therefore, in contrast to Hosokawa and Tomiyama ⁵ and Liu, ³⁰ who conducted two-phase bubbly flow experiments in vertical pipes of diameter $D=25\mathrm{mm}$ and $D=38\mathrm{mm}$ , respectively, Shawkat et al. ²⁰ analysed turbulent air-water bubbly flow in a vertical pipe of a larger diameter $D=200\mathrm{mm}$ .

The pipe length was $L=9560\mathrm{mm}$ . Measurements were performed in a cross section at $L_{\mathrm{m}}=8400\mathrm{mm}$ $(L_{\mathrm{m}}/D=42)$ above the pipe inlet.

For mixing of air and water, a showerhead injector was used. The injector had $550$ holes of $1\mathrm{mm}$ diameter. After the injector, the air-water mixture went through a honeycomb flow straightener. The goal was to reduce bubble swirls and improve bubble distribution.

During the experiments, the water temperature was kept within a range of $24.5\pm 0.1°\mathrm{C}$ . For the liquid turbulence characteristics measurements, hot film anemometry was used. By using a single hot film probe (TSI 1210-60W) and X-hot film probe (Dantec 55R63) liquid velocity and turbulent stresses were measured. A dual optical probe was used for the measurement of bubble characteristics. The vertical and horizontal distances between the probe tips were $1.16\mathrm{mm}$ and $1\mathrm{mm}$ . Relative uncertainty for the liquid velocity $U_{\mathrm{L}}$ is within $\pm (6-9)\mathrm{\%}$ , for turbulent stresses up to $13\mathrm{\%}$ and for void fraction ${\alpha }_{\mathrm{G}}$ $14\mathrm{\%}$ .

In this experimental database radial distributions of liquid velocity $U_{\mathrm{L}}$ , gas velocity $U_{\mathrm{G}}$ and void fraction ${\alpha }_{\mathrm{G}}$ are given.

Numerical simulation set-up

To perform numerical simulations of two-phase bubbly flow in a vertical pipe, solver twoPhaseEulerFoam of OpenFOAM v1906 is modified. These modifications include implementations of the following:

• proposed relation for increase in the drag coefficient due to shear, given by equation (12)

• lift force, equation (14)

• wall lubrication force, equations (16) and (17)

• BIT source terms, equations (26) and (27).

The computational domain, initial and boundary conditions and physical parameters of phases are explained in the continuation of this section.

Computational domain

Since in all analysed experiments the flow was assumed as axisymmetric, in order to reduce the number of computational cells, a quarter-pipe geometry is used in all simulations. The sketch of the computational domain is shown in Figure 2.

Figure 2.

Sketch of a computational domain. A quarter-pipe is used in simulations.

The computational domain is bounded by five surfaces: inlet and outlet in the horizontal plane, pipe wall and two vertical symmetry planes.

Values of pipe diameters $D$ , pipe lengths $L$ and distances from the inlet to measurement section $L_{\mathrm{m}}$ for experimental databases of Hosokawa and Tomiyama, ⁵ Liu ³⁰ and Shawkat et al. ²⁰ are presented in Table 1.

The computational domain is discretised using O-grid block topology. This provides a high level of control of important mesh quality parameters like orthogonality, skewness and stretching of the cells. A typical example of this mesh in a pipe cross-section is shown in Figure 3. Mesh is generated using the blockMesh utility of OpenFOAM.

Figure 3.

Computational mesh of a pipe cross section.

Since the LRN model is used, the value of non-dimensional distance $y^{+}$ of the adjacent wall cells is set to be less than 1. Edges of the cells are graded towards the wall, while in the flow direction cells have uniform edge lengths.

Mesh refinement study

Suitable discretisation of computational domain concerning the optimal number of cells is determined by grid sensitivity studies for each experimental database. Numerical calculations are performed on meshes with a different number of cells. Results are compared and optimal mesh is chosen.

As an example, a sensitivity study for one of the experiments of Hosokawa and Tomiyama ⁵ (case H1, see Table 1) is shown in Figure 4. Numerical calculations are performed on three meshes, named Mesh 1 (low refinement), Mesh 2 (medium refinement) and Mesh 3 (high refinement). The number of cells is presented in Table 2.

Figure 4.

The mesh sensitivity test for case H1 (see Table 1) is performed on three meshes: Mesh 1 with 165,800 cells (), Mesh 2 with 320,000 cells (), and Mesh 3 with 409,400 cells (······). The results are shown for: (a) relative velocity, (b) void fraction and (c) liquid turbulent kinetic energy.

Table 2.

Cell numbers of meshes used in mesh sensitivity tests for cases H1, L1 and S1 (see Table 1).

Case	Low resolution	Medium resolution	High resolution
H1	165,800	320,000	409,400
L1	210,200	425,000	563,800
S1	368,000	760,400	908,400

It is noticeable that the radial distribution of relative velocity and void fraction does not depend on the mesh density (Figure 4(a) and (b)). However, the kinetic energy profile (Figure 4(c)) obtained using Mesh 1 differs from the profile obtained using Mesh 2. No additional differences are found from the medium (Mesh 2) to the high refinement (Mesh 3) resolution.

It is concluded that the results obtained using Mesh 2 do not depend on the mesh density. Similar analyses were performed for each case listed in Table 1.

Table 2 provides data on meshes corresponding to cases H1, L1 and S1, of Hosokawa and Tomiyama, ⁵ Liu ³⁰ and Shawkat et al., ²⁰ respectively (see Table 1). Meshes with medium refinement are found as optimal computational meshes.

Initial and boundary conditions

Uniform phase velocities are set as the boundary condition at the inlet surface. These velocities are calculated as the ratio of phase volumetric fluxes and the corresponding phase volume fraction, $U_j=J_j/{\alpha }_j$ ( $j$  =  $\mathrm{L}$ , $\mathrm{G}$ ). Values of liquid turbulent kinetic energy $k_{\mathrm{L}}$ and turbulent dissipation rate ${\epsilon }_{\mathrm{L}}$ in inlet section are estimated by the following expressions: $k_{\mathrm{L}}=3/2(I|U_{\mathrm{L}}|{)}^2$ and ${\epsilon }_{\mathrm{L}}=(C_{\mu }^{0.75}{k}_{\mathrm{L}}^{1.5})/D$ , where $I=0.1$ is turbulence intensity and $C_{\mu }=0.09$ is the model constant. The uniform value of the void fraction is defined at the inlet.

At the outlet surface, a zero-gradient boundary condition is defined for $k_{\mathrm{L}}$ , ${\epsilon }_{\mathrm{L}}$ , void fraction ${\alpha }_{\mathrm{G}}$ and both phase velocities $U_{\mathrm{L}}$ and $U_{\mathrm{G}}$ .

At the pipe wall, a no-slip boundary condition is applied for both phase velocities and fixed zero value for the void fraction. For turbulent kinetic energy and dissipation zero values are prescribed at wall (actually, values of ${10}^{-10}{\mathrm{m}}^2/{\mathrm{s}}^3$ in order to avoid dividing with zero).

Initial conditions inside the computational domain for all quantities are identical to their values in the inlet section.

Physical parameters of phases

Liquid flow is calculated as incompressible. Gas is treated as an ideal gas. Phase densities, phase viscosities and surface tension are prescribed based on their values at atmospheric pressure and temperatures from experiments.

Mean void fraction and mean bubble diameter are defined in papers of Hosokawa and Tomiyama ⁵ and Liu, ³⁰ while for experimental case of Shawkat et al. ²⁰ these quantities are calculated by integration of experimental results:

$$\begin{array}{c}〈{\alpha }_{\mathrm{G}}〉=\frac{2}{R^2}{\int }_0^R{\alpha }_{\mathrm{G}}(r)r\mathrm{d}r,\hfill \\ 〈D_{\mathrm{B}}〉=\frac{2}{R^2}{\int }_0^R{D}_{\mathrm{B}}(r)r\mathrm{d}r.\hfill \end{array}$$

(28)

Numerical simulations are conducted with constant mean bubble diameter. Spherical bubbles are assumed.

Volumetric fluxes of gas and liquid, mean bubble diameters and void fractions for experimental cases examined are specified in Table 1.

Solution method

The pseudo-steady state solution is obtained using PIMPLE algorithm for pressure-velocity coupling. Time step size is set up to be adaptive, based on the condition of maximum Courant number equal to $0.25$ for each time step. Simulations are run for multiple flow-through times to obtain converged solutions. For cases of Hosokawa and Tomiyama ⁵ and Liu, ³⁰ it is necessary to run simulations for $6$  s of physical time, and for cases of Shawkat et al. ²⁰ (due to lower flow through time) it is necessary to run simulations for $38$  s of physical time.

Spatial and time discretisation schemes are of the second order.

Results and discussion

In order to determine the effects of the new proposed correlation for the drag coefficient of a single bubble in a shear flow, given by equation (12), the results of numerical simulations conducted with this new correlation are compared to the results of numerical simulations conducted with the drag coefficient for a single bubble in a uniform stream, given by equation (5). The remaining conditions for these simulations are the same, as described in previous sections. Results of both numerical simulations are tested against the experimental data of Hosokawa and Tomiyama, ⁵ Liu ³⁰ and Shawcat et al. ²⁰

Hosokawa and Tomiyama

The first column of Figure 5 shows radial profiles of relative velocity calculated in numerical simulations along with experimental data of Hosokawa and Tomiyama. ⁵

Figure 5.

Radial distribution of relative velocity $U_{\mathrm{R}}$ , void fraction ${\alpha }_{\mathrm{G}}$ and liquid turbulent kinetic energy $k_{\mathrm{L}}$ : () experimental data from Hosokawa and Tomiyama ⁵ ; () numerical predictions when drag coefficient is calculated for a single bubble in a uniform stream, equation (5); () numerical predictions with included new correlation for the drag coefficient, equation (12). The first three diagrams (a), (b) and (c) represents case H1, the diagrams (d), (e) and (f) represents case H2, the diagrams (g), (h) and (i) case H3 and last three diagrams (j), (k) and (l) represents case H4.

Numerical simulations with a drag coefficient for a single bubble in a uniform stream, given by equation (5), for all cases examined (H1, H2, H3 and H4) predict uniform relative velocity up to $r/R\approx 0.98$ , followed by small velocity increase near the wall and then a steep decrease to a zero relative velocity at the wall. For cases H1 and H2 (Figure 5(a) and (d)), the agreement of the relative velocity obtained from numerical simulations obtained using equation (5) and experimental results is moderately good for $r/R<0.8$ , while the relative velocity for $r/R>0.8$ is overpredicted. The discrepancy between experimental and numerical results obtained using equation (5) increases with the increase of liquid volumetric flux, as shown for cases H3 and H4 (Figure 5(g) and (j)). For these two cases, this discrepancy is pronounced for $r/R>0.6$ , where is the high influence of shear.

In contrast, if the drag coefficient is calculated for a single bubble in a shear flow, given by equation (12), numerical simulation results for relative velocity compare very well against corresponding experimental results, for all cases examined (Figure 5(a), (d), (g) and (j)). With the increase of liquid volumetric flux (for the same pipe diameter), there is an increase in the liquid shear rate. That further increases the drag coefficient of a single bubble (see equation (12)) and reduces the relative velocity, assuming that the drag force on a bubble is constant. This explanation is in accordance with experimental data: a stronger relative velocity reduction is observed for cases H3 and H4 (Figure 5(g) and (j), respectively) with higher liquid volumetric fluxes compared to cases H1 and H2 (Figure 5(a) and (d), respectively), with lower liquid volumetric fluxes.

This shows that in this particular case, under all the above conditions of two-phase bubbly flow, it is necessary to take into account the influence of shear flow on the drag coefficient (described by the new equation (12)) in order to correctly predict the relative velocity profile.

Radial distributions of void fractions obtained from numerical simulations and experimental results of Hosokawa and Tomiyama ⁵ are shown in the second column of Figure 5. According to experimental data, all cases examined have pronounced wall-peak void fraction profiles, except case H2 (Figure 5(e)) that has an approximately uniform void fraction up to $r/R<0.8$ and from that point towards the wall void fraction decreases. Numerical simulations predict wall-peaked void fraction profiles for all cases examined (cases H1, H2, H3 and H4). According to these cases, it can be concluded that the increase in a drag coefficient of a single bubble due to shear (compared to the drag coefficient of a bubble in the uniform stream) leads to a very small increase in the void fraction in the core region of a pipe, while this increase in the drag coefficient does not influence the radial distribution of void fraction for $r/R>0.8$ . For $r/R<0.7$ , in case H1 (Figure 5(b)) numerical simulation results of void fraction distribution compare very well against experimental results, whereas for case H2 (Figure 5(e)) numerical simulation underpredicts void fraction. For $r/R>0.7$ , for cases H1 and H2 there is a small overprediction of the void fractions computed in numerical simulations, compared to experimental data. Regarding void fraction distributions, deviations of numerical simulation results for higher liquid volumetric fluxes (cases H3 and H4, Figure 5(h) and (k), respectively) are more pronounced compared to lower liquid volumetric fluxes (cases H1 and H2). For cases H3 and H4, values of void fraction distribution obtained from numerical simulation are higher than experimental results for $r/R<0.6$ in case H3 and $r/R<0.8$ in case H4. In contrast, for $r/R>0.6$ in case H3 and $r/R>0.8$ in case H4, numerical simulations underpredict values of void fraction. The reason for this discrepancy could be in the interfacial force models or BIT models that demand more investigations in future studies. Also, in numerical simulations monodisperse approximation is used, while in real experiment, there is a certain distribution of bubble sizes.

The third column in Figure 5 represents the radial distribution of liquid turbulent kinetic energy obtained from numerical simulations and experimental data of Hosokawa and Tomiyama. ⁵ The overall agreement between experimental and numerical simulation results is good. For smaller liquid volumetric flux, cases H1 and H2 (Figure 5(c) and (f), respectively) there is a small underprediction of numerical simulation results for $r/R<0.8$ . However, for higher liquid volumetric flux, cases H3 and H4 (Figure 5(i) and (l), respectively), there is an overprediction of numerical simulation results for $r/R<0.8$ . For all cases examined, there is an underprediction of numerical results of liquid kinetic energy peak near the wall. This is in agreement with the study of Parekh and Rzehak. ⁸ They suggest that this underprediction could be a shortcoming of RANS models in general.

Considering that results of numerical simulations conducted with drag coefficient for a single bubble in uniform stream, equation (5), almost coincide with results of numerical simulations conducted with drag coefficient for a single bubble in shear flow equation (12), it can be concluded that an increase in the drag coefficient due to shear flow does not affect liquid turbulent kinetic energy.

According to the case of Hosokawa and Tomiyama, ⁵ it can be concluded that the biggest influence of a drag coefficient of a single bubble in a shear flow (compared to a drag coefficient of a single bubble in a uniform stream) is achieved on the relative velocity of air and water. This conclusion is checked when the flow occurs in pipes of larger diameter in experiments of Liu ³⁰ and Shawkat et al. ²⁰

Liu

In contrast to the experimental database of Hosokawa and Tomiyama ⁵ where relative velocities of phases are given, in the experimental database of Liu ³⁰ both air and water velocity profiles for all experiments are available. Therefore, the influence of the drag coefficient of bubbles, given by equations (5) and (12), on both the liquid and gas phase, can be analysed separately.

Liquid velocity profiles from numerical simulations and experimental database of Liu ³⁰ are plotted in the first row of Figure 6. It can be seen that there is an excellent agreement of numerical simulation results conducted with a drag coefficient for a single bubble in a uniform flow, given by equation (5) and results of numerical simulations conducted with a drag coefficient for a single bubble in the shear flow, given by equation (12). The reason for this behaviour is that due to small values of mean void fractions $〈{\alpha }_{\mathrm{G}}〉$ (in the range from $1.91$ to $4.11\mathrm{\%}$ in cases examined), the magnitudes of inertial forces of water are not affected if a drag coefficient of bubbles is calculated for a single bubble in the uniform stream or a single bubble in the shear flow. Therefore, there is no modification of the liquid velocity profile.

Figure 6.

Radial distribution of liquid phase velocity $U_{\mathrm{L}}$ , gas phase velocity $U_{\mathrm{G}}$ , void fraction ${\alpha }_{\mathrm{G}}$ and liquid turbulent kinetic energy $k_{\mathrm{L}}$ : () experimental data from Liu ³⁰ ; () numerical predictions when drag coefficient is calculated for a single bubble in a uniform stream, equation (5); () numerical predictions with included new correlation for the drag coefficient, equation (12). The diagrams in the first column (a), (d), (g) and (j) represents case L1, the diagrams (b), (e), (h) and (k) represents case L15, and the diagrams (c), (f), (i) and (l) represents case L29.

For case L1 shown in Figure 6(a), liquid velocity profile from experimental data agree very well with liquid velocity profile calculated in a numerical simulation. However, for case L15, shown in Figure 6(b), there is a small overprediction of numerical simulation results for $r/R<0.6$ and small underprediction of these results for $r/R>0.6$ . With the increase of liquid volumetric flux, these discrepancies are more pronounced, which is confirmed in case L29 shown in Figure 6(c).

Regarding gas velocity profiles, for cases L1, L15 and L29, shown in Figure 6(d) to (f), respectively, results of numerical simulations conducted with a drag coefficient for a single bubble in the uniform stream, equation (5), agree very well with results of numerical simulations conducted with a drag coefficient for a single bubble in the shear flow, equation (12). For all cases examined, only a small difference between these results is visible for $r/R>0.8$ .

For case L1 (Figure 6(d)), the gas velocity profile obtained from numerical simulations is slightly overpredicted for $r/R<0.8$ compared to the corresponding experimental results, while the agreement of numerical simulation results and experimental measurements for $r/R>0.8$ is very good. For case L15 (Figure 6(e)), the gas velocity calculated in numerical simulations is slightly larger than the gas velocity from experimental data for $r/R<0.5$ . In contrast, for $r/R>0.5$ gas velocity calculated in numerical simulation is smaller than the gas velocity measured in the experiment. This velocity underestimation increases towards the pipe wall. For case L29 (Figure 6(f)), there is an agreement between experimental and numerical simulation results around pipe axis (for $r/R<0.2$ ). However, for $r/R>0.2$ , experimental results of gas velocity are larger than the gas velocity calculated in numerical simulations. Similar to the case L15 (Figure 6(e)), this discrepancy increases towards the pipe wall. It can be concluded that increase in liquid volumetric flux increases discrepancy between gas velocity profiles calculated in numerical simulations and the corresponding experimental results.

Unlike the water velocity profile, the air velocity profile shows the dependence on the choice of the drag coefficient modelling. Taking into account the influence of shear flow on the drag coefficient has an impact on the air velocity profile, but in this case, it is small.

For cases L1, L15 and L29, it can be seen that profiles of void fraction (shown in Figure 6(g)–(i)), as well as profiles of liquid turbulent kinetic energy (shown in Figure 6(j)–(l)), obtained from numerical simulations conducted with a drag coefficient for a single bubble in a uniform stream agree very well with corresponding numerical simulation results performed with a drag coefficient for a single bubble in a shear flow.

For case L1, (Figure 6(g)), the agreement of radial void fraction calculated in the numerical simulation and experimental results is moderately good. However, with the increase of liquid volumetric flux, cases L15 and L29 (shown in Figure 6(h) and (i), respectively), numerical simulation results of void fraction become larger than the corresponding experimental values for $r/R<0.8$ . At the same time, the void fraction peak near the pipe wall becomes underpredicted by numerical simulation (for both cases L15 and L29).

Experimental data of the radial distribution of liquid turbulent kinetic energy for case L1, shown in Figure 6(j), are slightly larger and more uniform across the cross-section than the results of liquid turbulent kinetic energy obtained from numerical simulations. For case L15, shown in Figure 6(k), the agreement of experimental results of liquid turbulent kinetic energy and corresponding results of numerical simulations is moderately good. However, further increase in liquid volumetric flux increases discrepancies between experimental and numerical simulation results, which is seen in Figure 6(l). As discussed in the results of Hosokawa and Tomiyama, ⁵ the reason for this discrepancy could be in the interfacial force models and/or in BIT models. Also, in numerical simulations monodisperse approximation has been used.

Shawcat et al

The first column of Figure 7 shows the liquid velocity profiles from the experiment of Shawcat et al. ²⁰ and corresponding liquid velocity profiles calculated in numerical simulations. In these figures, there is no difference in liquid velocity profiles if numerical simulations are conducted with the drag coefficient for a single bubble in a uniform stream (equation (5)) or with a drag coefficient for a single bubble in a shear flow (equation (12). The reason for this is that due to small values of mean void fractions $〈{\alpha }_{\mathrm{G}}〉$ ( $2.30$ and $3.55\mathrm{\%}$ in cases examined), large magnitudes of inertial forces of water are not affected by small changes of drag forces of bubbles. This behaviour is also described in the analysis of experimental and numerical simulation results of cases from Liu. ³⁰

Figure 7.

Radial distribution of liquid phase velocity $U_{\mathrm{L}}$ , gas phase velocity $U_{\mathrm{G}}$ and void fraction ${\alpha }_{\mathrm{G}}$ : () experimental data from Shawkat et al. ²⁰ ; () numerical predictions when drag coefficient is calculated for a single bubble in a uniform stream, equation (5); () numerical predictions with included new correlation for the drag coefficient, equation (12). The diagrams (a), (b) and (c) represents the results for case S1, while the diagrams (d), (e) and (f) represents the results for case S2.

Experimental results show that liquid velocity distribution is approximately uniform for $r/R<0.6$ in Figure 7(a) and for $r/R<0.4$ in Figure 7(d)). Towards the pipe wall, liquid velocity decreases slightly for both cases.

The second column of Figure 7 shows the radial distribution of air velocity. It is seen that air velocity profiles obtained in numerical calculations with a drag coefficient for a single bubble in the uniform flow, equation (5), and with a drag coefficient for a single bubble in the shear flow equation (12) are almost identical. The pipe diameter used in Shawkat et al. ²⁰ is $D=200\mathrm{mm}$ , which is larger than the pipe diameter $D=38\mathrm{mm}$ used in experiments of Liu. ²⁰ The increase in the pipe diameter, at similar liquid bulk velocity values, decreases the liquid velocity gradient. Considering that a drag coefficient of a single bubble in a shear flow depends on a liquid velocity gradient, which is seen from equations (11) and (12), it is concluded that the influence of the new correlation for the drag coefficient of a single bubble in a shear flow on air velocity profiles decreases as a liquid velocity gradient decreases. This is why the influence of correlation (12) on air velocity profiles is visible in cases of Liu ³⁰ (Figure 6(d)–(f)), while it is not visible in the cases of Shawkat et al. (Figure 7(b) and (e)).

Results of numerical computations overpredict liquid velocity in almost the entire cross-section, shown in Figure 7(a) and (d). In contrast, numerical computations predict slightly lower air velocities relative to the experimentally measured values, shown in Figure 7(b) and (e).

The third column of Figure 7 shows the radial distribution of the void fraction from experimental results and void fraction calculated in numerical simulations. Experimental results for the void fraction show an approximately uniform radial distribution with a small peak at $r/R=0.85$ . Towards the wall, void fraction slowly decreases. However, numerical calculations overpredict void fraction in the entire domain. This is more pronounced in the case S2 (Figure 7(f)), for higher gas volumetric flux. Numerical results show that the peak is at a radial distance $r/R=0.97$ and after that, void fraction decreases sharply. No values of turbulence kinetic energy are available in this database.

Besides imperfections in interfacial forces and/or in BIT models, as well as due to monodisperse approximation used in numerical simulations, an additional source for disagreement between experimental and numerical results in Figure 7 could be the difference between liquid and gas volumetric fluxes obtained from flowmeters and by integration of experimental measurements obtained using optical and hot-film probes in the measurement section. These differences for cases of Shawkat et al. ²⁰ go up to around $10\mathrm{\%}$ , while for the cases of Liu ³⁰ they go up to $1\mathrm{\%}$ . These errors cannot be calculated for cases of Hosokawa and Tomiyama, ⁵ since only relative velocity is available in the experimental database. According to Figure 7(b) and (e), probes overestimate air velocity along the complete cross-section. In numerical simulations phase velocities calculated from flowmeters are specified as velocity boundary conditions at the inlet. These velocities are shown in Table 1.

Conclusions

Upward turbulent bubbly flows in vertical pipes are studied in the Euler-Euler framework of OpenFOAM. In almost all previous numerical studies of these flows available in the literature, high-Reynolds-number (HRN) models were used. In this paper, the turbulent bubbly flow is analysed using the low-Reynolds-number (LRN) $k$ - $\epsilon$ model.

To account for the increase in the drag coefficient of a single bubble due to shear (compared to drag coefficient of a single bubble in uniform flow), a relation of Legendre and Magnaudet ²¹ is examined. It is found that the drag coefficient defined by Legendre and Magnaudet, ²¹ equation (10), goes to infinity towards the pipe wall, with an increase especially pronounced in cells near the pipe wall. That causes numerical instabilities and therefore this correlation is not applicable with LRN models.

To consider the increase in the drag coefficient of a single bubble due to shear, but to provide stability of numerical simulations performed with LRN models, a new correlation for the drag coefficient (equation (12)) is proposed. The new correlation is based on linear regression with experimental data of Hosokawa and Tomiyama. ⁵

To examine the effects of the proposed correlation for the drag coefficient for a single bubble in a shear flow (equation (12)), the results of numerical simulations conducted with it are compared to the results of numerical simulations conducted with the drag coefficient for a single bubble in a uniform stream (equation (5)). The remaining conditions for these two types of simulations are the same. Besides the drag force, other interfacial forces taken into account by these numerical simulations are: lift force of Shaver and Podowski, ²⁵ turbulent dispersion force of Burns et al., ²⁶ wall lubrication force of Lubchenko et al. ⁶ and virtual mass force. The bubble induced turbulence model of Colombo and Fairweather ²⁹ is used in numerical simulations.

The numerical simulation results are compared against experimental databases of Hosokawa and Tomiyama, ⁵ Liu ³⁰ and Shawkat et al. ²⁰

A comparison of numerical simulation results with all three experimental databases shows that the choice of correlation for calculating the drag coefficient, whether it takes into account the influence of shear flow (equation (12)) or not (equation (5)), has no significant effect on the profiles of liquid turbulent kinetic energy $k$ and void fraction ${\alpha }_{\mathrm{G}}$ .

The biggest influence of the increase of the drag coefficient of a single bubble in shear flow (equation (12)), compared to the drag coefficient of a single bubble in uniform flow (equation (5)), is achieved on the relative velocity of air and water, in the case of Hosokawa and Tomiyama. ⁵ In this case, without implementing a new correlation for the drag coefficient, equation (12), numerical simulations cannot correctly predict the relative velocity profile.

In the case of Hosokawa and Tomiyama ⁵ relative velocity profiles are available, while in the other two cases of Liu ³⁰ and Shawkat et al., ²⁰ velocity profiles of both phases are presented.

Analysing the results of Liu, ³⁰ it is concluded that the air velocity profile depends on the correlation choice for the drag coefficient (equation (5) or equation (12)), while this choice does not affect the water velocity profile.

As expected, in the case of Shawkat et al. ²⁰ the water velocity profiles do not depend on whether equation (5) or equation (12) is used in the calculation. However, in this case, the air velocity profile also does not depend on the drag coefficient modelling. As discussed in the previous chapter, the influence of the new correlation equation (12) on the air velocity profile decreases as the liquid velocity gradient decreases. This is often the case in larger diameter pipes.

Based on all the above, it is recommended to use a new correlation for the drag coefficient equation (12) which considers the influence of shear flow and which is compatible with LRN models. For flows through pipes of relatively small diameter (e.g. Hosokawa and Tomiyama ⁵ and Liu ³⁰ ), the velocity gradient has higher values. In such cases, the use of equation (12) will lead to a better prediction of the air velocity profile (and thus the relative velocity profile). When the flow takes place in pipes of larger diameter (e.g. Shawkat et al. ²⁰ ), the velocity gradient usually has smaller values, so the equation (12) will not lead to a change in the calculation results.

In all calculations, velocity profiles are predicted with satisfactory accuracy. In some places, there are larger deviations in the $k$ and ${\alpha }_{\mathrm{G}}$ profiles compared with experimental results. Different things can lead to that. In numerical simulations, bubbles are represented as monodisperse. In reality, the bubble sizes always differ by some amount. The difference in bubble size can influence flow dynamics. To capture this effect, in future studies, population balance simulations will be conducted. Also, to enable more accurate numerical studies of two-phase air-bubbly flows, additional improvements of interfacial forces models and BIT models are necessary.

The numerical solver created in this study is freely available at the Github: https://github.com/darkoradenkovic/EulerEulerSimulationsInOpenFoam.git.

The developed solver can further be used in the development of more advanced wall functions in the context of turbulent wall-bounded bubbly flows or in turbulent bubbly flow studies where fine resolution of flow quantities is necessary.