Modification of flow transition for a vortex ring by a bubble released at the axis

Abstract

The interactions between a vortex ring and a gas bubble released at the axis are studied numerically, which shed light on understanding more complicated bubble–turbulence coupling. We fix the Reynolds number at $Re_{\tau }=7500$ and consider various Weber numbers in the range of $We=130\unicode{x2013}870$. We find that the translating speed of the vortex ring is substantially lower than the case of the vortex ring without a bubble. It is explained with two different mechanisms, depending on the Weber number. In the low-$We$ range, the reduction of translating speed of the ring is due to the capture of bubbles into the ring core, leading to significant changes in the vorticity distribution within the core. In the high-$We$ range, the repeatedly generated secondary vortex rings perturb the primary one, which bring about an earlier flow transition, thereby reducing the translating speed of the vortex ring. On the other hand, the evolution of a gas bubble is also affected by the presence of the vortex ring. In the low-$We$ range, we observe binary breakup of the bubble after it is captured by the primary vortex ring. In contrast, in the high-$We$ range, it is interesting to find that the bubble experiences sequentially stretching, spreading and breakup stages. In the high-$We$ range, the numbers of smaller bubbles predicted by the classical Rayleigh–Plateau instability of a stretched cylindrical bubble agree well with our numerical simulations. Consistent with the previous experiments, this number keeps unchanged at 16 as $We$ further increases. An additional comparison is made between two higher Reynolds numbers, indicating that the finer eddies in a vortical field with a higher Reynolds number tend to tear the bubble into more fragments.

1. Introduction

Bubbly flows occur in a wide range of applications, from chemical reactors to ship hydrodynamics. The interaction between turbulent eddies and bubbles determines the resulting bubble size distribution and, the other way around, the bubbles play an important role in the turbulent properties of the bubbly flows. It has been well established that the local flow pattern can determine the deformation and breakup of the bubbles, depending mainly on the Weber number, $We=\rho \epsilon ^{2/3}R^{5/3}_b/ \sigma$, representing the ratio between the turbulent stresses acting on the surface of the bubble and the confining stresses due to surface tension, where $\epsilon$ is the turbulent dissipation rate, $R_b$ is the bubble radius and $\sigma$ is the surface tension (Martínez-Bazán 2015). The early works of Kolmogorov (1949) and Hinze (1955) proposed two mechanisms for bubble breakup, as reviewed by Martínez-Bazán (2015). The inertial mechanism leads to bubble breakup when the Weber number is high, when the pressure fluctuations acting on the bubble surface are sufficiently large to overcome the confining surface tension forces (Martinez-Bazan, Montanes & Lasheras 1999a,b). The second mechanism, at low Weber number, is a resonance between the oscillation of the bubble and the turbulent flow, which leads to an accumulation of surface energy due to the bubble interaction with a sequence of turbulent eddies until the bubble breaks (Risso & Fabre 1998; Lauterborn & Kurz 2010). A more general view of the bubble deformation and breakup problem can be found in the review by Risso (2000), or a more recent review of the direct numerical simulation of turbulent flows laden with droplets and bubbles (Elghobashi 2019). Recently, experimental results were obtained regarding the shape of a bubble placed in a horizontal solid body rotational flow, reporting that a resonance between the driving frequency and the eigenfrequency of the bubble leads to a binary breakup of the bubble (Rodgar et al. 2021).

We know that the reduction of frictional drag through gas or bubble injection has been a hot topic (Ceccio 2010; Murai 2014), which leads to extensive studies of the interactions between bubbles and turbulent boundary layers. It has been reported that a few volume percent (${\leq }4\,\%$) of bubbles can reduce the overall drag by up to $40\,\%$ and beyond (Van den Berg et al. 2005; van den Berg et al. 2007; Elbing et al. 2008,  2013; van Gils et al. 2013). However, the exact physical mechanisms responsible for frictional drag reduction using bubbles in turbulent boundary layers are still not completely understood (Murai 2014). Generally speaking, drag reduction by bubbles can either be caused by direction modification of fluid properties like density and viscosity (L'vov et al. 2005), or through more complicated interactions between bubbles and turbulent structures within the boundary layer (Ferrante & Elghobashi 2004; Lu, Fernández & Tryggvason 2005). The latter has attracted more attention due to the rich physics in this two-way coupled system. Apparently, the coherent structures present in turbulent flow can be modified by the presence of bubbles (Ferrante & Elghobashi 2004), and on the other hand, the strong shear in the boundary layer can deform the bubbles which in turn modifies the turbulent flow (Spandan, Verzicco & Lohse 2018). Lu et al. (2005) reported that the presence of bubbles could suppress streamwise vortices and enstrophy, therefore reducing the wall shear stress. It has been shown experimentally in turbulent Taylor–Couette flows that the bubble deformability is crucial in the drag reduction mechanism, which modifies the lift force coefficient of the lift acting on the bubble and thus the bubble distribution in the flow (van Gils et al. 2013; Verschoof et al. 2016), which has also been confirmed by numerical simulations (Spandan et al. 2018).

There have been relatively few investigations of the interaction of bubbles with a single vortex structure. Jha & Govardhan (2015) studied experimentally the interaction of a single bubble with a single vortex ring in water, covering a range of Weber numbers between 0.155 and 19 according to the definition $We_{\tau }=\rho _l (\varGamma _0/2 {\rm \pi}R_0)^2 / (\sigma / R_b )$. Here, $\rho _l$ is the water density, $\varGamma _0$ is the initial vortex ring circulation, $R_0$ is the ring radius, $R_b$ is the bubble radius and $\sigma$ is the surface tension. Another definition, $We=0.87\rho _l(\varGamma _0/2{\rm \pi} a)^2/(\sigma /D_b)$ with $D_b=2R_b$ and $a$ being the vortex core radius, gives a range of 3 to 406. In their experiments, a gas bubble was initially released outside the vortex ring, and the bubble was smaller than the vortex ring. They found that the bubble was attracted towards the vortex due to the radial pressure gradient induced by the ring, and got engulfed into the vortex core. Once inside the core, the bubble was elongated due to the azimuthal pressure gradient caused by its presence, until it eventually broke. They also reported that, in the low-Weber-number case (lower than a critical value), the vortex core was fragmented when the bubble was trapped, with a simultaneous reduction in enstrophy and ring convection velocity. Similar processes of a bubble being trapped into a vortex ring as well as its subsequent azimuthal elongation and breakup have also been observed in more recent research (Biswas & Govardhan 2020), with both a single bubble and multiple bubbles considered. With a larger bubble, which is comparable to the vortex ring, experiments were performed with a major focus on the bubble breakup parameters, such as breakup frequency, number of daughter bubbles resulting from the breakup and the size distribution (Zednikova et al. 2019).

Although the two-way coupling between bubble and vortex ring has been considered in the previous studies mentioned above (Jha & Govardhan 2015; Zednikova et al. 2019; Biswas & Govardhan 2020), their major focus was still on the deformation and breakup of bubbles. It would be of interest to study the effect of bubbles on the characteristic properties of the vortex ring to extrapolate the results to more complex turbulent bubbly flows, as addressed by Martínez-Bazán (2015). It is known that three-dimensional direct numerical simulations can give more details of the flow structure, particularly the breakup and reconnection of vortex filaments related to flow transition. In our present work, we study numerically the idealised case of the interaction of a single bubble with a single vortex ring, at a fixed Reynolds number $Re_\varGamma =7500$, with varying Weber numbers. We are interested in both the bubble dynamics perspective and the modified transition process induced in the vortex ring due to the interaction.

2. Formulation of the problem

2.1. Problem set-up

The initial flow field is a quasi-Gaussian or the Lamb–Oseen vortex (Fabre, Sipp & Jacquin 2006; Mao & Hussain 2017). Non-dimensionalising the lengths by the ring radius $R_0$, the vortex profile is centred at $(x, \eta )=(0, 1)$, and its non-dimensional velocity components in the $(x, \eta, \phi )$ coordinates (see figure 1) are written as

(2.1)\begin{gather} u_{x}=\frac{\varGamma_0}{2{\rm \pi}}\frac{1-\eta}{x^2+(\eta-1)^2} ( 1- {\rm e}^{-(x^2+(\eta-1)^2)/{a^2}} ), \end{gather}

(2.2)\begin{gather}u_{\eta}=\frac{\varGamma_0}{2{\rm \pi}}\frac{x}{x^2+(\eta-1)^2}( 1- {\rm e}^{-(x^2+(\eta-1)^2)/{a^2}}), \end{gather}

(2.3)\begin{gather}u_{\phi}=0, \end{gather}

where $\varGamma _0$ is the initial circulation and $a$ is the initial radius of the vortex core. Here, $a/R_0=0.2$ is applied throughout this work. The initial bubble is spherical, with a size of $R_b/R_0=0.284$ and a distance $0.2R_0$ from the ring centre. We note that the induced velocity has an axial component of $1/{2{\rm \pi} } \approx 0.159$ at the centre of the vortex ring, i.e. $u_x \approx 0.159$ at $(x, \eta )=(0, 0)$. To initialise the vortex ring from being perfectly circular, we perturb the local radius $\eta$ by adding a small parameter $\zeta \ll 1$ multiplied by the sum of a set of $N$ Fourier modes, each with unit amplitude and random phase, as adopted in a previous study (Archer, Thomas & Coleman 2008). After initialisation, the flow is relaxed quickly to the real solution by solving the Navier–Stokes equations. The velocity in the cylindrical coordinate system, written in (2.1) to (2.3), can be easily transformed to the Cartesian coordinate system used in our three-dimensional simulations.

Figure 1. A vortex ring and a gas bubble in the cylindrical coordinate system, where $x, \eta$ and $\phi$ are the axial, radial and azimuthal coordinates respectively, and the origin $O$ is at the centre of the ring. The gas bubble is released initially at the $x$ axis with a short distance $0.2R_0$ from $O$, as illustrated in this figure.

2.2. Numerical method

The flow involving two immiscible phases can be solved using a single-fluid continuum approach. We solve the time-dependent Navier–Stokes equations in incompressible form under a finite-volume framework, considering the translation of a liquid vortex ring and the deformation and translation of a gas bubble induced by the vortex ring. The differential forms of the governing equations are

(2.4)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}

(2.5)\begin{gather}\frac{\partial\rho\boldsymbol{u}}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{uu})={-}\boldsymbol{\nabla} p +\rho\boldsymbol{g}+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau} +\boldsymbol{F}_{\sigma}, \end{gather}

where $\boldsymbol {u}$ is the velocity vector, $p$ is pressure, $\boldsymbol {g}$ is the gravitational acceleration pointing in the negative direction of the $x$ axis, $\boldsymbol {\tau }$ represents the deviatoric stress tensor, $\boldsymbol {F}_{\sigma }$ is included to represent the surface tension of the liquid–gas interface and $\rho$ and $\mu$ are the mixture density and dynamic viscosity of the fluid, respectively, defined as

(2.6)\begin{gather} \rho = \varphi \rho_l + (1-\varphi)\rho_g, \end{gather}

(2.7)\begin{gather}\mu = \varphi \mu_l + (1-\varphi)\mu_g. \end{gather}

The subscripts $l$ and $g$ represent respectively the liquid and gas phases. The density ratio $\rho _l /\rho _g$ is set to be 1000, approximating that between water and air. Here $\varphi$ is the volume fraction of liquid that has a value of one in the liquid and zero in the gas. In the volume-of-fluid method, $\varphi$ represents volume fractions of one of the two fluids and varies smoothly between 0 and 1. The interface can be defined as the surface where the value of $\varphi$ is equal to 0.5. The volume fraction $\varphi$ is evolved using the following advection equation:

(2.8)\begin{equation} \frac{\partial\varphi}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\varphi \boldsymbol{u})=0.\end{equation}

An algorithm called isoAdvector (Roenby, Bredmose & Jasak 2016) based on the volume-of-fluid method is used to reconstruct geometrically and advect the phase interfaces. This algorithm is satisfactory in terms of volume conservation, boundedness, surface sharpness and efficiency. The surface tension term $\boldsymbol {F}_{\sigma }$ in (2.5), at any point on the interface, assuming the interfacial tension $\sigma$ is constant, can be obtained as

(2.9)\begin{equation} \boldsymbol{F}_{\sigma}={\sigma\kappa}\boldsymbol{\nabla}\varphi, \end{equation}

where $\kappa$ is the curvature of the interface calculated by the continuum surface force method (Brackbill, Kothe & Zemach 1992).

To provide satisfactory accuracy in space for the evolution of the vortex ring, as well as that of the gas bubble, several levels of refinements are applied to an initially uniform Cartesian grid. The mesh cells in a box, which encloses the vortex ring and its later developed flow structure, have a length of $0.01R_0$ in all three dimensions. In a smaller box, which encloses the gas bubble, the cell length is $0.005R_0$. Therefore, 80 cells are distributed along the bubble diameter. The whole computational domain has approximately 35 million cells. To ensure spacial accuracy as the flow evolves, we move the computational domain axially with the speed at which the vortex ring translates. The integral form of the governing (conservation) equations defined in a moving volume $V$ bounded by a closed surface $S$ can be found in our previous study (Deng & Caulfield 2018). As the volume $V$ is no longer fixed in space, its motion is captured by the motion of its bounding surface $S$ at the boundary velocity $\boldsymbol {u}_s$.

The space discretisations for (2.5) are second-order upwind for the convection terms and central differences for the Laplacian terms. The temporal term is discretised by the implicit Crank–Nicolson method. The pressure–velocity coupling in incompressible flow simulation is obtained using the pressure-implicit with splitting of operator scheme (Ferziger, Perić & Street 2002). At the outer boundary of the computational domain (a cube with a side length of $100R$) we set the normal components of the velocity gradients to be zero, and the pressure (excluding the hydrostatic pressure) to be zero. To ensure time-discretisation independence, we choose the size of time step, meeting the requirement $Co=\delta t |\boldsymbol {U}|/\delta x <1$ for all cells, in which $|\boldsymbol {U}|$ is the magnitude of the velocity through a cell and $\delta x$ is the cell width in the direction of the velocity.

2.3. Non-dimensional parameters

Since the vortex ring at $Re_\varGamma =7500$ has been extensively studied by both experimental measurements (Weigand & Gharib 1994) and numerical simulations (Bergdorf, Koumoutsakos & Leonard 2007), we fix this Reynolds number for most simulations. Here, $Re_{\varGamma }$ is defined as

(2.10)\begin{equation} Re_\varGamma=\rho_l\frac{\varGamma_0}{\mu_l}, \end{equation}

where $\mu _l$ is the dynamic viscosity of the liquid. The Weber number, which is a measure of the relative importance of the fluid's inertia compared to its surface tension, is defined as

(2.11)\begin{equation} We=\frac{0.87\rho_l(\varGamma_0/2{\rm \pi} a)^2}{(\sigma /D_b)}, \end{equation}

following a previous study (Oweis et al. 2005), which may be thought of as the ratio of the pressure difference (${\rm \Delta} P=0.87\rho (\varGamma _0/2{\rm \pi} a)^2$) between the vortex core and the far field to the Laplace pressure ($\sigma /D_b$) for a spherical bubble. In this study, $We$ varies in the range of 130.5–870. We also note that for a freely rising gas bubble in quiescent liquid, it can be characterised by two parameters, the Eötvös number and the Galileo number, defined respectively as

(2.12a,b)\begin{equation} Eo=\frac{\rho_{l}gD_b^{2}}{\sigma} \quad {\rm{and}}\quad Ga=\frac{\rho_{l}^{2}gD_b^{3}}{\mu_{l}^{2}},\end{equation}

where $g$ is the magnitude of gravitational acceleration and $D_b$ is the diameter for an equivalent spherical bubble (here $D_b=2R_b$). According to the currently chosen parameters, $Ga$ is fixed at $43.2$, and $Eo$ varies from $0.00048$ to $0.016$, corresponding to a perfectly spherical bubble when it is driven only by buoyancy force (Tripathi, Sahu & Govindarajan 2015). A direct comparison between the $We$ and $Eo$ ranges indicates that the flow induced by the vortex ring plays a profoundly more significant role than the buoyancy force.

3. Transition of the vortex ring at $Re_\varGamma =7500$

For a translating vortex ring, its main instability is referred to as the Tsai–Widnall instability, or azimuthal instability, which leads to an extensive bending of the initially circular vortex lines. Following a previous study (Archer et al. 2008), in three-dimensional isosurface plots, we distinguish between the region of intense vorticity at the core centre, which is called the ‘inner core’, and the surrounding outer-core region of lower vorticity, which we call ‘halo’ vorticity, as shown in figure 2. It is seen that after the formation of a number of standing waves around the core azimuth, due to a narrow band of independently growing linear instability modes, the flow enters into a nonlinear transition phase. At the early stage, the constructive interference between different modes causes a noticeable ‘lop-sidedness’ to the wave growth and associated core displacement. As shown in figure 2(b), the inner-core displacement becomes appreciable and the halo vorticity rolls up into an interwoven mesh of secondary structure, which consists of a series of loops encompassing the inner core. We note that neighbouring loops are of alternating signed vorticity, and two loops wrap around each azimuthal wave (see figure 2a), consistent with previous experimental observations (Dazin, Dupont & Stanislas 2006b) and simulations (Bergdorf et al. 2007; Archer et al. 2008).

Figure 2. Three-dimensional isosurface visualisations of the secondary structures at $t^*=130$ (a,b) and $t^*=160$ (c,d), where $t^*=t\varGamma _0/R_0^{2}$. In (a,c), the dark-grey surfaces correspond to the inner-core regions $|\omega |R_0^2/\varGamma _0=2.8$; the red and blue surfaces correspond to $\omega _x R_0^2/\varGamma _0=0.85$ and $\omega _x R_0^2/\varGamma _0=-0.85$ respectively, visualising the secondary structure. In (b,c), the light-grey surfaces correspond to $QR_0^4/\varGamma _0^2=0.2$.

At later times, the stationary coherent vortical structure found during the laminar and transitional phases is superseded by the swirling of vorticity filaments, as seen in figure 2(c,d). The core region is no longer stationary, but bends and twists with time, with vorticity filaments similar to the secondary structure generated continually. The vorticity filaments circulate around the core and gradually pass out of the vortex bubble and into the wake as a stream of vorticity filaments and hairpin vortices.

Transforming the simulated results back to the cylindrical coordinate system, we can investigate the global behaviour of flow evolution more conveniently. The circulation can be obtained by

(3.1)\begin{equation} \varGamma = \int_{{-}L}^{{+}L} \int_{0}^{\infty} \langle \omega_\phi \rangle \, {\rm d}\eta \,{\rm d}\kern0.06em x, \end{equation}

where $\langle \omega _\phi \rangle$ represents the mean azimuthal vorticity. The centre of the vortex ring in the ($x,\eta$) space is evaluated by

(3.2a,b)\begin{equation} x_o = \int_{{-}L}^{{+}L} \int_{0}^{\infty} x \langle \omega_\phi \rangle \, {\rm d}\eta \, {\rm d}\kern0.06em x,\quad \eta_o = \int_{{-}L}^{{+}L} \int_{0}^{\infty} \eta \langle \omega_\phi \rangle \, {\rm d}\eta \, {\rm d}\kern0.06em x. \end{equation}

To quantitatively study the characteristics of azimuthal instabilities, we perform the Fourier analysis of the axial vorticity $\omega _x(x,\eta,\phi )$. For each $\eta \unicode{x2013}\phi$ plane, the azimuthal component decomposition of $\omega _x$ is written as

(3.3)\begin{equation} \omega_k (x,\eta,\phi,t)=C_k(x,\eta,t) \cos (k\phi)+S_k(x,\eta,t) \sin (k\phi),\end{equation}

where $C_k$ and $S_k$ represent the Fourier coefficients defined by

(3.4)\begin{equation} \{ C_k(x,\eta,t),S_k(x,\eta,t) \} =\frac{1}{\rm \pi} \int_{0}^{2{\rm \pi}}\omega_x (x,\eta,\phi,t) \{ \cos(k\phi), \sin (k\phi) \} \, {\rm d}\phi.\end{equation}

Based on this decomposition, a measure of the amplitude of the azimuthal component of $\omega _x$ for each wave number on the axial plane is given by

(3.5)\begin{equation} A_k(x,t)=\left( \int_{0}^{2{\rm \pi}} \int_{0}^{\infty} {\omega_k^2 (x,\eta,\phi,t)} \eta \, {\rm d}\eta \, {\rm d} \phi \right)^{{1}/{2}}. \end{equation}

To investigate the development of instability, we measure the first 30 modes at three different $x$ slices, as shown in figure 3(a–d). It is clearly shown that the streamwise vorticity, $\omega _x$, appears with the onset of the Widnall instability, which is apparent in the excitation of the linear unstable mode $k=8$ (Bergdorf et al. 2007). At the early stage, $\omega _x$ is relatively weak from the pattern of $\omega _x$, as shown at $t^*=120$. With the rapid evolution of the secondary instability, energy is transferred to the second and third harmonics $k=16$ and $k=24$, as shown in figure 3(c) for $t^*=140$. With the vortices breaking into small-scale ones, the dominant mode and its harmonics decay rapidly and all the modes exhibit the same order of magnitude in figure 3(d), indicating the flow transition to a turbulent state. These observations agree well with the previous findings in both experiments (Dazin, Dupont & Stanislas 2006a; Dazin et al. 2006b) and numerical simulations (Bergdorf et al. 2007). To further confirm the turbulent state as the nonlinear instability fully develops, we present the energy $E_k$ for different azimuthal modes in figure 4(a), which can be calculated by

(3.6)\begin{equation} E_k(x,t)=\int_{0}^{2{\rm \pi}} \int_{0}^{\infty} \frac{1}{2}|{\boldsymbol{u}_k^{\prime}(x,\eta,\phi,t)}|^2 \eta \, {\rm d}\eta \, {\rm d} \phi,\end{equation}

where $\boldsymbol {u}_k^{\prime }$ is the azimuthal component decomposition of the fluctuation velocity, i.e. the flow velocity subtracting the average translating speed of the vortex ring. It is identified that the law of energy decay follows a characteristic $k^{-5/3}$ law, indicating that the vortex ring has broken into a turbulent state. Alternatively, we create a structured $256\times 256 \times 256$ isotropic mesh, enclosing the vortex ring, and map the flow velocities onto this mesh. Integrating over all wavenumbers $\boldsymbol {k}_2$ of magnitude $|\boldsymbol {k}_2|=k$ for the trace of a velocity-spectrum tensor, we obtain the energy-spectrum function $E_{k2}$, as shown in figure 4(b) (Pope 2000). Here, we use $k_2$ to distinguish from the azimuthal mode number $k$. The inertial subrange characterised by the well-known Kolmogorov $k^{-5/3}$ law is seen more clearly for the higher Reynolds number $Re_{\tau }=15\ 000$, which spans down to a length scale around $0.11R_0$. We should point out that at $Re_{\tau }=15\ 000$ the vortex ring is still in the laminar-transitional regime, without being fully turbulent yet, as noted previously (Foronda-Trillo et al. 2021).

Figure 3. (a–d) Development of the first 30 harmonics. The three lines with different colours in the spectra represent different measuring planes $x= x_0-0.5a$, $x_0$ and $x_0+0.5a$ (referring to (3.2a,b)). (e–h) The corresponding $\omega _x$ structures, where the initial vortex ring is marked by a circle.

Figure 4. (a) Development of azimuthal modal energies at three sequential times: $t^*=140$, $t^*=160$ and $t^*=200$, for $Re_{\tau }=7500$. The energy decay for $t^*=200$ follows a characteristic $k^{-5/3}$ law, which is denoted by a dash-dotted line. (b) The traditional Kolmogorov energy spectra for $Re_{\tau }=7500$ and $Re_{\tau }=15\ 000$.

4. Interaction of a vortex ring with a bubble

We now investigate how a bubble affects the evolution of a translating vortex ring and, vice versa, the bubble deformation caused by the strain of the vortical flow. The bubble is released from the position shown in figure 1 at time $t^*=10$, when the vortex ring has relaxed to a divergence-free flow. We present in figure 5 the normalised horizontal position of the vortex ring as it interacts with the bubble for different $We$ cases. The motion of the same ring in the absence of the bubble has also been included as the baseline case. The translation speed at any given instant can be easily obtained as the slope of the position–time data shown. For the vortex ring alone, we know that the formula of translation or convective speed was found as

(4.1)\begin{equation} U_c = \frac{\varGamma_0}{4{\rm \pi} R_0} \left[ \ln \frac{8R_0}{a} + C \right],\end{equation}

where $C$ is a function of the shape of the vorticity distribution and $a$ is a measure of the core radius. For a uniform distribution $C=-1/4$, and for a Gaussian core vorticity distribution $C\simeq -0.558$ (Saffman 1970). For a translating vortex ring, the core spreads over time to be of order $a \sim a_0 + (\nu t)^{1/2} \sim a_0 + (\nu \varGamma )^{1/2}R_0$, where $a_0$ is the initial value of the core radius, indicating that the decay of the translation speed is dominated by viscous diffusion. Fukumoto & Moffatt (2008) estimated the translation speed of a vortex ring in a viscous fluid medium as

(4.2)\begin{equation} U_c = \frac{\varGamma_0}{4{\rm \pi} R_0} \left[ \ln \left( \frac{4R_0}{\sqrt{\nu t}} \right) -0.558 - 3.6716 \frac{\nu t}{R_0^2} \right].\end{equation}

In figure 5, we see that the distance travelled by the vortex ring interacting with a bubble is less than that by the non-interacting vortex ring for all Weber numbers. It is interesting to find that the difference between the ring–bubble interacting case and the baseline ring increases as $We$ is decreased in the low-$We$ cases, i.e. $We=130, 174$ and $217$, implying an increased contribution from the surface tension forces. For example, at $We=130$, there is a significant decline in the translation speed of the ring at early times when it starts to interact with the bubble (see the inset of figure 5), which is more pronounced at later times (see the main plot of figure 5). In contrast, we cannot find a notable monotonic behaviour in the high-$We$ range. This difference between low- and high-$We$ ranges implies a different ring–bubble interacting mechanism, which we discuss in more detail by examining their flow structures.

Figure 5. The horizontal position ($x/R_0$) of the vortex ring as a function of time. The results with different Weber numbers are divided into two groups: data for the low $We$ values are shown with solid lines of different colours and those for the high $We$ values are shown with symbols. The black solid line represents the motion of the same ring without bubbles. In the inset, the data in a smaller time range are shown.

We also plot the histories of the vortex translational velocity in figure 6, in which the decaying speed provided by the empirical formula (4.2) is also presented for reference. In figure 6(a), in the low-$We$ range, we observe that the translational speed of the vortex ring for $We=130$ decays more rapidly, and eventually to a much lower level than the other cases. A comparison between figures 6(a) and 6(b) shows that the low-$We$ range exhibits stronger fluctuations of the translational speed at the early stage, due to the oscillating motion (accelerating and decelerating process) of the bubble along the vortex ring axis before being trapped into the vortex core.

Figure 6. Histories of vortex translational speed for (a) the low-$We$ range and (b) the high-$We$ range. In (a), the black line is the Saffman formula (4.2), which has been shifted to match the numerical results. In the inset of (a), the stepwise decay of the translational velocity is observed for the case without bubbles.

To quantify the strength of secondary flow structure due to the azimuthal symmetry breaking as well as the flow transition to turbulence, we measure the total enstrophy according to the streamwise vorticity, expressed as

(4.3)\begin{equation} \varOmega_x=\frac{1}{2}\int \omega_x^2 \, {\rm d} V. \end{equation}

We show the evolution of $\varOmega _x$ in figure 7. It is seen that the rapid growth of enstrophy is earlier for the ring–bubble interacting cases than that in the absence of a bubble, indicating the earlier appearance of secondary flow structure. In the low-$We$ range, the lowest Weber number, i.e. $We=130$, presents the largest growth rate of enstrophy and reaches the peak earlier than the rest, as shown in figure 7(a). In contrast, we cannot find a strong correlation between Weber number and the evolution of enstrophy in the high-$We$ range, shown in figure 7(b). However, it is interesting to notice that there is a bump around $t^*\simeq 20$ as marked by the dashed circle in figure 7(b), which is due to the separation of the boundary layer around the bubble when it is accelerated or decelerated by the vortex, as addressed in a previous study (Foronda-Trillo et al. 2021). The previous study has also addressed that this detached vorticity is responsible for triggering the vortex instability. Clearly, the larger the Weber number, the more intense is this secondary vortex, as the bubble deforms more intensely.

Figure 7. Evolution of enstrophy for the secondary vorticity filaments wrapping around the primary vortex ring in the flow field for (a) the low-$We$ range and (b) the high-$We$ range. The case without a bubble is also included in both panels.

As implied in figure 7, the interaction between the vortex ring and bubble is complicated and correlated to the surface tension forces of the bubble or equivalently the Weber number ($We$). To distinguish between the low- and high-$We$ ranges, we select two representative cases from each, and present their evolution of flow structure and bubble deformation, as shown in figures 8 and 10 for $We=130$ and $We=348$ respectively. In figure 8, at the early stage, we see that the bubble translates with the vortex ring, while staying at the centre of the ring, when the whole flow structure is still symmetric. As seen at $t^*=15$, a small secondary vortex ring is induced by the bubble, which propagates backwards through the primary vortex ring (see $t^*=25$), accompanied by the generation of another secondary ring. The first secondary ring translates further downstream, and eventually disappears in the wake (see $t^*=45$). The repeatedly generated secondary rings perturb the primary one, leading to its azimuthal symmetry breaking, as shown at $t^*=45$ and $t^*=55$. The primary vortex ring then exhibits the same unstable wavenumber, resembling that shown in figure 2, i.e. $k=8$, for the vortex ring in the absence of bubbles. Then, the bubble moves towards the ring core and is captured by it, due to the low pressure within the core (see $t^*=65$). During this process of bubble entertainment, we also notice an overshoot of the bubble, as seen at $t^*=75$. After that, as the bubble re-approaches the ring core, the bubble breaks into two, as seen at $t^*=85$. Here, we note that $We=130$ corresponds to $We_\tau =3$, defined as $We_{\tau }=\rho (\varGamma _0/2 {\rm \pi}R_0)^2 / (\sigma / R_b )$. In a previous study (Foronda-Trillo et al. 2021), with a three-dimensional simulation, similar binary breakup of the bubble was also reported at $We_\tau =0.75$ with $R_b/R_0=1/3$, while at a different Reynolds number $Re_\tau =15\ 000$ and the bubble moves in the opposite direction of a translating vortex ring. It is easy to identify that the number of small bubbles is $n=12$ after breakup.

Figure 8. Time sequence of (a) top view and (b) side view images, presenting the interaction of a bubble with a vortex ring at $We=130$. The light-grey surfaces correspond to $QR_0^4/\varGamma _0^2=0.2$. The blue surfaces correspond to the volume fraction of liquid $\varphi =0.5$, representing the translation and breakup of the bubble. Note that the vortex ring translates towards the outside direction of the paper in (a), while from the left to right in (b).

Shown in figure 9 are the magnitudes of the first 30 modes at different time instants for $We=130$, corresponding to figure 8. They resemble the development of the first 30 harmonics without bubbles (see figure 3), while with an earlier transition to turbulent states. The overshoot of the bubble induced by its inertial migration is clearly seen in figure 9(g), resulting in symmetric distributions of the streamwise vorticity regarding the shooting direction.

Figure 9. (a–d) Development of the first 30 harmonics for $We=130$. The spectra are measured on the plane $x=x_0$ (referring to (3.2a,b)). (e–h) The corresponding $\omega _x$ structures, where the initial vortex ring is marked by a circle.

In figure 10, we present the evolution of interacting vortex ring and bubble for $We=348$. The corresponding development of the first 30 harmonics for this case is shown in figure 11. Instead of deviating from the ring centre, as shown in figure 8 for the low-$We$ range, the high-$We$ cases exhibit a significantly different behaviour. Attracted by the low pressure in the ring core, the bubble is stretched radially to a very thin ring, as shown at $t^*=15$. This bubble ring continues to expand in the radial direction, until the surface tension forces cannot keep the integrity of the bubble ring. It breaks into multiple fragments, or small bubbles, as shown at $t^*=20$, where a further expansion of the fragmented bubble ring is seen, as the small bubble reaches the outermost region of the vortex ring. This is followed by a contraction of the bubbles towards the ring centre. After that, the bubbles sit in the centre of the vortex ring core (see from $t^*=25$ to $t^*=65$). As seen at these time instants, the vortex ring encloses a set of smaller bubbles of unequal sizes, which does not further break up, resembling that from a previous experimental visualisation (see figure 1f in Jha & Govardhan (2015)) though with a different set-up. A direct comparison between figures 8 and 10 shows that the secondary vortex filaments become apparent at $t^*=65$ for $We=130$, while at $t^*=85$ for $We=348$, corresponding respectively to their rising stages of the enstrophy histories in figure 7. In figure 10, though much weaker, a secondary vortex ring induced by the spreading bubble can be observed, which stays behind as the small bubbles are entrained into the primary vortex ring, as shown at $t^*=20$.

Figure 10. Time sequence of (a) top view and (b) side view images, presenting the interaction of a bubble with a vortex ring at $We=348$. The light-grey surfaces correspond to $QR_0^4/\varGamma _0^2=0.2$. The blue surfaces correspond to the volume fraction of liquid $\varphi =0.5$, representing the translation and breakup of the bubble. Note that the vortex ring translates towards the outside direction of the paper in (a), while from the left to right in (b).

Figure 11. (a–d) Development of the first 30 harmonics for $We=348$. The spectra are measured on the plane $x=x_0$ (referring to (3.2a,b)). (e–h) The corresponding $\omega _x$ structures, where the initial vortex ring is marked by a circle.

To understand more clearly the stretching, spreading and subsequent breakup process of the bubble, we present in figure 12 the time sequence of the early bubble spreading stage with smaller time intervals. It is interesting to find that the bubble spreads to a flat disk before the thin bubble ring is formed, which gets flatter and flatter until it is torn from the centre, as seen at $t^*=13.5$. The ring bubble breaks into multiple fragments as its radius increases beyond a critical value (between $t^*=16.5$ and $t^*=18$). Here, the critical radius for the bubble ring is identified as $r_c=0.76 R_0$. We also note that the small bubbles move along a spiralling path around the vortex core during the capture process, before they reside in the core centre. Another phenomenon worthy of attention is that at around $t^*=22.5$ when pairs of fine streamwise vortical filaments induced by the small bubbles are observed. It could potentially trigger the azimuthal instability mode of the wavenumber $k=12$, depending on the number of small bubbles. However, this mode of perturbation decays quickly after the bubble capture process completes, with only weak vortex legs remaining in the wake, as seen from $t^*=25$ to $t^*=65$ in figure 10. Instead, the mode $k=8$ instability retakes the dominance, and the primary vortex ring undergoes the secondary transition resembling that without bubbles.

Figure 12. (a–j) Time sequence of vortex–bubble interaction in the early stage, with smaller time intervals, at $We=348$. Note that the light-grey surfaces correspond to $QR_0^4/\varGamma _0^2=0.1$, and the blue surfaces correspond to the volume fraction of liquid $\varphi =0.5$. The vortex ring translates towards the outside direction of the paper.

After examining the vortex–bubble interacting process, it is easy to identify the number of smaller bubbles for each Weber number case. As exhibited in figure 12, there is a fixed, distinguishable number of smaller bubbles before the vortex ring transits further to turbulence. Shown in figure 13 is the number of smaller bubbles as a function of $We$. For the three low-$We$ cases, the number is fixed at 2, while it increases with $We$ in the high-$We$ range, reaching 16 at $We=391$. This maximum number of 16 agrees coincidentally with a previous experimental observation (Jha & Govardhan 2015). The number of smaller bubbles keeps at 16 as $We$ further increases, up to the highest Weber number $We=870$ investigated in the current study. From the stretching and spreading process of the bubble observed in figure 12, we conjecture that the classical Rayleigh–Plateau instability is a possible mechanism determining the number of smaller bubbles after breakup, which is a simplified model examining a stretched cylindrical bubble. In this model, the wavelength corresponding to the maximum growth rate for the instability is about 4.5 times the diameter of the cylindrical bubble. In the current study, first, we can identify the maximum radius of the bubble ring, and thereby its azimuthal length $L_b$, before its breakup, as demonstrated in figure 12. Then we can calculate the tube diameter $D_t$ based on the gas volume $V_g$, which does not change due to the conservation of the gas volume, and $L_b$. Lastly, we can estimate the number of smaller bubbles by $n=L_b/(4.5*D_t$). To compare with our numerical simulations, we include the estimated bubble numbers in figure 13. Clearly, the visualised numbers of broken-off bubbles obtained in our simulations agree well with that estimated by the Rayleigh–Plateau instability theory. It is therefore easy to understand that the maximum number of smaller bubbles is limited to 16, due to the fact that the maximum radius of the bubble ring is always smaller than that of the primary vortex ring.

Figure 13. Variation of the number of bubbles with $We$ after complete breakup of the bubble. The number of Rayleigh–Plateau instability waves that fit on the cylindrical bubble is also shown.

In addition, we present the development of modal energies for $We=130$ and $We=348$ respectively in figures 14 and 15. It is not surprising that the energy decay follows a characteristic $k^{-5/3}$ law as the flow develops to fully turbulent states. However, comparing with the case without bubbles (see figure 4) and that with only two broken small bubbles (see figure 14), the $We=348$ case with 12 small bubbles distributed azimuthally along the perturbed vortex ring presents distinctly different spectra, as shown in figure 15, where a clear $k^{-1}$ law is identified. We understand that the $-1$ power region of the turbulence spectrum is a constant energy region responsible for energy production. Here, the $-1$ power region implies the continuous energy injection caused by the two-phase interactions at the early stage of turbulence. However, since it overlaps the inertial range as the turbulence fully develops (see the $-5/3$ power at $t^*=160$), this energy-containing region cannot be distinguished at that time. It will be of interest to simulate at higher Reynolds numbers; then it will be possible to observe these two regions simultaneously. We leave it as an open question.

Figure 14. Development of modal energies for $We=130$ at sequential times: $t^*=68, t^*=85$, $t^*=140$ and $t^*=160$. The energy decay for $t^*=160$ follows a characteristic $k^{-5/3}$ law, which is denoted by a dash-dotted line.

Figure 15. Development of modal energies for $We=348$ at sequential times: $t^*=100$, $t^*=105, t^*=120, t^*=140$ and $t^*=160$. The energy decay for the late time follows a characteristic $k^{-5/3}$ law, while for the time instants $t^*=100$ and $t^*=105$ it follows a characteristic $k^{-1}$, denoted by two dash-dotted lines.

5. Discussions on the effects of Reynolds number

In the previous experimental study, the only variable was the piston velocity, resulting in different ring radii, core radii and circulation strengths of the vortex rings (Jha & Govardhan 2015). Therefore, the Weber number always varies with the Reynolds number. In this study, we investigate these two key parameters, Weber number and Reynolds number, separately. Specifically, we investigate the effects of increasing the Reynolds number while keeping the Weber number fixed at $We=217$, to determine how the capture and deformation processes of a bubble are affected.

To expedite the capture process, we release the bubble slightly off-axis. Unlike the oscillatory and translational motion of the bubble along the ring axis observed in figure 8, the bubble immediately enters the ring core, as seen in figure 16. Subsequently, the bubble expands in the azimuthal direction along the core, which is consistent with previous visualisation experiments that released the bubble from outside the ring (Jha & Govardhan 2015). Below a certain critical Reynolds number, the elongated bubble undergoes binary breakup, resulting in two fragments. It is interesting to note that this binary breakup has also been observed for lower Reynolds numbers at the same Weber number of $We=217$, as we have discussed above. At higher Reynolds numbers, such as $Re_{\tau }=30\,000$, the bubble initially breaks into three fragments (see $t^*=30$ in figure 16b). The middle fragment then undergoes a binary breakup, coalescence and re-breakup process, resulting in the formation of four daughter bubbles. In both cases, the bubble is captured in the right core of the vortex ring, causing the right core to fragment while the left core remains intact, as the bubble does not reach it.

Figure 16. Time sequence of vortex–bubble interaction in the early stage, with smaller time intervals for two higher Reynolds numbers (a) $Re_{\tau }=25\,000$ and (b) $Re_{\tau }=30\,000$, at a fixed Weber number $We=217$. Note that the light-grey surfaces correspond to $QR_0^4/\varGamma _0^2=0.1$, and the blue surfaces correspond to the volume fraction of liquid $\varphi =0.5$. The vortex ring translates towards the outside direction of the paper. See also supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.517 for the case in (b).

As a turbulent flow is multiscale in nature and bubbles present in this flow are subject to breaking due to collision with multiple eddies, the dynamics of the breakup process becomes more complicated. As illustrated in figure 16, a direct comparison between two Reynolds numbers reveals that the finer eddies in a vortical field with a higher Reynolds number tend to tear the bubble into more fragments.

6. Concluding remarks

Differentiating from previous experiments (Jha & Govardhan 2015), in which a bubble was released from outside the ring, while following closely their way of analysis, we study numerically the interactions between a vortex ring and a gas bubble released at the axis. As addressed in a previous ‘Focus on Fluids’ paper (Martínez-Bazán 2015), which made a brief review and discussion of that experimental work, such an elementary system can represent an idealised scenario providing a framework to shed light on understanding more complicated bubble–turbulence coupling. In the current study, we fix the Reynolds number at $Re_{\tau }=7500$ and consider various Weber numbers in the range of $We=130\unicode{x2013}870$.

First, we investigate the role that bubbles play in the vortex ring dynamics. Clearly shown is that the translating speed of the vortex ring is substantially lower than that of the base of the vortex ring without a bubble. By looking at the enstrophy or volume integral of the square of the streamwise vorticity, we find that the rapid growth of enstrophy is earlier for the ring–bubble interaction cases than that without bubbles, indicating the earlier breaking of azimuthal symmetry, or flow transition. This reduction of ring translating speed, or the earlier flow transition, can be explained with two different mechanisms, depending on the Weber number. For the low-$We$ cases, consistent with that revealed in a previous study (Jha & Govardhan 2015), the reduction of convection speed of the ring is due to the capture of bubbles into the ring core, leading to significant changes in the vorticity distribution within the core. However, in the high-$We$ cases, the enstrophy rises before the bubbles are captured into the ring core. Instead, the repeatedly generated secondary vortex rings perturb the primary one, which bring about an earlier flow transition.

On the other hand, we examine the evolution of gas bubbles affected by the strain of the vortical flow. For relatively low Weber numbers, we observe binary breakup of the bubble after it is captured by the primary vortex ring. In contrast, for high Weber numbers, it is interesting to find that the bubble experiences sequentially stretching, spreading and breakup stages. Shortly after release, the bubble spreads to a flat disk, which gets flatter and flatter, until it is torn from the centre, forming a bubble ring. This bubble ring continues to expand in the radial direction, until reaching a critical radius, when the surface tension cannot keep the integrity of the bubble ring. After that, it breaks into multiple fragments, or smaller bubbles. The number of these smaller bubbles is distinguishable. For the three low-$We$ cases, the number is fixed at 2, while it increases from 4 to 16 with $We$ in the high-$We$ range. Furthermore, we find that the numbers of smaller bubbles predicated by the classical Rayleigh–Plateau instability of a stretched cylindrical bubble agree well with our numerical simulations in the high-$We$ range. Consistent with previous experiments (Jha & Govardhan 2015), this number reaches 16 as $We$ is sufficiently large.

Additionally, we also investigate the effects of Reynolds number. A direct comparison between two Reynolds numbers reveals that the finer eddies in a vortical field with a higher Reynolds number tend to tear the bubble into more fragments.

A future extension of the present work can be devoted to the cases with even lower Weber numbers, when the breakup of the bubble will not be observed. The critical Weber number for bubble breakup can be identified. Moreover, it would be interesting to examine the effects of gravity, as well as that of the initial bubble position, as a direct comparison with previous experiments. In that scenario, the initial accelerating process of the bubble driven by the buoyancy force will play an important role in the bubble–ring interactions.