# Bubble deformation by a turbulent flow

Stéphane Perrard1 Aliénor Rivière1,2 Wouter Mostert2 Luc Deike2,3 1LPENS,
Département de Physique, Ecole Normale Supérieure, PSL University, 75005
Paris, France
2Department of Mechanical and Aerospace Engineering, Princeton University
3 Princeton Environmental institute, Princeton University

(; revised ?; accepted ?. - To be entered by editorial office)

###### Abstract

We investigate the modes of deformation of an initially spherical bubble
immersed in an homogeneous and isotropic turbulent background flow. We perform
direct numerical simulations of the two-phase Navier-Stokes equations and
resolving the bubble deformation in an homogeneous and isotropic turbulent
flow at various Weber number (the ratio of turbulent and surface tension
forces), and discuss a theoretical framework for the bubble deformation in a
turbulent flow using a spherical harmonic decomposition. We derive at the
linear order an harmonic equation for each mode of bubble deformation, with a
forcing term originating from the background turbulent fluctuations, formally
relating the bubble deformation and the turbulent velocity fluctuations. The
growth of the total surface deformation and growth of each individual modes is
computed from the direct numerical simulations using an appropriate Voronoi
decomposition of the bubble surface. We show that two successive temporal
regimes occur: the first regime corresponds to deformations driven only by
inertial forces, with a linear growth of the interface deformation with time,
in good agreement with the theoretical predictions, whereas the second regime
results from a balance between inertial forces and surface tension, the
transition time between the two regimes being given by the period of the first
Rayleigh mode of bubble oscillation. We discuss how our theoretical framework
can be used to relate the bubble lifetime with the turbulence statistics and
eventually show that at high Weber number, bubble lifetime and break-up
characteristics can be deduced from the statistics of turbulent fluctuations
at the bubble scale.

##  1 Introduction

The interaction of turbulent flow with a free surface occurs in numerous
physical systems, from ocean waves forced by the turbulent wind (Phillips,
1957; Perrard et al., 2019), river surface patterns driven by underwater
turbulence (Peregrine, 1976; Brocchini & Peregrine, 2001), as well as drops
and bubbles in turbulent flow (Balachandar & Eaton, 2010; Elghobashi, 2019;
Mathai et al., 2020).

Deformation and break-up of bubbles and droplets in a turbulent flow control
exchanges of heat, mass and momentum in numerous natural and engineering
processes, from bubble-mediated gas exchange at the ocean-atmosphere surface
(Deike & Melville, 2018; Reichl & Deike, 2020), chemical reactors (Risso,
2018), the fall of rain drops (Villermaux & Bossa, 2009), the evaporation of
sea spray in turbulent boundary layer (Veron, 2015), dynamics of water
droplets in clouds (Balachandar & Eaton, 2010) as well as industrial liquid
atomization and fragmentation (Eggers & Villermaux, 2008).

The dynamics of a single bubble evolving freely in low-viscosity liquid at
rest has been extensively studied, with analytical and experimentalf results
describing bubble oscillations (Prosperetti, 1980; Miller & Scriven, 1968),
rise velocity (Moore, 1965; Maxworthy et al., 1996) and path instability
(Magnaudet & Eames, 2000; Ern et al., 2012). The rise dynamics are altered
once bubbles are close enough to each other and can interact through their
wakes (Harper, 1970; Yuan & Prosperetti, 1994), which at high void fraction
can interact through their wakes, a process that leads to collective effects
and a complex bubble induced turbulence dynamic (Lance & Bataille, 1991;
Risso, 2018).

The deformation dynamics and potential break-up of a bubble depend primarily
on the ability of the surrounding fluid to deform the bubble against surface
tension forces. This defines the Weber number, comparing inertial forces
generated by the turbulent carrier flow and the capillary cohesive forces.
Considering the velocity fluctuation at the bubble diameter scale $d_{0}$,
formalized by the longitudinal velocity increment $\delta u({\bf
d_{0}})=u_{L}({\bf r},t)-u_{L}({\bf r}+{\bf d_{0}},t)$, where $u_{L}$ is the
velocity component along the direction of $\bf d_{0}$, the turbulent Weber
number is defined as $\textrm{We}=\rho_{\ell}\langle\delta u(d_{0})^{2}\rangle
d_{0}/\gamma$ (Hinze, 1955; Risso & Fabre, 1998) with $\rho_{\ell}$ the
density of water, $\gamma$ the air-water surface tension and $\langle\rangle$
the average over the flow configurations. In an homogeneous and isotropic
turbulent flow, the velocity fluctuations at the bubble scale $\delta
u(d_{0})^{2}$ can be related to the mean dissipation rate of energy $\epsilon$
using Kolmogorov (1941) theory $\langle\delta u(d_{0})^{2}\rangle=C(\epsilon
d_{0})^{2/3}$ for $d_{0}$ in the inertial range. Experimental studies have
observed $C\in[2,2.2]$ depending on Reynolds number (Cowen & Variano, 2008;
Pope, 2000). We chose $C=2$ for consistency with Risso & Fabre (1998), and the
Weber number eventually writes:

$\displaystyle\textrm{We}=\frac{2\rho_{\ell}\epsilon^{2/3}d_{0}^{5/3}}{\gamma},$
(1)

for a bubble of diameter $d_{0}$ immersed in an homogeneous and isotropic
turbulent flow. Experimental studies of bubble dynamics in turbulence have
identified a critical Weber number $\textrm{We}_{c}$ of order unity, above
which the occurence of break-up becomes statistically dominant. The value of
$\textrm{We}_{c}$ reported from laboratory experiments varies among authors in
the range [1,5], corresponding to variation in experimental conditions, which
introduces other flow parameters such as large scale shear or spatial
variations in the dissipation rate (Martinez-Bazan et al., 1999; Andersson &
Andersson, 2006; Ravelet et al., 2011; Vejražka et al., 2018). Note also that
the critical Weber number is defined statistically, and can be influenced by
the temporal and spatial windows of observation of the bubbles. In this
article, we consider $\textrm{We}_{c}=3$, in accordance with our numerical
dataset, obtained from ensemble of simulations. Two mechanisms driving the
deformation and break-up have been discussed (Martinez-Bazan et al., 1999;
Andersson & Andersson, 2006; Ravelet et al., 2011; Vejražka et al., 2018),
either direct strong action of an eddy at the scale of the bubble leading to
large deformation and break-up, or resonance mechanism between deformation
caused by weaker eddies and oscillations of the bubble (Risso & Fabre, 1998).
Experimental studies have identified an oscillatory response of bubbles in
turbulence associated to the second eigenmode (Risso & Fabre, 1998; Ravelet et
al., 2011).

In such configurations, the gas-liquid interface is surrounded by stochastic
turbulent flow, characterized by perturbations of various strengths at various
scales. In this paper, we present a theoretical framework of bubble
oscillations to describe the temporal evolution of the bubble interface
deformation immersed in a turbulent background flow. To do so, we consider
linearized deformations and we neglect any feedback of the interface on the
statistics of the turbulent background flow. Our framework has two main
inspirations. On one hand, it extends the spherical harmonic decomposition of
bubble deformation (Prosperetti, 1980) to the presence of external forcing. On
the other hand, it is analogous to a liquid-gas interface interpreted as a
collection of harmonic oscillators which are forced by a turbulent background
flow. This is similar to the case of wind wave generation by pressure
fluctuations in a turbulent boundary layer, as proposed by Perrard et al.
(2019) and building on earlier work by Phillips (1957). However, as we will
see, the solutions describing the growth of the bubble interface deformation
lies in a different dynamical regime, in which a large separation of time
scale between stochastic growth and saturation is not fulfilled, in parts
because the turbulent forcing occurs in the dense liquid, imposing rapid
response from the gas phase, contrary to forcing in the turbulent air of much
weaker inertia inducing deformation at the surface of the dense liquid phase.
The utility of the present theoretical framework hence lies mainly on our
ability to test numerically the relation between surface deformation and
turbulent background flow statistics.

We derive an equation linking the growth of the bubble deformation decomposed
in spherical harmonics to the turbulence statistics. Doing so, we identify
different regimes. At short time, for $t\ll t_{2},t_{c}$, where $t_{2}$ is the
capillary time describing the oscillatory response of the bubble, and $t_{c}$
is the eddy turn-over time at the scale of the bubble, the bubble response is
independent of surface tension (Weber number) and the total bubble deformation
$\zeta_{\Omega}$ follows a linear scaling $\zeta_{\Omega}\sim t$ . The
amplitudes of the spherical harmonics coefficients also follow a linear growth
in time, while the response is dominated by the modes 2 of oscillation. At
intermediate time scale $t_{2}<t<t_{c}$, bubble deformation are described by a
sub-linear regime, where the transition time is a function of surface tension
(Weber number). The theory is completed by direct numerical simulations of
bubble deformation and break-up in a turbulent flow, resolving the full two-
phase air-water Navier-Stokes equations. We compute the bubble deformation, as
well as the evolution of the various dynamical modes of deformation. Very good
agreement between the simulations and theoretical predictions are obtained,
and we retrieve the different regimes predicted theoretically. Finally at long
time, the bubble deformation either saturates or the bubble breaks, both
regimes not being described by our linear theory but captured in the
simulations. We eventually link the bubble deformation growth to the
statistical properties of the velocity increments. Doing so, we provide an
estimate of the bubble lifetime statistics, which is parametrized by the
Reynolds number and the Weber number.

The paper is organized as follows, in §2, we present the theory, while in §3
we present the simulations and comparison with the theoretical prediction.
Conclusions and discussions are presented in §4.

##  2 Bubble deformation in a turbulent flow

###  2.1 Setting the scene

We consider an initially spherical gas bubble of radius $R_{0}=d_{0}/2$,
density $\rho_{g}$ and viscosity $\mu_{g}$, immersed in a liquid turbulent
background flow with velocity $\bm{u}(\bm{r},t)$, and of liquid density
$\rho_{\ell}$, viscosity $\mu_{l}$, and surface tension $\gamma$ between the
gas and the liquid phase. The origin of space is set at the centre of mass of
the bubble. Focusing on the early stage of bubble deformation, we limit our
description to bubble surface that can be parametrized in spherical
coordinates $(r,\theta,\varphi)$ by the single-valued distance ${\bf
R}=(R_{0}(t)+\zeta(\theta,\varphi,t)){\bf\hat{r}}$ to the center (Figure 1).
We assume the gas bubble density to be much smaller than the liquid density,
$\rho\equiv\rho_{\ell}\gg\rho_{g}$. In this limit of vanishing bubble inertia,
we neglect the feedback of the bubble interface motion on the velocity and
pressure statistics of the turbulent background flow, valid at low air
fraction in water (Risso, 2018).

We decompose the velocity field $\bm{u}(\bm{r},t)$ into:

${\bf u}={\bf u}_{T}+{\bf u}_{b},$ (2)

where ${\bf u}_{T}$ is the turbulent background flow in the absence of bubble,
and ${\bf u}_{b}$ is the flow associated to the interface motion. The
statistics of the turbulent background flow is not specified a priori. To
characterize fluctuating quantities, we introduce the operation
$\langle\rangle$ for the average over realisations, and the primes (′)
designates the root mean squared value, e.g. $u^{\prime}=\sqrt{\langle
u^{2}\rangle}$.

In the limit of large bubble with respect to Kolmogorov length
$\eta_{K}=(\nu^{3}/\epsilon)^{1/4}$, ($R_{0}\gg\eta_{K}$), the viscous effects
on bubble deformation can be first neglected, the bubble Reynolds number for
low viscosity liquid such as water being relatively high (typically a few
hundred) (Prosperetti, 1980; Magnaudet & Eames, 2000). We then seek ${\bf
u}_{b}$ under the form of a potential flow ${\bf u}_{b}={\bf\nabla}\Phi$. The
velocity field involved being much smaller than the speed of sound in each
phase, we also assume incompressible flow in both phases. The following
calculation is a first order description of the bubble deformation in a
turbulent flow.

###  2.2 Spherical decomposition

Figure 1: Bubble deformation in a turbulent flow : sketch & definitions.

To decompose the velocity field and the surface bubble into a basis adapted to
the spherical geometry, we introduce the spherical harmonic functions
${H}_{\ell,m}(\theta,\varphi)$ :

${H}_{\ell,m}(\theta,\varphi)=\sqrt{\frac{(2\ell+1)(\ell-m)!}{(\ell+m)!}}P^{m}_{\ell}(\cos{\theta})e^{im\varphi}$
(3)

with $\ell\in\mathbf{N}$ and $m\in\left[-\ell,\ell\right]$, $\theta$ and
$\varphi$ are the spherical coordinates and $P^{m}_{\ell}$ the associated
Legendre polynomial. The real form ${Y}_{\ell,m}(\theta,\varphi)$ of the
spherical harmonics writes :

$Y_{\ell}^{m}=\left\\{\begin{array}[]{ll}\sqrt{2}(-1)^{m}\mathbf{Im}(H_{\ell}^{-m})&\mbox{if
}m<0\\\ H_{\ell}^{0}&\mbox{if }m=0\\\
\sqrt{2}(-1)^{m}\mathbf{Re}(H_{\ell}^{m})&\mbox{if }m>0.\\\
\end{array}\right.$ (4)

The real spherical harmonics $Y_{\ell}^{m}$ verifies
$\forall(\ell,m)\in\mathbb{N}\times
[-\ell,\ell],\nabla^{2}r^{-\ell}Y_{\ell}^{m}=0$ for $r>0$. In this basis, the
expressions of bubble shape $R$ and velocity potential $\Phi$ write:

$\displaystyle R(\theta,\varphi,t)$ $\displaystyle=$ $\displaystyle
R_{0}(t)+\sum_{\ell=2}^{+\infty}\sum_{m=-\ell}^{\ell}a_{\ell,m}(t)Y_{\ell}^{m}(\theta,\varphi)$
(5) $\displaystyle\Phi(r,\theta,\varphi,t)$ $\displaystyle=$
$\displaystyle\sum_{\ell=2}^{+\infty}\sum_{m=-\ell}^{\ell}\frac{c_{\ell,m}(t)}{r^{\ell+1}}Y_{\ell}^{m}(\theta,\varphi).$
(6)

where $R_{0}$ is the mean bubble radius and $\Phi$ verifies
$\nabla^{2}\Phi=0$. The coefficients $a_{\ell,m}(t)$ and $c_{\ell,m}(t)$ are
stochastic variables which depend on the background velocity flow
configuration. The deformation will be linked to the radial velocity
component:

$\displaystyle{{\bf u}}\cdot{\bf\hat{r}}=\left(\frac{d{\bf
r_{c}}}{dt}+\bf\tilde{u}_{T}\right)\cdot{\bf\hat{r}}+\partial_{r}\Phi,$ (7)

where $d/dt$ refers to the total derivative with respect to time, $\bf r_{c}$
is the position of the bubble center of mass and $\bf\tilde{u}_{T}$ the
turbulent velocity fluctuations in the bubble frame of reference. We introduce
the spherical harmonic decomposition of $\bf\tilde{u}_{T}$ at the interface in
${\bf r}=\bm{R}$ :

$\displaystyle{{\tilde{\bf u}}_{T}}\cdot{\bf\hat{r}}_{|{\bf
r}=\bm{R}}=\sum_{\ell=2}^{+\infty}\sum_{m=-\ell}^{\ell}u_{\ell,m}(t)Y_{\ell}^{m}(\theta,\varphi),$
(8)

The mode zero coefficient $u_{0,0}$ vanishes from incompressibility. The
statistics of $u_{\ell,m}$ are a priori unknown : an evaluation of their
statistics will be provided using direct numerical simulations in section 3 in
the absence of bubble.

Using the continuity of normal stress in $\bm{r}=\bm{R}$ in the limit of
vanishing viscosity, we also decompose the pressure field $P_{\ell}$ in the
liquid phase at the interface (${\bf r}=\bm{R}$) in the spherical harmonics
base under the form:

$\displaystyle P_{\ell}(\theta,\varphi,t)$ $\displaystyle=$ $\displaystyle
P_{0}(t)-\gamma\kappa(\theta,\varphi,t)+\rho_{\ell}\sum_{\ell=2}^{+\infty}\sum_{m=-\ell}^{\ell}\beta_{\ell,m}\left(\frac{R_{0}}{r}\right)^{\ell+1}Y_{\ell}^{m}(\theta,\varphi),$
(9)

where $P_{0}(t)$ is the pressure inside the bubble, $\kappa$ is the sum of the
two principal curvatures, and the coefficients ${\beta}_{\ell,m}$ encode the
turbulent pressure fluctuations associated with the turbulent background flow.
Their statistics will be similarly estimated in the DNS section 3 in the
absence of bubble.

###  2.3 Forced oscillator solution

In this section, we look for an equation for the coefficient $a_{\ell,m}$ at
the linear order in $a_{\ell,m},u_{\ell,m}$ and ${\beta}_{\ell,m}$, which
holds for small deformations $\zeta\ll R_{0}$. The bubble surface
$R=R_{0}+\zeta$ and velocity field ${\bf u}$ are related throughout the
kinematic boundary condition

$\frac{d\zeta}{dt}={\tilde{{\bf u}}}({\bf R},t)\cdot{\bf\hat{r}}.$ (10)

For small deformation ($\zeta/R_{0}\ll 1$), the convective terms
$u_{\theta}/r\partial_{\theta}\zeta$ and
$u_{\varphi}/r\partial_{\varphi}\zeta$ are negligible compared to $u$ and
$d\zeta/dt$ is approximated by $\partial\zeta/\partial t$. In the limit of
incompressible flow, mass conservation imposes:

$\frac{d}{dt}\iint\textrm{d}\Omega\left(R_{0}(t)+\zeta(\theta,\varphi,t)\right)^{3}=0,$
(11)

where $d\Omega=\sin(\theta)d\theta d\varphi$ is the elementary solid angle. In
the limit of small deformation $\zeta\ll R_{0}$, we have $\dot{R}_{0}=0$.
Using the harmonic decomposition (Eq. 5, 6) in Eq. 10 and incompressible
limit, we get at the first order in $a_{\ell,m}/R_{0}$,

$\displaystyle\dot{a}_{\ell,m}=\frac{\ell}{R_{0}^{\ell+2}}c_{\ell,m}+u_{\ell,m},$
(12)

where overdots ($\dot{}$) corresponds to differentiation with respect to time.
Combining Eq. 10 with the momentum equation
$\partial_{t}u_{r}+1/2\partial_{r}(u_{r}^{2})=-1/\rho\partial_{r}p$ integrated
from $+\infty$ to $r=R(\theta,\varphi,t)$, we get :

$\displaystyle\left[\partial_{t}\Phi+\frac{1}{2}(\partial_{r}\Phi)^{2}\right]_{+\infty}^{R_{0}}=\frac{1}{\rho}(P_{\infty}-P_{\ell}),$
(13)

where Using the harmonics decomposition of $\Phi$ and the pressure $P$
combined with Eqs. 13, 12 yields an equation at the first order in
$a_{\ell,m}$ :

$\displaystyle\ddot{a}_{\ell,m}+(\ell-1)(\ell+1)(\ell+2)\omega_{0}^{2}a_{\ell,m}=\dot{u}_{\ell,m}+\ell\frac{1}{R_{0}}\beta_{\ell,m},$
(14)

where the curvature $\kappa$ has been obtained at the first order in
$a_{\ell,m}$ using $\kappa\ell=(\ell-1)(\ell+1)(\ell+2)a_{\ell,m}/R_{0}^{2}$
111https://ciks.cbt.nist.gov/garbocz/paper134/node4.html. We have introduced
$\omega_{0}^{2}=\gamma/(R_{0}^{3}\rho)$ the typical frequency associated with
bubble oscillation and
$\omega_{\ell}=\omega_{0}\sqrt{(\ell-1)(\ell+1)(\ell+2)}$, with in particular
$\omega_{2}=2\sqrt{3}\omega_{0}$. We also introduce the characteristic time
$t_{\ell}=1/\omega_{\ell}$ of mode $\ell$.

Eq. 14 provides a description of the bubble deformation due to turbulent
stresses under the set of hypothesis we have made, (i) absence of feedback of
the surface deformation $\zeta$ on the velocity and pressure statistics
$u_{\ell,m}$ and ${\beta}_{\ell,m}$, (ii) small deformation $\zeta/R_{0}\ll
1$, (iii) negligible viscous dissipation, i.e.
$\nu_{\ell}/R_{0}^{2}\ll\omega_{2}$. This assumption of negligible viscous
dissipation is coherent with Risso & Fabre (1998) discussion on the viscous
relaxation time scale $t_{\nu}=d_{0}^{2}/(80\nu)$. Eq. 14 corresponds to an
harmonic oscillator with a source term given by the pressure and velocity
fluctuations at the interface. The l.h.s is identical to Prosperetti (1980),
while the r.h.s derives from the forcing by the velocity acceleration and the
background pressure fluctuations.

Solving Eq. 14 requires initial conditions on $a_{\ell,m}$ and
$\dot{a}_{\ell,m}$, which are set by the initial bubble shape and the
kinematic boundary condition at $t=0$. Starting from an initially spherical
bubble at rest, we have :

$\displaystyle R(\theta,\varphi,0)=R_{0};\text{ and
}\dot{R}(\theta,\varphi,0)=u_{r}(R_{0},\theta,\varphi,0),$ (15)

which gives the initial conditions for $a_{\ell,m}$:

$\displaystyle a_{\ell,m}(0)=0;\text{ and }\dot{a}_{\ell,m}(0)$
$\displaystyle=$ $\displaystyle u_{\ell,m},$ (16)

###  2.4 Link with surface deformation

We seek for an estimate of the average global surface deformation. We
introduce the root mean squared $\zeta_{\Omega}$ of the bubble interface from
an integration over all solid angles:

$\zeta_{\Omega}^{2}=\frac{1}{4\pi}\iint\textrm{d}\Omega\zeta^{2}.$ (17)

According to the orthogonality of spherical harmonics
$\frac{1}{4\pi}\iint\textrm{d}\Omega
Y_{\ell}^{m}Y_{\ell^{\prime}}^{m^{\prime}}=\delta_{\ell,\ell^{\prime}}\delta_{m,m^{\prime}}$,
$\zeta_{\Omega}$ can be expressed as

$\zeta_{\Omega}^{2}(t)=\sum_{\ell=2}^{+\infty}\sum_{m=-\ell}^{\ell}a_{\ell,m}(t)^{2}.$
(18)

The surface deformation is hence directly related to the solutions of Eq. 14.
$a_{\ell,m}$ being stochastic variables, we look for its temporal behaviour as
well as temporal behavior of its ensemble averaged value $\langle
a_{\ell,m}^{2}\rangle$.

###  2.5 Controlling time scales

Figure 2: Ratio of characteristic bubble oscillation time scales with the
turbulence time scale of the problem as a function of the Weber number. The
turbulence time scale is the eddy turn over time at the bubble scale,
$t_{c}=\epsilon^{-1/3}d_{0}^{2/3}$, the capillary oscillation time is
$t_{0}=1/\omega_{0}=1/\sqrt{\gamma/(R_{0}^{3}\rho)}$ and the time associated
to mode 2 is $t_{2}=1/\omega_{2}=1/(2\sqrt{3}\omega_{0})$. Symbols correspond
to the DNS cases in section 3.

The solutions of Eq. 14 depends on the various time scales that are inherently
associated to an noise-excited oscillator. For a bubble evolving in a
turbulent flow, previous studies have clearly identified one controlling time
scale, the natural bubble oscillation response occurring at a capillary time
scale (Prosperetti, 1980; Risso & Fabre, 1998; Ravelet et al., 2011). Another
relevant time scale lies in the turbulence correlation time of the
fluctuations at the bubble scale also discussed by several authors (Risso &
Fabre, 1998; Martinez-Bazan et al., 1999; Ravelet et al., 2011). Depending on
the Weber number, a final time scale is eventually given either by the bubble
lifetime or a saturation time.

The modes $\ell=2$ carry most of the surface energy as they correspond to the
largest available spatial scale. Therefore, the typical characteristic
amplitude and time scales of $\zeta_{\Omega}$ are set by the modes $\ell=2$ of
oscillation. The first relevant time scale is the oscillator reduced periods
$t_{2}=1/\omega_{2}$ of the second mode $\ell=2$ of oscillation, which has
been proven to dominate the deformation (Risso & Fabre, 1998). On a time
significantly shorter than the reduced period $\omega_{2}t<1$, the response
will be purely inertial and independent of surface tension (Risso & Fabre,
1998).

The second characteristic time scale arises from the turbulence temporal
correlations ${\tau}_{\ell,m}$ and ${\tilde{\tau}}_{\ell,m}$ of respectively
the coefficients $u_{\ell,m}$ and ${\beta}_{\ell,m}$ of the spherical
decomposition of the velocity and pressure fluctuations. They can be defined
from the autocorrelation functions of velocity and pressure respectively :

$\displaystyle{\tau}_{\ell,m}$ $\displaystyle=$
$\displaystyle\frac{1}{u_{\ell,m}^{2}}\int_{-\infty}^{+\infty}\textrm{d}s~{}\langle
u_{\ell,m}(t)\cdot u_{\ell,m}(t+s)\rangle.$ (19)
$\displaystyle{\tilde{\tau}}_{\ell,m}$ $\displaystyle=$
$\displaystyle\frac{1}{{\beta}_{\ell,m}^{2}}\int_{-\infty}^{+\infty}\textrm{d}s~{}\langle{\beta}_{\ell,m}(t)\cdot{\beta}_{\ell,m}(t+s)\rangle.$
(20)

In the absence of consensus theory on temporal correlation functions in
turbulence (Pope, 2000), we estimate the correlation times of the dominant
modes $\tau_{2,m}$ and $\tilde{\tau}_{2,m}$ from the eddy turn over time at
the bubble scale $t_{c}=\epsilon^{-1/3}d_{0}^{2/3}$.The ratio
$\tilde{\tau}_{2,m}/t_{c}$, however, may depend on the detailed statistics of
the turbulent fluctuations. The correlation times $\tau_{\ell,m}$ and
$\tilde{\tau}_{\ell,m}$ of the higher order modes ($\ell>2$) are expected to
be smaller, as they correspond to smaller scales. In the following, we neglect
the temporal decoherence of higher order modes ($\ell>2$) as we expect the
main deformation to originate from the modes $\ell=2$ of deformation. Note
that for higher modes, $\ell>2$, the capillary time 1$/\omega_{\ell}$ is
expected to be shorter than the corresponding high wavenumber turbulent time
$\tau_{\ell}$.

Figure 2 shows the ratio of the typical capillary time scales ($t_{0}$ and
$t_{2}$) with the eddy turn-over time scale $t_{c}$ as a function of Weber
number. Below the critical value $\textrm{We}<\textrm{We}_{c}$ with
$\textrm{We}_{c}=3$, we have $t_{2}\ll t_{c}$, and this separation persists at
moderate We while being reduced. The crossover between these two time scales,
$t_{2}=t_{c}$ occurs at higher We number, exactly $\textrm{We}=192$ (and
$\omega_{2}=\sqrt{192}/(t_{c}\sqrt{We})$). For reference, we provide typical
time scales in table 1, corresponding to experimental conditions of bubbles
evolving in turbulent flow created from underwater turbulent jets (Ruth et
al., 2019; Risso & Fabre, 1998; Martinez-Bazan et al., 1999; Vejražka et al.,
2018). Note that these experimental systems can present large variations in
turbulence levels within the flow, and significant mean shear. A complete
discussion on the challenges of experimental systems to study bubble
deformation and break-up in homogeneous and isotropic turbulence can be found
in Masuk et al. (2019). The long time evolution of bubble deformation depends
eventually on the bubble stability. The stability being dependent on the
deformation, we postpone the discussion on either the saturation or the break-
up time to the end of the direct numerical simulations section §3. The
following sections derive asymptotic solutions associated to each temporal
regimes.

Reference |  $\epsilon$ (m2s-3) |  $L$ (mm) |  $d$ (mm) |  $t_{2}$ (ms) |  $t_{c}$ (ms)
---|---|---|---|---|---
Ruth et al. (2019) |  0.5-1 |  13 |  2-6 |  1-6 |  20-40
Risso & Fabre (1998) |  0.1-2 |  21 |  4-16 |  2-40 |  50-140
Ravelet et al. (2011) |  1 |  10 |  9 |  10 |  50
Martinez-Bazan et al. (1999) |  3-100 |  48 - 115 |  3 |  2 |  5-35
Vejražka et al. (2018) |  0.1-10 |  4-6 |  2-6 |  1-6 |  15-140
Masuk et al. (2019) |  0.16 |  45-60 |  2-6 |  1-6 |  15-30
Table 1: Estimated experimental time scales and length scales for typical
laboratory experiments of bubble deformation and break-up in turbulence. All
experiments are made with air bubbles in water. Dissipation rate $\epsilon$
and bubble diameter $d$ were directly provided by all authors. An estimate of
the integral length scale $L$ were provided by most of the authors, we infer
the values for Martinez-Bazan et al. (1999) using the usual statistical
properties of a turbulent jet (Hussein et al., 1994). The oscillation time is
given by $t_{2}=\frac{1}{2\sqrt{3}\sqrt{\gamma/(R_{0}^{3}\rho)}}$ and the eddy
turn over time at the bubble scale is given by
$t_{c}=\epsilon^{-1/3}d_{0}^{2/3}$.

###  2.6 Early time solution: linear regime

Starting from a spherical shape at $t=0$, the first instants verify
$\omega_{\ell}t\ll 1$ and $t\ll t_{c}$. In this limit, we can consider
$u_{\ell,m}$ and ${\beta}_{\ell,m}$ as independent of time, and neglect the
surface tension term. From eq. 14 and its initial conditions, we get by double
integration over time:

$\displaystyle a_{\ell,m}=u_{\ell,m}t+\frac{1}{2R_{0}}{\beta}_{\ell,m}t^{2}.$
(21)

The relative influence of each term will depend on the characteristic time
scale $2R_{0}u_{\ell,m}/{\beta}_{\ell,m}$, of order $t_{c}$ as it also
corresponds to a characteristic time scale of the turbulence at the bubble
scale $R_{0}$. For $t\ll t_{c}$, the pressure contribution is then negligible
and the dynamics of each mode $a_{\ell,m}$ reduces to a simple advection from
the background flow. Performing ensemble average operation yields:

$\displaystyle\langle a_{\ell,m}^{2}\rangle$ $\displaystyle=$
$\displaystyle\langle u_{\ell,m}^{2}\rangle t^{2},$ (22)

At short time ($\omega_{2}t\ll 1$, $t\ll t_{c}$), the evolution of the
$a_{\ell,m}$ is then expected to be independent of surface tension and
pressure fluctuations, and this first linear growth regime should hold for all
Weber numbers. Note that this limit is similar to the frozen turbulence regime
identified for bubble pinch-off in turbulence, in which most of the pinching
dynamics occur on a fraction of the eddy turn over time at the scale of the
neck (Ruth et al., 2019). In both cases, the influence of the turbulent
background flow reduces to the set of the initial conditions.

Later on, the linear growth regime is limited in time by either the bubble
natural oscillations at $1/\omega_{2}$ or by the correlation time of the
turbulent flow $t_{c}$, depending on the Weber number as shown in fig. 2. For
$\textrm{We}\gg 1$, when $t_{c}<t_{2}$, linear regime occurs until $t\approx
R_{0}/u_{\ell,m}^{\prime}\approx t_{c}$, which corresponds to
$\zeta_{\Omega}/R_{0}\approx 1$, a deformation amplitude sufficient to trigger
bubble breaking. $t_{c}$ then coincides with the bubble lifetime, and no later
time evolution is expected.

###  2.7 Forced oscillator regime

We consider $\textrm{We}<100$, in which the linear growth is first limited by
surface tension. For $t<t_{c}$, we can still consider $u_{\ell,m}$ and
${\beta}_{\ell,m}$ as independent of time. We have :

$\displaystyle\ddot{a}_{\ell,m}+\omega_{\ell}^{2}a_{\ell,m}$ $\displaystyle=$
$\displaystyle\ell\frac{1}{R_{0}}\beta_{\ell,m}$ (23)
$\displaystyle\dot{a}_{\ell,m}(0)$ $\displaystyle=$ $\displaystyle u_{\ell,m}$
(24)

Which admits a solution of the form :

$\displaystyle
a_{\ell,m}=\frac{u_{\ell,m}}{\omega_{\ell}}\sin(\omega_{\ell}t)+\frac{\ell{\beta}_{\ell,m}}{R_{0}\omega_{\ell}^{2}}(1-\cos(\omega_{\ell}t))$
(25)

The evolution of $\langle a_{\ell,m}^{2}\rangle$ at short time is given by :

$\displaystyle\langle a_{\ell,m}^{2}\rangle=\langle u_{\ell,m}^{2}\rangle
t^{2}+\frac{\ell\langle{\beta}_{\ell,m}u_{\ell,m}\rangle}{R_{0}}t^{3}-\frac{1}{3}\omega_{\ell}^{2}\langle
u_{\ell,m}^{2}\rangle
t^{4}+\frac{\ell^{2}\langle{\beta}_{\ell,m}^{2}\rangle}{R_{0}^{2}}\frac{t^{4}}{4}+O(t^{5}),$
(26)

However, the terms in $\langle{\beta}_{\ell,m}u_{\ell,m}\rangle/R_{0}$ and
$\langle{\beta}_{\ell,m}^{2}\rangle/R_{0}^{2}$ are respectively of order
$\langle u_{\ell,m}^{2}\rangle/t_{c}$ and $\langle
u_{\ell,m}^{2}\rangle/t_{c}^{2}$, their contribution can be neglected for
$t\ll t_{c}$. The first relevant non linear term in the range
$1/\omega_{\ell}<t<t_{c}$ is then the term $-1/3\langle u_{\ell,m}^{2}\rangle
t^{4}$, and the equation for $a_{\ell,m}$ becomes :

$\displaystyle\langle a_{\ell,m}^{2}\rangle=\langle u_{\ell,m}^{2}\rangle
t^{2}-\frac{1}{3}\omega_{\ell}^{2}\langle u_{\ell,m}^{2}\rangle
t^{4}+O(t^{5}),$ (27)

which corresponds to a saturation of the linear growth by capillary forces on
a time scale $1/\omega_{\ell}$.

###  2.8 Analogy with the Phillips regime

In the limit of small Weber number $\textrm{We}<\textrm{We}_{c}$, there will
be a later time evolution for $t>t_{c}$ at which the pressure and velocity
fluctuations decorrelate. Eq. 14 becomes a Langevin equation with a forcing
term characterized by its correlation time $t_{c}$ and its root mean squared
amplitude $A$ :

$\displaystyle\langle
A^{2}\rangle=\sum_{\ell=2}^{+\infty}\sum_{m=-\ell}^{\ell}\langle{\beta}_{\ell,m}^{2}\rangle+\langle\dot{u}_{\ell,m}^{2}\rangle.$
(28)

For $t\gg t_{c}$ we have by analogy with a Langevin equation :

$\displaystyle{\zeta_{\Omega}}^{2}=R_{0}^{2}\langle
a_{2}^{2}\rangle=\frac{4\langle
A^{2}\rangle}{R_{0}^{2}}\frac{tt_{c}}{\omega_{2}^{2}},$ (29)

which is the analogous for bubble deformation of the Phillips regime for wind-
wave generation (Phillips, 1957; Perrard et al., 2019):

$\displaystyle\langle\zeta_{\Omega}^{2}\rangle=\frac{\langle p^{2}\rangle
t}{2\sqrt{2}\rho_{\ell}^{2}gU_{c}}$ (30)

with $U_{c}=1/(kt_{c})$ and $\omega^{2}=gk$. In practice, this regime may not
be observed, as the saturation by non linear effects should occur on a time
scale comparable with $t_{c}$. The main difference between the Phillips regime
and the current case is the ratio of density between the pressure fluctuations
and the liquid inertia. In the Phillips case, the forcing term is in the air
phase, whereas the inertia is associated to the density of the liquid phase
$\rho_{\ell}$. The energy transfer per oscillation period is then much weaker
than the wave energy, and a slow accumulation of energy per period can be
observed. In the bubble case, the forcing being located in the liquid phase,
the pressure fluctuations and the wave energy are of the same order of
magnitude, and an oscillation of amplitude unity is usually reached in only
few periods of oscillation.

##  3 Direct Numerical simulation of bubble deformation in a turbulent flow

###  3.1 Numerical methods: the Basilisk flow solver

We perform direct numerical simulations of the three-dimensional,
incompressible Navier-Stokes equations, either with a single phase (the
turbulence precursor simulation) or with two-phases (air bubble and turbulent
water) with surface tension, using the free software Basilisk
222http://basilisk.fr/ (Popinet, 2009, 2018). We use a spatial adaptive octree
grid allowing to save computational time while resolving the different length
scales of the problem and a momentum conserving scheme and the interface is
reconstructed by a sharp geometric Volume of Fluid (VOF) method (Popinet,
2009, 2018). The solver has been extensively described in recent publications
(Popinet, 2015, 2018; van Hooft et al., 2018; Fuster & Popinet, 2018; Mostert
& Deike, 2020), and its accuracy has been largely validated on complex
multiphase flow, including bubble dynamics (Cano-Lozano et al., 2016), bubble
bursting (Deike et al., 2018; Lai et al., 2018; Berny et al., 2020), and wave
breaking (Deike et al., 2016; Mostert & Deike, 2020; Mostert et al., 2020). We
do not consider the effect of gravity in this work. The turbulent two-phase
simulations of bubble deformation in turbulence is presented below.

###  3.2 Preparation: creation of the turbulence and insertion of the bubble

The turbulent flow is generated by adding to the Navier-Stokes equation a
volumetric forcing term locally proportional to the velocity field $f{\bf u}$,
following the approach of Rosales & Meneveau (2005), previously implemented
and provided as an example in the Basilisk
library333http://basilisk.fr/src/examples/isotropic.c. Rosales & Meneveau
(2005) showed that such forcing in every point of the real space leads to a
well characterized homogeneous and isotropic turbulent flow with properties
similar to those obtained with a spectral code and forcing. Such approach has
also been used by Loisy & Naso (2017) to study rising bubbles in a turbulent
flow.

We consider a cubic box of size L with periodic boundary conditions on each
side. Adaptive mesh refinement is used on the velocity field, and the maximum
level of refinement $N$ allows comparison with fixed grid resolution having
$2^{N}$ grid points in each direction. The turbulent flow is generated for
increasing resolutions with $N$ going from 6 to 8, corresponding to equivalent
$64^{3}$ to $256^{3}$ grid size on a fixed grid. This resolution is modest but
will be increased around the interface when we inject the bubble in the flow,
and is chosen to keep the computational cost reasonable for each simulation as
we aim to perform a large number of simulations to access statistical
properties of the bubble deformation dynamics.

Figure 3a shows the time evolution of the kinetic energy in the precursor
simulation,
$K=\frac{1}{V}\iiint\frac{1}{2}\rho_{\ell}u(\vec{x},t)^{2}\mathrm{d}V$. After
a short transient, injection and dissipation of energy eventually balances on
average and we obtain a statistically stationary, homogeneous and isotropic
turbulence. The statistically stationary state is reached after approximately
15 eddy turn over time at the integral scale $\tau=u^{\prime 2}/\epsilon$,
where $\epsilon$ is the asymptotic dissipation rate and $u^{\prime}$ the
asymptotic root mean squared velocity. The associated Taylor Reynolds number
is $Re_{{\lambda}}=\frac{2K}{3\mu}\sqrt{\frac{15\mu}{\epsilon}}\approx 38$
which is a typical value for current two-phase simulations of turbulent flow
(Loisy & Naso, 2017; Elghobashi, 2019). Convergence of the statistical
properties of the flow with respect to numerical resolution is achieved and is
shown with simulations performed at lower maximum refinement. Note that while
the flow is statistically stationary and equivalent for the various
resolutions, because of the chaotic nature of the flow and the forcing each
realization presents differences one from each other.

Figure 3: (a) Kinetic energy as a function of time in the precursor
simulation, for increasing numerical resolution, respectively level N=6, 7 and
8. A statistically stationary state is reached after 15$\tau$ and numerical
convergence is reached for level N=7, with a corresponding turbulent Reynolds
number Reλ=38. A bubble will be immersed in this statistically stationary
flow, using different initial times as initial conditions. (b) Second order
structure function $D_{LL}$ and $D_{NN}$ in the longitudinal and transverse
directions respectively, compensated by the homogeneous and isotropic
turbulence scaling $(r\epsilon)^{-2/3}$ and $D_{LL}=3/4D_{NN}$. Kolmogorov
theory is superimposed in red dashed line. (c) Visualization of a snapshot of
the instantaneous radial velocity component of the turbulent flow evaluated on
a sphere of radius $R/\eta_{k}=32$ that a bubble of corresponding size will
face.

Figure 3b shows the statistical properties of the turbulent flow once the
stationary state is reached. We characterize the fluctuations using the second
order structure functions in the longitudinal $D_{LL}(d)$ and in the
transverse direction $D_{NN}(d)$, defined as :

$\displaystyle D_{LL}(d)$ $\displaystyle=$
$\displaystyle\frac{1}{3}\sum_{i}\langle\left(u_{i}(\bm{r},t)-u_{i}(\bm{r}+d{\bf\hat{r}_{i}},t)\right)^{2}\rangle$
(31) $\displaystyle D_{NN}(d)$ $\displaystyle=$
$\displaystyle\frac{1}{6}\sum_{i\neq
j}\langle\left(u_{i}(\bm{r},t)-u_{i}(\bm{r}+d{\bf\hat{r}_{j}},t)\right)^{2}\rangle,$
(32)

for homogeneous and isotropic flows, with ${\bf\hat{r}_{i}}$ the unit vector
along the $i$ direction. The transverse structure function is compensated by
its scaling for an homogeneous and isotropic flow $(d\epsilon)^{2/3}$, and we
indeed observe a plateau value close to $C=2$ (Pope, 2000). The relation
$D_{LL}=3/4D_{NN}$ is also verified by representing the compensated
longitudinal structure function $4/3D_{LL}(d)(d\epsilon)^{-2/3}$. The inertial
range is obviously quite limited due to the relatively coarse resolution, but
the turbulent flow at the scale of the bubble to be injected is reasonable,
and the bubble radius lies within the inertial range. The quantity of interest
for the bubble deformation is the spherical velocity increments, which are
defined on a sphere of radius $R$ by :

$\displaystyle\delta u_{S}(R,\theta,\varphi)$ $\displaystyle=$
$\displaystyle{\tilde{{\bf u}}}({\bf R},t)\cdot{\bf\hat{r}}$ (33)
$\displaystyle\tilde{{\bf u}}({\bf R})$ $\displaystyle=$ $\displaystyle{\bf
u}_{T}-\frac{1}{4\pi}\iint\textrm{d}\Omega{\bf u}_{T}({\bf
R},\theta,\varphi,t).$ (34)

From homogeneity and isotropy, the statistical properties of $\delta u_{S}(R)$
shall only depend on $r$. Figure 3c shows the velocity increment $\delta
u_{S}$ over a sphere of radius $R/\eta_{k}=32$, displaying an example of the
broad range of forcing scales the bubble sees and feels in the flow. The
ensemble average value $\delta u^{\prime}_{S}(R)=\sqrt{\langle(\delta
u_{S}(R))^{2}\rangle}$ compensated by $(\epsilon d)^{2/3}$ is represented in
red in figure 3b. We observe $\delta u^{\prime}_{S}(R)=D_{NN}(d/2)$ for
$Re_{\lambda}=38$ in the entire inertial range. The statistical properties of
the spherical increment and their link with surface deformations will be
further discussed in section 3.

Once the statistically stationary regime is reached, the temporal recording of
the velocity field is stored from $20$ to $60\tau$. Different instants are
used as initial times for numerical simulations of bubble deformation and
break-up. A central sphere of radius $R_{0}$, diameter $d_{0}$ and density
$\rho_{g}=1/850\rho_{\ell}$ is placed in the periodic box, and the flow in the
inner phase is initially set to zero. The bubble diameter is located in the
inertial range. For $Re_{\lambda}=38$ we have $d_{0}/\eta_{K}=17.6$,
$d_{0}/\lambda=1.49$ and $d_{0}/L=0.13$ where $\lambda$ is the Taylor
microscale $\lambda=\sqrt{15\nu u^{\prime 2}/\epsilon}$ for homogeneous and
isotropic turbulence and $L$ is the box size. The level of refinement around
the interface is $N=9$ for the majority of the numerical runs, and convergence
tests at level $N=10$ have been performed. The results on bubble deformation
by the turbulent flow presented here are independent of the resolution between
refinement levels $N=9$ and $N=10$ as shown in the Appendix. These resolutions
correspond to 70 and 140 points across the initial bubble diameter, which is
comparable to the resolution successfully used to resolve rising motion and
path instability of bubbles using the same numerical methods (Cano-Lozano et
al., 2016).

In a turbulent flow, the critical Weber number, or break-up threshold, is
defined in a statistical sense, therefore the probability of breaking does not
vanish for Weber number immediately smaller than
$\textrm{We}=\textrm{We}_{c}$. To determine the critical Weber number, we
perform ensemble of simulations by using different initial times from the
turbulence precursor simulation (each precursor typically spaced by
$1-3t_{c}$). We do not observe any breakup at $\textrm{We}=1.5$ while running
the simulations up to $20t_{c}$ for more than 10 runs with different initial
conditions. For $\textrm{We}=3$, we observe a little over 50% of break-ups
after $20t_{c}$ for an ensemble of 35 runs with different initial conditions.
It is clear that given the broad lifetime statistics of a bubble in these
conditions close to stability, the percentage of cases that do not break might
change if the simulations are run for longer times. For $\textrm{We}=6$, and
above, all bubbles break within a few $t_{c}$. As such we consider
$\textrm{We}_{c}=3$ as the critical Weber number in our configuration and this
value is within the variation of experimental measurements discussed in the
literature.

Figure 4: Snapshots of a bubble for one run at $\textrm{We}=3$, with interface
in white and each background plane showing one component of the velocity. The
bubble is injected in the center of the turbulent domain at $t=0$ with no
velocity and starts to deform quickly. During the first eddy turnover time, no
break-up event occurs, but we observe strong and erratic deformations over
time as well as advection by the turbulent flow.

The Weber number, eq. 1, is defined considering the temporal average of the
turbulence dissipation rate, and we vary the Weber number by changing the
surface tension. Starting from a sphere, we study the growth of the surface
deformation until either the saturation of bubble deformation for stable
bubbles or the breaking for unstable cases. We perform simulations for a wide
range of Weber number, $1.5<We<45$. The critical Weber number is found at
$\textrm{We}_{c}=3$, for which about half of the runs do not exhibit bubble
break-up after 20$\tau$. We will discuss simulations where for small Weber
numbers ($\textrm{We}\leq 1.5$) no break-up occurs, while for large Weber
number ($\textrm{We}\geq 3$) we will analyze the dynamics before the first
break-up. We perform an ensemble of simulations of 10 to 15 runs for each We
number, using various stored initial conditions from the precursor flow field,
leading to a total of about 80 simulations.

Figure 4 shows an example of a bubble evolving in the turbulent flow field in
the vicinity of the instability threshold at We=3, during the first eddy turn
over time at the bubble scale $t_{c}$. The bubble is initially spherical and
quickly deforms for $t<0.5t_{c}$, and then starts to exhibit erratic
oscillations. These deformations do not lead to break-up, and occur together
with advection by the turbulent background flow.

###  3.3 Computing the bubble deformation through a Voronoi decomposition

The interface sampling points $\bf r_{i}$ are by essence non uniformly spread,
as more points are dynamically added in regions of larger curvature due to the
adaptive algorithm. To compute averages accurately over the bubble interface,
we then construct a surface scheme using a spherical Voronoi decomposition
(Voronoi, 1908). The bubble surface is partitioned into regions close to each
of the points. For each region located around the point ${\bf r}_{i}$, we
associate a region area $A_{i}$ and a corresponding solid angle
$\Omega_{i}=A_{i}/r_{i}^{2}$. We evaluate the surface integral of any function
$f$ by :

$\displaystyle\frac{1}{4\pi}\iint\textrm{d}\Omega f({\bf
r})=\frac{1}{4\pi}\sum_{i}f({{\bf r}}_{i})\Omega_{i},$ (35)

where $\sum_{i}\Omega_{i}=4\pi$. The computation of the region areas $A_{i}$
and the corresponding solid angles $\Omega_{i}$ is illustrated in figure 5.
The $\Omega_{i}$ are computed in three steps. We first project all interface
points on the unit sphere (step 1). We then compute the region locations using
a spherical Voronoi algorithm based on a robust Delaunay triangulation (Caroli
et al., 2009) (step 2). The area $A_{i}$ of each convex polygon is eventually
estimated using the shoelace formula (step 3). On the unit sphere, the $A_{i}$
and the $\Omega_{i}$ are identical. An example of our Voronoi decomposition is
shown in Figure 5b on the unit sphere, and in Figure 5c after projection on
the initial bubble shape. The numerical error on the total solid angle
$\sum_{i}\Omega_{i}$ is less than 0.1$\%$.

Figure 5: Illustration of the interface tesselation using spherical Voronoi
diagram. a) Bubble deformed by the turbulent background (one case at
$\textrm{We}=3,t/t_{c}=0.22$). The output data points are unevenly distributed
at the interface. The bubble surface in grey levels has been reconstructed
using a two dimensional spline interpolation. b) Projection of each point on
the unit sphere, and computation of a spherical Voronoi diagram. Each polygon
corresponds to the region associated to a single data point. c) Projection of
the Voronoi diagram onto the initial shape bubble using the spline
interpolation. The weight of each datapoint is given by the associated polygon
area.

To analyse the interface dynamics using Eq. 14, we introduce the spherical
harmonics $Y_{\ell}^{m}(\theta,\varphi)$. The numerical estimate of the
harmonic coefficient $a_{\ell,m}$ is obtained by a sum over the interface
points ${\bf r}_{i}$ weighted by their local solid angle $\Omega_{i}$ :

$\displaystyle a_{\ell,m}=\sum_{i}\Omega_{i}\zeta({\bf
r}_{i})Y_{\ell}^{m}(\theta_{i},\varphi_{i})$ (36)

where $(r_{i},\theta_{i},\varphi_{i})$ are the spherical coordinates of ${\bf
r}_{i}$. We have first checked the orthonormal property of the spherical
harmonics $\frac{1}{4\pi}\iint\textrm{d}\Omega Y_{\ell}^{m}\cdot
Y_{(\ell^{\prime},m^{\prime})}=\delta_{\ell,\ell^{\prime}}\delta_{m,m^{\prime}}$
for $\ell\in[0,10]$ and all corresponding $m$ values on a set of 900 points,
typical of a surface bubble sampling, and randomly spread on a unit sphere.
The typical error is 0.1% for $\ell<5$, and increases for higher $\ell$
values. The harmonic decomposition has been tested on synthetic shapes, i.e. a
set of 900 points randomly spread on an interface of known spherical
decomposition. The relative error on the harmonic coefficient is about $0.5\%$
of the largest non-zero harmonic coefficient for $\ell\in[1,5]$. In practice,
the amplitude decreasing rapidly with $\ell$, an accurate estimate of the
harmonics coefficient is limited to $\ell\leq 5$ for the interface deformation
and around $\ell\sim 10$ for the velocity field.

The bubble center position evolves with time, the interface being advected by
the turbulent background flow. To analyse specifically the surface
deformations, we look for the centre position ${\bf r}_{c}$ of the bubble
frame of reference, for which $a_{1,m}=0$ for $m\in\\{-1,0,1\\}$, as defined
in section 2. The centre position ${\bf r}_{c}$ is computed recursively as
follow. At each step, we compute the Voronoi diagram and the associated
spherical harmonic functions $Y_{1,m}$ using the regions area $A_{i}$. Then,
since each function $Y_{1,m}$ presents a symmetry of revolution, a one
dimensional gradient descent on each function is sufficient to find the
position center that minimizes each coefficient $a_{1,m}$. For bubble
interfaces with single valued radial distance $r_{i}(\theta,\varphi)$, the
gradient descent indeed converges to a centre position ${\bf r}_{c}$ for which
the three mode 1 coefficients $a_{1,m}$ vanish.

###  3.4 Deformation dynamics for stable conditions: bubble deformation and
temporal evolution of spherical modes

Figure 6: (a) Total deformation $\zeta_{\Omega}/R_{0}$ as a function of time,
together with the rms amplitude of the mode 2 oscillations. (b) Energy
$a_{\ell}^{2}=\sum_{m=-l}^{m=l}a_{l,m}^{2}$ of the spherical harmonics for
modes $\ell=1$ to $\ell=5$ for one run at $\textrm{We}=3$. Mode 1 energy
vanishes by definition. Mode 2 appears as the most energetic mode, in
agreement with experimental literature, with oscillation close to
$f_{2}=\omega_{2}/(2\pi)$. As $\ell$ increases, the energy contained decreases
sharply. (c) Temporal evolution of the amplitude coefficients $a_{2,m}$ for
mode $\ell=2$, and $m=-2$ to $m=+2$. All modes $m$ oscillate with similar
amplitude and frequency.

We first describe the temporal evolution of the coefficients $a_{\ell,m}$ for
a single run at We=3, as an illustration of the analysis performed on each
simulation. Several modes $a_{\ell,m}$ of same $\ell$ and different $m$ values
are associated to the same oscillation frequency (see Eq. 14), therefore we
introduce the global coefficient $a_{\ell}$ describing the energy contained in
each mode $\ell$ and defined by :

$\displaystyle a_{\ell}^{2}=\sum_{m=-\ell}^{\ell}a_{\ell,m}^{2},$ (37)

where $a_{\ell}$ is then positive by convention. The global surface
deformation is obtained by a sum over all modes $\ell$. We recall the
expression of the root mean squared deformation $\zeta_{\Omega}$:
$\zeta_{\Omega}^{2}(t)=\sum_{\ell=2}^{+\infty}\sum_{m=-\ell}^{\ell}a_{\ell,m}^{2}(t)$.
Figure 6a shows the time evolution of $\zeta_{\Omega}/R_{0}$ together with the
total amplitude of the mode 2, $a_{2}/R_{0}$. We observe a rapid linear rise,
consistent with the prediction for $t\ll t_{2},t_{c}$, before oscillations and
saturation of the total deformation around $\zeta_{\Omega}/R_{0}\approx 0.15$.
As shown in the figure, most of the deformation comes from the mode $\ell=2$.
This first observation is in agreement with various experimental studies on
large bubbles immersed in a turbulent flow, which have reported the dominance
of mode 2 deformation in bubble dynamics (Risso & Fabre, 1998; Ravelet et al.,
2011). In practice, a sum over the first three modes $\ell=2,3,4$ estimates
the amplitude deformation $\zeta_{\Omega}$ within less than 2% of error as
long as $\zeta_{\Omega}/R_{0}<1$ for all cases. It confirms the predominance
of the first modes of oscillation, and validates the spherical decomposition
approach. In the following, we focus on the modes $\ell=2,3$ and 4 which
contribute mostly to the surface deformation.

Figure 7: Bubble deformation and break-up as time progresses for the same
initial conditions but increasing Weber number. Time increases from left to
right while We increases from top to bottom. At We=3, this particular
condition does not lead to break-up. At We=6, break-up occurs for
$t/t_{c}\approx 1.5$, while break-up occurs earlier for We=15, at
$t/t_{c}\approx 0.8$. The initial stage of deformation, for $t<0.5t_{c}$
appears very similar at all Weber numbers.

Figure 6b shows the temporal evolution of $a_{\ell}^{2}$ for the first few
spherical harmonics $\ell\in[1,4]$ for a single run. The modes $\ell=1$ have
zero amplitude from the choice of the center position, and the modes $\ell=2$
dominate the surface energy deformation. The energy contained in modes
$\ell=2$ oscillates with time, at a frequency close to $2\omega_{2}$ and the
energies in the higher order modes $\ell=3,4$ and $5$ remain smaller at all
times. Figure 6c shows the temporal evolution of each individual harmonic
coefficient $a_{2,m}$ for different $m$ components. All 5 coefficients appear
to oscillate with comparable amplitude and with frequencies close to
$\omega_{2}$.

###  3.5 Deformation dynamics for increasing Weber number and unstable
conditions

For one particular configuration of the turbulent flow, we investigate the
role of the Weber number on the deformation by changing only the value of the
surface tension. Figure 7 shows snapshots of the bubbles for increasing Weber
number We = 3, 6 and 15, for the same initial turbulent conditions. For all
cases, we observe early growth of surface deformation, while the turbulent
background flow stays identical for all runs and independent of surface
tension, validating the no-feedback hypothesis introduced in section 2. For
all We, the short time deformation under the same turbulence conditions appear
identical in agreement with the theoretical description of the linear regime
for $t\ll t_{c}$. For stable conditions (here We=3), the deformation rapidly
saturates as previously described. For unstable conditions (here We=6 and 15),
the deformations eventually lead to break-up at various time. The highest
Weber number (We=15) displays a break-up event relatively early at
$t/t_{c}\approx 0.8$, while for conditions closer to the stability threshold,
break-up occurs at $t/t_{c}>1$.

A quantitative description of the Weber number influence for one particular
run is given in Figure 8. We compute the total deformation $\zeta_{\Omega}$
starting from the same initial condition of the flow configuration, and
increasing values of the Weber number We = 1.5, 3, 6, 15, 30 and 45. Figure 8
shows the deformation $\zeta_{\Omega}/R_{0}$ as a function of time $t/t_{c}$.
We observe a universal rapid linear growth, followed by a saturation that
depends on the Weber number. The rapid linear growth is independent of the
Weber number for $t/t_{c}<0.1$, as expected from Eq. 21. As time passes, the
curves diverge from each other, the lower Weber number curves exhibits earlier
sub-linear growth. The saturation leads to long time oscillations at the low
Weber numbers namely We=1.5, 3 and 6 during several eddy turnover times, with
increasing amplitude as the Weber number increases. For larger Weber number
(here We=15 and 45), the bubble experiences break-up during the first eddy
turn over time, when surface deformation reaches $\zeta_{\Omega}/R_{0}\approx
0.74$.

Figure 8: (a) Root mean square surface deformation $\zeta_{\Omega}$ as a
function of time $t/t_{c}$ for the same initial turbulent flow condition and
increasing Weber number We = 1.5, 3, 6, 15, 30. (b) Zoom-in on the early-stage
dynamics. For $t/t_{c}<0.1$, a universal linear increase of the deformation is
observed at all We number. Earlier saturation time is reached for higher
surface tension forces (lower We), leading to lower values of the saturated
deformation amplitude.

###  3.6 Ensemble-averaged deformations at one Weber number: linear growth and
saturated regime.

Figure 9: a) Individual realisations of $\zeta_{\Omega}$ at We = 3, as a
function of the dimensionless time $t/t_{c}$ in transparent colors. Ensemble
average $\langle\zeta_{\Omega}\rangle$ of the surface deformation is
superimposed (black full symbolds), together with the ensemble average mode 2
amplitude (red full symbol). b) Coefficient $a_{2}$ as a function of
dimensionless time $t/t_{c}$ for $\textrm{We}=3$ (individual realisations in
transparent colors). Ensemble average $\langle a_{2}\rangle$ is superimposed
(full red symbol). Fit by $A(1-e^{-t/t_{sat}})$ of the ensemble average
$a_{2}$ c) Coefficient $a_{3}$ as a function of dimensionless time $t/t_{c}$
for $\textrm{We}=3$ and equivalent exponential fit (individual realisations
are transparent colors and green full triangle is ensemble average $\langle
a_{3}\rangle$).

To obtain ensemble average quantities, we analysed between 10 to 15 runs for
each Weber number, starting from different turbulent flow configurations,
which are obtained from different times of the precursor simulations. For each
run, we perform the analysis described above, from the Voronoi decomposition
to the computation of the surface deformation $\zeta_{\Omega}$ and the first
harmonics coefficients $a_{2},a_{3}$ and $a_{4}$.

Figure 9a) shows the temporal evolution of each individual realisation of
$\zeta_{\Omega}/R_{0}$ performed at $\textrm{We}=3$, as a function of the
dimensionless time $t/t_{c}$. The ensemble average
$\langle\zeta_{\Omega}\rangle/R_{0}$ is superimposed in black, and follows a
linear increase, as predicted by eq. 22. The mode 2 amplitude $a_{2}/R_{0}$ is
superimposed in figure 9a) (${\color[rgb]{1,0,0}\square}$), confirming in a
statistical sense the predominance of the mode 2 of oscillations. The modes
$a_{2}/R_{0}$ and $a_{3}/R_{0}$ follow the same trend, as shown in figure 9b)
and c), each coloured curve being an individual realisation. Again, at short
time the linear regime predicted by eq. 22 is observed.

We observe a transition to a sub-linear regime for
$\langle\zeta_{\Omega}\rangle$, $\langle a_{2}\rangle$ and $\langle
a_{3}\rangle$ around $0.1<t<0.3$. This transition is well described by the
theoretical expression derived in section 2 (eq. 27) i.e. $\propto
t\sqrt{(1-(\kappa t)^{2})}$ with $\kappa=\omega_{2}/\sqrt{3}$ for
$\langle\zeta_{\Omega}\rangle$ or $\langle a_{2}\rangle$, and
$\kappa=\omega_{3}/\sqrt{3}$ for $\langle a_{3}\rangle$. The prefactor is
related to the turbulence velocity field statistics $u_{\ell,m}$, which will
be evaluated in an upcoming section. In order to evaluate the saturation with
time, we describe the temporal evolution of $\langle\zeta_{\Omega}\rangle$,
$\langle a_{2}\rangle$ and $\langle a_{3}\rangle$ by the empirical expression
$a_{\ell}^{\infty}\left(1-e^{-t/t_{sat}}\right)$ with two fit parameters
$a_{\ell}^{\infty}$ and $t_{sat}$. The early time evolution is proportional to
$a_{\ell}^{\infty}/t_{sat}$, while the characteristic time of saturation is
given by $t_{sat}$. Eventually, we perform for each ensemble average mode
$\langle a_{\ell}\rangle$ a linear fit $a_{\ell}=s_{\ell}t$ in the early time
evolution ($t<t_{sat}/2$) by a simple linear model, which gives an accurate
measure of the growth velocity $s_{\ell}$.

###  3.7 Ensemble averaged deformations for increasing Weber number

Figure 10: a) Ensemble-averaged deformation $\langle\zeta_{\Omega}\rangle$ as
a function of time for increasing We number. For the lower We, a fit by
$A(1-e^{-t/t_{sat}})$ captures the transition from linear growth to
saturation, with saturation values increasing with We. Open symbols are used
once half of the bubbles in the ensemble have broken. b) Slope at the origin
as a function of the Weber number for $\langle\zeta_{\Omega}\rangle$, $\langle
a_{2}\rangle$, $\langle a_{3}\rangle$, $\langle a_{4}\rangle$. The slope at
the origin for $\langle\zeta_{\Omega}\rangle$ is independent of We. c)
Amplitude of saturation of the global deformation and modes as a function of
We.

Figure 10 summarizes the surface deformation growth for increasing We number.
Figure 10a) shows the evolution of $\langle\zeta_{\Omega}\rangle$ averaged
over at least 10 runs by Weber number, as a function of the Weber number. The
evolution is shown in solid symbols until half of the bubbles have broken, and
in open symbols later. At early time, we show the existence of an universal
regime, with a linear increase in time of the deformations independent of the
Weber number We. Similar results are observed when considering the spherical
harmonics amplitude $\langle a_{2}\rangle$, $\langle a_{3}\rangle$ and
$\langle a_{4}\rangle$. The dimensionless growth velocity
$s_{\ell}t_{c}/R_{0}$ at short time, obtained from the linear fits is shown in
figure 10b) as a function of the Weber number for
$\langle\zeta_{\Omega}\rangle/R_{0}$ ($\circ$), $\langle a_{2}\rangle$ (red
square), $\langle a_{3}\rangle$ (green triangle) and $\langle a_{4}\rangle$
(blue pentagon). For the low Weber numbers We=1.5 and We = 3, only the modes
that are significantly above the noise level have been represented, excluding
modes 3 and 4 at We = 1.5, and mode 4 at We = 3. The growth velocity
$s_{\ell}$ depends on the mode $\ell$, but is independent of Weber number, as
expected from Eq. 14 and 22.

The departure from linear growth that occurs at a time increasing with Weber
number is compatible with the prediction given that the second mode reduced
period $t_{2}$. Eventually, the saturation can be measured using the
exponential fit for the lowest Weber number (We= 1.5, 3 and 6) where a
saturation is clearly visible. For higher Weber number, we consider the value
of the linear fitting model, evaluated at the time of the first break-up. The
saturation time for $\textrm{We}<10$ and the earliest break-up value for
$\textrm{We}>10$ is shown in figure 10c). We observe a continuous growth with
Weber number, premising the likelihood of bubble break-up. At higher Weber
number, we also remark that the relative amplitude of modes 3 and modes 4 are
increasing for both the growth velocity $s_{\ell}$ and the saturation
$a_{\ell}^{\infty}$. At high We number, the higher order modes of deformation
in the global dynamics becomes more important and is likely to influence
break-up geometry that the bubble will experience. Note that the saturation at
intermediate Weber number would require to be treated theoretically by
considering the role of non-linear and viscous effects. Note also that for
$\zeta_{\Omega}/R_{0}>0.25$, a significant part of the bubbles have broken,
and the mean deformation observed when considering all break-up events is
$\zeta_{\Omega}/R_{0}\approx 0.74$, which is compatible with the observations
on areas of deformation at break-up made by Risso & Fabre (1998).

Having access to the full velocity flow, we can now relate the surface
deformation growth parametrized by $a_{\ell}$ and $s_{\ell}$ to the statistics
of the surrounding turbulent flow.

###  3.8 Statistics of velocity increments on a sphere predict the linear
growth of deformation

Figure 11: (a) Ensemble average of spherical harmonics growth rate $\langle
s_{\ell}\rangle R_{0}/t_{c}$ for $\ell=2,3$ and 4, together with the
corresponding spherical harmonics amplitude $\langle u_{2}\rangle t_{c}/R_{0}$
of the turbulent velocity fluctuations at the bubble scale shown as red
circle. Black squares are $\langle s_{\ell}\rangle t_{c}/R_{0}$ averaged for
all We, while the colored-squares correspond to the increasing We number.
Theoretical prediction $a_{\ell}=u_{\ell}t$ from section 2 is well verified,
with velocity fluctuations providing a quantitative prediction of the bubble
deformation. The interface deformation statistics at early times can thus be
extracted from velocity fluctuations at the bubble scale. (b) Probability
distribution of the velocity fluctuations $\langle u_{2}\rangle$ corresponding
to the statistics of velocity fluctuations on the bubble sphere responsible
for mode 2 deformation. A fit by a $\Gamma$-distribution of expression
$f(x,k,p)=x^{k-1}e^{-x/p}/\Gamma(k)p^{k}$ is superimposed (red dashed line)
with $k=7.5$ and $p=0.135$.

From the theory developed in section 2, the pre-factor of the linear growth
regime of the spherical harmonics mode amplitude shall be related to the
statistics of the turbulent flow. Indeed, the growth of the spherical
harmonics coefficients $a_{\ell,m}$ are linked to their counterpart in the
turbulent fluctuations, ${u}_{\ell,m}$ and ${\beta}_{\ell,m}$ in equations 14,
22 and 27. It bears significant importance from a practical point of view, as
measuring velocity statistics in an experimental or natural turbulent flow is
relatively accessible while measuring deformation properties on a bubble
interface is much more challenging. We recall the description of the spherical
increments, $\delta u_{S}$, which can be decomposed in the spherical harmonics
base, $\delta
u_{S}=\sum_{\ell=2}^{+\infty}\sum_{m=-\ell}^{\ell}{u}_{\ell,m}Y_{\ell}^{m}(\theta,\varphi)$,
where ${u}_{\ell,m}$ depends on the radius $r$ of the sphere.

Using the velocity flow in the DNS, we computed the ensemble average $\langle
u_{\ell}\rangle$ of the harmonic coefficients of $\delta u_{S}(r)$ for
$r=R_{0}$, using the same procedure as for the $a_{\ell}$. The dimensionless
coefficients $\langle u_{\ell}\rangle t_{c}/R_{0}$ are shown in figure 11a),
for $\ell=2,3$ and $4$ and correspond to the intensity of velocity
fluctuations at the bubble scale for specific modes $\ell$. The corresponding
coefficients $s_{\ell}t_{c}/R_{0}$ obtained from the processing of surface
deformations are superimposed in red and correspond to the intensity of
deformations of the bubbles for the corresponding modes $\ell$. For modes
$\ell=2$, 3 and 4, we find a quantitative agreement between $\langle
a_{\ell,m}\rangle$ and $\langle u_{\ell,m}\rangle t$, as expected from section
2, eq. 22. For $\ell=4$, we observe a slight difference $s_{4}>u_{4}$, which
can be attributed to the limit of resolution of the harmonics coefficient
computation on the bubble deformation for higher $\ell$ values.

Beyond the equality of the ensemble average values, the prediction made in
section 2 shall be valid for each individual realisation. This equality can be
used to infer the full statistics of $a_{\ell,m}$ at short time, using the
statistics of the $u_{\ell,m}$ in the background flow.

We= |  1.5 |  3 |  6 |  15 |  30 |  45 |  Reλ = |  38
---|---|---|---|---|---|---|---|---
$\langle s_{2}\rangle R_{0}/t_{c}$ |  1.06 |  1 |  1.03 |  0.93 |  0.98 |  0.95 |  $\langle u_{2}\rangle t_{c}/R_{0}$ |  0.95
$\langle s_{3}\rangle R_{0}/t_{c}$ |  0.56 |  0.55 |  0.54 |  0.44 |  0.51 |  0.64 |  $\langle u_{3}\rangle t_{c}/R_{0}$ |  0.50
$\langle s_{4}\rangle R_{0}/t_{c}$ |  0.44 |  0.35 |  0.34 |  0.31 |  0.33 |  0.35 |  $\langle u_{4}\rangle t_{c}/R_{0}$ |  0.26
$s_{2}^{\prime}R_{0}/t_{c}$ |  0.41 |  0.36 |  0.33 |  0.32 |  0.33 |  0.27 |  $u_{2}^{\prime}t_{c}/R_{0}$ |  0.34
$s_{3}^{\prime}R_{0}/t_{c}$ |  0.23 |  0.24 |  0.23 |  0.15 |  0.23 |  0.27 |  $u_{3}^{\prime}t_{c}/R_{0}$ |  0.17
$s_{4}^{\prime}R_{0}/t_{c}$ |  0.08 |  0.06 |  0.08 |  0.08 |  0.10 |  0.09 |  $u_{4}^{\prime}t_{c}/R_{0}$ |  0.09
Table 2: Values of the coefficients shown in figure 10 and 11. Ensemble
averaged of spherical harmonics growth rate $\langle s_{\ell}\rangle
R_{0}/t_{c}$ for $\ell=2,3$ and 4. The corresponding spherical harmonics
amplitude $\langle u_{2}\rangle t_{c}/R_{0}$ of the turbulent velocity
fluctuations at the bubble scale are also provided, computed for Reλ=38.

Figure 11b) shows the probability density function of the dimensionless mode 2
$u_{2}t_{c}/R_{0}$ at the bubble scale $R_{0}$. In contrary to the two-points
velocity increment $D_{LL}(R_{0})$, the pdf of $u_{2}$ in the inertial range
does not exhibit large tails, and 95% of $u_{2}$ values lie in the range
$[0.5\langle u_{2}\rangle,2\langle u_{2}\rangle]$. The difference between the
large tails of the pdf of $D_{LL}$ at $d/\lambda=1.5$ and the short tails of
$u_{2}$ can be attributed to the spatial average operation on the sphere,
which smooths out all the intermittent structures at a scale smaller than the
bubble, and suggests that the flow intermittency has a limited influence on
the bubble deformation in the inertial range. Considering a positive definite
quantity, a fit of the probability density function by a $\Gamma$-distribution
of expression $f(x,k,p)=x^{k-1}e^{-x/p}/\Gamma(k)p^{k}$ gives $k=7.5$ and
$p=0.135$.

##  4 Implications for bubble lifetime statistics

Figure 12: (a) Estimated probability distribution function of bubble lifetime
in the limit of large Weber number, obtained from
$N(T/t_{c})=N(\zeta_{c}/u_{2}t_{c})$ (eq. 4.1), and valid at high Weber
number, i.e. $\textrm{We}\approx>30$, using $\zeta_{c}=0.74R_{0}$. Mean and
rms values are given in table 3. (b) Inverse of the bubble lifetime $t_{c}/T$,
and its standard deviation, similar to the bubble break-up frequency, as a
function of the We number. The mean values for 10 to 20 simulations for each
Weber number are used, and the error-bar corresponds to the standard deviation
in the lifetime. Dashed line corresponds to the model proposed by Martinez-
Bazan et al. (1999, 2010),
$t_{c}/T\propto\sqrt{1-\textrm{We}_{c}/\textrm{We}}$. At high We number, the
bubble lifetime can be predicted from the turbulence statistics, and given by
$\zeta_{c}/u_{2}t_{c}$

Given the linear growth of deformation with time and a growth rate given by
the mode $u_{2}$, the pdf of $u_{2}$ can be used to evaluate the bubble
lifetime distribution at high Weber number, when surface tension becomes
negligible. In the limit of inertial break up, with correlated velocity
fluctuations The distribution of lifetime $N(T/t_{c})$ at high Weber number is
given by :

$\displaystyle
N(T/t_{c})=N\left(\frac{\langle\zeta_{c}\rangle}{u_{2}t_{c}}\right),$ (38)

where $\langle\zeta_{c}\rangle$ is the average critical deformation of
$\zeta_{\Omega}$ at break up. $\zeta_{c}$ cannot be computed from the Voronoi
decomposition of the interface, since the radius becomes multivalued at high
deformation. To evaluate the critical deformation $\zeta_{c}$, we approximate
the shape at break-up by an ellipsoid of the same volume than the initial
bubble and whose longest axe corresponds to the maximum distance between two
points on the bubble interface. From the large Weber number data
($\textrm{We}\geq 15$), we find $\langle\zeta_{c}\rangle=0.74$. Figure 12a
shows the distribution $N(T/t_{c})$ obtained from eq. 4.1, for a bubble of
size $d_{0}/\lambda=1.5$ using the statistics of the mode $u_{2}$ of the
harmonic decomposition of $u_{T}$ for a Taylor Reynolds number
$Re_{\lambda}=38$ (the turbulence considered in the present paper). The lower
lifetime bound is 0.1$t_{c}$, while the distribution shows exponential tales
for large values of the bubble lifetime. We also compute the same predicted
lifetime statistics for a turbulent flow obtained similarly but with higher
Reynolds number, $Re_{\lambda}=76$, and test values of $d_{0}/\lambda=1.5$
(red line on fig 12a) and $d_{0}/\lambda=3$ and observe only little changes in
the distribution. This suggests that the lifetime statistics is fairly
insensitive to the Reynolds number for bubbles within the inertial range
($d_{0}/\lambda>1$).

The average bubble lifetime can be computed from the numerical runs for each
Weber number from the break-up time of the initial bubble. Figure 12b shows
the inverse of the bubble lifetime $t_{c}/T$, similar to the break-up
frequency introduced in the literature, for increasing Weber number. The
errorbars correspond to the rms values $T_{rms}/t_{c}$. We observe that the
bubble lifetime decreases with Weber number, and reaches a value independent
of Weber number for $\textrm{We}\gg\textrm{We}_{c}$, here above We=30. We
observe that our data can be reasonably well be described by the model
proposed by Martinez-Bazan et al. (1999, 2010) (dashed line),
$t_{c}/T=C_{g}\sqrt{1-\textrm{We}_{c}/\textrm{We}}$., with the saturation
value at high Weber number being the only adjusted parameter. We use
$\textrm{We}_{c}=3$ and $C_{g}=1.5$ is fitted to our numerical data, whereas
Martinez-Bazan et al. (1999, 2010) initial data used $C_{g}=0.673$ and
$\textrm{We}_{c}=1$ to fit their data at Weber numbers far from the threshold
$\textrm{We}_{c}$.

|  We = 3 |  We = 6 |  We = 15 |  We = 30 |  We = 45 |  Re${}_{\lambda}=38$
---|---|---|---|---|---|---
$\\#$ of elements |  45 (13) |  39 (14) |  20 (8) |  20 (7) |  20 (7) |
$T/t_{c}$ |  3.01 |  2.15 |  1.06 |  0.90 |  0.88 |  $\vline~{}\langle\zeta_{c}\rangle/(\langle u_{2}\rangle t_{c})=0.91$
$T_{rms}/t_{c}$ |  0.86 |  1.18 |  0.65 |  0.45 |  0.40 |  $\vline~{}\langle\zeta_{c}\rangle/(u_{2}^{\prime}t_{c})=0.41$
$\langle\zeta_{c}\rangle/R_{0}$ |  0.25 |  0.5 |  0.68 |  0.85 |  0.69 |
Table 3: Bubble life-time obtained as a function of Weber number, obtained
from the DNS and from the turbulence statistics. The mean $T$ and rms
$T_{rms}$ values of the bubble life-time obtained from ensembles of 20 to 45
simulations at increasing Weber number (number of elements indicated in first
row). We also provide the computed mean deformation at break-up
$\langle\zeta\rangle$. The number in parentheses corresponds to the number of
elements used to compute the deformation at break-up $\zeta_{c}$ from the
bubble deformation simulations. On the right of the table, we show the
predicted mean and rms lifetime from the turbulence statistics,
$\zeta_{c}/(\langle u_{2}\rangle t_{c})$, for a constant deformation threshold
$\zeta_{c}/R_{0}=0.74$, for Re${}_{\lambda}=38$. Excellent agreement between
the numerical values and the ones predicted by the turbulence statistics is
observed at high Weber number.

Table 3 summarizes the bubble lifetime and the estimated lifetime from
velocity statistics. The mean and standard deviation values of the inferred
lifetime from the turbulence statistics can be compared with the bubble
lifetime obtained from the direct numerical simulations for ensembles of 20 to
45 simulations. We observe excellent agreement between the inferred and
simulated mean lifetime, as well as for the rms value of the distribution,
considering a constant break-up threshold of $\zeta_{c}=0.74R_{0}$, inferred
from the bubble deformation at break-up for high Weber number cases
($\textrm{We}\geq 15$). The full distribution
$N(T/t_{c})=N\left(\langle\zeta_{c}\rangle/(u_{2}t_{c})\right)$ could hence be
a good estimate of bubble lifetime statistics in the limit of high Weber
number.

##  5 Conclusion

We have presented a theoretical framework for bubble deformation in a
turbulent flow, by performing a spherical harmonics decomposition of the
bubble deformation and deriving a general forced oscillator equation for these
spherical harmonics mode, where each mode is forced by the corresponding
turbulent fluctuations mode. We identify various regimes in the time evolution
of the deformation, in particular a short time scale regime where deformations
grow linearly with time, the pre-factor being given by the strength of the
turbulent flow.

We perform direct numerical simulations of bubbles evolving in homogeneous and
isotropic turbulent flow at increasing We number and verify the theoretical
predictions. We observe the linear regime of deformation growth for individual
events and ensemble averaged quantities and show that the pre-factor of the
linear growth is indeed given by the statistics of the turbulence deformation
at the bubble scale. We observe in the simulations that the level of
deformation saturates, with a saturation time and level of saturation that
depends on the Weber number, with lower Weber number saturating earlier and at
lower deformation. At low to intermediate Weber number, the eigen-mode 2
dominates the deformation. At lower We number, a much broader lifetime
distribution is observed, with break-up that can occur at much later time, as
multiple eddies participate in the deformation and break-up, as described in
Risso & Fabre (1998). The broad lifetime distribution can be interpreted as a
consequence of a stochastic process with threshold, which gives rise to long
oscillation before a fluctuation of higher amplitude lead to break-up. Further
stochastic analysis of the surface deformation fluctuations around the average
saturation, accounting for non-linear and viscous effects could provide a
quantitative prediction of the lifetime distribution at intermediate Weber
number.

At high Weber number, break-up occurs within a single eddy turn-over time, and
such close to the linear regime of deformation. These findings have
significant implications: the ability to predict bubble deformation and break-
up time from the turbulence statistics is now possible, by extracting the
turbulence statistics and computing spherical harmonics decomposition of the
corresponding mode 2 and 3. According to our theory and as shown in the final
section, the statistics of such modes directly provide the bubble lifetime, if
one assumes a maximum deformation threshold. Our data typically suggest
$\zeta_{\Omega}/R_{0}\approx 0.74$. The deformation also shows significant
increase of higher frequency modes, which could be related at high Weber
number to the appearance of much smaller child bubbles. More elaborate model
would require a better estimate of the deformation statistics at break-up.

Experimental work probing the full three dimensional turbulence, while
challenging, could further test these results on bubble deformation, bubble
lifetime and break-up in turbulence, together with simulations at higher
turbulent Reynolds number investigating the role of the ratio between the
turbulent length scale and the bubble deformation dynamics. Similarly, this
numerical framework could be used to explore the child bubble size
distribution resulting from bubble break-up and its relationship with the
turbulence statistics.

## Appendix : Convergence test

We present a convergence test on the interface resolution, by running
simulations for maximum interface refinement levels $N=9$ and $N=10$, using
the same precursor simulation. As discussed in section 3, these resolutions
correspond to 70 and 140 grid points per diameter, while the resolution of the
turbulent flow remains the same, and correspond to the $Re_{\lambda}=38$ flow
with a velocity refinement of level 7. The interface deformation
$\zeta_{\Omega}$ as a function of time is shown in figure 12, for increasing
Weber number. We observe quasi-identical results between the maximum levels 9
and 10, for all cases, for times of the order of the eddy turn-over time. This
shows that the interface deformation is correctly resolved at the chosen
resolution, and the results presented in the paper are independent of
numerical resolution. It is important to remember that the increase in
resolution here only applies to the interface resolution, while the resolution
on the turbulent velocity field remains the same. However, we have
independently verified that the statistical properties of the turbulent flow
are also converged at the resolution we are using (see figure 3).

Figure 13: Convergence test on $\zeta_{\Omega}$ as a function of dimensionless
time $t/t_{c}$, for We=3,6,15,30,45. Level $N=9$ (blue) and level $N=10$
(red). Excellent convergence agreement is observed between the two numerical
resolutions.

###### Acknowledgements.

We thank Stéphane Popinet for scientific discussion and development of the
Basilisk package. We thank the three anonymous reviewers whose comments have
helped improve the quality of the manuscript. This work was supported by the
NSF CAREER award 1844932 to L.D., and the American Chemical Society Petroleum
Research Fund Grant 59697-DNI9 to L.D. A.R. was supported by an International
Fund grant from Princeton University to L.D. S.P. and A.R. were supported by
the Labex ENS-ICFP. We would like to acknowledge high-performance computing
support from Cheyenne (doi:10.5065/D6RX99HX) provided by NCAR’s Computational
and Information Systems Laboratory, sponsored by the National Science
Foundation. Computations were also performed on the Princeton supercomputer
Tiger2, as well as on Stampede, through XSEDE allocations to L.D. and W.M.,
XSEDE is an NSF funded program 1548562.

Declaration of Interests. The authors report no conflict of interest.

## References

  * Andersson & Andersson (2006) Andersson, Ronnie & Andersson, Bengt 2006 On the breakup of fluid particles in turbulent flows. AIChE Journal 52 (6), 2020–2030.
  * Balachandar & Eaton (2010) Balachandar, S & Eaton, John K 2010 Turbulent dispersed multiphase flow. Annual review of fluid mechanics 42, 111–133.
  * Berny et al. (2020) Berny, Alexis, Deike, Luc, Séon, Thomas & Popinet, Stéphane 2020 Role of all jet drops in mass transfer from bursting bubbles. Physical Review Fluids 5 (3), 033605.
  * Brocchini & Peregrine (2001) Brocchini, M. & Peregrine, D. H. 2001 The dynamics of strong turbulence at free surfaces. part 1. description. J. Fluid Mech. 449, 225–254.
  * Cano-Lozano et al. (2016) Cano-Lozano, José Carlos, Martinez-Bazan, Carlos, Magnaudet, Jacques & Tchoufag, Joël 2016 Paths and wakes of deformable nearly spheroidal rising bubbles close to the transition to path instability. Physical Review Fluids 1 (5), 053604.
  * Caroli et al. (2009) Caroli, M., de Castro, P. M. M., Loriot, S., Rouiller, O., Teillaud, M. & Wormser, C. 2009 Robust and efficient delaunay triangulations ofpoints on or close to a sphere. Research Report RR-7004. INRIA.
  * Cowen & Variano (2008) Cowen, E. A. & Variano, E. A. 2008 A random-jet-stirred turbulence tank. J. Fluid Mech. 604, 1–32.
  * Deike et al. (2018) Deike, L., Ghabache, E., Liger-Belair, G., Das, A. K., Zaleski, S., Popinet, S. & Seon, T. 2018 The dynamics of jets produced by bursting bubbles. Phys. Rev. Fluids 3.
  * Deike & Melville (2018) Deike, Luc & Melville, W Kendall 2018 Gas transfer by breaking waves. Geophysical Research Letters 45 (19), 10–482.
  * Deike et al. (2016) Deike, Luc, Melville, W. Kendall & Popinet, Stéphane 2016 Air entrainment and bubble statistics in breaking waves. Journal of Fluid Mechanics 801, 91–129.
  * Eggers & Villermaux (2008) Eggers, Jens & Villermaux, Emmanuel 2008 Physics of liquid jets. Reports on progress in physics 71 (3), 036601.
  * Elghobashi (2019) Elghobashi, Said 2019 Direct numerical simulation of turbulent flows laden with droplets or bubbles. Annual Review of Fluid Mechanics 51 (1), 217–244.
  * Ern et al. (2012) Ern, Patricia, Risso, Frédéric, Fabre, David & Magnaudet, Jacques 2012 Wake-induced oscillatory paths of bodies freely rising or falling in fluids. Annual Review of Fluid Mechanics 44, 97–121.
  * Fuster & Popinet (2018) Fuster, Daniel & Popinet, Stéphane 2018 An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension. Journal of Computational Physics 374, 752–768.
  * Harper (1970) Harper, JF 1970 On bubbles rising in line at large reynolds numbers. Journal of Fluid Mechanics 41 (4), 751–758.
  * Hinze (1955) Hinze, J.O. 1955 Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes. AIChE Journal 1 (3), 289–295.
  * van Hooft et al. (2018) van Hooft, J Antoon, Popinet, Stéphane, van Heerwaarden, Chiel C, van der Linden, Steven JA, de Roode, Stephan R & van de Wiel, Bas JH 2018 Towards adaptive grids for atmospheric boundary-layer simulations. Boundary-layer meteorology 167 (3), 421–443.
  * Hussein et al. (1994) Hussein, H. J., Capp, S. P. & George, W. K. 1994 Velocity measurements in high-reynolds number, momentum-conserving, axisymmetric, turbulent jets. J. Fluid Mech. 258, 31–75.
  * Kolmogorov (1941) Kolmogorov, A. N. 1941 The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers. Dokl. Akad. Nauk SSSR 30:301.
  * Lai et al. (2018) Lai, Ching-Yao, Eggers, Jens & Deike, Luc 2018 Bubble bursting: Universal cavity and jet profiles. Phys. Rev. Lett. 121, 144501\.
  * Lance & Bataille (1991) Lance, M & Bataille, J 1991 Turbulence in the liquid phase of a uniform bubbly air–water flow. Journal of fluid mechanics 222, 95–118.
  * Loisy & Naso (2017) Loisy, Aurore & Naso, Aurore 2017 Interaction between a large buoyant bubble and turbulence. Phys. Rev. Fluids 2, 014606.
  * Magnaudet & Eames (2000) Magnaudet, Jacques & Eames, Ian 2000 The motion of high-reynolds-number bubbles in inhomogeneous flows. Annual Review of Fluid Mechanics 32 (1), 659–708.
  * Martinez-Bazan et al. (1999) Martinez-Bazan, C., Montanes, J.L. & Lasheras, J.C. 1999 On the breakup of an air bubble injected into a fully developed turbulent flow. part 1. breakup frequency. J. Fluid Mech. 401, 157–182.
  * Martinez-Bazan et al. (2010) Martinez-Bazan, Carlos, Rodriguez-Rodriguez, Javier, Deane, GB, Montañes, JL & Lasheras, JC 2010 Considerations on bubble fragmentation models. Journal of Fluid Mechanics 661, 159–177.
  * Masuk et al. (2019) Masuk, Ashik Ullah Mohammad, Salibindla, Ashwanth, Tan, Shiyong & Ni, Rui 2019 V-onset (vertical octagonal noncorrosive stirred energetic turbulence): A vertical water tunnel with a large energy dissipation rate to study bubble/droplet deformation and breakup in strong turbulence. Review of Scientific Instruments 90 (8), 085105.
  * Mathai et al. (2020) Mathai, Varghese, Lohse, Detlef & Sun, Chao 2020 Bubble and buoyant particle laden turbulent flows. Annu. Rev. Condens. Matter Phys 11.
  * Maxworthy et al. (1996) Maxworthy, T, Gnann, C, Kürten, M & Durst, F 1996 Experiments on the rise of air bubbles in clean viscous liquids. Journal of fluid mechanics 321, 421–441.
  * Miller & Scriven (1968) Miller, CA & Scriven, LE 1968 The oscillations of a fluid droplet immersed in another fluid. Journal of fluid mechanics 32 (3), 417–435.
  * Moore (1965) Moore, DW 1965 The velocity of rise of distorted gas bubbles in a liquid of small viscosity. Journal of Fluid Mechanics 23 (4), 749–766.
  * Mostert & Deike (2020) Mostert, W & Deike, L 2020 Inertial energy dissipation in shallow-water breaking waves. Journal of Fluid Mechanics 890.
  * Mostert et al. (2020) Mostert, W, Popinet, S & Deike, L 2020 Transition to turbulence in ocean breaking waves. Submitted .
  * Peregrine (1976) Peregrine, D. H. 1976 Interaction of water waves and currents. Advances in Applied Mechanics. 16, 9–117.
  * Perrard et al. (2019) Perrard, S., Lozano-Durán, A., Rabaud, M., Benzaquen, M. & Moisy, F. 2019 Turbulent windprint on a liquid surface. J. Fluid Mech 873, 1020–1054.
  * Phillips (1957) Phillips, O. M. 1957 On the generation of waves by turbulent wind. J. Fluid Mech. 2 (05), 417–445.
  * Pope (2000) Pope, S. B. 2000 Turbulent Flows. Cambridge University Press.
  * Popinet (2009) Popinet, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. Journal of Computational Physics 228, 5838–5866.
  * Popinet (2015) Popinet, Stéphane 2015 A quadtree-adaptive multigrid solver for the serre–green–naghdi equations. Journal of Computational Physics 302, 336 – 358.
  * Popinet (2018) Popinet, Stéphane 2018 Numerical models of surface tension. Annual Review of Fluid Mechanics 50, 49–75.
  * Prosperetti (1980) Prosperetti, A. 1980 Free oscillations of drops and bubbles: the initial-value problem. J. Fluid Mech. 100, 333–347.
  * Ravelet et al. (2011) Ravelet, F., Colin, C. & Risso, F. 2011 On the dynamics and breakup of a bubble immersed in a turbulent flow. Phys. Fluids 23, 103301.
  * Reichl & Deike (2020) Reichl, Brandon G & Deike, Luc 2020 Contribution of sea-state dependent bubbles to air-sea carbon dioxide fluxes. Geophysical Research Letters p. e2020GL087267.
  * Risso (2018) Risso, Frédéric 2018 Agitation, mixing, and transfers induced by bubbles. Annual Review of Fluid Mechanics 50, 25–48.
  * Risso & Fabre (1998) Risso, Frédéric & Fabre, Jean 1998 Oscillations and breakup of a bubble immersed in a turbulent field. Journal of Fluid Mechanics 372, 323–355.
  * Rosales & Meneveau (2005) Rosales, Carlos & Meneveau, Charles 2005 Linear forcing in numerical simulations of isotropic turbulence: Physical space implementations and convergence properties. Physics of Fluids 17 (9), 095106.
  * Ruth et al. (2019) Ruth, D., Mostert, W., Perrard, S. & Deike, L. 2019 Bubble pinch-off in turbulence. PNAS .
  * Vejražka et al. (2018) Vejražka, Jiří, Zedníková, Mária & Stanovskỳ, Petr 2018 Experiments on breakup of bubbles in a turbulent flow. AIChE Journal 64 (2), 740–757.
  * Veron (2015) Veron, F. 2015 Ocean spray. Ann. Rev. Fluid Mech. 47, 507–538.
  * Villermaux & Bossa (2009) Villermaux, Emmanuel & Bossa, Benjamin 2009 Single-drop fragmentation determines size distribution of raindrops. Nature Physics 5 (9), 697.
  * Voronoi (1908) Voronoi, G. 1908 Nouvelles applications des paramètres continus à la théorie des formes quadratiques. J. Reine Angew. Math. 133, 97–178.
  * Yuan & Prosperetti (1994) Yuan, H & Prosperetti, A 1994 On the in-line motion of two spherical bubbles in a viscous fluid. Journal of Fluid Mechanics 278, 325–349.