Details of the experimental setup.

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

### Keywords

- bubble column CFD simulation
- breakup model
- bubble shape variations
- interfacial area
- mass transfer coefficient

## 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 *deq,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

where *i*-th group,

The source term for the *i*-th group,

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

where *i*-th class fraction of the total volume fraction and *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.

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 *deq,i* to construct the collision tube, a nominal diameter, *dV*, that approximately represents the size of the projected area of the bubble, is defined in a bounded range given by expression (4),

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.

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

where

Here, *ξ* =*λ* /*deq,i* is the non-dimensional size of eddies that may contribute to the breakage of bubbles with size *d*_{i.} The breakage probability function

where

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

where the breakage volume fraction is given by

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.

For an air-water system under atmospheric pressure and room temperature conditions, the size boundary to categorize between spherical and ellipsoidal bubbles represented by *deq,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 *deq*. For bubbles with ellipsoidal shapes, by assuming an oblate type of ellipsoid, the surface area can be calculated using Eq. (11),

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

Here,

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

where *dc* is determined to be 17.3 mm for the air-water system. For a single spherical-cap bubble, the wake angle ^{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:

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

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),

where *RC,j* and *RC,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

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

where

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

where **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

where

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)

The Sato coefficient

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

where

where *Reb* is the bubble Reynolds number given by

Here, *Mo* is the Morton number defined by *J* is determined by the piecewise function calculated using the empirical expression (30)

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

where *Eo* is the Eötvös number and

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

where

where *Eo’* is the modified Eötvös number based on the maximum horizontal dimension of the deformable bubble, *dh*, as defined and given, respectively, by Eqs. (34) and (35)

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),

where *CVM* 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 |

The mesh setup is illustrated in Figure 4. Grid 2 consists of

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.

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

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

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 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_{10} 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.

### 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, *kLa*, estimated based on the improved breakup model for each bubble class is presented in Figure 12.

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

where *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)

*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)

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.

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