2.1. Set-up

Figure 1(a) shows the schematic of the experimental set-up. The vertical channel is made of transparent acrylic plastic. Water is driven upward in the vertical channel by a pump. The dimensions of the horizontal section of the channel are $40\times 400$ mm and the channel height is 2.5 m. We define the frame of reference ($\boldsymbol {e}_x, \boldsymbol {e}_y, \boldsymbol {e}_z$) such that $\boldsymbol {e}_x$ and $\boldsymbol {e}_y$ are aligned with the spanwise and vertical directions of the channel, respectively. Monodisperse bubbles of 1.0 mm diameter are injected from a bubble generator attached in the channel at a section 1.0 m above the lower end of the channel. The gas flow rate and the duty cycle of the gas injection are dynamically controlled. Further details of the experimental set-up are provided in Appendix A.

Figure 1. (a) Schematic of the experimental set-up of high-speed imaging. (b) Representative image of bubbles captured in the test section.

We capture the images of bubbles rising near the wall using a high-speed camera (Photoron SA-2). The test section of the camera is square with a dimension of $60\times 60$ mm and placed at the middle of the front surface of the channel. The edges of the section are parallel to the vertical and spanwise directions of the channel. The vertical distance between the centre of the test section and the bubble generator, $h$, can be varied. Throughout the study, the camera continuously captures images with a frame rate of 1000 f.p.s. with a resolution of $2048\times 2048$ pixels. A representative image is shown in figure 1(b). The liquid flow rate is set such that the bulk Reynolds number is $Re_B=5100$ without gas flow and 1-Pentanol is added to the water with a concentration of 20 p.p.m. to prevent the coalescence of bubbles. Due to the contamination by surfactant, the drag of the bubbles is slightly increased compared with a clean spherical bubble with the same size. In the present setting, the mean vertical velocity of isolated bubbles rising near the wall is $205\ \textrm {mm} \ \textrm {s}^{-1}$, which corresponds to $Re=205$. The detailed statistics of the velocity for single and pairwise bubbles are addressed in the following sections.

2.2. Configuration of pairwise bubbles

We use the set-up to capture the configurations of pairwise bubbles. The test section is placed at $h=0.3$ m. To minimize the disturbance in the flow field caused by the injection of air as well as to maintain a low void fraction, we intermittently injected the gas, for a period of 0.04 s with a duty cycle of 25. By this operation, sparse bubbles are periodically released into the channel without forming clusters. From the high-speed images, we extract isolated pairs of bubbles with a condition that no bubble is present within 10 diameters from the centres of both bubbles. The bubbles captured in the test section are nearly spherical throughout the measurement. The mean diameter of the sampled bubbles is 1.0 mm with a variation of 4.9 %.

As shown in figure 2, we define the configuration of pairwise bubbles, namely bubble-$i$ and bubble-$j$, in $x$–$y$–$z$ Cartesian coordinates with its origin located at the middle point of the line connecting the centres of the bubbles, $O_i$ and $O_j$. The separation distance between the centres is $O_iO_j=d$. We define the non-dimensional separation distance as $S=d/R$. The pairwise angle $\theta _{ij}$ is defined between the $x$-axis and the line connecting the centres of the bubbles such that $-{\rm \pi} /2\leqslant \theta _{ij}\leqslant {\rm \pi}/2$. Additionally, we define the $\xi$–$\eta$–$z$ Cartesian coordinate system, which shares its origin with the $x$–$y$–$z$ coordinate system. The $\xi$-axis and the $\eta$-axis are respectively parallel and orthogonal to the line connecting the centres of the bubbles. The $x$–$y$–$z$ coordinates and the $\xi$–$\eta$–$z$ coordinates are associated with the two frames of reference, ($\boldsymbol {e}_x, \boldsymbol {e}_y, \boldsymbol {e}_z$) and ($\boldsymbol {e}_{\xi }, \boldsymbol {e}_{\eta }, \boldsymbol {e}_z$), and $\boldsymbol {U}_i$ and $\boldsymbol {U}_j$ represent the velocities of bubble-$i$ and bubble-$j$, respectively. Gravity is taken as downward along the $y$-axis as $\boldsymbol {g}=-g\boldsymbol {e}_y$. In the images, the bubbles smoothly move near the wall without violent collisions with the wall. As mentioned, the migrations of the bubbles toward the wall can be explained by the shear-induced lift force from the mean flow. Therefore, we can assume that the wall-normal velocity of the bubble, $U_z$, is zero. This is the key ansatz made based on this experiment that will be used to simplify our model that will be introduced in the following sections.

Figure 2. Configuration of a pair of bubbles.

To analyse the angle configuration of the pairwise bubbles and its dependence on the inter-bubble distance, we define the conditional pair probability distribution function (C-PPDF) as

(2.1)\begin{equation} \varGamma_{C}(S, \theta) =\frac{\varOmega_{C}}{N_p^{2}}\sum_{k=1}^{N_p}\sum_{l=1}^{N_p}\delta(S-S_k)\delta(\theta-|\theta_l|), \end{equation}

where $N_p$ is the number of pairs sampled, $\varOmega _{C}$ is a normalization constant and $H$ is the Heaviside step function; $S_k$ and $\theta _l$ are the inter-bubble distance of the $k$th pair and the pairwise angle of the $l$th pair, respectively. For convenience, we normalize the C-PPDF at each $S$ to obtain the conditional angular pair distribution function (C-APDF) as

(2.2)\begin{equation} G_{C}(\theta)=\frac{\varGamma_{C}(S,\theta)}{\displaystyle\int \varGamma_{C}(S,\theta)\,\textrm{d}\theta}, \end{equation}

which corresponds to the probability of observing pairs with a pairwise angle of $\theta$, among those with an inter-bubble distance of $S$. To obtain $\varGamma _{C}$ and $G_{C}$ from discrete bubbles, similar to an approach of Bunner & Tryggvason (Reference Bunner and Tryggvason2002), we approximate $\delta (S)$ and $\delta (\theta )$ by averaging over a circular strip with a thickness of $\Delta S$ with a radius of $S$, and a circular section with a central angle of $\Delta \theta$, respectively. We use $\Delta r=1.0$ and $\Delta \theta ={\rm \pi} /10$.

Figure 3 shows the obtained C-APDFs at various $S$. For pairs at $5/2\leqslant S<9/2$, more than half of the pairs take the pairwise angle within $0<\theta <{\rm \pi} /10$. For pairs at $11/2< S<17/2$, the probability is more uniform about $\theta$, but peaks appear at $2{\rm \pi} /5<\theta <9{\rm \pi} /10$. Thus, pairs with a short inter-bubble distance tend to take the side-by-side configurations while those with a long inter-bubble distance tend to take nearly in-line, oblique configurations. The peak decays as $S$ becomes large to approach $10$. These biases in the pairwise angle indicate that the motions of bubbles are subjected to the inter-bubble interactions.

Figure 3. Conditional angular pair distribution function for isolated pairs obtained from the experimental images.

Figure 4 shows the trajectories of representative bubbles captured in the images. The trajectories are not perfectly straight along $\boldsymbol {e}_y$ perpendicular to gravity but fluctuated. The fluctuations are observed for bubbles regardless of their proximity to others, indicating that the bubbles are subjected to turbulent agitation.

Figure 4. Trajectories of representative bubbles captured in the test section.

3.1. Formulation

We model the translational motions of the Lagrangian point bubbles shown in figure 2. This formulation extends that reported by Maeda et al. (Reference Maeda, Date, Sugiyama, Takagi and Matsumoto2013) to model the effect of turbulent agitation. Throughout the modelling, we assume that bubbles are spherical. This assumption is based on the experimental observation and is also consistent with the aforementioned numerical and experimental studies on the motion of bubbles rising at $Re=O(100)$. The motion of a spherical bubble immersed in a fluid can be formulated as a point mass undergoing a translational motion (Auton, Hunt & Prud'Homme Reference Auton, Hunt and Prud'Homme1988; Magnaudet & Eames Reference Magnaudet and Eames2000; Merle, Legendre & Magnaudet Reference Merle, Legendre and Magnaudet2005). From this perspective, the equation of motion of an isolated spherical bubble of radius $R$ rising under buoyancy at velocity $\boldsymbol {U}$ in an incompressible viscous fluid can be expressed as

(3.1)\begin{equation} (\rho_\textrm{{B}}+C_\textrm{{M}}\rho_\textrm{{L}}){V}_\textrm{{B}}\frac{\textrm{d}\boldsymbol{U}}{\textrm{d}t} ={-}{\frac{3}{8}}\rho_\textrm{{L}}\frac{{V}_\textrm{{B}}}{{R}}{C}_\textrm{{D}}|\boldsymbol{U}|\boldsymbol{U} +(\rho_\textrm{{B}-} \rho_\textrm{{L}}){V}_\textrm{{B}}\boldsymbol{g}, \end{equation}

where $\rho _{\textrm {B}}$ and $\rho _{\textrm {L}}$ are the densities of the bubble and the fluid, respectively, $C_{\textrm {D}}$ is the drag coefficient, $C_{\textrm {M}}$ is the added-mass coefficient and $V_{\textrm {B}}$ is the bubble volume. The history force is neglected compared with the quasi-steady drag since we consider regimes in that the acceleration is not large compared with $\boldsymbol {U}/2R$ (Merle et al. Reference Merle, Legendre and Magnaudet2005). Considering that $\rho _b\ll \rho _L$ and $C_{\rm M}=1/2$ (Batchelor Reference Batchelor1967), the equation can be simplified as

(3.2)\begin{equation} \frac{1}{2}\rho_\textrm{{L}}{V}_\textrm{{B}}\frac{\textrm{d}\boldsymbol{U}}{\textrm{d}t} ={-}{\frac{3}{8}}\rho_\textrm{{L}}\frac{{V}_\textrm{{B}}}{{R}}{C}_\textrm{{D}}|\boldsymbol{U}|\boldsymbol{U} - \rho_\textrm{{L}}{V}_\textrm{{B}}\boldsymbol{g}, \end{equation}

where the first and second terms on the right-hand side correspond to the drag and the buoyancy force, respectively. For a spherical bubble moving in an incompressible viscous fluid, Levich (Reference Levich1962) derived ${C}_\textrm {{D}}=Re/{48}$ by considering the global kinetic energy balance in viscous potential flow. Moore (Reference Moore1963) derived ${C}_\textrm {{D}}=Re/{48}({1}-M_{\infty }/Re^{1/2})$, where the correction $M_{\infty }=2.211\dots$ models the effect of the viscous diffusion in the boundary layer at the bubble surface and the viscous wake formed behind the bubble formed due to the separation of the layer.

For later convenience, we decompose the drag coefficient into the potential and viscous contributions as ${C}_\textrm {{D}}={C}_{\textrm {Dpot}}+{C}_{\textrm {Dvis}}$, where ${C}_{\textrm {Dpot}}={48}/Re$ and ${C}_{\textrm {Dvis}}=M_{\infty }/{48}Re^{3/2}$. For the present pairwise bubbles, we consider additional terms that represent the interaction force and turbulent agitation. The equation of motion for bubble-$i$ is expressed as

(3.3)\begin{equation} \frac{1}{2}\rho_\textrm{{L}}{V}_\textrm{{B}}\frac{\textrm{d}\boldsymbol{U}_{i}}{\textrm{d}t} ={-}{\frac{3}{8}}\rho_\textrm{{L}}\frac{{V}_\textrm{{B}}}{{R}}{C}_\textrm{{D}}|\boldsymbol{U}_{i}|\boldsymbol{U}_{i} - \rho_\textrm{{L}}{V}_{{B}}\boldsymbol{g} + \boldsymbol{F}_{\textrm{{Int}}{ij}} + \boldsymbol{F}_{Wi},\end{equation}

where $\boldsymbol {F}_{\textrm {Int}ij}$ represents the interaction force exerted on bubble-$i$ due to the presence of bubble-$j$, and $\boldsymbol {F}_{Wi}$ represents the fluctuation due to turbulent agitation. Note that the added-mass coefficient of pairwise bubbles under mutual interactions can deviate from $1/2$. This deviation is small for the regime considered in our study. As a reference, for bubbles rising side by side at $Re=O(100)$, it has been shown that $C_{\rm M}$ of both bubbles increases from $1/2$ with a decrease in the inter-bubble distance. Meanwhile, the factor of this increase is $O(0.1)$ at $S=2.5$ and rapidly decays with inter-bubble distance like $S^{-3}$ (Legendre et al. Reference Legendre, Magnaudet and Mougin2003). We consider that the correction to $C_{\rm M}$ is negligible for our analysis and use $C_{\rm M}=1/2$ for simplicity.

We decompose $\boldsymbol {F}_{\textrm {Int}ij}$ into the drag and lift components as

(3.4)\begin{equation} \boldsymbol{F}_{\textrm{{Int}}{ij}}=\boldsymbol{F}_{\textrm{{D}}_{\textrm{Int}}{ij}}+\boldsymbol{F}_{{\textrm{L}_{\textrm{Int}}}{ij}}. \end{equation}

We further decompose $\boldsymbol {F}_{{\textrm {D}_{\textrm {Int}}}{ij}}$ and $\boldsymbol {F}_{{\textrm {L}_{\textrm {Int}}}{ij}}$ into the contribution of the potential flow (potential interaction) and that of the viscosity (viscous interaction) as

(3.5)$$\begin{gather} \boldsymbol{F}_{{\textrm{D}_{\textrm{Int}}}{ij}}=\boldsymbol{F}_{{\textrm{D}_{\textrm{Int}}}{ij}\textrm{pot}}+\boldsymbol{F}_{{\textrm{D}_{\textrm{Int}}}{ij} \textrm{vis}}, \end{gather}$$

(3.6)$$\begin{gather}\boldsymbol{F}_{{\textrm{L}_{\textrm{Int}}}{ij}}=\boldsymbol{F}_{{\textrm{L}_{\textrm{Int}}}{ij}\textrm{pot}}+\boldsymbol{F}_{{\textrm{L}_{\textrm{Int}}}{ij}\textrm{vis}}. \end{gather}$$

We normalize each term on the right-hand-sides of these equations by $1/2\rho |\boldsymbol {U}_i|^{2}{\rm \pi} R^{2}$ to obtain the corresponding potential and viscous contributions of the drag and lift coefficients as ${C}_{{\textrm {D}_{\textrm {Int}}}{ij} \textrm {pot}}$, ${C}_{{\textrm {L}_{\textrm {Int}}}{ij} \textrm {pot}}$, ${C}_{{\textrm {D}_{\textrm {Int}}}{ij} \textrm {vis}}$ and ${C}_{{\textrm {L}_{\textrm {Int}}}{ij} \textrm {vis}}$.

To further formulate the forces acting on the interacting bubbles, we invoke two key assumptions. First, based on the experimental observation, we assume that the bubbles are trapped by the wall and the wall-normal motions are negligible, to the accuracy of our interest. To this end, we assume for simplicity that $\boldsymbol {F}_{\textrm {{Int}}}\boldsymbol {\cdot }\boldsymbol {e}_z=\boldsymbol {F}_{W}\boldsymbol {\cdot }\boldsymbol {e}_z\approx 0$. Second, we neglect the variation of the drag coefficient due to the presence of the wall. Shi et al. (Reference Shi, Rzehak, Lucas and Magnaudet2020) has recently quantified the increase of the drag coefficient of isolated spherical bubbles moving in a linear shear flow near a plane wall, with a relative shear rate of up to 0.5. The relative drag increase identified in their study is as small as 0.1 for a bubble moving at $Re=O(100)$ with its surface a quarter-radius away from the wall, compared to the bubble moving in an irrotational, unbounded flow. The consideration of such a slight increase in the drag may not drastically improve the overall accuracy of the model prediction, among other uncertainties. Based on this reason, in the following sections we express the forces based on formulations for bubbles in an unbounded flow.

3.2. Potential interaction

We denote the velocity potential of fluid surrounding the bubbles as $\phi$. To induce ${C}_{{\textrm {D}_{\textrm {Int}}}{ij}\textrm {pot}}$ and ${C}_{{\textrm {L}_{\textrm {Int}}}{ij}\textrm {pot}}$, we employ a bi-polar decomposition of $\phi$ on the ($\xi , \eta , z$) coordinate (Van Wijngaarden Reference Van Wijngaarden1976; Kok Reference Kok1993)

(3.7)\begin{equation} \phi=\phi_i+\phi_j, \end{equation}

where $\phi _i$ and $\phi _j$ are the velocity potentials expanded from $O_i$ and $O_j$, respectively. We define the spherical coordinate systems, $r_i$–$\psi _i$–$\varphi$ and $r_j$–$\psi _j$–$\varphi$, with their origins located at $O_i$ and $O_j$, respectively. Using Legendre polynomials and associated Legendre polynomials, $\phi _i$ and $\phi _j$ can be expressed as

(3.8)$$\begin{gather} \phi_i=\sum_{n=1}^{\infty}\left(\frac{R}{r_i}\right)^{n+1}\{g^{(1)}_{in}P_n(\cos{\psi_i})+(g^{(2)}_{in}\cos{\varphi}+g^{(3)}_{in}\sin{\varphi})P^{1}_n(\cos{\psi_i})\}, \end{gather}$$

(3.9)$$\begin{gather}\phi_j=\sum_{n=1}^{\infty}\left(\frac{R}{r_j}\right)^{n+1}\{g^{(1)}_{jn}P_n(\cos{\psi_j})+(g^{(2)}_{jn} \cos{\varphi}+g^{(3)}_{jn}\sin{\psi})P^{1}_n(\cos \psi_j)\}, \end{gather}$$

where $g^{(k)}_{in}$, and $g^{(k)}_{jn}$ ($k=1, 2, 3$) are coefficients determined by the boundary condition at the surface of the bubbles. For convenience, we decompose the interaction force $\boldsymbol {F}_{{\textrm {Int}}{ij}{\textrm {pot}}}$ into components corresponding to the forces along the $\xi$- and $\eta$-axes as $\boldsymbol {F}_{\textrm {{Int}}{ij}\textrm {{pot}}}={\boldsymbol {F}}_{\xi _\textrm {{Int}}ij \textrm {pot}}+{\boldsymbol {F}}_{\eta _{\textrm {Int}}ij \textrm {pot}}$. Here, ${\boldsymbol {F}}_{\xi _{\textrm {Int}}ij \textrm {pot}}$ and ${\boldsymbol {F}}_{\eta _{\textrm {Int}}ij \textrm {pot}}$ are expressed as polynomials in terms of $S$

(3.10)$$\begin{gather} {\boldsymbol{F}}_{\xi_\textrm{Int}ij \textrm{pot}}=f_{\xi_\textrm{Int}ij\textrm{pot}}^{({-}4)}S^{{-}4}+f_{\xi_\textrm{Int}ij\textrm{pot}}^{({-}7)} S^{{-}7}+\cdots, \end{gather}$$

(3.11)$$\begin{gather}{\boldsymbol{F}}_{{\eta_\textrm{Int}ij \textrm{pot}}}=f_{\eta_\textrm{Int}ij\textrm{pot}}^{({-}4)}S^{{-}4}+f_{\eta_\textrm{Int}ij \textrm{pot}}^{({-}7)}S^{{-}7}+\cdots, \end{gather}$$

where

(3.12)$$\begin{gather} f_{\xi_\textrm{Int}ij\textrm{pot}}^{({-}4)}=\rho{\rm \pi} R^{2}({-6{U}^{2}_{j\xi}+3{U}_{i\eta}{U}_{j\eta}}), \end{gather}$$

(3.13)$$\begin{gather}f_{\xi_\textrm{Int}ij\textrm{pot}}^{({-}7)}=\rho{\rm \pi} R^{2} \left({-}6U^{2}_{i\xi}+12U_{i\xi}U_{j\xi}+6U^{2}_{j\xi}+\tfrac{3}{2}{U^{2}_{i\eta}}+\tfrac{3}{2}{U^{2}_{j\eta}}\right), \end{gather}$$

(3.14)$$\begin{gather}f_{\eta_\textrm{Int}ij\textrm{pot}}^{({-}4)}=\rho{\rm \pi} R^{2} 3{U}_{j\xi}({U}_{i\eta}+{U}_{j\eta}), \end{gather}$$

(3.15)$$\begin{gather}f_{\eta_\textrm{Int}ij\textrm{pot}}^{({-}7)}=\rho{\rm \pi} R^{2}\left\{\left(-\tfrac{9}{2}U_{i\xi}+3U_{j\xi}\right)U_{i\eta}-\tfrac{3}{2}U_{j\xi}U_{j\eta}\right\}. \end{gather}$$

Normalizing ${\boldsymbol {F}}_{\xi _\textrm {Int}ij \textrm {pot}}$ and ${\boldsymbol {F}}_{\xi _\textrm {Int}ij \textrm {pot}}$ by $1/2\rho |\boldsymbol {U}_i|^{2}{\rm \pi} R^{2}$, we define the dimensionless coefficients ${C}_{\xi _\textrm {Int}ij \textrm {pot}}$ and ${C}_{\eta _\textrm {Int}ij \textrm {pot}}$ expressed as

(3.16)$$\begin{gather} {C}_{\xi_\textrm{Int}ij \textrm{pot}}=c_{\xi_\textrm{Int}ij\textrm{pot}}^{({-}4)}S^{{-}4}+c_{\xi_\textrm{Int}ij\textrm{pot}}^{({-}7)} S^{{-}7}+\cdots, \end{gather}$$

(3.17)$$\begin{gather}{C}_{{\eta_\textrm{Int}ij \textrm{pot}}}=c_{\eta_\textrm{Int}ij\textrm{pot}}^{({-}4)}S^{{-}4}+c_{\eta_\textrm{Int}ij\textrm{pot}}^{({-}7)}S^{{-}7}+\cdots, \end{gather}$$

where

(3.18)$$\begin{gather} c_{\xi_\textrm{Int}ij\textrm{pot}}^{({-}4)}=\frac{2f_{\xi_\textrm{Int}ij\textrm{pot}}^{({-}4)}}{\rho{\rm \pi} R^{2}|\boldsymbol{U}_i|^{2}}=\frac{-12{U}^{2}_{j\xi}+6{U}_{i\eta}{U}_{j\eta}}{|\boldsymbol{U}_{i}|^{2}}, \end{gather}$$

(3.19)$$\begin{gather}c_{\xi_\textrm{Int}ij\textrm{pot}}^{({-}7)}=\frac{2f_{\xi_\textrm{Int}ij\textrm{pot}}^{({-}7)}}{\rho{\rm \pi} R^{2}|\boldsymbol{U}_i|^{2}}=\frac{-12U^{2}_{i\xi}+24U_{i\xi}U_{j\xi}+12U^{2}_{j\xi}+3{U^{2}_{i\eta}}+3{U^{2}_{j\eta}}}{|\boldsymbol{U}_{i}|^{2}}, \end{gather}$$

(3.20)$$\begin{gather}c_{\eta_\textrm{Int}ij\textrm{pot}}^{({-}4)}=\frac{2f_{\eta_\textrm{Int}ij\textrm{pot}}^{({-}4)}}{\rho{\rm \pi} R^{2}|\boldsymbol{U}_i|^{2}}=\frac{6{U}_{j\xi}({U}_{i\eta}+{U}_{j\eta})}{|\boldsymbol{U}_{i}|^{2}}, \end{gather}$$

(3.21)$$\begin{gather}c_{\eta_\textrm{Int}ij\textrm{pot}}^{({-}7)}=\frac{2f_{\eta_\textrm{Int}ij\textrm{pot}}^{({-}7)}}{\rho{\rm \pi} R^{2}|\boldsymbol{U}_i|^{2}}=\frac{({-}9U_{i\xi}+6U_{j\xi})U_{i\eta}-3U_{j\xi}U_{j\eta}}{|\boldsymbol{U}_{i}|^{2}}. \end{gather}$$

Using the relations between the $x$–$y$–$z$ coordinates and the $\xi$–$\eta$–$z$ coordinates, we obtain

(3.22)$$\begin{gather} {C}_{{\rm D}_\textrm{Int}ij \textrm{pot}}={-}{C}_{\xi_\textrm{Int}ij \textrm{pot}}\cos{\theta}+{C}_{\eta_\textrm{Int}ij \textrm{pot}}\sin{\theta}, \end{gather}$$

(3.23)$$\begin{gather}{C}_{{\rm L}_\textrm{Int}ij \textrm{pot}}={-}{C}_{\xi_\textrm{Int}ij \textrm{pot}}\sin{\theta}+ {C}_{\eta_\textrm{Int}ij \textrm{pot}}\cos{\theta}. \end{gather}$$

Further details of the formulations are provided in Appendix B.

3.3. Viscous interaction

To express the viscous contributions of the lift and drag forces due to the inter-bubble interaction, ${\boldsymbol {F}}_{{\rm D}_\textrm {Int}ij \textrm {vis}}$ and ${\boldsymbol {F}}_{{\rm L}_\textrm {Int}ij \textrm {vis}}$, we follow an analytical model introduced by Hallez & Legendre (Reference Hallez and Legendre2011) (hereafter, we will recall it as the HL model). The HL model describes forces acting on pairwise spherical bubbles under mutual interactions rising at the equal velocity at $50\leqslant Re\leqslant 500$ in an incompressible viscous fluid quiescent at infinity. In particular, from the model we extract the viscous contributions of the interaction forces induced by the wake formed behind the leading bubble on the trailing bubble as well as those induced by the viscous diffusion in the boundary layers at the surface of bubbles. Although the velocities of the two bubbles in our model are not necessarily equal, they are rising upward under buoyancy at similar velocities: $\boldsymbol {U}_i\simeq \boldsymbol {U}_j\simeq {U}_{iy}\boldsymbol {e}_y$. Based on the similarity, we adopt the HL model to express ${\boldsymbol {F}}_{{\rm D}_\textrm {Int} ij \textrm {vis}}$ and ${\boldsymbol {F}}_{{\rm L}_\textrm {Int}ij \textrm {vis}}$ as

(3.24)$$\begin{gather} \boldsymbol{F}_{{\textrm{D}_{\textrm{Int}}}{ij}\textrm{vis}}\approx{-}{\frac{3}{8}}\rho_{{\textrm{L}}}\frac{{V}_{{\textrm{B}}}}{{R}}{C}_{{\textrm{D}}_{\textrm{Int}}{ij}\textrm{vis}}| \boldsymbol{U}_{{i}}|^{2}\boldsymbol{e}_{y}, \end{gather}$$

(3.25)$$\begin{gather}\boldsymbol{F}_{{\textrm{L}_{\textrm{Int}}}{ij}\textrm{vis}}\approx{\frac{3}{8}}\rho_{{\textrm{L}}}\frac{{V}_{{\textrm{B}}}}{{R}}{C}_{{\textrm{L}}_{\textrm{Int}}{ij} \textrm{vis}}|\boldsymbol{U}_{{i}}|^{2}\boldsymbol{e}_{x}, \end{gather}$$

where $C_{{\textrm {D}_{\textrm {Int}} ij\textrm {vis}}}$ and $C_{{\textrm {D}_{\textrm {Int}}ij\textrm {vis}}}$ are the viscous contributions of the drag and lift coefficient of bubble-$i$ in the presence of bubble-$j$

(3.26)$$\begin{gather} {C}_{{\textrm{D}}_{\textrm{Int}}{ij}\textrm{vis}}=\left\{\begin{array}{@{}ll} \displaystyle -M_{\infty}\dfrac{48}{Re^{{3}/{2}}}S^{{-}3}, & 0\leqslant\theta\leqslant{\rm \pi},\\ \displaystyle -(M^{*}_{2}(1+S^{3})-M_{\infty})\dfrac{48}{Re^{{3}/{2}}}, & -{\rm \pi}\leqslant\theta\leqslant 0. \end{array} \right. \end{gather}$$

(3.27)$$\begin{gather}C_{{\textrm{L}_{\textrm{Int}} ij\textrm{vis}}}=\left\{\begin{array}{@{}ll} \displaystyle -\dfrac{240}{ReS^{4}}\cos{\theta}, & 0\leqslant\theta\leqslant{\rm \pi},\\ \displaystyle -\dfrac{240}{ReS^{4}}\cos{\theta}-\dfrac{2}{3}C_{{\textrm{L}\,\textrm{wake}}}\beta Re\dfrac{\cos{\theta}}{\sin{\theta}}(1-\tilde{u}_{{\textrm{wake}}})\tilde{u}_{{rot}}, & -{\rm \pi}\leqslant\theta\leqslant 0. \end{array} \right. \end{gather}$$

When the pairwise angle is at $0\leqslant \theta \leqslant {\rm \pi}$, bubble-$i$ is at the leading position against bubble-$j$ along the $y$-axis, and viscous interaction for bubble-$i$ comes only from the viscous diffusion in the boundary layer of bubble-$j$. On the other hand, when $-{\rm \pi} \leqslant \theta \leqslant 0$, bubble-$i$ is located at the trailing position against bubble-$j$ and subjected also to the influence of the wake formed behind bubble-$j$. In the wake region, dragged by the leading bubble, the distribution of the magnitude of the fluid's velocity along the direction of the bubble's motion takes a Gaussian profile which decays asymptotically to zero in regions distant from the wake (Batchelor Reference Batchelor1967). As discussed in detail by Yuan & Prosperetti (Reference Yuan and Prosperetti1994) and Harper (Reference Harper1997), due to this local velocity distribution within the wake region, the drag coefficient of the trailing bubble in the wake region against the background flow field is effectively reduced (slip-stream effect). The coefficient $M_2^{*}$ in relation (3.26) accounts for this drag reduction. Moreover, due to the non-zero vorticity distribution in the wake, the trailing bubble is subjected to the shear-induced lift force in the direction outward from the wake. The magnitude of this shear-induced lift force is expressed as $C_{{\rm L}_{\textrm {Single}}} \boldsymbol {u}\times (\boldsymbol {\nabla \times {u}})$, where $C_{{\rm L}_{\textrm {Single}}}$ is the lift coefficient of a single isolated bubble immersed in an uniform shear flow. Auton (Reference Auton1987) analytically induced $C_{{\rm L}_{\textrm {Single}}}=1/2$ in the limit of weak shear. Legendre & Magnaudet (Reference Legendre and Magnaudet1998) extended the result as $C_{{\rm L}_{\textrm {Single}}}=(Re+16/2(Re+29))$ for $Re\geqslant 5$. The HL model modifies this expression to define the lift, which appears in the term that includes $C_{{\textrm {L}\,\textrm {wake}}}$. This wake-induced lift force plays a key role in the pairwise dynamics considered in the present study. In the following sections, in order to assess the effects of the viscous contribution to the interaction force, we simulate cases without including viscous interactions by setting ${C}_{{\textrm {D}}_{\textrm {Int}}{ij}\textrm {vis}}={C}_{{\textrm {L}}_{\textrm {Int}}{ij}\textrm {vis}}=0$, and compare results obtained with the complete model.

3.4. Turbulent agitation

The random motions of particles in turbulent flow field can be modelled as Brownian motion (Taylor Reference Taylor1922). Meanwhile, the fluctuations of the velocities of the present pairwise bubbles with a finite inter-bubble distance may be correlated, since the turbulent agitation that drives the fluctuations may have spatial structures at the scale of the inter-bubble distance. For pairwise particles that move on the streamlines of the background flow field (tracer particles), their velocity correlation corresponds to the Lagrangian velocity correlation of the background single-phase turbulence (Batchelor Reference Batchelor1953; Durbin Reference Durbin1980; Thomson Reference Thomson1990; Sawford Reference Sawford2001). On the contrary, the present bubbles have relative velocities (slip velocity) against the carrier fluid, and the resulting trajectories do not follow the streamlines of the fluid that would otherwise be present in the space filled by the bubbles. The correlation of the fluctuating velocities of bubbles therefore may not correspond to the Lagrangian velocity correlation of single-phase turbulence. Systematic derivation of such a correlation is not a simple task. The slip velocity results from the inter-phase momentum transfer at the bubble interface, and computation of this momentum transfer requires modelling of the thin boundary layer at the surface of the bubbles. The state of the boundary layer can be influenced by surfactant molecules as well as dynamically change due to interactions with turbulent eddies on the small scale. The consideration of such a complex dynamics is a subject of high-fidelity simulation and, as mentioned, is out of the scope of the present study.

Instead, we model the effect of the agitation through spatio-temporal stochastic forcing of bubbles, $\boldsymbol {F}_{Wi}$, based on an ansatz that it is Gaussian in time and coloured in space. We also crudely assume that the spatial correlation is isotropic on the $x\text {--}y$ plane. As a result, the spatial component of the correlation can be represented by a scalar function about the inter-bubble distance. Rather than treating the equation (3.3) as a stochastic differential equation, we empirically express $\boldsymbol {F}_{Wi}$ as a pseudo-stochastic process

(3.28)\begin{equation} \boldsymbol{F}_{W_i} \approx \sum^{Ns}_k C_{k}\boldsymbol{\chi}_k(\mathrm{sign}(i-j)d_{ij})\sum^{Nt}_l D_{k,l}\zeta_{k,l}(t),\end{equation}

where $\boldsymbol {\chi }$ and $\zeta$ are the spatial and temporal basis functions, respectively. Here, $C$ and $D$ are weighting coefficients that respectively define the spatial and temporal spectra of the noise. In the present study, we employ phase-randomized Fourier basis functions for both $\boldsymbol {\xi }$ and $\zeta$. This pseudo-stochastic representation of the agitation extends the expression of Kraichnan (Reference Kraichnan1970) that was used to represent the random velocity field in isotropic turbulence. In the present study, the parameters are obtained such that simulations reproduce the ensemble-averaged statistics in the experiment a posteriori, including the velocity power spectral density (PSD) and the mean square displacement (MSD) of isolated bubbles, as well as the spatial correlation of the fluctuating velocities of pairwise bubbles; $(N_s, N_t)=(10^{3}, 10^{2})$ is found sufficient to obtain statistical convergence. The comparisons of the experimental and numerical statistics are discussed in the following section. Further details on the stochastic representation and the parameters are provided in Appendix C. Throughout the study, (3.3) is integrated to simulate the motion of bubbles, using the standard fourth-order Runge–Kutta scheme (RK4) with a sufficiently small time step size.

4.1. Single isolated bubble

For isolated bubbles, the turbulent agitation is independent of the spatial coordinates of bubbles, within the assumption that the response of bubbles to the agitation is spatially homogeneous on the wall. To characterize the fluctuating velocity from the experimental data as well as to validate the temporal component of the stochastic forcing, we address the temporal statistics of single isolated bubbles.

Figure 5 compares the ensemble-averaged PSD of the fluctuating velocity of isolated bubbles obtained from the experiment and the simulation. Both results are obtained from the ensemble average of $O(10^{2})$ distinct realizations. The fluctuating velocity is defined as the mean-subtracted velocity of each bubble within a temporal window of 0.15 s. The overall profile of the modelled PSD agrees with the experimental one, although the experimental PSD is less smooth. The PSDs present a sharp cutoff at $f_c\approx 50$ Hz, which corresponds to the maximum frequency of the agitation.

Figure 5. Ensemble-averaged PSD of the velocity fluctuation of isolated bubbles.

Figure 6 compares the spanwise MSDs of the same isolated bubbles obtained in the experiment and the simulation. The overall profiles of the MSD agree with each other, except that the modelled MSD is slightly shifted upward due to the faster initial growth. For both plots, the MSD presents a faster growth near the origin. The slope then becomes linear at $t\approx 0.05$ s. The profiles of the MSD agree with the theory of turbulent diffusion (Taylor Reference Taylor1922). The initial growth corresponds to the ballistic translation of the bubbles at the time scale below that of the agitation, $1/f_c\approx 0.02$ s. The linear growth at $t>1/f_c$ corresponds to the Brownian diffusion. The initial discrepancy between the two plots could be due to the fluctuation of the (mean) drag coefficient against the sudden acceleration, which is not formulated in the present model. Nevertheless, the slope of the linear regime agrees very well with the experiment, and the difference between the two profiles is constant over time. For the purpose of reduced-order modelling in the present study, we consider that the difference is admissible for the aimed at accuracy of the model prediction.

Figure 6. Ensemble-averaged mean square spanwise displacement of isolated bubbles.

Figures 7(a) and 7(b) respectively show the evolution of the fluctuating velocity of five representative bubbles which were randomly chosen from the experiment and that from the simulation. The profiles of the two figures qualitatively agree well, indicating that the modelled velocity fluctuation reproduces the experimental one.

Figure 7. Evolution of the spanwise velocity of representative bubbles in (a) experiment and (b) simulation, within the sampling time window.

Figures 8(a) and 8(b) respectively compare the probability density functions (PDFs) of fluctuating components of the spanwise and streamwise velocities obtained from the model and the experiment. The modelled and experimental PDFs agree with each other in both figures. As expected from the modelling ansatz that the stochastic forcing is white in time, the modelled PDFs are excellently fitted by a normal distribution. Meanwhile, the experimental PDFs are non-Gaussian; their peaks at $v=0$ are greater and their tails are heavier, compared with the modelled PDFs. The heavy tail may be associated with the intermittency of turbulent fluctuations in the career fluid. These high-kurtosis features are more pronounced in the PDF of the streamwise velocity. The difference in the profiles of the spanwise and streamwise PDFs indicates the spatial anisotropy of the agitation, which is not considered in the present model. Nevertheless, similar to the argument made about the MSD, we consider that the discrepancies between the modelled and experimental PDFs are admissible for the purpose of the present study.

Figure 8. Probability density functions of the (a) spanwise velocity and (b) streamwise velocity of a single isolated bubbles.

4.2. Spatial correlation of the agitation

In order to obtain the spatial correlation of the agitation, we compute the Lagrangian two-point correlation of the fluctuating velocities of pairs of bubbles. To isolate the effect of potential interaction on the correlation, we exclude pairs with their inter-bubble distance smaller than 5$S$ at any instance on their trajectories.

Figure 9 shows the measured and modelled correlations, $R_{u'u'}=(R_{U_x'U_x'}+R_{U_y'U_y'})/2$ as well as a fit with a sinc function, $\mathrm {sinc}(Rk_cS)$. Interestingly, the obtained correlation shows reasonable agreement with the sinc function with $k_c=2.05\ \textrm {Hz}^{-1}\ \textrm {rad}$. Although other candidate functions can be use to fit to the data, this agreement suggests that the spectrum of the spatial correlation of the data can simply be approximated by a top-hat function with a cutoff wavenumber of $k_c$. The cutoff wavenumber indicates that the minimum spatial length scale of the correlation is $2{\rm \pi} /k_c\approx 3$ mm. It should again be emphasized that this integral length scale is empirical and unique to the fluctuating velocities of the bubbles, and not immediately associated with the spatial correlation of the background (single-phase) turbulent flow field. In the meantime, from an a priori point of view, the motions of the bubbles are expected to be insensitive to the velocity fluctuations whose structural length scale is of the order of or smaller than the bubble size, $O(1)$ mm. The agreement of the a priori and a posteriori integral length scales of the correlation provides additional confidence in the validity of the present model.

Figure 9. Second-order correlation function of the fluctuating velocities of pairwise bubbles, $R_{u'u'}$.

5.1. Effect of viscous wake

In order to assess the preferential configurations captured in the experiment, we simulate the motions of pairwise bubbles. We first address pairs in initially near in-line configurations without agitation. We track the motion of the trailing bubble in coordinates whose origin is placed at the centre of the trailing bubble, $(\Delta x, \Delta y) =(x_j-x, y_j-y)$. The trailing bubble is initially released at $(\Delta x_0, \Delta y_0):\Delta x_0\in [\varepsilon , 2], \Delta y_0\in [-5R, 0]$. We terminate simulations when the bubbles collide. We set the initial velocity of both bubbles as $(U_x,U_y)=(0,0.311)\ \textrm {ms}^{-1}$. The corresponding Reynolds number is $Re=311$. This velocity is the terminal rise velocity of an isolated single spherical bubble with a diameter of 1.0 mm in quiescent water, based on (3.3). We simulate cases with and without viscous interaction, and compare results.

Figure 10(a–c) shows the relative trajectories of the trailing bubble for representative pairs with various initial configurations. In figure 10(a), cases with $\Delta x_0=\varepsilon$ are shown in that the bubbles are initially in line with the small perturbation. Without viscous interaction, in all the cases the trailing bubble travels away from the leading bubble, nearly vertically due to the repulsion induced by potential flow. With viscous interaction, in the case with $\Delta y_0=4R$, the trailing bubble first travels downward and then translates in the positive $x$ direction. The bubble subsequently trails a curved path to horizontally approach the leading bubble. In the cases with $\Delta y_0=6R$ and $8R$, the trailing bubble horizontally translates in the positive $x$ direction, and then draws a similar curved path to approach the leading bubble. The initial downward motion of the trailing bubble with $\Delta y_0=4R$ indicates the dominance of the potential-induced repulsion. Figure 10(b) shows cases with $\Delta x_0=R$. Without viscous interaction, for all cases, the trailing bubble first travels diagonally downward to increase the inter-bubble distance, and then draws curved paths to approach the leading bubble. The radius of the curved path increases with the initial inter-bubble distance, unlike those in figure 10(a). With viscous interaction, the bubble initially takes a more horizontal path and then trails the curved path. For the same initial condition, the radius of the curved path becomes smaller with viscous interaction. Compared with the cases shown in figure 10(a), the difference in the trajectories due to viscous interaction is less pronounced. Figure 10(c) shows cases with $x_0=2R$. All trajectories of the trailing bubble are similar to those in 10(b); the bubbles move diagonally downward and then trail curved paths to approach the leading bubble. Viscous interaction makes almost no difference in this figure.

Figure 10. Representative trajectories of the trailing bubble relative to the leading bubble, with various initial configurations without agitation. (a–c) Present cases with distinct values of initial horizontal distances: $\Delta x_0=({a})\varepsilon$, $(b)R$ and $(c)2R$. In each panel, cases with three distinct values of lateral distance are shown: $\Delta y=({a})4R$, $(b)6R$ and $(c)8R$. For each initial condition, results (-) with and (–) without viscous interactions are plotted in the coordinates fixed on the centre of the leading bubble. Units are normalized by $R$.

Overall, the relative trajectories shown in figure 10 indicate the strong influence of the viscous wake on the pairwise dynamics. For in-line pairs for which the trailing bubble is initially at $\Delta x<2R$, the trailing bubble is pushed out from the wake region and takes significantly shorter paths to become side by side with the leading bubble. The results also support the theory that the in-line configuration is unstable, as predicted by Harper (Reference Harper1970) and Hallez & Legendre (Reference Hallez and Legendre2011). Note that when the pair is initially in a perfectly in-line configuration ($\Delta x_0 = 0$), the bubbles reach the equilibrium state such that $(\Delta x, \Delta y)=(0, -6.8R)$, with a Reynolds number of $Re_{eq}=326$. The equilibrium distance and $Re_{eq}$ correspond to those derived in previous studies (Yuan & Prosperetti Reference Yuan and Prosperetti1994; Harper Reference Harper1997; Hallez & Legendre Reference Hallez and Legendre2011).

To assess the effects of the viscous interaction on the velocity of the trailing bubble, figures 11(a) and 11(b) respectively show the trajectories with its initial positions located at $(\Delta x_0, \Delta y_0)=(\varepsilon ,-6R)$ and $(2R,-6R)$, with markers placed every $5.0\times 10^{-3}$ s. In the former initial condition, the trailing bubble is initially immersed in the wake region, while in the latter condition it is out of the wake region. In figure 11(a), the markers indicate that, during the initial horizontal motion, the trailing bubble attains its maximum velocity of $O(0.1)\ \textrm {ms}^{-1}$, relative to the leading bubble. The motion becomes slower on the subsequent curved part of the trajectory, and then it is re-accelerated as the inter-bubble distance becomes small. The initial velocity in figure 11(b) is smaller by an order of magnitude than that in figure 11(a) up to $\Delta x_0\approx 4R$, while the velocities on the later part of the trajectory are similar. The difference in the initial velocity indicates the large amplitude of the wake-induced lift force compared with the force due to potential interaction. In experiments, it is unlikely that pairs of bubbles are observed in the in-line configuration with the aforementioned equilibrium distance since small disturbances are always present which can immediately trigger this symmetry breaking. This implication agrees with the instability of the in-line configuration observed in the experiments by Katz & Maneveau (Reference Katz and Maneveau1996), Sanada et al. (Reference Sanada, Watanabe, Fukano and Kariyasaki2005) and Kusuno et al. (Reference Kusuno, Yamamoto and Sanada2019).

Figure 11. Trajectories of the trailing bubble relative to the leading bubble with initial configurations of $(\Delta x_0,\Delta y_0)= ({a})(\varepsilon ,6R)$ and $(b)(2R,-6R)$, with viscous interaction. Markers are placed every 5.0 ms. The total travel time of the bubbles is (a) 2.73 s and (b) 8.0 s.

To assess the effect of viscous interaction on the total travel time of the trailing bubble until clustering with the leading bubble, in figure 12(a,b) we respectively show contours of the travel time simulated with and without viscous interaction, when released at various relative coordinates. The agitation was not modelled in these simulations. For both cases, the travel time tends to increase with both the pairwise angle and the inter-bubble distance. With viscous interaction, the travel time is greater for initial configurations slightly out of line than that in perfectly in-line configurations. The maximum travel time in the domain is $O(10^{3})$ s at $(\Delta x,\Delta y)=(2R,-10R)$. Without viscous interaction, the travel time monotonically increases with the pairwise angle for any given inter-bubble distance. Figure 12(c) shows a contour of the difference between the cases with and without viscous interaction, normalized by the travel time in the case with viscous interaction. The region with the difference greater than unity is confined to the wake region, $\Delta x:\Delta x<2R$. On the $\Delta y$-axis, with decreasing $\Delta y$, the difference increases toward $\Delta y_0 \approx -6R$ to take a peak value of $O(10^{3})$, and then decays. In the far field, the difference is expected to converge to zero as the interaction force becomes infinitesimally small. Overall, results in figure 12 highlight the significant acceleration of the side-by-side clustering of in-line pairs by the wake-induced lift force acting on the trailing bubble.

Figure 12. Contours of the time required for the trailing bubble to catch up with the leading bubble released in the domain without agitation. Results are shown simulated (a) with viscous interaction, $T_{t,vis}$, and (b) without viscous interaction, $T_{t,invis}$. Their normalized difference, $|T_{t,invis}-T_{t,vis}|/T_{t,vis}$, is shown in (c).

5.2. Effect of turbulent agitation

We analyse the effects of the modelled agitation on the pairwise dynamics. We first address the fluctuation in the travel time of the trailing bubble during the wake-induced translation. Figure 13(a) shows the box plot of the travel time required for the trailing bubble to reach $\Delta x=2R$, when it is released at $(\Delta x, \Delta y)=(\varepsilon ,\Delta y_0)$, where $\Delta y_0 \in [-8.5R, -2.5R]$. At all $\Delta y_0$, the travel time is of the order of $O(10^{-2})$ s. This time scale is as small as the minimum time scale of the fluctuation, $1/f_c$, and the bubble travels on a ballistic path during the simulation. The median travel time is slightly smaller than that of the case without agitation at all $\Delta y_0$. The median of the travel time slightly decreases with decreasing $\Delta y_0$ up to $\Delta y_0=-4.5 R$ and then becomes nearly constant at smaller $\Delta y_0$. The fluctuation of the travel time is spread within 80 % of the mean travel time at all $\Delta y_0$. Figure 13(b) shows the probability distribution of the travel time for the case with $\Delta y_0=-2.5R$. The distribution presents a skewness of 0.85. This positive weak skewness can be explained by the growth of the potential drag $C_{D\xi ,int}$ with increasing horizontal velocity of the trailing bubble. Overall, results shown in figure 13 indicate that the time scale of the wake-induced translation is small and not influenced by the agitation, regardless of the inter-bubble distance considered.

Figure 13. (a) Box plots of the time required for the trailing bubble to reach $\Delta x=2R$, when released at $\Delta x_0=\varepsilon$ and $\Delta y_0:\Delta y_0 \in [-8.5R, -2.5R]$. Results are obtained with agitation modelled; ($-\circ -$) denotes results without agitation. (b) Probability distribution of $T_s$ for $\Delta y_0=-8.5R$.

We also point out that the order of the travel time of the trailing bubble identified here, $T_w$, agrees with the result of the recent DNS of pairwise bubbles rising in viscous liquid reported by Zhang et al. (Reference Zhang, Ni and Magnaudet2020). Zhang et al. (Reference Zhang, Ni and Magnaudet2020) characterizes the dynamics of pairwise, in-line bubbles using the Galilei number (Ga) and the Bond number (Bo). The value of these non-dimensional parameters for the pairwise bubbles involved in figure 13 in the present study corresponds to $(Ga, Bo) = (35, 0.03)$ based on the radius of the bubble, and this value nearly matches the parameter group $(Ga, Bo) = (30, 0.02\text {--}0.05)$ involved in figure 15 of Zhang et al. (Reference Zhang, Ni and Magnaudet2020). An estimation can be made from this figure that the travel time to a transverse distance $\Delta S = 2R$ of the trailing bubble is $\hat {T}_w=O(10)$ in non-dimensional units after the breaking of the in-line configuration. This non-dimensional travel time corresponds to a physical time of $T_w=\sqrt {R/g}\hat {T}_w = O(10^{-2})$ s. The estimated order of the travel time agrees with the result identified in the present study. This agreement can provide additional confidence to the accuracy of our modelling.

Next, we address the fluctuation during the trailing bubble travelling on curved paths driven by potential interaction, for oblique pairs out of the influence of the wake. Figure 14(a) shows the box plot of the time required for the trailing bubble to catch up with the leading bubble, when the bubbles are released at $\Delta x_0=2R$, with various values of $\Delta y_0$. In these conditions, the trailing bubble is subjected almost only to the potential-induced attraction and the agitation. The travel time without agitation scales like $|\Delta y_0|^{4}$. This scaling agrees with the theoretical prediction of (3.3), based on the leading order of the potential interaction, $S^{-4}$. The magnitude of the fluctuation (box height) significantly increases with decreasing $\Delta y_0$. The median values of the travel time in the case with agitation do not significantly deviate from that without agitation. The mean travel time, on the other hand, positively deviates from the travel time without agitation with a greater magnitude compared with the median. The difference between the median and the mean also increases with decreasing $\Delta y_0$. The increase in the deviation between the median and the mean suggests the growth of the skewness of the fluctuation. Figure 14(b) shows the probability distribution of the travel time at $\Delta y_0=-5R$, up to the third quantile. As predicted, the distribution is highly non-Gaussian; it spreads over greater than two decades and its profile is skewed even on the semi-logarithmic axis.

Figure 14. (a) Box plots of the time required for the trailing bubble to catch up with the leading bubble, when it is released at the coordinate of $\Delta x_0=2R$ and $\Delta y_0:\Delta y_0 \in [2.5R, 6.5S]$; ($-\circ -$) denotes results without agitation, ($-\diamond -$) denotes the box mean, ($--$) denotes $T_t\sim \Delta y_0^{4}$. (b) Probability distribution of $T_s$ for $\Delta y_0=-5R$, up to the third quantile.

To gain further insights into the fluctuation in the travel time, figure 15 shows representative trajectories of the trailing bubble relative to the leading bubble simulated with viscous interaction and agitation, with various initial configurations. While fluctuating, the trajectories are overall similar to the cases without agitation. The magnitude of the fluctuations, however, apparently becomes greater for pairs with a larger initial inter-bubble distance, suggesting that the fluctuation in the travel time results from that of the trajectories. Compared with the pairwise bubbles without agitation, these results also suggest that the travel time for near in-line pairs to become side by side in turbulence can fluctuate by orders of magnitude. The enhancement of the fluctuation in both the travel time and the trajectory against the increase in the inter-bubble distance can be explained by the growth of the relative magnitude of the agitation compared with the potential flow attraction; the amplitude of the agitation is homogeneous, while the amplitude of potential interaction rapidly decays with increasing the inter-bubble distance. The significant acceleration of the side-by-side clustering by the additive fluctuation may be counter-intuitive, compared with the deceleration. Our interpretation is that the acceleration is due to the highly path-dependent nature of the travel time. With a finite probability, the agitation can drive the trailing bubble closer to the leading bubble. In that case, the trailing bubble moves to an inner curved trajectory with a smaller radius on which the bubble travels faster. For pairs at a shallow angle, fluctuations can also drive the trailing bubble upper relative to the leading bubble, and this motion can trigger side-by-side potential attraction earlier than the case without agitation, resulting in faster clustering.

Figure 15. Representative trajectories of the trailing bubble relative to the leading bubble, with viscous interaction and agitation. The initial positions of the trailing bubble are $(\Delta x_0, \Delta y_0)= (a)(2R, -2R)$, $(b)(2R, -3R)$ and $(c)(2R, -4R)$. The trajectories in each panel are obtained from independent realizations.

The time scales of the wake-induced symmetry breaking and the potential-driven side-by-side clustering identified in this section, and their comparisons with the travel time of bubbles in the experimental imaging, $O(1)$ s, explain the preferential configurations captured in the images. The time scale of the initial symmetry breaking is $O(10^{-2})$ s, regardless of the agitation. The in-line configuration is therefore likely to be broken before reaching the experimental test section. Pairs initially at $S<5$ reach side by side and those at $S>5$ remain in the near in-line configurations, within the time scale of $O(1)$ s. Therefore, assuming that the distribution of the pairwise angle is initially uniform, the probability of observing the near in-line pairs is increased and that of the perfect-inline pairs is decreased, for pairs at $S>5$ in the test section. On the other hand, pairs initially at $S<5$ take side-by-side configurations before they reach the test section.

6.1. Onset of dilute cluster formation in the experiment

To further discuss the implications of the pairwise dynamics for the mechanism of the formation of larger bubble clusters, we capture the onset of clustering of multiple bubbles in the channel using the same experimental set-up. The bulk flow condition and the surfactant concentration follow § 2. This time the air is continuously supplied to realize an average void fraction of 0.1 %. At this void fraction, the average inter-bubble distance becomes approximately $6R$, assuming that all bubbles are trapped by the wall. To capture the evolution of the spatial distributions of bubbles, we captured instantaneous images of bubbles at three different test sections located at various heights from the bubble inlet: $h=0.1$, $0.3$ and $0.5$ m. The times required for the bubbles to reach the centre of each test section from the inlet are $t=0.5$ s, 1.0 s and 2.5 s, respectively.

Figure 16(a–c) shows the images of bubbles obtained at the test sections. To quantify the state of clusters, we compute the pair probability distribution function (PPDF)

(6.1)\begin{equation} \varGamma(S, \theta) =\frac{\varOmega}{N_b(N_b-1)}\sum_{i=1}^{N_b}\sum_{j=1, j\neq i}^{N_b}\delta(S-|S_{ij}|)\delta(\theta-|\theta_{ij}|), \end{equation}

where $N_b$ is the total number of bubbles in each image and $\varOmega$ is normalization constant. We approximate the delta functions following the approach used to obtain the C-PPDF. Notice that the physical interpretation of the PPDF is distinct from the C-PPDF, since to obtain the PPDF all bubbles are sampled from each image without conditioning. To obtain the statistical mean, we average $\varGamma$ using over $10^{3}$ images at each test section without overlap. For convenience, from the PPDF we compute the residual of the angular pair distribution function (APDF) normalized at each $S$ after subtracting unity

(6.2)\begin{equation} \Delta G(\theta)=\frac{\varGamma(S,\theta) \displaystyle\int \,\textrm{d}\theta}{\displaystyle\int\varGamma(S,\theta) \, \textrm{d}\theta}-1, \end{equation}

which represents the bias of the probability of observing pairs with a configuration angle $\theta$ at each $S$. If the distribution is uniform at all angles, $\Delta G(\theta )=0$. Figures 16(d)–16( f) respectively show $\Delta G(\theta )$ for $S=5.0$, 7.5, 10 and 12.5, which are obtained from the images captured in the corresponding test sections in figure 16(a–c). In figure 16(d), $\Delta G$ largely varies with $S$. For $S=5.0$, $\Delta G$ takes its peak with values of around 0.2 at $\theta =0$; $\Delta G$ decays with $\theta$ and becomes negative at $\theta :\theta >0.2{\rm \pi}$. For larger $S$, the distribution varies with $h$. Small peaks are observed at $\theta =0$ for $S=7.5$ and 10. For $S=12.5$, $\Delta G$ presents a weak positive correlation with $\theta$. These plots indicate that pairs with a short inter-bubble distance ($S<7.5$) tend to become side by side, while specific structures are not observed at larger scales. On the other hand, in panels (e,f) of figure 16, the profiles of the plots for various values of $S$ are similar. In all plots in both panels, $\Delta G$ takes its peak values at $\theta =0$ and almost steadily decays with $\theta$; $\Delta G$ takes negative values at $\theta >0.2$ for all plots. The magnitude of the slope decays with increasing $S$. These profiles in figure 16(e, f) indicate that bubbles tend to be aligned horizontally regardless of their mutual distance, at $h=0.3$ and 0.5 m. Moreover, the resemblance of figure 16(e, f), and their clear difference from figure 16(d) indicate that the evolution of the horizontal clustering reaches its stationary state before bubbles reach $h=0.3$ m after $h=0.1$ m, therefore at $t:0.5< t<1.5$ s. The angle dependence of $\Delta G$ and its evolution in figure 16(d–f) agree with the visual orientation of bubbles in figure 16(a–c). In figure 16(a), many pairs are in side-by-side configurations, while larger clusters are not present. In figure 16(b), on the other hand, we clearly observe long, horizontal chains of bubbles. Similar clusters are observed in figure 16(c).

Figure 16. (a–c) Images of the bubbles captured in the test sections located at different heights. (d–f) Angular pair distribution function of the bubbles obtained from the sequence of images the sections corresponding to (a–c). Symbols: $\circ$, $S=5.0$; $\diamond$, $S=7.5$; $\square$, $S=10.0$; $\times$, $S=12.5R$.

6.2. Probability of pairwise clustering

The clustering of multiple bubbles involves many-body interactions, and implications from the pairwise interaction are clearly limited. Nevertheless, potential interaction is short range and the effects of the nearest neighbours are expected to be the most dominant. Therefore, the formation of dilute clusters is plausibly controlled by a finite number of local inter-bubble interactions, and its time scale can be associated with that of pairwise clustering.

To assess the possible effects of viscous interaction and the agitation as well as the spatial correlation of the agitation on the onset of cluster formation, we simulate the dynamics of pairwise bubbles with various initial configurations during the clustering time scale identified in the experiment, $O(1)$ s, with various model settings as summarized in table 1. In these sets, we address combinations of the interaction models with and without viscous effects, and agitation models with and without spatial correlation. In the cases without the agitation (set B, C, E and F), the initial position of the trailing bubble is defined at the centre of $16\times 20$ grids that uniformly discretize the domain of $\Delta x\in [0,3.2]$ and $\Delta y\in [-4,0]$. A total of 500 simulations are conducted for the set of $(\Delta x, \Delta y)$ at each grid cell to obtain the probability of side-by-side clustering occurring within time $t$: $p_c(t,\Delta x,\Delta y)$. In the cases with the agitation, $48\times 60$ grids are used for the same domain. A single simulation is conducted at each grid cell.

Table 1. Summary of the model settings for the simulation of the probability of side-by-side clustering.

For each case in table 1, figure 17(a–f) shows contours of the probability at $t=1.5$ s, $p_{T}(t=1.5~s,\Delta x,\Delta y)$. Figure 17(a,d) shows almost no difference, suggesting the small influence of the wake on the probability, in the case without agitation. Due to the absence of agitation, the probability is binary (0 or 1) changing across the contour line. In the plots of the cases with agitation, the contours change with a finite gradient due to the smooth change of the probability. The contour lines of 0.5 and 0.9 are nearly common between these plots. The lines of 0.5 in both plots are also similar to those in figure 17(a,b). Interestingly, the contour lines of 0.1 present differences among the contours. In figure 17(b,c), the contour follows the wake profile, indicating that the probability of initially in-line pairs reaching a side-by-side position within 1.5 s is sufficiently small. On the other hand, in figure 17(b), the contour line of 0.1 is pushed down below $\Delta y =-7R$. In figure 17(c) the corresponding contour is not even present. These changes indicate that the probability for in-line bubbles is significantly increased by the presence of the wake. The difference between figure 17(b,c), and that between figure 17(e, f) is small. This result indicates that the spatial correlation of the agitation has a small influence on the probability.

Figure 17. Contours of the probability for pairs of bubbles to reach the side-by-side position within 1.5 s, in the case that the trailing bubble is initially placed at each given coordinate relative to the leading bubble in the domain. Panels (a), (b) and (c) are obtained without agitation, with spatially uncorrelated agitation and with spatially correlated agitation, respectively, without viscous interaction. Panels (d–f) are counterparts of (a–c) with viscous interaction.

Figure 18 shows the evolution of the conditional probability of side-by-side clustering to occur for in-line configurations of pairs. The conditional probability is defined as the probability field integrated over the wake region at each instant

(6.3)\begin{equation} p_{T,c}(t)=\frac{1}{\varOmega_w}\int_{\varOmega_w}p_c(t,\Delta x,\Delta y)\,\textrm{d}(\Delta x)\,\textrm{d}(\Delta y), \end{equation}

where $\varOmega _w$ is the wake region defined as $\varOmega _w:\Delta x\in [0,2R], \Delta y\in [-8R,-2.25R]$. For all cases, the probability monotonically increases with time, and the slopes of the plots grow up to $t\approx 0.5$ s and then decay. The increase is because pairs with larger inter-bubble distance reach side by side at a later time. Without agitation, regardless of viscous interaction, the probabilities reach approximately 0.075 at $t=1.5$ s (cases A and D). With agitation, viscous interaction makes a difference to the probability evolution. Without viscous interaction and spatially uncorrelated agitation (case B), the probability reaches approximately 0.15 at $t=1.5$ s. With the use of spatially correlated agitation (case C), the probability becomes slightly larger at all $t$. With viscous interaction and spatially uncorrelated agitation (case E), the probability reaches approximately 0.3 at $t=1.5$. The spatial correlation of the agitation likewise increases the probability with a slight magnitude (case F). Overall, compared with the cases without agitation, the probability is increased by a factor of 2 by modelling the agitation, and by a factor of 4 by modelling both the agitation and viscous interaction, regardless of the spatial correlation of the agitation, at all $t$. These results indicate that the combination of viscous interaction and turbulence can significantly enhance the side-by-side clustering of initially in-line pairs, at the time scale of the dilute cluster formation in the experiment.

Figure 18. Evolution of the probability of side-by-side clustering in the case that the trailing bubble is initially released in the wake region, obtained from the sets shown in table 1.

6.3. Implication for the formation of larger clusters

Viscous wakes and potential interactions are common for spherical bubbles rising at $Re>1$, and the symmetry breaking of the in-line configurations as well as side-by-side clustering of bubbles that were identified in the previous sections could be observed in various regimes of bubbly flows. The potential-induced attraction that leads to the formation of clusters perpendicular to gravity was predicted by previous models of three-dimensional bubbly dispersion in potential flows (Sangani & Didwania Reference Sangani and Didwania1993; Smereka Reference Smereka1993; Yurkovetsky & Brady Reference Yurkovetsky and Brady1996; Spelt & Sangani Reference Spelt and Sangani1998). Such strong clustering has not been observed in experiments. The discrepancy between the theory and experiments has been associated with fluctuations induced by bubbles that are not included in these theories, including path instability due to deformations of bubbles and the wake-induced disturbances in the flow (Risso Reference Risso2018). On the contrary, the present results suggest that, for bubbles trapped near the wall, the viscous wake may accelerate the clustering of bubbles. Bouche et al. (Reference Bouche, Roig, Risso and Billet2012) pointed out that the viscous wake is largely attenuated when bubbles rise in a thin gap. In the present experiment, the depth of the channel is much greater than the size of bubbles and bubbles are not subjected to such constraints. Given the agreement between the experiment and the model, we consider that the effects of the channel geometry may not drastically influence the physics of pairwise clustering identified in the previous sections, within the validity of the model.

At a bulk void fraction of $O(1)$ % or greater, denser and larger clusters can appear, compared with the chain-like clusters shown in figure 16. These clusters rise faster than isolated bubbles, and catch up with leading bubbles to merge and grow vertically as well as breakup (Takagi & Matsumoto Reference Takagi and Matsumoto2011). Detailed mechanisms of inter-bubble interactions that lead to this dynamics remain uncovered. In this regime, bubble clusters not only interact with the background turbulent flow field but also can alter the flow structures through two-way coupling. The present modelling assumption may not hold that the modelled agitation is homogeneous and isotropic on the plane where bubbles reside. Nevertheless, during the onset of clustering from a state of sparse random dispersion, the present approach could still be of use to predict the time scale of clustering.