Modelling of Bubbly Flow in Bubble Column Reactors with an Improved Breakup Kernel Accounting for Bubble Shape Variations

Abstract

Bubble shapes have been assumed to be spherical in the currently available breakup models such as the one developed by Luo and Svendsen (1996). This particular breakup model has been widely accepted and implemented into CFD modelling of gas-liquid two-phase flows. However, simulation results obtained based on this model usually yield unreliable predictions about the breakage of very small bubbles. The incorporation of bubble shape variation into breakup models has rarely been documented in the study but the bubble shape plays an important role when considering the interactions with the surrounding turbulent eddies in turbulent bubbly flows, especially when the effects of bubble deformation, distortion and bubble internal pressure change are considered during the events of eddy-bubble collision. Thus, the assumption of spherical bubbles seems to be no longer appropriate in reflecting this phenomenon. This study proposes and implements an improved bubble breakup model, which accounts for the variation of bubble shapes when solving the population balance equations for CFD simulation of gas-liquid two-phase flows in bubble columns.

1. Introduction

Previous CFD studies have employed the assumption of a unified bubble diameter, which can generate reliable predictions if the bubble size distribution is very narrow. However, numerical modelling of gas-liquid two-phase flow behaviors should also take into account scenarios where wide bubble size distributions and eddy/bubble-bubble interactions exist. These are very influential factors in the calculation of the gas-liquid interfacial area, which in turn affects the prediction of the mass and heat transfer between the two phases. By solving the population balance equations (PBEs) during the numerical simulation, the bubble size distribution can be derived directly, while the behaviour of the eddy/bubble-bubble interactions can be reflected within coalescence and breakup models.

For the process of bubble breakup, Coulaloglou and Tavlarides [1] assumed that the breakup process would only occur if the energy from turbulent eddies acting on the fluid particle was more than the surface energy it contains. Prince and Blanch [2] acknowledged that bubble breakup is caused by eddy-bubble collision and proposed that bubble breakup can only be induced by eddies with approximately the same characteristic length. For instance, eddies at a much larger length scale transports the bubbles without causing any breakups. Luo and Svendsen [3] described the bubble breakup process by considering both the length scale and the amount of energy contained in the arriving eddies. The minimum length scale of eddies that are responsible for the breakup process is equivalent to 11.4 times the Kolmogorov scale. The critical probability of bubble breakup is related to the ratio of surface energy increase of bubbles after breakup to the mean turbulent kinetic energy of the colliding eddy. Therefore, very small eddies do not contain sufficient energy to cause the bubble breakup process. Lehr et al. [4] proposed a slightly different breakup mechanism from Luo and Svendsen [3] by considering the minimum length scale of eddies to be determined by the size of the smaller bubble after breakup. They also specified that the breakup process is dependent on the inertial force of the arriving eddy and the interfacial force of the bubble. Based on the results of Luo and Svendsen [3] and Lehr et al. [4], Wang et al. [5] proposed an energy constraint and capillary constraint criteria for the breakup model. The energy constraint requires the eddy energy to be greater than or equal to the increase of surface energy of bubbles after the breakage has occurred. The capillary constraint requires the dynamic pressure of the eddy to exceed the capillary pressure of the bubble. The use of these two breakup criteria has restricted the occurrence of breakage that generates unphysically small daughter bubbles and demonstrated more reliable results than that of Luo and Svendsen [3]. Similar ideas to those of Wang et al. [5] have also been adopted by Zhao and Ge [6], Andersson and Andersson [7] and Liao et al. [8]. A more concise breakup constraint of energy density increase was proposed by Han et al. [9]. The constraint of energy density increase involves only one term, which is the energy density itself, to represent what was originally expressed by two terms: capillary pressure and surface energy. It was shown that the energy density increase during the entire breakup process should not exceed the energy density of the parent bubble.

Incorporation of a bubble shape variation into the breakup model has rarely been documented in the open literature. Therefore, the aim of this study is to consider the influence of bubble shape variation on the bubble breakage process in bubble column flows. A breakage model accounting for the variation of bubble shapes, coupled with the breakage criterion of energy density increase, is proposed here.

2. Mathematical modelling

2.1. Bubble size distribution

The bubble size distribution is determined by employing the population balance model with a consideration of bubble coalescence and breakup. Bubbles are divided into several size groups with different diameters specified by the parameter d_eq,i and an equivalent phase with the Sauter mean diameter to represent the bubble classes. In this study, 16 bubble classes with diameters ranging from 1 to 32 mm are applied based on the geometric discretization method where $V_i=2V_{i-1}$. The population balance equation is expressed by Eq. (1)

$$\frac{\partial n_i}{\partial t}+\nabla \cdot \left({{\displaystyle \overrightarrow{v}}}_i\cdot n_i\right)=S_i$$E1

where $n_i$ is the number density for i-th group, ${\overrightarrow{v}}_i$ is the mass average velocity vector and $S_i$ is the source term.

The source term for the i-th group, $S_i$, can be thought of as the birth and death of bubbles due to coalescence and breakup, respectively. The expression for this particular term is given by Eq. (2)

$$\begin{array}{ll}{S}_i\hfill & =B_{\mathit{\text{coalescence}},i}-D_{\mathit{\text{coalescence}},i}+B_{\mathit{\text{breakup}},i}-D_{\mathit{\text{breakup}},i}\\ & =\sum _{d_{\mathit{eq},j}=d_{\mathit{eq},\text{min}}}^{d_{\mathit{eq},i}/2}{\mathrm{\Omega }}_C\left(d_{\mathit{eq},j}:d_{\mathit{eq},i}-d_{\mathit{eq},j}\right)-\sum _{d_{\mathit{eq},j}}^{d_{\mathit{eq},\text{max}}-d_{\mathit{eq},i}}{\mathrm{\Omega }}_C\left(d_{\mathit{eq},j}:d_{\mathit{eq},i}\right)+\sum _{d_j=d_i}^{d_{\max }}{\mathrm{\Omega }}_B\left(d_{\mathit{eq},j}:d_{\mathit{eq},i}\right)-{\mathrm{\Omega }}_B\left(d_{\mathit{eq},i}\right)\hfill \end{array}$$E2

The local gas volume fraction can be calculated using Eq. (3),

$${\alpha }_g{f}_i=n_i{V}_i$$E3

where $f_i$ is the i-th class fraction of the total volume fraction and $V_i$ is the volume for the i-th class.

To describe the coalescence between two bubbles, the coalescence kernel proposed by Luo [10] was utilized in this study. As this is not the main concern of this work, further details of the coalescence kernel can be found in Luo’s paper.

The breakup model proposed in this study is based on the work of Luo and Svendsen [3]. Several improvements have been introduced in this study to produce a more realistic breakup model. In Luo and Svendsen’s model, the shape of breakage bubbles was assumed to be spherical. However, the experimental studies and statistical results, such as Grace et al. [11] and Tomiyama [12], have found that bubbles exist in various shapes and the dynamics of bubble motion strongly depend on the shape of the bubbles. For example, Figure 1 shows the experimentally recorded breakup process of a spherical-cap bubble found in an operating bubble column used in an ongoing research project funded by the Natural Science Foundation of China (NSFC). The spherical-cap bubble is colliding with a bombarding eddy that was generated as a consequence of shedding eddy from the preceding bubbles. The spherical-cap bubble then becomes deformed and distorted and finally breaks into two ellipsoidal bubbles. This phenomenon may lead to two major implications. Firstly, the shed eddies that interact with subsequently formed bubbles are mainly induced by the presence of preceding bubbles. These shed eddies dissipate mainly due to the viscous influence and they will decay downstream in a slightly short distance. Thus, these eddies will have the size the same order as the preceding bubbles. This kind of bubble-induced turbulence may exhibit different dynamic behaviour as can be distinguished from the typical Kolmogorov −5/3 law on the turbulence kinetic energy spectrum. It should be pointed here with caution that more fundamental investigations are required to reveal the interactions between the eddy generated by bubble-induced turbulence and the bubbles, and the impact of this interaction on the bubble breakage process. Secondly, although the bubble shape has been assumed to be spherical in the previous studies for the simplification of models, the variation of bubble shapes could potentially become a critical factor for better prediction of the bubbly flow characteristics of the gas phase in CFD simulations, because the type of geometrical shape has a strong impact on the surface energy of bubbles and interfacial area.

Figure 1.

Time sequence photos of the breakup of a rising bubble in a 150-mm diameter cylindrical bubble column (Ug = 0.02 m/s; total duration: 0.0366 s).

From experimental observations, bubble shapes can be classified into different types. Thus, the effects of different bubble shapes are taken into account in this study. However, due to the uncertainty of the spatial orientation of the bubbles during their movement, the determination of the contact angle of the bombarding eddy is very difficult but this needs to be tackled as the contact angle will directly affect the projection/sweep area of the eddy-bubble collision tube. On the contrary, if the bubble that is about to breakup is assumed to be spherical, the projection/sweep area of the collision tube will be consistent no matter which direction of the bombarding eddy comes from. Instead of using the original bubble size d_eq,i to construct the collision tube, a nominal diameter, d_V, that approximately represents the size of the projected area of the bubble, is defined in a bounded range given by expression (4),

$$c\le d_V\le a$$E4

where c and a are the length of the short axis and long axis, respectively. The new imaginary collision tube is presented in Figure 2.

Figure 2.

Diagram showing an eddy entering a collision tube and moving through it with a mean velocity.

The breakup rate for one individual parent bubble that forms into two daughter bubbles can be calculated using Eq. (5),

$${\mathrm{\Omega }}_B={\int }_{{\lambda }_{\min }}^d{\omega }_B^T{p}_B\mathit{d\lambda }$$E5

where ${\omega }_B^T$ is the collision probability density, which can be estimated from Luo and Swendsen [3], as defined by Eq. (6)

$${\omega }_B^T\left(\xi \right)=0.923\left(1-{\alpha }_g\right){\left(\epsilon d_{\mathit{eq},i}\right)}^{1/3}{n}_i\frac{\left(d_{V,s}/d_{\mathit{eq},i}+\xi \right)\left(d_{V,l}/d_{\mathit{eq},i}+\xi \right)}{d_{\mathit{eq},i}^2{\xi }^{11/3}}$$E6

Here, ξ =λ /d_eq,i is the non-dimensional size of eddies that may contribute to the breakage of bubbles with size d_i. The breakage probability function $p_B$ used by Luo and Svendsen [3] is given by Eq. (7),

$$p_B=\text{exp}\left(-\frac{e_s}{{\displaystyle \overline{e}}}\right)$$E7

where $\overline{e}$ is the mean turbulent kinetic energy for eddies of size $\lambda$ and $e_s$ is the increase in surface energy of bubbles after breakage. The mean turbulent kinetic energy can be determined by Eq. (8)

$$\overline{e}={\rho }_l\frac{\pi }{6}{\lambda }^3\frac{{\overline{u}}_{\lambda }^2}{2}=\frac{\mathit{\pi \beta }}{12}{\rho }_l{\left(\epsilon d_{\mathit{eq},i}\right)}^{2/3}{d}_{\mathit{eq},i}^3{\xi }^{11/3}$$E8

By assuming the bubbles before and after breakage have deformed shapes with an equivalent diameter, when the parent bubble of size d_eq,i breaks into two bubbles of size d_eq,j and (d_eq,i³-d_eq,j³)^1/3, the increase in surface energy can be estimated using Eq. (9),

$$e_s\left(d_{\mathit{eq},i},d_{\mathit{eq},j}\right)=\sigma \cdot \pi d_{\mathit{eq},i}^2\left[f_V^{2/3}+{\left(1-f_V\right)}^{2/3}-1\right]$$E9

where the breakage volume fraction is given by $f_V=d_{\mathit{eq},j}^3/d_{\mathit{eq},i}^3$. Since the effects of different shapes of bubbles are now taken into account, Eq. (9) can be rewritten in a general form in terms of the surface area of the bubbles, $S$, as defined by Eq. (10)

$$e_s=\sigma \cdot \left(S_{j,1}+S_{j,2}-S_i\right)$$E10

Although there have been some recent developments on the instability of bubble shapes, such as the studies by Cano-Lozano et al. [13], Zhou and Dusek [14] and Tripathi et al. [15], there is no consensuses on concise definitions on bubble shapes and bubble shape model. Therefore, a more commonly accepted statistical model of bubble shapes by Tomiyama [12] has been employed in this study. In addition, the lift model described by Tomiyama [12] has been adopted as it has been well implemented by different commercial CFD packages. According to the criteria proposed by Tomiyama [12] and Tomiyama et al. [16], there are three main types of bubble shapes that should be considered in the bubble columns of this study. These shapes include spherical, ellipsoidal and spherical-capped bubbles. These three types of bubble shapes may also be considered for modelling gas-liquid two-phase flow or gas-liquid-solid three-phase flow in bubble columns with similar scales that operate at similar conditions to what is applied in this work. The details of these three types of bubble shapes and their potential breakage scenarios are depicted in Figure 3.

Figure 3.

Classification of the three types of bubble shapes and the possible breakage scenarios.

For an air-water system under atmospheric pressure and room temperature conditions, the size boundary to categorize between spherical and ellipsoidal bubbles represented by d_eq,1 is roughly 1.16 mm for the pure system and approximately 1.36 mm for a slightly contaminated system. It is very important to point out that the volumes of ellipsoidal bubbles and spherical-cap bubbles should be equal to the volumes of their equivalent spherical bubbles with diameter d_eq. For bubbles with ellipsoidal shapes, by assuming an oblate type of ellipsoid, the surface area can be calculated using Eq. (11),

$$S_{\mathit{\text{ellipsoid}}}=\frac{\pi }{2}{d}_{\mathit{eq}}^2{E}^{1/3}\left(1+\frac{1}{2E\sqrt{E^2-1}}\text{ln}\left(2E^2-1+2E\sqrt{E^2-1}\right)\right)$$E11

where the aspect ratio E can be expressed using the empirical correlation described by Wellek et al. [17], which is given by Eq. (12)

$$E=a/b=1+0.163{\mathit{Eo}}^{0.757}$$E12

Here, $\mathit{Eo}$ is the Eötvös number as defined by Eq. (13)

$$\mathit{Eo}=\frac{g\left({\rho }_l-{\rho }_g\right)d_{\mathit{eq}}^2}{\sigma }$$E13

The size boundary to divide between ellipsoidal and spherical-cap bubbles represented by d_C is estimated using Eq. (14),

$$d_C=\sqrt{\raisebox{1ex}{$40\sigma $}\!\left/ \!\raisebox{-1ex}{$g\left({\rho }_L-{\rho }_g\right)$}\right.}$$E14

where d_c is determined to be 17.3 mm for the air-water system. For a single spherical-cap bubble, the wake angle ${\theta }_W$ is assumed to be 50^o based on the work of Tomiyama [12]. As the volume of the spherical-cap bubble is equivalent to the volume of the spherical bubble, Eq. (15) can be formulated as follows:

$$R_S^3=\frac{\pi d_{\mathit{eq}}^3/6}{{\left(1-\text{cos}{\theta }_W\right)}^2-{\left(1-\text{cos}{\theta }_W\right)}^3/3}$$E15

The curved surface area for the front edge can be calculated from the following relationship given by Eq. (16):

$$S_{\mathit{Cap}}=2\pi R^2\left(1-\text{cos}{\theta }_W\right)$$E16

The experimental observations of Davenport et al. [18] and Landel et al. [19] have clearly indicated that the rear surface of a single spherical-cap bubble follows a constantly oscillating lenticular shape, resulting from external perturbations acting on the rear surface. This lenticular shaped rear surface can be considered to be essentially flat, and the surface energy increase required to breakup the surface can be neglected based on the consideration that when any arriving eddies bombard the flat surface, the energy resulting from the surface tension force action will be far smaller than the kinetic energy exuded by the turbulent eddies. It should be noted with caution that these are rough approximations, and more complicated crown bubble systems are not considered in this work. The influence of the variation of bubble shapes on the increase in surface energy is further illustrated in Figure 7.

The breakup model proposed by Luo and Svendsen [3] only considered the surface energy requirement for breakup events but it should be noted that bubble breakage may also be subjected to the pressure head difference of the bubble and its surrounding eddies, especially when the breakage volume fraction is small. Therefore, on the basis of the interaction force balance proposed by Lehr et al. [4], the pressure energy requirement also needs to be considered as a competitive breakup mechanism. This can be imposed as a constraint. The same idea has been adopted by Zhao and Ge [6], Liao et al. [8] and Guo et al. [20]. The pressure energy requirement can be expressed using Eq. (17),

$$e_P=\frac{\sigma }{\text{min}\left(R_{C,j},R_{C,k}\right)}\cdot \frac{\pi {\left(\text{min}\left(d_{\mathit{eq},j},d_{\mathit{eq},k}\right)\right)}^3}{6}$$E17

where R_C,j and R_C,k are the equivalent radius of curvature of daughter bubbles. The theoretical prediction of the surface energy and pressure energy requirement is shown in Figure 6.

As pointed out by Han et al. [21], from a volume-based energy perspective, the surface energy density of the parent bubble should always exceed the maximum value of the energy density increase during the entire breakup process. This is an important breakup criterion that has been adopted in this study and concurrently relates the size of the parent bubble to the sizes of the daughter bubbles. This restricts the generation of very small bubbles from the breakup process because the energy densities of the daughter bubbles will tend towards infinity when their sizes tend to zero. The energy density criterion can be expressed by Eq. (18) if it is coupled with the variation of bubble shapes

E18

The breakup frequency can be obtained by substituting Eqs. (6)–(17) into Eq. (5), which results in Eq. (19),

$$\begin{array}{l}{\mathrm{\Omega }}_B=\left\{\begin{array}{l}0.923\left(1-{\alpha }_g\right)n_i{\left(\epsilon /d_{\mathit{eq},i}^2\right)}^{1/3}\cdot {\int }_{{\xi }_{\min }}^1\frac{\left(d_{V,s}/d_{\mathit{eq},i}+\xi \right)\left(d_{V,l}/d_{\mathit{eq},i}+\xi \right)}{{\xi }^{11/3}}\text{exp}\left(-\frac{12\sigma \left(S_j+S_k-S_i\right)}{{\mathit{\pi \beta \rho }}_l{\epsilon }^{2/3}{\xi }^{11/3}{d}_{\mathit{eq},i}^{11/3}}\right)\mathit{d\xi },\\ \mathit{\text{when}}\frac{6\sigma \left(S_j+S_k-S_i\right)}{\pi d_{\mathit{eq},i}^3}\ge \frac{\sigma }{\text{min}\left({\mathit{Rc}}_j,{\mathit{Rc}}_k\right)}\\ 0.923\left(1-{\alpha }_g\right)n_i{\left(\epsilon /d_{\mathit{eq},i}^2\right)}^{1/3}\cdot {\int }_{{\xi }_{\min }}^1\frac{\left(d_{V,s}/d_{\mathit{eq},i}+\xi \right)\left(d_{V,l}/d_{\mathit{eq},i}+\xi \right)}{{\xi }^{11/3}}\text{exp}\left(-\frac{2\sigma {\left(\text{min}\left(d_{\mathit{eq},j},d_{\mathit{eq},k}\right)\right)}^3}{\text{min}\left({\mathit{Rc}}_j,{\mathit{Rc}}_k\right){\mathit{\beta \rho }}_l{\epsilon }^{2/3}{\xi }^{11/3}{d}_{\mathit{eq},i}^{11/3}}\right)\mathit{d\xi },\\ \mathit{\text{when}}\frac{6\sigma \left(S_j+S_k-S_i\right)}{\pi d_{\mathit{eq},i}^3}<\frac{\sigma }{\text{min}\left({\mathit{Rc}}_j,{\mathit{Rc}}_k\right)}\end{array}\right.\end{array}$$E19

where ${\xi }_{\min }$ is the minimum breakage volume fraction that is able to satisfy the energy density criterion shown in Eq. (18).

2.2. Governing Eqs

A three-dimensional (3D) transient CFD model is employed in this work to simulate the local hydrodynamics of the gas-liquid two-phase bubble column. An Eulerian-Eulerian approach is adopted in order to describe the flow behaviors for both phases, that is, water as the continuous phase and air as the dispersed phase.

The mass and momentum balance equations are given by Eqs. (20) and (21), respectively,

$$\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\rho }_k{\alpha }_k\right)\mathbf{+}\mathbf{\nabla }\mathbf{\bullet }\left({\rho }_k{\alpha }_k{\mathit{u}}_k\right)\mathbf{=}\mathbf{0}$$E20

$$\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\rho }_k{\alpha }_k{\mathit{u}}_k\right)\mathbf{+}\mathbf{\nabla }\mathbf{\bullet }\left({\rho }_k{\alpha }_k{\mathit{u}}_k{\mathit{u}}_k\right)\mathbf{=}\mathbf{-}{\alpha }_k\mathbf{\nabla }p\mathbf{+}\mathbf{\nabla }\mathbf{\bullet }{\mathit{\tau }}_k+{\alpha }_k{\rho }_k\mathit{g}+{\mathit{F}}_k$$E21

where ${\rho }_k$, ${\alpha }_k$, u_k, τ_k and F_k represent the density, volume fraction, velocity vector, viscous stress tensor and the interphase momentum exchange term for the k (liquid or gas) phase, respectively. The sum of the volume fractions for both phases is equal to 1.

A modified $k˜\epsilon$ turbulence model with the consideration of bubble-induced turbulence by Sato and Sekoguchi [22] is used for turbulence closure. The turbulent kinetic energy $k_l$ and dissipation rate ${\epsilon }_l$ are computed using Eqs. (22) and (23),

$$\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\rho }_i{\alpha }_i{k}_i\right)+\mathbf{\nabla }\bullet \left({\rho }_i{\alpha }_i{k}_i{\mathit{u}}_i\right)=\mathbf{\nabla }\bullet \left[{\alpha }_{\mathit{i}}\left({\mathit{\mu }}_{\mathit{i}}+\frac{{\mathit{\mu }}_{\mathit{eff},\mathit{i}}}{{\mathit{\sigma }}_{\mathit{k}}}\right)\mathbf{\nabla }{k}_{\mathit{i}}\right]+{\alpha }_i\left(G_{k,i}-{\mathrm{\rho }}_i{\epsilon }_i\right)$$E22

$$\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\rho }_i{\alpha }_i{\epsilon }_i\right)\mathbf{+}\mathbf{\nabla }\mathbf{\bullet }\left({\rho }_i{\alpha }_i{\mathit{\epsilon }}_i{\mathit{u}}_i\right)\mathbf{=}\mathbf{\nabla }\mathbf{\bullet }\left[{\alpha }_{\mathit{i}}\left({\mathit{\mu }}_{\mathit{i}}\mathbf{+}\frac{{\mathit{\mu }}_{\mathit{eff}\mathbf{,}\mathit{i}}}{{\mathit{\sigma }}_{\mathit{\epsilon }}}\right)\mathbf{\nabla }{\epsilon }_{\mathit{i}}\right]+{\alpha }_i\frac{{\epsilon }_i}{k_i}\left(C_{1\epsilon }{G}_{k,i}-C_{2\epsilon }{\mathrm{\rho }}_i{\epsilon }_i\right)$$E23

where $G_{k,l}$ is the production of turbulent kinetic energy and ${\mu }_{t,l}$ is the turbulent viscosity. In this work, the standard $k˜\epsilon$ model constants used are $C_{\mu }$ = 0.09, $C_{1\epsilon }$ = 1.44, $C_{2\epsilon }$ = 1.92, ${\sigma }_k$ = 1.0, ${\sigma }_{\epsilon }$ = 1.3.

The effective viscosity is composed of the contributions of turbulent viscosity and an extra term considering the effect of bubble-induced turbulence and is defined by Eq. (24)

$${\mathrm{\mu }}_{\mathit{eff},l}={\rho }_l{C}_{\mathrm{\mu }}\frac{k_l^2}{{\epsilon }_l}+{\rho }_l{C}_{\mathrm{\mu },\mathit{BIT}}{\alpha }_g{d}_b\left|{\mathit{u}}_g-{\mathit{u}}_l\right|$$E24

The Sato coefficient $C_{\mu ,\mathit{BIT}}=0.6$ is adopted according to the study [22].

2.3. Interphase momentum transfer

In this study, drag force, lift force and added mass force are considered as the main interactions between the continuous liquid phase and the dispersed gas phase. The drag force is calculated using Eq. (25),

$${\mathit{F}}_D=\frac{3}{4}\frac{C_D}{d_{\mathit{eq}}}{\rho }_l{\alpha }_g\left|{\mathit{u}}_g-{\mathit{u}}_l\right|\left({\mathit{u}}_g-{\mathit{u}}_l\right)$$E25

where $C_D$ is the drag coefficient, which can be obtained from the model by Grace et al. [11]. The Grace model is well suited for gas-liquid flows in which the bubbles exhibit a range of shapes, such as sphere, ellipsoid and spherical-cap. However, instead of comparing the values of drag coefficients in the original Grace model, the drag coefficient can be applied directly into the present model as the variation of bubble shapes has been taken into account. The drag coefficients for the different shapes of bubbles are calculated using Eqs. (26)–(28),

$$C_{D,\mathit{\text{sphere}}}=\left\{\begin{array}{l}24/{\mathit{Re}}_b{\mathit{Re}}_b<0.01\\ 24\left(1+0.15{\mathit{Re}}_b^{0.687}\right)/{\mathit{Re}}_b{\mathit{Re}}_b\ge 0.01\end{array}\right.$$E26

$$C_{D,\mathit{cap}}=\frac{8}{3}$$E27

$$C_{D,\mathit{\text{ellipse}}}=\frac{4}{3}\frac{{\mathit{gd}}_{\mathit{eq}}}{U_t^2}\frac{\left({\rho }_l-{\rho }_g\right)}{{\rho }_l}$$E28

where Re_b is the bubble Reynolds number given by ${\mathit{Re}}_b=\frac{{\rho }_l\left|{\mathit{u}}_{\mathit{g}}-{\mathit{u}}_l\right|d_{\mathit{eq}}}{{\mu }_l}$ and $U_t$ is the terminal velocity calculated using Eq. (29),

$$U_t=\frac{{\mu }_l}{{\rho }_{l}d}{\mathit{Mo}}^{-0.149}\left(J-0.857\right)$$E29

Here, Mo is the Morton number defined by $\mathit{Mo}=\frac{{\mu }_l^{4}g\left({\rho }_l-{\rho }_g\right)}{{\rho }_l^2{\sigma }^3}$ and J is determined by the piecewise function calculated using the empirical expression (30)

$$J=\left\{\begin{array}{l}0.94H^{0.757}2<H<59.3\\ 3.42H^{0.441}H>59.3\end{array}\right.$$E30

H in expression (30) is defined by Eq. (31),

$$H=\frac{4}{3}{\mathit{\text{EoMo}}}^{-0.149}{\left(\frac{{\mu }_l}{{\mu }_{\mathit{ref}}}\right)}^{-0.14}$$E31

where Eo is the Eötvös number and ${\mu }_{\mathit{ref}}=0.0009\mathit{kg}/\left(m\cdot s\right)$.

The lift force acting perpendicular to the direction of relative motion of the two phases can be calculated by using Eq. (32)

$${\mathit{F}}_{\mathit{\text{Lift}}}=C_L{\rho }_l{\alpha }_g\left({\mathit{u}}_g-{\mathit{u}}_l\right)\times \left(\nabla \times {\mathit{u}}_l\right)$$E32

where $C_L$ is the lift coefficient and is estimated by the Tomiyama lift force correlation [12], as described by the following empirical relation (33),

$$C_L=\left\{\begin{array}{l}\text{min}\left[0.288\text{tanh}\left(0.121{\mathit{Re}}_b\right),f\left({\mathit{Eo}}^{\text{'}}\right)\right]{\mathit{Eo}}^{\text{'}}\le 4\\ f\left({\mathit{Eo}}^{\text{'}}\right)4<{\mathit{Eo}}^{\text{'}}<10\\ -0.29{\mathit{Eo}}^{\text{'}}>10\end{array}\right.$$E33

where $f\left({\mathit{Eo}}^{\text{'}}\right)=0.00105\mathit{Eo}{\text{'}}^3-0.0159\mathit{Eo}{\text{'}}^2-0.0204{\mathit{Eo}}^{\text{'}}+0.474$. Eo’ is the modified Eötvös number based on the maximum horizontal dimension of the deformable bubble, d_h, as defined and given, respectively, by Eqs. (34) and (35)

$${\mathit{Eo}}^{\text{'}}=\frac{g\left({\rho }_l-{\rho }_g\right)d_h^2}{\sigma }$$E34

$$d_h=d{\left(1+0.163{\mathit{Eo}}^{0.757}\right)}^{1/3}$$E35

The virtual mass force is also significant when the gas phase density is smaller than the liquid phase density. The estimation of the virtual mass force due to the deformation of bubbles is one of the unresolved issues that require further investigation. With the caution, the virtual mass force is still calculated using Eq. (36),

$${\mathit{F}}_{\mathit{VM}}=C_{\mathit{VM}}{\rho }_l{\alpha }_g\left(\frac{d{\mathit{u}}_{\mathit{l}}}{\mathit{dt}}-\frac{d{\mathit{u}}_{\mathit{g}}}{\mathit{dt}}\right)$$E36

where C_VM is the virtual mass coefficient defined as 0.5 in this study.

2.4. Numerical modelling

To validate the influence of variations in bubble shapes considered in the breakup model, numerical simulations have been carried out for the air-water bubble columns used by Kulkarni et al. [23] and Camarasa et al. [24] denoted by Case 1 and Case 2, respectively, in Table 1.

	Diameter (m)	Height (m)	Superficial velocity (m/s)	Static liquid height (m)
Case 1	0.15	0.8	0.0382	0.65
Case 2	0.1	2	0.0606	0.9

Table 1.

Details of the experimental setup.

The mesh setup is illustrated in Figure 4. Grid 2 consists of $20\left(r\right)\times 40\left(\theta \right)\times 100\left(z\right)$ nodes in the radial, circumferential and axial directions, respectively. The grid independence was tested in a coarser Grid 1 of $16\left(r\right)\times 32\left(\theta \right)\times 80\left(z\right)$ nodes and a refined Grid 3 of $26\left(r\right)\times 48\left(\theta \right)\times 126\left(z\right)$ nodes, in which case the total number of cells is doubled gradually. The grid independence test for these three setups has yielded similar results quantitatively, even though the overall trend of overprediction occurred for all three grids, as shown in Figure 5. Grid 2 was chosen and used in subsequent simulations to investigate the effects of the improved breakup model.

Figure 4.

Horizontal cross section and front view of the mesh setup for the main body of the bubble column.

Figure 5.

Comparison of the simulated gas holdup profile to the data reported in Camarasa et al. [24] with three different grid configurations.

ANSYS Fluent 3D pressure-based solver is employed in CFD-PBM modelling. The time step is set to be 0.001 s for all simulations, which is considered to be sufficient for illustrating the time-averaged characteristics of the flow fields by carrying out the data-sampling statistics for typically 120 s after the quasi-steady state has been achieved. The improved breakup model is integrated into the simulations by using the user define function (UDF). At the inlet boundary, the volume fraction of gas phase is set to be 1. The treatment of the inlet velocity is different from using a constant superficial gas velocity, but a normally distributed velocity profile is applied by using the model proposed by Shi et al. [25], which can accurately reflect the experimental conditions employed in the study by Camarasa et al. [24]. Further information about the reasons, theoretical basis and the effects of using the inlet model can be found in their published work. The outlet boundary is set to be a pressure outlet at the top. No-slip conditions are applied for both the liquid and gas phases at the bubble column wall.

3. Results and discussion

3.1. Effect of deformed bubble shape variations on the pressure and surface energy required for bubble breakage

To illustrate the influence of pressure energy control breakup, theoretical predictions of the surface energy and the pressure energy requirements for the breakage of ellipsoidal and spherical-capped bubble are shown, respectively, in Figure 6. It can be clearly seen from Figure 6 that the energy requirement for ellipsoid bubble shifts from pressure energy to surface energy with an increase in the breakup volume fraction. This may be attributed to a higher dynamic pressure being required inside a smaller bubble for resisting the surrounding eddy pressure in order to sustain its own existence. However, the spherical bubble requires most of the surface energy for its breakage. This may mainly be due to the contribution of the large front surface of spherical-capped bubbles.

Figure 6.

Two competitive control mechanisms of the breakage of ellipsoidal bubbles (d_eq,i = 16 mm) and spherical-capped bubbles (d_eq,i = 32 mm).

The surface energy requirement for bubble breakage in Figure 6 has taken into account the bubble shape variations. To further illustrate the significance of considering the variation of bubble shapes, a theoretical comparison of the increase in surface energy for the breakage of the original spherical bubbles versus various shapes of bubbles has been shown in Figure 7. The generation of spherical bubbles due to eddy collision with large ellipsoidal or spherical-capped bubble is not covered, as the breakage volume fraction will be far smaller than 0.05. The generation of small spherical bubbles occurs more frequently due to the interaction of the shed eddies with the bubble skirt. This phenomenon was concisely described and explained by numerical modelling work carried out by Fu and Ishii [26]. It is shown in Figure 7 that the maximum increase in surface energy for ellipsoidal bubbles and spherical-capped bubbles is different. As binary breakage is assumed, a large ellipsoidal bubble breaks into two smaller ellipsoidal bubbles in most cases. The maximum increase in surface energy is demonstrated when equal-size breakage occurs, which suggests that the parent ellipsoidal bubble has been through a large deformation process itself. However, the spherical-capped bubble can break into different combinations of daughter bubble types, including one ellipsoidal and one spherical-capped bubble, two ellipsoidal bubbles, or two spherical-capped bubbles. The maximum increase in surface energy for the breakage of a spherical-capped parent bubble is found with the largest volume fraction of ellipsoidal daughter bubble. This result coincides with the existing experimental observations: the ellipsoidal bubble has a more stable structure that is able to resist bombarding eddies from both the front and the rear, whereas the spherical-capped bubble can only resist eddies hitting from the front but is easily and rapidly ruptured by eddies hitting from the rear.

Figure 7.

Normalized increase in the surface energy for the breakage of original spherical bubbles and various shapes of bubbles.

Figure 8 compares the time-averaged gas holdup predicted by the original breakup model and the improved breakup model. It can be found that the improved breakup model has achieved results very similar to the experimental data at the core region of the column (r/R < 0.6), while underestimation is shown near the column wall for both the original breakup model and the improved breakup model. Since the standard $k˜\epsilon$ turbulence model is still applied in this study, the underestimation of gas holdup may be due to the slight poor prediction of the turbulence dissipation rate. The issue of underestimation on the gas holdup distribution has also been addressed by Chen et al. [27], in which case the breakup rate was artificially increased by a factor of 10 to obtain a “better” agreement with the experimental data.

Figure 8.

Comparison of the original breakup model and the improved breakup model for the prediction of the gas holdup profile in the radial direction for Case 1.

Figure 9 shows the radial distribution of the time-averaged turbulence dissipation rate for Case 1. The turbulence dissipation rate distribution predicted by the standard $k˜\epsilon$ model is smaller than the result obtained by the RNG $k˜\epsilon$ model. This is because the RNG $k˜\epsilon$ model has a specific contribution from the local strain rate as the correction to the turbulence dissipation rate. The tendency of the standard $k˜\epsilon$ model to underestimate the turbulence dissipation rate can also be seen in the studies carried out by Laborde-Boutet et al. [28], Chen [29] and Jakobsen et al. [30]. As a result, the standard $k˜\epsilon$ model is insufficient to properly estimate the turbulence dissipation rate in the regions with rapidly strained flows, which most likely corresponds to the near wall region in the bubble columns. It can be seen from Eq. (19) that the breakup rate ${\mathrm{\Omega }}_B˜{\epsilon }^{1/3}\cdot \text{exp}\left(-{\epsilon }^{-2/3}\right)$, which is at least equivalent to the dissipation rate ε of the order of −1/3. Therefore, the equilibrium state of bubble coalescence and breakup phenomena cannot be reasonably addressed with an inaccurate estimation of the turbulence dissipation rate and inevitably affect the predictions of gas holdup. Also, as the predicted coalescence rate is about one order of magnitude higher than the predicted breakup rate, the bubble coalescence and breakup phenomena cannot be reasonably addressed under this scenario and will inevitably affect the predictions of gas holdup. In addition, as pointed out by Jakobsen et al. [30], despite the accuracy of calculating the local turbulence dissipation rate from the $k˜\epsilon$ turbulence model, this turbulence dissipation rate merely represents a fit of a turbulence length scale to single-phase pipe flow data. Therefore, the contribution of turbulence eddies that are induced by the bubbles has not being included. More importantly, the mechanism of bubble breakage caused by the interactions of bubble-induced turbulence eddies with the subsequent bubbles, which may be dominant in the core region of the bubble column, cannot be revealed through the breakage kernels that are very sensitive to the turbulence dissipation rate.

Figure 9.

Radial distribution of time-averaged turbulence dissipation rate for Case 1.

Figure 10 shows the radial distribution of time-averaged gas holdup at different cross sections in the axial direction. The results are obtained by using the improved breakup model. It can be seen clearly from Figure 10 that the predicted time-averaged gas holdup in the fully developed region (H/D > 5) has achieved self-preserving characteristics regardless of the axial positions. It appears that the inlet conditions have a weak influence on this self-preserving nature in the bubble columns, which is a result concurring with some previous experimental findings [31, 32].

Figure 10.

Radial distribution of time-averaged gas holdup at different axial positions.

Figure 11 presents the unit volume-based interfacial area in the bulk region for each bubble class. Due to the large differences in size from the smallest to the largest bubble class, the y-axis is shown in a log₁₀ scale. Interfacial area is a key parameter that largely affects the prediction of heat and mass transfer between gas and liquid phase in the bubble columns. Although the differences in the simulated interfacial area between the improved breakup model and the original breakup model are not significant when the bubble size is relatively small, the influence of the bubble shapes is gradually reflected when the shape of the bubbles transforms from ellipsoid to spherical-cap, resulting in much larger interfacial areas for spherical-capped bubbles.

Figure 11.

Comparison of the simulated interfacial area in the bubble column for Case 2.

3.2. Effect of deformed bubble shape variations on the interfacial mass transfer across bubble surfaces

The interfacial area obtained by the improved breakup model is based on the statistical model of bubble shapes. The results will be slightly different when a more realistic model, which considers the dynamic deformations that occur during bubble motions, is implanted into the simulations. Indeed, the current results have implied that assuming all bubbles to be of a spherical shape may lead to significant underestimation of the interfacial area and hence affect the predictions of the heat and mass transfer rate when chemical reactions are considered in the bubble column reactors. To further address this issue, the volumetric mass transfer coefficient, k_La, estimated based on the improved breakup model for each bubble class is presented in Figure 12.

Figure 12.

Comparison of the volumetric mass transfer coefficient for each bubble class.

The convective mass transfer film coefficient can be defined by Eq. (37)

$$k_L=\frac{\overline{D}}{d}\mathit{Sh}$$E37

where $\overline{D}$ is mass diffusivity, d is the bubble diameter and Sh is the Sherwood number. The Sherwood number represents the ratio of the convective mass transfer to the rate of diffusive mass transfer. It can be determined by using the Frossling equation described by Eq. (38)

$$\mathit{Sh}=2+0.552{\text{Re}}^{\frac{1}{2}}{\mathit{Sc}}^{\frac{1}{3}}$$E38

Re in Eq. (38) is the bubble Reynolds number and Sc is the Schmidt number. The Schmidt number is the ratio of momentum diffusivity to mass diffusivity, defined by Eq. (39)

$$\mathit{Sc}=\frac{\upsilon }{\overline{D}}$$E39

According to the analogy between heat and mass transport phenomena, a similar method can be applied to calculate the Nusselt number by simply replacing the Schmidt number with the Prandtl number. By doing so, the ratio of convective heat transfer to conductive heat transfer can be characterized.

It is observed that the volumetric mass transfer coefficient is greatly increased due to the contribution of ellipsoidal and spherical-capped bubbles. However, the peak value obtained based on the improved breakup model may be attributed to the predicted number density of the corresponding bubble class. As illustrated in Figure 7, the improved breakup model requires a higher increase in surface energy at the boundary of ellipsoidal and spherical-capped bubbles, which makes the smallest spherical-capped bubbles more difficult to break. The results for this particular bubble class may not be a good reflection of the physical phenomenon in reality, but the overall enhancement of the mass transfer coefficient is still very significant. The predictions on the overall mass transfer coefficient are shown in Figure 13. Figure 14 presents the local mass transfer coefficient at different cross sections along the height of the bubble column. It can be seen from Figure 14 that the mass transfer rate estimations based on Luo and Svendsen model and the improved breakup model are obviously very different. The results based on the Luo and Svendsen model may imply that the mass transfer is mainly associated with the regions where the larger Sauter mean bubble diameter has been predicted. The results based on the improved breakup model suggest that the mass transfer is more uniformly distributed, in which case the enhanced overall mass transfer estimation comes from the statistical sum of the contributions of each bubble class.

Figure 13.

Volumetric mass transfer coefficient predicted using both the improved breakup model and Luo and Svendsen model [3].

Figure 14.

Distribution of estimated volumetric mass transfer coefficient at different cross sections in the bubble column for Case 2. (a) Luo and Svendsen model; (b) improved breakup model. (from top to bottom: H = 0.6, 0.5, 0.4, 0.3 and 0.2 m.).

4. Conclusion

In this study, an improved breakup model has been proposed based on the breakup model by Luo and Svendsen [3]. This improved breakup model takes into account the variation of bubble shapes in bubble columns, which include spherical, ellipsoid and spherical-cap shaped bubbles. In addition, the model considers the pressure energy controlled breakup coupled with modified breakage criteria. The simulation results demonstrate an overall agreement with the experimental data reported in the open literature. The difference between the surface energy and the pressure energy requirements for forming various daughter bubbles has been illustrated. The energy density constraint has been applied to prevent overestimating the breakage rate of small bubbles. This study on the dynamic behaviour of various bubble shapes could potentially lead to a more comprehensive understanding of the mass and heat transfer characteristics of multiphase flows in the bubble column.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 91534118). Weibin Shi would also like to acknowledge the Ph.D. scholarship of the International Doctoral Innovation Centre (IDIC) of University of Nottingham Ningbo China and the support of EPSRC (Grant no. EP/G037345/1).

This chapter is an extension of a conference paper that was presented at the 13th International Conference on Heat Transfer, Fluid Mechanics and Thermodynamics (HEFAT2017) and was nominated for invitation into the HEFAT2017 Special Issue of Heat Transfer Engineering based on its designation as a high-quality paper of relevance to the modelling of fluids based systems.

Nomenclature

a	long half axis length of a ellipse, m
c	short half axis length of a ellipse, m
CD	effective drag coefficient for a bubble around a swarm, dimensionless
D	bubble column diameter, m
D¯	mass diffusivity, m2/s
d	bubble diameter, m
deq	equivalent bubble diameter, m
dV	length of virtual axis, m
Eo	Eötvös number, dimensionless
ē	mean turbulence kinetic energy, kg·m2/s2
es	increase in surface energy, kg·m2/s2
FD	drag force, N/m3
FLift	lift force, N/m3
FVM	virtual mass force, N/m3
fV	breakage volume fraction, dimensionless
g	gravity acceleration, m/s2
H	distance from the bottom surface, m
k	turbulence kinetic energy, m2/s2
kL	convective mass transfer film coefficient, m/s
Mo	Morton number, dimensionless
n	number density per unit volume, m−3
t	time, s
Rc	radius of curvature, m
Re	Reynolds number, dimensionless
S	surface area, m2
Sh	Sherwood number, dimensionless
Sc	Schmidt number, dimensionless
U	superficial velocity, m/s
Ut	terminal velocity, m/s
ūλ	mean velocity of turbulence eddies, m/s
u	velocity vector, m/s
V	volume, m3
α	phase volume fraction, gas holdup
ε	turbulence dissipation rate, m2/s3
λ	characteristic length scale of eddy, m
μ	molecular dynamic viscosity, Pa·s
μeff	effective turbulence dynamic viscosity, Pa·s
υ	kinematic viscosity, m2/s
ρ	fluid density, kg/m3
σ	surface tension, N/m
τ	shear stress, Pa
b	bubble
g	gas
i	i-th class bubble
j/k	daughter bubble
l	liquid