# Instability of a swirling bubble ring

Yonghui Xu, Ivan Delbende, Daniel Fuster, Maurice Rossi

► To cite this version:

Yonghui Xu, Ivan Delbende, Daniel Fuster, Maurice Rossi. Instability of a swirling bubble ring. Comptes Rendus Mécanique, 2020, 348 (6-7), pp.519-535. 10.5802/crmeca.22. hal-02993273

HAL Id: hal-02993273

https://hal.science/hal-02993273

Submitted on 6 Nov 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Comptes Rendus
## Mécanique

Yonghui Xu, Ivan Delbende, Daniel Fuster and Maurice Rossi

**Instability of a swirling bubble ring**

Volume 1, issue 0 (0000), p. 000-000.

<https://doi.org/10.5802/crmeca.22>

© Académie des sciences, Paris and the authors, 0000.

Some rights reserved.

<figure><img src="image_1.png" /></figure>

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Thematic Issue / *Fascicule thématique*

# Instability of a swirling bubble ring

Yonghui Xu<sup>a</sup>, Ivan Delbende<sup>b, c</sup>, Daniel Fuster<sup>d</sup> and Maurice Rossi<sup>*</sup>,<sup>d</sup>

<sup>a</sup> Sorbonne Université, UMR 7190, Institut Jean Le Rond d'Alembert, 75005 Paris, France

<sup>b</sup> Sorbonne Université, UFR d'Ingénierie, 4 place Jussieu, 75005 Paris, France

<sup>c</sup> LIMSI, CNRS, Université Paris-Saclay, rue du Belvédère, 91405 Orsay, France

<sup>d</sup> CNRS, UMR 7190, Institut Jean Le Rond d'Alembert, 75005 Paris, France.

*E-mails*: yoghui.xu@sorbonne-universite.fr (Y. Xu), Ivan.Delbende@sorbonne-universite.fr (I. Delbende), daniel.fuster@sorbonne-universite.fr (D. Fuster), maurice.rossi@sorbonne-universite.fr (M. Rossi).

**Abstract.** A toroidal bubble or a cylindrical gas jet are known to be subjected to the Rayleigh-Plateau instability. Air bubble rings produced by beluga whales and dolphins however are observed that remain stable for long times. In the present work, we analyse the generation of such toroidal bubbles *via* numerical simulations, in particular how the process depends on surface tension. Their stability properties are then briefly analysed. For the estimated Reynolds and Weber numbers relative to the bubbles produced by these animals, the presence of a vortex inside and around the bubble is found to strongly stabilize the Rayleigh-Plateau instability.

**Keywords.** Instability, Vortex ring, Two-phase flow, Direct numerical simulation.

Manuscript received 30th June 2020, accepted 6th July 2020.

## 1. Introduction

This paper is written in *memoriam* Y. Couder. Yves was a colleague of deep thoughts and a humanist. In everyday life, he was prone to enthusiasm and passion for knowledge. He was an extraordinary observer and many of his breakthroughs originated from careful observations. As a modest homage, this paper presents a topic that takes its origin from an observation on Youtube! In his turbulence experiments [1, 2], Yves used bubbles to visualize vortices. The present work concerns a topic somewhat related: the generation and dynamics of an air bubble in the form of toroidal rings. It was some kind of surprise to see dolphins [3] or beluga whales [4] generate bubble rings from their mouth or blowholes and thereafter play for a long period of time with rings in which any sign of instability was absent.

<footnote>* Corresponding author.</footnote>
Indeed an axisymmetric columnar jet of fluid 1 embedded in a fluid 2, when slightly disturbed, is known to break into small droplets or bubbles. This subject was first studied by Lord Rayleigh [5, 6]. He computed how infinitesimal perturbations (more specifically normal modes of axial *k* and azimuthal *m* wavenumbers) evolve in time and obtained analytically the complex eigenvalue *s*(*k*, *m*) = *s*<sub>*r*</sub> + *is*<sub>*i*</sub>, and precisely *s*<sub>*r*</sub> their growth rate as a function of wavenumbers *k* and *m*. For an inviscid two-phase column of radius *a*<sub>0</sub>, axisymmetric perturbations *m* = 0 are the most dangerous ones and the dimensionless eigenvalue *s̄* is a function of the dimensionless wavenumber *k̄*:

$$ \bar{s}^2 \left( 1 + \frac{\rho^{(2)} K_0(\bar{k}) I'_0(\bar{k})}{\rho^{(1)} K_1(\bar{k}) I_0(\bar{k})} \right) = \bar{k} (1 - \bar{k}^2) \frac{I'_0(\bar{k})}{I_0(\bar{k})} \quad \text{with } \bar{s} \equiv \sqrt{\frac{a_0^3 \rho^{(1)}}{\sigma}} \text{, } \bar{k} \equiv ka_0. \quad (1) $$

For $0 \le \bar{k} \le 1$, $\bar{s}$ is real with the two possible signs, meaning that the jet is unstable. For $\bar{k} \ge 1$, $\bar{s} = \text{is}_i$ is purely imaginary indicating a stable jet. This result can be explained by the theory of minimal energy: a multiphase system tends to possess a minimum capillary energy at rest. This instability known as Rayleigh-Plateau instability is due to a capillary mechanism. It breaks the jet into several drops or bubbles when $\bar{k} \le 1$ so that the interface area is minimized. In order to stabilize Rayleigh-Plateau instability, an added rotation is an interesting possibility for applications, e.g. liquid atomization, spray generation, combustion processes. Previous works on the effect of rotation are those of [7–11]. In our case, we introduce this aspect to understand the video of dolphins playing with ring bubbles!

The present paper covers the three stages observed during the dolphin's video inside three separate sections. First, Section 2 introduces the generation of bubble rings through blowhole or mouth. To do so, one considers an idealized axisymmetric nozzle containing initially a gas bubble which is rapidly expelled. This part uses direct numerical simulation. A second stage in which the toroidal bubble is stretched is analyzed in Section 3 by analytical means. As an end result, this provides the vorticity profile reached after the stretching period. The third stage is discussed in Section 4 in which we study the stability of this new profile. The ring is replaced by a straight columnar bubble embedded in a vortex and a standard linear stability method is employed to investigate how superimposed infinitesimal perturbations evolve in time.

## 2. Two-phase ring generation

A bubble ring may appear in various ways: a first method consists in generating tiny gas bubbles in advance in a region which is thereafter crossed by a vortex ring produced away from this location [12]. The tiny bubbles migrate because of centrifugal force inside the vorticity ring core and then coalesce to form a toroidal gas bubble transported by the pre-formed vortex ring. In a second method, the vortex is generated at the same spot where tiny bubbles are located [13]. An other method consists in initializing a unique gas bubble inside a nozzle. Ejecting the fluid then causes the bubble to deform and break, yielding a toroidal bubble at the same time and location in which vorticity rolls up, e.g. near the trailing edge of the nozzle. This latter case is probably pertinent for vortex ring bubbles produced through the mouth or blowholes by belugas or dolphins. The presence of the tongue and its motion probably facilitate the toroidal bubble formation as mentioned and demonstrated by a scuba diver [14] but it is shown here not to be an essential feature.

We study below the toroidal bubble production through the mouth or blowholes in an idealized and simplest fashion without introducing the effect of the tongue. The geometry is an axisymmetric nozzle of radius $R_0$ and of thickness $2b_0 = R_0/4$ with an edge at 26.5° (see Figure 1a). The nozzle edge enables to fix the region of boundary layer separation and ends at $x = 2.8R_0$. Initially the two fluids are at rest and inside the tube, an air bubble of axial length $L_B$ is located
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Scheme of the numerical nozzle. (b) The whole numerical domain. In the yellow region, the mesh size is $\Delta_c = 0.125R_0$ and in the red region, $\Delta_c = 0.0078R_0$. (c) The three initial cases: the flow is at rest and a gas bubble (in red) of length $L_B$ is contained in the nozzle. From top to bottom: $L_B = 0.6R_0$, $1.2R_0$ and $2.4R_0$.</figcaption></figure>

within the axial interval $x \in [2.6R_0 - L_B, 2.6R_0]$ (see Figure 1c). The numerical domain (see Figure 1b) has been checked to be large enough to avoid boundary effects during the time of simulation: it is a square box $(x, r) \in [0, 16R_0] \times [0, 16R_0]$. Outflow boundary condition (stress free) are imposed at the downstream domain boundary $x = 16R_0$, impenetrability and slip condition on the lateral wall $r = 16R_0$ and no-slip everywhere else except for the region $x = 0$ and $r \le R_0$ where the velocity is imposed. The air bubble is suddenly expelled from the nozzle and injected inside the ambient water. This occurs *together with* the phenomenon of roll-up of the boundary layer vorticity. The velocity field far inside the nozzle (at $x = 0$) mimics this rapid ejection

$$u_x(x=0, r, t) = U_0 \operatorname{erf}\left(\frac{t}{\tau}\right) \operatorname{erf}(\eta), \quad u_r = 0 \quad \text{with} \quad \eta = \frac{R_0 - r}{\sqrt{v^{(w)}t}} \ge 0, \qquad (2)$$

where

$$\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} \exp(-s^2) \, ds$$

stands for the error function. The center velocity is changing from rest to a constant value $U_0$ within a time scale $\tau$. For this problem we assume that $\tau \ll R_0/U_0$: the center velocity reaches its constant value rapidly compared to the time of roll-up itself. In our simulations, this parameter is fixed at $\bar{\tau} = \tau U_0/R_0 = 0.1$. Moreover the flow is almost spatially uniform inside the tube except near the tube border $r = R_0$ where a boundary layer is present to ensure no-slip at the wall. This solution is close to what would be observed in a long tube subjected to a sudden pressure gradient. Kinematical viscosity of air and water are taken to be respectively (dimensional values are given in SI units) $v^{(a)} \sim 1.8 \times 10^{-5}$ and $v^{(w)} \sim 9 \times 10^{-7}$ and air and water density respectively $\rho^{(a)} = 1$, $\rho^{(w)} = 10^3$. Finally, air-water surface tension is assumed to be $\sigma = 73.4 \times 10^{-3}$.

Based on the nozzle radius $R_0$ as characteristic length scale and $U_0$ as the velocity scale, as well as water density and viscosity, this flow is defined by several dimensionless numbers: a *Reynolds* number, a *Weber* number and the relative length $r_B$ of the bubble:

$$Re = \frac{U_0 R_0}{v^{(w)}}, \quad We = \frac{\rho^{(w)} U_0^2 R_0}{\sigma}, \quad r_B = \frac{L_B}{R_0}. \qquad (3)$$

Note that $Re/We^{1/2} = (\rho^{(w)}\sigma R_0)^{1/2}/\mu^{(w)}$ is independent on velocity and is equal to $Re/We^{1/2} \approx 10^4 R_0^{1/2}$ when $R_0$ is expressed in SI units. The boundary layer thickness $\delta_0$ verifies $\delta_0/R_0 = \sqrt{\bar{\tau}/Re}$, where a nondimensional time variable $\bar{\tau} = U_0t/R_0$ is introduced. This thickness remains small for advective times much less than the Reynolds number $\bar{\tau} \ll Re$. The largest vorticity
<table><caption>**Table 1.** Estimated typical Reynolds and Weber numbers for scuba divers, beluga whales or dolphins</caption><thead><tr><th>SI units</th><th>R<sub>0</sub> = 0.02, U<sub>0</sub> = 0.1</th><th>R<sub>0</sub> = 0.1, U<sub>0</sub> = 0.1</th><th>R<sub>0</sub> = 0.02, U<sub>0</sub> = 0.3</th><th>R<sub>0</sub> = 0.1, U<sub>0</sub> = 0.3</th></tr></thead><tbody><tr><td>Re</td><td>2222</td><td>11111</td><td>6666</td><td>22222</td></tr><tr><td>We</td><td>2.72</td><td>13.6</td><td>24.5</td><td>120</td></tr></tbody></table>

occurs at the tube wall. In the bubble, for times such that $\bar{t}$ is of order one, the dimensionless vorticity maximum is of order

$$ \frac{R_0}{U_0} \max(\omega) \sim -\frac{2}{\sqrt{\pi}} \frac{\sqrt{Re}}{\sqrt{V^{(a)}/V^{(w)}}} \sim -\sqrt{\frac{Re}{15.7}} \quad (4) $$

For $Re = 5000$, the maximum is close to the value 17.8. Using dimensionless quantities, the entrance condition reads

$$ u_x = \text{erf}\left(\frac{\bar{t}}{\bar{t}}\right) \text{erf}(\eta), \quad u_r = 0, \quad \text{with} \quad \eta = \frac{1-r}{\sqrt{\bar{t}/Re}}. \quad (5) $$

The axisymmetric flow is described by the velocity fields $u_i(x, r, t)$ with $i = 1, 2$ related to the pseudo-scalar vorticity component $\omega_\theta(x, r, t)$, the plane $(x, r)$ being oriented by the normal unit vector $e_z$. In the present work, we focus on the transient period once the ejection has started: we hence study how the roll-up is modified by surface tension. For the formation process of the ring, axisymmetry is an adequate hypothesis but later on, one expects three-dimensional instabilities to appear on the starting bubble ring for a monophasic [15, 16] as for a two-phase vortex [15]. The numerical code we use in the present work is the Basilisk code [17]. It has been largely tested and used to simulate two-phase flow problems. One problematic aspect is the motion of the triple line between bubble, solid and liquid, present during the initial ejection stage. We will not go into numerical details here but let us simply mention *en passant* that Basilisk numerically mimics the presence near the triple line of a Navier condition, i.e. the presence of a slip length condition which is used in the phenomenology of triple line motion. Here the slip length is half the mesh size near the solid boundary, namely $\frac{1}{2}\Delta_c = 0.0039 R_0$.

The present work is not intended to be an extensive investigation of the parameter space of this problem, which in itself would cover several articles. We simply provide here some evidence about the generation of a bubble ring. Table 1 gives estimated Reynolds and Weber numbers for the toroidal bubble produced by scuba divers, beluga whales or dolphins. As a consequence, we use a realistic Reynolds number which is still affordable by numerical simulations running on a standard workstation. For most simulations, the Reynolds number is set at $Re = 5000$. In the monophasic case, the standard vortex is formed by vorticity sheet roll-up. In the two-phase situation, this mechanism is still active but the bubble expelled from the tube may break up in the same time interval. More importantly, the bubble contains and is surrounded by positive vorticity transported from the boundary layer at the nozzle wall but also from positive and negative vorticity produced at the interface. This interfacial vorticity (related in particular to capillary waves and bubble retraction) is shed in the bulk in the form of vortices. At later times, such vortices may induce a dynamics on the bubble as do the vortices generated by the nozzle wall vorticity roll-up.

The first set of simulations corresponds to a small Weber number $We = 3$. For a length $r_B = 0.6$ (Figure 2), the toroidal vortex is not formed but a bubble remains with vortices inside. This means that when surface tension is too strong, inertial effects are too weak to deform the bubble and thus to generate a toroidal bubble. The second simulation (Figure 3) corresponds to a larger width $r_B = 2.4$. In this instance, the bubble ring is again not formed but the starting period is quite complicated with phase of disconnection followed by reconnection. In this process, multiple
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Simulations at $Re = 5000$, $We = 3$ and $r_B = 0.6$. From left to right, the snapshots correspond to times $\bar{t} = 2, 3.2$ and $6.4$. The upper half displays phases and interfaces (red is gas and blue water) and the lower half the vorticity field. The vorticity value is colored with a maximum scale of $\omega = 3.75$, which is one fifth of the absolute vorticity maximum at the wall boundary layer (see (4)).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Simulations at $Re = 5000$, $We = 3$ and $r_B = 2.4$. From left to right, first row: $\bar{t} = 4$, $6, 8$; second row: $\bar{t} = 8.4, 10$ and $16$. Visualizations are similar to Figure 2.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Simulations at $Re = 5000$, $We = 100$ and $r_B = 0.6$. From left to right: $\bar{t} = 2.8, 4.8, 7.2$ and $15.2$. Visualizations are similar to Figure 2.</figcaption></figure>

vortices shed at the interface tend to favor this reconnection process. For larger bubbles, the starting vortex is shed inside the air bubble. As a conclusion, it can be checked that for small Weber numbers (e.g. $We = 3$), the bubble that is formed is not toroidal!

The second set displayed on Figures 4–5 corresponds to a large Weber number $We = 100$ for $r_B = 0.6$ and $r_B = 2.4$ respectively. In the first case ($r_B = 0.6$, Figure 4), the bubble becomes a thin layer near the centre and a rim forms near $r = R_0$ with the vorticity roll-up ($\bar{t} = 2.8$). The layer breaks thereafter into small fragments since surface tension is too weak and a toroidal bubble is formed near the edge attached to the nozzle edge ($\bar{t} = 4.8$). The flow then detaches the toroidal bubble ($\bar{t} = 7.2$) and the system evolves through capillary waves dynamics and interactions with
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Simulations at $Re = 5000$, $We = 100$ and $r_B = 2.4$. From left to right, first row: $\bar{t} = 4.8, 7.2, 8, 11.2$; second row: $\bar{t} = 13.6$ and $19.6$. Visualizations are similar to Figure 2.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Simulations at $Re = 5000$, $We = 100$ and $r_B = 8$. From left to right: $\bar{t} = 14, 20, 24$ and $40$. Visualizations are similar to Figure 2.</figcaption></figure>

vortices ($\bar{t} = 15.2$). For a larger $r_B$ (here $r_B = 2.4$, Figure 5), the boundary layer is not thin but the roll-up of the vortex sheet is still located near $r = R_0$ ($\bar{t} = 4.8$). At the center, the two interfaces get nearby ($\bar{t} = 7.2$) and then reconnect forming a unique toroidal bubble ($\bar{t} = 8$). Later the bubble is stretched ($\bar{t} = 11.2$) before being detached ($\bar{t} = 13.6$). The result is a larger toroidal bubble ($\bar{t} = 19.6$). The process is identical for $r_B = 8$ (Figure 6). For even larger values of $r_B$, the bubble does not break and the vorticity roll-up occurs inside the gas bubble. The largest toroidal bubble appears for a value of $r_B$ which is clearly dependent on the Weber number.

The intermediate case $We = 20$ is an in-between situation (Figures 7–9) in which a fine balance between the classical roll-up of vorticity and the bubble break-up occurs. For $r_B = 2.4$ (Figure 7), the same process occurs as at larger $We$: a toroidal bubble is formed, first attached to the edge ($\bar{t} = 9.6$) then detached ($\bar{t} = 14.4$). Thereafter its cross-section becomes almost circular ($\bar{t} = 20, \bar{t} = 24$). For $r_B = 1.2$ (Figure 8), the bubble is detached ($\bar{t} = 5.6$) before being broken ($\bar{t} = 9.6$). This is due to the presence of a vortex that was previously shed by the boundary layer. More generally, vorticity production at the interface makes the process quite complex since the toroidal bubble interacts with the vorticity generated at its interface beforehand, or generated by the wall boundary layer (see Figure 8 at times $\bar{t} = 9.6, 20$ and $30$). In many simulations for $We = 20$, the chaotic interaction between such vortices and the toroidal bubble brings this bubble towards the symmetry axis, merging it into a simple bubble. This process is seen here on a simulation for $r_B = 0.6$ (Figure 9): the toroidal bubble is formed at $\bar{t} = 4.8$ but its reconnects towards the center at $\bar{t} = 8$. This was also observed for $r_B = 2.4$ by slightly changing the Navier condition for the triple line.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Simulations at $Re = 5000$, $We = 20$ with $r_B = 2.4$. From left to right, first row: $\bar{t} = 4.4, 7.2, 9.6$; second row: $\bar{t} = 14.4, 20$ and 24. Visualizations are similar to Figure 2.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Simulations at $Re = 5000$, $We = 20$ with $r_B = 1.2$. From left to right, first row: $\bar{t} = 4.4, 5.6, 8, 9.6$; second row: $\bar{t} = 20, 30$. Visualizations are similar to Figure 2.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Simulations at $Re = 5000$, $We = 20$ with $r_B = 0.6$. From left to right: $\bar{t} = 2, 3.6, 4.8$ and 8. Visualizations are similar to Figure 2.</figcaption></figure>

It is worth mentioning that modifying the Reynolds number, as illustrated in Figure 10, and/or the position of the bubble may change some dynamical features but does not modify the overall picture described above. When the toroidal bubble is formed in the ejection problem, the structure of the vorticity field is complex inside the bubble and around it (see Figure 11). This is also seen on vorticity profiles (see Figure 12). One may define a bubble-ring total length $L_B$ which is of order $2\pi R_0$, as well as the size $a_B$ of the toroidal bubble section S and the mean circulation $\Gamma_B$ inside the bubble:

$$a_B = \sqrt{\frac{I_1}{\pi}} R_0 \quad \text{and} \quad \Gamma_B = I_2 U_0 R_0, \quad \text{with} \quad I_1 = \int_S dx dr, \quad I_2 = \int_S \omega dx dr. \quad (6)$$
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Simulations with $Re = 10000$ and $r_B = 2.4$. From left to right, top row: $We = 3$ at $\bar{t} = 8, 10, 14.6$; bottom row: $We = 100$ at $\bar{t} = 8, 10$ and $14.6$. Visualizations are similar to Figure 2.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Simulations with $Re = 5000$, $r_B = 4$. Vorticity field near the bubble for (a) $We = 20$ at $\bar{t} = 20$; (b) $We = 100$ at $\bar{t} = 26.4$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Profiles of vorticity $\omega_{\theta}(x, r)$ at $Re = 5000$, $We = 100$, $r_B = 4$ and $\bar{t} = 17.6$ (corresponding to the second case of Figure 11) on two orthogonal lines passing through the center of section S: (a) along the x direction, (b) along the r direction. The red dots are located inside the gas bubble, blue dots are outside.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Swirl number $q_s$ as a function of nondimensional time $\tilde{t}$. Simulations at $Re = 5000$, $We = 100$ with (a) $r_B = 4$, (b) $r_B = 5$.</figcaption></figure>

Based on such quantities, one computes the mean vorticity inside the bubble $\Gamma_B/(\pi a_B^2)$. For a uniform vorticity inside a bubble of radius $a_B$, the azimuthal velocity at the interface is equal to $\Gamma_B/(2\pi a_B)$. A dimensionless parameter called the swirl number $q_s$ can thus be defined as the ratio of this azimuthal velocity and the characteristic capillary velocity $U_{cap} \equiv [\sigma/(\rho^{(w)}a_B)]^{1/2}$:

$$q_s = \frac{\Gamma_B}{2\pi a_B U_{cap}}. \qquad (7)$$

After a transient period, this swirl seems to fluctuate around an average value, as illustrated in Figure 13.

### 3. The action of stretching on the bubble ring

In the videos [3, 4], dolphins or beluga whales play with the vortex rings and are able by a motion of their head to increase the bubble length. This is typically known as stretching mechanism of vortex tubes. Starting from a vorticity structure similar to the field obtained numerically at the end of the simulations in the previous section (that is those corresponding to vorticity fields in Figure 11), one should model the effect of the dolphin head. In order to remain simple, we analyse the effect of this stretching mechanism by a simplified model which enables us to use analytical methods. An air column of radius $a_B$ replaces the two-phase bubble ring with a vorticity field vanishing away from the bubble and the effect of stretching for the ring is introduced in the associated columnar bubble via an axial axisymmetric unsteady stretching

$$\mathbf{u} = [\gamma(t)x, u_y(y, z, t) - \frac{1}{2}\gamma(t)y, u_z(y, z, t) - \frac{1}{2}\gamma(t)z]. \qquad (8)$$

The vorticity structure at $t = 0$ is characterized by a unique vorticity component $\omega_x(y, z)$ along $x$ which is identical to $\omega_\theta(r, z)$ obtained at the end of previous simulations (see for instance Figure 11). The stretching occurs during a finite period of time between time $t = 0$ and $t = T_s$. The stretched solution is again characterized by the unique vorticity component $\omega_x(y, z, t)$. This field satisfies the governing equation

$$\frac{D\omega_x}{Dt} = \gamma(t)\omega_x + v\Delta_{2D}\omega_x, \qquad (9)$$

where $\Delta_{2D}$ stands for the two-dimensional Laplacian and $D \cdot /Dt$ for the material derivative:

$$\Delta_{2D} \equiv \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}, \quad \frac{D}{Dt} \equiv \frac{\partial}{\partial t} - \frac{1}{2}\gamma(t) \left[ y \frac{\partial}{\partial y} + z \frac{\partial}{\partial z} \right] + u_y \frac{\partial}{\partial y} + u_z \frac{\partial}{\partial z}. \qquad (10)$$
Across any point at the interface, the velocity field is continuous i.e. $u_y^{(1)} = u_y^{(2)}$ and $u_z^{(1)} = u_z^{(2)}$, and the tangential stress is continuous as well. Let us now introduce the following change of variables [18]

$$ \tau = \int_{0}^{t} S(t') dt', \quad \text{with } S(t) = \exp\left[\int_{0}^{t} \gamma(t') dt'\right], \quad \chi = \sqrt{S(t)}y, \quad \eta = \sqrt{S(t)}z, \quad (11) $$

as well as a rescaling of vorticity and velocity fields

$$ \tilde{\omega}_x(\chi, \eta, \tau) = \omega_x/S(t), \quad \tilde{\omega}_y(\chi, \eta, \tau) = u_y/\sqrt{S(t)}, \quad \tilde{u}_z(\chi, \eta, \tau) = u_z/\sqrt{S(t)}. \quad (12) $$

After such changes, the field $[\tilde{u}_y(\chi, \eta, \tau), \tilde{u}_z(\chi, \eta, \tau)]$ has the dynamics of an *unstretched solution* with the same initial condition:

$$ \frac{D\tilde{\omega}_x}{D\tau} = \nu\tilde{\Delta}_{2D}\tilde{\omega}_x, \quad \text{with } \tilde{\Delta}_{2D} = \frac{\partial^2}{\partial\chi^2} + \frac{\partial^2}{\partial\eta^2} \quad \text{and} \quad \frac{D}{D\tau} = \frac{\partial}{\partial\tau} + \tilde{u}_y\frac{\partial}{\partial\chi} + \tilde{u}_z\frac{\partial}{\partial\eta}, \quad (13) $$

associated to equivalent boundary conditions.

This trick allows one to obtain the effect of stretching on a vortex profile in two steps: first, one solves the pure two-dimensional advection-diffusion equation (13) starting with the velocity field at $\tau = t = 0$ and ending the simulation at $\tau = \tau_S = \int_0^{t=T_S} S(t') dt'$; second, one applies the scalings (11)–(12) back to the initial variables pertinent to describe the physical state. In the videos [3], the bubble ring length $L_B$ increases in a substantial manner: assume this increase to be by a factor of $F$ and $\gamma$ to be constant, this means that $F = \exp(\gamma T_S)$ and thus $\tau_S = \int_0^{t=T_S} S(t') dt'$ yields $(F-1)/\gamma$. Hence, the ratio $\tau_S/T_S = (F-1)/\ln F$ is large: for a factor $F = 10$, its value is around 4. The first stage of the simulation is thus a pure two-dimensional advection diffusion equation applied during a time $\tau_S$ larger than $T_S$. In the monophasic case, differential rotation in such an equation is known to accelerate diffusion [19], so that the vorticity profile rapidly tends towards a smoother axisymmetric one. More generally, two-dimensional vortical flows localized in a finite region are known to decay towards an axisymmetric monopole or, in case the total circulation is zero, towards a dipole [20, 21]. The present case is slightly different because of the presence of two phases. Since $v^{(a)}/v^{(w)} = 20$ is large, vorticity actually diffuses much more rapidly inside the bubble than outside. This implies that vorticity tends to become uniform inside the bubble. Vorticity outside becomes smoother than initially and the interface becomes axisymmetric by damping of capillary waves. This evolution is partially seen on Figure 14. In order to go back to the vortex solution at $t = T_S$ (second step), one applies the transformation inverse of (11)–(12). The stretching transformation leaves the vorticity uniform in the bubble, the section of which is now circular with radius $a_{new} \ll a(t=0)$. Indeed, by conservation of mass, $a_{new}/a(t=0)$ is of the order of $1/\sqrt{F}$. If we assume that the circulation $\Gamma_B$ inside the bubble is not changing much, this stretched vortex also possesses a new swirl number

$$ q_S^{\text{new}} \sim \sqrt{\frac{a(t=0)}{a_{\text{new}}}} q_S(t=0) \sim \sqrt{F} q_S(t=0). \quad (14) $$

The swirl number is thus expected to increase when the bubble ring length increases.

### 4. Instability of a columnar interface in the presence of surface tension and rotation

After the stretching phase, a toroidal gas bubble inside a vortex with an almost circular section is generated. In order to understand its stability, we proceed again by studying the same simplified problem: an air column of radius $a_0 = a_{new}$ replaces a two-phase bubble ring with a vorticity field vanishing away from the bubble. This simplified problem is much more tractable since the base flow is homogeneous in the axial direction $x$. Whenever necessary, the notation $Q^{(p)}$ is explicitly used to represent the quantity $Q$ in the inner phase $p=1$ or the outer phase $p=2$. The jump of a field $Q$ across a point $x$ of the interface is denoted by $[[Q]] = Q^{(1)}(x, t) - Q^{(2)}(x, t)$.
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Evolution of the vorticity field $\tilde{\omega}_x(\chi, \eta, \tau)$ near the bubble, starting from the initial condition corresponding to the second case in Figure 11 and evolving through the two-dimensional advection-diffusion equation (13), represented (left to right and top to bottom) at times $\tau = 0, 8, 16, 24, 32$ and $40$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Growth rate $\bar{s}$ as a function of wavenumber $\bar{k}$ for axisymmetric modes $m = 0$. (a) Inviscid column of liquid 1 surrounded by a liquid 2 at various density ratios $\rho^{(2)}/\rho^{(1)} = 10^3, 10^2, 10, 1, 10^{-1}$ and $10^{-3}$. (b) Viscous air column surrounded by water at various Ohnesorge numbers $Oh$ between $10^{-3}$ and $10^{-1}$ with $\mu^{(2)}/\mu^{(1)} = 50$. The inviscid curve almost coincides with that of case $Oh = 10^{-3}$.</figcaption></figure>

The theoretical relation (1) indicates that the presence of an external liquid does not change the stability criterion for various density ratios $\rho^{(2)}/\rho^{(1)}$ (Figure 15a). Here, the density ratio is $\rho^{(w)}/\rho^{(a)} = 10^3$, which simply makes the system of a columnar fluid 1 more stable: while the Laplace pressure pushes fluid 1 outward, fluid 2 will oppose a force, which slows down the instability and shift the most unstable wavenumber $\bar{k}$ to a smaller value [22]. The significance of dynamical viscosity $\mu^{(1)}$ (fixing the ratio $\mu^{(2)}/\mu^{(1)}$ to the value $\mu^{(w)}/\mu^{(a)} \approx 50$) compared to surface tension $\sigma$ is described by the dimensionless Ohnesorge number $Oh = \mu^{(1)}/(\rho^{(1)}\sigma a_0)^{1/2}$, a larger Ohnesorge number indicating an increased influence of viscosity. Since $\mu^{(1)}/(\rho^{(1)}\sigma)^{1/2} \approx 6.64 \times 10^{-5}$, $Oh \approx 6.64 \times 10^{-4}$ for a radius $a_0 = 1$ cm. Figure 15b presents the classical dispersion
relation for an air column surrounded by water when viscosity is taken into account. Again it is found that viscosity makes the system less unstable and shifts the most amplified wavelength to a longer on from $Oh = 10^{-2}$ on. The added inertia and viscosity are hence incapable to bring the system from an unstable to a stable configuration.

We assume here that an additional rotational motion stabilizes the toroidal bubble generated by dolphins: with an appropriate azimuthal speed, the centrifugal force stabilizes the Rayleigh-Plateau instability. The stability analysis for a swirling jet has been already studied by Hocking [7]. In the case of an inviscid fluid column with no external fluid and subjected to rotation of constant angular rotation rate $\Omega_0$ (or equivalently of constant vorticity $2\Omega_0$), a stability criterion for azimuthal wavenumber $m$ at $k = 0$ was obtained using the Weber number based on the velocity $\Omega_0 a_0$, i.e. $We = \rho^{(1)} a_0^3 \Omega_0^2 / \sigma$ (not to be confused with the Weber number used in (3) which was based on $U_0$ and $R_0$). This criterion implies that a mode with $k = 0$ and $m \neq 0$ is stable when $We \le m(m+1)$. The planar mode ($k = 0, m = 1$) being a displacement mode, it is neutral and stability for ($k = 0, m \neq 0$) modes is ensured when $We \le We_c = 6$. In follow-on studies, Gillis [10] provided a general stability criterion for three-dimensional disturbances on a viscous fluid column with no external fluid

$$We \le (ka_0)^2 + m^2 - 1. \qquad (15)$$

Weidman [11] showed how the dominant azimuthal mode for a two-phase axisymmetric rotating system depends on the Reynolds number $Re$. Finally the linear stability of a uniformly rotating viscous liquid column has been investigated by Kubitschek and Weidman [8] and validated experimentally [9]: the dominating mode depends on $Re$ and on the rotation speed. In the present work, we study how perturbations evolve when superimposed on an axisymmetric air column in an infinite ambient fluid, namely water. We consider the case where surface tension, viscosity, density contrast and centrifugal force act on the instability growth rate, and in which the angular velocity profile $U_\theta(r)/r$ is not uniform outside the bubble. To the best of our knowledge, this stability problem has not yet been performed.

### 4.1. The governing equation for a columnar flow with rotation

We adopt the model of an infinitely long two-phase capillary jet of radius $a_0$. The inner and outer fluid are both incompressible, immiscible and viscous. Using now radius $a_0$, $[\rho^{(1)} a_0^3 / \sigma]^{1/2}$ and $\rho^{(1)}$ as characteristic scales for length, time and density, we write the Navier-Stokes equations inside each fluid $p = 1$ or 2:

$$\bar{\nabla} \cdot \bar{\mathbf{u}}^{(p)} = 0 \qquad (16)$$

$$\bar{\rho}^{(p)} \frac{D}{D\bar{t}} \bar{\mathbf{u}}^{(p)} = -\bar{\nabla} \bar{\rho}^{(p)} + \frac{\bar{\mu}^{(p)}}{Re_B} \nabla^2 \bar{\mathbf{u}}^{(p)}, \quad \frac{1}{Re_B} = Oh = \frac{\mu^{(1)}}{\sqrt{\rho^{(1)} \sigma a_0}} \qquad (17)$$

where the Reynolds number $Re_B$ of the vortex is the inverse of the Ohnesorge number. Ratios $\bar{\mu}^{(p)} = \mu^{(p)} / \mu^{(1)}$, and $\bar{\rho}^{(p)} = \rho^{(p)} / \rho^{(1)}$ are fixed: $\bar{\mu}^{(2)} = 50$ and $\bar{\rho}^{(2)} = 10^3$ and by definition $\bar{\mu}^{(1)} = \bar{\rho}^{(1)} = 1$. In the following, we use only dimensionless quantities except when specified and the bar notation is assumed. Since the inner fluid $p = 1$ is air and the outer fluid $p = 2$ is water, it is clear that

$$Re_B = \sqrt{\frac{1}{\rho^{(2)}}} \sqrt{\frac{a_0}{R_0}} \frac{Re}{\sqrt{We}} \approx 300 \sqrt{a_0} \qquad (18)$$

when $a_0$ is the bubble radius in SI units. When $a_0 \sim 1$ cm, this implies $Re_B \sim 30$, and when $a_0 \sim 10$ cm, this implies $Re_B \sim 100$. Alternatively, if one assumes that $a_0 / R_0 \sim 0.1$, for simulations $Re = 10000$ and $We = 100$, one gets $Re_B = 10$; for simulations $Re = 10000$ and $We = 1$, one gets $Re_B = 100$.
The Navier-Stokes equations are written in cylindrical coordinates. The radial position of the interface between the two fluids can be described by a function $ζ$ in cylindrical coordinates: $r = ζ(θ, x, t)$. This interface is characterized by an outward normal $n$ and two tangential vectors $t_x$ and $t_θ$:

$$ t_x = \frac{\left(\frac{\partial \zeta}{\partial x}, 0, 1\right)}{\sqrt{1 + \left(\frac{\partial \zeta}{\partial x}\right)^2}}, \quad n = \frac{\left(1, -\frac{1}{\zeta} \frac{\partial \zeta}{\partial \theta}, -\frac{\partial \zeta}{\partial x}\right)}{\sqrt{1 + \left(\frac{\partial \zeta}{\partial x}\right)^2 + \left(\frac{1}{\zeta} \frac{\partial \zeta}{\partial \theta}\right)^2}}, \quad t_\theta = t_x \times n, \quad (19) $$

forming an orthonormal basis ($t_x, n, t_θ$). The velocity is continuous across the interface $[[u]] = 0$. The normal velocity must be compatible with the interface motion

$$ \frac{D\zeta}{Dt} = v_r(\zeta, \theta, x, t). \quad (20) $$

The dynamic conditions express that normal stress and tangential shear must be balanced by the surface tension effect:

$$ \left[ \left( -p + 2 \frac{\mu}{Re_B} \mathbf{n} \cdot \mathbf{e} \cdot \mathbf{n} \right) \right] = \frac{1}{R_1} + \frac{1}{R_2}, \quad \left[ \left( 2\mu t_\theta \cdot \mathbf{e} \cdot \mathbf{n} \right) \right] = 0 \quad \text{with} \quad e_{ij}^{(p)} = \frac{1}{2} \left( \frac{\partial u_i^{(p)}}{\partial x_j} + \frac{\partial u_j^{(p)}}{\partial x_i} \right). \quad (21) $$

## 4.2. The rotating base state

Here, we mainly study a basic state $u = W_θ(r, t)e_θ$ which possesses the feature described in Section 3: it is characterized by a non-uniform angular velocity $W_θ/r$ and a non-uniform axial vorticity $ω = [W_θ/r + \partial W_θ/\partial r]e_x$ in phase 2 only, and the interface is circular $ζ(x, θ, t) = 1$. If one neglects the unsteadiness due to bulk viscosity, such profile is a possible solution. In physical terms, this approximation is pertinent since instability occurs a priori much faster than the diffusion of the base state. To be mathematically correct, one introduces a body force to counterbalance the bulk diffusion for the basic state. However viscous effects are still present in the boundary conditions at $r = 1$:

$$ W_{\theta}^{(1)}(1) = W_{\theta}^{(2)}(1), \quad (22) $$

$$ \mu^{(1)} \left( \frac{\partial W_{\theta}^{(1)}}{\partial r}(1) - W_{\theta}^{(1)}(1) \right) = \mu^{(2)} \left( \frac{\partial W_{\theta}^{(2)}}{\partial r}(1) - W_{\theta}^{(2)}(1) \right). \quad (23) $$

This latter condition can be rewritten as well as

$$ \omega_x^{(2)}(1) = \frac{\mu^{(1)}}{\mu^{(2)}} \omega_x^{(1)}(1) - 2 \left( \frac{\mu^{(1)}}{\mu^{(2)}} - 1 \right) W_{\theta}^{(1)}(1). \quad (24) $$

In each phase ($p$), the basic pressure $\Pi(r)$ should satisfy equation $ρW_θ^2(r)/r = \partial\Pi/\partial r$ and a jump condition at $r = 1$. This leads to a pressure of the following form:

$$ \Pi^{(1)}(r) = \rho^{(1)} \int_{0}^{r} \frac{[W_{\theta}^{(1)}(r')]^2}{r'} dr', \quad \Pi^{(2)}(r) = \Pi^{(1)}(1) + \rho^{(2)} \int_{1}^{r} \frac{[W_{\theta}^{(2)}(r')]^2}{r'} dr' - 1. \quad (25) $$

The vorticity profile is assumed to be uniform in fluid 1 that is

$$ W_{\theta}^{(1)}(r) = qr, \quad \omega_{x}^{(1)} = 2q \quad \text{for } r < 1, $$

where the dimensionless parameter $q$ is equal to the ratio of the dimensional azimuthal velocity at the interface with respect to the characteristic capillary velocity $[\sigma/(\rho^{(1)}a_0)]^{1/2}$ — this quantity is equal to $0.27/a_0^{1/2}$ when $a_0$ is expressed in SI units. Note that

$$ q \equiv \sqrt{\frac{\rho^{(1)}}{\rho^{(2)}}} q_S = \sqrt{\frac{\rho^{(a)}}{\rho^{(w)}}} q_S. \quad (26) $$
<figure><img src="image_17.png" /><figcaption>**Figure 16.** (a) Normalized vorticity $\omega_x(r)/(2q)$ and (b) normalized velocity $W_\theta(r)/q$ profiles for the base state with $\rho^{(2)}/\rho^{(1)} = 10^3$, $\mu^{(2)}/\mu^{(1)} = 50$ for $\delta = 0.3$ (left), $\delta = 1$ (center) or $\delta = 2$ (right).</figcaption></figure>

Vorticity in phase 1 being uniform, Equation (24) imposes that $\omega_x^{(2)}(1) = 2q$. Apart from this constraint, let us assume a vorticity profile $\omega^{(2)}(r) = \omega_x^{(2)}(1)f(r)$ in phase 2, where function $f$ is chosen as $f(r) = \exp[-(r-1)^2/\delta^2]$ and decays in the radial direction, thus generating a vortex ribbon of radial size $\delta$ in phase 2 surrounding the interface. It is easily seen that

$$W_{\theta}^{(2)}(r) = \frac{q}{r} + 2\frac{q}{r}\int_{1}^{r} r'f(r')dr', \quad \omega_{x}^{(2)}(r) = 2qf(r) \quad \text{for } r > 1.$$

This expression, which can be written as well as

$$W_{\theta}^{(2)}(r) = \frac{q}{r} \left[ 1 + \delta^2 \left( 1 - e^{-(r-1)^2/\delta^2} \right) + \delta \sqrt{\pi} \operatorname{erf} \left( \frac{r-1}{\delta} \right) \right] \quad \text{for } r > 1, \quad (27)$$

is displayed in Figure 16.

### 4.3. The linear instability results

We study the linear stability of the above base flow using standard but cumbersome methods. Indeed, the problem is an eigenvalue problem $\bar{s}(\bar{k}, m)$ for given axial $\bar{k}$ and azimuthal $m$ wavenumbers, that contains viscous diffusion in the bulk as well as on the interface. A code has been specifically developed, that solves coupled Orr–Sommerfeld equations in cylindrical coordinates for the two phases. For obvious reasons, details on the implementation will be published elsewhere, and only results for the air-water columnar bubble are presented here.

Figure 17 displays the growthrate $\bar{s}$ as a function of wavenumber $\bar{k}$ for the axisymmetric mode $m=0$, which is known to be the most unstable mode for $q=0$. The presence of rotation is found to stabilize axisymmetric modes, and for $q$ larger than a critical value $q_c \approx 0.03$, waves become stable for all wavelengths. This observation holds irrespective of the Reynolds number and $\delta$ within the studied interval $\delta \in [0.3, 2]$. Actually the structure in the corona at $r > 1$ seems to play a minor role for the axisymmetric mode.
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Growth rate $\bar{s}$ as a function of wavenumber $\bar{k}$ for axisymmetric mode $m = 0$ at various swirl numbers $q = 0, 0.01, 0.02, 0.03$. Reynolds number is (a) $Re_B = 10$, (b) $Re_B = 100$. Viscous air column surrounded by water with $\delta = 0.3$ (left) or $\delta = 2$ (right), with $\rho^{(2)}/\rho^{(1)} = 10^3$ and $\mu^{(2)}/\mu^{(1)} = 50$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Growth rate $\bar{s}$ as a function of wavenumber $\bar{k}$ for a viscous column of air surrounded by water at $\delta = 0.3$, $Re_B = 10$ and various swirl $q = 0.2, 0.4, 0.6$. (a) $m = \pm 1$; (b) $m = \pm 2$; (c) $m = \pm 3$. Positive $m$: symbols; negative $m$: bold line.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Same as Figure 18 but for $\delta = 1$ ($Re_B = 10$).</figcaption></figure>

Figures 18–21 provide linear instability results for non axisymmetric modes $m \pm 1, \pm 2$ and $\pm 3$. The first observation is that, at given $\bar{k}$ and $|m|$ values, the growth rate is found almost the same
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** Same as Figure 18 but for $\delta = 0.3$ and $Re_B = 100$.</figcaption></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 21.** Same as Figure 18 but for $\delta = 1$ and $Re_B = 100$.</figcaption></figure>

for $m > 0$ and $m < 0$. The equality exactly holds when $q = 0$ for symmetry reasons, and the fact that it holds also roughly for $q \neq 0$ can be explained by the uniform vorticity inside the bubble. Interestingly, rotation does not influence mode $|m| = 1$, that remains stable. For sufficiently low Reynolds number and small $\delta$, other azimuthal modes $m \neq 0$ for small $k$ can be destabilized. For instance, the mode $m = 2$ is found unstable for $Re_B = 10$ and $\delta = 0.3$ (Figure 18b). The growth rate decreases as the Reynolds number increases (Figure 20b). This is reminiscent of planar modes $k = 0$ in Ref. [8]. Note that for larger values of $\delta$, this mode is restabilized: for $\delta = 1$, all modes tested are stabilized by rotation, this conclusion holding irrespective of the Reynolds number (Figures 19 and 21).

Finally, going back to the bubble ring, the axisymmetric Rayleigh-Plateau modes are stabilized by rotation. Indeed, from the simulations in Section 2, the pertinent values of $q$ (Figure 13 and (26)) are found to be $q \sim 10^{-3/2} \times 30 \sim 0.94$. When a stretching factor is applied, $q$ generally increases (Equations (14) and (26)) and is thus larger than the stability threshold of $q_c \approx 0.03$.

**5. Conclusions**

This paper covers three different aspects of the dynamics of toroidal bubbles such as those produced by scuba divers, beluga whales or dolphins. During the initial production of such bubble, vorticity trailing from the nozzle or from the opening rolls up into a toroidal vortex ring able to trap gas in its core; to do so, it has to overcome surface tension which would favor the formation of a single spherical bubble. Simulations show that vortices shed from the nozzle or generated at the bubble interface itself can influence the formation dynamics in a non trivial way. Once formed, the bubble ring may be stretched by an appropriate flow so as to increase its length. A theoretical model where curvature effects are overlooked shows that, during the stretching process, the vorticity in the gas tends to become uniform, the cross-section tends to circular while the core radius is highly reduced and the swirl is enhanced. This latter finding is important since a linear stability study of a gas columnar vortex surrounded by water eventually shows that the dominant axisymmetric Plateau instability is stabilized as swirl increases above a
well-defined critical value. Put all together, these pieces of investigation give some serious clues
to physically understand the surprising stability of swirling bubble rings.

References

[1] Y. Couder, "Observation expérimentale de la turbulence bidimensionnelle dans un film liquide mince", *C. R. Acad. Sci. Paris II* **297** (1983), p. 641-645.

[2] O. Cadot, S. Douady, Y. Couder, "Characterization of the low pressure filaments in three-dimensional turbulent shear flow", *Phys. Fluids* **7** (1995), p. 630-646.

[3] Youtube, last seen 2020. https://www.youtube.com/watch?v=2m6ie3MVIAw.

[4] Youtube, last seen 2020. https://www.youtube.com/watch?v=Iiuv2ELttaQ.

[5] Lord Rayleigh, "Xvi. On the instability of a cylinder of viscous liquid under capillary force", *Lond. Edinburgh Dublin Phil. Mag. J. Sci.* **34** (1892), no. 207, p. 145-154.

[6] Lord Rayleigh, "Xix. On the instability of cylindrical fluid surfaces", *Lond. Edinburgh Dublin Phil. Mag. J. Sci.* **34** (1892), no. 207, p. 177-180.

[7] L. M. Hocking, D. H. Michael, "The stability of a column of rotating liquid", *Mathematica* **6** (1959), no. 1, p. 25-32.

[8] J. P. Kubitschek, P. D. Weidman, "The effect of viscosity on the stability of a uniformly rotating liquid column in zero gravity", *J. Fluid Mech.* **572** (2007), p. 261-286.

[9] J. P. Kubitschek, P. D. Weidman, "Helical instability of a rotating viscous liquid jet", *Phys. Fluids* **19** (2007), no. 11, 114108.

[10] J. Gillis, B. Kaufman, "The stability of a rotating viscous jet", *Q. Appl. Math.* **19** (1962), no. 4, p. 301-308.

[11] P. D. Weidman, M. Goto, A. Fridberg, "On the instability of inviscid, rigidly rotating immiscible fluids in zero gravity", *Z. Angew. Math. Phys.* **48** (1997), no. 6, p. 921-950.

[12] Youtube, last seen 2020. https://www.youtube.com/watch?v=5KdsVsiZ6bU.

[13] Youtube, last seen 2020. https://www.youtube.com/watch?v=4mbu2ueMe2E.

[14] Youtube, last seen 2020. https://youtu.be/vFbN31EuMvc.

[15] P. G. Saffman, *Vortex Dynamics*, Cambridge University Press, Cambridge, 1992.

[16] T. T. Lim, T. B. Nickels, "Vortex rings", in *Fluid Vortices* (S. I. Green, ed.), Fluid Mechanics and Its Applications, vol. 30, Springer, Dordrecht, 1995, p. 95-153.

[17] S. Popinet, "A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations", *J. Comput. Phys.* **302** (2015), p. 336-358.

[18] J. D. Gibbon, A. S. Fokas, C. R. Doering, "Dynamically stretched vortices as solution of the 3D Navier-Stokes equations", *Physica D* **132** (1999), p. 497-510.

[19] A. J. Bernoff, J. F. Lingevitch, "Rapid relaxation of an axisymmetric vortex", *Phys. Fluids* **6** (1994), no. 11, p. 3717-3723.

[20] J. C. McWilliams, "The emergence of isolated coherent vortices in turbulent flow", *J. Fluid Mech.* **146** (1984), p. 21-43.

[21] R. Robert, J. Sommeria, "Statistical equilibrium states for two-dimensional flows", *J. Fluid Mech.* **229** (1991), p. 291-310.

[22] J. Eggers, E. Villermaux, "Physics of liquid jets", *Rep. Prog. Phys.* **71** (2008), p. 16-18.

[23] S. Popinet, Basilisk, http://basilisk.fr.

<figure><img src="image_27.png" /></figure>