# Bubble breakup reduced to a one-dimensional nonlinear oscillator

Aliénor Rivière, Laurent Duchemin, Christophe Josserand, Stéphane Perrard

► To cite this version:

Aliénor Rivière, Laurent Duchemin, Christophe Josserand, Stéphane Perrard. Bubble breakup reduced to a one-dimensional nonlinear oscillator. Physical Review Fluids, 2023, 8 (9), pp.094004. 10.1103/PhysRevFluids.8.094004. hal-04232278

HAL Id: hal-04232278

https://hal.science/hal-04232278v1

Submitted on 18 Jan 2024

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Bubble break-up reduced to a 1D non-linear oscillator

Aliénor Rivière¹,*, Laurent Duchemin¹, Christophe Josserand², and Stéphane Perrard¹

¹PMMH, CNRS, ESPCI Paris, Université PSL, Sorbonne Université, Université de Paris, 75005, Paris, France and

²Laboratoire d'Hydrodynamique (LadHyX), UMR 7646 CNRS-Ecole Polytechnique, IP Paris, 91128 Palaiseau, France

(Dated: January 18, 2024)

Breaking dynamics of bubbles in turbulence produce a wide range of bubble sizes, which mediates gas transfer, in particular at the ocean/atmosphere interface. At the scales close to the stability limit of bubbles torn away by inertial forces, a typical geometry that induces bubble break-up is the uni-axial straining flow. In this configuration, the bubble shapes and their limit of stability have been studied theoretically and numerically near their equilibrium. Using numerical simulations, we investigate the bubble dynamics and break-up in such flows, starting from initial shapes far from equilibrium. We show that the break-up threshold is significantly smaller than the previous linear predictions and evidence that the break-up threshold depends on both the Reynolds number at the bubble size, and the initial bubble shape (ellipsoids). To rationalize the bubble dynamics and the observed thresholds, we propose a reduced model for the oblate/prolate oscillations (second Rayleigh mode) based on an effective potential that depends on the control parameters and the initial bubble shape. Our model successfully reproduces bubble oscillations, the maximal deformation below the threshold and the bubble lifetime above the threshold.

## I. INTRODUCTION

The evolution of bubbles and droplets in turbulent flows has important fundamental and practical applications. Bubbles drive low solubility gas exchanges, such as CO₂, at the ocean-atmosphere interface [1–4] and play a major role in nuclear reactors and chemical reactions [5]. Due to the inherent complexity of turbulent flows, identifying the key flow ingredients leading to break-up remains challenging both experimentally and numerically [6–8]. In inertial flows, bubble fate is primarily controlled by the ratio between the inertial forces and the capillary forces, namely the Weber number, We. A Weber number of order unity separates stable (low We) from unstable (large We) bubbles. The effect of viscosity, on its side is quantified by the Reynolds number that balances inertia with viscous effects. Experimental measurements of bubble break-up in turbulence have reported a broad range of critical Weber number, questioning the nature of the transition and the threshold definition. In simplified flow geometries, one can perform rigorous stability analysis, to understand the physical mechanism at play. Stagnation point flows, for instance, have been studied to model bubble deformations and break-up in turbulence [9, 10]. By investigating the stationary shapes and their linear stability, it has been shown that below a critical Weber number, $We_c^S$, a stable and an unstable stationary solutions coexist. At $We_c^S$, the two solutions merge and no stationary solution remains beyond $We_c^S$: any bubble will surely break [11–14]. However, this transition, called a saddle node bifurcation, only defines an upper bound for the critical Weber number $We_c$ that separates breaking from non breaking bubbles in a given experimental or numerical setup. In sub-critical transitions, the knowledge of $We_c^S$ is insufficient to predict the dynamics in realistic conditions since finite amplitude perturbation can lead to a state change well below the critical value of global stability loss. Such transitions have been evidenced for instance in parallel flows [15], open flows [16] or in spatially extended systems [17] such as dissipative solitons [18, 19], and in viscous drop break-up [20, 21]. A dynamical description of bubble deformations far from the stable states is therefore still lacking.

The relevance of stationary extensional flows for bubble break-up in turbulence has been brought back into the spotlight by recent experimental studies, showing that the turbulence is frozen during the break-up process [22], and that extensional flows are among the relevant geometries for bubble break-up [23]. In this paper, we thus investigate numerically the dynamics of a bubble in a uni-axial straining flow, starting from initial shapes far from their equilibrium position. We demonstrate that the initial condition strongly affects the bubble fate and evidence the sub-critical nature of the transition to break-up. We then characterize the effective critical Weber number, as a function of Reynolds number and initial bubble shapes (spheroids). We also show that the whole dynamics can be reduced to a peculiar non-linear 1D oscillator whose parameters depend on the initial bubble shape, the Weber number and the Reynolds number. Eventually, we successfully reproduce the bubble temporal deformations and the behavior close to the critical Weber number.

<footnote>* alienor.riviere@espci.fr</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. a) Scheme of a bubble at the center of a uni-axial straining flow. (Oz) is the axis of symmetry. Arrows show typical streamlines in the absence of bubble. b) Enlargement around a bubble experiencing its maximum deformation at Re = 400 and We = 7.3 ≈ We_c (Re = 400). The blue line denotes the bubble interface, black lines, isocontours of the stream function.</figcaption></figure>

## II. NUMERICAL SET-UP

We inject a spherical bubble of diameter D into a uniaxial straining flow $u(z, r) = Eze_z - \frac{1}{2}Ere_r$, where E is the typical shear amplitude, $e_z$ and $e_r$ are unit vectors of the axi-symmetric coordinate system $(z, r)$ (see figure 1). Density and viscosity ratios are 850 and 55 respectively, close to air-water ratios. Physically, this numerical experiment corresponds to a bubble quickly transported or suddenly submitted to a new straining region. Both phases are assumed incompressible and non condensible. The bubble dynamics is controlled by two dimensionless numbers, the Weber number

$$We = \frac{\rho E^2 D^3}{\gamma}, \qquad (1)$$

which compares inertia and capillary forces, and the Reynolds number

$$Re = \frac{ED^2}{\nu}, \qquad (2)$$

which compares inertia and viscous force, with $\rho$ and $\nu$ the liquid density and the kinematic viscosity respectively and $\gamma$ the surface tension between the two phases. One could use alternatively to the Reynolds number, the Ohnesorge number $\text{Oh} = \nu\sqrt{\rho/(D\gamma)}$ which thus does not include the straining rate but only the material and geometrical properties of the bubble. While the Ohnesorge number plays an important role in the final stage of a fluid break-up [24, 25], it is more pertinent in our configuration to use dimensionless numbers that involve the inertia of the flow. We perform numerical simulations of the Navier-Stokes equations for two incompressible fluids, using the free open source software Basilisk [26–29]. The interface is described using a geometric Volume Of Fluid method combined with a second-order sharp interface reconstruction and we use Adaptive Meshgrid Refinement (AMR). We consider a half domain to enforce left-right symmetry: the numerical domain is a square of size $L = 10R_0$, ie $z \in [0, L]$ and $r \in [0, L]$, with $R_0$ the initial bubble radius. We perform the numerical simulations in two steps, with a minimal grid size of $L/2^9$ and $L/2^{10}$ respectively, corresponding to 51.2 and 102.4 points per bubble equivalent radius. The rational of these numerical resolutions and the numerical convergence study can be found in Supplemental Material [30]. First we create a stationary stagnation point flow without a bubble, called precursor, using Dirichlet boundary conditions for the velocity field at $r = L$ (inflow) and Neumann condition at $z = L$ (outflow), and conversely for the pressure. Starting from a fluid at rest, the velocity field converges to the straining flow, with an error on the total kinetic energy density smaller than 1% of the theoretical value, $11\rho E^2 L^2/48$. Note that, since the straining flow is vorticity free, we can use the same setup to simulate both viscous and inviscid flows. We extract the final stationary state and inject a bubble of radius $R_0 = L/10$ at the stagnation point ($r = z = 0$), by changing only density and viscosity. The boundary conditions at the bubble interface are initially not fulfilled, nevertheless the code adapts in a few time
<figure><img src="image_3.png" /><figcaption>FIG. 2. a) Critical Weber number versus Re (crosses) with an error inferior to $10^{-2}$. The solid black line is the inviscid value $We_c^\infty$. The blue curve has expression $We_c^\infty \exp(-100/(We_c^\infty Re))$. Open circles and dashed line are $We_c^S(Re)$ and its inviscid limit as computed by Sierra-Ausin et al. [14] and Miksis [11] respectively. The inset shows the viscous correction to $We_c^\infty$. The dotted line follows $100/Re$. b) Critical Weber number versus the ellipsoid shape parameter $a_0$, for inviscid simulations. The blue curve is a polynomial fit of degree 2 with a maximum at $a_0 = 1$. Initial (black) and critical (red) shapes are represented for every $a_0$.</figcaption></figure>

steps to restore a solution, which corresponds to a much shorter time than all physical time scales considered in the following. This bubble injection method has been successfully used to study bubble dynamics in other flows [31, 32]. For each Reynolds number, ranging from 10 to 800, we vary the Weber number. We also perform series of inviscid simulations.

Figure 1b shows an enlargement around a bubble at $Re = 400$ and $We = 7.3$, experiencing its maximum deformation. Bubble interface is the blue line. The flow field, visualized by the isocontours of the stream function (black lines), smoothly goes around the bubble. We also evidence in this figure a recirculation air flow inside the bubble.

### III. AN INITIAL VALUE PROBLEM.

At low Weber number, the bubble first elongates and then relaxes to its equilibrium shape, either via damped oscillations (for $Re$ typically larger than 100), or monotonic relaxation ($Re < 50$). For sufficiently large $We$ the bubble elongates along the $z$ direction and breaks. We denote by $We_c$ the critical Weber number which separates breaking from non breaking configurations. We measure $We_c$ as a function of $Re$ using a bisection method. The result is shown in figure 2a. The critical Weber number converges to the inviscid value (solid line), $We_c^\infty$ at large Reynolds number, with a viscous correction, $We_c^\infty - We_c$, following $1/Re$ (inset plot). We observe that reducing $Re$ for a fixed $We$ allows the bubble to pass from a stable to an unstable configuration. Indeed, viscosity plays a destabilizing role through viscous shear at the interface. The static critical Weber number $We_c^S$ has previously been investigating by several authors from quasi-static deformations [11, 12] or linear stability analysis [14]. The recent computation of
<figure><img src="image_4.png" /><figcaption>FIG. 3. a) & b) Several typical temporal evolution of $x$ at Re = 400. Simulation data are in color, the model (4) is superimposed (black dotted line). a) An initially spherical bubble at low We. b) Three evolutions close to $We_c(a_0)$, for a stable sphere (blue line), an unstable sphere (grey line) and a stable ellipsoid (orange line, $a_0 = 0.8$).</figcaption></figure>

$We_c^S$ as a function of the Reynolds number from [14] is shown in figure 2a (open circles) together with the inviscid limit from [11] (dashed line). The critical Weber number $We_c$ we measure is significantly smaller than $We_c^S$. Indeed, in practice, the threshold $We_c^S$ would be observed for quasi-static deformations of bubbles, henceforth neglecting inertial effects. Starting from an initially spherical bubble, inertia cannot in fact be neglected. Above the break-up threshold, $We_c < We < We_c^S$, there still exists a stable shape surrounded by a finite basin of attraction, but the initial condition, i.e. the initial shape which is deformed compared to the stable shape, leads to the escape from this basin, and therefore, to break-up. The observed break-up transition is henceforth a sub-critical bifurcation. For such bifurcations, the response to an initial finite perturbation is dramatic, and the dynamics cannot be investigated using only linear stability analysis [16]. For viscous suspended drops in extensional flows at low Reynolds number, similar sub-critical break-ups have been evidenced experimentally and numerically [20, 21, 33, 34]. However, for bubbles, to the best of our knowledge, the sub-critical nature of the transition has not yet been reported.

To test the sensitivity to the initial conditions, we also consider ellipsoids of revolution of the same volume, $4/3\pi R_0^3$, with a local radius $R(\theta)$, in a inviscid flow (see figure 1 for the definition of $\theta$). The semi axis $a_0 = R(\pi/2)/R_0$ sets the whole initial shape from volume conservation, with prolate shapes corresponding to $a_0 < 1$ and oblate shapes to $a_0 > 1$.

The critical Weber number (fig. 2b) dramatically depends on the initial bubble shape, as expected for a sub-critical transition. The critical Weber number is maximum for the sphere, and decreases for both oblate and prolate shapes as the additional surface energy takes part in the break-up process. Near the maximum in $a_0 = 1$, we expect a quadratic dependency of $We_c$ with the distance to the sphere $|1 - a_0|$, as shown by the parabolic fit (blue line) in figure 2b.

To quantify the bubble deformation dynamics, we introduce the second Rayleigh mode of oscillation [35]:

$$x = 2 \int_{0}^{\pi/2} \frac{R(\theta)}{R_0} Y_2^0(\cos \theta) \sin \theta d\theta \quad (3)$$

where $Y_2^0$ is the spherical harmonics of principal number $\ell = 2$ and secondary number $m = 0$, corresponding to oblate-prolate oscillations at an angular frequency $\omega_2 = (12\gamma/(\rho R_0^3))^{1/2}$. The modes $\ell = 2$ are known to capture accurately most of bubble deformations [31, 36, 37]. Here, due to symmetries, only the mode $m = 0$ is present.

Figure 3 illustrates the various dynamics of mode 2. Far from $We_c$ the mode amplitude $x$ exhibits damped oscillations and converges to a finite value $x_\infty$, corresponding to a non spherical stable shape, as shown in 3a for an initially spherical bubble at Re = 400. The same behaviour is also observed for different initial conditions. Figure 3b illustrates the dynamics near the critical threshold. Slightly below the critical Weber number $We_c$, for an initially spherical bubble (blue curve), the amplitude first approaches a plateau with a maximum value, $x_{max}$ close to its critical value $x_c$ and eventually converges to a stable shape. Just above the threshold (grey curve), the dynamics is initially indistinguishable from the stable case, until the amplitude grows exponentially and finally decays abruptly, right before break-up. For a different initial shape, as illustrated in figure 3b with $a_0 = 0.8$ (orange curve), we observe the same behaviour, however the critical deformation $x_c$ at threshold increases. These curves are symptomatic of a sub-critical transition with a stable and an unstable equilibrium positions, in which the stability depends on both the control parameters (We and Re) and the initial conditions.
<figure><img src="image_5.png" /><figcaption>FIG. 4. a) Evolution of the potential V, defined in (4), with We for initially spherical bubbles at Re = 400. The range of explored x values have a more intense color. Stable equilibrium positions are denoted by red diamonds (red dashed lines in a and b). The maximum values, $x_{max}$, are denoted by black stars (same in 3a, 3b). b) Evolution of V with Re for a fixed of We = 5.</figcaption></figure>

IV. REDUCED MODEL

We assume that the whole bubble dynamics can be described by a damped non linear oscillator for x of the form:

$$ \ddot{x} + \lambda \dot{x} = -\nabla V(x, \dots), \qquad (4) $$

where $V(x, \dots)$ is an effective potential that may depend on all control parameters (We, Re and $a_0$). Time is made dimensionless using the mode angular frequency $\omega_2$. $\lambda = 20\sqrt{2/3}Oh$, with $Oh = \sqrt{We}/Re$, the Ohnesorge number, is the theoretical linear damping factor as computed by Lamb [38]. This theoretical expression perfectly captures dissipation in our complete dataset. The case of a harmonic potential was investigated by Kang and Leal [39]. Here, we look for a stationary polynomial potential, V, of degree 3, the minimum degree allowing to have two equilibrium positions. To do so we minimize, for every simulation, the mismatch between left and right hand side of equation (4), that is to say the quantity:

$$ I = \int (\ddot{x} + \lambda \dot{x} - p_0 - p_1 x - p_2 x^2)^2 dt \qquad (5) $$

by optimizing $p_0$, $p_1$ and $p_2$. This model perfectly describes the temporal evolution of x, both far from We$_c$ (see black dotted line in figure 3a) and close to We$_c$ (figure 3b). The effective potential V depends on both We, as illustrated in figure 4a for Re = 400, where blue curves correspond to We → 0 and yellow curves to We → We$_c$, as well as on Re (see figure 4b) and $a_0$ (data not shown). For a given Re, as Weber increases, the stable equilibrium (red diamonds) shifts to the right, in agreement with the literature [11, 14, 39]. Concomitantly, as We increases, the initial velocity, $\dot{x}_0 \propto \sqrt{We}$, increases and the energy barrier decreases, leading to the critical case where $x_{max} = x_c$ and We = We$_c$ (yellow curve). Increasing We further more would lead to the merging of the two equilibrium positions at We = We$_c^S$, corresponding to the global stability loss.

Conversely, for a fixed We, as Re decreases, the equilibrium position shifts to the positive values and the energy barrier decreases, as illustrated in figure 4b for We = 5. Both effects are consistent with the destabilizing role of viscosity which leads to a decrease of We$_c$ with Re, as can be seen on figure 2.
<figure><img src="image_6.png" /><figcaption>FIG. 5. Evolution of the coefficients of the three parameters defined in (5) with Re and We for initially spherical bubbles. Colored circles correspond to finite Re simulations and black crosses to inviscid simulations. a) Constant coefficient $p_0$ as a function of We. $p_0$ evolves linearly with We. b) Evolution of $p_0$ with Re. The dotted black line is the average inviscid value. We recover a more efficient forcing at small Re, compatible with a $1/Re$ scaling (solid black line). c) Evolution of the linear coefficient $p_1$ with Re. The black dotted line is the prediction of the pulsation from Kang and Leal [39] found from a linear development. d) Evolution of the quadratic coefficient $p_2$ with We.</figcaption></figure>

Figures 5 show the evolution of the three coefficients $p_0$, $p_1$ and $p_2$ of equation (5) with We and Re, for initially spherical bubbles. Circles denote finite Reynolds number simulations while black crosses are for inviscid simulations. For all Re, the constant forcing $p_0$, depends linearly on We (figure 5a), as was found theoretically in the inviscid case by Kang and Leal [39]. However, there was no theoretical prediction for the Re-dependency. Figure 5b shows that $p_0$ decreases with Re, in agreement with the destabilizing effect of viscosity and converges to its inviscid value (black dotted line). The shape is compatible with a $1/Re$ decay (solid black line), reminiscent of the Reynolds dependency of $We_c$. The linear coefficient $p_1$ is found to always be negative (figure 5c): the linear restoring force is positive, and $-p_1$ is the oscillator angular frequency. We found a weak dependency of $p_1$ with Re from 100 to 800. For the lowest Reynolds number $Re = 50$, the value of $p_1$ is under-determined, as the oscillations are overdamped. As We increases, the angular frequency $|p_1|$ decreases, as can be visualized on figures 3. For small We, the dynamics is well approximated by a linear oscillator, and we recover the theoretical prediction from Kang and Leal [39] (black dotted line). They show, by developing the dynamics around the equilibrium position, that the linear dependency of the pulsation with We comes from the coupling of mode 2 with the mode 4. The nonlinear coefficient $p_2$ is found to be always positive, in agreement with the existence of an unstable equilibrium position for $x$. $p_2$ is of order 1, and exhibits non monotonic evolution with the Weber number (figure 5d). However, $p_2$ values may be underdetermined
<figure><img src="image_7.png" /><figcaption>FIG. 6. a) Maximal deformation as a function of the distance to We$_{c}$. Finite Re simulations are denoted with circles, inviscid ellipsoids by triangles and inviscid spheres by diamonds. As We → We$_{c}$, the maximal deformation converges to its critical value D$_{c}$. Inset plot: Rescaled D$_{max}$ for the spheres, with two parameters that depend on Re: D$_{c}$ and α. In b) the black dotted line denotes the inviscid value for spheres (and correspond to a$_{0}$ = 1 in c). b) & c) Evolution of D$_{c}$ with Re and a$_{0}$ respectively. d) & e) Similar plots for the evolution of the slope with Re and a$_{0}$.</figcaption></figure>

for small and moderate values of We, for which the dynamics is mostly linear. We did not find a clear dependency of p$_{2}$ with Re.

V. DYNAMICS CLOSE TO THE CRITICAL POINT.

We now use this effective potential V description for the dynamics of x, to quantify the evolution of the maximal deformation as We → We$_{c}$. For the sake of simplicity, we consider the limit of negligible dissipation, in which energy is conserved:

$$
\frac{1}{2} \dot{x}_0^2 + V(x_0, \text{We}) = \frac{1}{2} \dot{x}^2 + V(x, \text{We}), \qquad (6)
$$

With the initial condition $\dot{x}_0 = \sqrt{2We}$. Without any loss of generality we set $V(x_0, We) = 0$ for all We. We define
$\max x = x_{\max}$. At this point $\dot{x}_{\max} = 0$, so that $x_{\max}$ is solution of:

$$
We = V(x_{\text{max}}, We). \tag{7}
$$
<figure><img src="image_8.png" /><figcaption>FIG. 7. a) Dimensionless lifetime, ET, as a function of the distance to We$_{c}$. Finite Re simulations are denoted with circles, inviscid ellipsoids by triangles and inviscid spheres by diamonds. As We → We$_{c}$, the lifetime diverges logarithmically. Inset plot: Rescaled ET for the spheres, with two parameters that depend on Re: ET(2We) and β. In b & c the black dotted line denotes the inviscid value for spheres (and correspond to a$_{0}$ = 1 in c & e). b) & c) Evolution of ET(2We) with Re and a$_{0}$ respectively. d) & e) Similar plots for the evolution of the slope β with Re and a$_{0}$.</figcaption></figure>

Since $\partial_x V|_{x_c} = 0$, developing (7) around the unstable position $x_c$, at the lowest orders in We$_{c}$ – We and $x_c - x_{max}$ gives:

$$ \mathrm{We}_c + (\mathrm{We} - \mathrm{We}_c) = V_c + (\mathrm{We} - \mathrm{We}_c)\partial_{\mathrm{We}}V|_c + \frac{1}{2}(x_c - x_{\mathrm{max}})^2\partial_{xx}V|_c \quad (8) $$

where $V_c = V(x_c, We_c)$ and $|_c = |_{x_c, We_c}$. Since, by definition, $We_c = V_c$ we get:

$$ (x_c - x_{\max})^2 = A(We_c - We) \quad (9) $$

with $A = 2(\partial_{We}V|_c - 1)/\partial_{xx}V|_c$. Figures 4a and 4b show that $\partial_{xx}V|_c$ is always negative, indeed $x_c$ is a maximum, and that $\partial_{We}V|_c < 0$. This ensures that $A > 0$ so that $x_c - x_{max}$ writes:

$$ x_c - x_{\max} = \sqrt{A}(We_c - We) \quad (10) $$

To give a simpler description of the shapes close to the critical point, we introduce the deformation parameter $D = 1 - R(\pi/2, t)/R(0, t)$ with $D < 0$ for oblate shapes, $D > 0$ for prolate shapes and $D \to 1$ for an infinitely long gas filament along $z$. Figure 6 shows $D$ for both the inviscid ellipsoids (triangles) and the spheres at finite Re (circles).
For $We \to 0$, the bubble is insensitive to the surrounding flow and $D_{max} \to D_0 = 1 - a_0^3$ when $a_0 \le 1$. Conversely, for $We \to We_c$ the maximum deformation converges to a critical value $D_c$, which depends on both Re and $a_0$. We find that both the initial bubble deformation and the Reynolds number increases the critical deformation $D_c$, with a slight difference between oblate and prolate shapes for the same distance to the sphere $|1 - a_0|$. The total deformation being mainly given by the amplitude $x$, we expect $D$ to follow:

$$D_c - D_{max} = \alpha \sqrt{1 - \frac{We}{We_c}}, \qquad (11)$$

as We approaches $We_c$. For each $Re$ and $a_0$ values, we fit the two parameters $\alpha$ and $D_c$ of equation (11). All data sets collapse on the same master curve as represented in the inset of figure 6a for sphere, showing that equation (11) also holds for non conservative systems (finite Re). Figures 6b and 6c show the evolution of $D_c$ with Re and $a_0$ respectively, while 6d and 6e present the evolution of the slope. In both 6b and 6d, the black dotted line (which corresponds to $a_0 = 1$ in 6c and 6e) represents inviscid values. $D_c$ increases with Re and decreases with the distance to the sphere. These dependencies are reminiscent of the evolution of $We_c$ with both Re and $a_0$. Indeed, since larger deformations need to be reached in order to break for larger Re or distance to the sphere, $We_c$ increases. The deformation $D_c$ varies strongly with $a_0$, which implies that bubble faith is highly dependent on bubble history. The critical shapes show that bubbles are more deformed at criticality when the initial shape is not spherical (see figure 2b), revealing the importance of inertial effect in the breakup process. We also conclude that there is no absolute maximum deformation after which the bubble breaks.

Similar developments can be performed to model the lifetime slightly above the critical Weber number. Indeed, when $We \to We_c$, the bubble lifetime is dominated by the time spent close to the unstable shape. We show in the Supplemental Material [30] that this time can be expressed as:

$$ET = ET(2We) - \beta \log(We/We_c - 1) \qquad (12)$$

where $ET(2We)$ and $\beta$ are two constants, and $We > We_c$. Figure 7a shows the dimensionless lifetime, $ET$, as a function of the distance to the critical point. In the limit of large Weber number, for all cases, the lifetime converges to the advection time $1/E$. Near $We = We_c$, the lifetime diverges logarithmically. After adjusting the two constants of equation (12) for each dataset, all the data collapse onto a single curve, as shown on the inset plot of figure 7a for initially spherical bubbles. Figures 7b and 7c show that $ET(2We) \approx 1$ a value that is independent on both Re and $a_0$. On the contrary, the slope $\beta$ increases as the initial shape gets away from the sphere (see figure 7e) and slightly increases with Re (7d).

# VI. CONCLUSION AND PERSPECTIVES

In this paper, we evidence that bubbles can break in a uniaxial straining flow, even when there still exists a stable equilibrium position. The threshold at which break-up occurs depends on both the Reynolds number and the initial bubble shape. Since, in real configurations bubbles dynamics are rarely quasi-static, these results have practical important consequences: history matters. The critical Weber number at which bubbles break should always be considered together with a set of initial conditions or at least understood in a statistical sense. In turbulent flows for instance, the probability that a bubble encounters a large pressure or velocity fluctuation that breaks it depends on its size, but all bubbles can break.

Taking advantage of the dynamical system approach, we show that bubble dynamics can be described by a simple one dimensional oscillator which depends on We, Re and on the initial bubble shape. This model successfully captures the maximum deformation and the lifetime close to critical conditions.

In a turbulent flow, a bubble will be immersed in a succession of various flow geometries of random duration. The relevance of the uniaxial strain flow has been shown by the experimental work of Masuk et al. [23], who measured the relative orientation between the bubble principal axis of deformation and the velocity gradient tensor at the bubble scale (see figure 5 of [23]). From their measurements, two main flow geometries were identified: the slip of a bubble with respect to the surrounding flow and the straining flow that elongates (respectively compresses) the bubble.

The persistence time of each flow geometry can be estimated, considering turbulent scaling laws for bubbles within the inertial range. A lower bound of the typical correlation time of velocity fluctuations at the bubble scale is given by the eulerian correlation time, namely the eddy turn-over time, $t_c(d) \sim \epsilon^{-1/3}d^{2/3}$, where $\epsilon$ is the energy dissipation rate. This time can be compared to the bubble capillary period $T_2 = \pi/(2\sqrt{6})\sqrt{\rho d^3/\gamma}$ of the dynamics. We have the relation $t_c/T_2 \propto 4\sqrt{3}/\pi We_t^{-1/2}$, where $We_t$ is the Weber number at the bubble scale, usually defined by $We_t = 2\rho\epsilon^{2/3}d^{5/3}/\gamma$. In practice, as long as $We_t < 4.8$, $t_c/T_2 > 1$, and the flow is correlated over more than one period of oscillation and can be considered as frozen. In this regime, corresponding to all bubbles near the critical stability threshold, the dynamics described in this article may hold. One main perspective of this work will then be
to model turbulence as a succession of stationary uni-axial straining flows of given orientation and strain rate. One could then evaluate the bubble deformation experienced in a given flow configuration, and iterate the process using a newly found initial shape and a randomly picked flow configuration until the dynamical critical Weber number is exceeded. This approach could be used to provide a versatile statistical framework for bubble break-up.

[1] W. K. Melville, The role of surface-wave breaking in air-sea interaction, Annual Review of Fluid Mechanics **28**, 279 (1996).

[2] G. B. Deane and M. D. Stokes, Scale dependence of bubble creation mechanisms in breaking waves, Nature **418**, 839 (2002).

[3] L. Deike and W. K. Melville, Gas transfer by breaking waves, Geophysical Research Letters **45**, 10,482 (2018), https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2018GL078758.

[4] L. Deike, Mass transfer at the ocean-atmosphere interface: the role of wave breaking, droplets, and bubbles, Annual Review of Fluid Mechanics **54**, 191 (2022).

[5] H. A. Jakobsen, *Chemical reactor modeling* (Springer, 2008).

[6] J. Vejrazka, M. Zednikova, and P. Stanovsky, Experiments on breakup of bubbles in a turbulent flow, A.I.Ch.E. Journal **1** (2017).

[7] S. Elghobashi, Direct numerical simulation of turbulent flows laden with droplets or bubbles, Annual Review of Fluid Mechanics **51**, 217 (2019), https://doi.org/10.1146/annurev-fluid-010518-040401.

[8] A. U. M. Masuk, A. K. Salibindla, and R. Ni, Simultaneous measurements of deforming hinze-scale bubbles with surrounding turbulence, Journal of Fluid Mechanics **910** (2021).

[9] J. Rodríguez-Rodríguez, J. Gordillo, and C. Martínez-Bazán, Breakup time and morphology of drops and bubbles in a high-reynolds-number flow, Journal of Fluid Mechanics **548**, 69 (2006).

[10] A. Revuelta, J. Rodriguez-Rodriguez, and C. Martínez-Bazán, Bubble break-up in a straining flow at finite reynolds numbers, Journal of Fluid Mechanics **551**, 175 (2006).

[11] M. J. Miksis, A bubble in an axially symmetric shear flow, The Physics of Fluids **24**, 1229 (1981).

[12] G. Ryskin and L. Leal, Numerical solution of free-boundary problems in fluid mechanics. part 3. bubble deformation in an axisymmetric straining flow, Journal of Fluid Mechanics **148**, 37 (1984).

[13] I. Kang and L. Leal, Numerical solution of axisymmetric, unsteady free-boundary problems at finite reynolds number. i. finite-difference scheme and its application to the deformation of a bubble in a uniaxial straining flow, The Physics of fluids **30**, 1929 (1987).

[14] J. Sierra-Ausin, P. Bonnefis, A. Tirri, D. Fabre, and J. Magnaudet, Dynamics of a gas bubble in a straining flow: Deformation, oscillations, self-propulsion, Phys. Rev. Fluids **7**, 113603 (2022).

[15] T. S. Chen and D. D. Joseph, Subcritical bifurcation of plane poiseuille flow, J. Fluid Mech. **58**, 337 (1973).

[16] J.-M. Chomaz, Global instabilities in spatially developing flows: non-normality and nonlinearity, Ann. Rev. of Fluid Mech. **37**, 357 (2005).

[17] P. Coullet, L. Gil, and D. Repaux, Defects and subcritical bifurcations, Phys. Rev. Lett. **62**, 2957 (1989).

[18] S. Fauve and O. Thual, Solitary waves generated by subcritical instabilities in dissipative systems, Phys. Rev. Lett. **64** (1990).

[19] L. Gordillo and M. A. Garcia-Nustes, Dissipation-driven behavior of nonpropagating hydrodynamic solitons under confinement, Phys. Rev. Lett. **112** (2014).

[20] H. A. Stone, B. Bentley, and L. Leal, An experimental study of transient effects in the breakup of viscous drops, Journal of Fluid Mechanics **173**, 131 (1986).

[21] H. A. Stone and L. G. Leal, Relaxation and breakup of an initially extended drop in an otherwise quiescent fluid, Journal of Fluid Mechanics **198**, 399 (1989).

[22] D. Ruth, W. Mostert, S. Perrard, and L. Deike, Bubble pinch-off in turbulence, PNAS (2019).

[23] A. U. M. Masuk, A. K. Salibindla, and R. Ni, The orientational dynamics of deformable finite-sized bubbles in turbulence, Journal of Fluid Mechanics **915**, A79 (2021).

[24] J. Eggers, Nonlinear dynamics and breakup of free-surface flows, Rev. Mod. Phys. **69**, 865 (1997).

[25] J. Burton, R. Waldrep, and P. Taborek, Scaling and instabilities in bubble pinch-off, Phys. Rev. Lett. **94** (2005).

[26] http://basilisk.fr, http://basilisk.fr.

[27] S. Popinet, Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries, J. Comput. Phys. **190**, 572 (2003).

[28] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, Journal of Computational Physics **228**, 5838 (2009).

[29] M. O. Abu-Al-Saud, S. Popinet, and H. A. Tchelepi, A conservative and well-balanced surface tension model, Journal of Computational Physics (2018).

[30] See Supplemental Material at for the numerical convergence analysis with and without the bubble, and for the derivation of equation (12). The Supplemental Material also contains Refs. [40].

[31] S. Perrard, A. Rivière, W. Mostert, and L. Deike, Bubble deformation by a turbulent flow, Journal of Fluid Mechanics **920**, A15 (2021).

[32] A. Rivière, W. Mostert, S. Perrard, and L. Deike, Sub-hinze scale bubble production in turbulent bubble break-up, Journal
of Fluid Mechanics **917**, A40 (2021).

[33] J. Rallison, A numerical study of the deformation and burst of a viscous drop in general shear flows, Journal of Fluid Mechanics **109**, 465 (1981).

[34] G. Gallino, T. M. Schneider, and F. m. c. Gallaire, Edge states control droplet breakup in subcritical extensional flows, Phys. Rev. Fluids **3**, 073603 (2018).

[35] L. Rayleigh *et al.*, On the capillary phenomena of jets, Proc. R. Soc. London **29**, 71 (1879).

[36] F. Risso and J. Fabre, Oscillations and breakup of a bubble immersed in a turbulent field, Journal of Fluid Mechanics **372**, 323 (1998).

[37] I. Kang and L. Leal, Bubble dynamics in time-periodic straining flows, Journal of Fluid Mechanics **218**, 41 (1990).

[38] H. Lamb, Hydrodynamics, Aufl., Cambridge: Univ. Press 1879-1932 **427** (1932).

[39] I. Kang and L. Leal, Small-amplitude perturbations of shape for a nearly spherical bubble in an inviscid straining flow (steady shapes and oscillatory motion), Journal of Fluid Mechanics **187**, 231 (1988).

[40] J. A. Van Hooft, S. Popinet, C. C. Van Heerwaarden, S. J. Van der Linden, S. R. De Roode, and B. J. Van de Wiel, Towards adaptive grids for atmospheric boundary-layer simulations, Boundary-layer meteorology **167**, 421 (2018).