2.1. Bubble oscillation in a fluid at rest

We consider an incompressible, initially spherical gas bubble of radius $R_0=d_0/2$, density $\rho _g$ and viscosity $\mu _g$, immersed in a liquid at rest of liquid density $\rho _\ell$, viscosity $\mu _l$ and surface tension between the gas and the liquid phase $\gamma$. The origin of space is set at the centre of mass of the bubble.

The chosen frame of reference is then non-inertial. Note that we do not study here the motion of the centre of mass but aim at a description of the bubble deformation. Focusing on the early stage, we limit our description to a bubble surface that can be parametrized in spherical coordinates $(r,\theta ,\varphi )$ by the single-valued distance ${\boldsymbol R} = (R_0(t) + \zeta (\theta ,\varphi ,t)) \hat {{\boldsymbol r}}$ to the centre (figure 1). We aim for a phenomenological set of equations for the modes of deformation of the bubble in the presence of a turbulent background flow, inspired by the forced oscillator model from Risso & Fabre (Reference Risso and Fabre1998) and Lalanne et al. (Reference Lalanne, Masbernat and Risso2019).

Figure 1.

Bubble deformation in a turbulent flow: sketch and definitions.

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

(2.1)\begin{equation} {H}_{\ell,m}(\theta,\varphi) = \sqrt{\frac{(2\ell+1)(\ell-m)!}{(\ell+m)!}} P^m_\ell(\cos{\theta}) \,\textrm{e}^{\textrm{i}m\varphi}, \end{equation}

with $\ell \in \boldsymbol {N}$ and $m \in [-\ell , \ell ]$, and $P^m_\ell$ is the associated Legendre polynomial. The real form ${Y}_{\ell ,m}(\theta ,\varphi )$ of the spherical harmonics reads as

(2.2)\begin{equation} Y_{\ell}^{m} = \begin{cases} \sqrt{2}({-}1)^m \textbf{Im}(H_{\ell}^{{-}m}) & \mbox{if } m<0, \\ H_\ell^0 & \mbox{if } m=0, \\ \sqrt{2}({-}1)^m \textbf{Re}(H_{\ell}^{m}) & \mbox{if } m>0. \end{cases} \end{equation}

The real spherical harmonics $Y_{\ell }^{m}$ verifies $\forall \ell \in \mathbb {N},\ \forall m \in [-\ell ,\ell ], \nabla ^2 r^{-\ell } Y_{\ell }^{m} = 0$ for $r>0$. In this basis, the expression of the bubble surface $R(\theta ,\varphi ,t)$ reads as

(2.3)\begin{equation} R(\theta,\varphi,t) = R_0(t) + \sum_{\ell=2}^{+\infty} \sum_{m={-}\ell}^\ell a_{\ell,m}(t) Y_{\ell}^{m}(\theta,\varphi), \end{equation}

where $R_0(t)$ is the mean bubble radius. In the following, we neglect the radius variation associated with gas compressibility and we assume that $R_0$ is independent of time. In the absence of any background flow, at the linear order in the deformation amplitude $\zeta /R_0$ and small viscous dissipation in the bubble boundary layer, the coefficients $a_{\ell ,m}$ are solutions of a set of oscillator equations of the form (Prosperetti Reference Prosperetti1980)

(2.4)\begin{equation} \ddot a_{\ell,m} + 2 \beta_\ell \dot a_{\ell,m} + (\ell-1)(\ell+1)(\ell+2) \omega_0^2 a_{\ell,m} = 0, \end{equation}

where ${\cdot}$ denotes the time derivative, $\omega _0 = \sqrt {\gamma /(R_0^3\rho _\ell )}$ is the typical angular frequency associated with bubble oscillation and $\beta _\ell$ is the damping coefficient of mode $\ell$. The expression of $\beta _2$ is summarized in Lalanne et al. (Reference Lalanne, Masbernat and Risso2019) for low viscosity, and the general formula can be found in Miller & Scriven (Reference Miller and Scriven1968). We introduce $\omega _\ell = \omega _0 \sqrt {(\ell -1)(\ell +1)(\ell +2)}$ the angular frequency associated to mode $\ell$ and the associated characteristic time $t_\ell = 1/\omega _\ell$ of mode $\ell$. We have, in particular, $\omega _2 = 2 \sqrt {3} \omega _0$.

2.2. Bubble oscillation in a turbulent background flow

We consider that the bubble is immersed in a liquid animated by turbulent motion, with a velocity $\boldsymbol u(\boldsymbol r,t)$ and a pressure field $p(\boldsymbol r,t)$. In the general case, modelling the action of the turbulent background flow on a bubble is particularly complex. A two-way coupling between the velocity fluctuations near the bubble interface and the interface deformation may take place. The turbulent flow deforms the bubble under the action of normal and tangential stresses, and the turbulent flow statistics can be modified in the immediate vicinity of the bubble interface. In order to construct a phenomenological model for the modes of deformation of the bubble, we follow the approach of Risso & Fabre (Reference Risso and Fabre1998) and Lalanne et al. (Reference Lalanne, Masbernat and Risso2019), which considers a dynamical equation for the second mode of deformation $a_2$, assumed to be the main mode excited by the flow. The amplitude $a_2(t)$ then represents the oscillation of the bubble between an oblate and a prolate shape.

We introduce two adaptations of the formulation of Risso & Fabre (Reference Risso and Fabre1998) and Lalanne et al. (Reference Lalanne, Masbernat and Risso2019). (i) We consider that the turbulent flow can excite other modes $(\ell ,m)$ of deformation. The contribution of these higher-order modes $(\ell > 2)$ could be non-negligible for high Weber numbers ($We \gg We_c$). (ii) We aim to estimate the turbulent excitation in terms of pressure and velocity statistics evaluated on a sphere of the same radius, in the bulk of the turbulent flow. We perform a linear expansion in the amplitude of the modes of deformation, and we neglect the effect of the interface on the turbulent statistics in its immediate vicinity. The decoupling approach was first motivated by the high-Weber-number case, in which capillary forces are small compared with inertial forces. For a lower Weber number ($We \approx We_c$), our approach is expected to still hold at short time compared with the response time $1/\omega _\ell$. We also consider a large dissipative time compared with the bubble period of oscillation $(\omega _\ell /\beta _\ell \gg 1)$, such that we neglect the viscous dissipation. For all these reasons, the following approach is rather phenomenological and by no means exact. The validity of our approach is further discussed by comparison with direct numerical simulations of a single bubble of initial spherical shape deformed by a homogeneous and isotropic turbulent flow.

At the linear order of deformation ($a_{\ell ,m}/R_0 \ll 1$), we assume that the bubble surface can still be described by a set of oscillatory equations, adding a forcing term that represents the turbulent background flow

(2.5)\begin{equation} \ddot a_{\ell,m} + \omega_\ell^2 a_{\ell,m} = {F}_{\ell,m}, \end{equation}

where ${F}_{\ell ,m}$ is a function of $\boldsymbol u(\boldsymbol R,t)$, $\textrm {d}\boldsymbol u/\textrm {d} t(\boldsymbol R,t)$ and $p(\boldsymbol R,t)$, respectively, the velocity field, the acceleration field and the pressure field evaluated at the surface of the bubble. We neglect the shear contribution and, looking at the early-time dynamics, we neglect viscous dissipation. To specify the function ${F}_{\ell ,m}$, we introduce the spherical harmonic decomposition of $u_r(R_0 \hat {{\boldsymbol r}},t) = {\boldsymbol u} \boldsymbol{\cdot} {\hat {{\boldsymbol r}}}$ and $p(R_0 \hat {{\boldsymbol r}},t)$ evaluated on a sphere of radius $R_0$,

(2.6)$$\begin{gather} u_r(R_0 \hat{{\boldsymbol r}},t) = \sum_{\ell=2}^{+\infty} \sum_{m={-}\ell}^\ell {u}_{\ell,m}(t) Y_{\ell}^{m}(\theta,\varphi), \end{gather}$$

(2.7)$$\begin{gather}p(R_0 \hat{{\boldsymbol r}},t) = \rho_\ell \sum_{\ell=2}^{+\infty} \sum_{m={-}\ell}^\ell {\rm \pi}_{\ell,m}(t) Y_{\ell}^{m}(\theta,\varphi), \end{gather}$$

where $\hat {{\boldsymbol r}}$ is the unit vector along $\boldsymbol r$. We then approximate the fields $\boldsymbol u(\boldsymbol R,t)$, $\textrm {d}\boldsymbol u/\textrm {d} t(\boldsymbol R,t)$ and $p(\boldsymbol R,t)$ by their value on an undeformed bubble, whose interface position is $R_0 \hat {{\boldsymbol r}}$. Doing so, we neglect the contribution of cross-terms of the type $a_{\ell ,m} u_{\ell ,m}$ and terms of $a_{\ell ,m}$, $u_{\ell ,m}$, ${\rm \pi} _{\ell ,m}$ of quadratic and higher order. In this limit, the coefficients of the spherical harmonic decomposition of $\dot u_r(R_0 \hat {{\boldsymbol r}},t)$ become ${\dot u}_{\ell ,m}$ as $R_0 \hat {{\boldsymbol r}}$ is independent of time. We also neglect the coupling between two different spherical harmonic modes $(\ell ,m)$ and $(\ell ',m')$ which arises at least at the order $a_{\ell ,m} u_{(\ell ',m')}$ or $a_{\ell ,m} a_{(\ell ',m')}$. For $a_{\ell ,m}/R_0 \ll 1$, we then assume that ${F}_{\ell ,m}$ is a linear function of $u_{\ell ,m}, {\dot u}_{\ell ,m}$ and ${{\rm \pi} }_{\ell ,m}$, and we postulate the following form:

(2.8)\begin{equation} \ddot a_{\ell,m} + \omega_\ell^2 a_{\ell,m} = {\dot u}_{\ell,m} + \frac{C(\ell)}{R_0} {\rm \pi}_{\ell,m}. \end{equation}

Here $C$ is a function of $\ell$. The kinematic boundary conditions at $t=0$ give $R(\theta ,\varphi ,0) = R_0$ and $\dot R(\theta ,\varphi ,0) = u_r(R_0,\theta ,\varphi ,0)$, which gives the initial conditions for $a_{\ell ,m}$,

(2.9a,b)\begin{equation} a_{\ell,m}(0) = 0; \quad \text{and}\quad \dot a_{\ell,m}(0) = u_{\ell,m}(0). \end{equation}

2.3. Link with surface deformation

We seek for an estimate of the average global surface deformation. We introduce the root mean squared $\zeta _\varOmega$ of the bubble interface from an integration over all solid angles,

(2.10)\begin{equation} \zeta_\varOmega^2 = \frac{1}{4{\rm \pi}} \iint \textrm{d}\varOmega \zeta^2. \end{equation}

According to the orthogonality of spherical harmonics $({1}/{4{\rm \pi} }) \iint \textrm {d}\varOmega Y_{\ell }^{m} Y_{\ell '}^{m'} = \delta _{\ell ,\ell '} \delta _{m,m'}$, $\zeta _\varOmega$ can be expressed as

(2.11)\begin{equation} \zeta_\varOmega^2(t) = \sum_{\ell=2}^{+\infty} \sum_{m={-}\ell}^\ell a_{\ell,m}(t)^2. \end{equation}

The surface deformation is then related to the solutions of (2.8). The coefficients $a_{\ell ,m}$ being stochastic variables, we look for their temporal behaviour as well as temporal behaviour of their ensemble-averaged values $\langle a_{\ell ,m}^2 \rangle$.

2.4. Controlling time scales

The solutions of (2.8) depend on the various time scales that are inherently associated to a 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 (Risso & Fabre Reference Risso and Fabre1998; Ravelet et al. Reference Ravelet, Colin and Risso2011). Another relevant time scale lies in the turbulence correlation time of the fluctuations at the bubble scale also discussed by several authors (Risso & Fabre Reference Risso and Fabre1998; Martinez-Bazan et al. Reference Martinez-Bazan, Montanes and Lasheras1999; Ravelet et al. Reference Ravelet, Colin and Risso2011). Depending on the Weber number, another time scale is eventually given either by the bubble lifetime or the 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 _\varOmega$ are set by the modes $\ell =2$ of oscillation. The first relevant time scale is the oscillator reduced period $t_{2}=1/\omega _2$ of the second mode $\ell =2$ of oscillation, which has been proven to dominate the deformation (Risso & Fabre Reference Risso and Fabre1998). 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 Reference Risso and Fabre1998).

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 ${{\rm \pi} }_{\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, i.e.

(2.12)$$\begin{gather} {\tau}_{\ell,m} = \frac{1}{\langle u_{\ell,m}^2 \rangle}\int_{-\infty}^{+\infty} \textrm{d}s \,\langle u_{\ell,m}(t) \, u_{\ell,m}(t+s) \rangle. \end{gather}$$

(2.13)$$\begin{gather}{\tilde\tau}_{\ell,m} = \frac{1}{\langle {\rm \pi}_{\ell,m}^2 \rangle} \int_{-\infty}^{+\infty} \textrm{d}s \,\langle {\rm \pi}_{\ell,m}(t) \, {\rm \pi}_{\ell,m}(t+s) \rangle. \end{gather}$$

The estimate of ${\tau }_{\ell ,m}$ presents some difficulties, as there is no clear consensus on the time scales in turbulence. The velocity fluctuations we consider here are neither associated with the Eulerian frame of reference nor a local Lagrangian framework. To provide an estimate of the time scale $\tau _{2,m}$ of the dominant mode of deformation, we consider the time scales $t_E$ and $t_L$, respectively, of the Eulerian velocity increments and the Lagrangian velocity increments at the bubble scale. The ratio $t_E/t_L$ scales as $(Re_\lambda )^{1/2}$ from the random sweeping hypothesis (Tennekes Reference Tennekes1975), such that $t_E$ and $t_L$ can respectively be considered as lower and upper bounds of the correlation time $\tau _{2,m}$. The frozen Taylor hypothesis will henceforth be considered as valid for $t< t_E$, and we estimate $\tilde \tau _{2,m}$ from the eddy turnover time at the bubble scale $t_c = \epsilon ^{-1/3} d_0^{2/3}$. The correlation times $\tau _{\ell ,m}$ and $\tilde \tau _{\ell ,m}$ of the higher-order modes ($\ell > 2$) are expected to scale as $\ell ^{-2/3}$ whereas the capillary time 1$/\omega _\ell$ scales as $\ell ^{-3/2}$. The criterion $\omega _2 \tau _2 < 1$ therefore implies that $\omega _\ell \tau _\ell < 1$ for higher values of $\ell$.

Figure 2 shows the ratio of the typical capillary time scales ($t_0$ and $t_2$) with the eddy turnover time scale $t_c$ as a function of Weber number. Below the critical value ${We}<{We}_c$ with ${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 Weber number, namely ${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 (Risso & Fabre Reference Risso and Fabre1998; Martinez-Bazan et al. Reference Martinez-Bazan, Montanes and Lasheras1999; Vejražka et al. Reference Vejražka, Zedníková and Stanovskỳ2018; Ruth et al. Reference Ruth, Mostert, Perrard and Deike2019). 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. (Reference Masuk, Salibindla, Tan and Ni2019). 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.

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 turnover 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 _\ell )}$ and the time associated to mode 2 is $t_2=1/\omega _2 = 1/(2 \sqrt {3} \omega _0)$. Symbols correspond to the direct numerical simulation cases in § 3.

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. Estimates of the integral length scale $L$ were provided by most of the authors, we infer the values for Martinez-Bazan et al. (Reference Martinez-Bazan, Montanes and Lasheras1999) using the usual statistical properties of a turbulent jet (Hussein, Capp & George Reference Hussein, Capp and George1994). The oscillation time is given by $t_2={1}/{2 \sqrt {3}\sqrt {\gamma /(R_0^3\rho _\ell )}}$ and the eddy turnover time at the bubble scale is given by $t_c = \epsilon ^{-1/3} d_0^{2/3}$.

We note that the compressible mode of deformation for a bubble or diameter $2R$ in a infinite domain of water oscillates at the Minnaert frequency $f_M = 1/(2{\rm \pi} R) (3 \varGamma p_0/\rho _\ell )^{1/2}$ (Minnaert Reference Minnaert1933), where $\varGamma$ is the polytropic coefficient and $p_0$ the bubble pressure. When considering millimetric air bubbles in water, the Minnaert frequency $f_M$ ranges from $[0.5,10]$ kHz for sizes in the range $[10,0.5]$ mm (Deane & Stokes Reference Deane and Stokes2002). The acoustic frequency $f_M$ can be compared with the Lamb oscillation frequency $1/t_2=\omega _2/(2 {\rm \pi})$ (Lamb Reference Lamb1995) which is driven by surface tension forces and characteristic of the large deformation and break-up dynamics. For an air bubble in water of millimetric size, this second mode of oscillation $1/t_2$ ranges between [0.025,1] kHz, so that the ratio $f_M t_2 > 10$ and air compressibility can be neglected when analysing large deformations leading to break-up.

2.5. 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 ${{\rm \pi} }_{\ell ,m}$ as independent of time, and neglect the surface tension term. From (2.8) and its initial conditions, we get by double integration over time,

(2.14)\begin{equation} a_{\ell,m} = u_{\ell,m} t + \frac{C(\ell)}{2 R_0} {\rm \pi}_{\ell,m} t^2. \end{equation}

The relative influence of each term will depend on the characteristic time scale $2 R_0 u_{\ell ,m}/{{\rm \pi} }_{\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}$ reduce to a simple advection by the background flow. Performing an ensemble-average operation yields

(2.15)\begin{equation} \langle a_{\ell,m}^2 \rangle = \langle u_{\ell,m}^2 \rangle t^2. \end{equation}

At short time ($\omega _2 t \ll 1$, $t \ll t_c$), the evolution of $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 turnover time at the scale of the neck (Ruth et al. Reference Ruth, Mostert, Perrard and Deike2019). In both cases, the influence of the turbulent background flow reduces to the set of 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 figure 2, or by viscous effects which are not investigated here. For ${We} \gg 1$, when $t_c < t_2$, the linear regime occurs until $t \approx R_0/u_{\ell ,m}' \approx t_c$, which corresponds to $\zeta _\varOmega /R_0 \approx 1$, a deformation amplitude sufficient to trigger bubble breaking. Here $t_c$ then coincides with the bubble lifetime, and no later-time evolution is expected.

2.6. Forced oscillator regime

We consider ${We} < 100$, in which the linear growth is first limited by surface tension ($\omega _2 t_c >1$). For $t < t_c$, we still consider that $u_{\ell ,m}$ and ${{\rm \pi} }_{\ell ,m}$ are independent of time. We have

(2.16)$$\begin{gather} \ddot a_{\ell,m} + \omega_\ell^2 a_{\ell,m} = \dot u_{\ell,m} + \frac{C(\ell)}{R_0} {\rm \pi}_{\ell,m}, \end{gather}$$

(2.17)$$\begin{gather}\dot a_{\ell,m}(0) = u_{\ell,m}, \end{gather}$$

which admits a solution of the form

(2.18)\begin{equation} a_{\ell,m} = \frac{u_{\ell,m}}{\omega_\ell} \sin(\omega_\ell t) + \frac{C(\ell) {\rm \pi}_{\ell,m}}{R_0 \omega_\ell^2} (1 - \cos(\omega_\ell t)). \end{equation}

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

(2.19)\begin{align} \langle a_{\ell,m}^2 \rangle =\langle u_{\ell,m}^2 \rangle t^2 + \frac{C(\ell) \langle {\rm \pi}_{\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{C(\ell)^2 \langle {\rm \pi}_{\ell,m}^2 \rangle}{R_0^2} \frac{t^4}{4} + O(t^5). \end{align}

However, since the terms in $\langle {\rm \pi}_{\ell ,m} u_{\ell ,m} \rangle /R_0$ and $\langle {\rm \pi}_{\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 nonlinear 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

(2.20)\begin{equation} \langle a_{\ell,m}^2 \rangle = \langle u_{\ell,m}^2 \rangle t^2 - \tfrac{1}{3} \omega_\ell^2 \langle u_{\ell,m}^2 \rangle t^4 + O(t^5), \end{equation}

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

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 (http://basilisk.fr/) (Popinet Reference Popinet2009, Reference Popinet2018). We use a spatial adaptive octree grid allowing us to save computational time while resolving the different length scales of the problem and a momentum conserving scheme. The interface is reconstructed by a sharp geometric volume of fluid method (Popinet Reference Popinet2009, Reference Popinet2018). The solver has been extensively described in recent publications (Popinet Reference Popinet2015; Fuster & Popinet Reference Fuster and Popinet2018; Popinet Reference Popinet2018; van Hooft et al. Reference van Hooft, Popinet, van Heerwaarden, van der Linden, de Roode and van de Wiel2018; Mostert & Deike Reference Mostert and Deike2020), and its accuracy has been largely validated on complex multiphase flow, including bubble dynamics (Cano-Lozano et al. Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016), bubble bursting (Deike et al. Reference Deike, Ghabache, Liger-Belair, Das, Zaleski, Popinet and Seon2018; Lai, Eggers & Deike Reference Lai, Eggers and Deike2018; Berny et al. Reference Berny, Deike, Séon and Popinet2020) and wave breaking (Deike, Melville & Popinet Reference Deike, Melville and Popinet2016; Mostert & Deike Reference Mostert and Deike2020; Mostert, Popinet & Deike Reference Mostert, Popinet and Deike2021). We do not consider the effect of gravity in this work. The turbulent two-phase simulations of bubble deformation in turbulence are 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 {\boldsymbol u}$, following the approach of Rosales & Meneveau (Reference Rosales and Meneveau2005), previously implemented and provided as an example in the Basilisk library (http://basilisk.fr/src/examples/isotropic.c). Rosales & Meneveau (Reference Rosales and Meneveau2005) 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 an approach has also been used by Loisy & Naso (Reference Loisy and Naso2017) 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 a 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 an 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 assess statistical properties of the bubble deformation dynamics.

Figure 3(a) shows the time evolution of the kinetic energy in the precursor simulation, $K = ({1}/{V})\iiint \frac {1}{2}\rho _\ell u(\boldsymbol {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 turnover times at the integral scale $\tau = u'^2 / \epsilon$, where $\epsilon$ is the asymptotic dissipation rate and $u'$ the asymptotic root mean squared velocity. The associated Taylor Reynolds number is $Re_{{\lambda }} = ({2 K}/{3 \mu })\sqrt {{15 \mu }/{\epsilon }}\approx 38$, which is a typical value for current two-phase simulations of turbulent flow (Loisy & Naso Reference Loisy and Naso2017; Elghobashi Reference Elghobashi2019). 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{_{\lambda }}=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/4 D_{NN}$. Kolmogorov theory is superimposed in the 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 3(b) 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

(3.1)$$\begin{gather} D_{LL}(d) = \tfrac{1}{3} \sum_i \langle \left(u_i(\boldsymbol r,t) - u_i(\boldsymbol r+ d \hat{{\boldsymbol r}}_{{\boldsymbol i}},t) \right)^2 \rangle, \end{gather}$$

(3.2)$$\begin{gather}D_{NN}(d) = \tfrac{1}{6} \sum_{i \neq j} \langle \left(u_i(\boldsymbol r,t) - u_i(\boldsymbol r+ d \hat{{\boldsymbol r}}_{{\boldsymbol j}},t) \right)^2 \rangle, \end{gather}$$

for homogeneous and isotropic flows, with $\hat {{\boldsymbol r}}_{{\boldsymbol i}}$ the unit vector along the $i$ direction. The transverse structure function is compensated by its scaling for a homogeneous and isotropic flow $(d\epsilon )^{2/3}$, and we indeed observe a plateau value close to $C=2$ (Pope Reference Pope2000). The relation $D_{LL} = 3/4 D_{NN}$ is also verified by representing the compensated longitudinal structure function $4/3 D_{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

(3.3)$$\begin{gather} \delta u_S(R,\theta,\varphi) = {\tilde {{\boldsymbol u}}}({\boldsymbol R},t) \cdot \hat{{\boldsymbol r}}, \end{gather}$$

(3.4)$$\begin{gather}\tilde {{\boldsymbol u}}({\boldsymbol R}) = {\boldsymbol u} - \frac{1}{4{\rm \pi}} \iint \textrm{d}\varOmega {\boldsymbol u}({\boldsymbol R},\theta,\varphi,t). \end{gather}$$

From homogeneity and isotropy, the statistical properties of $\delta u_S(R)$ shall only depend on $r$. Figure 3(c) 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'_S(R) =\sqrt { \langle (\delta u_S(R))^2 \rangle }$ compensated by $(\epsilon d)^{2/3}$ is represented in red in figure 3(b). We observe $\delta u'_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 § 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 $\eta _K = (\nu ^3/\epsilon )^{1/4}$ is the Kolmogorov length scale, $\lambda$ is the Taylor microscale $\lambda = \sqrt {15 \nu u'^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. Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016).

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 a Weber number immediately smaller than ${We} = {We}_c$. To determine the critical Weber number, we perform an 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 break-up at ${We} = 1.5$ while running the simulations up to $20t_c$ for more than 10 runs with different initial conditions. For ${We}=3$, we observe break-ups a little over 50 % of the time by $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 ${We}=6$, and above, all bubbles break within a few $t_c$. As such, we consider ${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.

Note that we have verified that the order of magnitude of the pressure fluctuations at the bubble surface are too small to induce compressibility effects in the bubble's dynamical response, hence validating the incompressible framework to study millimetric bubbles in water turbulence.

The Weber number, (1.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 numbers, $1.5 < We < 45$. The critical Weber number is found at ${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 (${We} \leq 1.5$), no break-up occurs, while for large Weber number, (${We} \geq 3$), we will analyse the dynamics before the first break-up. We perform an ensemble of simulations of 10 to 15 runs for each Weber 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 turnover time at the bubble scale $t_c$. The bubble is initially spherical and quickly deforms for $t<0.5 t_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.

Figure 4.

Snapshots of a bubble for one run at ${We} = 3$, with the interface in white and each background plane showing one component of the velocity. The bubble is injected in the centre 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.

3.3. Computing the bubble deformation through a Voronoi decomposition

The interface sampling points $\boldsymbol 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 Reference Voronoi1908). The bubble surface is partitioned into regions close to each of the points.

For each region located around the point ${\boldsymbol r}_i$, we associate a region area $A_i$ and a corresponding solid angle $\varOmega _i = A_i/r_i^2$. We evaluate the surface integral of any function $f$ by

(3.5)\begin{equation} \frac{1}{4{\rm \pi}} \iint \textrm{d}\varOmega f({\boldsymbol r}) = \frac{1}{4 {\rm \pi}}\sum_i f({{\boldsymbol r}}_i) \varOmega_i, \end{equation}

where $\sum _i \varOmega _i = 4 {\rm \pi}$. The computation of the region areas $A_i$ and the corresponding solid angles $\varOmega _i$ is illustrated in figure 5. The $\varOmega _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. Reference Caroli, de Castro, Loriot, Rouiller, Teillaud and Wormser2009) (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 $\varOmega _i$ are identical. An example of our Voronoi decomposition is shown in figure 5(b) on the unit sphere, and in figure 5(c) after projection on the initial bubble shape. The numerical error on the total solid angle $\sum _i \varOmega _i$ is less than 0.1 %.

Figure 5.

Illustration of the interface tesselation using a spherical Voronoi diagram. (a) Bubble deformed by the turbulent background (one case at ${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 data point is given by the associated polygon area.

To analyse the interface dynamics using (2.8), 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 ${\boldsymbol r}_i$ weighted by their local solid angle $\varOmega _i$,

(3.6)\begin{equation} a_{\ell,m} = \sum_i \varOmega_i \zeta({\boldsymbol r}_i) Y_{\ell}^{m} (\theta_i,\varphi_i), \end{equation}

where $(r_i, \theta _i,\varphi _i)$ are the spherical coordinates of ${\boldsymbol r}_i$. We have first checked the orthonormal property of the spherical harmonics $({1}/{4{\rm \pi} }) \iint \textrm {d}\varOmega \,Y_{\ell }^{m} \cdot Y_{(\ell ',m')} = \delta _{\ell ,\ell '} \delta _{m,m'}$ 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, since the amplitude decreases 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 centre 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 ${\boldsymbol r}_c$ of the bubble frame of reference, for which $a_{1,m} = 0$ for $m \in \{-1,0,1\}$, as defined in § 2. The centre position ${\boldsymbol r}_c$ is computed recursively as follows. 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 centre that minimizes each coefficient $a_{1,m}$. For bubble interfaces with a single-valued radial distance $r_i(\theta ,\varphi )$, the gradient descent indeed converges to a centre position ${\boldsymbol 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

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 the same $\ell$ and different $m$ values are associated to the same oscillation frequency (see (2.8)); therefore, we introduce the global coefficient $a_\ell$ describing the energy contained in each mode $\ell$, defined by

(3.7)\begin{equation} a_\ell^2 = \sum_{m={-}\ell}^{\ell} a_{\ell,m}^2, \end{equation}

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 _\varOmega$: $\zeta _\varOmega ^2(t) = \sum _{\ell =2}^{+\infty } \sum _{m=-\ell }^\ell a_{\ell ,m}^2(t)$. Figure 6(a) shows the time evolution of $\zeta _\varOmega /R_0$ together with the total amplitude of 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 _\varOmega /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 Reference Risso and Fabre1998; Ravelet et al. Reference Ravelet, Colin and Risso2011). In practice, a sum over the first three modes $\ell = 2,3,4$ estimates the amplitude deformation $\zeta _\varOmega$ within less than 2 % error as long as $\zeta _\varOmega /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 6.

(a) Total deformation $\zeta _\varOmega /R_0$ as a function of time, together with the root mean square (r.m.s.) 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 =4$ for one run at ${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 {\rm \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.

Figure 6(b) 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 centre 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 6(c) shows the temporal evolution of each individual harmonic coefficient $a_{2,m}$ for different $m$ components. All five 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 numbers, ${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 § 2. For all ${We}$, the short time deformation under the same turbulence conditions appears 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 times. 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$.

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.

A quantitative description of the Weber number influence for one particular run is given in figure 8. We compute the total deformation $\zeta _\varOmega$ 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 _\varOmega /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 (2.14). As time passes, the curves diverge from each other, and the lower Weber number curves exhibit earlier sublinear 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 numbers (here ${We}=15$ and 45), the bubble experiences break-up during the first eddy turnover time, when surface deformation reaches $\zeta _\varOmega /R_0 \approx 0.74$. Note that the deformation at break-up is estimated from a fit of the bubble interface by an ellipsoidal shape at the break-up time, as the Voronoi decomposition fails at large deformation, significantly before break-up.

Figure 8.

(a) Root mean square surface deformation $\zeta _\varOmega$ 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 Weber numbers. 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

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 _\varOmega$ and the first harmonics coefficients $a_2, a_3$ and $a_4$.

Figure 9(a) shows the temporal evolution of each individual realisation of $\zeta _\varOmega /R_0$ performed at ${We} = 3$, as a function of the dimensionless time $t/t_c$. The ensemble average $\langle \zeta _\varOmega \rangle /R_0$ is superimposed in black, and follows a linear increase, as predicted by (2.15). The mode 2 amplitude $a_2/R_0$ is superimposed in figure 9(a) (${\square }$, red), confirming in a statistical sense the predominance of the mode 2 oscillations. The modes $a_2/R_0$ and $a_3/R_0$ follow the same trend, as shown in figures 9(b) and 9(c), each coloured curve being an individual realisation. Again, at short time the linear regime predicted by (2.15) is observed.

Figure 9.

(a) Individual realisations of $\zeta _\varOmega$ at ${We} = 3$, as a function of the dimensionless time $t/t_c$ in transparent colours. Ensemble average $\langle \zeta _\varOmega \rangle$ of the surface deformation is superimposed (black full symbols), 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 ${We}=3$ (individual realisations in transparent colours). Ensemble average $\langle a_2\rangle$ is superimposed (full red symbol). Fit by $A (1-\textrm {e}^{-t/t_{sat}})$ of the ensemble average $a_2$ in shown by the dashed lines. (c) Coefficient $a_3$ as a function of dimensionless time $t/t_c$ for ${We}=3$ and the equivalent exponential fit (individual realisations are transparent colours and the green full triangle is ensemble average $\langle a_3\rangle$).

We observe a transition to a sublinear regime for $\langle \zeta _\varOmega \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 § 2 (2.20), i.e. $\propto t \sqrt {(1 - (\kappa t)^2)}$, with $\kappa = \omega _2/\sqrt {3}$ for $\langle \zeta _\varOmega \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 _\varOmega \rangle$, $\langle a_2 \rangle$ and $\langle a_3 \rangle$ by the empirical expression $a_\ell ^\infty (1 - \textrm {e}^{-t/t_{sat}} )$ 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 summarizes the surface deformation growth for increasing Weber number. Figure 10(a) shows the evolution of $\langle \zeta _\varOmega \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 a universal regime, with a linear increase in time of the deformations independent of the Weber number. 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 10(b) as a function of the Weber number for $\langle \zeta _\varOmega \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 (2.8) and (2.15).

Figure 10.

(a) Ensemble-averaged deformation $\langle \zeta _\varOmega \rangle$ as a function of time for increasing Weber number. For the lower ${We}$, a fit by $A (1-\textrm {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 _\varOmega \rangle$, $\langle a_2 \rangle$, $\langle a_3 \rangle$, $\langle a_4 \rangle$. The slope at the origin for $\langle \zeta _\varOmega \rangle$ is independent of $We$. (c) Amplitude of saturation of the global deformation and modes as a function of $We$.

The departure from linear growth that occurs at a time increasing with Weber number is compatible with the prediction given 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 numbers, we consider the value of the linear fitting model, evaluated at the time of the first break-up. The saturation time for ${We} < 10$ and the earliest break-up value for ${We} > 10$ is shown in figure 10(c). We observe a continuous growth with Weber number, premising the likelihood of bubble break-up. At higher Weber numbers, we also remark that the relative amplitudes of modes 3 and 4 are increasing for both the growth velocity $s_\ell$ and the saturation $a_\ell ^\infty$. At high Weber numbers, the higher-order modes of deformation in the global dynamics become more important and are likely to influence the break-up geometry that the bubble will experience.

Note that treating the saturation at an intermediate Weber number would require considering the role of nonlinear and viscous effects. Note also that for $\zeta _\varOmega /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 _\varOmega /R_0\approx 0.74$, which is compatible with the observations on areas of deformation at break-up made by Risso & Fabre (Reference Risso and Fabre1998).

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

From the theory developed in § 2, the prefactor 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}$ is linked to their counterpart in the turbulent fluctuations, ${u}_{\ell ,m}$ and ${{\rm \pi} }_{\ell ,m}$ in (2.8), (2.15) and (2.20). 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 direct numerical simulations, 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 11(a), and in table 2 as red stars for $\ell = 2, 3$ and $4$ and correspond to the intensity of velocity fluctuations at the bubble scale for specific modes $\ell$. The coefficients $s_\ell t_c/R_0$ obtained from the processing of surface deformations are superimposed in black circles 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 § 2, (2.15). For $\ell = 4$, we observe a slight difference for $\ell = 4$ with $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.

Figure 11.

(a) Ensemble average of the spherical harmonics growth rate $\langle s_\ell \rangle R_0/t_c$ for $\ell =2, 3$ and 4, together with the corresponding spherical harmonics flow coefficient $\langle u_\ell \rangle t_c/R_0$ of the turbulent velocity fluctuations at the bubble scale, shown as red stars. Black circles are $\langle s_\ell \rangle t_c/R_0$ averaged for all ${We}$, while the coloured squares correspond to the Weber number scale (see colourbar). The theoretical prediction $a_\ell = u_\ell t$ from § 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 $\varGamma$-distribution of expression $f(x,k,p)=x^{k-1} \,\textrm {e}^{-x/p}/\varGamma (k) p^k$ is superimposed (red dashed line) with $k=7.5$ and $p=0.135$.

Table 2.

Values of the coefficients shown in figures 10 and 11. Ensemble average 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_\ell \rangle t_c/R_0$ of the turbulent velocity fluctuations at the bubble scale is also provided, computed for $Re_\lambda =38$.

Beyond the equality of the ensemble-average values, the prediction made in § 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.

Figure 11(b) shows the probability density function (p.d.f.) 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 p.d.f. 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 p.d.f. 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 smoothes 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 p.d.f. by a $\varGamma$-distribution of expression $f(x,k,p)=x^{k-1} \,\textrm {e}^{-x/p}/\varGamma (k) p^k$ gives $k=7.5$ and $p=0.135$.