# Spectral analysis of dispersed multiphase flows in the presence of fluid interfaces

Gabriel Ramirez, Alan Burlot, Rémi Zamansky, Guillaume Bois, Frédéric Risso

► To cite this version:

Gabriel Ramirez, Alan Burlot, Rémi Zamansky, Guillaume Bois, Frédéric Risso. Spectral analysis of dispersed multiphase flows in the presence of fluid interfaces. International Journal of Multiphase Flow, 2023, 177, pp.104860. 10.1016/j.ijmultiphaseflow.2023.104860. hal-04295313v4

HAL Id: hal-04295313

https://hal.science/hal-04295313v4

Submitted on 14 Oct 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Spectral analysis of dispersed multiphase flows in the presence of fluid interfaces

Gabriel Ramirez<sup>a,b</sup>, Alan Burlot<sup>b</sup>, Rémi Zamansky<sup>a</sup>, Guillaume Bois<sup>b</sup>, Frédéric Risso<sup>a,*</sup>

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

## ARTICLE INFO

For the Special 50th Anniversary Issue of IJMF

**Keywords:**

- Dispersed flows

- Bubbles

- Droplets

- Turbulence

- Body-induced-agitation

- Spectral analysis

## ABSTRACT

Spectral analysis of dispersed two-phase flows is highly desirable to reveal the interplay of the various flow scales, much larger or much smaller than the size of the dispersed bodies. This is a challenging task as the matching conditions at the body interfaces generate singularities in the fields describing the two-phase mixture. The nature of these singularities and their consequences on the spectra are theoretically analyzed for bubble or droplet flows. Results of direct numerical simulations are reported and spatial spectra of the mixture velocity, the flow forces and their power are examined. The regular part of the spectral densities of energy production, dissipation and transfers between scales are separated from their singular part. The resulting spectral energy balance, free of the footprint of the singularities, is found in agreement with coarse-grained simulations where the interfaces are filtered out before solving the Navier-Stokes equations. These results pave the way for the spectral analysis of more complex turbulent dispersed flows.

## 1. Introduction

The dynamics of dispersed multiphase flows are controlled by the interactions between a population of bodies – either droplets, bubbles or solid particles – and a fluid phase in which they are randomly distributed. The bodies have a great impact on flow fluctuations either by modulating a pre-existing shear-induced turbulence (Balachandar and Eaton, 2010; Brandt and Coletti, 2022) or by directly inducing agitation through their motion relative to the fluid (Risso, 2018). Spectral analysis is a powerful tool to understand turbulence, since two-point correlations give access to the flow structure. Applied to multiphase flows, it should reveal the interplay between scales and lead to a better modeling. In particular, our objective is to determine the spectral density of the power of all the forces that contribute to the energy balance, distinguishing between rate of production, dissipation and transfer between scales. However, the question of how to achieve this in the presence of numerous sharp interfaces between the carrier and dispersed phases is not trivial.

Since the pioneering work of Lance and Bataille (1991), many experimental works have reported one-dimensional spectra of the liquid phase velocity in bubbly flows. Different methods were used to deal with the interruptions by the bubbles: removing the parts of the signal belonging to the gas phase and filling the gaps by a smooth function (Lance and Bataille, 1991), retaining only parts of the signals between two bubbles which are not interrupted (Martínez Mercado

et al., 2010; Mendez-Diaz et al., 2013; Roghair et al., 2011; Prakash et al., 2016; Alméras et al., 2017), measuring the flow just behind a rising bubble swarm (Riboux et al., 2010). A few numerical studies have proposed spectral analysis of the velocity field of the carrying phase in order to make direct comparisons with experiments (Roghair et al., 2011). However, spectra obtained from Direct Numerical Simulations (DNS) of dispersed two-phase flows are generally calculated by considering the entire flow field without distinguishing between carrying and dispersed phases (Tryggvason et al., 2002; Lucci et al., 2010; Dodd and Ferrante, 2016; Pandey et al., 2020; Innocenti et al., 2021; Pandey et al., 2022; Cialesi-Esposito et al., 2022).

The question of whether the spectral analysis should take into account the entire flow mixture or only the carrying phase deserves to be discussed. In this work, we focus on numerical simulations where the flow is known everywhere and consider the case of a statistically homogeneous dispersed flow in a large periodic domain, which is well suited to a spectral description in the k-wavenumber domain.

We examine first the option of building spectra of the continuous phase only, and immediately set aside the method of filling the region occupied by the dispersed phase by an arbitrary smooth field. Considering a region of the flow that is never crossed by an interface ensures that the field of any physical quantity is smooth. However, the surface of the bodies are now boundaries from which momentum and

<footnote><sup>a</sup> Institut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse and CNRS, Toulouse, France</footnote>

<footnote><sup>b</sup> Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides, 91191, Gif-sur-Yvette, France</footnote>

<footnote>* Corresponding author.
E-mail address: frisso@imft.fr (F. Risso).</footnote>
energy are supplied to the system under study. All the local statistical quantities, such as velocity variance, energy dissipation rate, pressure gradient or inertial forces, depend on the distance to those boundaries. In such a heterogeneous flow, an average spatial spectrum $S(k)$ cannot represent a meaningful spectral density of any of these fields. This would make no more sense than constructing a spatial spectrum from the velocity along a line perpendicular to a wall in a channel flow. An attempt to overcome this limitation has been done by Freund and Ferrante (2019) who analyzed DNS results of a droplet laden flow. They used wavelets, a tool developed to deal with inhomogeneous fields by allowing the construction of spectra, $S(x, k)$, that depend also on the location $x$. For each wavenumber, the domain was decomposed into three regions in contact but without intersections: an intermediary region $I$ crossed by the interfaces, a region $C$ fully located in the continuous phase and a region $D$ fully located in the droplets. Vector $x$ was only used to determine to which region belongs a given point and three spectra, $S_I(k)$, $S_C(k)$ and $S_D(k)$, one for each region, were calculated. Since the spectra depend only on $k$, they still face the inconsistency of describing highly heterogeneous fields. A more sophisticated use of wavelets preserving the double dependency in $x$ and $k$ could lead to a relevant description, but has yet to be developed for the study of dispersed multiphase flows. In this work, we stay with the conclusion that, in the context of a homogenized description based on spectra depending only on the wavenumber, it is not relevant to consider the flow field of the carrying phase only.

The remaining option is to consider the flow of the entire mixture. Since it is statistically homogeneous, the use of the Fourier transform is appropriate. However, the fields of the physical quantities under investigation experience singularities at the interfaces between the phases. These singularities are due to the presence of jumps across the interface of some physical properties, such as density and viscosity, and to surface tension. They are of several kinds, such as the discontinuities of the velocity derivatives or the Dirac delta function that describes the pressure gradient. These singularities may have a strong footprint on the spectra. In particular, they can generate oscillations, which were very early reported on the velocity spectra of dispersed two phase flows with large bubbles (Tryggvason et al., 2002) or solid particles (Lucci et al., 2010). These oscillations are associated with the Gibbs phenomenon, familiar in signal processing, and well exemplified by the sine cardinal function, which is the Fourier transform of a gate. Their signature also contains a power-law decay at large wavenumbers, which must not be confused with a dynamical phenomenon (Risso, 2011). The spectra of any flow-field quantity contains a regular part that describes the smooth variations in the bulk phases and a singular part that accounts for interfacial singularities.

In what follows, we focus on droplets or bubbles in a fluid, which means that the dynamics of the two phases are described by the Navier-Stokes equations and the interface by an interfacial tension. A similar analysis could be carried out for solid particles in a fluid, but it is out of the scope of the present work. Our objective is to examine the consequences of the presence of interfacial singularities on spectra. The singular part of a spectrum is sometimes qualified as spurious, since it complicates the physical interpretation. However, in the context of the analysis of the mixture fields, it is a part of the mathematical solution of the physical problem. The possibility of separating the regular and singular part is a central question that motivates this work.

This paper is organized as follows. Section 2 presents the dynamics equations of the mixture flow and examines the nature of the different interfacial singularities. Section 3 analyzes the spectral signatures of basic singularities caused by droplet or bubble interfaces. Section 4 reports results of direct numerical simulations and discusses the spectral density of the velocity, of the terms of the Navier-Stokes equations and of the terms of the energy balance. Section 5 compares with results of coarse-grained simulations, where interfacial singularities are filtered out before solving the Navier-Stokes equations. Section 6 summarizes the main findings and concludes.

<figure><img src="image_4.png" /><figcaption>Fig. 1. Local coordinates nearby an interface.</figcaption></figure>

## 2. Equations of motion of a two-phase mixture and field singularities

In this section, we examine the nature of the singularities encountered by the fields describing the flow, which are due to the presence of fluid interfaces.

We decided to analyze the spectral distribution of the entire mixture. For this reason, we write the equations of motion of the fluids as a single set of partial differential equations for the velocity $u$ and pressure $P$, which is valid throughout both phases. This means that the matching of the dynamic conditions between the two phases at the interfaces have to be inserted within the Navier-Stokes equations for the mixture. Because that matching involves jumps of stresses through the interfaces, the use of generalized functions (Gelfand and Shilov, 1964), known as distributions, is required. Note that this approach is known as one-fluid sharp interface methods in the context of numerical simulations of two-phase flow (Kataoka, 1986; Lalanne et al., 2015).

We consider a population of fluid particles, either drops or bubbles, dispersed in another fluid. The two fluids are immiscible and form a two-phase mixture, without phase change. The interfaces between the phases are massless, of zero thickness and entirely characterized by a surface tension $\sigma$. The density and viscosity of the two fluids are different, but constant in each phase. In the following, the carrier phase, dispersed phase and interfaces are distinguished by the subscripts $c$, $d$ and $s$, respectively. We introduce the indicator function, $\chi$, which is equal to one in the dispersed phase and to zero in the carrying phase. The density $\rho$ and the viscosity $\mu$ of the mixture can be written:

$$ \rho = \rho_c + (\rho_d - \rho_c)\chi, \qquad (1) $$

$$ \mu = \mu_c + (\mu_d - \mu_c)\chi. \qquad (2) $$

In the context of continuum mechanics, in the absence of phase change, the velocity field is continuous. The mass conservation is

$$ \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = D_t \rho + \rho \nabla \cdot \mathbf{u} = 0. \qquad (3) $$

Because $\rho$ is constant in each phase and no fluid particle crosses interfaces, the material derivative ($D_t = \partial_t + \mathbf{u} \cdot \nabla$) of the density is zero and the velocity field is divergence-free: $\nabla \cdot \mathbf{u} = 0$. This result has further consequences for the differentiability of $\mathbf{u}$. Let us consider a portion of the interface and name $x_{\parallel 1}$ and $x_{\parallel 2}$ the local coordinates along the interface, and $x_{\perp}$ the coordinate orthogonal to it (Fig. 1). Since $\mathbf{u}$ is continuous across the interface, the four strain-rate components $\partial_{\parallel \xi} u_{\eta}$ (with $\xi$ and $\eta$ equal to 1 or 2) are equal on both sides of the interface and thus have no discontinuity. Then, considering that $\nabla \cdot \mathbf{u} = 0$, we get that the normal strain rate $\partial_{\perp} u_{\perp} = -\partial_{\parallel 1} u_{\parallel 1} - \partial_{\parallel 2} u_{\parallel 2}$ is continuous too. Thus, the only derivatives of $\mathbf{u}$ that can present a discontinuity are the shear-rate components $\partial_{\perp} u_{\parallel \xi}$, the regularity of which is prescribed by the shear-stress matching across the interface and will be discussed later.
<table><caption>**Table 1**
Interfacial discontinuities (Eq. (7)) involved in each terms of the momentum equations (Eq. (4)). A check mark means that the corresponding discontinuity-type is present.</caption><thead><tr><th></th><th>Regular part</th><th>p ≥ 1<br>(continuous)</th><th>p = 0<br>Heaviside type</th><th>p = -1<br>Dirac type</th></tr></thead><tbody><tr><td>Buoyancy: ρg</td><td>✓</td><td></td><td>✓</td><td></td></tr><tr><td>Inertia: -ρD<sub>l</sub>u</td><td>✓</td><td></td><td>✓</td><td></td></tr><tr><td>Interfacial force: F<sub>σ</sub></td><td></td><td></td><td></td><td>✓</td></tr><tr><td>Pressure gradient: -∇P</td><td>✓</td><td>?</td><td>✓</td><td>✓</td></tr><tr><td>Viscous diffusion: ∇·(2μS)</td><td>✓</td><td>?</td><td>✓</td><td>✓<br>if μ<sub>c</sub> ≠ μ<sub>d</sub><br>or ∇<sup>S</sup>σ ≠ 0</td></tr></tbody></table>

The momentum equations are formally written as the classic Navier–Stokes equations with a non-constant viscosity in which surface tension is taken into account by adding an interfacial force:

$$-\rho D_l \mathbf{u} - \nabla P + \rho \mathbf{g} + \nabla \cdot (2\mu \mathbf{S}) + \mathbf{F}_{\sigma} = 0, \quad (4)$$

where $S = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$ is the strain-rate tensor. The interfacial force is written as

$$\mathbf{F}^{\sigma} = -\sigma \nabla^{S} \mathbf{n} \nabla \chi + \nabla^{S} \sigma \int_{\mathbf{x}_{s}} \delta_{3D}(\|\mathbf{x} - \mathbf{x}_{s}\|) dS, \quad (5)$$

where $\mathbf{n}$ is unit normal to the interface pointing outward from the dispersed phase, $\nabla^S$ is the two-dimensional gradient operator on the interface, $\mathbf{x}_s$ is any point on the interface $S$ and $\delta_{3D}(\|\mathbf{x} - \mathbf{x}_s\|) = \delta(x - x_s)\delta(y - y_s)\delta(z - z_s)$ with $\delta$ the Dirac delta function. The surface divergence $\nabla^S \mathbf{n}$ of the normal vector is twice the mean curvature $\kappa$ of the interface and the gradient of the indicator function can also be written in terms of the Dirac delta function as

$$\nabla \chi = -\mathbf{n} \int_{\mathbf{x}_s} \delta_{3D}(\|\mathbf{x} - \mathbf{x}_s\|) dS. \quad (6)$$

The first term of Eq. (5) corresponds to the interfacial jump of the normal stress due to the capillary pressure: $-2\sigma\kappa$. The second term accounts for the interfacial jump of the shear stress due to the Marangoni effect: $\partial_{\|\xi}\sigma$.

Let us examine the nature of the singularities of various fields of interest. A field $V$ can be decomposed in a regular singularity-free part $V_{reg}$ and interfacial jumps $[V_{irreg}^{(p)}]$ of its successive derivatives of order $p$:

$$V(\mathbf{x}) = V_{reg}(\mathbf{x}) + \sum_{p \ge P_{min}}^{P_{max}} [V_{irreg}^{(p)}(\mathbf{x} = \mathbf{x}_s)], \quad \text{with } p_{min} \ge -1. \quad (7)$$

The integer $p_{min}$ is the order of the lowest derivative that is discontinuous, which characterizes the strongest singularity of $V$. The sum may continue to infinity or stop at a finite value $p_{max}$, which then characterizes the less sharp singularity. At this point, it is worth mentioning that the regularity of the solutions of the Navier–Stokes equations is a very hard mathematical issue, which is still an open question even in the case of a single-phase flow (Onsager, 1949; Duchon and Robert, 2000; Dubrulle, 2019). As the aim of the present work is to examine the practical consequences of the presence of fluid interfaces on the spectra of the fields, we shall concentrate on the lower values of $p$, which correspond to the more severe singularities with the strongest consequences. Note that $p_{min} = 1$ corresponds to a field that is continuous whereas its derivative experiences a jump across the interface. Handling values of $p$ that are less than unity means that we consider the word “derivative” in a broader sense. Indeed, we have to deal with the cases where the function itself is discontinuous ($p_{min} = 0$) or is the derivative of a discontinuous function ($p_{min} = -1$). For instance, $\chi$, as well as $\rho$ and $\mu$, is a pure step function of Heaviside type, characterized by $p_{min} = p_{max} = 0$.

Table 1 lists the type of discontinuities which are expected to be involved in the various terms of the momentum equations, according

to the decomposition defined by Eq. (7). For the balance to be satisfied, every singularity appearing in one term of the momentum equations must cancel out with singularities of the same order in the other terms. In other words, the sum of the terms in each column must be zero. We examine now the volume forces of the momentum equations one by one.

The buoyancy force $\rho g$ is the product of the density and the acceleration of gravity; it is thus of pure Heaviside type. Check marks are therefore present only in the columns corresponding to the regular part and the $p = 0$ singularity.

Since a fluid particle never crosses an interface and the velocity is continuous, its acceleration $D_l u$ is expected to be regular. However, because it involves the density, the nature of the singularity of the inertial force $\rho D_l u$ is the same as that of the buoyancy force.

Eqs. (5)–(6) show that $F^\sigma$ is of pure Dirac type and associated with $p_{min} = p_{max} = -1$. A check mark is therefore present only in the $p = -1$ column. Note that this is the only term without a regular part, as it is not defined outside the interfaces.

The pressure gradient has generally both Heaviside and Dirac singularities. This is easy to show by considering the elementary case of an interface of constant $\sigma$ between two fluids at rest in a constant gravity field: for example a droplet pending at the tip of a capillary tube. Eq. (4) simplifies to

$$\nabla P = -\rho g - 2\sigma\kappa \nabla \chi, \quad (8)$$

where the pressure gradient has to balance both the Heaviside-type buoyancy force and the Dirac-type interfacial force. The first corresponds to the jump in the hydrostatic pressure: $\partial P_d / \partial z - \partial P_c / \partial z = -(\rho_d - \rho_c)g$. The second accounts for the Laplace pressure jump: $P_d - P_c = -2\kappa\sigma$.

As regards the viscous force, $\nabla·(2\mu\mathbf{S})$, the fact that the viscosity, which is discontinuous, appears under a differential operator immediately suggests the presence of a Dirac-type singularity. However, it is interesting to analyze its significance by distinguishing again between the tangential and normal directions, and by considering its two possible origins, namely: a viscosity jump or a Marangoni stress.

As we have shown before, $\partial_\perp u_\perp$ is continuous across the interface. Therefore, the jump of the normal viscous stress through the interface is equal to $2(\mu_d - \mu_c)\partial_\perp u_\perp$. This means that the existence of a viscosity difference between the phases is sufficient to generate a finite jump in the viscous stress, which leads the viscous volume force in the momentum equation to include a Dirac-type contribution. The matching of the normal stresses at interface is written

$$2(\mu_d - \mu_c)\partial_\perp u_\perp + (P_d - P_c) + 2\kappa\sigma = 0, \quad (9)$$

which shows that even in the case of vanishing surface tension, a jump in viscosity is sufficient to generate a pressure jump. On the other hand, when $\nabla^S \sigma \neq 0$, the Marangoni stress is responsible for a jump in the tangential shear stress. Finally, it thus turns out that the viscous volume force contains a Dirac-type singularity ($p_{min} = -1$) if either $\mu_d - \mu_c \neq 0$ or $\nabla^S \sigma \neq 0$. Since the pressure gradient involves a Heaviside-type discontinuity ($p = 0$), we can expect that the viscous force also does so to satisfy the equilibrium of the normal stresses at the interface under flow conditions. Then, as we have no argument for deciding whether or not the second (and higher) derivative of the viscous stress should be the same on both sides of the interface, we have added question marks in the $p \ge 1$ columns of pressure gradient and viscous force.

To end this section, we conclude with the velocity field. We have already shown that all partial derivatives of $u$ are continuous except the shear-rate components $\partial_\perp u_\perp$. Therefore, if $\mu_d - \mu_c = 0$ and $\nabla^S \sigma = 0$, the continuity of the shear stress ensures that of the shear-rate, so all the velocity derivatives are continuous and $u$ is characterized by $p_{min} \ge 2$. Otherwise, in the general case, $p_{min} = 1$.
### 3. Spectral signatures of basic singularities caused by droplet or bubble interfaces

In this section we examine the spectral signature of singularities of various orders $p$ encountered in multiphase flows. For a more general discussion of spectral analysis of signal with singularities, the reader is referred to Broer and Takens (1993), Bacry et al. (1993), Pikovsky et al. (1995). The energy density spectrum $E_V$ of a scalar field $V(x)$ is defined as follows. First, we take the three dimensional Fourier transform of $V(x)$:

$$ \hat{V}(\mathbf{k}) = \int e^{i\mathbf{k} \cdot \mathbf{x}} V(\mathbf{x}) d\mathbf{x}. \qquad (10) $$

Second, we multiply it by its complex conjugate ($*$), sum over shells $\Lambda_{iw}$ of constant $\|\mathbf{k}\|$ and take the average:

$$ E_V(k) = \left\langle \int_{\|\mathbf{w}\|=k} \hat{V}(\mathbf{w}).\hat{V}^*(\mathbf{w})d\Lambda_{iw} \right\rangle, \qquad (11) $$

(depending on whether $V$ is a scalar or a vector, the operation “.” is a multiplication of numbers or a scalar product).

#### 3.1. Analytical results for spherical interfaces

We introduce the indicator function $H_{sphere}$ of a ball (interior of a sphere) as a reference Heaviside-type field, and the indicator function $\delta_{sphere}$ of a sphere (surface of a ball) as a Dirac-type field. These fields are respectively representative of the density and the interfacial force of a two-phase mixture of spherical droplets. Their Fourier transforms are

$$ \hat{H}_{\text{sphere}}(\mathbf{k}) = \frac{4\pi}{\|\mathbf{k}\|^3} (\sin(R\|\mathbf{k}\|) - R\|\mathbf{k}\| \cos(R\|\mathbf{k}\|)) \qquad (12) $$

and

$$ \hat{\delta}_{\text{sphere}}(\mathbf{k}) = \frac{4\pi R}{\|\mathbf{k}\|} \sin(R\|\mathbf{k}\|) \qquad (13) $$

where $R$ is the sphere radius (Gelfand and Shilov, 1964). Their spectra are

$$ E_{H_{\text{sphere}}}(k) = \frac{(4\pi)^2}{k^4} (\sin(Rk) - Rk \cos(Rk))^2, \qquad (14) $$

and

$$ E_{\delta_{\text{sphere}}}(k) = (4\pi R)^2 \sin^2(Rk). \qquad (15) $$

The spectra of the two cases present oscillations of period $k_b = \pi/R$, corresponding to $\sin^2(Rk)$. On the one hand, the oscillations of $E_{\delta_{sphere}}$ keep constant phase and amplitude for all $k$. On the other hand, $E_{H_{sphere}}$ behaves as $R^6 k^2$ at small $k$ and as $k^{-2}\cos^2(Rk)$ at large $k$, which means it cancels out for both $k=0$ and $k=\infty$.

#### 3.2. Numerical results for spherical interfaces

We have numerically computed the spectra of basic singular scalar fields $V_{pmin}$, characterized by a value $p_{min}$ from -1 to 2, defined as follows:

$$ V_{-1}(\mathcal{X}) = \delta(\mathcal{X}), \qquad (16) $$

$$ V_0(\mathcal{X}) = \mathcal{H}(\mathcal{X}), \qquad (17) $$

$$ V_1(\mathcal{X}) = (\mathcal{X}^2 - 1)((J_1/2 + 1)\mathcal{X}^2 - 1)\mathcal{H}(\mathcal{X}), \qquad (18) $$

$$ V_2(\mathcal{X}) = (\mathcal{X}^2 - 1)^2((J_2/8 - 1)\mathcal{X}^2 + 1)\mathcal{H}(\mathcal{X}), \qquad (19) $$

where $H(\xi)$ is the classic Heaviside function (equal to zero for $\xi < 0$ and to one for $\xi > 0$), $\mathcal{X}$ is the normalized distance from the surface of the sphere: $\mathcal{X} = 1 - \|\mathbf{x} - \mathbf{x}_c\|/R$. $V_{-1}$ and $V_0$ correspond to $\delta_{sphere}$ and $H_{sphere}$. Then, it is easy to check that $V_1$ and $V_2$ have the expected level of singularity provided $J_1$ and $J_2$ are finite. Indeed, both are regular inside the sphere and zero outside, while at the interfaces: $V_1 = 0$ and $V'_1$ undergoes a jump $J_1$; $V_2 = V'_2 = 0$ and $V''_2$ undergoes a jump $J_2$.

(In contrast with $\delta_{sphere}$ and $H_{sphere}$, $V_1$ and $V_2$ do not correspond to a pure singularity of order $p=1$ or $2$, which means that $p_{max} \neq p_{min}$.) These four fields have been mapped on a three-dimensional regular grid of spacing $\Delta = R/25$, which implies that numerical interfaces have a non zero thickness. Then, the spectrum of each of them is computed numerically by making use of the discrete Fourier transform.

Fig. 2 shows the spectra of $V_{-1}$, $V_0$, $V_1$ and $V_2$ as functions of the wavenumber $k$ normalized by $k_b$. The numerical spectrum of $\delta_{sphere}$ is in agreement with the analytical result (Eq. (15)), except at large $k$ where the numerical approximation of the delta function becomes rough. However, the numerical spectrum of $H_{sphere}$ remains accurate beyond $k = 10k_b$, so we can trust as well the numerical spectra of $V_1$ and $V_2$ in the considered range of wavenumbers.

The conclusions obtained from the analytical expressions of $E_{\delta_{sphere}}$ and $E_{H_{sphere}}$ can be generalized to larger values of $p_{min}$. All spectra show oscillations of period $k_b$ and a final power-law decay. At wavenumbers much smaller than $1/R$, they all show a similar pattern and behave as $k^2$ as $k$ tends towards zero. At large $k$, the spectrum behaves as $k^{-2p-2} \sin^2(R\|\mathbf{k}\|+\phi)$. At this stage, two main conclusions can be drawn. On the one hand, the decay is controlled by the order of the singularity: the more regular the function, the lower $p_{min}$, the faster the decay. On the other hand, the period of the oscillations and the wavenumber where the singularity spectrum is maximum are controlled by the droplet size, whatever the nature of the discontinuity.

The power-law decay generated by an interface singularity of a field should not be confused with the power-law subrange exhibited by a regular physical field, such as the inertial $k^{-5/3}$ of single-phase flow turbulence or the $k^{-3}$ of bubble-induced agitation (Lance and Bataille, 1991; Risso, 2018). Such physical subranges always have a cutoff at a certain wavenumber: the Kolmogorov microscale in the case of turbulence (Pope, 2000), the size of smaller bubble disturbances in the model of bubble-induced spectrum by Risso (2011). A singularity power-law never stops and always exceeds the regular part of a physical signal at large wavelengths.

At scales close to the size of a droplet, the question of which one dominates in a physical signal between the regular part or the singularity has no general answer and must be examined in each specific situation.

#### 3.3. Effects of deviation from sphericity and of interface thickness

In real situations, droplets or bubbles are often not spherical. Also, spatial resolution is finite, which means that a fluid interface can never be described by a surface of zero thickness. In particular, numerical simulations based on the sharp-interface formulation have to deal with finite-size meshes that contain the two phases. In this section, we examine the consequences on the spectrum of a singularity of the drop non-sphericity, defined by an aspect ratio $a/b$, and of a finite interface thickness, $\epsilon_s$. In what follows, $R$ is defined as the radius of a sphere of the same volume and all calculations have been performed on the same mesh as that used in the previous section.

We begin with the case of the Heaviside singularity ($p_{min} = 0$). Fig. 3 shows numerical spectra $E_{V_0}$ obtained by considering an oblate ellipsoidal interface with $a/b = 1.35$ and 2, which is still sharp ($\epsilon_s = 0$), together with the theoretical result (Eq. (15)) and the numerical result for a sphere (Fig. 2). Deformation clearly dampens oscillations, and the greater $a/b$, the lower the oscillation amplitude. This is easily understood by noting that, for an ellipsoid, the distance $2\beta$ between two diametrically opposed points on the interface varies between $2a$ and $2b$ depending on the direction considered. The three-dimensional spectrum therefore mixes oscillations of various periods $\pi/\beta$, which cannot remain in phase as $k$ increases, causing them to cancel each other out. Fig. 3 also shows the spectrum computed by considering spherical interfaces of random radii uniformly distributed between $R$ and $3R$. One sees that a similar damping of the oscillations occurs when one considers a population of drops of different sizes, with the
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Numerical spectra of three-dimensional scalar fields of various singularity orders: $p_{min} = -1$ (red), $p_{min} = 0$ (green), $p_{min} = 1$ with $J_1 = 1$ (blue), $p_{min} = 2$ with $J_2 = 10$ (orange). The dashed straight lines correspond to $k^{-2-2p_{min}}$, with the exponent written on right side. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Numerical spectra of the indicator function of oblate ellipsoidal balls ($p_{min} = 0$, Heaviside type) of single size with various aspect ratios $a/b$ or spherical balls with random sizes. The dashed line corresponds to $k^{-2p_{min}}$.</figcaption></figure>

variation of $R$ from drop to drop substituting for the variation of $2\beta$ with direction. In any case, this process has no effect on the baseline of the spectrum, the decays of which remains the same.

Fig. 4 presents numerical $E_{V_0}$ obtained by considering a spherical interface of finite thickness. The interface has been smoothly extended over a thickness $\epsilon_s$ by replacing the Heaviside function in Eq. (17) by an error function: $V_0(\chi) = 1 - \{\operatorname{erf}[(\chi - 1)/\epsilon_s]\}/2$. For $\epsilon_s = R/100$, the interface thickness is below the grid resolution $\Delta = R/25$, and the numerical result is the same as that of Fig. 3, and in agreement with the theory. The effect of $\epsilon_s$ becomes visible when it becomes larger than $\Delta$. A finite thickness affects the small scales by increasing the final decay of the spectrum, without having any effect on the oscillations. However, it has no impact on the scales that are much larger than $\epsilon_s$. This means that a very strong localized gradient of a regular field, such as that of the velocity in a large-Reynolds number boundary layer of thickness $\epsilon_{BL}$ around a drop, would generate a similar spectral signature as a velocity jump at scales much larger than $\epsilon_{BL}$. We can also wonder whether singularities associated with intermittency (Le Berre et al., 2023) could also affect turbulence spectra in a similar way.

Now, we examine the case of the Dirac singularity ($p_{min} = -1$). The Dirac delta function in Eq. (16) is replaced by finite boxcar function

of height $\epsilon_s$ and width $1/\epsilon_s$. Fig. 5 shows the numerical spectra for spherical interfaces of various thicknesses, while Fig. 6 shows the same results for oblate ellipsoids of aspect ratio $a/b = 2$. The conclusions are the same as for the Heaviside singularity. Increasing the interface thickness causes a faster decay of the spectrum at large wavenumber. Increasing the deformation generates a damping of the oscillations.

### 3.4. Spectral signatures of singularities in physical cases

We showed in the previous section (Eq. 2) that a finite jump in the fluid properties (density or viscosity) across the interface or a non-zero interfacial tension causes singularities in the fields that characterize the flow. In the present section, we examined the signature of basic singularities of various orders corresponding to fluid interfaces of droplets. For spherical drops having all the same radius $R$, the spectrum of the interface singularities is maximum close to $k_b = \pi/R$, shows oscillations of period $k_b$ and decays at small scales as $k^{-2\rho^{-2}}$. For polydisperse or deformed drops the oscillations may be absent, while the baseline of the singularity spectrum is still present.

A real physical field is more complex, since it includes a contribution from the regular part that accounts for the dynamics in the
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Numerical spectra of the indicator function of spherical balls ($p_{min} = 0$, Heaviside Type) with smooth interface of various thicknesses $ε_s$. The dashed line corresponds to $k^{-2p_{min}}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Numerical spectra of the indicator function of spheres ($p_{min} = -1$, Dirac type) of various thicknesses $ε_s$. The dashed line corresponds to $k^{-2p_{min}}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Numerical spectra of the indicator function of an oblate ellipsoids ($p_{min} = -1$, Dirac type, $a/b = 2$) with a smooth interface of various thicknesses $ε_s$. The dashed line corresponds to $k^{-2p_{min}}$.</figcaption></figure>
bulk phases, in addition to singular contributions that account for the interfaces. Is it possible to distinguish between the regular and the singular parts? First of all, we may try to attribute the presence of oscillations as $sin^2(kR + \phi)$ to the singular part. This first approach has two limitations. On the one hand, for drops moving at large Reynolds number relative to the carrying phase, oscillations may be due to the presence of a dynamic boundary layer at the drop surface. On the other hand, spectra calculated from DNS of highly deformed drops with a broad size distribution immersed in a turbulent field by Crialesi-Esposito et al. (2022) show no oscillations, for the reason explained above. Secondly, we can rely on small scales to detect the singular part. In principle, the singular part undergoes a never-ending power-law decay as $k$ approaches infinity. It must therefore eventually emerge from the regular part, which has a physical cutoff, corresponding to the Kolmogorov microscale in a turbulent flow. However, the dissipative range is most of the time not resolved, which means that we need to detect the singular contribution at larger scales.

At this stage, it is therefore difficult to anticipate a general method for distinguishing the regular part from the singular part, which would probably have to be done on a case-by-case basis. Furthermore, even if such a distinction can be made, the question of its physical relevance remains open.

## 4. Direct numerical simulations of bubble-induced agitation

In the previous section, we examined the theoretical consequences of basic jump conditions on the spectrum of some reference fields. We are now interested in analyzing their practical effects on the spectrum of physical fields, such as the mixture velocity and the volume forces involved in the momentum balance. Since the spectrum of interface discontinuities is significant for wavelengths of the order of the drop size $R$ and smaller, it has negligible consequences when most of the energy is contained at scales much larger than $R$, as in the case of droplets or bubbles smaller than the Kolmogorov microscale in a turbulent field. On the other hand, the case of a swarm of rising bubbles in an otherwise quiescent fluid is particularly impacted, since the energy of the fluctuations is generated by the bubbles at scales close to $R$.

The fluid velocity in both phases can hardly be measured in experiments. Most of the time, only the velocity field of the outer phase is measured in two-phase dispersed flows. In addition, the flow in the close vicinity of bubbles or droplets is often removed because of spurious measurements. On the other hand, DNS using sharp interface methods provide exhaustive information on all physical fields at any location. In what follows, we analyze in detail the results of DNS of a homogeneous swarm of bubbles rising at high Reynolds number.

### 4.1. Presentation of the simulations

Bubbles having all the same size are initially homogeneously distributed over a triperiodic cubic domain and then rise under the action of buoyancy. After an initial transient, the flow statistics cease to evolve in time. All subsequent results are obtained in this statistically steady state. The simulations are carried out with the open-source TRUST/TrioCFD code, which uses a front-tracking method to solve the continuity and Navier-Stokes equations in the sharp-interface form (Eqs. (3), (4)) across the entire domain, including both the gas and liquid phases. The code has been described in detail and validated for high-Reynolds-number bubbly flow with parameters similar to those of the present work in du Cluzeau et al. (2019), du Cluzeau et al. (2022).

The physical parameters are: $\mu_c = \mu_d = 3.73 \times 10^{-4}$ Pa s, $\rho_c = 1.1713 \times 10^3$ kg/m$^3$, $\rho_d = 87.545$ kg/m$^3$, $\sigma = 18.05 \times 10^{-3}$ N/m, $d = 2R = 10^{-3}$ m, $g = 9.81$ m/s$^2$. This gives an Archimedes number $Ar = \frac{\rho_c d \sqrt{(1-\rho_d/\rho_c)g d}}{\mu_c} \approx 300$ and a Bond number $Bo = (\rho_c - \rho_d)g d^2/\sigma \approx 0.6$, corresponding to slightly deformed bubbles in a flow regime dominated by inertia. The density ratio, $\rho_d/\rho_c \approx 1/13$, and viscosity ratio, $\mu_d/\mu_c = 1$, are outside the range expected for bubbly flow at common room temperature and

<figure><img src="image_10.png" /><figcaption>Fig. 7. Snapshot showing the vertical velocity in a vertical plane ($\alpha = 6\%$, refined mesh).</figcaption></figure>

pressure. Note also that surface tension is kept constant. These conditions have been chosen to simplify the entanglement of the various singularities in the analysis of the spectrum of the various terms. In particular, the theoretical analysis suggests that a viscosity ratio equal to unity should ensure a more regular velocity and, together with a constant $\sigma$, to a more regular viscous force.

The simulations are performed on a 432³ regular mesh in a cubic domain of side $L = 31.7 \times 10^{-3}$ m or $L = 15.85 \times 10^{-3}$ m. In the larger domain, the mesh-grid spacing is thus $\Delta = 73.5 \times 10^{-6}$ m, while in the smaller domain it is $\Delta = 36.77 \times 10^{-6}$ m.

Two gas volume fractions have been considered: $\alpha = 0.03$ and $0.06$. A snapshot of the flow is shown in Fig. 7. The case $\alpha = 0.03$ has been solved on the larger domain ($L = 32d$) with the lower resolution ($\Delta = d/14$). It gives the following overall statistics for the two-phase mixture: average bubble relative velocity $u_r = 0.17$ m/s, half the variance of the mixture velocity $K_m = 5.4 \times 10^{-2} u_r^2$ and dissipation rate per unit volume $\epsilon_m = 7.45 \times 10^{-3} \rho_c u_r^2/d$. The average bubble Reynolds number, Weber number and aspect ratio are: $Re = \frac{\rho_c d u_r}{\mu_c} = 528$, $We = \frac{\rho_c d u_r^2}{\sigma} = 1.83$ and $a/b = 1.39$.

The case $\alpha = 0.06$ has been solved on both the larger domain ($L = 32d$, $\Delta = d/14$) and the smaller domain with a better resolution ($L = 16d$, $\Delta = d/27$). The lower resolution gives: $u_r = 0.150$ m/s, $K_m = 11.2 \times 10^{-2} u_r^2$, $\epsilon_m = 18.5 \times 10^{-3} \rho_c u_r^3/d$, $Re = 470$, $We = 1.45$ and $a/b = 1.33$. The greater resolution gives: $u_r = 0.147$ m/s, $K_m = 11.8 \times 10^{-2} u_r^2$, $\epsilon_m = 22.0 \times 10^{-3} \rho_c u_r^3/d$, $Re = 462$, $We = 1.40$ and $a/b = 1.37$. Rigorously, even our finest resolution is not enough to exactly capture the boundary layer around the bubbles at such a Reynolds number (Innocenti et al., 2021), which explains small differences between the two resolutions. In the literature, such a limitation is encountered in many numerical simulations involving a large number of dispersed particles. This is not a serious limitation for the aim of the present work, which is to analyze the spectral signature of interface discontinuities, as we shall see later when comparing the results of our coarse and refined meshes.
<figure><img src="image_11.png" /><figcaption>Fig. 8. Normalized spectrum of the density (DNS results). The straight line shows the asymptotic evolution $k^{-2\rho_{min}-2}$, with $\rho_{min} = 0$.</figcaption></figure>

## 4.2. Spectrum of density and velocity

In a numerical simulation, even in the case where the mesh grid is fitted to follow the interface, the fields are described with a finite resolution $\Delta$. The consequences for the spectrum have been discussed in Section 3.3 where an interface of finite thickness has been considered. A coarser resolution leads to a faster decrease in the spectral footprint of the singularities. In addition, with the front tracking method used here, the two-dimensional mesh that describes the interface differs from the three-dimensional one on which is discretized the bulk flow. In a cell containing both phases, the values of the density and viscosity have to be interpolated between those of the two phases, as it is also the case with common other sharp-interface methods, such as VOF or Level-Set. Here the density is determined as the arithmetic mean of the densities of the two phases weighted by their volume fraction within the cell considered. Since the viscosity is assumed constant, no interpolation is required. The interfacial force is computed from the curvature of the two-dimensional interface mesh. Then, $F_\sigma$ is distributed over the cells of the three-dimensional mesh that are located at a distance lower than $\Delta$ from the interface. Therefore, the Dirac delta function is smoothed over a width of approximately $2\Delta$.

The effect of a region where the mixture is described by a fluid with intermediate properties can hardly be anticipated, but it may retribute the singularities between the different terms of the momentum balance in a way that probably depends on the numerical scheme.

Let us consider first the consequences of the numerical discretization on the spectra of two basic fields: the density $\rho$, which is only impacted by the finite resolution, and the velocity $u$, which is affected by both the finite resolution and by the way the continuity and Navier-Stokes equations are solved. Note that the spectra are normalized by using $k_b$ and their variance $\langle(\cdot)^2\rangle$, in order to ease comparison between the two volume fractions.

Fig. 8 shows the spectrum of the density for the three different simulations and the analytical solution $E_{H_{sphere}}$ (Eq. (14)). Whatever the volume fraction or the resolution, the spectrum shows oscillations as $\sin^2(Rk)$ in agreement with the analytical solution. Oscillation amplitudes are smaller due to departure from sphericity. Cases $\alpha = 3\%$ and $\alpha = 6\%$ at the coarser resolution match at all wavenumbers and agree at large scales with the case at finer resolution. The coarser cases show a decay in agreement with the theoretical $k^{-2}$ up to $k = 3k_b$, while the finer case follows it up to $6k_b$, showing that halving the mesh spacing

actually halves the scale from which the interface begins to be seen as less sharp than it should be. Regarding the density, the numerical spectrum behaves as predicted by the theory.

Fig. 9 shows the spectrum $E_u$ of the velocity. The spectra of the three cases are very similar, featuring oscillations with a period $k_b$. We note only differences at very large scales due to finite domain size, and at very small scales close to the Nyquist wavenumber. This indicates that the velocity spectrum is very robust to changes in resolution. The spectrum reaches a maximum near $k_b$, followed by a $k^{-3}$ subrange, which is known as a signature of the bubble-induced agitation (Risso, 2018). Then, it shows a steeper decay as $k^{-4}$ for $k > 3k_b$, followed by an even steeper one as $k^{-n}$, with $4 < n < 5$, for $k > 6k_b$. Under present conditions, with no viscosity jump or interfacial tension gradient, the predictions for $u$ are $p_{min} = 2$, which corresponds to a final decay as $k^{-6}$. The fact that $u$ is less regular than expected is probably a consequence of numerical approximations in the cells crossed by the interface.

## 4.3. Spectrum of the terms of the momentum equations

We examine now the spectra of the five terms involved in the momentum balance (4) from the DNS at $\alpha = 6\%$ with either the coarse or the refined resolution: buoyancy $(\rho - \rho_m)g$, interfacial force $F_\sigma$, pressure gradient $-\nabla P$, inertia $-\rho D_i u$ and viscous diffusion $\nabla \cdot (2\mu S)$. Note that $\rho_m = (1-\alpha)\rho_c + \alpha\rho_d$ is the average mixture density and we have chosen to include the mean hydrostatic pressure gradient $\rho_m g$ in the buoyancy term so that it takes into account all the external forces acting on the two-phase system. In this section, all spectra are normalized using $k_b$ and the variance of $(\rho - \rho_m)g$. It is worth recalling that the spectrum is quadratic in the quantity considered, so that the balance of the five terms does not lead to a simple constraint on their spectra. When only two terms are in equilibrium, their spectrum is the same, but it is not as simple for more terms.

Since $g$ is constant, the spectrum of $(\rho - \rho_m)g$ (Fig. 10) is similar to that of $\rho$ (Fig. 8) and leads to the same conclusions.

Fig. 11 presents the spectra of $F_\sigma$ and $-\nabla P$. At large scales, for wavenumbers smaller than $\approx 2k_b$, both are independent of the mesh refinement. At small scales, the Dirac-type singularity introduced by interfacial tension gives a constant intensity. We indeed observe a plateau of the baseline before a cut-off due to the finite resolution, which appears at a smaller $k$ at lower resolution. Still, at small scales,
<figure><img src="image_12.png" /><figcaption>Fig. 9. Normalized spectrum of the velocity (DNS results). Note that the expected final decay is $k^{-6}$ ($p_{min} = 2$) for the case without viscosity jump.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Normalized spectrum $E_{(\rho-\rho_m)/\rho}$ of buoyancy force at $\alpha = 6\%$ (DNS results).</figcaption></figure>

the spectra of $F_\sigma$ and $-\nabla P$ are identical, which confirms that the Dirac-type singularities cancel out entirely between the interfacial force and the pressure gradients, in agreement with the equilibrium of normal stresses at the interface. The theoretical spectra $E_{H\text{sphere}}$ and $E_{\delta\text{sphere}}$ are also reported. It turns out that the oscillations of the pressure gradient are in phase with $E_{H\text{sphere}}$ at large scales and with $E_{\delta\text{sphere}}$ at small scales, which underlines the fact that pressure undergoes both a jump in its value and its derivative across the interface. So far, the spectra of the forces are in agreement with the theoretical predictions of Section 3.

The spectra of the terms that explicitly involve the velocity, $-\rho D_t u$ and $\nabla \cdot (2\mu S)$, show a somewhat different picture (Fig. 12). Both are independent of the mesh refinement for $k/k_b < 3$ for $-\rho D_t u$ and $k/k_b < 4$ for $\nabla \cdot (2\mu S)$. The fact that they show no oscillations while showing no decay in this range recalls the footprint of a Dirac-type singularity. As demonstrated earlier, these terms should not present such a singularity. That spectral behavior at small scales is therefore a consequence of the approximation done by using a discrete description of the interface. Even if we are not able to relate the spectra of these two quantities to the detail of the numerical scheme, it is interesting to note that the

small-scale plateaus reached by $E_{\rho D_t u}$ and $E_{\nabla \cdot (2\mu S)}$ are lower than those of $E_{\nabla P}$ and $E_{F_\sigma}$, suggesting that a part of the Dirac-like singularity of the latter two is redirected to the former two due to the approximate numerical treatment near the interface.

In summary, this section leads to the following conclusions. The spectral signature of the singularities is as predicted by theory for buoyancy, pressure gradient and interfacial force, but different for inertia and viscous force due to numerical approximations in the cells crossed by the interface. However, all spectra are remarkably independent of the mesh spacing, the effect of which is only significant at the small scales where the effect of singularities becomes dominant.

### 4.4. Energy budget in the spectral domain

The spectral energy balance of the kinetic energy of the two-phase mixture can be obtained in the same way as for a single-phase flow. The Fourier transform (Eq. (10)) of each term in the Navier–Stokes equation (Eq. (4)) is calculated, and contracted with the complex conjugate $\hat{u}^*$
<figure><img src="image_14.png" /><figcaption>Fig. 11. Normalized spectra $E_{v,P}$ of pressure gradient, and $E_{F_u}$ of interfacial force at $\alpha = 6\%$ (DNS results).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Normalized spectra $E_{\rho,D,u}$ of inertia, and $E_{v,(2\mu S)}$ of viscous force at $\alpha = 6\%$ (DNS results).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. Normalized spectral density of the terms of the energy balance for the refined grid at $\alpha = 6\%$ (DNS results).</figcaption></figure>
of the Fourier transform of the velocity:

$$ \underbrace{+\mu \nabla^2 \mathbf{u} \cdot \hat{\mathbf{u}}^*}_{\text{Dissipation}} + \underbrace{(\rho - \rho_m) \mathbf{g} \cdot \hat{\mathbf{u}}^*}_{\text{Production}} - \underbrace{-\rho D_t \mathbf{u} \cdot \hat{\mathbf{u}}^*}_{\text{Inertia}} + \underbrace{-\nabla P \cdot \hat{\mathbf{u}}^*}_{\text{Pressure}} + \underbrace{\mathbf{F}_\sigma \cdot \hat{\mathbf{u}}^*}_{\text{Interfacial tension}} = 0. \quad (20) $$

Each term is a complex-valued field in the wavenumber vector k. We take its real part and integrate it over shells of constant k = ||k|| (Eq. (11)) to finally obtain a real-valued function of k. This means we will not be examining flow anisotropy in what follows.

Note that according to the Plancherel's theorem, each term T of the momentum equation satisfies the relation: $\int T(\mathbf{x})\mathbf{u}(\mathbf{x})d\mathbf{x}^3 = \int T(\mathbf{k})\mathbf{u}(\mathbf{k})^* d\mathbf{k}^3$, which means that the integral over the wavenumbers of each term of the spectral balance (Eq. (20)) is equal to the average power of T. In particular, the first term can be expressed with the velocity spectrum as $\mu k^2 E_u$, and, its integral over k, noted $-\epsilon$, is equal to the dissipation rate of kinetic energy. It can thus be interpreted as the spectral density of the dissipation rate.

Since we are considering a statistical steady state, the integral over k of the second term is equal to $\epsilon$ and corresponds to the rate of energy supplied to the system by external forces. The buoyancy term therefore represents the spectral density of the production rate of kinetic energy.

The integral over k of each of the three remaining terms is zero, which means that inertia, pressure and interfacial forces do not contribute to the total amount of power, their role being to transfer energy between scales. Note that it is only true in steady state.

Before examining the DNS results, it is useful to comment on a specificity of the inertia term when the densities of the two phases are different, as it was misinterpreted by one of the few previous works that dealt with such a situation (Pandey et al., 2020). The inertia term can be split into two parts as $-\rho D_t \mathbf{u} \cdot \hat{\mathbf{u}}^* = -\rho \partial_t \mathbf{u} \cdot \hat{\mathbf{u}}^* - \rho \mathbf{u} \cdot \nabla \mathbf{u} \cdot \hat{\mathbf{u}}^*$. When density is constant, $\rho$ can be taken out of the Fourier-transform operator and the first term of the right-hand-side can be re-written $\partial_t(\hat{\mathbf{u}} \cdot \hat{\mathbf{u}}^*)$, which is zero at steady state. On the other hand, when $\rho$ is not constant, $-\rho \partial_t \mathbf{u} \cdot \hat{\mathbf{u}}^*$ does not vanishes at steady state and contributes to the transfer of energy between the scales. It is important to stress that while $-\rho D_t \mathbf{u} \cdot \hat{\mathbf{u}}^*$ is Galilean invariant (for $||k|| > 0$), the two parts of its decomposition are not. This is not a problem when $\rho$ is constant, since we can always choose the reference frame where the average velocity is zero. However, it means that this decomposition is irrelevant in the present case involving two phases with different densities moving at different velocities. Pandey et al. (2020) adopted a different perspective. They decided to decompose the fluid acceleration term into a Eulerian time derivative and a convective part in a way so that the former term vanishes in steady state. For each term T of the Navier-Stokes equations, they defined its spectral density as: $\hat{T} \cdot \hat{\mathbf{u}}^* + \frac{\hat{T}^*}{\rho} \rho \hat{\mathbf{u}}$.

The Eulerian time-derivative term is now written $\partial_t(\hat{\mathbf{u}} \cdot \hat{\mathbf{u}}^*)$, which is indeed zero in steady state. However, this comes at a high price, as the physical meaning of terms involving $T/\rho$ is unclear. Moreover, the division by $\rho$ strengthens singularities and their effect on the spectrum. What is more, their mathematical relevance is uncertain, in particular regarding $F^\sigma/\rho$, which does not represent the interfacial stress jump and involves a ratio between the Dirac and Heaviside distributions. We therefore see no advantage in adopting this approach and use the classic method defined by Eq. (20).

Fig. 13 presents the spectral energy budget (Eq. (20)), where all spectra have been normalized using $k_b$ and $\epsilon$. Oscillations with a period $k_b$ of significant amplitudes are visible on all terms, indicating that the singular parts are far from being negligible.

The production term shows a large positive and almost symmetric peak around $k = 0.75k_b$. Beyond $k = 1.5k_b$, it becomes a secondary term dominated by oscillations of decaying amplitude. Between $k = 0$ and $1.5k_b$, production is balanced by dissipation, inertia and interfacial tension. In contrast with the turbulence of a single-phase flow, transfer here is ensured to a comparable extent by inertia and interfacial force,

while the role of the pressure gradient remains negligible. However, the inertial and interfacial contributions exhibit very strong oscillations, causing them to change sign and making it difficult to understand them separately. The physical interpretation of the balance is simplified if inertia, interfacial force and pressure gradient are added together to construct the total transfer term, represented by the gray line in the figure. This makes clear the existence of two regimes, both corresponding to an energy transfer from large to small scales. Before the production peak ($0 \le k \le 0.75k_b$), the energy supplied by the work of buoyancy is mainly balanced by the transfer term, which transports energy to smaller scales. After the dissipation peak ($k \ge k_b$), the energy supplied from larger scales by the transfer term is dissipated. This picture is close to single-phase flow turbulence, with the notable difference that there is no scale separation between the peak of production and the peak of dissipation, which explains the absence of a $k^{-5/3}$ inertial subrange. It should also be noted that there is no subrange where production and dissipation are in equilibrium, in contrast to what was postulated by Lance and Bataille (1991) to explain the existence of the $k^{-3}$ subrange. The same conclusion was recently reached by Zamansky et al. (2024) who showed that the $k^{-3}$ subrange corresponds to a constant shear rate imposed by the bubble wakes and may be explained by a mechanism of return to isotropy that takes place once energy production by bubbles has ceased.

The total transfer possesses another property, which is of major interest for the purpose of the present work. It has no oscillations, indicating that the singularities of its three components cancel each other out, leaving only the regular part. Since the sum of the production and dissipation terms balances the transfer term, their singularities are the opposite of each other. It is therefore sufficient to identify the singular part of one of them to obtain two distinct, regular and singular spectral energy balances.

Since buoyancy acts on the bubble scale, it does not supply energy at scales that are significantly smaller than R. Consequently, the regular part of the production spectrum must decay very fast after $k_b$. The DNS spectrum plotted in log scale in Fig. 14a shows that the decay after the peak is interrupted by strong oscillations on a power-law decaying baseline, which can be attributed to the singular part. The regular part of the production spectrum can thus be obtained by low-pass filtering of the DNS spectrum. This has been done in Fig. 14a, where the blue line shows the production spectrum multiplied by a sharp filter, $\exp(-(k/k_c)^4)$, with $k_c = 1.3k_b$. (Changing the cutoff wavenumber $k_c$, slightly changes the result in the vicinity of $k_c$, but has no effect beyond). The singular part (gray dashed line), obtained as the difference between the total spectrum and the regular part, combines several types of singularities since its oscillations slowly evolves from a $\cos^2(kR)$ behavior to a $\sin^2(kR)$ behavior, and its power decay from a slope $-3$ to $-4$.

In contrast with production, the dissipation is only expected to experience a viscous cutoff at scales much smaller than the bubbles, so there is no easy way to distinguish between its singular and regular parts. Nevertheless, as noted above, its singular part is the opposite of that of the production spectrum and its regular part is obtained by subtracting it from the total dissipation spectrum. Fig. 14b shows a log plot of the total dissipation spectrum, as well as its regular and singular parts obtained in that way. The regular part shows a clear $k^{-1}$ evolution in the range from $k = 0.8k_b$ (just after the production peak) to $2-3k_b$.

Provided the pressure gradient, inertia and interfacial contributions are considered all together in a total transfer term, we eventually obtain two separate energy budgets, one for the regular part, the other for the singular part. The production (Fig. 14a) and the dissipation (Fig. 14b) involve both a regular and a singular part, while the total transfer (Fig. 14b) is regular.

Finally, the singular and regular parts of the velocity spectrum are obtained by dividing the corresponding part of the dissipation spectrum by $\mu k^2$. The regular part shows a $k^{-3}$ evolution in the region where the dissipation spectrum shows a $k^{-1}$ behavior.
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 14. Normalized spectral density of the terms of the energy balance and energy spectrum. DNS at $\alpha = 6\%$ and CGS at $\alpha = 7.5\%$. (a): Production, (b): Dissipation, (c): Transfers, (d): Energy.</figcaption></figure>

## 5. Comparison with coarse-grained simulations

In this work, the Navier-Stokes equations have been directly solved without additional modeling and the effect of the singularities due to the presence of the interfaces have been filtered a posteriori out from the numerical results. An opposite approach is to filter the flow in the region close to the interfaces before solving the Navier-Stokes equations. In this case, all fields are regular. Such coarse-grained simulations (CGS) were introduced by Riboux et al. (2013) to model the agitation generated by a flow through an array of fixed bubbles and extended to simulate freely moving bubbles by Le Roy De Bonneville et al. (2021). The case of a homogeneous swarm of rising bubbles has been extensively investigated by this method in Zamansky et al. (2024), where a detailed analysis of the spectral energy budget was provided. In this section, we compare these CGS spectra to the DNS spectra presented in the previous section.

Here we briefly present the CGS simulations, focusing on elements that are useful for understanding the comparison with the DNS. A detailed description of the method is available in Le Roy De Bonneville et al. (2021), Zamansky et al. (2024). The carrying phase fills the entire domain without interruption. Its dynamics is described by the continuity and Navier-Stokes equations for an incompressible fluid of constant density and viscosity. The presence of the bubbles is accounted for in the Navier-Stokes equations by a forcing term $\mathbf{f}_{b \to f}$,

$$ \mathbf{f}_{b \to f}(\mathbf{x}, t) = - \sum_{b=1}^{N_b} \mathbf{F}_{f \to b} G(\mathbf{x} - \mathbf{x}_b(t)), \quad (21) $$

where a Gaussian kernel,

$$ G(\xi) = \frac{1}{(2\pi\sigma_G^2)^{3/2}} \exp(-\frac{\|\xi\|^2}{2\sigma_G^2}), \quad (22) $$

is used to distribute the force $-\mathbf{F}_{f \to b}$ exerted on the fluid by each bubble $b$. The location $\mathbf{x}_b$ of each bubble is calculated by solving Newton's second law, where the force $\mathbf{F}_{f \to b}$ exerted by the fluid included drag,

added-mass, buoyancy and Tchen forces, calculated by using classic expressions involving both the bubble velocity and the fluid velocity at the bubble location. The main difference between the DNS and CGS approaches lies in the forcing added to the Navier-Stokes equations. In DNS, $\mathbf{F}_\sigma$ is localized at the interfaces, while, in CGS, $\mathbf{f}_{b \to f}$ is spread over a region of scale $\sigma_G$.

The CGS simulations have been carried out in a triperiodic cubic domain of dimension $L = 70d$, on a regular mesh of spacing $\Delta = d/15$. The spatial resolution is thus the same as that of the coarse DNS mesh, but the momentum supplied by the bubble is here filtered at a scale of the order of that of the bubble, since $2\sigma_G = 0.28d$. The physical parameters are those of air bubbles of diameter $d = 2.5$ mm rising in water. The results reported here correspond to a volume faction $\alpha = 7.5\%$, and a bubble Reynolds number $Re = 760$.

The spectra are determined from the fluid pressure and velocity fields computed by the CGS. The spectral densities of energy ($E_\mu$) and dissipation ($\mu\nabla^2\mathbf{u}\cdot\hat{\mathbf{u}}^*$) are obtained exactly as in the DNS. The buoyancy force is applied to the bubbles, which then transfer it to the fluid through the momentum forcing $\mathbf{f}_{b \to f}$. The production spectrum is thus calculated from the work of the coupling force: $\hat{\mathbf{f}}_{b \to f} \cdot \hat{\mathbf{u}}^*$. The spectral transfers are, as in a single-phase flow, ensured by inertia ($-\rho\hat{D}_i\mathbf{u}\cdot\hat{\mathbf{u}}^*$) and pressure ($-\nabla\hat{P}\cdot\hat{\mathbf{u}}^*$), with a vanishing contribution of the pressure in homogeneous flow condition.

Normalized CGS spectra have been reported in Fig. 14 (red lines). Since $\rho$ and $\mu$ are constant and the coupling force $\hat{\mathbf{f}}_{b \to f}$ is smooth, CGS spectra do not have a singular part. The agreement between CGS spectra and the regular part of DNS spectra is remarkable and constitutes a cross-validation of both. On the one hand, the spreading of the momentum transferred from the bubble to the fluid over a region of comparable size to that of the bubble is confirmed to be significantly relevant, not only at scales larger than the filtering threshold, but also at significantly smaller scales. On the other hand, this gives us confidence in the method used to separate the regular and singular parts of the DNS spectra.
## 6. Conclusion

The spectral analysis of dispersed two-phase flows is complicated by the presence of interfaces between the two phases, across which some quantities are discontinuous. This difficulty can be circumvented by considering subdomains containing only one of the two phases, but at the price of the severe limitation: transfers of momentum and energy between phases cannot be described, making it impossible to construct a spectral balance. Analyzing fields extending over the entire domain is therefore a much more fruitful approach to understand the physical mechanisms underlying the dynamics of turbulent fluctuations.

In the case of bubbles or droplets transported by a carrying phase, the combination of jumps of density, viscosity and stresses across the interface leads the various fields characterizing the dynamics of flow mixture to experience singularities, which can be ordered from the less to the most regular by means of integer *p*: Dirac delta function ($p = -1$), Heaviside step function ($p = 0$), discontinuity of the first derivatives ($p = 1$), discontinuity of the second derivatives ($p = 2$)... It is important to stress that a given field has generally various interfacial singularities. For example, the pressure gradient generally contains a Dirac singularity caused by interfacial tension ($\kappa\sigma$) and a Heaviside singularity caused by the jump in the gravity force ($\rho g$). A field can be characterized by its most severe singularity, $p_{min}$. The velocity field is generally continuous but its derivatives are not, leading to $p_{min} = 1$.

The spectrum of any field combines a regular part, which represents the smooth evolution outside the interface, and a singular part, which reflects interfacial singularities. At large wavenumbers, the spectrum of a singularity of order $p$, follows a power-law decay as $-2p - 2$. Therefore there exists a wavenumber above which, the spectrum of a physical field is dominated by its singular part, and it ends up by being dominated by its most severe singularity, of order $p_{min}$. The wavenumbers where the singular part begins to become significant depends on the amplitude of the regular part relative to that of the singular part, and the way it evolves with $k$. When the interfaces are the surface of droplets or bubbles, the singular part of the spectrum reaches a maximum at a wavelength close to the droplet scale and follows a power-law at smaller scales. When the droplets are spheres of the same radius $R$, the spectrum of the singularities shows oscillations as $\sin^2(Rk + \phi)$, of period $k_b = \pi/R$. These oscillations vanish when the droplets are deformed or their size distributions is broad.

In a dispersed two-phase flow, the singular part of the spectra is therefore expected to be important for wavenumbers larger than $k_b$. For instance, it is insignificant on the major part of the turbulent spectra in the case of a turbulent flow laden with droplets of size smaller than the Kolmogorov micro-scale. However, it cannot be neglected when the droplet size is in the energetic wavelength range. In this work, we have investigated the case where its impact is the most important, a swarm of bubbles rising at a high Reynolds number, since the fluctuating energy is produced in a range around $k_b$. In addition, we have considered almost spherical bubbles of a single diameter $d$, for which the oscillations of the spectrum at a period $k_b$ are a visible signature of the singular part. We also chose the same viscosity for the two phases to simplify the analysis.

A homogeneous rising swarm, at volume fractions $\alpha = 3\%$ and $6\%$ and Reynolds number $Re = 500$, has been computed by Direct Numerical Simulations using a front tracking of the interfaces, with a mesh grid spacing $\Delta = d/14$ or $d/27$. The conclusions reached are independent of $\alpha$ and $\Delta$. The numerical spectra of the velocity and of all the terms of the momentum equations clearly show the signature of singularities of order $p = -1$ (Dirac type), $p = 0$ (Heaviside type) and $p \ge 1$ (continuous fields). Buoyancy, pressure gradient and interfacial forces show the singularities predicted by the theory. However, the two terms that explicitly imply velocity, namely inertia and viscous forces, are less regular than they should be. This indicates that, even though the momentum balance is locally satisfied everywhere, the numerical treatment in the meshes that are crossed by an interface, redistributes

the singularities between the different terms, leading to a less regular velocity field.

The spectral power densities of all the forces contributing to the energy balance were determined, in the same way as usually done in turbulence studies. This budget is made of five terms. The buoyancy term corresponds to the production of kinetic energy. The viscous term corresponds to the dissipation of kinetic energy into heat. The integral over the wavenumbers of the production term is the opposite of that of the dissipation term and equal to the dissipation rate of energy $\epsilon$. The inertia, interfacial tension and pressure-gradient terms transfer energy between scales without contributing to the total amount of power. It turns out that the singular part cancels out when we add the three transfer terms together, leading to a regular total transfer term. Then, taking advantage of the fact that the regular part of the production term decreases very sharply beyond $k_b$, its singular and regular parts can be separated. Since the singular part of the dissipation term must balance the singular part of the production term, we end up with two separate budgets, one for the regular part, one for the singular part. These results have been compared with those of Coarse-Grained Simulations, in which the momentum transfer between the bubbles and the fluid are filtered before the Navier-Stokes equations are solved. The CGS spectra are found to be in remarkable agreement with the regular part of the DNS spectra, validating both approaches to obtain the regular part of the spectra.

This work proves that it is relevant to analyze DNS spectra of bubbly flows computed by considering the fields over the entire domain. In the case of homogeneous bubble-induced agitation, it has been possible to separate the regular and singular parts of the spectral energy budget thanks to the two following properties: (1) the singular parts of the transfer terms cancels out when we add them, and (2) the production term decay very sharply after $k_b$. The question now is under which conditions these two properties are valid. In a future work, we shall address the case of coupled forced isotropic homogeneous turbulence and bubble-induced agitation.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgment

This work was granted access to the HPC resources of TGCC under the allocation 2023-A0142B07712 made by GENCI.

## References

Alméras, E., Mathai, Varghese, Lohse, Detlef, Toschi, F., 2017. Experimental investigation of the turbulence induced by a bubble swarm rising within incident turbulence. J. Fluid Mech. 825, 1091-1112. http://dx.doi.org/10.1017/jfm.2017.410.

Bacry, E., Muzy, J.F., Arnéodo, A., 1993. Singularity spectrum of fractal signals from wavelet analysis: Exact results. J. Stat. Phys. 70 (3-4), 635-674. http://dx.doi.org/10.1007/BF01053588.

Balachandar, S., Eaton, J.K., 2010. Turbulent Dispersed Multiphase Flow. Ann. Rev. Fluid Mech. 42 (1), 111-133. http://dx.doi.org/10.1146/annurev.fluid.010908.165243, Vol 43.

Brandt, L., Coletti, F., 2022. Particle-Laden turbulence: Progress and perspectives. Annu. Rev. Fluid Mech. 54 (1), 159-189. http://dx.doi.org/10.1146/annurev-fluid-030121-021103, ISSN (1545) 0066-4189-4479.

Broer, H., Takens, F., 1993. Mixed spectra and rotational symmetry. Arch. Ration. Mech. Anal. 124 (1), 13-42. http://dx.doi.org/10.1007/BF00392202.

du Cluzeau, A., Bois, G., Leoni, N., Toutant, A., 2022. Analysis and modeling of bubble-induced agitation from direct numerical simulation of homogeneous bubbly flows. Phys. Rev. Fluids 7 (4), 044604. http://dx.doi.org/10.1103/PhysRevFluids.7.044604.
Crialesi-Esposito, M., Rosti, M.E., Chibbaro, S., Brandt, L., 2022. Modulation of homogeneous and isotropic turbulence in emulsions. J. Fluid Mech. 940, A19. http://dx.doi.org/10.1017/jfm.2022.179.

Dodd, M.S., Ferrante, A., 2016. On the interaction of Taylor length scale size droplets and isotropic turbulence. J. Fluid Mech. 806, 356–412. http://dx.doi.org/10.1017/jfm.2016.550.

du Cluzeau, A., Bois, G., Toutant, A., 2019. Analysis and modelling of Reynolds stresses in turbulent bubbly up-flows from direct numerical simulations. J. Fluid Mech. 866, 132–168. http://dx.doi.org/10.1017/jfm.2019.100.

Dubrulle, B., 2019. Beyond kolmogorov cascades. J. Fluid Mech. 687, P1–63. http://dx.doi.org/10.1017/jfm.2019.98.

Duchon, J., Robert, R., 2000. Inertial energy dissipation for weak solutions of incompressible Euler and Navier-Stokes equations. Nonlinearity 13 (1), 249.

Freund, A., Ferrante, A., 2019. Wavelet-spectral analysis of droplet-laden isotropic turbulence. J. Fluid Mech. 875, 914–928. http://dx.doi.org/10.1017/jfm.2019.515.

Gelfand, I.M., Shilov, G.E., 1964. Generalized Functions Vol 1 Properties and Operations. Academic Press - New York and London.

Innocenti, A., Jaccod, A., Popinet, S., Chibbaro, S., 2021. Direct numerical simulation of bubble-induced turbulence. J. Fluid Mech. 918, http://dx.doi.org/10.1017/jfm.2021.288.

Kataoka, I., 1986. Local instant formulation of two-phase flow. Int. J. Multiph. Flow 12 (5), 745–758. http://dx.doi.org/10.1016/0301-9322(86)90049-2.

Lalanne, B., Villegas Rueda, L., Tanguy, S., Risso, F., 2015. On the computation of viscous terms for incompressible two-phase flows with Level Set/Ghost fluid method. J. Comput. Phys. 301, 289–307. http://dx.doi.org/10.1016/j.jcp.2015.08.036.

Lance, M., Bataille, J., 1991. Turbulence in the liquid-phase of a uniform bubbly air water-flow. J. Fluid Mech. 222 (-1), 95–118.

Le Berre, M., Lehner, T., Pomeau, Y., 2023. Singularities in turbulent flows: How to observe them? Physica D 443, 133563.

Le Roy De Bonneville, F., Zamansky, R., Risso, F., Boulin, A., Haquet, J.-F., 2021. Numerical simulations of the agitation generated by coarse-grained bubbles moving at large Reynolds number. J. Fluid Mech. 926, A20. http://dx.doi.org/10.1017/jfm.2021.670.

Lucci, F., Ferrante, A., Elghobashi, S., 2010. Modulation of isotropic turbulence by particles of Taylor length-scale size. J. Fluid Mech. 650, 5–55. http://dx.doi.org/10.1017/S0022112009994022.

Martínez Mercado, J., Chehata, G.D., Van Gils, D., Sun, C., Lohse, D., 2010. On bubble clustering and energy spectra in pseudo-turbulence. J. Fluid Mech. 650, 287–306.

Mendez-Diaz, S., Serrano-García, J.C., Zenit, R., Hernández-Cordero, J.A., 2013. Power spectral distributions of pseudo-turbulent bubbly flows. Phys. Fluids 25 (4), 043303.

Onsager, L., 1949. Statistical hydrodynamics. Il Nuovo Cimento 6 (2), 279–287.

Pandey, V., Mitra, D., Perlekar, P., 2022. Turbulence modulation in buoyancy-driven bubbly flows. J. Fluid Mech. 932, A19. http://dx.doi.org/10.1017/jfm.2021.942, ISSN (1469) 0022-1120-7645.

Pandey, V., Ramadugu, P., Perlekar, R., 2020. Liquid velocity fluctuations and energy spectra in three-dimensional buoyancy-driven bubbly flows. J. Fluid Mech. 884, R6. http://dx.doi.org/10.1017/jfm.2019.991.

Pikovsky, A.S., Zaks, M.A., Feudel, U., Kurths, J., 1995. Singular continuous spectra in dissipative dynamics. Phys. Rev. E 52 (1), 285–296. http://dx.doi.org/10.1103/PhysRevE.52.285.

Pope, S.B., 2000. Turbulent Flows. Cambridge University Press.

Prakash, V.N., Martínez Mercado, J., van Wijngaarden, L., Mancilla, E., Tagawa, Y., Lohse, D., Sun, C., 2016. Energy spectra in turbulent bubbly flows. J. Fluid Mech. 791, 174–190.

Riboux, G., Legendre, D., Risso, 2013. A model of bubble-induced turbulence based on large-scale wake interactions. J. Fluid Mech. 719, 362–387. http://dx.doi.org/10.1017/jfm.2013.12.

Riboux, G., Risso, F., Legendre, D., 2010. Experimental characterization of the agitation generated by bubbles rising at high Reynolds number. J. Fluid Mech. 643, 509–539.

Risso, F., 2011. Theoretical model for $k^{-3}$ spectra in dispersed multiphase flows. Phys. Fluids 23 (1), 011701.

Risso, F., 2018. Agitation, mixing, and transfers induced by bubbles. Annu. Rev. Fluid Mech. 50 (1), 25–48, Vol 43.

Roghair, I., Mercado, J.M., Annaland, M.Van.Sint., Kuipers, H., Sun, C., Lohse, D., 2011. Energy spectra and bubble velocity distributions in pseudo-turbulence: Numerical simulations vs. experiments. Int. J. Multiph. Flow 37 (9), 1093–1098.

Tryggvason, G., Bunner, Bernard, Tryggvason, Gretar, 2002. Dynamics of homogeneous bubbly flows Part 2. Velocity fluctuations. J. Fluid Mech. 466, 53–84. http://dx.doi.org/10.1017/S0022112002001180.

Zamansky, R., Le Roy De Bonneville, F., Risso, F., 2024. Turbulence induced by a swarm of rising bubbles from coarse-grained simulations. J. Fluid Mech. 984, A68. http://dx.doi.org/10.1017/jfm.2024.230.