An experimental study on the multiscale properties of turbulence in bubble-laden flows

Abstract

The properties of bubble-laden turbulent flows at different scales are investigated experimentally, focusing on the flow kinetic energy, energy transfer and extreme events. The experiments employed particle shadow velocimetry measurements to measure the flow in a column generated by a homogeneous bubble swarm rising in water, for two different bubble diameters ($2.7$ mm and $3.9$ mm) and moderate gas volume fractions ($0.26\,\%\sim 1.31\,\%$). The two velocity components were measured at high resolution, and used to construct structure functions up to twelfth order for separations spanning the small to large scales in the flow. Concerning the flow anisotropy, the velocity structure functions are found to differ for separations in the vertical and horizontal directions of the flow, and the cases with smaller bubbles are the most anisotropic, with a dependence on void fraction. The degree of anisotropy is shown to increase as the order of the structure functions is increased, showing that extreme events in the flow are the most anisotropic. Our results show that the average energy transfer with the horizontal velocity component is downscale, just as for the three-dimensional single-phase turbulence. However, the energy transfer associated with the vertical component of the fluid velocity is upscale. The probability density functions of the velocity increments reveal that extreme values become more probable with decreasing Reynolds number, the opposite of the behaviour in single-phase turbulence. We visualize those extreme events and find that regions of intense small-scale velocity increments occur near the turbulent/non-turbulent interface at the boundary of the bubble wake.

1. Introduction

Particle-laden flows are a central topic in fluid mechanics and omnipresent in nature and technology (Prosperetti & Tryggvason 2009; Balachandar & Eaton 2010). While great attention has been given to investigating flows containing suspensions of small inertial, heavy particles (Maxey 1987; Bec et al. 2006, 2007; Fox 2014; Gustavsson & Mehlig 2016; Ireland, Bragg & Collins 2016a,b; Hogendoorn & Poelma 2018; Dou et al. 2018a,b; Petersen, Baker & Coletti 2019; Tom & Bragg 2019; Berk & Coletti 2021), there has been less of a focus on bubbly turbulent flows, partly due to the increased complexity associated with performing experiments or simulations for such flows (Lohse 2018). Among the various topics relevant to bubbly flows, bubble-induced turbulence (BIT) is an important area for investigation, both for its own fundamental importance and also for understanding bubble motion (Magnaudet & Eames 2000; Loisy & Naso 2017; Mathai et al. 2018), deformation (Masuk et al. 2021a; Perrard et al. 2021), coalescence/breakup (Liao & Lucas 2010; Rivière et al. 2021), clustering mechanisms (Zenit, Koch & Sangani 2001; Maeda et al. 2021) and mixing processes (Alméras et al. 2015, 2019) in bubbly flows.

Recently, attention has been given to investigating single-point turbulence statistics in bubble-laden flows. Readers are referred to Risso (2018) and Mathai, Lohse & Sun (2020) for detailed reviews. For flows with low to moderate gas void fraction ($\alpha <5\,\%$), we summarize the following characteristics for homogeneous bubble swarms rising either within a background quiescent or weakly turbulent carrier liquid: (i) the liquid velocity fluctuations are highly anisotropic, with much larger fluctuations in the direction of the mean bubble motion (vertical direction in standard coordinates) (Mudde 2005; Lu & Tryggvason 2013; Ma et al. 2020b; Liao & Ma 2022); (ii) the probability density functions (PDFs) of all fluctuating velocity components are non-Gaussian, and the PDF of the vertical velocity fluctuation is strongly positively skewed, while the other two directions have symmetric PDFs (Riboux, Risso & Legendre 2010; Roghair et al. 2011; Riboux, Legendre & Risso 2013; Bouche et al. 2014; Lai & Socolofsky 2019); (iii) bubbles wakes introduce additional turbulence and enhance turbulence dissipation rates in the vicinity of the bubble surface (Santarelli, Roussel & Fröhlich 2016; du Cluzeau, Bois & Toutant 2019; Masuk, Salibindla & Ni 2021b); (iv) modulation of the liquid mean velocity profile due to interphase momentum transfer, resulting in modifications to the background shear-induced turbulence (Lu & Tryggvason 2008; Bragg et al. 2021).

Studies that have explored the multiscale properties of bubble-laden turbulence have mainly focused on investigating modifications to the energy spectra of the liquid velocity fluctuations due to the bubbles. Lance & Bataille (1991) were the first to find a power law scaling with a slope of approximately $-3$ for the vertical velocity fluctuation in a bubble-laden turbulent channel flows using hot-film anemometry. This scaling that emerges in BIT-dominated flows is in contrast to the classical $-5/3$ scaling that appears in many single-phase turbulent flows (Pope 2000). This $-3$ scaling was also confirmed by subsequent studies for all the components of the fluctuating fluid velocity, and has been observed both in experiments and direct numerical simulations (DNS) (Riboux et al. 2010; Mendez-Diaz et al. 2013; Pandey, Ramadugu & Perlekar 2020; Innocenti et al. 2021) for varying bubble properties.

Beyond the behaviour of the energy spectra, knowledge about the multiscale properties of bubble-laden turbulence is quite limited. A pioneering study on this topic is that by Rensen, Luther & Lohse (2005), who performed hot-film anemometry measurements in a bubble-laden water tunnel and found an increase of the second-order structure function for the two-phase case compared with the single-phase case for the same bulk Reynolds number, and that this increase was more pronounced at the small scales than the large scales. Moreover, they considered the PDFs of the velocity increments and used extended self-similarity (Benzi et al. 1993) to show that the flow intermittency is enhanced by the bubbles. Furthermore, they argued that, once the bubbles are present in the flow, the dependence of the flow properties on the actual bubble concentration is weak for the range they investigated ($0.5\,\%\leq \alpha \leq 2.9\,\%$). Similar behaviour was also observed later in Biferale et al. (2012) and Ma et al. (2021) when comparing the small-scale properties of bubble-laden and unladen turbulent flows.

To gain more insight into the multiscale energetics of bubble-laden turbulence, Pandey et al. (2020) computed the scale-by-scale average energy budget equation in Fourier space using results from an interface-resolved DNS with several tens of bubbles rising in an initially quiescent flow. They showed that on average there is a downscale energy transfer, just as occurs for the single-phase turbulence in three dimensions (Alexakis & Biferale 2018). A similar finding was also reported by two more recent DNS studies (Innocenti et al. 2021; Ma et al. 2021). An issue with our previous study (Ma et al. 2021), however, is that the conclusion was drawn based on a one-dimensional dataset, which only allowed us to construct the velocity increments for separations in the spanwise direction of the bubble-laden turbulent channel flow.

Another important point to be quantified is how the bubbles influence the anisotropy of the turbulent flow across the scales. Thanks to significant research efforts, single-phase multiscale anisotropy is understood in considerable detail (Sreenivasan & Antonia 1997; Biferale & Procaccia 2005). While phenomenological turbulence theories postulate a return to isotropy at small scales (Kolmogorov 1941; Frisch 1995), experimental and numerical data have shown persistent small-scale anisotropy (Shen & Warhaft 2000; Ouellette et al. 2006; Pumir, Xu & Siggia 2016), especially when considering high-order structure functions (Kurien & Sreenivasan 2000; Carter & Coletti 2017). In contrast to single-phase turbulence, where anisotropy is usually injected into the flow at the large scales (Chang, Bewley & Bodenschatz 2012), bubbles inject anisotropy into the flow at the scale of their diameter/wake, and this often corresponds to the small scales of the turbulence. As a result of this, bubble-laden turbulent flows can exhibit much stronger anisotropy at the small scales of the flow compared with their single-phase counterparts with the same bulk Reynolds number (Ma et al. 2021). This behaviour was quantified in our recent study by developing a new method based on the barycentric map approach (Banerjee et al. 2007), and the results also revealed that the bubble Reynolds number is the key factor responsible for governing the flow anisotropy, whereas the void fraction does not seem to play an important role, at least for the void fractions considered.

To advance the understanding of the multiscale properties of bubble-laden turbulence, in this paper we present an experimental study based on measurements of millimetre-sized air bubbles with (approximately) fixed shape/size rising in a vertical column of water, using high-resolution particle shadow velocimetry (PSV). The water is deliberately contaminated to allow for the no-slip condition to be satisfied on the bubble surface (Elghobashi 2019), similar to what is usually assumed in DNS studies. The bubbles are injected uniformly at a low to moderate gas void fractions to minimize the large-scale velocity fluctuations in the flow (Harteveld, Mudde & van den Akker 2003), and we focus on a region of the flow sufficiently far away from the walls of the column where the one-point statistics from both phases are approximately homogeneous across the measurement window. Unlike our previous study that used a one-dimensional DNS dataset (Ma et al. 2021), PSV measures two-dimensional, two-component velocity fields, providing access to both longitudinal and transverse structure functions associated with separations along two directions. Moreover, we record a large number of uncorrelated velocity fields, so that structure functions up to twelfth order can be computed. One of the ways the present work advances that in Ma et al. (2021) is that the resolution of the DNS used in Ma et al. (2021) was insufficient to resolve extreme fluctuations in the flow, and therefore only structure functions up to order four were considered. Here, since experimental data are used there is no question regarding whether the physics of the flow is properly captured in the measured quantities.

The rest of this paper is organized as follows. In § 2, we introduce the experimental set-up and the measurement techniques. We then first present the single-point statistics for both phases in § 3. The multipoint results are divided into three parts, namely, anisotropy in § 4, energy transfer in § 5 and intermittency in § 6.

2. Experimental set-up and measurement techniques

2.1. Experimental facility

The experimental data required for the present study have only recently become available thanks to the advance of the PSV measurement technique for two-phase flows (Hessenkemper & Ziegenhein 2018). We used this method in our experiments, which were conducted in Helmholtz-Zentrum Dresden – Rossendorf, Germany. The experimental set-up consists of a rectangular column (depth $50\ \mathrm {mm}$ and width $112.5\ \mathrm {mm}$) made of acrylic glass which is filled with tap water to a height of $1100\ \mathrm {mm}$ (figure 1). The temperature in the laboratory was kept constant at $20\,^{\circ }\mathrm {C}$ and the density and kinematic viscosity of the water are assumed to be constant with values $\rho _l=998\ \mathrm {kg}\ \mathrm {m}^{-3}$ and $\nu =1\times 10^{-6}\ \mathrm {m}^{2}\ \mathrm {s}^{-1}$, respectively.

Figure 1. Sketch of the bubble column used in the experiments (note that in the actual experiment, the number of bubbles in the column is $O(10^{3})$). The sketch is not to scale; the channel depth is many times larger than the bubble diameter. Top right shows a representation of the shallow depth of field seen from the side view and bottom right shows the sparger arrangement.

To suppress bubble deformation, $1000\ \mathrm {ppm}$ 1-Pentanol is added to the flow. As demonstrated by Tagawa, Takagi & Matsumoto (2014), 1-Pentanol at such a high concentration leads to a rapid full contamination of the bubbles. The thereby induced immobilization of the otherwise mobile bubble surface results in a nearly no-slip condition at the gas–liquid interphase (Takagi & Matsumoto 2011; Manikantan & Squires 2020). A reduction of the bubble rise velocity can also be clearly observed in the experiment by adding this amount of 1-Pentanol, which is explained by the Marangoni effect (Takagi & Matsumoto 2011). Another feature of 1-Pentanol is that the adsorbed surfactant inhibits bubble coalescence and break-up, so that in each experiment the bubbles are mono-disperse with a fixed bubble size.

Air bubbles are injected through several spargers that are inserted into $11$ holes that have been drilled into the bottom of the column (see bottom right of figure 1 for the sparger configuration). All the spargers are removable, hence, both the gas fraction and the bubble size can be varied by removing or replacing different spargers. In the present study, either all $11$ spargers or the $8$ outer spargers are used to ensure a homogenous bubble distribution at the measurement height. Here, we use two different spargers with the inner diameters $0.2\ \mathrm {mm}$ and $0.6\ \mathrm {mm}$, corresponding to constant gas flow rates of $0.04\ \mathrm {l}\ \mathrm {min}^{-1}$ and $0.1 \mathrm {l}\ \mathrm {min}^{-1}$ per sparger, respectively (our preliminary tests show that these flow rates lead to a stable production of bubbles of constant size). In summary, we consider two different bubble sizes, and for each bubble size we consider two gas void fractions. These four mono-dispersed cases are labelled as SmLess, SmMore, LaLess and LaMore in table 1 (Sm/La for smaller/larger bubbles and More/Less for higher/lower gas void fraction), including some basic characteristic dimensionless numbers for the bubbles. It should be noted that SmLess and LaLess do not have the same gas void fraction, nor do SmMore and LaMore. The naming convention is instead to distinguish two cases with the same sized bubbles but different void fractions, e.g. SmLess and SmMore correspond to cases with the same (small) sized bubbles, but one with greater void fraction than the other, and similarly for LaLess and LaMore.

Table 1. Selected basic statistics of the gas phase for the four investigated cases. Here, $\alpha _p$ is the averaged gas void fraction, $d_p$ the equivalent bubble diameter, $\chi$ the aspect ratio, $Ga\equiv \sqrt {|{\rm \pi} _\rho -1|gd_{p}^{3}}/\nu$ the Galileo number. The values of $Re_p$, the bubble Reynolds number, and $C_D$, the drag coefficient, are based on $d_p$ and the bubble to fluid relative velocity obtained from the experiment. Here, $\varDelta$ is one PSV grid, i.e. the smallest spatial resolution of the present experiment.

2.2. Measurement techniques

2.2.1. Liquid velocity measurement

Although it is possible to use conventional particle image velocimetry (PIV) with a laser sheet to measure the liquid velocity in bubbly flows, such a side-wise high intensity illumination can result in an inhomogeneous illumination due to unwanted lateral shadows of the bubbles as well as strong light scattering and reflection at the gas–liquid interfaces (Bröder & Sommerfeld 2007; Tropea 2011; Ziegenhein & Lucas 2016). To circumvent these problems, the very recently developed PSV method for dispersed two-phase flows (Hessenkemper & Ziegenhein 2018) is performed here. This method, which has been introduced and validated for single-phase flows in Estevadeordal & Goss (2005) and Goss & Estevadeordal (2006), employs volumetric direct in-line illumination with e.g. LED-backlights for the region of interest, due to which scattering effects are greatly reduced and no lateral bubble shadows occur. By using a shallow depth of field (DoF), sharp tracer particle shadows positioned inside the DoF region can be identified and the particle displacement is evaluated in a PIV-like manner. Hence sometimes this method is also called ‘PIV with LED’, or ‘particle shadow image velocimetry’ (Estevadeordal & Goss 2005). The top right of figure 1 shows a representation of such a shallow DoF based on the side view of the column. This kind of thin DoF is also frequently used in laser-based $\mu$PIV measurements, since a laser-sheet thickness below $0.5\ \mathrm {mm}$ is hard to achieve, while the DoF can be adjusted to even smaller depth expansions (Wereley & Meinhart 2010). With the present set-up, the effective DoF (i.e. depth of correlation (DoC) in the terminology of measurement) is calculated using (2) in Bourdon, Olsen & Gorby (2006) as $\mathrm {DoC}\simeq 370\ \mathrm {\mu }\mathrm {m}$, which is approximately one seventh of the smaller bubble diameter (table 1).

The liquid velocity measurements are performed along the $x_1$–$x_2$ symmetry plane at the centre of the depth ($x_3$). In the following, we refer to $x_1$ and $x_2$ as the vertical and horizontal directions, respectively. The measurement height is $x_1=0.65\ \mathrm {m}$ based on the centre of the field of view (FOV) (figure 1) to ensure that bubbles entering the measurement region have already lost any memory of the way they were injected. The flow is seeded with $10\ \mathrm {\mu }\mathrm {m}$ hollow glass spheres (HGS) – Dantec, and illuminated by a $200\ \mathrm {W}$ LED-lamp, which consists of $160$ small LED arrays arranged in a circular area of $100\ \mathrm {mm}$ diameter and located directly facing the high-speed camera in the $x_3$-direction. Since the HGS particles are hydrophilic, almost no flotation due to the added surfactant occurs. To fully capture the small-scale fluctuations of the flow, the Stokes number of the seeding particle, defined as $St=\tau _d\tau _f$ must be much smaller than unity, where $\tau _d$ is the particle response time and $\tau _f$ is the characteristic time scale of the flow. In the current experiments the highest $St$ is $O(10^{-3})$ for $\tau _f$ based on the estimated Kolmogorov time scale (table 2), whose estimate is explained in § 3. Image pairs for the velocity determination are acquired with a MotionPro Y3 high-speed CMOS camera, with a frame rate of $900$ f.p.s. for the smaller bubble cases and 1100 f.p.s. for the larger bubble cases.

Table 2. Selected basic statistics of the liquid phase for the four investigated flow configurations.

Before the velocity interrogation is performed, the comparatively large bubble shadows are masked and only interrogation areas outside the mask are considered. This is done by first applying a Gaussian and a median filter on the image to blur the small tracer particles, and then the much darker bubble shadows are masked (see figure 2) using a grey value threshold of 110. Velocity fields are generated using a multipass-refinement procedure, with $3$ refinement steps that result in a final interrogation window size of $32 \times 32$ pixels and $50\,\%$ overlap. Measurements are obtained by mounting a macro-lens (Samyang) with a focal length of $100\ \mathrm {mm}$ and an f-stop of $2.8$. This provides a FOV of $13.9\ \mathrm {mm}\,({H_1}) \times 19.5 \mathrm {mm}\,({H_2})$, with a pixel size of $14.7\ \mathrm {\mu }$m. Since this pixel size is even slightly larger than the seeding particle size, the so-called particle sharpening step (Hessenkemper & Ziegenhein 2018) is used, which artificially increases the particle size due to a swelling to the neighbouring pixels. Hence, an effective particle size of approximately $3\text{--}4$ pixels is achieved. In the resulting liquid velocity fields, the spatial resolution $\varDelta$ (PSV grid) is approximately $0.23\ \mathrm {mm}$ ($d_p/\varDelta \approx 11{\sim}17$, see table 1 for different cases). For more details on the PSV method used in the present measurements, including the whole PIV-like post-processing steps, we refer the readers to Hessenkemper & Ziegenhein (2018).

Figure 2. Instantaneous realization of velocity vector over the FOV in the case SmMore with a detected in-focus bubble (sharp interface) in the middle left. The horizontal/vertical dashed lines are examples, indicating the representative points along the lines where the structure functions were computed.

In figure 2 a typical instantaneous FOV with the bubbles and the resulting liquid velocity vector field are shown for the case SmMore. In this snapshot, while there is only one in-focus bubble (identified by the sharp interface), with strong wake entrainment identified in the velocity vector field, the other bubbles are blurred, i.e. out of focus. Tracer particles at these locations cannot be detected, so that no liquid velocity fields are available there.

To obtain converged high-order statistics, for each case considered $60\ 000$ velocity fields were recorded with an image pair acquisition rate of $1\ \mathrm {Hz}$. With $60\times 84$ interrogation windows in each velocity field, this yields in total ${\approx}3\times 10^{9}$ data points for each case.

2.2.2. Bubble statistics measurement

The bubble statistics and properties are evaluated with a separate set of measurements. For these measurements, the flow configurations are kept exactly the same with the exception that no tracer particles are added to the flow and that a larger FOV of $40\ \mathrm {mm} \times 80\ \mathrm {mm}$ (width $\times$ height) than with PSV for the liquid phase is investigated to obtain sufficient bubble statistics. Here, we have evaluated $500$ such image pairs for each case considered.

A machine learning procedure is used to detect and intersect overlapping bubbles in the images (Hessenkemper, Ziegenhein & Lucas 2022). The underlying neural network uses a U-Net architecture (Ronneberger, Fischer & Brox 2015) and is trained to find the intersection boundary of a bubble in front of another bubble. Using these intersections, the bubbles are split and an ellipse is fitted to the contour of each bubble. The ellipse parameters are then used to calculate the volume of an axisymmetric spheroid. Since the gas fraction is comparatively small, the number of overlapping bubbles is also small and up to now unevaluated errors connected to this are estimated to approximately $5\,\%$. Due to the contamination of the bubbles, deviations from an ellipse due to irregular (wobbling) bubble surfaces are only minor. Example images for the two different bubble sizes with fitted ellipses are given in figure 3 for all the cases. The bubble size is then calculated using the volume-equivalent bubble diameter of a spheroid as $d_p=(d_{maj}^{2}d_{min})^{1/3}$, where $d_{maj}$ and $d_{min}$ are the lengths of the major and minor axes of the fitted ellipse, respectively. For the smaller bubbles, the shape is close to a sphere with a small aspect ratio $\chi =d_{maj}/d_{min}=1.12$ and we obtain $d_p=2.7\ \mathrm {mm}$. In contrast with the larger bubbles, we have $d_p=3.8{\sim}3.9\ \mathrm {mm}$ and $\chi =1.3$. Here, ${\rm \Delta} \rho$ is the density difference between the carrier phase and the dispersed phase and $\sigma$ is the surface tensor coefficient. The range of the bubble Reynolds number, $Re_p$, in our experiments is $500{\sim}900$ based on $d_p$ and the bubble to liquid relative velocity. Note that, although no net liquid flow through a cross-section ($x_2$–$x_3$ plane) is present (i.e. the bulk velocity of the liquid is zero), for bubble columns, the flow exhibits a large-scale circulation in a time-averaged sense with upward flow in the centre and downward flow close to the walls (Mudde 2005). This results in a non-negligible upward flow in the measurement area of the present PSV (shown later in § 3), so that the relative velocity in the FOV is in general not equal to the bubble terminal rise velocity. Table 1 lists all the relevant parameters for the bubble properties.

Figure 3. Example images of the bubbles with fitted ellipses for an arbitrary instant: (a) case SmLess, (b) case SmMore, (c) case LaLess and (d) case LaMore (figures are cut from a region of the FOV for the gas phase for the corresponding case).

A shallow DoF is used for the measurements conducted to obtain bubble statistics, but it is an order of magnitude larger than that for the liquid velocity measurements due to the different optical system for these two sets of experiments. To measure the bubbles in the $x_3$-centre region, where the liquid velocity measurements were also taken, the average grey value derivative along the bubble contour is calculated. A threshold for this average grey value derivative is used so that only bubbles in the $x_3$-centre region are considered, and the depth of this $x_3$-centre region is approximately $5\ \mathrm {mm}$. Hence, the evaluated depth for the gas phase is an order of magnitude larger than for the liquid velocity, which is, however, necessary due to the size of the bubbles. The number of bubbles that fall into this region ($40\ \mathrm {mm} \times 80\ \mathrm {mm} \times 5\ \mathrm {mm}$) is between $2500$ and $5000$ summed up over the considered $500$ image pairs depending on the case being considered, and this corresponds to an average void fraction range of $0.26\,\%$ to $1.31\,\%$ (table 1). The bubble velocity was determined by calculating the distance translated by individual bubble centroids, using a standard particle tracking velocimetry algorithm that matches the closest neighbouring bubble in a small search window in the presumed direction of bubble motion (Bröder & Sommerfeld 2007).

3. Basic flow characterization

Before exploring the properties of the flow at different scales, we begin with some basic characteristics of the flow as quantified by single-point statistics of both phases. The instantaneous velocity $\tilde {\boldsymbol {u}}(\boldsymbol {x},t)$ is decomposed into a ensemble-averaged part $\boldsymbol {U}(\boldsymbol {x},t)$ and a fluctuating part $\boldsymbol {u}(\boldsymbol {x},t)$, with associated components $\tilde {u}_i=U_i+u_i$.

For all the cases investigated, it is observed that the local void fraction (figure 4a), the vertical gas/liquid mean velocity (figure 4b) and the liquid fluctuating velocity profiles (figure 4c,d) in the observation region (FOV of either liquid or gas) of the column are constant to a very good approximation, implying statistical homogeneity of the flow in the observation region. As mentioned earlier, for most cases there is a slight upwards mean flow for the liquid phase in the FOV.

Figure 4. One-point statistics along the horizontal axis of FOV for the four considered cases: (a) gas void fraction, (b) liquid/gas vertical velocity and liquid fluctuating velocity in (c) vertical and (d) horizontal components.

According to figure 4, the root-mean-square value of the vertical velocity fluctuation $u_1^{rms}$ is significantly larger than the horizontal component $u_2^{rms}$ for each case, indicating the large-scale anisotropy in the flow due to the preferential direction of the bubble motion. The ratio of the vertical to horizontal velocity fluctuation $u_1^{rms}/u_2^{rms}$ is in the range $1.4$ to $1.6$ for the four cases investigated, which is close to the previous studies of Riboux et al. (2010) and Ma et al. (2021) for bubbles with sizes in the range $1$ to $3\ \mathrm {mm}$. For both components, we find the velocity fluctuations increase with increasing $\alpha$ – from SmLess to LaMore.

Following Ma et al. (2021), we define a Reynolds number $Re_{H_2}\equiv u^{\ast } H_2/\nu$ which indicates the range of scales in the turbulent bubbly flows. Here, $u^{\ast }\equiv \sqrt {(2/3)k_{{FOV}}}$, and $k_{{FOV}}$ is the weighted turbulent kinetic energy (TKE), $k_{{FOV}}=((u_1^{rms})^{2}+2(u_2^{rms})^{2})/2$, assuming axisymmetry of the flow in the FOV about the vertical direction, and averaged over the FOV of the liquid phase. Note that for our experiments the bulk Reynolds number is zero, since there is no averaged net liquid flow when averaged over the entire flow cross-section. Moreover, we are interested in the properties of the fluctuating velocity field, and hence it is more appropriate to consider a Reynolds number based on $u^{\ast }$.

In figure 5 we plot $Re_{H_2}$ vs large-scale anisotropy ratio $u_1^{rms}/u_2^{rms}$ (averaged over the FOV of liquid). The figure shows that $Re_{H_2}$ increases in the order of SmLess, SmMore, LaLess to LaMore, which corresponds to larger bubbles and higher void fraction. Since the Reynolds number of the flow is usually understood to be related to the range of excited scales of motion in the system, this result implies that in the aforementioned sequence the range of excited scales also increases in the flow, hence, the flow becomes increasingly multiscale. Furthermore, we find that for large scales whose velocities are characterized by $u_1^{rms}$ and $u_2^{rms}$, the smaller bubbles produce more anisotropy in the flow than the larger bubbles, as indicated by a larger ratio of $u_1^{rms}/u_2^{rms}$ for the cases SmLess and SmMore. This is in very close agreement with our previous study based on DNS data of bubble-laden turbulent channel flow driven by a vertical pressure gradient (Ma et al. 2021).

Figure 5. Reynolds number, $Re_{H_2}$, plotted vs large-scale anisotropy ratio, $u_1^{rms}/u_2^{rms}$.

To quantify the PSV resolution with respect to the Kolmogorov scale $\eta$, we estimate the mean dissipation $\epsilon$ based on an algebraic relation derived by Ma, Lucas & Bragg (2020a) for BIT-dominated flows

(3.1)\begin{equation} \epsilon=S_k/(1-\alpha) , \end{equation}

where the interfacial term $S_k$ for the TKE transport equation is adopted from Ma et al. (2017) and Ma (2017)

(3.2)\begin{equation} S_{k}=\min(0.18Re_{p}^{0.23},1)\boldsymbol{F}_{D}\boldsymbol{\cdot}\boldsymbol{u}_r . \end{equation}

Here, $\boldsymbol {u_r}$ is the mean slip velocity between the bubble and the liquid, and $\boldsymbol {F}_{D}=(3/4d_p)C_D\alpha \|\boldsymbol {u}_r\|\boldsymbol {u}_r$ is the drag force on the bubbles averaged over the FOV of the liquid. We can then compute the Kolmogorov length scale $\eta \equiv (\nu ^{3}/\epsilon )^{1/4}$, as well the time scale $\tau _\eta \equiv (\nu /\epsilon )^{1/2}$. Both parameters and the ratio of $\varDelta /\eta$ (ranges between $1.8$ and $2.8$) are given in table 2 for all the cases. The current PSV spatial resolution and the size of the liquid FOV are comparable to the recent experimental study of Carter & Coletti (2017) on single-phase turbulence using standard PIV (see their small FOV/high resolution measurement).

In figure 6 we plot the PDFs of the liquid velocity fluctuations for both directions, normalized by their standard deviations. Due to the large quantity of velocity fields recorded from the experiment the PDFs are well converged with tails extending to extreme values, with values of the PDF spanning six orders of magnitude, which is much greater than previous experiments for bubbly flows (Riboux et al. 2010; Alméras et al. 2017; Lai & Socolofsky 2019). Both the vertical and horizontal velocity PDFs are in good quantitative agreement with the previous experimental studies just mentioned. While PDFs of the horizontal velocity fluctuations are symmetric and non-Gaussian, for the vertical component, the PDFs are strongly positively skewed for all the cases. This positive asymmetry originates from the wake entrainment as visualized in figure 2 for the region directly behind the in-focus bubble, which leads to a larger probability of upward fluctuations. The asymmetry of the PDFs is highest for the case SmLess with the lowest gas void fraction and gradually reduces at larger $Re_{H_2}$. In our four cases, $Re_{H_2}$ increases with increasing $\alpha$, and a similar trend was also reported in figure 5(a) of Prakash et al. (2016). Furthermore, for the PDF of the horizontal velocity component, SmLess is the furthest from Gaussian, associated with stronger large-scale intermittency than the other cases. In many aspects, those results are very different from those of homogeneous isotropic turbulence for single-phase flows, for which the PDFs of the velocity fluctuations are almost Gaussian (Gotoh, Fukayama & Nakano 2002).

Figure 6. Normalized PDFs of liquid velocity fluctuations: (a) the vertical component and (b) the horizontal component.

4. Turbulence anisotropy quantified using structure functions

While Kolmogorov's theory assumes local isotropy at the small scales of turbulent flows (Kolmogorov 1941, K41 for brevity), many experimental and DNS results reveal persistent anisotropy because the process of a return to isotropy at progressively smaller scales can be very slow (Pumir & Shraiman 1995; Shen & Warhaft 2000; Ouellette et al. 2006; Carter & Coletti 2017). The most systematic approach for characterizing anisotropy at different scales is based on the use of irreducible representations of the SO(3) group (Arad, L'vov & Procaccia 1999; Biferale & Procaccia 2005). However, this method requires information on from the complete three-dimensional flow fields, which is often not available. Due to this practical difficulty, many studies on turbulence anisotropy are based on the velocity structure function tensor, and then characterizing anisotropy based on how this quantity varies in different flow directions. The $n$th-order structure function is defined as

(4.1)\begin{equation} \boldsymbol{D}_n(\boldsymbol{r},t)\equiv\langle {\rm \Delta} \boldsymbol{u}^{n}(\boldsymbol{x},\boldsymbol{r},t) \rangle , \end{equation}

where ${\rm \Delta} \boldsymbol {u}(\boldsymbol {x},\boldsymbol {r},t)\equiv \boldsymbol {u}(\boldsymbol {x}+\boldsymbol {r},t)-\boldsymbol {u}(\boldsymbol {x},t)$ is the fluid velocity increment, and $\langle \cdot \rangle$ denotes an ensemble average. The calculation of liquid velocity increments in a bubbly flow is, however, delicate, since the liquid velocity is not defined at points occupied by a bubble. To overcome this non-continuous velocity signal challenge, we use the method proposed in Ma et al. (2017) – store only liquid velocity data along the horizontal/vertical PSV grid lines whenever the entire line is free from bubbles in the considered FOV (see the two dashed lines in figure 2 for example). Based on the $60\ 000$ velocity fields that were recorded for each case, we were able to extract data along $500\ 000$ to $1\ 500\ 000$ such vertical/horizontal lines, depending on the case. Note that Freund & Ferrante (2019) employed wavelet transforms to handle discontinuities of the velocity field at the bubble interface, and it would be interesting in future work to see how this could be used in analysing experimental data as an alternative to only extracting data along lines free of bubbles.

4.1. Second-order structure function

We first consider the second-order structure function, whose components are

(4.2)\begin{equation} D_2^{ij}(\boldsymbol{x},\boldsymbol{r},t)\equiv\langle {\rm \Delta} u_i(\boldsymbol{x},\boldsymbol{r},t){\rm \Delta} u_j(\boldsymbol{x},\boldsymbol{r},t) \rangle . \end{equation}

Hereafter, we suppress the space $\boldsymbol {x}$ and time $t$ arguments since we are considering a flow which is statistically homogeneous and stationary flows over the FOV. The PSV technique provides access to data associated with separations along two directions, namely, the vertical separation $\boldsymbol {r}=r_1\boldsymbol {e}_1 (r_1\equiv \|\boldsymbol {r}\|)$ and the horizontal separation $\boldsymbol {r}=r_2\boldsymbol {e}_2 (r_2\equiv \|\boldsymbol {r}\|)$. Hence, we are able to compute the four contributions

(4.3)$$\begin{gather} D^{L}_2(r_1)=D_2^{11}(r_1) , \end{gather}$$

(4.4)$$\begin{gather}D^{L}_2(r_2)=D_2^{22}(r_2) , \end{gather}$$

(4.5)$$\begin{gather}D^{T}_2(r_1)=D_2^{22}(r_1) , \end{gather}$$

(4.6)$$\begin{gather}D^{T}_2(r_2)=D_2^{11}(r_2) , \end{gather}$$

based on the Cartesian coordinate system depicted in figure 1. For an incompressible and isotropic flow, the following relation holds for the transverse structure function:

(4.7)\begin{equation} D^{T}_{2,iso}(r_\gamma)=D^{L}_{2}(r_\gamma)+\frac{r_\gamma}{2} \frac{\partial}{\partial{r_\gamma}}D^{L}_{2}(r_\gamma) , \end{equation}

where no summation over $\gamma$ is implied.

Figure 7 shows all the measured components of the second-order structure function, as well as $D^{T}_{2,iso}(r_\gamma )$ obtained from (4.7) for the representative case SmLess in figure 7(a,d). The results show that the values of the structure functions increase in the order SmLess, SmMore, LaLess and LaMore, which corresponds to increasing $Re_{H_2}$, gas void fraction and/or bubble Reynolds number. This relationship holds for all of the components computed and across all scales. A similar trend was also reported for all three diagonal components of the second-order structure function based on separations in the spanwise direction of a bubble-laden turbulent channel flow in our recent study (Ma et al. 2021). The second-order structure function is also related to the energy spectrum, which is often considered in analysing bubble laden flows. In figure 8 we plot the energy spectra corresponding to $u_1$ for horizontal wavenumbers, and the results show the behaviour $E_{11}\sim \kappa ^{-3}$ for all the cases, as observed in previous studies which was discussed in § 1. We also note that the $-3$ scaling spans wavelengths in the approximate range $1\times 10^{-3}$ to $5\times 10^{-3}\ {\rm m}$. If the bubble wake is approximately statistically axisymmetric due to the bubble shape, and if the $-3$ range is associated with the bubble wake, then since the upper value $5\times 10^{-3}\ {\rm m}$ is ten times smaller than the depth of the column, the wake should be only weakly affected by the walls of the column in the $x_3$ direction, for bubbles in the FOV. However, this is only an estimate, and it is possible that the finite column depth relative to the bubble size could play some role in restricting the bubble wake development, a point for exploration in future work.

Figure 7. Measured second-order transverse (a,d) and longitudinal (b,c) structure functions, with separations along the horizontal (a,b) and the vertical (c,d) directions. Note that $D^{T}_{2,iso}(r_i)$ calculated with (4.7) is shown for SmLess in (a,d).

Figure 8. One-dimensional energy spectra of the vertical component of the fluctuating liquid velocity, where the wavevector is in the horizontal direction. The two vertical dashed lines denote the wavelength $\kappa ^{-1}=d_p$ for smaller and larger bubbles, respectively.

At the large scales of a homogeneous flow both $(1/2)D^{11}_2(r_\gamma )/\langle u_1u_1\rangle$ and $(1/2)D^{22}_2(r_\gamma )/\langle u_2u_2\rangle$ approach unity. Our data for these normalized quantities are shown in figure 20 in Appendix A and reveal that, along the separation $r_2$, the quantities converge to unity when $r_2\rightarrow H_2$. However, while $(1/2)D^{22}_2(r_1\rightarrow H_1)/\langle u_2u_2\rangle$ approaches unity, $(1/2)D^{11}_2(r_1\rightarrow H_1)/\langle u_1u_1\rangle \approx 0.8$ for all the cases. This implies that the FOV is not large enough in the vertical direction to resolve the integral length scale of the flow in this direction, while it is for the horizontal direction. However, this may be simply due to the fact that $H_1< H_2$ for the FOV.

For the SmLess case, departures from $D^{11}_2(r_2)=D^{11,iso}_2(r_2)$ are larger at the large scales, reflecting the large-scale anisotropy characterized by e.g. $u^{{rms}}_1/u^{{rms}}_2$ in table 2. At for smallest values of $r_2$, $D^{11}_2(r_2)\approx D^{11,iso}_2(r_2)$. However, for separations along the vertical direction there are strong departures from $D^{22}_2(r_1)=D^{22,iso}_2(r_1)$ at all scales (figure 7d). We obtained similar results (not shown) for the other three cases. Interestingly, such behaviour was also observed in Carter & Coletti (2017) for single-phase turbulence generated by jet stirring (with zero mean flow) at a similar range of $u^{{rms}}_1/u^{{rms}}_2$.

For the second-order structure functions, one of other ways to quantify anisotropy is by the ratios of components such as $D^{L}_2(r_1)/D^{L}_2(r_2)$ and $D^{T}_2(r_2)/D^{T}_2(r_1)$, which would be equal to unity for an isotropic flow. The results in figure 9 show that both ratios monotonically decrease for decreasing separation. While the ratio of the longitudinal structure functions is very close to unity at small scales, the ratio of transverse structure functions departs more strongly from unity. By comparing the different cases, the results indicate that the smaller bubbles generate stronger anisotropy in the flow than the larger bubbles across the scales, which is in agreement with the DNS results in Ma et al. (2021) and also the results of large-scale anisotropy in § 3.

Figure 9. Ratio of longitudinal (a) and transverse (b) structure functions in different separation directions for all the cases. In (a,b) the horizontal lines indicate the value of unity and the two vertical dashed lines show $r=d_p$ for smaller and larger bubbles, respectively.

In general, the differing behaviour of the second-order structure function in the longitudinal and transverse directions shows that both need to be considered in order to fully consider the turbulence anisotropy across scales. This conclusion based on our bubble-laden flow is in agreement with that of Carter & Coletti (2017) for single-phase turbulence.

4.2. High-order structure function

Next, we consider measures of the anisotropy based on even-order structure functions up to order twelve by considering the ratios $D^{L}_n(r_1)/D^{L}_n(r_2)$ and $D^{T}_n(r_2)/D^{T}_n(r_1)$ which would be equal to unity for an isotropic flow. Using higher-order structure functions allows for a characterization of how anisotropic the large fluctuations in the system are. The results in figures 10 and 11 show that the data are statistically well converged up to order 8. We also performed a convergence check by comparing the higher-order moments computed only using half of the data, and using the full data. The results showed that for structure functions up to eighth order the two results were almost identical, indicating that our results are statistically converged due to the very large number of data points obtained in the experiments. For orders 10 and 12, the results are not fully converged. However, although there is considerable noise, the general trend of the results seems clear. The results show that the ratios $D^{L}_n(r_1)/D^{L}_n(r_2)$ and $D^{T}_n(r_2)/D^{T}_n(r_1)$ deviate increasingly strongly from unity as $n$ is increased, with the results for $n\geq 6$ reaching values very far from unity across most of the flow scales. This shows that extreme fluctuations in the flow are more anisotropic than the ‘typical’ fluctuations characterized by the $n=2$ results. Similar behaviour has also been observed for single-phase turbulence, showing that higher-order moments are the most anisotropic (Kurien & Sreenivasan 2000; Warhaft & Shen 2002). For the longitudinal directions, the results show that, just as for the $n=2$ results, the cases with smaller bubbles exhibit greater anisotropy than those with larger bubbles. However, for the transverse direction the LaLess case is by far the most anisotropic. We are not aware of a good explanation for this differing dependence of the longitudinal and transverse anisotropy on the properties of the bubbles.

Figure 10. Ratio of $n$th-even-order longitudinal structure functions in different separation directions for all the cases ($n=2, 4, 6, 8, 10, 12$). The grey horizontal line in each plot indicates the isotropic value of unity. Note that (a) presents the same value as figure 9(a), just with a different label of $y$-axis.

Figure 11. Ratio of $n$th-even-order transverse structure functions in different separation directions for all the cases ($n=2, 4, 6, 8, 10, 12$). Note that (a) presents the same value as figure 9(b), just with a different label of $y$-axis.

5. Energy transfer

We now turn to considering the mean energy transfer in the flow, by analysing the third-order structure function as well as the related inter-scale energy transfer term (Hill 2001; Alexakis & Biferale 2018)

(5.1)\begin{equation} \mathcal{F}(\boldsymbol{r})=\left.\sum_{\gamma=1}^{3}\mathcal{F}_\gamma(\boldsymbol{r})\equiv \sum_{\gamma=1}^{3} \sum_{\gamma=1}^{3}\frac{\partial}{\partial \xi_\gamma}D^{\gamma ii}_3(\boldsymbol{\xi})\right\vert_{\boldsymbol{\xi}=\boldsymbol{r}},\end{equation}

that appears in the Kármán–Howarth-type equation governing $D^{ii}_2(\boldsymbol {r})$, where $D^{\gamma ii}_3(\boldsymbol {r})\equiv \langle {\rm \Delta} {u_\gamma }(\boldsymbol {r}) {\rm \Delta} u_i(\boldsymbol {r}){\rm \Delta} u_i(\boldsymbol {r})\rangle$. Note that the divergence term in (5.1) has been written using the variable $\boldsymbol {\xi }$ to emphasize that this term involves the divergence evaluated at $\boldsymbol {r}$.

The sign of $\mathcal {F}$ indicates the direction of the nonlinear energy transfer between the scales of the flow, and in three-dimensional single-phase turbulence, we have $\mathcal {F}<0$ in the inertial and dissipation ranges, corresponding to a downscale transfer of energy from large to small scales on average. With our two-dimensional PSV data we are not able to construct the full quantity $\mathcal {F}$, however, we are able to consider some of the contributions $\mathcal {F}_\gamma$.

In figure 12 we plot the components $\mathcal {F}_1(r_1)=d(D^{111}_3+D^{122}_3)/dr_1$ and $\mathcal {F}_2(r_2)=d(D^{222}_3+D^{211}_3)/dr_2$. The peaks in the magnitudes of both $\mathcal {F}_1$ and $\mathcal {F}_2$ increase in the sequence SmLess, SmMore, LaLess and LaMore, which corresponds to higher $Re_{H_2}$, $\alpha _p$ and/or $Re_p$. This is in part simply due to the fact that the kinetic energy in the flow increases in this sequence also, and with more kinetic energy at the large scales, there is more available to be passed downscale. Concerning $\mathcal {F}_1$, the results show that at scales $r>O(d_p)$, $\mathcal {F}_1<0$ corresponding to energy moving downscale, while at scales $r< O(d_p)$, $\mathcal {F}_1>0$, corresponding to energy moving upscale. In contrast to $\mathcal {F}_1$, the contribution $\mathcal {F}_2$ shows an opposite trend, with $\mathcal {F}_2<0$ for at smaller scales indicating a downscale energy transfer. The magnitude of $\mathcal {F}_2$ is, however, in general much larger than $\mathcal {F}_1$, implying a higher energy transfer rate in the horizontal direction of the flow. At the low Reynolds number of the flow, there is no inertial range over which either $\mathcal {F}_1$ or $\mathcal {F}_2$ are constant, indicating that the nonlinear energy transfer is not in the form of a constant flux cascade.

Figure 12. Contribution $\mathcal {F}_1$ (a) and $\mathcal {F}_2$ (b) to the nonlinear energy transfer term for all the cases. In (a,b) the two vertical dashed lines show $r=d_p$ for smaller and larger bubbles, respectively.

Our previous study (Ma et al. 2021) explored bubble-laden turbulent channel flows and also observed that there is a downscale energy transfer for $\mathcal {F}_3(r_3)$ (where in that paper the 3 direction corresponded to the spanwise direction of the channel flow). However, in that study we could not compute either $\mathcal {F}_1$ (streamwise) or $\mathcal {F}_2$ (wall normal) due to the one-dimensional dataset being used. Although the type of flow presently under consideration is considerably different from a pressure driven bubble-laden turbulent channel flow, the results nevertheless indicate that the nonlinear energy transfer behaviour may be very different in the different directions of the flow.

These results therefore show that the direction of the energy flux can differ from one direction to another, with $\mathcal {F}_1$ revealing an upscale energy transfer at smaller scales, while $\mathcal {F}_2$ reveals the opposite. Therefore, although the total energy transfer may be downscale in bubbly turbulent flows (Pandey et al. 2020), some directions of the flow can exhibit an upscale transfer.

The energy transfer term $\mathcal {F}$ depends on the third-order structure function, and its behaviour may be further understood by considering the associated skewness of the velocity increments. Figure 13 shows $D^{L}_3/(D^{L}_2)^{3/2}$ (the skewness of the longitudinal velocity increment) for the two separations $r_1$ and $r_2$. The results show that when $D^{L}_3(r_2)$ is plotted in this normalized form, the results for all of the cases almost collapse, and $D^{L}_3(r_2)/(D^{L}_2(r_2))^{3/2}$ is not far from the value $-0.4$ that occurs for single-phase isotropic turbulence (Davidson 2015) for $r\rightarrow 0$. The negative longitudinal skewness along the separation perpendicular to the bubble motion direction was also found in Ma et al. (2021) for some of the bubble-laden cases considered, however, for other cases this skewness was actually positive. This is very different from our results where all the cases considered have very similar negative values of the skewness. It is not clear as to whether this is due to the different kinds of flows being studied here and in Ma et al. (2021), or perhaps also because of differences in the parameter space of the bubbles being explored. For vertical separations $r_1$, $D^{L}_3(r_1)/(D^{L}_2(r_1))^{3/2}$ also almost collapses for all the cases but is positive at the small scales with a similar magnitude $\approx 0.4$ to the results for $D^{L}_3(r_2)/(D^{L}_2(r_2))^{3/2}$. This again indicates strong anisotropy for the odd-order structure functions, just as for the even-order structure functions discussed in § 4.

Figure 13. Normalized third-order longitudinal structure functions (skewness), with separations along the horizontal (a) and the vertical (b) directions. The horizontal lines in (a,b) indicate the value of zero.

Similar plots are obtained for the transverse structure functions $D^{T}_3(r_2)$ and $D^{T}_3(r_1)$. The corresponding skewness values were found to be close to zero in both separations directions for all the cases (not shown here).

6. Extreme fluctuations in the flow

Having considered the flow anisotropy and nonlinear energy transfer, we now turn to consider the extreme fluctuations in the flow. Single-phase turbulence at high Reynolds number exhibits fluctuations in space and time of its small-scale quantities that are orders of magnitude larger than the average values, a phenomenon referred to as small-scale intermittency (Frisch 1995; Sreenivasan & Antonia 1997). Such extreme fluctuations are critical for many processes in both nature and engineering (Sapsis 2021), and their investigation continues to be an area of active research (Yeung, Zhai & Sreenivasan 2015; Buaria et al. 2019). However, very few studies have explored the topic of intermittency in bubble-laden turbulent flows (see the related works introduced in § 1). In what follows, we focus on extreme fluctuations of the velocity increments ${\rm \Delta} \boldsymbol {u}$ in the bubble-laden flow we are considering.

6.1. Probability density functions

Figure 14 shows the instantaneous velocity increments ${\rm \Delta} u_i(r_2=\varDelta )$, normalized by their standard deviations $\sigma _{{\rm \Delta} u_i}$ along a line (without bubbles) in the direction of $x_2$ within the FOV. The plots reveal intermittent fluctuations of ${\rm \Delta} u_i(r_2=\varDelta )$, with relatively small regions where ${\rm \Delta} u_i(r_2=\varDelta )/\sigma _{{\rm \Delta} u_i}$ takes on large values, between which there are relatively small fluctuations of ${\rm \Delta} u_i(r_2=\varDelta )/\sigma _{{\rm \Delta} u_i}$. We will return later to considering how these intermittent large fluctuations of ${\rm \Delta} u_i(r_2=\varDelta )/\sigma _{{\rm \Delta} u_i}$ are associated with the bubble wakes in the flow.

Figure 14. Instantaneous velocity increments normalized by standard deviations: (a) ${\rm \Delta} u_1(r_2=\varDelta )/\sigma _{{\rm \Delta} u_1}$ and (b) ${\rm \Delta} u_2(r_2=\varDelta )/\sigma _{{\rm \Delta} u_2}$, along an arbitrary horizontal line for an instant, when the line is free of bubbles for each considered case. In (a,b) the curves for SmMore, LaLess and LaMore are shifted upward by $7, 14$ and $21$, respectively, for a better visual representation.

To statistically quantify the extreme events in the flow, the PDFs of the velocity increments can be considered, and these are shown in figure 15, normalized by the standard deviation, for separations equal to one PSV grid $(r=\varDelta )$. All of the results in figure 15 show that the velocity increments have strongly non-Gaussian PDFs, consistent with the observations based on figure 14. While the PDFs of the transverse velocity increments in panels (a,d) are almost symmetric, the PDFs of the longitudinal velocity increment are negatively skewed for horizontal separations $r_2$, panel (b), and positively skewed for vertical separations $r_1$, panel (c), consistent with the skewness results in figure 13. The results also show that the PDFs become increasingly non-Gaussian in the order of LaMore, LaLess, SmMore, SmLess, which corresponds to the order of decreasing Reynolds number. Therefore, surprisingly, while non-Gaussianity of the PDFs of the velocity increments becomes stronger as Reynolds number is increased in single-phase turbulence (Frisch 1995), the opposite occurs for the bubble-laden turbulent flow considered here.

Figure 15. Normalized PDFs of the transverse (a,d) and longitudinal (b,c) velocity increments, with the separation $r=\varDelta$ along the horizontal (a,b) and the vertical (c,d) directions for all four cases.

In turbulent flows, intermittency is characterized not simply by non-Gaussianity of the small-scale flow quantities, but also by the non-Gaussianity depending on scale (as well as Reynolds number) (Frisch 1995). To consider intermittency in the flow we therefore consider how the normalized PDFs of the velocity increments change with scale. Figure 16 shows the PDFs of the transverse (a,d) and longitudinal (b,c) velocity increments for various separation distances $r_\gamma =\varDelta, 10\varDelta,20\varDelta, 40\varDelta, H_\gamma$ and for the representative case SmMore. The general behaviour is very similar to that observed for single-phase turbulence (Ishihara, Gotoh & Kaneda 2009): as the scale is decreased, the PDFs of both the transverse (a,d) and longitudinal (b,c) velocity increments deviate more and more from a Gaussian distribution with the tails becoming increasingly heavy, and the PDFs becoming increasingly peaked. For the PDFs of longitudinal velocity increments are negatively skewed for ${\rm \Delta} u_2(r_2)$ and positively skewed for ${\rm \Delta} u_1(r_1)$, but the skewness becomes progressively weaker as the scale increases, and for the largest scale captured in the FOV, i.e. $r_\gamma =H_\gamma$, we find almost symmetric PDFs.

Figure 16. Normalized PDFs of the transverse (a,d) and longitudinal (b,c) velocity increments in the case SmMore for various separation distances $r_i=\varDelta, 10\varDelta,20\varDelta, 40\varDelta, H_i$, where $\varDelta$ is one PSV grid. (a,b) Separations along the horizontal direction and (c,d) separations along the vertical direction.

We note that, while the PDFs of the velocity increments should asymptote to the single-point PDFs of the velocity in the limit $r\rightarrow \infty$, our data for the PDF of ${\rm \Delta} u_1(r_\gamma =H_\gamma )$ are very different (compare the left tails of the PDFs) from single-point PDF of $u_1$ shown in figure 6(a). This is probably simply due to the fact that, as already noted earlier, the size of the FOV is not large enough for this asymptotic regime to be observed. The aforementioned PDFs are, however, so different that the flow behaviour must undergo an interesting transition where negative fluctuations of the vertical velocity begin to be suppressed (but not positive fluctuations) as the scale continues to be increased. This point will be explored in a future work by employing a larger FOV to observe this asymptotic regime.

6.2. Flatness of velocity increment

The intermittency and non-Gaussianity can be quantified using the flatness measures defined as $D^{T}_4(r_\gamma )/(D^{T}_2(r_\gamma ))^{2}$ (figure 17a,d) and $D^{L}_4(r_\gamma )/(D^{L}_2(r_\gamma ))^{2}$ (figure 17b,c). The first observation from the results in figure 17 is that all four bubble-laden cases show a rather similar behaviour as $r$ increases, slowly reducing from values of up to 40 towards the Gaussian value of 3. However, there is still considerable deviation from 3 even at the largest scale $r_\gamma /H_\gamma =1$, consistent with the PDF results for ${\rm \Delta} \boldsymbol {u}(r_\gamma =H_\gamma )$ in figure 16. We also find that, for both the transverse or longitudinal flatness measures and along either separation direction, the flatness is largest for the case with smaller $Re_{H_2}$, reflecting the same trend as the PDFs of the velocity increments at the single scale $\varDelta$ (figure 15). The values and behaviour of the flatness are similar to those observed in Ma et al. (2021) for the three bubble-laden cases with $\alpha _p=2.14\,\%$, although for the present flow the flatness values are larger than those in Ma et al. (2021). For separations $r_2$, the transverse flatness values are larger than the longitudinal ones, similar to that observed in isotropic turbulence (Li & Meneveau 2005; Ishihara et al. 2009). In contrast, for separations $r_1$, the transverse flatness values are slightly smaller than the longitudinal ones.

Figure 17. Normalized fourth-order transverse (a,d) and longitudinal (b,c) structure functions, corresponding to the flatness of the velocity increments. The horizontal lines in (a–d) indicate the Gaussian value of $3$ for the flatness.

6.3. Flow structures associated with extreme fluctuations

We now use flow visualization to understand which regions of the flow are associated with the extreme velocity increments. Figure 18(c,d) gives two examples of such flow visualization extracted from a randomly chosen (but representative) FOV result for ${\rm \Delta} u_1(r_2=\varDelta )/\sigma _{u_1}$ for the cases SmMore and LaMore, respectively. Each grid in the figure corresponds to an interrogation window of the present PSV. For convenience, we also plot in (a,b) the corresponding original velocity vector fields and denote in each case the in-focus bubbles. The images show that for both cases the regions where the velocity increment takes on extreme values are located in the outer-most region of the individual bubble wakes, while the fluctuations inside the wakes are relatively weak (see c,d). At a result of this it appears that the small-scale intermittency in the present flow already characterized is quite different from the small-scale intermittency that occurs in single-phase turbulence that arises due to nonlinear self-amplification of the velocity gradients (Tsinober 2001; Li & Meneveau 2005; Buaria et al. 2019; Buaria, Bodenschatz & Pumir 2020). In particular, the intermittency in the present flow is more similar to what is often called external intermittency (Townsend 1949), which arises near the interface between two regions of the flow with significantly differing behaviour, e.g. such as at the turbulent/non-turbulent interface in turbulent jets (da Silva et al. 2014; Gauding et al. 2021). In the bubble-laden flow we are considering, the intermittency arises at the interface between the edge of the bubble wake and the surrounding relatively quiescent flow.

Figure 18. Snapshot of the original velocity vector (a,b) and the intensity distributions of normalized velocity increment ${\rm \Delta} u_1(r_2=\varDelta )/\sigma _{{\rm \Delta} u_1}$ (c,d). Here, (a,c) are from the same instant based on the SmMore case and (b,d) are from the same instant based on the LaMore case. The in-focus bubbles are denoted in (a,b).

One further observation on figure 18(c,d) is that the extreme values are larger for the smaller bubble wakes than the larger bubble wakes (note the different legend scales in c,d). This point is made clearer in figure 19(a,b) where instead $|{\rm \Delta} u_1(r_2=\varDelta ) |\geq 8\sigma _{{\rm \Delta} u_1}$ is plotted. Looking at the in-focus bubble on the right-hand side for both cases (where a more complete wake structure is visible than for the left in-focus bubble), we find that the regions of extreme values are larger for the SmMore case. By increasing the threshold to $|{\rm \Delta} u_1(r_2=\varDelta ) |\geq 16\sigma _{{\rm \Delta} u_1}$ (figure 19c,d), we see that there are still some locations with extreme values satisfying this threshold for the SmMore, while there are none for the LaMore case. Although this observation is only for one snapshot, we have conducted the same analysis for many different snapshots and virtually all of them show the same behaviour. These observations are consistent with the PDF results in figure 15 that show that cases with smaller bubbles (lower $Re_{H_2}$) have heavier tails than the cases with larger bubbles (higher $Re_{H_2}$).

Figure 19. Snapshots indicating the regions with extreme values, where $|{\rm \Delta} u_1(r_2=\varDelta ) |\geq 8\sigma _{{\rm \Delta} u_1}$ (a,b) and $|{\rm \Delta} u_1(r_2=\varDelta ) |\geq 16\sigma _{{\rm \Delta} u_1}$ (c,d) highlighted in red. Panels (a,c) are from the same instant from SmMore and (b,d) are from the same instant from LaMore. The in-focus bubbles are denoted in each panel.

An explanation for why the intermittency is largest for the cases with smaller bubbles and lower volume fractions was given in Ma et al. (2021). When the volume fraction is not too large (but still large enough for the bubbles to make a statistically significant contribution to the flow properties), there are relatively few regions in the flow where turbulence is produced due to the bubble wakes, meaning that turbulence in the flow is very patchy and therefore intermittent. As the volume fraction increases, the regions of the flow influenced by the bubble wakes becomes less rare, and hence intermittency in the flow reduces. This then explains why the SmLess case is the most intermittent, and the LaMore case the least intermittent. Moreover, since the flow we are considering has no background turbulence (i.e. the only turbulence in the flow is the result of the bubble wakes), then this may explain why the intermittency we see is higher than that in Ma et al. (2021) where the flow considered had background turbulence, so that turbulent activity in the flow was not restricted to regions close to the bubbles.

7. Conclusions

In this paper we have presented an analysis of the multiscale properties of a bubble-laden turbulent flow, based on experimental data of a flow in a vertical column with bubble swarms rising in water. The experiment takes advantage of a recently developed PSV technique for particle/bubble flows and provides the first comprehensive dataset for computing multipoint measurements (without having to invoke Taylor's hypothesis, in contrast to, e.g. Rensen et al. 2005) of flows laden with finite-sized bubbles. By acquiring a very large number of snapshots of the velocity field we are able to compute structure functions of up to twelfth order, and PDFs that resolve the heavy tails associated with extreme fluctuations in the flow.

The results show that the level of anisotropy in the flow produced by the rising bubbles is strong in general, and not negligible at any scale in the flow. Moreover, the results show that (i) the differing behaviour of the second-order longitudinal and transverse structure functions when measured for separations in different directions shows that both velocity components and separation directions of the two-dimensional data need to be considered in order to fully characterize the anisotropy of the flow; (ii) the bubble size and void fraction are both important parameters determining the amount of anisotropy in the flow; (iii) higher-order structure functions reveal greater anisotropy across the scales of the flow, such that extreme fluctuations in the flow are the most anisotropic. Since the PSV measurements capture two-dimensional data, we were able to consider the energy transfer between scales for two separation directions in the flow. The results revealed a downscale energy transfer on average for horizontal separations, but an upscale energy transfer on average in the vertical direction. However, the horizontal energy transfer was much stronger than that in the vertical direction.

We also investigated extreme events in the flow by considering the normalized PDFs of the velocity increments in the flow. The results showed that the probability of extreme fluctuations increases with decreasing scale, just as in single-phase turbulence. However, the results also showed that for a given scale the probability of extreme events decreased with increasing Reynolds number, contrary to what occurs in single-phase turbulence. To explore the origin of these extreme fluctuations in the bubble-laden flows, we visualized regions of extreme small-scale velocity increments in the FOV and observed that they are typically located at the boundary of the wakes produced by the bubbles. For the cases with smaller bubbles and lower void fractions, vast regions outside of the bubbles wakes exhibit weak fluctuations, and so this combined with the extreme fluctuations at the bubble wake boundaries leads to strong intermittency. For larger bubbles which produce larger flow Reynolds numbers, and with larger void fractions, the wake regions become less rare in the flow and hence the flow is less intermittent than the former case, even though the Reynolds number is higher. Furthermore, the extreme values were also observed to reach larger values (compared with the standard deviation) for the smaller bubbles, which again causes the smaller bubble cases to exhibit greater intermittency than the larger bubble cases, in addition to the effect arising from the fraction of the flow modified by the bubble wakes.

Finally, it is worth pointing out two limitations with the present PSV data. First, the FOV is not large enough to resolve the integral length scales of the flow, and therefore, as discussed in § 6.1, we were not able to fully resolve the transition in the properties of the bubble-laden flow as the scale approaches the integral length scale. Second, the bubbles considered here all have approximately fixed shape, mimicking in nature small bubbles in a contaminated flow. We are currently performing new experiments with deformable bubbles and differing FOVs in order to address these two points and will report the results in future work.