*Two-Fluid Model*

The two-fluid model for the flow of the gas bubbles in the Newtonian fluid with the Euler–Euler model was employed. Both phases were considered as a continuum and could be defined with the help of differential equations. A thin surface separated both the fluids and at the interphase, and consequently a jump continuity condition could be employed [44]. The bubbles were distributed homogenously in the flow and considered to be of the same spherical shape. It was assumed that no mass was transferred between the two phases and isothermal conditions were maintained. The density of the liquid was constant, whereas for the gas-phase it depended on pressure, *p*. Bubble redispersion was neglected and there was low gas holdup. With these assumptions, continuity equations for both the liquid and bubble phase were of the form [45]

$$\nabla.(\mathcal{U}\_l) = 0,\tag{2}$$

$$
\frac{
\partial
\left(
\varepsilon\_{\mathcal{S}} \rho\_{\mathcal{S}}
\right)
}{
\partial t
} + \nabla \left(
\varepsilon\_{\mathcal{S}} \rho\_{\mathcal{S}} \mathcal{U}\_{\mathcal{S}}
\right) = 0,
\tag{3}
$$

where Ul and Ub were the liquid and bubble velocities, respectively, whereas ε<sup>g</sup> and ρ<sup>g</sup> represented the void space and density of the gaseous phase with the closure relation constraint ε<sup>l</sup> + ε<sup>g</sup> = 1. For low gas <sup>ε</sup><sup>l</sup> <sup>≈</sup> 1. Using ideal gas laws, <sup>ρ</sup><sup>g</sup> <sup>=</sup> <sup>p</sup> RT0 , ρ<sup>l</sup> ≈ constant, and because isothermal conditions were taken into account T0 remained constant.

The momentum transfer for both the liquid and gaseous phase [46] were defined as

$$
\rho\_l \left( \frac{\partial \mathcal{U}\_l}{\partial t} + \mathcal{U}\_l \frac{\partial \mathcal{U}\_l}{\partial X} + \mathcal{V}\_l \frac{\partial \mathcal{U}\_l}{\partial Y} \right) = -\frac{\partial P}{\partial X} + \mu \left( \frac{\partial^2 \mathcal{U}\_l}{\partial X^2} + \frac{\partial^2 \mathcal{U}\_l}{\partial Y^2} \right) + \rho\_l \mathcal{g} - \varepsilon\_{\mathcal{S}} \rho\_l \mathcal{g}\_{\mathcal{V}} \tag{4}
$$

$$
\rho \eta \left( \frac{\partial V\_I}{\partial t} + \mathcal{U}\_I \frac{\partial V\_I}{\partial X} + V\_I \frac{\partial V\_I}{\partial Y} \right) = -\frac{\partial P}{\partial Y} + \mu \left( \frac{\partial^2 V\_I}{\partial X^2} + \frac{\partial^2 V\_I}{\partial Y^2} \right), \tag{5}
$$

$$\varepsilon\_{\mathcal{S}}\rho\_{\mathcal{S}}\Big(\frac{\partial \mathcal{U}\_{\mathcal{S}}}{\partial t} + \mathcal{U}\_{\mathcal{S}}\frac{\partial \mathcal{U}\_{\mathcal{S}}}{\partial X} + \mathcal{V}\_{\mathcal{S}}\frac{\partial \mathcal{U}\_{\mathcal{S}}}{\partial Y}\Big) = -\varepsilon\_{\mathcal{S}}\frac{\partial \mathcal{P}}{\partial X} + \mu\varepsilon\_{\mathcal{S}}\Big(\frac{\partial^{2} \mathcal{U}\_{\mathcal{S}}}{\partial X^{2}} + \frac{\partial^{2} \mathcal{U}\_{\mathcal{S}}}{\partial Y^{2}}\Big) + \rho\_{\mathcal{S}}g\varepsilon\_{\mathcal{S}} - F\_{\text{int},\text{x}}\tag{6}$$

$$
\varepsilon\_{\mathcal{S}} \rho\_{\mathcal{S}} \left( \frac{\partial V\_{\mathcal{S}}}{\partial t} + \mathcal{U}\_{\mathcal{S}} \frac{\partial V\_{\mathcal{S}}}{\partial X} + V\_{\mathcal{S}} \frac{\partial V\_{\mathcal{S}}}{\partial Y} \right) = -\varepsilon\_{\mathcal{S}} \frac{\partial P}{\partial Y} + \mu \varepsilon\_{\mathcal{S}} \left( \frac{\partial^2 V\_{\mathcal{S}}}{\partial X^2} + \frac{\partial^2 V\_{\mathcal{S}}}{\partial Y^2} \right) - F\_{\text{int},y} \tag{7}
$$

where Fint is the interaction force between the liquid and gaseous phases. These forces were able to be divided into three parts—i) force due to drag, which is incorporated in the uniform flow, ii) added mass force due to the acceleration of the bubble, and iii) lift force for the spherical bubble [47]. Mathematically this could be written as

$$F\_{\rm int} = F\_d + F\_{am} + F\_{l\prime} \tag{8}$$

$$F\_d = -\mathbb{C}\_w V\_b (\mathcal{U}\_b - \mathcal{U}\_l)\_\prime \tag{9}$$

$$F\_{am} = -\mathcal{C}\_{am} V\_{b} \rho\_{l} \left(\frac{d\mathcal{U}\_{slip}}{dt}\right) \tag{10}$$

$$F\_I = -\mathbb{C}\_I V\_b \rho\_l (\mathcal{U}\_b - \mathcal{U}\_l) \times (\nabla \times \mathcal{U}\_l),\tag{11}$$

in which Cw was a constant involving bubble diameter and drag force, Vb was the volume of the bubble, Cam was a constant related to the volume fraction of liquid accelerated with the bubble, Uslip was slip velocity defined as *Uslip* = *Ub* − *Ul*, and *Cl* was the lift coefficient. Finally, *Ub* was the speed of a single bubble and it was related to gas velocity by *Ug* = *Ub* + *Udri f t*, where *Udri f t* is the average drift velocity, which was proportional to the ratio of the special change in gas void fraction to total void fraction of gas.

To transform from a fixed frame to wave frame, Lorentz transformations were employed.

$$X = \mathbf{x} \ast \mathbf{-} \text{ct}, \; \mathbf{Y} = \mathbf{y} \ast, \; \mathcal{U} = \mathbf{u} \ast - \text{c}, \; \mathcal{V} = \mathbf{v} \ast, \mathbf{p} \ast \begin{pmatrix} \mathbf{x} \ast \; \mathbf{y} \ast \end{pmatrix} = \mathbf{P}(\mathbf{X}, \mathbf{Y}, \mathbf{t}). \tag{12}$$

Equations (2)-(11) became

$$\nabla.(\mu\_l\*) = 0,\tag{13}$$

$$\varepsilon \frac{\partial \left( \varepsilon\_{\mathcal{S}} \rho\_{\mathcal{S}} \right)}{\partial \mathbf{x}^\*} + \nabla \left( \varepsilon\_{\mathcal{S}} \rho\_{\mathcal{S}} \mu\_{\mathcal{S}} \mathbf{\hat{z}} \right) = 0,\tag{14}$$

$$
\rho\_l \left( \mu\_l \ast \frac{\partial \boldsymbol{u}\_{l^\bullet}}{\partial \mathbf{x}^\bullet} + \boldsymbol{v}\_l \ast \frac{\partial \boldsymbol{u}\_{l^\bullet}}{\partial \mathbf{y}^\bullet} \right) = -\frac{\partial p}{\partial \mathbf{x}^\bullet} + \mu \left( \frac{\partial^2 \boldsymbol{u}\_{l^\bullet}}{\partial \mathbf{x}^\bullet} + \frac{\partial^2 \boldsymbol{u}\_{l^\bullet}}{\partial \mathbf{y}^\bullet} \right) + (1 - \boldsymbol{\varepsilon}\_{\boldsymbol{\mathcal{g}}}) \rho\_l \mathbf{g}\_{l^\bullet} \tag{15}
$$

$$\rho\_l \rho\_l \left( u\_l \ast \frac{\partial \upsilon\_{l^\bullet}}{\partial \mathbf{x}^\bullet} + \upsilon\_l \ast \frac{\partial \upsilon\_{l^\bullet}}{\partial y^\bullet} \right) = -\frac{\partial p}{\partial y^\bullet} + \mu \left( \frac{\partial^2 \upsilon\_{l^\bullet}}{\partial \mathbf{x}^\bullet} + \frac{\partial^2 \upsilon\_{l^\bullet}}{\partial y^\bullet} \right) \tag{16}$$

$$\varepsilon\_{\mathcal{S}} \rho\_{\mathcal{S}} \Big( \mu\_{\mathcal{S}} \* \frac{\partial \mathfrak{u}\_{\mathcal{S}} \*}{\partial \mathbf{x}^{\bullet}} + \upsilon\_{\mathcal{S}} \* \frac{\partial \mathfrak{u}\_{\mathcal{S}} \*}{\partial \mathfrak{y}^{\bullet}} \Big) = -\varepsilon\_{\mathcal{S}} \frac{\partial p}{\partial \mathbf{x}^{\bullet}} + \mu \varepsilon\_{\mathcal{S}} \Big( \frac{\partial^{2} \mathfrak{u}\_{\mathcal{S}} \*}{\partial \mathbf{x}^{\bullet}} + \frac{\partial^{2} \mathfrak{u}\_{\mathcal{S}} \*}{\partial \mathfrak{y}^{\bullet}} \Big) + \rho\_{\mathcal{S}} g \varepsilon\_{\mathcal{S}} - f\_{\text{int}, \mathsf{x}} \tag{17}$$

$$
\varepsilon\_{\mathcal{S}} \rho\_{\mathcal{S}} \Big( u\_{\mathcal{S}} \ast \frac{\partial v\_{\mathcal{S}} \ast}{\partial \mathbf{x}^{\star}} + v\_{\mathcal{S}} \ast \frac{\partial v\_{\mathcal{S}} \ast}{\partial y^{\star}} \Big) = -\varepsilon\_{\mathcal{S}} \frac{\partial p^{\star}}{\partial y^{\star}} + \mu \varepsilon\_{\mathcal{S}} \Big( \frac{\partial^{2} v\_{\mathcal{S}} \ast}{\partial \mathbf{x}^{\star 2}} + \frac{\partial^{2} v\_{\mathcal{S}} \ast}{\partial y^{\star 2}} \Big) - f\_{\text{int}, y} \tag{18}
$$

$$f\_{\rm int} = f\_d + f\_{\rm am} + f\_{l\_2} \tag{19}$$

$$f\_d = -\mathbb{C}\_{\text{w}} V\_b (\mu\_b \ast -\mu\_l \ast)\_\prime \tag{20}$$

$$F\_{\rm am} = -\mathbb{C}\_{\rm anu} V\_{b} \rho\_{l} c \Big( \frac{\partial u\_{b} \*}{\partial \mathbf{x}} - \frac{\partial u\_{l} \*}{\partial \mathbf{x}} + (u\_{b} \* \cdot \nabla) u\_{b} \* - (u\_{l} \* \cdot \nabla) u\_{l} \* \Big) \tag{21}$$

$$F\_l = -\mathbb{C}\_l V\_b \rho\_l (\boldsymbol{\mu}\_b \ast -\boldsymbol{\mu}\_l \ast) \times (\nabla \times \boldsymbol{\mu}\_l \ast),\tag{22}$$

and by employing the dimensionless quantities were

 $\mathbf{x} = \frac{\mathbf{x}^\*}{\lambda}, \mathbf{y} = \frac{\mathbf{y}^\*}{\sigma}, \mathbf{u}\_{l, \emptyset} = \frac{\mathbf{u}\_{l, \emptyset} \mathbf{x}}{\varepsilon}, \mathbf{v}\_{l, \emptyset} = \frac{\mathbf{v}\_{l, \emptyset} \mathbf{x}}{\varepsilon}, \delta = \frac{\mathbf{r}\_{l, \emptyset} \mathbf{R}}{\varepsilon}, \mathbf{R} = \frac{\rho \mathbf{c} \mathbf{u}}{\mu}, p = \frac{\rho \mathbf{c} \mathbf{u}}{\mu \lambda \mathbf{c}} P$   $\mathbf{E} = \frac{a^2}{\sigma} \rho \mathbf{j}$ ,  $\mathbf{R}\_{\mathbf{b} \mathbf{c}} = \frac{a \mathbf{c}}{\mu} \rho\_l V\_{\mathbf{b} \mathbf{v}} A = \frac{\mu \mathbf{c}}{\sigma}$ .

By imposing Equation (23) on Equations (12)-(22) and considering Re and δ approach zero, the corresponding expressions resulted in

$$
\varepsilon\_{\mathcal{S}} = \frac{\varepsilon\_1}{\rho\_{\mathcal{S}} u\_{\mathcal{b}}},
\tag{24}
$$

$$\frac{\partial p}{\partial \mathbf{x}} = \frac{\partial^2 u\_l}{\partial y^2} + \left(1 - \varepsilon\_{\mathcal{S}}\right) \frac{Eo}{A},\tag{25}$$

$$A\rho\_l \frac{\partial \mathcal{p}}{\partial \mathbf{x}} = Eo \left(\rho\_\mathcal{g} + \frac{V\_b \mathcal{C}\_w}{\mathcal{g}}\right) - R\_{bc} A \rho\_l \left(\mathcal{C}\_l \mathcal{W}\_s \frac{\partial u\_l}{\partial z} + \mathcal{C}\_{nm} \mathcal{W}\_s \frac{\partial (u\_b - u\_l)}{\partial y}\right) \tag{26}$$

along with the consequential nondimensional boundary limitations

$$
\mu\_l(-h) = 0,\\
\mu\_l(h) = 0,\\
\varepsilon\_\S(0) = 1.\tag{27}
$$

#### **3. Mathematical Solutions and Results**

To compute the solutions of complicated nonlinear coupled equations, a powerful and efficient technique called the homotopy perturbation method (HPM) was used for finding analytic solutions. HPM is a powerful method which works even without the need of a linearization process [48–50]. It tends to reduce the nonlinear equations into a system of linear equations and generates an asymptotic solution. To serve the purpose, the initial guess was formed as

$$
\mu\_{l,0} = \frac{1}{2} \left( -2 - h^2 P + Py^2 \right) \\
\text{and } \mu\_{b,0} = 1 + Py. \tag{28}
$$

The linear operators were

$$L\_1 = \frac{d^2 u\_l}{dy^2}, \text{ and } L\_2 = \frac{du\_b}{dy}. \tag{29}$$

From Equations (22)-(23)

$$Eq1(\mathbf{x}, \mathbf{y}, \mathbf{q}) = P(\mathbf{x}, \mathbf{q}) - \frac{\partial^2 \mu\_l}{\partial y^2}(\mathbf{x}, y, \mathbf{q}) - \left(1 - \frac{c\_1}{\rho\_l \mu\_b(\mathbf{x}, \mathbf{y}, \mathbf{q})}\right) \frac{Eo}{A} \tag{30}$$

$$\begin{aligned} \text{Eq2(x,y,q)} &= -A\rho\_l P(\mathbf{x,q}) + Eo(\rho\_\mathcal{S}(\mathbf{x,y,q}) + \frac{V\_b \mathcal{L}\_w}{\mathcal{S}}) - \\ R\_{bc} A\rho\_l \big( \mathcal{C}\_l \mathcal{W}\_s \frac{\partial u\_l(\mathbf{x,y,q})}{\partial z} + \mathcal{C}\_{am} \mathcal{W}\_s \frac{\partial (u\_b(\mathbf{x,y,q}) - u\_l(\mathbf{x,y,q}))}{\partial y} \big) \end{aligned} \tag{31}$$

Constructing the homotopy that satisfy

$$H(\mathbf{U}, \mathbf{q}) = (1 - q) \Big( L\_1(\mathbf{U}) - L\_2(\mathbf{u}\_{l,0}) \Big) + q(Eq1(\mathbf{U})) = 0,\tag{32}$$

$$H(V, \mathbf{q}) = (1 - q) \Big( L\_2(V) - L\_2(\mathbf{u}\_{b,0}) \Big) + q(Eq2(V)) = 0,\tag{33}$$

with *q* ∈ (0, 1). When *q* = 0, the equation provided an initial guess, and for *q* = 1 the equation generated a required solution. The solution should be of the form

$$\mathcal{U}\mathcal{U} = \mathcal{U}\_0 + p\mathcal{U}\_1 + \dots,\\ \mathcal{V} = V\_0 + pV\_1 + \dots \tag{34}$$

Setting *p* = 1 the solution will be

$$\mathcal{U}I = \mathcal{U}I\_0 + \mathcal{U}\_1 + \dots, \mathcal{V} = V\_0 + V\_1 + \dots \tag{35}$$

The expressions achieved up to second order were

*uf* = *<sup>K</sup>*<sup>5</sup> + *<sup>K</sup>*<sup>6</sup> *<sup>y</sup>* + K7 *<sup>y</sup>*<sup>2</sup> + *<sup>K</sup>*<sup>8</sup> *<sup>y</sup>*<sup>3</sup> + *<sup>K</sup>*<sup>9</sup> *<sup>y</sup>*<sup>4</sup> <sup>−</sup> *<sup>K</sup>*<sup>10</sup> *<sup>y</sup>*5, K5 <sup>=</sup> <sup>−</sup><sup>1</sup> <sup>−</sup> <sup>h</sup>2P <sup>2</sup> <sup>−</sup> <sup>1</sup> <sup>12</sup>Eoh4P2 <sup>−</sup> <sup>1</sup> <sup>24</sup>CamEoh4PRbeWs <sup>+</sup> <sup>1</sup> <sup>24</sup>ClEoh4PRbeWs <sup>−</sup> CwEo2h4PVb 24gρ<sup>l</sup> , *K*<sup>6</sup> = <sup>1</sup> <sup>3</sup>Eo*h*2*<sup>P</sup>* <sup>−</sup> <sup>1</sup> <sup>6</sup>C*am*Eo*h*2*P*R*be*W*<sup>s</sup>* <sup>−</sup> <sup>C</sup>*w*Eo2*h*2V*<sup>b</sup>* <sup>3</sup>*g*ρ*<sup>l</sup>* <sup>−</sup> 3C*w*Eo2*h*4*P*V*<sup>b</sup>* <sup>40</sup>*g*ρ*<sup>l</sup>* , *K*<sup>7</sup> = *<sup>P</sup>* <sup>2</sup> , *<sup>K</sup>*<sup>8</sup> <sup>=</sup> <sup>−</sup>Eo*<sup>P</sup>* <sup>3</sup> <sup>+</sup> <sup>1</sup> <sup>6</sup>C*am*Eo*P*W*<sup>s</sup>* <sup>+</sup> <sup>C</sup>*w*Eo2V*<sup>b</sup>* <sup>3</sup>*g*ρ*<sup>l</sup>* <sup>+</sup> <sup>C</sup>*w*Eo2*h*2*P*V*<sup>b</sup>* <sup>12</sup>*g*ρ*<sup>l</sup>* , *<sup>K</sup>*<sup>9</sup> = Eo*P*<sup>2</sup> <sup>12</sup> <sup>+</sup> <sup>1</sup> <sup>24</sup>C*am*Eo*P*R*be*W*<sup>s</sup>* <sup>−</sup> <sup>1</sup> <sup>24</sup>C*l*Eo*P*R*be*W*<sup>s</sup>* <sup>+</sup> <sup>C</sup>*w*Eo2*P*V*<sup>b</sup>* <sup>24</sup>*g*ρ*<sup>l</sup>* , *<sup>K</sup>*<sup>10</sup> <sup>=</sup> <sup>C</sup>*w*Eo2*P*V*<sup>b</sup>* 120*g*ρ*<sup>l</sup>* (36)

*ub* = 1 + *K*<sup>11</sup> *y* + *K*<sup>12</sup> *y*<sup>2</sup> + *K*<sup>13</sup> *y*<sup>3</sup> + *K*<sup>14</sup> *y*4, *K*<sup>11</sup> = <sup>1</sup> 6*g*ρ*l* - 3C*w*EoV*<sup>b</sup>* - 4 + *h*2*P* (−<sup>2</sup> <sup>+</sup> *Cam*W*s*) <sup>+</sup> *gP*- <sup>18</sup> <sup>+</sup> <sup>C</sup>*l*Eo*h*2R*be*W*<sup>s</sup>* <sup>−</sup> <sup>C</sup>*am*- <sup>18</sup> <sup>+</sup> Eo*h*2R*be* W*<sup>s</sup>* + 6C*am*2W*<sup>s</sup>* 2 ρ*l K*<sup>12</sup> = <sup>1</sup> <sup>12</sup>*g*2ρ*<sup>l</sup>* 2 - 3C*w*2Eo2 - 4 + *h*2*P* V*b* <sup>2</sup> + C*w*Eo*gP*V*<sup>b</sup>* - −18 + Eo*h*<sup>2</sup> + 12*Cam*W*<sup>s</sup>* ρl + 6(*Cam* − C*l*)*g*2*P*R*be*W*s*(−2 + *Cam*W*s*)ρ*<sup>l</sup>* 2 K13 = <sup>1</sup> <sup>6</sup>*g*2ρ*<sup>l</sup>* 2 - Eo*P* - C*w*2EoV*<sup>b</sup>* <sup>2</sup> + C*<sup>w</sup> g*V*<sup>b</sup>* (2 + *Cam*(−1 + R*be*)W*<sup>s</sup>* − C*l*R*be*W*s*)ρl + (*Cam* − C*l*)*g*2R*be*W*s*ρ*<sup>l</sup>* 2 *<sup>K</sup>*<sup>14</sup> <sup>=</sup> <sup>−</sup> <sup>C</sup>*w*Eo2*P*V*<sup>b</sup>* (<sup>C</sup>*w*V*b*+*g*ρ*l*) <sup>24</sup>*g*2ρ*<sup>l</sup>* <sup>2</sup> . (37)

Here, *P* = <sup>∂</sup>*<sup>p</sup>* <sup>∂</sup>*<sup>x</sup>* . The flow rate in the fixed frame [51] was defined by

$$Q\_f = \int\_{-h}^{h} u\_f(x, y) dy\_r \tag{38}$$

$$Q\_{\mathcal{S}} = \int\_{-\mathrm{l}}^{\mathrm{l}} u\_{\mathcal{S}}(\mathbf{x}, \mathbf{y}) d\mathbf{y},\tag{39}$$

$$Q = Q\_f + Q\_\mathcal{S} = \int\_{-h}^{h} u\_f(x, y) + u\_\mathcal{S}(x, y) dy. \tag{40}$$

Solving the equation to get P in terms of Q and x was

$$\begin{split} \frac{\partial p}{\partial x} &= \frac{1}{4Eo\mathcal{H}^{5}\rho\_{l}} \Big( \mathcal{K}\_{1} \cdot h^{5} + \mathcal{K}\_{2} \cdot h^{3} - \sqrt{(\mathcal{K}\_{1} \cdot h^{5} + \mathcal{K}\_{2} \cdot h^{3})^{2} - (\mathcal{K}gh^{5} + \mathcal{K}\_{4}h^{6})} \Big) \\ \mathcal{K}\_{1} &= Eo(\mathcal{C}\_{l}\mathcal{g}\mathcal{R}\_{\text{hr}}\mathcal{W}\_{s}\rho\_{l} - \mathcal{C}\_{\text{w}}AEo\mathcal{V}\_{b} - \mathcal{C}\_{\text{aw}}\mathcal{g}\mathcal{R}\_{\text{hr}}\mathcal{W}\_{s}\rho\_{l}), \mathcal{K}\_{2} = -10\frac{\mathcal{G}\mathcal{P}}{A} \\ \mathcal{K}3 &= 120 \frac{\mathcal{E}\alpha\_{l}^{2}\rho\_{l}^{2}}{A}(1-Q), \mathcal{K}4 = 240 \frac{\mathcal{G}^{2}h\rho\_{l}^{2}\mathcal{E}o}{A} \end{split} \tag{41}$$

#### **4. Discussion**

To study the impact of numerous parameters, such as volume (*Vb*), Eotvos number (*Eo*), added mass coefficient (*Cam*), slip velocity (*Ws*), lift coefficient (*Cl*), and model coefficient (*Cw*), Figures 2–19 were plotted.

**Figure 3.** Pressure rise for several values of *Cam*.

**Figure 5.** Pressure rise for several values of *Ws*.

**Figure 6.** Fluid velocity curves for several values of *Cw*. Rbe = 10, *Vb* = 1.0, *Eo* = 1.0, *Cam* = 1.0, *Ws* = 10, *Cl* = 10.

**Figure 7.** Fluid velocity curves for several values of*Cl*. Rbe = 10, *Vb* = 1.0, *Eo* = 1.0, *Cam* = 1.0, *Ws* = 10, *Cw* = 10.

**Figure 8.** Fluid velocity curves for several values of *Eo*. Rbe = 10, *Vb* = 1.0, *Cl* = 10, *Cam* = 1.0, *Ws* = 10, and *Cw* = 1.0.

**Figure 9.** Fluid velocity curves for several values of *Vb*. Rbe = 10, *Eo* = 1.0, *Cl* = 10, *Cam* = 1.0, *Ws* = 10, and *Cw* = 1.0.

**Figure 10.** Fluid velocity curves for several values of*Cam*. Rbe = 10, *Eo* = 1.0, *Cl* = 10, *Vb* = 1.0, *Ws* = 10, and *Cw* = 1.0.

**Figure 11.** Fluid velocity curves for several values of Rbe. *Cam* = 1.0, *Eo* = 1.0, *Cl* = 10, *Vb* = 1.0, *Ws* = 10, and *Cw* = 1.0.

**Figure 12.** Fluid velocity curves for several values of *Ws*. *Cam* = 1.0, *Eo* = 1.0, *Cl* = 10, *Vb* = 1.0, Rbe = 10, and *Cw* = 1.0.

**Figure 13.** Bubble velocity curves for several values of *Cw*. Rbe = 10, *Vb* = 1.0, *Eo* = 1.0, *Cam* = 1.0, *Ws* = 10, and *Cl* = 10.

**Figure 14.** Bubble velocity curves for several values of *Cl*. Rbe = 10, *Vb* = 1.0, *Eo* = 1.0, *Cam* = 1.0, *Ws* = 10, and *Cw* = 1.0.

**Figure 15.** Bubble velocity curves for several values of *Eo*. Rbe = 10, *Vb* = 1.0, *Cl* = 10, *Cam* = 1.0, *Ws* = 10, and *Cw* = 1.0.

**Figure 16.** Bubble velocity curves for several values of *Vb*. Rbe = 10, *Eo* = 1.0, *Cl* = 10, *Cam* = 1.0, *Ws* = 10, and *Cw* = 1.0.

**Figure 17.** Bubble velocity curves for several values of *Cam*. Rbe = 10, *Eo* = 1.0,*Cl* = 10, *Vb* = 1.0, *Ws* = 10, and *Cw* = 1.0.

**Figure 18.** Bubble velocity curves for several values of Rbe. *Cam* = 1.0, *Eo* = 1.0,*Cl* = 10, *Vb* = 1.0, *Ws* = 10, and *Cw* = 1.0.

**Figure 19.** Bubble velocity curves for several values of *Ws*. Rbe = 10, *Vb* = 1.0, *Eo* = 1.0,*Cam* = 1.0, *Cw* = 10, and *Cl* = 10.
