Hostname: page-component-745bb68f8f-v2bm5 Total loading time: 0 Render date: 2025-02-07T03:14:49.998Z Has data issue: false hasContentIssue false

What happens to the vortex structures when the rising bubble transits from zigzag to spiral?

Published online by Cambridge University Press:  04 September 2017

Jie Zhang
Affiliation:
State Key Laboratory for Strength and Vibration of Mechanical Structures, School of Aerospace, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, PR China
Ming-Jiu Ni*
Affiliation:
State Key Laboratory for Strength and Vibration of Mechanical Structures, School of Aerospace, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, PR China School of Engineering Science, University of Chinese Academy of Sciences, Beijing 101408, PR China
*
Email address for correspondence: mjni@ucas.ac.cn

Abstract

It has been demonstrated by many experiments carried out over the last 60 years that in certain liquids a single millimetre-sized bubble will rise within an unstable path, which is sometimes observed to transit from zigzag to spiral. After performing several groups of direct numerical simulations, the present work gives a theoretical explanation to reveal the physical mechanism causing the transition, and the results are presented in two parts. In the first part, in which a freely rising bubble is simulated, equal-strength vortex pairs are observed to shed twice during a period of the pure zigzag path, and this type of motion is triggered by the amounts of streamwise vorticities accumulated on the bubble interface, when a critical value is reached. However, when the balance between the counter-rotating vortices is broken, an angular velocity is induced between the asymmetric vortex pairs, driving the bubble to rise in an opposite spiral path. Therefore, although there is no preference of the spiral direction as observed in experiments, it is actually determined by the sign of the stronger vortex thread. In the second part, external vertical magnetic fields are imposed onto the spirally rising bubble in order to further confirm the relations between the vortex structures and the unstable path patterns. As shown in our previous studies (Zhang & Ni, Phys. Fluids, vol. 26 (10), 2014, 102102), the strength of the double-threaded vortex pairs, as well as the imbalance between them, will be weakened under magnetic fields. Therefore, as the vortex pairs become more symmetric, the rotating radius of the spirally rising bubble is observed to decrease. We try to answer the question, put forward by Shew et al. (2005, Preprint, ENS, Lyon), ‘what caused the bubble to transit from zigzag to spiral naturally?’

Type
Papers
Copyright
© 2017 Cambridge University Press 

1 Introduction

When a millimetre-sized bubble is released freely in a liquid, the buoyant forces will drive the bubble to rise and the kinetic energy of the surrounding fluid is increased correspondingly. Depending on the physical properties of the liquids, the motion of the bubble may present unstable paths, such as zigzag or spiral, which have been summarized by Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a ) recently. Since Mougin & Magnaudet (Reference Mougin and Magnaudet2001) numerically verified that the path instability of the rising bubble was caused by the wake instability behind the bubble, the relations between the bubble paths and wake structures became clearer. However, at present several problems remain unclear, such as what drives the rising bubble to transit from zigzag to spiral naturally, as put forward by Shew, Poncet & Pinton (Reference Shew, Poncet and Pinton2005), Magnaudet & Eames (Reference Magnaudet and Eames2000) and Ern et al. (Reference Ern, Risso, Fabre and Magnaudet2012).

By using hyper-purified water, Duineveld (Reference Duineveld1995) confirmed for the first time that the bubble rising in water will present an oscillatory path when the bubble radius $R>0.91~\text{mm}$ , corresponding to a Reynolds number $Re\approx 660$ . Thereafter, more and more experimental studies were carried out to determine the relationships between the wake structures and the path instabilities (Lunde & Perkins Reference Lunde and Perkins1997; Brücker Reference Brücker1999; Ellingsen & Risso Reference Ellingsen and Risso2001). The spiral path was observed to associate with a continuous wake composed of double-threaded vortex pairs, while the zigzag motion was accompanied by the shedding of hairpin vorticities. However, most of the experiments suffered from contamination because of the added materials. Thus in order to keep the liquids highly purified while visualizing the wakes behind the bubble, a significant effort has been carried out over the last 15 years (De Vries, Biesheuvel & Van Wijngaarden Reference De Vries, Biesheuvel and Van Wijngaarden2002; Sanada, Shirota & Watanabe Reference Sanada, Shirota and Watanabe2007; Veldhuis, Biesheuvel & Van Wijngaarden Reference Veldhuis, Biesheuvel and Van Wijngaarden2008; Zenit & Magnaudet Reference Zenit and Magnaudet2009). Most of these observed a shedding process of double-threaded vortices at the rear of the bubble during the zigzag stage, while the vortex structures were more continuous and stable in the spiral stage. Moreover, by using different silicone oils, which do not need an ultrapure environment, Zenit & Magnaudet (Reference Zenit and Magnaudet2008) found that it was actually the aspect ratio $\unicode[STIX]{x1D712}$ rather than $Re$ that triggered the path instability.

Because of the difficult visualizing technique in highly purified water, numerical simulations were gradually employed to study the bubble dynamics by keeping the bubble with a frozen shape. Mougin & Magnaudet (Reference Mougin and Magnaudet2001) proved that the path instability was coupled with the unstable wake. They imposed manual perturbations to cause the bubble to transit from zigzag to spiral. Then Mougin & Magnaudet (Reference Mougin and Magnaudet2006) and Shew, Poncet & Pinton (Reference Shew, Poncet and Pinton2006) found the evolutions of the forces and torques on the bubble were different when it travelled within zigzag and spiral paths. Magnaudet & Mougin (Reference Magnaudet and Mougin2007) reported the critical value of $Re$ and $\unicode[STIX]{x1D712}$ for the appearance of wake instability by studying the uniform flow past an ellipsoidal bubble, as an extension of the numerical investigations by Leal (Reference Leal1989) and Blanco & Magnaudet (Reference Blanco and Magnaudet1995). On the basis of this, Cano-Lozano, Bohorquez & Martínez-Bazán (Reference Cano-Lozano, Bohorquez and Martínez-Bazán2013) performed a more precise study with realistic fore-and-aft asymmetric bubble shapes. Furthermore, by performing the linear stability analysis (Yang & Prosperetti Reference Yang and Prosperetti2007; Tchoufag, Magnaudet & Fabre Reference Tchoufag, Magnaudet and Fabre2013, Reference Tchoufag, Magnaudet and Fabre2014; Cano-Lozano et al. Reference Cano-Lozano, Tchoufag, Magnaudet and Martínez-Bazán2016b ), several unstable modes were identified for a rising bubble whether it was freely moving or set in uniform flows.

As far as the direct numerical simulations are concerned, in which the bubble is allowed to deform freely, the available studies are really rare, to the authors’ knowledge. Hua, Lin & Stene (Reference Hua, Lin and Stene2008) and Gaudlitz & Adams (Reference Gaudlitz and Adams2009) were able to show periodic vortex shedding during the zigzag motion, and the later work confirmed different vortex structures during the zigzag and spiral motions. Within the volume-of-fluid (VOF) framework, Tripathi, Sahu & Govindarajan (Reference Tripathi, Sahu and Govindarajan2015) and Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a ) confirmed that the double-threaded vortices were shed periodically in the zigzag stage while the vortex structures were more stable in the spiral stage. Zhang & Ni (Reference Zhang and Ni2014b ) and Zhang, Ni & Moreau (Reference Zhang, Ni and Moreau2016) observed distinct hairpin-like vortex structures when studying an Ar bubble rising in liquid metal GaInSn.

As summarized above, although the relation between the wake instability and the path instability is already recognized, it is still not clear why the bubble is sometimes observed to transit from a zigzag to a spiral. As put forward by Shew et al. (Reference Shew, Poncet and Pinton2005), ‘to unravel the causes of the transition from zigzag to spiral. Is one wake vortex stronger or, as we suggest, are the wake vortices simply unstable to rotation?’ Further, although it is believed that the accumulation of vorticities on the bubble surface will lead to vortex shedding, there is still no quantitative evidence about this. In the present study, we will answer the above questions by presenting the evolutions of the vortex structures during the zigzag and spiral motions. According to the experimental results by Zenit & Magnaudet (Reference Zenit and Magnaudet2008), path instability will even appear within low- $Re$ flows as long as the aspect ratio of the bubble exceeds a critical value. In the low- $Re$ regime, the wake structures are thought to be more elemental and regular than those under higher- $Re$ flows, where the vortex structures would be much more complicated. In addition, high- $Re$ flows contain more flow instabilities, so that it is difficult to identify the physical mechanism for the path transition. It is also pointed out by Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a ) that, for simulating high- $Re$ flows, capturing the boundary layer requires much thinner grids close to the bubble, and this is really time-consuming or otherwise the results are not reliable. As a consequence, the path instability, together with the wake instability, will be studied within lower- $Re$ flows in the present study. The mechanisms could also be used to understand the path transition from zigzag to spiral in high- $Re$ flows.

In particular, as presented in our previous study (Zhang & Ni Reference Zhang and Ni2014b ), an external vertical magnetic field (MF) will modify the double-threaded vortex structures markedly. Therefore, an MF can be used to further identify the relations between the unstable path patterns and the vortex structures. More information about the influence of a vertical MF on the vortex structures can be found in that paper. Here, the main conclusions are summarized briefly as follows. (1) Under vertical MFs, the twisted counter-rotating vortex threads will be more parallel with one another and the strengths of the vortex pairs are significantly weakened in the whole flow field. This influence is attributed to two aspects: the more spherical shape of the bubble and more uniform flow fields under vertical MFs. (2) Stronger vortex threads will be weakened more significantly because it is a typical influence of the MFs that the vorticities along with the MFs will be more homogeneous, as summarized by Moreau (Reference Moreau2013) and Sommeria & Moreau (Reference Sommeria and Moreau1982). After understanding this, we will further validate that it is rather the vortex structures, especially the relative strength between the double vortex threads, that determine the different path patterns of the rising bubble.

The paper is organized in the following manner. In § 2 we state the problem and in § 3 we describe the numerical methods, respectively. The numerical results are presented in § 4 and show how the vortex structures evolve when transiting from zigzag to spiral; moreover, the physics of the transition is discussed. Furthermore, external vertical MFs are imposed on the spirally rising bubble in order to validate the physical mechanisms of the transition. Finally, the conclusions are drawn in § 5.

2 Physical model

A spherical bubble is released at the bottom of a rectangular container where it will rise freely under gravity, as shown in figure 1. As presented, gravity is directed along the $z$ -direction, while the $x{-}y$ plane is the horizontal cross-section. The computational domain, with a size of $36R\times 36R\times 108R$ in the present study, is thought to be large enough for the bubble motion so that the boundary effect can be ignored.

Figure 1. Sketch of the physical model. If an external MF is imposed, it is along the gravity direction, given as $B_{z}$ .

Bubbly flows are thought to be incompressible so that they are governed by the Navier–Stokes equations

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}\right)=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}+\boldsymbol{F}_{s}+\boldsymbol{S}, & \displaystyle\end{eqnarray}$$
(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D64E}=\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}})$ is the viscous stress tensor, $\boldsymbol{F}_{s}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FF}_{s}\boldsymbol{n}$ stands for the surface tension force, and $\boldsymbol{S}=\unicode[STIX]{x1D70C}\boldsymbol{g}$ is the gravitational force. In addition, $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D707}$ represent the density and the dynamic viscosity of the continuous or disperse phase, respectively; $\unicode[STIX]{x1D70E}$ is the surface tension coefficient; $\unicode[STIX]{x1D705}$ is the interface curvature; $\boldsymbol{n}$ is the normal direction of the interface; and the force centralized on the interface is represented by the Dirac distribution function $\unicode[STIX]{x1D6FF}_{s}$ .

The following variables are introduced to rescale the equations:

(2.3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle L^{\ast }=\frac{L}{R},\quad \boldsymbol{u}^{\ast }=\frac{\boldsymbol{u}}{\sqrt{gR}},\quad t^{\ast }=\sqrt{\frac{g}{R}}t,\\ \displaystyle \unicode[STIX]{x1D70C}^{\ast }=\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}_{l}},\quad p^{\ast }=\frac{p}{\unicode[STIX]{x1D70C}_{l}gR},\quad \unicode[STIX]{x1D707}^{\ast }=\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D707}_{l}},\\ \displaystyle \unicode[STIX]{x1D705}^{\ast }=R\unicode[STIX]{x1D705},\quad \boldsymbol{g}^{\ast }=\frac{\boldsymbol{g}}{\Vert \boldsymbol{g}\Vert },\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where the subscript $l$ denotes the physical properties of the ambient liquid surrounding the bubble. Consequently, the Navier–Stokes equations (2.1) can be rewritten as

(2.4) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}^{\ast }}{\unicode[STIX]{x2202}t}+\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{\ast }=-\unicode[STIX]{x1D735}p^{\ast }+\frac{1}{Ga}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D707}^{\ast }(\unicode[STIX]{x1D735}\boldsymbol{u}^{\ast }+\unicode[STIX]{x1D735}\boldsymbol{u}^{\ast \text{T}}))+\frac{1}{Eo}\unicode[STIX]{x1D705}^{\ast }\unicode[STIX]{x1D6FF}_{s}^{\ast }\boldsymbol{n}+(\unicode[STIX]{x1D70C}^{\ast }-1)\boldsymbol{g}^{\ast }. & & \displaystyle\end{eqnarray}$$

Therefore, the dynamic behaviour of the rising bubble is actually determined by four dimensionless parameters, two of which are $\unicode[STIX]{x1D70C}^{\ast }$ and $\unicode[STIX]{x1D707}^{\ast }$ and the other two of which are defined as

(2.5a,b ) $$\begin{eqnarray}\displaystyle Ga=\frac{\unicode[STIX]{x1D70C}_{l}\sqrt{gR}R}{\unicode[STIX]{x1D707}_{l}},\quad Eo=\frac{\unicode[STIX]{x1D70C}_{l}gR^{2}}{\unicode[STIX]{x1D70E}}, & & \displaystyle\end{eqnarray}$$

where $Ga$ stands for the ratio of the gravitational force to the viscous force, and $Eo$ is the ratio of the gravitational force to the surface tension force. In fact, they are variant forms of the Reynolds number and Weber number ( $We$ ) obtained by replacing $u$ with the scaled velocity $\sqrt{gR}$ , given as

(2.6a,b ) $$\begin{eqnarray}\displaystyle Re=\frac{\unicode[STIX]{x1D70C}_{l}uR}{\unicode[STIX]{x1D707}_{l}},\quad We=\frac{\unicode[STIX]{x1D70C}_{l}u^{2}R}{\unicode[STIX]{x1D70E}}. & & \displaystyle\end{eqnarray}$$

In the present study, the density ratio and viscosity ratio between the gas phase and the liquid phase, namely $\unicode[STIX]{x1D70C}^{\ast }$ and $\unicode[STIX]{x1D707}^{\ast }$ , are fixed at $10^{-3}$ and $10^{-2}$ , respectively. Therefore, the bubble motions are entirely characterized by $Ga$ and $Eo$ according to (2.4). Another important dimensionless parameter employed to describe the flows is the Morton number ( $Mo$ ), defined as $Mo=Eo^{3}/Ga^{4}$ , which is only dependent on the physical properties of the liquid.

Moreover, when there are external MFs imposed on the electrically conducting flows, an electromagnetic force will be induced because of the motion of the conductive fluid, and it takes the form of $\boldsymbol{F}_{L}=N(\boldsymbol{J}\times \boldsymbol{B})$ added to the right-hand side of (2.1) and (2.4). In this formulation, $N$ is the dimensionless parameter indicating the ratio of the electromagnetic force to the inertial force, and it is defined as $N=\unicode[STIX]{x1D70E}_{e}B^{2}R/(\unicode[STIX]{x1D70C}\sqrt{gR})$ in the present study. Here, $\unicode[STIX]{x1D70E}_{e}$ is the electric conductivity of the fluid and $\boldsymbol{B}$ is the applied magnetic vector. Moreover, the induced current density, which is denoted $\boldsymbol{J}$ , can be calculated as $\boldsymbol{J}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}+\boldsymbol{u}\times \boldsymbol{B}$ , where $\unicode[STIX]{x1D711}$ is the induced electric potential. Owing to the conservative character of the current density, a divergence-free condition of $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{J}=0$ needs to be satisfied. Therefore, the solution of the electromagnetic field must solve an additional electric potential Poisson equation, giving as $\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D711})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\times \boldsymbol{B})$ .

3 Numerical methodology

The results to be discussed below are obtained by solving the three-dimensional Navier–Stokes equations in the entire fluid domain. In order to accomplish this, the open-source software of Gerris flow solver (Popinet Reference Popinet2009) is employed herein. The VOF method coupled with the adaptive refinement technique are implemented in the solver, leading it to be a very efficient tool in simulating multiphase flows. Recently, another two papers have been published that use Gerris to study the dynamic behaviour of single bubble motion, respectively by Tripathi et al. (Reference Tripathi, Sahu and Govindarajan2015) and Cano-Lozano et al. (Reference Cano-Lozano, Tchoufag, Magnaudet and Martínez-Bazán2016b ), and both of them show very reliable results. In the present study, the grid has a density of 32 cells per bubble radius close to the interphase, corresponding to the spatial resolution adopted by Tripathi et al. (Reference Tripathi, Sahu and Govindarajan2015), and the grid density in the wake region is also refined automatically into the smallest size of $R/16$ . Because we only consider the low- $Re$ flow regime, the grids are thought to be fine enough to capture the boundary layer and to resolve the vortex structures at the bubble rear, as validated by Cano-Lozano et al. (Reference Cano-Lozano, Tchoufag, Magnaudet and Martínez-Bazán2016b ). The equations and the discretization techniques implemented in Gerris are described in detail by Popinet (Reference Popinet2009) and will not be repeated here. Also, numerous validations are given by Tripathi et al. (Reference Tripathi, Sahu and Govindarajan2015) and Cano-Lozano et al. (Reference Cano-Lozano, Tchoufag, Magnaudet and Martínez-Bazán2016b ) too.

When it comes to solving the electromagnetic field, the details of the numerical methods, as well as the validation tests, can be found in our previous papers (Zhang & Ni Reference Zhang and Ni2014a ,Reference Zhang and Ni b ; Zhang et al. Reference Zhang, Ni and Moreau2016) and they will not be repeated here either. The methodology proves to be stable and accurate even when the ratio of the electric conductivity between the ambient fluid and bubble is as large as $10^{5}$ .

4 Numerical results

In the present study, different flow conditions are investigated with a limited regime of $Ga$ and $Eo$ , which are varied between 35.36 and 70.72 and between 0.7 and 3, respectively. However, in order to make the paper more compact, we will focus on the flow condition of $Ga=50$ while $Eo$ is varied between $1.0$ , $2.0$ and $3.0$ in the following study; more results are available in appendix A. We will demonstrate later that they come to similar conclusions when the bubble travels within a zigzag or spiral path. Under $Ga=50$ , the terminal $Re$ approaches $150$ , and the rising trajectories of the bubble centroid, with a three-dimensional view, are given in figure 2. As shown in the figures, different rising paths are observed as $Eo$ increases, from an absolutely zigzag path ( $Eo=1.0$ ) to a zigzag–spiral transition path ( $Eo=2.0$ ) and finally a perfect spiral path ( $Eo=3.0$ ).

It should be noted that our results are fully in accordance with the numerical results reported by Cano-Lozano et al. (Reference Cano-Lozano, Tchoufag, Magnaudet and Martínez-Bazán2016b ), who draw a phase diagram summarizing the different path patterns observed in their simulations, as presented in figure 22 of that paper. According to that diagram, the bubble motions in the present study (the parameters become $\{Ga,Eo\}=\{141.42,4\},\{141.42,8\},\{141.42,12\}$ if $D=2R$ is used as characteristic length) should fall into the planar zigzag, flattened spiral and spiral regions, respectively.

Figure 2. The different paths of the bubble ( $Ga=50,~\unicode[STIX]{x1D70C}^{\ast }=1/1000,~\unicode[STIX]{x1D707}^{\ast }=1/100$ ):  $Eo=1$ for zigzag (a), $Eo=2.0$ for zigzag–spiral transition (b) and $Eo=3.0$ for spiral (c).

4.1 Zigzag motion

Regarding the zigzag motion of $Eo=1.0$ , the transverse positions of the bubble centroid are drawn in figure 3(a) after a coordinate transformation, in order to make the new $x^{\prime }$ -axis orient along with the zigzag direction. As it is believed that such an oscillatory motion is connected with the dynamic behaviours of the streamwise vortex structures, the evolution of streamwise vorticity isocontours $\unicode[STIX]{x1D714}_{z}^{\ast }=\pm 0.5$ are presented in the smaller panels at four typical time points, when the bubble passes the furthermost and average positions respectively. In those smaller panels, the double-threaded wakes (one thread blocks the view of the other) are observed to be shed periodically from the bubble interface so that their signs are changed twice during one zigzag period, the frequency being nearly $f_{z}=1/11$ . This periodic shedding process is also presented by Cano-Lozano et al. (Reference Cano-Lozano, Tchoufag, Magnaudet and Martínez-Bazán2016b ); however, what we are concerned with here is the reason causing such vortex shedding.

Consequently, we plot the time histories of the norm of streamwise vorticities accumulated on the bubble surface, that is, $\Vert \unicode[STIX]{x1D70F}_{z}\Vert =\int _{S}\Vert \unicode[STIX]{x1D714}_{z}^{\ast }\Vert \,\text{d}S$ , where $S$ is the bubble interface. The result is shown in figure 3(b), where a periodic oscillation is also observed, and the frequency is twice that of the zigzag motion, being $f_{\unicode[STIX]{x1D714}}=1/5.5$ . Precisely, it is observed that, whenever the bubble passes the furthermost positions, such as A and C in figure 3(a), $\Vert \unicode[STIX]{x1D70F}_{z}\Vert$ decreases to the minimum value; in contrast, when the bubble passes the average positions, such as B and D in figure 3(a), $\Vert \unicode[STIX]{x1D70F}_{z}\Vert$ reaches the maximum value. This indicates that there is a critical value of $\Vert \unicode[STIX]{x1D70F}_{z}\Vert$ to trigger the vortex shedding. However, if we calculate the total vorticities ( $\Vert \unicode[STIX]{x1D714}^{\ast }\Vert =\sqrt{\unicode[STIX]{x1D714}_{x}^{2}+\unicode[STIX]{x1D714}_{y}^{2}+\unicode[STIX]{x1D714}_{z}^{2}}$  ) accumulated on the bubble surface, there are no such periodic oscillations observed, as shown in figure 3(c). Therefore, it is rather the accumulation of streamwise vorticities on the bubble surface that triggers the vortex shedding. Regarding this critical value, it is probably related with the curvature of the bubble interface, as presented by Magnaudet & Mougin (Reference Magnaudet and Mougin2007). In fact, for different flow conditions, we find the critical value of $\Vert \unicode[STIX]{x1D70F}_{z}\Vert$ that triggers the zigzag motion varies with $Ga$ and $Eo$ . In addition, with a given flow condition of constant $Ga$ and $Eo$ , the zigzag motion can be triggered earlier as long as the critical value of $\Vert \unicode[STIX]{x1D70F}_{z}\Vert$ is reached. For instance, when two bubbles rise side by side, the zigzag motion emerges earlier than in the single bubble case because $\Vert \unicode[STIX]{x1D70F}_{z}\Vert$ accumulates more quickly in the presence of another bubble; a more detailed analysis is presented by Zhang, Chen & Ni (Reference Zhang, Chen and Ni2017).

Figure 3. (a) Time histories of the horizontal positions of the zigzag bubble ( $Eo=1.0$ ), and the vortex structures within a zigzag period, respectively at $t^{\ast }=39$ , $42$ , $45$ and $48$ . In the vortex contour maps, the light colour (green) corresponds to $\unicode[STIX]{x1D714}_{z}^{\ast }=0.5$ while the dark colour (blue) is $\unicode[STIX]{x1D714}_{z}^{\ast }=-0.5$ . (b) The time histories of the norm of the streamwise vorticity accumulated on the bubble surface. (c) The time histories of the norm of the total vorticity accumulated on the bubble surface.

The mechanism to generate such a zigzag motion is shown in figure 4(a), which presents the top view of the counter-rotating vortices. In the figure, the solid line is the bubble interface while the dotted line corresponds to vortex pair II shown in figure 3(a). According to the theory of vortex pairs (Batchelor Reference Batchelor2000), when the vortex strengths between the counter-rotating vortices are equal, a horizontal velocity directing to the right will be induced, and under the conservation law of momentum, the bubble will be pushed by the liquid to the left. This theoretical interpretation is consistent with the numerical results. Nevertheless, it should be noted that, to induce such a straight-line motion, the vortex strength between the counter-rotating vortices should be equal. To quantitatively prove this, the isosurface of streamwise vorticity, $\unicode[STIX]{x1D714}_{z}^{\ast }=\pm 0.5$ , on a plane $4R$ downstream from the bubble centre is shown in figure 4(b). It is observed that the double-threaded vortex pairs are well matched in strength. Moreover, through calculation, the time history of the arithmetic integration of streamwise vorticities, being $\unicode[STIX]{x1D70F}_{plane}=\int \unicode[STIX]{x1D714}_{z}^{\ast }\,\text{d}S$ on that plane, is also shown in figure 4(b), whereas the value of zero indicates that, during the zigzag motion, the counter-rotating vortex pairs are perfectly symmetric. It should be noted that, in the following calculations, $\int \Vert \unicode[STIX]{x1D714}_{z}^{\ast }\Vert \,\text{d}S$ is very different from $\int \unicode[STIX]{x1D714}_{z}^{\ast }\,\text{d}S$ ; the latter distinguishes the vortices between positive values and negative values in order to investigate the balance between the positive and negative vortices, so we call it ‘arithmetic integration’.

Figure 4. (a) The sketch of the vortex pairs behind the zigzag bubble with $Eo=1.0$ . (b) The isosurface of $\unicode[STIX]{x1D714}_{z}^{\ast }=\pm 0.5$ on the plane of $4R$ downstream from the bubble centre. It is observed that they are well matched in strength. (c) The time series of the arithmetic integration of $\unicode[STIX]{x1D714}_{z}^{\ast }$ on the cut plane as $\unicode[STIX]{x1D70F}_{plane}=\int \unicode[STIX]{x1D714}_{z}^{\ast }\,\text{d}S$ , proving that the counter-rotating vortex pairs are perfectly symmetric.

Figure 5. (a) Time histories of the horizontal positions of the zigzag–spiral bubble ( $Eo=2.0$ ), and the evolution of the vortex structures during the transition stage from zigzag to spiral, respectively at $t^{\ast }=40$ , $42$ , $44$ and $46$ . In the contour map, the light colour corresponds to $\unicode[STIX]{x1D714}_{z}^{\ast }=0.5$ and the dark colour is $\unicode[STIX]{x1D714}_{z}^{\ast }=-0.5$ . (b) Time histories of $\Vert \unicode[STIX]{x1D714}_{z}^{\ast }\Vert$ accumulated on the bubble surface.

4.2 Transition from zigzagging motion to spiral motion

When the rising bubble transits from zigzag to spiral, in the case of $Eo=2$ , the time series of the transverse positions of the bubble centroid are plotted in figure 5(a). It is observed that, when $t^{\ast }<40$ , the bubble travels within a zigzag path, which develops to a spiral one after $t^{\ast }>40$ . Correspondingly, the evolution of $\Vert \unicode[STIX]{x1D70F}_{z}\Vert$ is also drawn in figure 5(b), in which the variation trend is different from that in the zigzag motion. The oscillation amplitude of $\Vert \unicode[STIX]{x1D70F}_{z}\Vert$ is much smaller, indicating that vortex shedding is not apparent in the spiral stage, and this is also observed in experiments (De Vries et al. Reference De Vries, Biesheuvel and Van Wijngaarden2002).

Moreover, our real concern is why this transition happens. The time evolution of the wake structures during the transition is also attached in figure 5(a), recording from $t^{\ast }=40$ to $46$ . In the isosurface contour maps, the light colour corresponds to $\unicode[STIX]{x1D714}_{z}^{\ast }=0.5$ while the dark colour is $\unicode[STIX]{x1D714}_{z}^{\ast }=-0.5$ . Through a close inspection, it is observed that the new vortex pairs, which are just shedding from the bubble surface, gradually lose symmetry and tend to twist with one another, finally leading the bubble to rotate within a spiral path. It should be noted that the positive vortex thread lies in the inner side while the negative one is at the outside. As a consequence, the asymmetry of the vortex pairs seems to be responsible for the path transition from zigzag to spiral. Moreover, as widely known, the amounts of vorticity accumulated on the bubble surface are directly dependent on the interface curvature (Magnaudet & Mougin Reference Magnaudet and Mougin2007). Consequently, the side views of the bubble, when the negative vortices and the positive vortices, respectively, are generated, are presented in figure 6 when it travels in a spiral and zigzag motion. It is obvious that, in the spiral stage, the curvature at the $y-$  side, where the positive vortices are shed from the bubble, is larger than that at the $y+$  side, where negative vortices are shed. Therefore, the positive vortex thread is gradually dominant over the negative one. In contrast, during the zigzag motion, as is shown in figure 6(b), the two sides are exactly symmetrical; thus the strengths between the double vortex threads are also equivalent.

Figure 6. The $y\pm$  side views of the bubble shape when it travels in (a) spiral motion (asymmetric shape) and (b) zigzag motion (symmetric shape), respectively.

Table 1. The top row shows the streamwise vorticity distributions on the bubble surface at different time periods. The bottom row shows the contour maps of the maximum interface curvature at corresponding times. The values of the contour levels are attached in the colour scale. It is observed that, during the zigzag period from $t^{\ast }=35$ to $37$ , both the interface curvatures and positive and negative vorticity distributions are symmetrical; however, in the transition stage from $t^{\ast }=41$ to $47$ , both symmetries are broken, so that a larger surface curvature results in stronger positive vortices; and after entering the spiral motion at $t^{\ast }>53$ , positive vortices completely dominate on the bubble surface because of the much larger local interface curvature.

To be more detailed on this issue, in the case $Eo=2$ , the distribution of the positive and negative streamwise vorticities on the bubble surface, as well as the interface curvature at corresponding times, are presented in table 1. It should be noted that we calculate the maximum of the absolute value of the principal curvatures, referred to as $\unicode[STIX]{x1D705}_{max}$ in table 1. In the diagrams, we would find that the positive (negative) vortex thread is shed at the right (left) side of the bubble. Several time periods are selected: In the first stage from $t^{\ast }=35$ to $37$ , the bubble still travels in the zigzag motion; it is obvious that positive and negative $\unicode[STIX]{x1D714}_{z}^{\ast }$ are symmetrically distributed at the two sides of the bubble, and correspondingly $\unicode[STIX]{x1D705}_{max}$ are also equal in both vortex regions. However, in the second period, when the bubble motion transits from zigzag to spiral, namely from $t^{\ast }=41$ to $47$ , the positive vortices start to prevail over the negative ones; in the meantime, the interface curvature also becomes larger at the positive vortex region. In the last stage after $t^{\ast }>53$ , when the spiral motion is finally generated, it is found that the positive vortices are much stronger than the negative ones, and the local interface curvature is also more pronounced. Therefore, it is quantitatively validated that, when the bubble shape is symmetrical at the two sides where positive and negative streamwise vortices are shed, the vorticity distributions on the bubble surface are also symmetrical; however, when the interface curvatures become asymmetric, the balance between the positive and negative vortices is also broken.

Figure 7. (a) Sketch of the motion of the counter-rotating vortex pairs when they are unequal in strength. (b) Top view of the vortex structures behind the spirally rising bubble. (c) Contour map of the streamwise vorticity in the cross-plane $4R$ downstream from the bubble centre. (d) Time evolution of the arithmetic summation of $\unicode[STIX]{x1D714}_{z}^{\ast }$ on the cross-plane. It is observed that the value jumps from zero to a positive value around $t^{\ast }=40$ , indicating that the inner positive vortex thread is dominant gradually in the transition stage.

Unlike the symmetric vortex pairs during the zigzag motion, the unequal vortex pairs will induce an angular velocity to drive them into rotation around an invariant vorticity centre, as shown by Leweke, Le Dizès & Williamson (Reference Leweke, Le Dizès and Williamson2016). This is depicted in figure 7(a). When it comes to the present problem, the top view of the vortex pairs is shown in figure 7(b), where the black line is the bubble interface. In figure 7(b), the inner positive vortex thread is stronger than the negative one, so that an anticlockwise angular velocity is thus induced to make the vortex pairs rotate with one another. To conserve momentum, the bubble will rotate in a clockwise direction so that a spiral motion is produced gradually. This conjecture is consistent with the result shown in figure 2(b), where a clockwise rotation of the rising bubble is observed. The isosurfaces of $\unicode[STIX]{x1D714}_{z}^{\ast }$ on the plane $4R$ downstream from the bubble centre are also plotted in figure 7(c). It is clear that the structures of the counter-rotating vortices are asymmetric, being very different from those during the zigzag motion. To further validate that the inner positive vortex thread is stronger than the negative one, the time history of the arithmetic integration of the streamwise vorticities, being $\unicode[STIX]{x1D70F}_{plane}=\int \unicode[STIX]{x1D714}_{z}^{\ast }\,\text{d}S$ on the cut plane, is presented in figure 7(d). From figure 7(d), it is observed that $\unicode[STIX]{x1D70F}_{plane}$ jumps from zero to a positive value around $t^{\ast }=40$ , indicating that the inner positive vortex thread dominates gradually in the transition stage, and the growth of the imbalance finally leads to a path transition from zigzag to spiral.

Another test case to further illustrate the physics of the path transition is presented in appendix B, the result of which is almost identical with the conclusion we draw here. The difference is that we impose perturbations on the bubble surface tension in order to produce asymmetries on the bubble shape as well as vortex pairs. As a consequence, the zigzagging bubble transits into the spiral stage.

Moreover, we notice that the asymmetric vortex structures in the spiral regime were also observed by Popinet (Reference Popinet2017), who is now developing a more powerful code called Basilisk by using numerical algorithms similar to those implemented in Gerris. Through an in-depth discussion with S. Popinet (2017, personal communication), the correctness of our results is further identified.

Therefore, we answer the question raised by Shew et al. (Reference Shew, Poncet and Pinton2005) and Ern et al. (Reference Ern, Risso, Fabre and Magnaudet2012). It is rather the asymmetric vortex pairs, which are produced by the asymmetric deformation of the bubble shape in the present case, that induce an angular velocity that drives the bubble in a spiral motion. This interpretation is consistent with the experimental observation by Brücker (Reference Brücker1999), who also suggested that the counter-rotating vortex pairs were no longer symmetric when the bubble transited from zigzag to spiral. In figure 12 of that paper, it is the same appearance as we observe above: the bubble rotates in the opposite direction against the inner thread of the vortex pairs. However, he did not give any further analysis or interpretations; besides, this observation did not receive additional attention in the following studies either.

Consequently, although there is no preference for the bubble to rotate in the clockwise or anticlockwise direction when rising in spiral motions, as summarized in the experimental studies, the relative strength between the double-threaded vortex pairs actually determines its rotation direction.

As mentioned earlier, more results of different flow conditions are presented in appendix A, as shown in tables 2 to 4, which correspond to the flow condition of $Ga=35.36$ , $Ga=42.43$ and $Ga=71.72$ , respectively. From the diagrams in the tables, it is observed that, at a given $Ga$ , the rising path of the single bubble transits from zigzag to spiral gradually as $Eo$ is increased, and this is not only qualitatively but also quantitatively consistent with the phase map given by Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a ) (their figure 22). By investigating the evolution of the vortex structures and the vortex strengths during different rising stages, we come to similar conclusions as drawn under $Ga=50$ , and detailed discussions can be found in appendix A.

4.3 Magnetic influence on the spiral motion

After detecting the physical mechanism for a path transition from zigzag to spiral, we intend to validate that, by diminishing the imbalance between the asymmetric vortex pairs, the radius of the spiral motion would decrease, or even return to a rectilinear motion. As introduced in § 1, it is found in our previous study (Zhang & Ni Reference Zhang and Ni2014b ) that the counter-rotating vortex pairs are weakened under the influence of vertical MFs, which will be used herein to provide further confirmation of our theoretical analysis. We select the flow condition of $Ga=50$ and $Eo=3$ for investigation, whereas a pure spiral motion is observed without MFs, as shown in figure 2(c).

After being released, the bubble rises freely without MFs till $t^{\ast }=42$ , so that the spiral motion has been generated already. At the moment $t^{\ast }=42$ , different vertical MFs, respectively $N_{z}=0.3$ , $N_{z}=0.5$ and $N_{z}=1.0$ , are imposed in order to study the influence of the magnetic strength. The subscript $z$ herein indicates the MFs is vertically directed.

To eliminate the influence of the varying terminal rising velocities on the bubble paths when imposing MFs, the time histories of the rising velocity under different $N_{z}$ are presented in figure 8. As found in our previous study (Zhang & Ni Reference Zhang and Ni2014b ), moderate vertical MFs will increase the terminal rising velocity of the bubble a little. This is also validated in figure 8, whereas a little larger terminal velocities are observed under MFs. Therefore, we do not think the variations of the rising velocities play a notable role if the bubble path is changed.

Figure 8. The time histories of the rising velocities under different MFs. It is observed that, except for the relatively lower velocity without MF, there is nearly no difference among other cases when $N_{z}>0$ .

After that, we calculate the arithmetic integration of $\unicode[STIX]{x1D714}_{z}^{\ast }$ integrated over the bubble interface before and after imposing vertical MFs, with $\unicode[STIX]{x1D70F}_{bubble}=\int _{S}\unicode[STIX]{x1D714}_{z}^{\ast }\,\text{d}S$ . The results are plotted in the top row of figure 9, where the red (green) lines depict the value before (after) $t^{\ast }=42$ . From figure 9, it is observed that, without MF, $\unicode[STIX]{x1D70F}_{bubble}$ is almost stabilized at $\unicode[STIX]{x1D70F}_{bubble}=3$ , indicating that there are more positive streamwise vorticities accumulated on the bubble interface. However, in the rest of the panels, while the MF is increased to $N_{z}=0.3$ , $N_{z}=0.5$ and $N_{z}=1.0$ , the value of $\unicode[STIX]{x1D70F}_{bubble}$ is observed to drop to $\unicode[STIX]{x1D70F}_{bubble}=2$ , $\unicode[STIX]{x1D70F}_{bubble}=1.5$ and $\unicode[STIX]{x1D70F}_{bubble}=0$ respectively. Therefore, under the influence of vertical MFs, the positive and negative streamwise vorticity are more symmetric on the bubble surface.

Figure 9. The evolution of the arithmetic summation of $\unicode[STIX]{x1D714}_{z}^{\ast }$ over the the bubble surface (ad) and the projection of the rising paths of the bubble (eh) under vertical MFs. The MFs are introduced at $t^{\ast }=42$ . It is observed that, with the increase of $N_{z}$ , the positive and negative streamwise vorticities are more symmetrical on the bubble surface, resulting in a decreased radius of spiral motion.

Subsequently, the projections of the rising paths are also presented in the bottom row of figure 9, where red (green) lines stand for paths before (after) vertical MFs being imposed. From the plots, it is clearly shown that, after imposing external MFs, the radius of the spiral motion is decreased gradually, complying with $r_{N_{z}=0}>r_{N_{z}=0.3}>r_{N_{z}=0.5}$ . This variation is caused by the more symmetrical counter-rotating vortex pairs. Moreover, in the case $N_{z}=1.0$ , the radius of the spiral motion is decreased to zero because the driving force for spiral motion vanishes as $\unicode[STIX]{x1D70F}_{bubble}=0$ . By comparing the time histories of $\unicode[STIX]{x1D70F}_{bubble}$ and the path under $N_{z}=1.0$ , it is observed that, even when $\unicode[STIX]{x1D70F}_{bubble}$ has already decreased to zero at $t^{\ast }=50$ , the bubble still travels within a weakened spiral motion, indicating that the variation of the spiral path lags behind the variation of the vortices; therefore, it is rather the vortex structures that determine the rising path.

On the basis of our theoretical analysis, after fully eliminating the imbalance between the asymmetric vortex pairs, there is another possibility for the path transition of the spirally rising bubble, that is recovering from spiral to zigzag again. In the subsequent paragraph we will interpret why this does not happen under the influence of vertical MFs.

Let us go back to § 4.1, where we demonstrate that, unless the accumulation of the streamwise vorticities on the bubble interface reaches a critical value, then the vortex shedding would happen to trigger the zigzag motion. As a consequence, the evolution of $\Vert \unicode[STIX]{x1D70F}_{z}\Vert =\int _{S}\Vert \unicode[STIX]{x1D714}_{z}^{\ast }\Vert \,\text{d}S$ , which depicts the norm of streamwise vorticities accumulated on the bubble surface when vertical MFs are imposed, are presented in figure 10.

Figure 10. Time evolution of the norm of $\Vert \unicode[STIX]{x1D714}_{z}^{\ast }\Vert$ accumulated on the the bubble surface when it rises under different vertical MFs. The MFs are introduced at $t^{\ast }=42$ . It is observed that $\Vert \unicode[STIX]{x1D70F}_{z}\Vert$ decreases with the intensification of the vertical MFs. In particular, in the case $N_{z}=1.0$ , where the vortex pairs are already symmetrical as shown in figure 9, it results in the rectilinear but not the zigzag path because $\Vert \unicode[STIX]{x1D714}_{z}^{\ast }\Vert$ is not large enough to trigger the vortex shedding.

Figure 11. The evolutions of the double-threaded vortex pairs with or without vertical MFs, respectively with $N=0$ , $N_{z}=0.5$ and $N_{z}=1.0$ . The time period is selected from $t^{\ast }=43$ , when the MFs are just imposed onto the spirally rising bubble. In the contour map, the light colour corresponds to $\unicode[STIX]{x1D714}_{z}^{\ast }=0.5$ while the dark colour is $\unicode[STIX]{x1D714}_{z}^{\ast }=-0.5$ . It is observed that in the case $N=0$ the vortex pairs twist with one another, and the structures are very stable. However, under $N_{z}=0.5$ , the twisted structures seem to be weakened, but never stand parallel with one another. When the magnetic strength increases to $N_{z}=1.0$ , the disintegration of the twisted vortex structures is quicker, and after $t^{\ast }>49$ , there is nearly no streamwise vorticities shed from the bubble surface.

In figure 10(b), it is observed that after entering the spiral stage, $\Vert \unicode[STIX]{x1D70F}_{z}\Vert$ almost stabilized around a value of $24$ ; however, after imposing vertical MFs, $\Vert \unicode[STIX]{x1D70F}_{z}\Vert$ decreases gradually as $N_{z}$ becomes larger. Let us focus on the case of $N_{z}=1.0$ , in which the vortex pairs are fully symmetric as shown in figure 9. It is found that in such a case $\Vert \unicode[STIX]{x1D70F}_{N_{z}=1.0}\Vert =2.5$ is too small to trigger the vortex shedding, so that the spiral motion becomes rectilinear directly.

After understanding how the radius of the spiral motion is decreased under the influence of vertical MFs, the evolutions of the double-threaded vortex pairs under $N=0$ , $N_{z}=0.5$ and $N_{z}=1.0$ are presented in figure 11, within a time period of $43<t^{\ast }<55$ . In the contour maps, the light colour corresponds to $\unicode[STIX]{x1D714}_{z}^{\ast }=0.5$ while the dark colour is $\unicode[STIX]{x1D714}_{z}^{\ast }=-0.5$ . From the figure, in the case $N=0$ , it is observed that the twisted vortex structures drive the bubble to rotate within a spiral motion. Meanwhile, the vortex pairs seem to be very stable so that no vortex shedding happens. However, in the case of $N_{z}=0.5$ , the twisted vortex structures are weakened gradually over time because the imbalance between the vortex pairs is attenuated, as explained in § 1. However, even so we do not observe the twisted vortex structures to break up completely, and the weakened twisted structures are stable again after $t^{\ast }>49$ . In the last case of $N_{z}=1.0$ , the evolution process of the vortex structures is similar, except that the disintegration of the twisted vortex structures is quicker. Ultimately, after $t^{\ast }>49$ , the vorticities accumulated on the bubble surface are too few to trigger any path instability, and the bubble finally rises in a rectilinear path.

5 Conclusion

The present work has detailed a mechanism explaining why and how the bubble travels from zigzag to spiral. Regarding rising within the zigzag motion, we find it is the accumulation of the streamwise vorticities on the bubble surface, rather than the total vorticities, that trigger the vortex shedding in a zigzag period. When the bubble transits from zigzag to spiral, previous studies always described the different behaviours of the bubble in the two travelling patterns, or different vortex structures behind the bubble. Our mechanism, which is based on the quantitative calculation of the counter-rotating vortex strengths, however, explains how different vortex structures are evolved, and how an angular velocity is induced between the asymmetrical vortex pairs to drive the bubble travelling from zigzag to spiral. Furthermore, by investigating the evolution of the interface curvature and vortex distribution on the bubble interface during the zigzag–spiral transition stage, we find the asymmetrically deforming bubble shape causes the unbalanced vortex distribution, where more vorticities are shed from the interface with a larger local curvature. To further validate this, in the spirally rising regime, we are able to produce asymmetric bubble shape and vortex shedding by imposing perturbations on the surface tension of the bubble, and a path transition from zigzag to spiral is correspondingly induced.

In the second part of the paper, by imposing external vertical MFs, which are used to diminish the asymmetry between the vortex pairs, the spiral motion is observed to be weakened, which finally transits into the rectilinear path under much stronger MFs. Nevertheless, we give detailed explanations about why a zigzag motion is not recovered; it is because not only symmetric vortex pairs are needed, but also strong enough vortex pairs are required to trigger vortex shedding. The influences of vertical MFs on vortex structures are also reported in detail in our previous study (Zhang & Ni Reference Zhang and Ni2014b ), and the present study is rather a particular case, which only focuses on the spiral stage of the bubble.

Furthermore, to make things simple, we only investigate the bubble motion in low- $Re$ flows, in which almost no other flow instability happens. Therefore, the imbalance between the vortex pairs can be ascribed to the asymmetric deformation of the bubble shape during the motion. However, when $Re$ becomes higher, other instabilities in the flow field would also lead the vortex pairs to be asymmetric, and thus spiral or even chaotic rising trajectories are observed.

Acknowledgements

This work has been supported by NSFC (51636009, 11502193), CAS (XDB22040201) and MOST (2013GB114000). The authors thank Professor J. Magnaudet and Professor S. Popinet for discussions about the results in this paper. We also thank the anonymous reviewers for their comments, which improved the quality of the paper dramatically.

Appendix A. More numerical simulations in broader parameter spaces

Besides $Ga=50$ , we have also considered other flow conditions in this study in order to provide more comprehensive results of the bubble motions. Their characteristic properties are presented in tables 24, which correspond to the flow conditions of $Ga=35.36$ , $Ga=42.43$ and $Ga=71.72$ , respectively (they are actually $Ga=100$ $120$ and $200$ by Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a ) if we use the diameter instead of the radius as the characteristic length). From the left row to the right, the top view of the rising path, the streamwise vorticity isocontours of $\unicode[STIX]{x1D714}_{z}^{\ast }=\pm 0.4$ and the evolution of $\unicode[STIX]{x1D70F}_{plane}$ on a plane $4R$ downstream from the bubble centre are displayed for each case. It should be noted that the bubbles start to rise at the initial horizontal position of $\{X,Y\}=\{1,0\}$ . As shown in the diagrams in the tables, with fixed $Ga$ , the rising path of the bubble is observed to transit from zigzag to spiral gradually by increasing $Eo$ , and this process of transition not only qualitatively, but also quantitatively, agrees with the phase map given by Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a ) (figure 22 in that paper) when we pay attention close to the transition borderline of the path instability. In particular, in the zigzag regime, i.e. $\{Ga,Eo\}=\{35.36,1.0\},\{42.43,1.5\}$ and $\{70.72,0.7\}$ , the vortex pairs are observed to shed from the bubble surface periodically and their relative strengths are also equal. On the other hand, in the spiral regime, i.e. $\{Ga,Eo\}=\{35.36,3.0\}$ and $\{70.72,3.0\}$ , the vortex pairs twist with one another so that the inner thread (negative one in $\{35.36,3.0\}$ and positive one in $\{70.72,3.0\}$ ) is stronger than the outside one through calculation. It is also observed that, if the inner vortex thread is positive (negative), the spiral path rotates in the clockwise (anticlockwise) direction. The characteristic behaviours in both regimes come to the same conclusion as we draw under $Ga=50$ in § 4.

Figure 12. Time evolution of the streamwise vorticity isocontours $\unicode[STIX]{x1D714}_{z}^{\ast }=\pm 0.4$ during the flattened spiral motion. The flow conditions are $Ga=35.36$ and $Eo=2.0$ . The time period is selected between $t^{\ast }=75$ and $83$ . It is observed that the vortex structures behave as the combined manner of the vortex shedding and the vortex twisting.

In particular, another typical path pattern belonging to the flattened spiralling regime has been identified, i.e.  $\{Ga,Eo\}=\{35.36,2.0\}$ and $\{42.43,2.0\}$ in the present study. This regime has already been reported by Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a ) but it is not observed in the flow condition of $Ga=50$ . According to Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a ), they find this is almost a long transient state that will eventually converge to either a zigzag or a spiral regime. However, this transient may be very long, so that in most of the experiments it is viewed as a stable state. According to our numerical results, it seems to be stable too, particularly in the case of $\{Ga,Eo\}=\{35.36,2.0\}$ , as shown in table 2. The vortices appear to be more or less the intermediate state between the parallel (zigzag) and twisted (spiral) vortex structures; and also, as displayed in the rightmost row, the inner vortex thread (negative one in both cases) is also proved to be a little stronger than the one at the outside. Consequently, the asymmetry between the vortex pairs diverts the bubble motion away from a pure zigzag path; however, it is not sufficient to generate a spiral motion. To shed more light on the relation between the flattened spiral motion and the vortex structures, the evolution of the streamwise vortices is presented in figure 12 during the time period $t^{\ast }=75{-}83$ . From the figure, the streamwise vortices evolve in a combined manner of the vortex shedding and the vortex twisting, that is, the vortex pairs are partially twisted with one another at the upper part of the threads while they are still periodically shed from the bubble. Therefore, the vortex threads in the flattened spiral stage have the concurrent structures within the zigzag and spiral motion. As a consequence of this particular vortex structure, the bubble travels within a flattened spiral trajectory as an intermediate state between zigzag and spiral motion.

Table 2. Characteristic properties of the single bubble motion under the flow condition of $Ga=35.36$ . From left to right: the top view of the rising path, the streamwise vorticity isocontours $\unicode[STIX]{x1D714}_{z}^{\ast }=\pm 0.4$ and the evolution of $\unicode[STIX]{x1D70F}_{plane}$ on a plane $4R$ downstream from the bubble centre. From top to bottom: $Eo=1.0$ , $Eo=2.0$ and $Eo=3.0$ . It is observed that, with increased $Eo$ , the rising path of the single bubble transits from zigzag to flattened spiral, and ultimately to spiral; meanwhile, the imbalance between the vortex pairs, given as $\unicode[STIX]{x1D70F}_{plane}$ , is also increased.

Appendix B. To trigger the spiral motion by imposing perturbations on the bubble shape

In this part, our objective is to further prove that through imposing perturbations on the bubble interface, by which the symmetries of the bubble shape and the vortex pairs are broken, the bubble motion would transit from zigzag to spiral. This is viewed as a supplement of § 4.2, in which the physical mechanisms for the path transition are discussed more deeply.

Table 3. Characteristic properties of the single bubble motion under the flow condition of $Ga=42.43$ . From top to bottom: $Eo=1.5$ and $Eo=2.0$ , corresponding to zigzag and flattened spiral motion, respectively. Other descriptions refer to those in table 2.

Figure 13. The characteristic properties of the bubble motion without and with imposing perturbations on the surface tension. From left to right: the projection of the rising path of the bubble and the evolution of the streamwise vorticity isocontours $\unicode[STIX]{x1D714}_{z}^{\ast }=\pm 0.4$ at different time periods; in addition, the streamwise vorticity distributions on the bubble interface at corresponding times are attached. It is observed, without perturbations, that the bubble rises within a zigzag path and the flow characteristics are symmetric; in contrast, such symmetries are broken when imposing the perturbations, and the bubble motion transits from zigzag to spiral gradually.

Table 4. Characteristic properties of the single bubble motion under the flow condition of $Ga=70.72$ . From top to bottom: $Eo=0.7$ and $Eo=3.0$ , corresponding to zigzag and spiral motion, respectively. Other descriptions refer to those in table 2.

The flow condition under investigation is $\{Ga,Eo\}=\{35.36,2.5\}$ , in which the bubble motion is thought to be spiral according to the phase map given by Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a ). Moreover, table 3 shows a path transition from zigzag to flattened spiral when $Eo$ increases from 1.5 to 2.0; therefore, a further increment to $Eo=2.5$ should lead to a spiral motion in sequence. In fact, within the computational time of $t^{\ast }=82$ , we just observe the bubble of $Eo=2.5$ to rise within a zigzag manner, as presented in figure 13 where ‘no perturbation’ is marked. For comparison, still in figure 13 where ‘perturbation’ is marked, small perturbations are imposed on the surface tension during the zigzag motion at $t^{\ast }=62$ for a short time period of $\unicode[STIX]{x0394}t^{\ast }=0.5$ , then the asymmetric bubble shape will be triggered, and, as a result, the double-threaded vortices will lose symmetry gradually. In the left column of figure 13, the projection of the bubble rising paths are displayed with and without perturbations, respectively, where the red (green) lines indicate the paths before (after) imposing the perturbations. In the right columns, the evolution of the streamwise vorticity isocontour $\unicode[STIX]{x1D714}_{z}^{\ast }=\pm 0.4$ is presented; the bubble shapes at the corresponding time points are also attached, so that the distribution of $\unicode[STIX]{x1D714}_{z}^{\ast }$ on the bubble interface are plotted.

Without imposing perturbations, it is observed that both the bubble shape and the vortex pairs retain symmetric structures throughout the computations; as a consequence, the bubble rises within a zigzag motion as we observe. The computation is stopped at $t^{\ast }=82$ due to the limitation of the computational resources. In contrast, after imposing perturbations on the surface tension during the numerical simulations, the bubble motion transits from zigzag to spiral gradually. To explore the physical mechanisms of such transition, we should still look into the evolutions of the bubble shape and vortex structures during this transition stage. It is observed that from $t^{\ast }=63$ to $72$ , the bubble shape and the distribution of $\unicode[STIX]{x1D714}^{\ast }$ on the interface develops to be asymmetric gradually, and the shedding double-threaded vortices tend to twist with one another. Subsequently, from $t^{\ast }=77$ to $92$ when the spiral motion is finally generated, the negative vorticities distributed on the bubble interface are dominant over the positive ones, and the shedding vortex pairs twist with one another completely. The evolution process of vorticity distribution during the transition stage is identical with table 1, and the only difference between them is which thread of the vortex pairs is stronger: positive one induces clockwise spiral motion while negative one induces anticlockwise motion.

Therefore, this part of the work contributes to strengthen the interpretations we give in § 4.2; it sheds more light on the relations between the bubble shapes, the vortex structures and the bubble motions.

References

Batchelor, G. K 2000 An Introduction to Fluid Dynamics. Cambridge University Press.Google Scholar
Blanco, A. & Magnaudet, J. 1995 The structure of the axisymmetric high-Reynolds number flow around an ellipsoidal bubble of fixed shape. Phys. Fluids 7 (6), 12651274.CrossRefGoogle Scholar
Brücker, C. 1999 Structure and dynamics of the wake of bubbles and its relevance for bubble interaction. Phys. Fluids 11 (7), 17811796.Google Scholar
Cano-Lozano, J. C., Bohorquez, P. & Martínez-Bazán, C. 2013 Wake instability of a fixed axisymmetric bubble of realistic shape. Intl J. Multiphase Flow 51, 1121.Google Scholar
Cano-Lozano, J. C., Martinez-Bazan, C., Magnaudet, J. & Tchoufag, J. 2016a Paths and wakes of deformable nearly spheroidal rising bubbles close to the transition to path instability. Phys. Rev. Fluids 1 (5), 053604.CrossRefGoogle Scholar
Cano-Lozano, J. C., Tchoufag, J., Magnaudet, J. & Martínez-Bazán, C. 2016b A global stability approach to wake and path instabilities of nearly oblate spheroidal rising bubbles. Phys. Fluids 28 (1), 014102.Google Scholar
De Vries, A. W. G., Biesheuvel, A. & Van Wijngaarden, L. 2002 Notes on the path and wake of a gas bubble rising in pure water. Intl J. Multiphase Flow 28 (11), 18231835.CrossRefGoogle Scholar
Duineveld, P. C. 1995 The rise velocity and shape of bubbles in pure water at high Reynolds number. J. Fluid Mech. 292, 325332.Google Scholar
Ellingsen, K. & Risso, F. 2001 On the rise of an ellipsoidal bubble in water: oscillatory paths and liquid-induced velocity. J. Fluid Mech. 440, 235268.CrossRefGoogle Scholar
Ern, P., Risso, F., Fabre, D. & Magnaudet, J. 2012 Wake-induced oscillatory paths of bodies freely rising or falling in fluids. Annu. Rev. Fluid Mech. 44, 97121.Google Scholar
Gaudlitz, D. & Adams, N. A. 2009 Numerical investigation of rising bubble wake and shape variations. Phys. Fluids 21 (12), 122102.CrossRefGoogle Scholar
Hua, J., Lin, P. & Stene, J. F. 2008 Numerical simulation of gas bubbles rising in viscous liquids at high Reynolds number. In Moving Interface Problems and Applications in Fluid Dynamics, Contemporary Mathematics, vol. 466, pp. 1734.Google Scholar
Leal, L. G. 1989 Vorticity transport and wake structure for bluff bodies at finite Reynolds number. Phys. Fluids A 1 (1), 124131.CrossRefGoogle Scholar
Leweke, T., Le Dizès, S. & Williamson, C. H. K. 2016 Dynamics and instabilities of vortex pairs. Annu. Rev. Fluid Mech. 48, 507541.Google Scholar
Lunde, K. & Perkins, R. J. 1997 Observations on wakes behind spheroidal bubbles and particles. In ASME FED Summer Meeting, Vancouver, Canada, Paper FEDSM, p. 141.Google Scholar
Magnaudet, J. & Eames, I. 2000 The motion of high-Reynolds-number bubbles in inhomogeneous flows. Annu. Rev. Fluid Mech. 32 (1), 659708.CrossRefGoogle Scholar
Magnaudet, J. & Mougin, G. 2007 Wake instability of a fixed spheroidal bubble. J. Fluid Mech. 572, 311337.Google Scholar
Moreau, R. J. 2013 Magnetohydrodynamics, vol. 3. Springer Science and Business Media.Google Scholar
Mougin, G. & Magnaudet, J. 2001 Path instability of a rising bubble. Phys. Rev. Lett. 88 (1), 014502.CrossRefGoogle ScholarPubMed
Mougin, G. & Magnaudet, J. 2006 Wake-induced forces and torques on a zigzagging/spiralling bubble. J. Fluid Mech. 567, 185194.Google Scholar
Popinet, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. J. Comput. Phys. 228 (16), 58385866.CrossRefGoogle Scholar
Sanada, T., Shirota, M. & Watanabe, M. 2007 Bubble wake visualization by using photochromic dye. Chem. Engng Sci. 62 (24), 72647273.Google Scholar
Shew, W. L., Poncet, S. & Pinton, J.-F.2005 Path instability and wake of a rising bubble Preprint, ENS, Lyon.Google Scholar
Shew, W. L., Poncet, S. & Pinton, J.-F. 2006 Force measurements on rising bubbles. J. Fluid Mech. 569, 5160.Google Scholar
Sommeria, J. & Moreau, R. 1982 Why, how, and when, MHD turbulence becomes two-dimensional. J. Fluid Mech. 118, 507518.Google Scholar
Tchoufag, J., Magnaudet, J. & Fabre, D. 2013 Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble. Phys. Fluids 25 (5), 054108.Google Scholar
Tchoufag, J., Magnaudet, J. & Fabre, D. 2014 Linear instability of the path of a freely rising spheroidal bubble. J. Fluid Mech. 751, R4.Google Scholar
Tripathi, M. K., Sahu, K. C. & Govindarajan, R. 2015 Dynamics of an initially spherical bubble rising in quiescent liquid. Nat. Commun. 6, 6268.Google Scholar
Veldhuis, C., Biesheuvel, A. & Van Wijngaarden, L. 2008 Shape oscillations on bubbles rising in clean and in tap water. Phys. Fluids 20 (4), 040705.Google Scholar
Yang, B. & Prosperetti, A. 2007 Linear stability of the flow past a spheroidal bubble. J. Fluid Mech. 582, 5378.Google Scholar
Zenit, R. & Magnaudet, J. 2008 Path instability of rising spheroidal air bubbles: a shape-controlled process. Phys. Fluids 20 (6), 061702.Google Scholar
Zenit, R. & Magnaudet, J. 2009 Measurements of the streamwise vorticity in the wake of an oscillating bubble. Intl J. Multiphase Flow 35 (2), 195203.Google Scholar
Zhang, J., Chen, L. & Ni, M.-J. 2017 The dynamic behaviors of the vortex interactions when two bubbles rising side by side in the viscous liquids. J. Fluid Mech.; (submitted).Google Scholar
Zhang, J. & Ni, M.-J. 2014a Direct simulation of multi-phase MHD flows on an unstructured Cartesian adaptive system. J. Comput. Phys. 270, 345365.Google Scholar
Zhang, J. & Ni, M.-J. 2014b Direct simulation of single bubble motion under vertical magnetic field: paths and wakes. Phys. Fluids 26 (10), 102102.Google Scholar
Zhang, J., Ni, M.-J. & Moreau, R. 2016 Rising motion of a single bubble through a liquid metal in the presence of a horizontal magnetic field. Phys. Fluids 28 (3), 032101.Google Scholar
Figure 0

Figure 1. Sketch of the physical model. If an external MF is imposed, it is along the gravity direction, given as $B_{z}$.

Figure 1

Figure 2. The different paths of the bubble ($Ga=50,~\unicode[STIX]{x1D70C}^{\ast }=1/1000,~\unicode[STIX]{x1D707}^{\ast }=1/100$): $Eo=1$ for zigzag (a), $Eo=2.0$ for zigzag–spiral transition (b) and $Eo=3.0$ for spiral (c).

Figure 2

Figure 3. (a) Time histories of the horizontal positions of the zigzag bubble ($Eo=1.0$), and the vortex structures within a zigzag period, respectively at $t^{\ast }=39$, $42$, $45$ and $48$. In the vortex contour maps, the light colour (green) corresponds to $\unicode[STIX]{x1D714}_{z}^{\ast }=0.5$ while the dark colour (blue) is $\unicode[STIX]{x1D714}_{z}^{\ast }=-0.5$. (b) The time histories of the norm of the streamwise vorticity accumulated on the bubble surface. (c) The time histories of the norm of the total vorticity accumulated on the bubble surface.

Figure 3

Figure 4. (a) The sketch of the vortex pairs behind the zigzag bubble with $Eo=1.0$. (b) The isosurface of $\unicode[STIX]{x1D714}_{z}^{\ast }=\pm 0.5$ on the plane of $4R$ downstream from the bubble centre. It is observed that they are well matched in strength. (c) The time series of the arithmetic integration of $\unicode[STIX]{x1D714}_{z}^{\ast }$ on the cut plane as $\unicode[STIX]{x1D70F}_{plane}=\int \unicode[STIX]{x1D714}_{z}^{\ast }\,\text{d}S$, proving that the counter-rotating vortex pairs are perfectly symmetric.

Figure 4

Figure 5. (a) Time histories of the horizontal positions of the zigzag–spiral bubble ($Eo=2.0$), and the evolution of the vortex structures during the transition stage from zigzag to spiral, respectively at $t^{\ast }=40$, $42$, $44$ and $46$. In the contour map, the light colour corresponds to $\unicode[STIX]{x1D714}_{z}^{\ast }=0.5$ and the dark colour is $\unicode[STIX]{x1D714}_{z}^{\ast }=-0.5$. (b) Time histories of $\Vert \unicode[STIX]{x1D714}_{z}^{\ast }\Vert$ accumulated on the bubble surface.

Figure 5

Figure 6. The $y\pm$ side views of the bubble shape when it travels in (a) spiral motion (asymmetric shape) and (b) zigzag motion (symmetric shape), respectively.

Figure 6

Table 1. The top row shows the streamwise vorticity distributions on the bubble surface at different time periods. The bottom row shows the contour maps of the maximum interface curvature at corresponding times. The values of the contour levels are attached in the colour scale. It is observed that, during the zigzag period from $t^{\ast }=35$ to $37$, both the interface curvatures and positive and negative vorticity distributions are symmetrical; however, in the transition stage from $t^{\ast }=41$ to $47$, both symmetries are broken, so that a larger surface curvature results in stronger positive vortices; and after entering the spiral motion at $t^{\ast }>53$, positive vortices completely dominate on the bubble surface because of the much larger local interface curvature.

Figure 7

Figure 7. (a) Sketch of the motion of the counter-rotating vortex pairs when they are unequal in strength. (b) Top view of the vortex structures behind the spirally rising bubble. (c) Contour map of the streamwise vorticity in the cross-plane $4R$ downstream from the bubble centre. (d) Time evolution of the arithmetic summation of $\unicode[STIX]{x1D714}_{z}^{\ast }$ on the cross-plane. It is observed that the value jumps from zero to a positive value around $t^{\ast }=40$, indicating that the inner positive vortex thread is dominant gradually in the transition stage.

Figure 8

Figure 8. The time histories of the rising velocities under different MFs. It is observed that, except for the relatively lower velocity without MF, there is nearly no difference among other cases when $N_{z}>0$.

Figure 9

Figure 9. The evolution of the arithmetic summation of $\unicode[STIX]{x1D714}_{z}^{\ast }$ over the the bubble surface (ad) and the projection of the rising paths of the bubble (eh) under vertical MFs. The MFs are introduced at $t^{\ast }=42$. It is observed that, with the increase of $N_{z}$, the positive and negative streamwise vorticities are more symmetrical on the bubble surface, resulting in a decreased radius of spiral motion.

Figure 10

Figure 10. Time evolution of the norm of $\Vert \unicode[STIX]{x1D714}_{z}^{\ast }\Vert$ accumulated on the the bubble surface when it rises under different vertical MFs. The MFs are introduced at $t^{\ast }=42$. It is observed that $\Vert \unicode[STIX]{x1D70F}_{z}\Vert$ decreases with the intensification of the vertical MFs. In particular, in the case $N_{z}=1.0$, where the vortex pairs are already symmetrical as shown in figure 9, it results in the rectilinear but not the zigzag path because $\Vert \unicode[STIX]{x1D714}_{z}^{\ast }\Vert$ is not large enough to trigger the vortex shedding.

Figure 11

Figure 11. The evolutions of the double-threaded vortex pairs with or without vertical MFs, respectively with $N=0$, $N_{z}=0.5$ and $N_{z}=1.0$. The time period is selected from $t^{\ast }=43$, when the MFs are just imposed onto the spirally rising bubble. In the contour map, the light colour corresponds to $\unicode[STIX]{x1D714}_{z}^{\ast }=0.5$ while the dark colour is $\unicode[STIX]{x1D714}_{z}^{\ast }=-0.5$. It is observed that in the case $N=0$ the vortex pairs twist with one another, and the structures are very stable. However, under $N_{z}=0.5$, the twisted structures seem to be weakened, but never stand parallel with one another. When the magnetic strength increases to $N_{z}=1.0$, the disintegration of the twisted vortex structures is quicker, and after $t^{\ast }>49$, there is nearly no streamwise vorticities shed from the bubble surface.

Figure 12

Figure 12. Time evolution of the streamwise vorticity isocontours $\unicode[STIX]{x1D714}_{z}^{\ast }=\pm 0.4$ during the flattened spiral motion. The flow conditions are $Ga=35.36$ and $Eo=2.0$. The time period is selected between $t^{\ast }=75$ and $83$. It is observed that the vortex structures behave as the combined manner of the vortex shedding and the vortex twisting.

Figure 13

Table 2. Characteristic properties of the single bubble motion under the flow condition of $Ga=35.36$. From left to right: the top view of the rising path, the streamwise vorticity isocontours $\unicode[STIX]{x1D714}_{z}^{\ast }=\pm 0.4$ and the evolution of $\unicode[STIX]{x1D70F}_{plane}$ on a plane $4R$ downstream from the bubble centre. From top to bottom: $Eo=1.0$, $Eo=2.0$ and $Eo=3.0$. It is observed that, with increased $Eo$, the rising path of the single bubble transits from zigzag to flattened spiral, and ultimately to spiral; meanwhile, the imbalance between the vortex pairs, given as $\unicode[STIX]{x1D70F}_{plane}$, is also increased.

Figure 14

Table 3. Characteristic properties of the single bubble motion under the flow condition of $Ga=42.43$. From top to bottom: $Eo=1.5$ and $Eo=2.0$, corresponding to zigzag and flattened spiral motion, respectively. Other descriptions refer to those in table 2.

Figure 15

Figure 13. The characteristic properties of the bubble motion without and with imposing perturbations on the surface tension. From left to right: the projection of the rising path of the bubble and the evolution of the streamwise vorticity isocontours $\unicode[STIX]{x1D714}_{z}^{\ast }=\pm 0.4$ at different time periods; in addition, the streamwise vorticity distributions on the bubble interface at corresponding times are attached. It is observed, without perturbations, that the bubble rises within a zigzag path and the flow characteristics are symmetric; in contrast, such symmetries are broken when imposing the perturbations, and the bubble motion transits from zigzag to spiral gradually.

Figure 16

Table 4. Characteristic properties of the single bubble motion under the flow condition of $Ga=70.72$. From top to bottom: $Eo=0.7$ and $Eo=3.0$, corresponding to zigzag and spiral motion, respectively. Other descriptions refer to those in table 2.