2.1. Bubble size distribution

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

66 Heat and Mass Transfer - Advances in Modelling and Experimental Study for Industrial Applications

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

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

within coalescence and breakup models.

should not exceed the energy density of the parent bubble.

density increase, is proposed here.

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 Vi ¼ 2Vi�1. The population balance equation is expressed by Eq. (1)

$$\frac{\partial n\_i}{\partial t} + \nabla \cdot \left(\overrightarrow{v}\_i \cdot n\_i\right) = \mathbf{S}\_i \tag{1}$$

where ni is the number density for i-th group, vi ! is the mass average velocity vector and Si is the source term.

The source term for the i-th group, Si, 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{split} \mathbf{S}\_{i} &= \mathbf{B}\_{\text{coupleence},i} - D\_{\text{coupleence},i} + B\_{\text{break}p,i} - D\_{\text{break}p,i} \\ &= \sum\_{d\_{\text{eq}}=d\_{\text{eq}}=d\_{\text{eq}}}^{d\_{\text{eq}}/2} \, \Omega\_{\mathbb{C}} \big( d\_{\text{eq},i} : d\_{\text{eq},i} - d\_{\text{eq},i} \big) - \sum\_{d\_{\text{eq}j}}^{d\_{\text{eq},\text{max}}-d\_{\text{eq},i}} \, \Omega\_{\mathbb{C}} \big( d\_{\text{eq},j} : d\_{\text{eq},i} \big) + \sum\_{d\_{j}=d\_{i}}^{d\_{\text{max}}} \, \Omega\_{\mathbb{B}} \big( d\_{\text{eq},j} : d\_{\text{eq},i} \big) - \, \Omega\_{\mathbb{B}}(d\_{\text{eq},j}) \end{split} \tag{2}$$

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

$$
\alpha\_{\mathcal{S}} f\_i = n\_i V\_i \tag{3}
$$

where f <sup>i</sup> is the i-th class fraction of the total volume fraction and Vi 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

collision tube, a nominal diameter, dV, that approximately represents the size of the projected

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

where c and a are the length of the short axis and long axis, respectively. The new imaginary

The breakup rate for one individual parent bubble that forms into two daughter bubbles can

<sup>B</sup> is the collision probability density, which can be estimated from Luo and Swendsen

ni

Here, ξ =λ /deq,i is the non-dimensional size of eddies that may contribute to the breakage of bubbles with size di. The breakage probability function pB used by Luo and Svendsen [3] is

pB <sup>¼</sup> exp � es

where e is the mean turbulent kinetic energy for eddies of size λ and es is the increase in surface energy of bubbles after breakage. The mean turbulent kinetic energy can be determined by Eq. (8)

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

e

ðd λmin ωT

� �<sup>1</sup>=<sup>3</sup>

Ω<sup>B</sup> ¼

� � εdeq,i

c ≤ dV ≤ a (4)

http://dx.doi.org/10.5772/intechopen.76448

69

<sup>B</sup>pB dλ (5)

� � (7)

<sup>ξ</sup><sup>11</sup>=<sup>3</sup> (6)

dV,s=deq,i <sup>þ</sup> <sup>ξ</sup> � � dV,l=deq,i <sup>þ</sup> <sup>ξ</sup> � � d2 eq,i

area of the bubble, is defined in a bounded range given by expression (4),

collision tube is presented in Figure 2.

be calculated using Eq. (5),

[3], as defined by Eq. (6)

given by Eq. (7),

ωT

<sup>B</sup>ð Þ¼ ξ 0:923 1 � α<sup>g</sup>

where ω<sup>T</sup>

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

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

$$c \le d\_V \le a \tag{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),

$$
\Omega\_B = \int\_{\lambda\_{\rm min}}^d \omega\_B^T p\_B \, d\lambda \tag{5}
$$

where ω<sup>T</sup> <sup>B</sup> is the collision probability density, which can be estimated from Luo and Swendsen [3], as defined by Eq. (6)

$$
\omega\_B^T(\xi) = 0.923 \left( 1 - \alpha\_\xi \right) \left( \varepsilon d\_{eq,i} \right)^{1/3} n\_i \frac{\left( d\_{V,s}/d\_{eq,i} + \xi \right) \left( d\_{V,l}/d\_{eq,i} + \xi \right)}{d\_{eq,i}^2 \xi^{11/3}} \tag{6}
$$

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

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

Figure 1. Time sequence photos of the breakup of a rising bubble in a 150-mm diameter cylindrical bubble column

68 Heat and Mass Transfer - Advances in Modelling and Experimental Study for Industrial Applications

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

energy of bubbles and interfacial area.

(Ug = 0.02 m/s; total duration: 0.0366 s).

$$p\_{\beta} = \exp\left(-\frac{e\_s}{\overline{e}}\right) \tag{7}$$

where e is the mean turbulent kinetic energy for eddies of size λ and es is the increase in surface energy of bubbles after breakage. The mean turbulent kinetic energy can be determined by Eq. (8)

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

$$\frac{\overline{\varepsilon} = \rho\_l \frac{\pi}{6} \lambda^3 \overline{u}\_{\lambda}^2}{2 = \frac{\pi \theta}{12} \rho\_l \left(\varepsilon d\_{eq,i}\right)^{2/3} d\_{eq,i}^3 \xi^{11/3}} \tag{8}$$

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

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

1

ffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>E</sup><sup>2</sup> � <sup>1</sup>

where the aspect ratio E can be expressed using the empirical correlation described by Wellek

g r<sup>l</sup> � r<sup>g</sup> � �d<sup>2</sup>

The size boundary to divide between ellipsoidal and spherical-cap bubbles represented by dC

dC <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>40</sup>σ=<sup>g</sup>ð Þ <sup>r</sup>L�r<sup>g</sup>

where dc is determined to be 17.3 mm for the air-water system. For a single spherical-cap bubble, the wake angle θ<sup>W</sup> is assumed to be 50o 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)

The curved surface area for the front edge can be calculated from the following relationship

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

eq=6

<sup>2</sup> � ð Þ <sup>1</sup> � cosθ<sup>W</sup>

3

<sup>p</sup> ln 2E<sup>2</sup> � <sup>1</sup> <sup>þ</sup> <sup>2</sup><sup>E</sup>

eq

ffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>E</sup><sup>2</sup> � <sup>1</sup>

http://dx.doi.org/10.5772/intechopen.76448

71

<sup>σ</sup> (13)

<sup>q</sup> (14)

ð Þ 1 � cosθ<sup>W</sup> (16)

<sup>=</sup><sup>3</sup> (15)

� � p ! (11)

<sup>E</sup> <sup>¼</sup> <sup>a</sup>=<sup>b</sup> <sup>¼</sup> <sup>1</sup> <sup>þ</sup> <sup>0</sup>:163Eo<sup>0</sup>:<sup>757</sup> (12)

2E

Eo ¼

<sup>S</sup> <sup>¼</sup> <sup>π</sup>d<sup>3</sup>

ð Þ 1 � cosθ<sup>W</sup>

SCap <sup>¼</sup> <sup>2</sup>πR<sup>2</sup>

area can be calculated using Eq. (11),

et al. [17], which is given by Eq. (12)

is estimated using Eq. (14),

can be formulated as follows:

given by Eq. (16):

Sellipsoid <sup>¼</sup> <sup>π</sup>

2 d2

Here, Eo is the Eötvös number as defined by Eq. (13)

R3

eqE<sup>1</sup>=<sup>3</sup> <sup>1</sup> <sup>þ</sup>

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,i 3 -deq, j 3 ) 1/3, the increase in surface energy can be estimated using Eq. (9),

$$\left(\varepsilon\_s \left(d\_{\text{eq},i}, d\_{\text{eq},j}\right)\right) = \sigma \cdot \pi d\_{\text{eq},i}^2 \left[f\_V^{2/3} + \left(1 - f\_V\right)^{2/3} - 1\right] \tag{9}$$

where the breakage volume fraction is given by <sup>f</sup> <sup>V</sup> <sup>¼</sup> <sup>d</sup><sup>3</sup> eq,j =d<sup>3</sup> eq,i . 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)

$$\mathbf{e}\_s = \sigma \cdot \left(\mathbf{S}\_{\mathbf{\hat{\cdot}},1} + \mathbf{S}\_{\mathbf{\hat{\cdot}},2} - \mathbf{S}\_{\mathbf{\hat{\cdot}}}\right) \tag{10}$$

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

$$S\_{\rm ellipsoid} = \frac{\pi}{2} d\_{eq}^2 E^{1/3} \left( 1 + \frac{1}{2E\sqrt{E^2 - 1}} \ln \left( 2E^2 - 1 + 2E\sqrt{E^2 - 1} \right) \right) \tag{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)

$$E = a/b = 1 + 0.163Eo^{0.757} \tag{12}$$

Here, Eo is the Eötvös number as defined by Eq. (13)

e ¼ r<sup>l</sup> π <sup>6</sup> <sup>λ</sup><sup>3</sup> u2 λ

<sup>12</sup> r<sup>l</sup> εdeq,i � �<sup>2</sup>=<sup>3</sup>

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

> eq,i f 2=3

shapes of bubbles are now taken into account, Eq. (9) can be rewritten in a general form in

es ¼ σ � Sj, <sup>1</sup> þ Sj,<sup>2</sup> � Si

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

d3 eq,i

<sup>V</sup> þ 1 � f <sup>V</sup>

eq,j =d<sup>3</sup> eq,i

� �<sup>2</sup>=<sup>3</sup> � <sup>1</sup> h i

� � (10)

<sup>ξ</sup><sup>11</sup>=<sup>3</sup> (8)

. Since the effects of different

3 -deq,

(9)

<sup>2</sup> <sup>¼</sup> πβ

70 Heat and Mass Transfer - Advances in Modelling and Experimental Study for Industrial Applications

1/3, the increase in surface energy can be estimated using Eq. (9),

� � <sup>¼</sup> <sup>σ</sup> � <sup>π</sup>d<sup>2</sup>

es deq,i; deq,j

where the breakage volume fraction is given by <sup>f</sup> <sup>V</sup> <sup>¼</sup> <sup>d</sup><sup>3</sup>

potential breakage scenarios are depicted in Figure 3.

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

terms of the surface area of the bubbles, S, as defined by Eq. (10)

j 3 )

$$E\sigma = \frac{\lg\left(\rho\_l - \rho\_\S\right)d\_{eq}^2}{\sigma} \tag{13}$$

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

$$\mathrm{d}\_{\mathbb{C}} = \sqrt{^{40\circ}|\_{\mathcal{S}} \, \_{\mathcal{S}} \, \_{\mathbb{C}}} \, \tag{14}$$

where dc is determined to be 17.3 mm for the air-water system. For a single spherical-cap bubble, the wake angle θ<sup>W</sup> is assumed to be 50o 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\_{eq}^3 / 6}{\left(1 - \cos \theta\_W\right)^2 - \left(1 - \cos \theta\_W\right)^3 / 3} \tag{15}$$

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

$$S\_{\rm Cap} = 2\pi R^2 (1 - \cos \theta\_W) \tag{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),

$$\varepsilon\_{P} = \frac{\sigma}{\min\{R\_{\mathbb{C},j}, R\_{\mathbb{C},k}\}} \cdot \frac{\pi\left(\min\{d\_{\text{eq},j}, d\_{\text{eq},k}\}\right)^3}{6} \tag{17}$$

2.2. Governing Eqs

ous phase and air as the dispersed phase.

∂ ∂t

tion rate ε<sup>l</sup> are computed using Eqs. (22) and (23),

∂ <sup>∂</sup><sup>t</sup> <sup>r</sup><sup>i</sup> αiki <sup>þ</sup> <sup>∇</sup><sup>∙</sup> <sup>r</sup><sup>i</sup>

<sup>þ</sup>∇<sup>∙</sup> <sup>r</sup><sup>i</sup>

2.3. Interphase momentum transfer

using Eq. (25),

∂ <sup>∂</sup><sup>t</sup> <sup>r</sup>iαiε<sup>i</sup>

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

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

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

where rk, αk, uk, τ<sup>k</sup> and F<sup>k</sup> represent the density, volume fraction, velocity vector, viscous stress tensor and the interphase momentum exchange term for the k (liquid or gas) phase, respec-

A modified k˜ε turbulence model with the consideration of bubble-induced turbulence by Sato and Sekoguchi [22] is used for turbulence closure. The turbulent kinetic energy kl and dissipa-

> þ μeff , <sup>i</sup> σε <sup>∇</sup>ε<sup>i</sup>

where Gk,l is the production of turbulent kinetic energy and μt,l is the turbulent viscosity. In this work, the standard k˜ε model constants used are C<sup>μ</sup> = 0.09, C<sup>1</sup><sup>ε</sup> = 1.44, C<sup>2</sup><sup>ε</sup> = 1.92, σ<sup>k</sup> = 1.0, σε = 1.3. The effective viscosity is composed of the contributions of turbulent viscosity and an extra

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

> α<sup>g</sup> u<sup>g</sup> � u<sup>l</sup>

ð Þþ rkα<sup>k</sup> ∇∙ð Þ¼ rkαku<sup>k</sup> 0 (20)

þ α<sup>i</sup> Gk,i � r<sup>i</sup>

C<sup>1</sup>εGk,i � C<sup>2</sup>εriε<sup>i</sup>

εi (22)

http://dx.doi.org/10.5772/intechopen.76448

73

(23)

(24)

(25)

ð Þþ rkαku<sup>k</sup> ∇∙ð Þ¼� rkαkuku<sup>k</sup> αk∇pþ∇∙τ<sup>k</sup> þ αkrkg þ F<sup>k</sup> (21)

μeff ,i σk <sup>∇</sup>k<sup>i</sup> 

> Cμ,BITαgdb u<sup>g</sup> � u<sup>l</sup>

> > u<sup>g</sup> � u<sup>l</sup>

þ α<sup>i</sup> εi ki

∂ ∂t

tively. The sum of the volume fractions for both phases is equal to 1.

αikiu<sup>i</sup>

μeff ,l ¼ r<sup>l</sup>

The Sato coefficient Cμ,BIT ¼ 0:6 is adopted according to the study [22].

<sup>F</sup><sup>D</sup> <sup>¼</sup> <sup>3</sup> 4 CD deq rl

αiεiu<sup>i</sup> <sup>¼</sup>∇<sup>∙</sup> <sup>α</sup><sup>i</sup> <sup>μ</sup><sup>i</sup>

<sup>¼</sup> <sup>∇</sup><sup>∙</sup> <sup>α</sup><sup>i</sup> <sup>μ</sup><sup>i</sup> <sup>þ</sup>

term considering the effect of bubble-induced turbulence and is defined by Eq. (24)

C<sup>μ</sup> k2 l εl þ r<sup>l</sup>

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

$$6\,\text{cS}\_{i}\Big/\mathfrak{m}\_{\text{eq},i}^{3} \geq 6\,\sigma\cdot\max(S\_{j}/\pi d\_{\text{eq},j}^{3}, S\_{k}/\pi d\_{\text{eq},k}^{3}) - 6\,\text{cS}\_{i}/\pi d\_{\text{eq},i}^{3}\tag{18}$$

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

Ω<sup>B</sup> ¼ 0:923 1 � α<sup>g</sup> � �ni ε=d<sup>2</sup> eq,i � �<sup>1</sup>=<sup>3</sup> � Ð 1 ξmin dV,s=deq,i <sup>þ</sup> <sup>ξ</sup> � � dV,l=deq,i <sup>þ</sup> <sup>ξ</sup> � � <sup>ξ</sup><sup>11</sup>=<sup>3</sup> exp � <sup>12</sup><sup>σ</sup> Sj <sup>þ</sup> Sk � Si � � πβr<sup>l</sup> ε<sup>2</sup>=<sup>3</sup>ξ<sup>11</sup>=<sup>3</sup> d 11=3 eq,i 0 @ 1 Adξ, when <sup>6</sup><sup>σ</sup> Sj <sup>þ</sup> Sk � Si � � πd<sup>3</sup> eq,i <sup>≥</sup> <sup>σ</sup> min Rcj; Rck � � 0:923 1 � α<sup>g</sup> � �ni ε=d<sup>2</sup> eq,i � �<sup>1</sup>=<sup>3</sup> � Ð 1 ξmin dV,s=deq,i <sup>þ</sup> <sup>ξ</sup> � � dV,l=deq,i <sup>þ</sup> <sup>ξ</sup> � � <sup>ξ</sup><sup>11</sup>=<sup>3</sup> exp � <sup>2</sup><sup>σ</sup> min deq,j; deq, <sup>k</sup> � � � � <sup>3</sup> min Rcj; Rck � �βr<sup>l</sup> ε<sup>2</sup>=<sup>3</sup>ξ<sup>11</sup>=<sup>3</sup> d 11=3 eq,i 0 @ 1 Adξ, when <sup>6</sup><sup>σ</sup> Sj <sup>þ</sup> Sk � Si � � πd<sup>3</sup> eq,i <sup>&</sup>lt; <sup>σ</sup> min Rcj;Rck � � 8 >>>>>>>>>>>>>>>>>>>< >>>>>>>>>>>>>>>>>>>: (19)

where ξmin is the minimum breakage volume fraction that is able to satisfy the energy density criterion shown in Eq. (18).

#### 2.2. Governing Eqs

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

72 Heat and Mass Transfer - Advances in Modelling and Experimental Study for Industrial Applications

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

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

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

when <sup>6</sup><sup>σ</sup> Sj <sup>þ</sup> Sk � Si � � πd<sup>3</sup> eq,i

when <sup>6</sup><sup>σ</sup> Sj <sup>þ</sup> Sk � Si � � πd<sup>3</sup> eq,i

where ξmin is the minimum breakage volume fraction that is able to satisfy the energy density

<sup>ξ</sup><sup>11</sup>=<sup>3</sup> exp � <sup>12</sup><sup>σ</sup> Sj <sup>þ</sup> Sk � Si

0 @

πβr<sup>l</sup>

<sup>≥</sup> <sup>σ</sup> min Rcj; Rck � �

<sup>&</sup>lt; <sup>σ</sup> min Rcj;Rck � �

<sup>ξ</sup><sup>11</sup>=<sup>3</sup> exp � <sup>2</sup><sup>σ</sup> min deq,j; deq, <sup>k</sup>

0 @ � �

1 Adξ,

� � � � <sup>3</sup>

ε<sup>2</sup>=<sup>3</sup>ξ<sup>11</sup>=<sup>3</sup> d 11=3 eq,i

1 Adξ,

ε<sup>2</sup>=<sup>3</sup>ξ<sup>11</sup>=<sup>3</sup> d 11=3 eq,i

min Rcj; Rck � �βr<sup>l</sup>

dV,s=deq,i <sup>þ</sup> <sup>ξ</sup> � � dV,l=deq,i <sup>þ</sup> <sup>ξ</sup> � �

dV,s=deq,i <sup>þ</sup> <sup>ξ</sup> � � dV,l=deq,i <sup>þ</sup> <sup>ξ</sup> � �

π min deq,j; deq, <sup>k</sup> � � � � <sup>3</sup>

<sup>6</sup> (17)

ð18Þ

(19)

increase in surface energy is further illustrated in Figure 7.

eP <sup>¼</sup> <sup>σ</sup>

min RC,j; RC, <sup>k</sup> � � �

requirement can be expressed using Eq. (17),

with the variation of bubble shapes

� �ni ε=d<sup>2</sup>

� �ni ε=d<sup>2</sup>

eq,i � �<sup>1</sup>=<sup>3</sup>

eq,i � �<sup>1</sup>=<sup>3</sup> � Ð 1 ξmin

� Ð 1 ξmin

0:923 1 � α<sup>g</sup>

0:923 1 � α<sup>g</sup>

criterion shown in Eq. (18).

in Eq. (19),

8

>>>>>>>>>>>>>>>>>>><

>>>>>>>>>>>>>>>>>>>:

Ω<sup>B</sup> ¼

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{\partial}{\partial t}(\rho\_k \alpha\_k) + \nabla \cdot (\rho\_k \alpha\_k \mathfrak{u}\_k) = \mathbf{0} \tag{20}$$

$$\frac{\partial}{\partial t}(\rho\_k \alpha\_k \mathfrak{u}\_k) + \nabla \cdot (\rho\_k \alpha\_k \mathfrak{u}\_k \mathfrak{u}\_k) = -\alpha\_k \nabla p + \nabla \cdot \mathfrak{T}\_k + \alpha\_k \rho\_k \mathfrak{g} + \mathcal{F}\_k \tag{21}$$

where rk, αk, uk, τ<sup>k</sup> and F<sup>k</sup> 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˜ε turbulence model with the consideration of bubble-induced turbulence by Sato and Sekoguchi [22] is used for turbulence closure. The turbulent kinetic energy kl and dissipation rate ε<sup>l</sup> are computed using Eqs. (22) and (23),

$$\frac{\partial}{\partial t} \left( \rho\_i a\_i k\_i \right) + \nabla \cdot \left( \rho\_i a\_i k\_i \mathfrak{u}\_i \right) = \nabla \cdot \left[ a\_i \left( \mu\_i + \frac{\mu\_{\text{eff},i}}{\sigma\_k} \right) \nabla k\_i \right] + a\_i \left( \mathcal{G}\_{k,i} - \rho\_i \varepsilon\_i \right) \tag{22}$$

$$\frac{\partial}{\partial t} \left( \rho\_i \alpha\_i \varepsilon\_i \right) + \nabla \cdot \left( \rho\_i \alpha\_i \varepsilon\_i \mu\_i \right) = \nabla \cdot \left[ \alpha\_i \left( \mu\_i + \frac{\mu\_{\varepsilon \mathcal{G},i}}{\sigma\_\varepsilon} \right) \nabla \varepsilon\_i \right] + \alpha\_i \frac{\varepsilon\_i}{k\_i} \left( \mathbb{C}\_{1\varepsilon} \mathbb{G}\_{k,i} - \mathbb{C}\_{2\varepsilon} \rho\_i \varepsilon\_i \right) \tag{23}$$

where Gk,l is the production of turbulent kinetic energy and μt,l is the turbulent viscosity. In this work, the standard k˜ε model constants used are C<sup>μ</sup> = 0.09, C<sup>1</sup><sup>ε</sup> = 1.44, C<sup>2</sup><sup>ε</sup> = 1.92, σ<sup>k</sup> = 1.0, σε = 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)

$$\left| \boldsymbol{\mu}\_{\rm eff,l} = \rho\_l \mathbf{C}\_{\mu} \frac{k\_l^2}{\varepsilon\_l} + \rho\_l \mathbf{C}\_{\mu, BIT} \alpha\_{\rm g} d\_b \left| \boldsymbol{\mu}\_{\rm g} - \boldsymbol{\mu}\_l \right| \tag{24}$$

The Sato coefficient Cμ,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),

$$F\_D = \frac{3}{4} \frac{C\_D}{d\_{cq}} \rho\_l \alpha\_\wp \left| \boldsymbol{\mu}\_\wp - \boldsymbol{\mu}\_l \right| \left( \boldsymbol{\mu}\_\wp - \boldsymbol{\mu}\_l \right) \tag{25}$$

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

$$\mathbf{C}\_{D,sphere} = \begin{cases} 24/\text{Re}\_b & \text{Re}\_b < 0.01\\ 24\left(1 + 0.15 \text{Re}\_b^{0.687}\right) / \text{Re}\_b & \text{Re}\_b \ge 0.01 \end{cases} \tag{26}$$

$$\mathbb{C}\_{D,cap} = \frac{8}{3} \tag{27}$$

CL ¼

given, respectively, by Eqs. (34) and (35)

force is still calculated using Eq. (36),

2.4. Numerical modelling

Table 1. Details of the experimental setup.

Table 1.

8 >>><

>>>:

where f Eo<sup>0</sup> � � <sup>¼</sup> <sup>0</sup>:00105Eo0<sup>3</sup> � <sup>0</sup>:0159Eo0<sup>2</sup> � <sup>0</sup>:0204Eo<sup>0</sup>

min 0:288tanh 0ð Þ :121Re<sup>b</sup> ; f Eo<sup>0</sup> h i � � Eo<sup>0</sup>

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

f Eo<sup>0</sup> � � <sup>4</sup> <sup>&</sup>lt; Eo<sup>0</sup>

�0:<sup>29</sup> Eo<sup>0</sup>

ber based on the maximum horizontal dimension of the deformable bubble, dh, as defined and

g r<sup>l</sup> � r<sup>g</sup> � �d<sup>2</sup>

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

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

The mesh setup is illustrated in Figure 4. Grid 2 consists of 20ð Þ� r 40ð Þ� θ 100ð Þz nodes in the radial, circumferential and axial directions, respectively. The grid independence was tested in a coarser Grid 1 of 16ð Þ� r 32ð Þ� θ 80ð Þz nodes and a refined Grid 3 of 26ð Þ� r 48ð Þ� θ 126ð Þz 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.

Case 1 0.15 0.8 0.0382 0.65 Case 2 0.1 2 0.0606 0.9

Diameter (m) Height (m) Superficial velocity (m/s) Static liquid height (m)

dul dt � <sup>d</sup>ug

h

Eo<sup>0</sup> ¼

FVM ¼ CVMrlα<sup>g</sup>

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

≤ 4

< 10

http://dx.doi.org/10.5772/intechopen.76448

(33)

75

> 10

þ 0:474. Eo' is the modified Eötvös num-

<sup>σ</sup> (34)

dt � � (36)

dh <sup>¼</sup> <sup>d</sup> <sup>1</sup> <sup>þ</sup> <sup>0</sup>:163Eo<sup>0</sup>:<sup>757</sup> � �<sup>1</sup>=<sup>3</sup> (35)

$$\mathbf{C}\_{D,ellipse} = \frac{4}{3} \frac{g d\_{eq}}{U\_t^2} \frac{\left(\rho\_l - \rho\_\mathcal{S}\right)}{\rho\_l} \tag{28}$$

where Reb is the bubble Reynolds number given by Re<sup>b</sup> <sup>¼</sup> <sup>r</sup>lj j ug�u<sup>l</sup> deq μl and Ut is the terminal velocity calculated using Eq. (29),

$$\mathcal{U}I\_t = \frac{\mu\_l}{\rho\_l d} M o^{-0.149} (I - 0.857) \tag{29}$$

Here, Mo is the Morton number defined by Mo <sup>¼</sup> <sup>μ</sup><sup>4</sup> <sup>l</sup> <sup>g</sup>ð Þ <sup>r</sup>l�r<sup>g</sup> r2 <sup>l</sup> <sup>σ</sup><sup>3</sup> and <sup>J</sup> is determined by the piecewise function calculated using the empirical expression (30)

$$J = \begin{cases} 0.94H^{0.757} & 2 < H < 59.3\\ 3.42H^{0.441} & H > 59.3 \end{cases} \tag{30}$$

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

$$H = \frac{4}{3} E o M o^{-0.149} \left(\frac{\mu\_l}{\mu\_{r\!f}}\right)^{-0.14} \tag{31}$$

where Eo is the Eötvös number and μref ¼ 0:0009 kg=ð Þ m � s .

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

$$F\_{L\circ\!\!\! } = \mathbb{C}\_{L} \rho\_{l} \alpha\_{\mathcal{g}} \left(\mathfrak{u}\_{\mathcal{g}} - \mathfrak{u}\_{l}\right) \times \left(\nabla \times \mathfrak{u}\_{l}\right) \tag{32}$$

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

Modelling of Bubbly Flow in Bubble Column Reactors with an Improved Breakup Kernel Accounting for Bubble… http://dx.doi.org/10.5772/intechopen.76448 75

$$\mathbf{C}\_{L} = \begin{cases} \min\left[0.288 \tanh(0.121 \text{Re}\_{b}), f\left(Eo^{'}\right)\right] & \text{Eo}^{'} \le 4\\ f\left(Eo^{'}\right) & 4 < Eo^{'} < 10\\ -0.29 & Eo^{'} > 10 \end{cases} \tag{33}$$

where f Eo<sup>0</sup> � � <sup>¼</sup> <sup>0</sup>:00105Eo0<sup>3</sup> � <sup>0</sup>:0159Eo0<sup>2</sup> � <sup>0</sup>:0204Eo<sup>0</sup> þ 0:474. 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)

$$E\sigma^{'} = \frac{g\left(\rho\_l - \rho\_{\mathcal{S}}\right)d\_h^2}{\sigma} \tag{34}$$

$$d\_{\rm h} = d \left( 1 + 0.163 Eo^{0.757} \right)^{1/3} \tag{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),

$$F\_{VM} = C\_{VM} \rho\_l \alpha\_\S \left(\frac{d\mathfrak{u}\_l}{dt} - \frac{d\mathfrak{u}\_\mathcal{g}}{dt}\right) \tag{36}$$

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

#### 2.4. Numerical modelling

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

> CD,sphere <sup>¼</sup> <sup>24</sup>=Re<sup>b</sup> Re<sup>b</sup> <sup>&</sup>lt; <sup>0</sup>:<sup>01</sup> 24 1 <sup>þ</sup> <sup>0</sup>:15Re<sup>0</sup>:<sup>687</sup>

> > CD, cap <sup>¼</sup> <sup>8</sup>

3 gdeq U2 t

<sup>J</sup> <sup>¼</sup> <sup>0</sup>:94H<sup>0</sup>:<sup>757</sup> <sup>2</sup> <sup>&</sup>lt; <sup>H</sup> <sup>&</sup>lt; <sup>59</sup>:<sup>3</sup> 3:42H<sup>0</sup>:<sup>441</sup> H > 59:3

EoMo�0:<sup>149</sup> <sup>μ</sup><sup>l</sup>

The lift force acting perpendicular to the direction of relative motion of the two phases can be

where CL is the lift coefficient and is estimated by the Tomiyama lift force correlation [12], as

FLift ¼ CLrlα<sup>g</sup> u<sup>g</sup> � u<sup>l</sup>

μref

!�0:<sup>14</sup>

CD,ellipse <sup>¼</sup> <sup>4</sup>

where Reb is the bubble Reynolds number given by Re<sup>b</sup> <sup>¼</sup> <sup>r</sup>lj j ug�u<sup>l</sup> deq

(

<sup>H</sup> <sup>¼</sup> <sup>4</sup> 3

where Eo is the Eötvös number and μref ¼ 0:0009 kg=ð Þ m � s .

described by the following empirical relation (33),

Ut <sup>¼</sup> <sup>μ</sup><sup>l</sup> rld

�

74 Heat and Mass Transfer - Advances in Modelling and Experimental Study for Industrial Applications

b

� �=Re<sup>b</sup> Re<sup>b</sup> ≥ 0:01

r<sup>l</sup> � r<sup>g</sup> � �

rl

<sup>l</sup> <sup>g</sup>ð Þ <sup>r</sup>l�r<sup>g</sup> r2

<sup>3</sup> (27)

μl

Mo�0:<sup>149</sup>ð Þ <sup>J</sup> � <sup>0</sup>:<sup>857</sup> (29)

� � � ð Þ <sup>∇</sup> � <sup>u</sup><sup>l</sup> (32)

<sup>l</sup> <sup>σ</sup><sup>3</sup> and <sup>J</sup> is determined by the piecewise

(26)

(28)

(30)

(31)

and Ut is the terminal

Eqs. (26)–(28),

velocity calculated using Eq. (29),

Here, Mo is the Morton number defined by Mo <sup>¼</sup> <sup>μ</sup><sup>4</sup>

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

calculated by using Eq. (32)

function calculated using the empirical expression (30)

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.

The mesh setup is illustrated in Figure 4. Grid 2 consists of 20ð Þ� r 40ð Þ� θ 100ð Þz nodes in the radial, circumferential and axial directions, respectively. The grid independence was tested in a coarser Grid 1 of 16ð Þ� r 32ð Þ� θ 80ð Þz nodes and a refined Grid 3 of 26ð Þ� r 48ð Þ� θ 126ð Þz 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.


Table 1. Details of the experimental setup.

76 Heat and Mass Transfer - Advances in Modelling and Experimental Study for Industrial Applications

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

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

http://dx.doi.org/10.5772/intechopen.76448

77

conditions are applied for both the liquid and gas phases at the bubble column wall.

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

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

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.

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

3. Results and discussion

required for bubble breakage

spherical-capped bubbles.

bubbles (deq,i = 32 mm).

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