Sage Journals: Discover world-class research

Abstract

Bubble–bubble interactions play important roles in the dynamic behaviours of multiple bubbles or bubble clouds in a vortical flow field. Based on the Rayleigh–Plesset equation and the modified Maxey–Riley equation of a single bubble, bubble–bubble interaction terms are derived and introduced for multiple bubbles. Thus, both the Rayleigh–Plesset and modified Maxey–Riley equations are improved by considering bubble–bubble interactions and then applied for the multiple bubbles entrainment into a stationary Gaussian vortex. Runge–Kutta fourth-order scheme is adopted to solve the coupled dynamic and kinematic equations and the convergence study has been conducted. Numerical result has also been compared and validated with the published experimental data. On this basis, the oscillation, trajectory and effects of different parameters of double-bubble and multi-bubble entrainment into Gaussian vortex have been studied and the results have been compared with those of the cases without bubble–bubble interactions. It indicates that bubble–bubble interactions influence the amplitudes and periods of bubble oscillations severely, but have small effects on bubble trajectories.

Keywords

Bubble–bubble interactions multiple bubbles bubbles entrainment Gaussian vortex

Introduction

Cavitation and gas bubbles entrainment into a vortical flow field is always of great importance, because of its either useful applications or undesirable effects, in chemical engineering,¹ hydraulic engineering,² marine and ocean engineering^3,4 and so on. One famous example in marine and ocean engineering field is propeller cavitation.⁵ It involves the kinematics and dynamics of multiple cavitation bubbles emerged at the finite hydrofoils. The vortex is generated from the interaction between the rotating propeller and the flow. The bubble nuclei translate and rotate along the vortex as well as oscillate and interact with each other. When the pressure at the bubble position induced by the vortex is lower than the local saturated vapour pressure, the bubble nucleus will grow explosively there. This process is usually called ‘cavitation inception’. Many studies^5,6 indicate that the vortex cavitation bubble clouds mentioned above are one of the main reasons for the white wake after the stern of a high-speed ship.

To reveal this complex phenomenon, researchers have made great efforts in experimental, analytical and numerical simulations. Main studies in the last century can be found in the review paper ‘Cavitation in vortical flows’ by Arndt.⁷ From the viewpoint of numerical modelling, work on cavitation and gas bubble entrainment into a vortical flow field can be broadly divided into two categories. In the first category, a small bubble nucleus is released and will be captured in the vortical flow field. Because just one single bubble is considered, it is relatively easy to track its motion and deformation induced by the vortex flow. Different methods were adopted to study the entrainment of the single bubble depending on whether the bubble was spherical or not. Hsiao and Pauley⁶ studied the cavitation inception in a tip vortex flow through the spherical Rayleigh–Plesset (RP) bubble dynamics equation coupled with the bubble motion equation. The tip vortex flow was obtained from a Reynolds-averaged Navier–Stokes (RANS) solver. On this basis, Hsiao et al.⁸ improved the spherical bubble dynamic model using the average pressure over the bubble surface to replace the pressure at the bubble centre. However, the bubble–vortex, vortex–vortex interaction and the non-spherical deformation of the bubble had not been considered in these researches. Hsiao and Chahine⁹ applied spherical and non-spherical bubble dynamic models to study the cavitation inception, deformation and the scaling of the bubble in a vortex flow. The spherical and non-spherical models were both modified, where the spherical model accounted for the slip velocity between the bubble and the background flow, while the non-spherical model was embedded in an unsteady RANS code. Hsiao and Chahine¹⁰ further extended this model to study the cavitation inception in a ducted propulsor and the vortex–vortex interaction was considered using Navier–Stokes (NS) computations. The work by Hsiao and Chahine et al. provided much insight on the single cavitation bubble captured by the vortex. Oweis et al.¹¹ applied direct numerical simulations and point-particle tracking method to study capture and inception of bubbles near a line vortex under non-cavitating and cavitating conditions. Numerical results were compared and validated with their experimental data from Oweis et al.¹² Choi et al.¹³ applied both axisymmetric and three-dimensional numerical models to examine the non-spherical growth, oscillation and collapse of a vortex cavitation bubble from a spherical nucleus. The numerical models, named as ‘DF_UNCLE’, were based on finite volume formulation to solve the NS equations including full bubble–vortex interaction. Numerical results were compared with the experimental study of Choi and Ceccio¹⁴ who recorded the deformation of individual cavitation vortex bubble produced by a laser-initiated nuclei passing through a region of pressure reduction and recovery. Ni et al.¹⁵ and Zhang and Ni¹⁶ used one-way coupled particle tracking method and the boundary element method to simulate the entire process of the motion of a bubble in a Gaussian line vortex field, including its release, capture, inception, deformation and splitting.

In the second category of the problem, multiple bubbles or bubble clouds are initially released in the vortical flow. The growth of the bubble number raises some challenging problems which a single bubble has not encountered. The most distinguishing one is the bubble will be influenced by the neighbouring bubbles, namely bubble–bubble interactions (BBI), which is short for BBI in this article. However, most previous researchers studied the multiple bubbles entrainment into vortex flow using the same method as the single bubble and just increased the number of equations directly. For example, Hsiao and Chahine¹⁷ simulated the acoustic pressure generated by cavitation inception in a tip vortex flow in water containing a realistic bubble nuclei size distribution through the spherical bubble dynamics model directly. Hsiao et al.⁵ considered the effect of gas diffusion across the bubble wall on the entrainment and dynamics of bubble nuclei population in the flow field of a propeller. Park et al.¹⁸ distributed multi-bubbles with prescribed size and position initially and studied the multi-bubble motions in a vortex flow of a hydrofoil by spherical RP model. It was found that cavitation inception of smaller nuclei is more sensitive to the change of cavitation number. Some other researchers focused on the bubble–vortex interactions rather than BBI. Finn et al.¹⁹ developed a discrete element model to simulate the bubble entrainment and interactions with two-dimensional vortical flows. It was found that the existence of the bubble resulted in the distortion of the vortex. Cihonski et al.^20,21 further extended the study to three-dimensional simulation of travelling vortex ring together with a few micro-bubbles using a volumetric coupling approach. However, in their model, they ignored bubble oscillations and just took bubbles as spheres with gas density.

On the other hand, BBI have been extensively studied in various flow fields, such as in infinite flow field,^22,23 near the walls,^24–26 in homogeneous bubbly flows,²⁷ in a sound field^28–31 and under given pressure waves.³² It is found that the coupling of bubbles may lead to a variety of phenomena that a single bubble can never exhibit, such as attraction and repulsion of pulsating bubbles,^24,28,33,34 phase delays³⁵ and even coalescence.^36,37 However, most previous studies on BBI in these flow fields focused on the oscillations of bubbles at their static balance positions, which denote the positions of stationary bubbles when the inner pressure equals to the fluid pressure outside, instead of the translations of bubbles, so the effect of BBI on the motion equations of moving bubbles has been seldom considered.

It seems that there has been little work which considers the BBI during multiple moving bubbles. On the other hand, inspired by the previous studies in other flow fields, BBI are expected to affect multiple bubbles entrainment into a vortical flow field, which is still not understood and needs to be studied further. These form the prime motivation of the present work and give its distinctive focus. Here, we adopt a stationary Gaussian vortex and distribute prescribed numbers of very small bubble nuclei inside it. RP equation and modified Maxey–Riley (MR) equation are both improved and solved simultaneously to obtain the oscillations and motions of bubbles. The numerical model is validated against the published experimental data. Influences of different parameters such as bubble radii and number are further considered.

The main motivation or contribution of this article from a modelling standpoint is the incorporation of additional terms in the RP equation and kinematic equation of bubble motion to account for BBI. It shows that the amplitude and period of bubble oscillation have been changed significantly relative to the un-modified equations. As a result, the vortex cavitation inception and motion of multi-bubble nuclei may be predicted more accurately by the model developed here, which would lead to a better way to control the propeller cavitation and improve the effectiveness of the propeller of high-speed vessels.

Mathematical formulation

Dynamic and kinematic models of multiple bubbles

Dynamic model

To study the dynamics of multiple bubbles, we first consider the oscillation of a spherical bubble in an infinite free flow field. It is assumed that the fluid is incompressible with constant density $ρ = 1000 kg / m^{3}$ and dynamic viscosity coefficient $μ = 1.0037 * 10^{- 3} Pa s$ . The NS equation in the spherical coordinate frame is³⁸

\frac{\partial u_{r}}{\partial t} + u_{r} \frac{\partial u_{r}}{\partial r} = - \frac{1}{ρ} \frac{\partial P}{\partial r} + υ (\frac{\partial^{2} u_{r}}{\partial r^{2}} + \frac{2}{r} \frac{\partial u_{r}}{\partial r} - \frac{2 u_{r}}{r^{2}})

(1)

where $u_{r}$ and P are the radial velocity and pressure, respectively; r is the radial coordinate to the bubble centre; t is time; and $υ$ is kinematic viscosity coefficient. According to the mass conservation equation of the incompressible fluid, $u_{r} (r, t) = \overset{\cdot}{R} (t) R (t)^{2} / r^{2}$ , where $R (t)$ is bubble radius, $\overset{\cdot}{R} (t)$ is the first-order derivative of $R (t)$ to time, namely, the normal velocity of the bubble surface. Substituting $u_{r}$ into equation (1), one has

- \frac{1}{ρ} \frac{\partial P}{\partial r} = \frac{\overset{\cdot\cdot}{R} R^{2} + 2 {\overset{\cdot}{R}}^{2} R}{r^{2}} - \frac{2 {\overset{\cdot}{R}}^{2} R^{4}}{r^{5}}

(2)

It should be noticed that the last term on the right-hand side of equation (1) becomes zero in equation (2). As discussed by Yang and Yeh,³⁹ the viscous term in Newton liquid to a spherical bubble is zero and the contribution of viscosity just comes from the boundary layer of the bubble, as shown in equation (4).

Integrating equation (2) from the bubble radius R to the infinity $\infty$ , one has

\frac{1}{ρ} (P_{l} - P_{\infty}) = \overset{\cdot\cdot}{R} R + \frac{3}{2} {\overset{\cdot}{R}}^{2}

(3)

where $P_{l}$ is the liquid pressure at the bubble surface and $P_{\infty}$ is the reference pressure at infinity. Equation (3) is just the famous RP equation.⁴⁰

Taking the surface tension and the viscous stress at the bubble surface into account, one can link the bubble inner pressure $P_{b}$ and $P_{l}$ as

P_{l} = P_{b} - \frac{4 μ \overset{\cdot}{R}}{R} - \frac{2 σ}{R} = P_{c} + P_{g 0} {(\frac{R_{0}}{R})}^{3 γ} - \frac{4 μ \overset{\cdot}{R}}{R} - \frac{2 σ}{R}

(4)

where $σ$ is the surface tension coefficient; the second equality is based on the assumption that the bubble contains liquid vapour and non-condensable gas, which follows a polytropic law. $P_{c}$ is the saturated vapour pressure of the condensable gas. $P_{g 0}$ and $R_{0}$ are the initial pressure and volume when the bubble is formed, $P_{g 0} = P_{l 0} + 2 σ / R_{0} - P_{c}$ for an initial stationary bubble with $P_{l 0}$ as the initial undisturbed liquid pressure. $γ$ is the ratio of specific heat of the gas. $σ$ , $P_{c}$ and $γ$ are taken as $0.072 N / m$ , $2338 Pa$ and 1.25, respectively, in this article.

Substituting equation (4) into equation (3) and considering the relative translational slip velocity between the bubble velocity $u_{b}$ and the carrying liquid velocity $u_{w}$ , one has⁴¹

R \overset{\cdot\cdot}{R} + \frac{3}{2} {\dot{R}}^{2} + 4 υ \frac{\dot{R}}{R} - \frac{{(u_{w} - u_{b})}^{2}}{4} = \frac{1}{ρ} [P_{c} + P_{g 0} {(\frac{R_{0}}{R})}^{3 γ} - \frac{2 σ}{R} - P_{w}]

(5)

where $P_{w}$ is the undisturbed vortex pressure and the last term on the left-hand side of equation (5) is due to the effect of relative translational slip velocity.⁸

Equation (5) can be used to study the oscillation of a single bubble in a carrying flow. However, when multiple bubbles are involved, BBI should be incorporated on the basis of equation (5). It is assumed any bubble is so small that the asphericity of bubble surface can be ignored. On the other hand, the bubble size is assumed to be much smaller than the distance between bubbles and the existence of a bubble will not change the sphericity of another one. By integrating equation (2), one can obtain the pressure wave emitted by a pulsating bubble j as⁴²

P_{j} = ρ (\frac{{\overset{\cdot\cdot}{R}}_{j} R_{j}^{2} + 2 {\overset{\cdot}{R}}_{j}^{2} R_{j}}{r_{j}} - \frac{{\overset{\cdot}{R}}_{j}^{2} R_{j}^{4}}{2 r_{j}^{4}}) = \frac{ρ}{r_{j}} \frac{d ({\overset{\cdot}{R}}_{j} R_{j}^{2})}{dt} + O (\frac{1}{r_{j}^{4}})

(6)

where $r_{j}$ is the distance measured from the centre of bubble j. The amplitude of the second term on the right-hand side of equation (6) is very small compared to the first one and can be neglected.⁴³

It is further assumed that bubble sizes and oscillation periods are much smaller than the length and time scale of the ambient flow, respectively. Thus, the presence and oscillation of the bubbles will not change the distribution of the ambient flow, which is also named as ‘one-way coupling’ Lagrangian point-particle model.²¹ The pressures induced by the neighbouring bubbles can be linearly added²² to the vortex flow pressure and thus the oscillation equation of bubble i can be obtained by

\begin{array}{l} R_{i} {\overset{\cdot\cdot}{R}}_{i} + \frac{3}{2} {\dot{R}}_{i}^{2} + 4 υ \frac{{\dot{R}}_{i}}{R_{i}} + \sum_{j = 1, j \neq i}^{N} \frac{1}{D_{i j}} \frac{d ({\dot{R}}_{j} R_{j}^{2})}{d t} \\ - \frac{{(u_{w} - u_{b i})}^{2}}{4} = \frac{1}{ρ} [P_{c} + P_{g 0 i} {(\frac{R_{0 i}}{R_{i}})}^{3 γ} - \frac{2 σ}{R_{i}} - P_{w}] \end{array}

(7)

where $D_{ij}$ is the distance between the centres of bubbles i and j. Equation (7) can also be derived using the perturbation theory and Bernoulli equation and the details can refer to Liang et al.³¹

Kinematic model

Based on the MR equation⁴⁴ of a small rigid sphere in a non-uniform flow, the motion equation of a small spherical bubble in an incompressible liquid in the Lagrangian reference frame can be written as

m_{b} \frac{d u_{b}}{dt} = F_{a} + F_{g} + F_{d} + F_{l} + F_{p}

(8)

where $m_{b} = ρ_{b} 4 π R^{3} / 3$ is the mass of the gas bubble with the gas density $ρ_{b}$ , added mass force $F_{a} = (2 / 3) π R^{3} ρ ((D u_{w} / Dt) - (d u_{b} / dt)) + 2 π R^{2} (dR / dt) ρ (u_{w} - u_{b})$ , net buoyancy force $F_{g} = - (4 / 3) π R^{3} (ρ - ρ_{b}) g$ , drag force $F_{d} = (1 / 2) ρ π R^{2} C_{D} (u_{w} - u_{b}) | u_{w} - u_{b} |$ , lift force $F_{l} = (1 / 2) ρ π R^{3} C_{L} (u_{w} - u_{b}) \times ω / χ$ and pressure gradient force $F_{p} = - (4 / 3) π R^{3} \nabla P_{w}$ , in which g is the gravity acceleration, $C_{D}$ and $C_{L}$ are the drag and lift coefficients, respectively, $χ = | ω | R / | u_{w} - u_{b} |$ is the dimensionless shear rate and $ω$ is the fluid vorticity of the undisturbed flow at bubble centre. In this article, $C_{D} = 24 (1 + 0.15 {Re}_{b}^{0.687}) / R e_{b}$ ⁴⁵ and $C_{L} = 0.59 χ^{0.25}$ ,⁴⁶ where the bubble Reynolds number is defined as $R e_{b} = 2 R | u_{w} - u_{b} | / υ$ . For the detailed derivation and discussion on these forces, see the study by Zhang and Ni.¹⁶

When multiple bubbles are considered, it is expected that the neighbouring bubbles will induce additional forces besides the forces above. The induced force $F_{ij}$ of bubble j to bubble i can be approximated by integral of induced pressure $P_{j}$ over the surface $S_{i}$

F_{ij} = - \int_{S_{i}} P_{j} n_{i} dS = - \int_{V_{i}} \nabla P_{j} dS \approx - V_{i} \nabla P_{j}

(9)

where $n_{i}$ is the normal vector of the interface of bubble i, pointing towards the liquid. Substituting equation (6) (or equation (2) directly) into equation (9) and omitting the higher order small amount, one has

\begin{array}{l} F_{i j} = - V_{i} e_{r} {\frac{\partial P_{j}}{\partial r} |}_{r = D_{i j}} \approx ρ \frac{4}{3} π R_{i}^{3} \frac{d ({\dot{R}}_{j} R_{j}^{2})}{d t} \frac{e_{r}}{D_{i j}^{2}} \\ = ρ \frac{4}{3} π R_{i}^{3} \frac{d ({\dot{R}}_{j} R_{j}^{2})}{d t} \frac{x_{i} - x_{j}}{{| x_{i} - x_{j} |}^{3}} \end{array}

(10)

where $e_{r}$ denotes the radial unit vector, pointing from the centre of bubble j to that of bubble i, and $x_{i}$ and $x_{j}$ are the coordinates of the centres of bubble i and j, respectively. $F_{ij}$ in equation (10) is also named as ‘secondary Bjerknes force’, while the pressure gradient force $F_{p}$ in equation (8) is also named as ‘primary Bjerknes force’.^29,47

Substituting equation (10) into equation (8) and considering gas density $ρ_{b} \approx ρ / 820$ , one can obtain the simplified bubble motion differential equation as

\begin{matrix} \frac{d u_{bi}}{dt} = & - g - \frac{3}{ρ} \nabla P_{w} + \frac{3}{4} \frac{C_{D}}{R_{i}} (u_{w} - u_{bi}) | u_{w} - u_{bi} | \\ + \frac{3}{4} C_{L} \frac{(u_{w} - u_{bi}) \times ω}{χ} \\ + \frac{3}{R} (u_{w} - u_{bi}) {\overset{\cdot}{R}}_{i} + \sum_{j = 1, j \neq i}^{N} 2 \frac{d ({\overset{\cdot}{R}}_{j} {R_{j}}^{2})}{dt} \frac{x_{i} - x_{j}}{{| x_{i} - x_{j} |}^{3}} \end{matrix}

(11)

Equation (11) can be adopted to calculate the translation of multiple bubbles in a given vortex flow. The Runge–Kutta fourth-order scheme is applied to solve equations (7) and (11) simultaneously. It should be noticed that although the induced velocity of bubble j may disturb the motion of bubble i also, this disturbance is $O [(R / D_{ij})^{2}]$ times smaller and can be neglected as discussed by Seo et al.²⁷

Gaussian vortex model

According to the hydrofoil theory, a finite aspect ratio hydrofoil generates a free vortex along the tip. Various free vortex models can be adopted corresponding to different kinds of foil shape, such as Gaussian vortex (also known as Lamb–Oseen vortex), Rankine vortex and Scully vortex. Gaussian vortex is adopted in this article.

Combining the Gaussian vortex flow and the incoming free-stream flow, one can obtain the velocity and pressure of the undisturbed carrying flow as

u_{w} = {\begin{matrix} u_{wr} = 0 \\ u_{w θ} = \frac{Γ}{2 π r'} (1 - e^{- ξ}) \\ u_{wz} = u_{\infty} \end{matrix}

(12)

\begin{array}{l} P_{w} = P_{\infty} + \int_{\infty}^{0} - \frac{ρ u_{w θ}^{2}}{r^{'}} d r^{'} = P_{\infty} - \frac{α ρ}{2} {(\frac{Γ}{2 π a_{c}})}^{2} \\ {ξ^{- 1} {(1 - e^{- ξ})}^{2} + 2 \int_{ξ}^{2 ξ} \frac{e^{- x}}{x} d x} \end{array}

(13)

where $r'$ is the radial distance from the measure point to the vortex line core; $Γ$ and $a_{c}$ are local vortex circulation and local core radius, respectively, depending on the foil shape; $u_{\infty}$ is the inflow velocity from infinity; $P_{\infty}$ is infinite hydrostatic pressure; and $ξ = α (r' / a_{c})^{2}$ with $α = 1.256$ . Corresponding vorticity and pressure gradient are

ω (r') = \nabla \times u_{w} = \frac{α Γ}{π a_{c}^{2}} e^{- ξ}

(14)

\nabla P_{w} = \frac{\partial P_{w} (r')}{\partial r'} = \frac{ρ u_{w θ}^{2}}{r'} = \frac{ρ Γ^{2}}{4 π^{2} {r'}^{3}} {(1 - e^{- ξ})}^{2}

(15)

It can be calculated that the maximum tangential velocity emerges at local core radius $a_{c}$ as $u_{w θ, max} = β Γ / 2 π a_{c}$ with $β = 1.256$ , while the minimum pressure emerges at vortex centreline as $P_{w, min} = P_{\infty} - η ρ (Γ / 2 π a_{c})^{2}$ with $η = 0.87$ . If the minimum pressure $P_{w, min}$ is smaller than saturated vapour pressure $P_{c}$ , the bubble is expected to expand at the vortex core explosively, which is also known as ‘cavitation inception’.⁴⁸ Here, we define the cavitation number as $σ_{\infty} = (P_{\infty} - P_{c}) / (0.5 ρ u_{\infty}^{2})$ . So cavitation inception will happen in theory when $σ_{\infty} ⩽ σ_{\infty, cri} = η Γ^{2} / 2 π^{2} a_{c}^{2} u_{\infty}^{2}$ .

Numerical results and discussions

Convergence study and comparison

Convergence study with time step

The present model is first verified through convergence study with time step. The parameters of the Gaussian vortex are same as those in Oweis et al.,¹² that is, $Γ = 0.290 m^{2} / s$ , $a_{c} = 5.6 mm$ , $u_{\infty} = 10 m / s$ and $σ_{\infty} = 1.4$ . One can predict that cavitation inception will not happen under this cavitation number based on $σ_{\infty, cri}$ defined above. In the flow field, the cylindrical coordinate system is established with the origin located at the local core of Gaussian vortex, as shown in Figure 1. In this section, the case of double bubbles named as bubble 1 and bubble 2 is considered. Their initial radii are $R_{10} = R_{20} = a_{c} / 50$ and initial coordinates are $r_{1} = 3 a_{c}$ , $θ_{1} = 0$ , $z_{1} = 0$ and $r_{2} = 2.6 a_{c}$ , $θ_{2} = 0.133 rad$ , $z_{2} = 0.4 a_{c}$ , respectively.

Figure 1.

Schematic of the problem.

The fixed time step $Δ t$ is decreased from $1 \times 10^{- 6} s$ to $1 \times 10^{- 7} s$ and then to $1 \times 10^{- 8} s$ . The radius oscillation of bubble 1 at different time steps is provided in Figure 2, where the time range is between 0.024 and 0.026 s for a clear view. It can be seen that the results between $1 \times 10^{- 7} s$ and $1 \times 10^{- 8} s$ are very close to each other, which shows good convergence. The results of bubble 2 are similar to those of bubble 1, which are not plotted here. In the following calculations, $Δ t$ is taken as $1 \times 10^{- 7} s$ .

Figure 2.

Radius oscillation of bubble 1 at different time steps.

Comparison with experimental data

The numerical model is further validated by comparison with experimental data. The experimental data on bubble rotating and captured by a vortex are rare, especially for multi-bubbles. Here, we adopted the photographic experimental data¹¹ of an initially static bubble releasing and moving in a line Gaussian vortex. The experimental parameters are $Γ = 0.290 m^{2} / s$ , $a_{c} = 5.6 mm$ , $u_{\infty} = 10 m / s$ , $σ_{\infty} = 3$ , $R_{0} = a_{c} / 6$ and $r = a_{c}$ . Numerical simulation adopts the same parameters and the initial pressure inside the bubble equals to the fluid pressure outside.

Figure 3 provides the comparison between numerical results and experimental data, where the circles come from the outlines of the bubble recorded in the photos by Oweis et al.¹¹ It can be seen that the position and volume of the numerical bubble agree with those of experiment in general. Because the numerical model is based on the spherical bubble theory, it can be understood that bubble shapes disagree. Besides, the capture time in experiment is 2.8 ms while the numerical one is 2.6 ms, which are close too. For more discussions, refer to Ni et al.¹⁵ The agreement validates the single-bubble model in the Gaussian vortex flow and the coupling effects of multi-bubbles will be studied on this basis.

Figure 3.

The comparison of bubble entrainment between experimental data (circles, from Oweis et al.¹¹) and numerical simulation (spheres).

Two-bubble case study

Based on the convergence study and validation with experimental data, the cases of double bubbles and multiple bubbles will be studied in this and next sections. Hereafter, the initial bubble radius ranged from $0.016 a_{c}$ to $0.024 a_{c}$ in our simulation, and the maximum bubble radius is no more than 1.6 times of the initial one, which are under ‘small bubble’ assumption of the ‘one-way coupling’ method.

Comparison between cases with and without BBI

The trajectory and oscillation of two bubbles are simulated in this section and the results are compared with those without BBI. The parameters are same with those in section ‘Convergence study with time step’.

Figures 4 and 5 show the radius oscillations of bubbles 1 and 2 in the cases with and without BBI, respectively. It is interesting to notice that the amplitude of bubble 1 is depressed while that of bubble 2 is intensified, compared to those without BBI. The average amplitude of bubble 1 is depressed by up to 49.0% while that of bubble 2 is intensified by up to 25.5%. As the bubble radii oscillate periodically, the amplitudes of them also show cyclical fluctuations. On the other hand, the oscillation periods of bubbles 1 and 2 increase by 0.6% and 2.2%, respectively, relative to the case without BBI. So BBI influence both the oscillation amplitude and period of bubbles distinctly. The changes are different from those of two static bubbles oscillating in a uniform pressure field, where two bubble oscillations are the same.³² Thus, we predict the pressure background or the bubble position affects the oscillation also. Actually, we have done a series of study by changing the initial relative positions of the two bubbles and found that the initial radial positions r affects the oscillation amplitude predominantly. The pressure varies with r, which induces the oscillation amplitude at a smaller r intensified but the one at a larger r depressed. However, it can be seen in Figures 4 and 5 that the bubble capture time and balance radius are almost same with those without BBI. It denotes that BBI may have little influence on the bubble trajectories, which will be validated in Figure 6.

Figure 4.

Radius oscillations of bubble 1 in cases with and without BBI: (a) within a long time and (b) within a short time.

Figure 5.

Radius oscillations of bubble 2 in cases with and without BBI: (a) within a long time and (b) within a short time.

Figure 6.

Trajectories of bubble 1 in the cases with and without BBI.

Figure 6 shows the comparison of trajectories of bubble 1 in the cases with and without BBI. It can be seen that bubble trajectory with BBI almost coincides with that of the bubble without BBI. This is because the BBI term $\sum_{j = 1, j \neq i}^{N} 2 (d ({\overset{\cdot}{R}}_{j} R_{j}^{2}) / dt) | x_{i} - x_{j} | / {| x_{i} - x_{j} |}^{3}$ and oscillation–motion interaction term $(3 / R_{i}) (u_{w} - u_{bi}) {\overset{\cdot}{R}}_{i}$ are both quite small compared with other forces in equation (11), especially with the drag and pressure gradient, which induces that BBI have little influence on the bubble trajectories.

Figure 7 shows trajectories and distance of two bubbles with BBI, where bubble positions corresponding to the largest distance in Figure 7(b) are marked in Figure 7(a). It can be seen that either bubble presents a kind of ‘reduced spiral’ motion. The largest distance happens when the two bubbles move to the two ends of the circumference relatively in their trajectories. When they approach to the vortex core or when they are captured, the two bubbles will move along the axis of the Gaussian vortex at carrying liquid velocity, so the distance between them remains constant as shown in Figure 7(b).

Figure 7.

Motion of bubbles with BBI: (a) trajectories of bubbles and (b) distance of bubbles.

Here, we will further investigate the influence of different parameters on BBI of double bubbles in a Gaussian vortex field, such as bubble radii and initial distance between two bubbles.

Initial radii of bubbles

The initial radii of the bubbles are important factors that influence bubble oscillations, as Li et al.³² showed in their work where the oscillation of double bubbles is driven by pressure fluctuation. As mentioned in section ‘Comparison with experimental data’, the initial radial position r affects bubble oscillations predominantly. In order to eliminate the influence of r, only z-coordinates of two bubbles are different here. The initial coordinates of the two bubbles are $r_{1} = 3 a_{c}$ , $θ_{1} = 0$ , $z_{1} = 0$ and $r_{2} = 3 a_{c}$ , $θ_{2} = 0$ , $z_{2} = 0.4 a_{c}$ , respectively. On this basis, we will keep the initial radius of bubble 1 as $R_{10} = a_{c} / 50$ constant and change the initial radius of bubble 2 from $0.8 R_{10}$ to $R_{10}$ and then to $1.2 R_{10}$ . The oscillation of bubble 1 will be studied.

Figure 8 shows the oscillation of bubble 1 at different ratios between the initial radii of the bubbles, where just a short time span is provided for a clearer view. It can be seen that initial radii ratio will influence the oscillation amplitude and period of the bubbles. For the oscillation period, it is observed that the period in the case of $R_{10} : R_{20} = 1 : 1$ is longest, which is $4.569 \times 10^{- 5} s$ , while the periods in the cases of $R_{10} : R_{20} = 1 : 0.8$ and $R_{10} : R_{20} = 1 : 1.2$ decrease by 2.4% and 3.2%, respectively. In terms of amplitude in Figure 8(b), it can be seen that the average amplitude in the case of $R_{10} : R_{20} = 1 : 0.8$ increases by up to 4.4% and that in the case of $R_{10} : R_{20} = 1 : 1.2$ decreases by up to 21.8% relative to the case of $R_{10} : R_{20} = 1 : 1$ . If the initial radius of the adjacent bubble decreases, the oscillation amplitude of the other bubble will be intensified periodically and vice versa.

Figure 8.

Oscillations of bubble 1 at different ratios of initial radii: (a) radius and (b) radius amplitude.

Initial distance between bubbles

The initial distance $D_{12}$ between the two bubbles will influence the strength of BBI. The initial bubble radii are taken as $R_{10} = R_{20} = a_{c} / 50$ and the initial coordinates of them are same as given in section ‘Initial radii of bubbles’.

Figure 9 shows the oscillation of the radius of bubble 1 with BBI at different $D_{12}$ . The period decreases from $4.569 \times 10^{- 5} s$ to $4.510 \times 10^{- 5} s$ and then to $4.491 \times 10^{- 5} s$ as $D_{12}$ increases from $0.4 a_{c}$ to $0.8 a_{c}$ and then to $1.2 a_{c}$ . If $D_{12}$ continues to increase, the period will approach to that in the case without BBI soon. As $D_{12}$ gets larger, the oscillation amplitude will get smaller, but the influence of the distance on amplitude is small relatively, with maximum rate of change only 0.8% and 1.2%, respectively.

Figure 9.

Oscillations of bubble 1 at different initial distances between two bubbles: (a) radius and (b) radius amplitude.

Multi-bubble case study

Based on the double-bubble case considered in previous sections, multi-bubble cases with different bubble numbers or distributions will be analysed in this section, which is more general in realistic situation.

Bubble number

Under different numbers of bubbles in the vortex field, bubble oscillations will change greatly. Different numbers of bubbles with same radii $a_{c} / 50$ are distributed along z axis, with $r = 3 a_{c}$ and $θ = 0$ . Their initial z-coordinates are $z_{1} = 0$ , $z_{2} = 0.4 a_{c}$ , $z_{3} = - 0.4 a_{c}$ , $z_{4} = 0.8 a_{c}$ and $z_{5} = - 0.8 a_{c}$ , respectively. For different cases, the corresponding bubble numbers will be selected.

Figure 10 shows radius oscillation of bubble 1 in multi-bubble case at different bubble numbers. It is observed that oscillation period of bubble 1 will increase by 1.6% and 3.8%, respectively, as bubble number rises from 2 to 3 and 5. Besides, the average amplitude of bubble 1 will also rise by up to 12.5% and by up to 16.5% correspondingly, respectively. It can be further seen in Figure 10(b) that as bubble number increases, the cyclical fluctuation of the oscillation amplitude will become more intense, with larger amplitude and smaller period.

Figure 10.

Oscillations of bubble 1 at different bubble numbers: (a) radius and (b) radius amplitude.

Bubble distribution

Three kinds of distribution of five bubbles with same radii $a_{c} / 50$ will be considered. They are distributed along r axis, z axis and the circle, respectively. In the first scenario, the initial radial coordinates of these five bubbles are $r_{1} = 3 a_{c}$ , $r_{2} = 2.6 a_{c}$ , $r_{3} = 2.2 a_{c}$ , $r_{4} = 3.4 a_{c}$ and $r_{5} = 3.8 a_{c}$ , respectively. In the second scenario, the initial axial coordinates are $z_{1} = 0$ , $z_{2} = 0.4 a_{c}$ , $z_{3} = 0.8 a_{c}$ , $z_{4} = - 0.4 a_{c}$ and $z_{5} = - 0.8 a_{c}$ . In the third scenario, the initial circumferential coordinates are $θ_{1} = 0$ , $θ_{2} = 0.4 π$ , $θ_{3} = 0.8 π$ , $θ_{4} = 1.2 π$ and $θ_{5} = 1.6 π$ , respectively.

Figure 11 shows the trajectories of these five bubbles at three kinds of distributions. It can be seen that all bubbles present kinds of ‘reduced spiral’ motion as mentioned in section ‘Comparison between cases with and without BBI’. As for the Figure 11(a), the bubble closer to the vortex core will be captured earlier. When bubble 1 and bubble 2 both approach to the vortex core, the distance between them becomes very small and the mathematical model does not work anymore. Thus, the calculation terminates then. Similar situation also happens for the case of distribution along the circle, where the calculation terminates when the distance between either two bubbles is smaller than a critical value, as shown in Figure 11(c).

Figure 11.

Bubble trajectories at three kinds of distributions: (a) r axis distribution, (b) z axis distribution and (c) circle distribution.

Conclusion

The dynamic behaviours of multiple bubbles entrainment into a vortical flow field have been investigated by solving the coupled improved RP and MR equations where BBI terms are introduced. On the basis of convergence study and comparison with experimental data, series of numerical simulation have been conducted. The results have also been compared with those without BBI. Some conclusions have been drawn as follows:

The oscillations of the bubbles are affected by BBI greatly, but the trajectories of the bubbles are hardly influenced.

When the bubbles have same initial radii, BBI is able to increase amplitude of the bubble with a smaller radial position by up to 25.5% while reduce that of the other one by up to 49.0% relative to the case without BBI over the parameter space studied. Meanwhile, it extends both the oscillation periods of the two bubbles by up to 2.2%.

When the bubbles have same radial positions, BBI is able to intensify the amplitude of the bigger bubble and depress that of the smaller bubble. In terms of the oscillation period, when the radii of two bubbles equal, the oscillation period reaches longest.

Both amplitudes and periods of the bubbles will rise with the increase of bubble number.

Footnotes

Academic Editor: Moran Wang

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by Lloyd’s Register Foundation through the joint centre involving University College London, Shanghai Jiaotong University and Harbin Engineering University, the National Natural Science Foundation of China (Nos 11302056, 11302057, 11472088 and 51579054) and the International Postdoctoral Exchange Fellowship Program (No. 20140068), to which the authors are most grateful. Lloyd’s Register Foundation helps to protect life and property by supporting engineering-related education, public engagement and the application of research.

References

Aoki

Hoque

Chanson

Physical modelling and similitude of air bubble entrainment at vertical circular plunging jets. Chem Eng Sci 2004; 56: 747–758.

Hoschek

Carrica

Weber

Bubble entrainment and distribution in a model spillway with application to total dissolved gas minimization. J Hydraul Eng: ASCE 2008; 134: 763–771.

Zhang

Liu

YL.

Improved three-dimensional bubble dynamics model based on boundary element method. J Comput Phys 2015; 294: 208–223.

Zhang

Cui

. Experimental study on bubble dynamics subject to buoyancy. J Fluid Mech 2015; 776: 137–160.

Hsiao

Jain

Chahine

. Effect of gas diffusion on bubble entrainment and dynamics around a propeller. In: Proceedings of the 26th symposium on naval hydrodynamics, Roman, 17–22 September 2006. Washington, DC: The National Academies Press.

Hsiao

Pauley

LL.

Study of vortex cavitation inception using Navier–Stokes computation and bubble dynamics model. J Fluid Eng: T ASME 1999; 121: 198–204.

Arndt

REA

. Cavitation in vortical flows. Annu Rev Fluid Mech 2002; 34: 143–175.

Hsiao

Chahine

Liu

HL.

Scaling effect on prediction of cavitation inception in a line vortex flow. J Fluid Eng: T ASME 2003; 125: 53–60.

Hsiao

Chahine

GL.

Prediction of tip vortex cavitation inception using coupled spherical and nonspherical bubble models and Navier–Stokes computations. J Mar Sci Technol 2004; 8: 99–108.

10.

Hsiao

Chahine

GL.

Numerical study of cavitation inception due to vortex/vortex interaction in a ducted propulsor. J Ship Res 2008; 52: 114–123.

11.

Oweis

Hout

Lyer

. Capture and inception of bubbles near line vortices. Phys Fluids 2005; 17: 105–118.

12.

Oweis

Choi

Ceccio

SL.

Dynamics and noise emission of laser induced cavitation bubbles in a vortical flow field. J Acoust Soc Am 2004; 115: 1049–1058.

13.

Choi

Hsiao

Chahine

. Growth, oscillation and collapse of vortex cavitation bubbles. J Fluid Mech 2009; 624: 255–279.

14.

Choi

Ceccio

SL.

Dynamics and noise emission of vortex cavitation bubbles. J Fluid Mech 2007; 575: 1–26.

15.

Zhang

Yao

. Numerical simulation of trajectory and deformation of bubble in tip vortex. Appl Math Mech 2012; 33: 701–716.

16.

Zhang

BY.

Influences of different forces on the bubble entrainment into a stationary Gaussian vortex. Sci China G: Phys Mech 2013; 56: 2162–2169.

17.

Hsiao

Chahine

GL.

Scaling of tip vortex cavitation inception noise with a bubble dynamics model accounting for nuclei size distribution. J Fluid Eng: T ASME 2005; 127: 55–65.

18.

Park

Seol

Choi

. Numerical prediction of tip vortex cavitation behavior and noise considering nuclei size and distribution. Appl Acoust 2009; 70: 674–680.

19.

Finn

Sham

Apte

SV.

Modeling and simulation of multiple bubble entrainment and interactions with two dimensional vortical flows. Phys Fluids 2011; 23: 301–320.

20.

Cihonski

Finn

Apte

. Modeling and simulation of multiple bubble entrainment and interactions with a traveling vortex ring. In: Proceedings of ASME 2012 fluids engineering division summer meeting, Rio Grande, PR, 8–12 July 2012. New York: ASME.

21.

Cihonski

Finn

Apte

SV.

Volume displacement effects during bubble entrainment in a travelling vortex ring. J Fluid Mech 2013; 721: 225–267.

22.

Ida

Multibubble cavitation inception. Phys Fluids 2009; 21: 113302.

23.

Chew

Klaseboer

Ohl

. Interaction of two differently sized oscillating bubbles in a free field. Phys Rev E 2011; 84: 1–10.

24.

Quinto-su

Ohl

Interaction between two laser-induced cavitation bubbles in a quasi-two-dimensional geometry. J Fluid Mech 2009; 633: 425–435.

25.

Zhang

Song

. Numerical simulation of bubble breakup phenomena in a narrow flow field. Appl Math Mech 2010; 31: 449–460.

26.

Zhang

Cui

Study on splitting of a toroidal bubble near a rigid boundary. Phys Fluids 2015; 27: 062102.

27.

Seo

Lele

Tryggvason

Investigation and modeling of bubble-bubble interaction effect in homogeneous bubbly flows. Phys Fluids 2010; 22: 063302–063319.

28.

Mettin

Akhatov

Parlitz

. Bjerknes forces between small cavitation bubbles in a strong acoustic held. Phys Rev E 1997; 56: 2924–2931.

29.

Mettin

Luther

Ohl

. Acoustic cavitation structures and simulations by a particle model. Ultrason Sonochem 1999; 6: 25–29.

30.

Formulation of multibubble cavitation. Phys Rev E 2011; 83: 066313.

31.

Liang

Chen

Shao

. Aspherical oscillation of two interacting bubbles in an ultrasound field. Chinese Phys Lett 2012; 29: 074701.

32.

Cai

Huai

. Interaction mechanism of double bubbles in hydrodynamic cavitation. J Therm Sci 2003; 22: 242–249.

33.

Oguza

Prosperetti

A generalization of the impulse and virial theorems with an application to bubble oscillation. J Fluid Mech 1990; 218: 143–162.

34.

Harkin

Kaper

Nadim

Coupled pulsation and translation of two gas bubbles in a liquid. J Fluid Mech 2001; 445: 377–411.

35.

Ida

Phase properties and interaction force of acoustically interacting bubbles: a complementary study of the transition frequency. Phys Fluids 2005; 17: 097107.

36.

Rungsiyaphornrat

Klaseboer

Khoo

. The merging of two gaseous bubbles with an application to underwater explosions. Comput Fluids 2003; 32: 1049–1074.

37.

Zhang

Yang

Huang

. Numerical simulation of column charge underwater explosion based on SPH and BEM combination. Comput Fluids 2013; 71: 169–178.

38.

Lamb

Hydrodynamics. 6th ed.Cambridge: Cambridge University Press, 1932 (reprinted by Dover, 1945).

39.

Yang

Yeh

HC.

Theoretical study of bubble dynamics in purely viscous fluids. AIChE J 1963; 12: 927–931.

40.

Rayleigh

On the pressure developed in liquid during the collapse of a spherical cavity. IDEA Fitness J 1917; 34: 94–98.

41.

Plesset

Prosperetti

Bubble dynamics and cavitation. Annu Rev Fluid Mech 1977; 9: 145–185.

42.

Prosperetti

Bubble phenomena in sound fields: part two. Ultrasonics 1984; 22: 115–124.

43.

Barbat

Ashgriz

Liu

GS.

Dynamics of two interacting bubbles in an acoustic field. J Fluid Mech 1999; 389: 137–168.

44.

Maxey

Riley

JJ.

Equation of motion for a small rigid sphere in a nonuniform flow. Phys Fluids 1983; 26: 883–889.

45.

Schiller

Naumann

A drag coefficient correlation. VDI Zeitung 1935; 77: 51–59.

46.

Sridhar

Katz

Drag and lift forces on microscopic bubbles entrained by a vortex. Phys Fluids 1995; 7: 389–399.

47.

Leighton

Walton

Pickworth

MJW

. Primary Bjerknes forces. Eur J Phys 1990; 11: 47.

48.

Brennen

CE.

Cavitation and bubble dynamics. Oxford: Oxford University Press, 1995, p.30.

Bubble–bubble interaction effects on dynamics of multiple bubbles in a vortical flow field

Abstract

Keywords

Introduction

Mathematical formulation

Dynamic and kinematic models of multiple bubbles

Dynamic model

Kinematic model

Gaussian vortex model

Numerical results and discussions

Convergence study and comparison

Convergence study with time step

Comparison with experimental data

Two-bubble case study

Comparison between cases with and without BBI

Initial radii of bubbles

Initial distance between bubbles

Multi-bubble case study

Bubble number

Bubble distribution

Conclusion

Footnotes

Declaration of conflicting interests

Funding

References