Scale-dependent anisotropy, energy transfer and intermittency in bubble-laden turbulent flows

Data from Direct Numerical Simulations of disperse bubbly flows in a vertical channel are used to study the effect of the bubbles on the carrier-phase turbulence. A new method is developed, based on the barycentric map approach, that allows to quantify the anisotropy and componentiality of the flow at any scale. Using this the bubbles are found to significantly enhance flow anisotropy at all scales compared with the unladen case, and for some bubble cases, very strong anisotropy persists down to the smallest flow scales. The strongest anisotropy observed was for the cases involving small bubbles. Concerning the inter-scale energy transfer, our results indicate that for the bubble-laden cases, the energy transfer is from large to small scales, just as for the unladen case. However, there is evidence of an upscale transfer when considering the transfer of energy associated with particular components of the velocity field. Although the direction of the energy transfer is the same with and without the bubbles, the transfer is much stronger for the bubble-laden cases, suggesting that the bubbles play a strong role in enhancing the activity of the nonlinear term in the flow. The normalized forms of the fourth and sixth-order structure functions are also considered, and reveal that the introduction of bubbles into the flow strongly enhances intermittency in the dissipation range, but suppresses it at larger scales. This strong enhancement of the dissipation scale intermittency has significant implications for understanding how the bubbles might modify the mixing properties of turbulent flows.


Introduction
Turbulence and multiphase flows are two of the most challenging topics in fluid mechanics and when combined they pose a formidable challenge, even in the dilute dispersed regime (Balachandar & Eaton 2010). The focus here is on liquid flows laden with disperse bubbles, which can be particularly challenging since the bubbles can strongly alter the liquid phase turbulence (Mudde 2005; Lohse 2018; Elghobashi 2019). In particular, the bubbles can modify the turbulence due to production effects arising from the bubble wakes (Riboux et al. 2010;Lai & Socolofsky 2019), enhanced local turbulent kinetic energy dissipation rates in the vicinity of the bubble surfaces (Santarelli et al. 2016;Masuk et al. 2021), and modulation of the liquid mean velocity profile due to interphase momentum transfer, resulting in an alteration of shear-induced turbulence (Lu & Tryggvason 2013;du Cluzeau et al. 2019;Cifani et al. 2020). Mathai et al. (2020) highlighted particular ways in which the classical scenario for single-phase turbulence, based on single-point statistical analysis, is modified due to the bubbles moving relative to the fluid. Turbulence arising from this relative motion is often referred to as bubble-induced turbulence (BIT) and its effects can be captured in the Reynolds-averaged Navier-Stokes (RANS) modelling framework through the inclusion of additional source terms in the relevant transport equations (Fox 2014;Joshi & Nandakumar 2015;Ma et al. 2020b).
While significant progress has been made in understanding and characterizing how the bubbles influence the single-point turbulence statistics of the liquid phase, less attention has been paid to the influence of the bubbles on the multiscale/multipoint flow statistics. Those that have considered this aspect have only focused on the kinetic energy spectrum of the liquid velocity fluctuations, with a key observation being that in BIT dominated flows, a power-law behavior for the energy spectrum arises with exponent −3 in both the wavenumber or frequency domain (Lance & Bataille 1991;Roghair et al. 2011;Mendez-Diaz et al. 2013;Ma et al. 2017). This behavior was also reported in Pandey et al. (2020), who investigated the energy budget equations in wavenumber space and explained the −3 slope as arising due to a balance between kinetic energy production due to the bubbles and viscous dissipation. In their analysis they evaluated for the first time (to our knowledge) the nonlinear scale-to-scale energy flux term for bubbly flows, showing a forward (downscale) energy cascade, just as also occurs for single-phase turbulence in three dimensions. An issue with their analysis, however, is that they included the values of the flow at grid points occupied by the bubbles when evaluating the statistics of the carrier phase, thereby contaminating the fluid statistics.
There have been very few studies exploring the multiscale properties of bubble-laden turbulent flows in physical space, e.g. using structure function analysis. Rensen et al. (2005) performed hot-film anemometry measurements in the Twente water tunnel and computed the longitudinal second-and fourth-order structure function with the aid of Taylor's hypothesis. They found an increase of the second-order structure function for the two-phase case compared with the singlephase case under the same bulk Reynolds number, and that this increase was more pronounced at the small scales than the large scales. Their fourth-order structure function results revealed an increase of the intermittency at the small scales of the flow when the flow contained bubbles, even for a relatively low gas void fraction (0.5%). Similar behavior was also observed in Biferale et al. (2012) when comparing the small-scale properties of boiling and non-boiling convective turbulent flows.
An important aspect yet to be quantified is how the bubbles affect the anisotropy of the flow at different scales. For single-phase turbulence, the energy containing scales in many flows such as those with shear, rotation, and buoyancy are anisotropic (Biferale & Procaccia 2005). Phenomenological theories of turbulence predict a return-to-isotropy at small enough scales (Kolmogorov 1941b;Frisch 1995;Sreenivasan & Antonia 1997). However, measurements have revealed persistent small-scale anisotropy (Pumir & Shraiman 1995;Shen & Warhaft 2000;Ouellette et al. 2006;Carter & Coletti 2017). In contrast to single-phase flow, where energy is often injected into the flow at large scales, bubbles can inject energy into the flow at the scale of their size, which usually corresponds to the small-scales of the turbulence. Since the bubbles have a preferential direction of motion due to buoyancy, this could lead to the injection of strong anisotropy into the flow at the small scales, leading to strong departures from the behavior of the single-phase case. The study of Pandey et al. (2020) was based on Fourier-space analysis with averaging over spherical shells in wavevector space, and so did not permit them to explore the anisotropy of the flow at different scales.
Another important point is that in Pandey et al. (2020) the flow had no background turbulence (i.e. all the turbulence was generated by the bubbles), and hence it was not possible to consider how the bubbles modify the turbulence compared with the single-phase case. In order to more fully understand how the bubbles modify the properties of the turbulence, it is desirable to consider a configuration in which the unladen flow is already turbulence, and then one can explore how the bubbles modify the properties of the turbulence when they are introduced.
In the present work we seek to advance the understanding of the properties of bubble-laden turbulent flows across its range of scales. To do this, data from Direct Numerical Simulation (DNS) of finite-size bubbles in a turbulent channel flow with a prescribed bulk Reynolds number is utilized, for different bubble sizes and for mono and bidisperse cases. A new method is developed based on the barycentric map (Banerjee et al. 2007), and this is used together with the DNS data to quantify the anisotropy of the bubble-laden turbulence flow across the range of its scales. By computing the structure functions of various orders, the direction-dependent liquid velocity fluctuations are also explored at different scales, as is the scale-to-scale energy transfer, and intermittency. These results provide new insights into the properties of bubble-laden turbulent flows, and how they differ from the single-phase counterpart at different scales in the flow.

Database
The DNS data we use is from the studies (Santarelli & Fröhlich 2015, 2016 that simulated the motion of many thousands of bubbles at low Eötvös number in a vertical turbulent channel flow. The bubbles are handled using an Immersed Boundary Method, and are modeled as rigid spherical objects with a no-slip condition enforced at their surface, representing the behaviour of air bubbles rising in contaminated water. Compared to other simulations of this type (see the related references in Mathai et al. 2020), these simulations are substantially closer to applications in that they involve a turbulent background flow, contaminated fluid, realistic density ratio, higher bubble Reynolds number, a much larger domain, and a much larger number of bubbles.
As shown in figure 1, the vertical flow takes place between two flat walls separated by the distance , and the size of the computational domain is × × = 4.41 × ×2.21 . Here, denotes the streamwise position coordinate, the wall-normal coordinate, and the spanwise coordinate, and the corresponding unit basis vectors are , , , respectively. The numerical grid employed has the same spacing Δ = /232 in all directions, resulting in 1024 × 232 × 512 grid points in the , , and directions, respectively. A no-slip condition was applied at the walls, and periodic boundary conditions were applied in the and directions. The gravitational force acts in the direction − , and the bulk velocity was kept constant by instantaneously adjusting a volume force, equivalent to a pressure gradient, thus imposing a desired bulk Reynolds number = / , where is the kinematic viscosity of the fluid. The DNS were all conducted with = 5263. The data used in this work were obtained for three monodisperse cases (SmMany, SmFew, LaMany) and one bi-disperse case labelled BiDisp, of the same void fraction as SmMany and LaMany with half the void fraction consisting of smaller bubbles and the other half of larger bubbles. Additionally, a single-phase simulation labelled Unladen was performed under the same conditions for comparison. Table 1 provides an overview of all cases with the corresponding labels.

Data processing
A standard way to analyze the multiscale properties of turbulence is to use the fluid velocity increments Δ ( , , ) ≡ ( + , )− ( , ), where ≡ − is the fluctuating fluid velocity, is the separation vector, and · denotes an ensemble average (estimated using appropriate space and time averages). The calculation of velocity increments in a bubble-laden flow is, however, delicate, since the phase boundaries can interrupt the fluid flow signal. To overcome this noncontinuous velocity signal challenge, different methods have been used in the literature, such as smoothing the discontinuities by a Gauss function (Lance & Bataille 1991); considering only intervals between bubbles where the velocity signal is continuous (Martínez et al. 2010;Mendez-Diaz et al. 2013;Roghair et al. 2011); and measuring the wake behind a rising swarm of bubbles, where there are no bubbles (Riboux et al. 2010). Here, we use a method ideally suited for interface-resolved DNS of disperse flows proposed in our previous study (Ma et al. 2017). In this method, the fluid velocity is recorded along grid lines in the spanwise direction whose wallnormal location lies within the centre region 0.474 < < 0.526 (this width corresponds to the smaller bubble diameter, see figure. 1), which is a sufficiently thin region for these lines to be considered statistically equivalent. For each line, data was recorded whenever the entire line was free from bubbles, and we recorded 1, 000, 000 instances of this for each case. With this method, we cannot compute Δ ( , , ) for arbitrary , but can compute Δ ( , 3 3 , ), allowing us to perform an extensive investigation into the scale-dependent properties of bubble-laden turbulent flows.

Reynolds number in bubble-laden turbulent flows
The range of scales in a turbulent flow is governed by its Reynolds number, and therefore before analyzing the properties of the bubble-laden turbulent channel flows at different scales we first consider the Reynolds numbers. As mentioned in §2, the bulk Reynolds number = / was kept fixed at = 5263. However, since we are interested in the properties of the fluctuating component of the velocity field, it is more informative for our purposes to consider a Reynolds number based on the fluctuating velocity field. To that end we consider the Reynolds number ≡ * / , where * ≡ √︁ (2/3) , and is the turbulent kinetic energy (TKE) ≡ (1/2) evaluated at the channel centre. In figure 2 we plot versus 2 , where 2 = is the second invariant of the Reynoldsstress anisotropy tensor (evaluated at the channel centre), = / −(2/3) , that quantifies the magnitude of the large-scale anisotropy in the flow. The results show that varies significantly across the cases. The SmFew case, is only slightly larger than the unladen case as the bubble size is small and the bulk void fraction is quite low. However, as and are increased, increases significantly, implying that as and are increased (at least over the range we consider), the range of excited scales in the flow also increases, and hence the flow becomes increasingly multiscale. The increase cannot continue indefinitely, however, since when and become sufficiently large, the problem becomes analogous to flow through a porous medium, for which the flow Reynolds number cannot be very large. This is also related to the fact that in such a regime, is no longer the relevant lengthscale in the flow Reynolds number, but rather the inter-bubble distance becomes the appropriate length scale.
In the considered range of bubble/channel parameters, figure 2 suggests a approximately linear relation between and 2 for the bubble-laden cases, with larger 2 corresponding to smaller . This behaviour may be understood in terms of the model for BIT from Ma et al. (2020a), according to which where is the interfacial term appearing in the TKE equation, ≡ with the mean slip velocity between the bubble and the fluid at the channel centre, 2 = 1.92, and is the void fraction in the channel centre. The algebraic expression in (3.1) was derived based on the fact that in BIT-dominated flows the dissipation term is balanced by the interfacial term in the TKE transport equation at the channel centre. We refer the reader to Ma et al. (2020a) for the detailed derivation.
Based on (3.1) we can re-express the velocity scale * by introducing the interfacial term adopted from Ma et al.
is the drag force on the bubbles at the channel centre. Using these results we obtain * = where = 1.12 if 0.18 0.23 < 1, and = 1 otherwise. This finally leads to In figure 3 we plot normalized by the factor √︁ /(1 − ) in order to test the prediction in (3.4). The results show that the relation ∝ √︁ /(1 − ) describes the data very well except for the SmFew case. This is due to the fact that for this case the void fraction is low (0.29%), so the influence of the background flow (ignored in (3.1)) is not negligible. Further data across a wider range of the parameter space is needed to more comprehensively validate (3.4) and to discover the extent of the range over which it holds.

Multiscale anisotropy and second-order structure function
A systematic approach for analyzing the multiscale anisotropy of a turbulent flow is to use the irreducible representations of the SO(3) group, which consists of projecting the multipoint turbulent correlation functions onto the space of spherical harmonics (Arad et al. 1999;Biferale & Toschi 2001;Biferale & Procaccia 2005). However, such an analysis requires information on the full three dimensional flow field, something that is usually not obtainable from experiments. Furthermore, as discussed in §2, phase boundaries interrupt the flow field in multiphase flows, introducing further challenges in applying this method.
Due to these challenges, many investigations on anisotropy in turbulent flows focus directly on the structure function tensor (which are essentially moments of the velocity increments), comparing the longitudinal and transverse components in order to discern the level of anisotropy in the flow. In single-phase flows, the main focus was on scrutinizing the postulate of local isotropy and its implications (Kolmogorov 1941b, K41 for brevity). In general, experiments and numerical simulations do not strictly confirm the convergence toward isotropy predicted within K41 theory as a function of the scale (Kurien & Sreenivasan 2000;Shen & Warhaft 2002;Dhruva et al. 1997). In the context of the single-point Reynolds stresses, the Lumley triangle (Lumley & Newman 1977) provided a powerful way to quantify and visualize anisotropy in the flow. It would be desirable to have something analogous to this for the structure functions, which would then provide a way to quantify and visualize anisotropy in the flow at different scales.

Second-order structure function and its anisotropy
Consider the second-order structure function (4.1) Hereafter we will suppress the time argument since we are focusing on statistically stationary flows. The Cartesian coordinate system chosen is depicted in figure 1, and as discussed in §2, our DNS data only allow us to compute the velocity increments for separations in the spanwise direction, i.e. = 3 3 (and ≡ = 3 ). In this case the longitudinal structure function is ( 3 ) = 33 ( 3 ), and for an incompressible, isotropic flow we would have  (2016) for the corresponding Reynolds normal stresses in the channel centre. Compared to the unladen case, the enhancement of the structure function level is in the sequence SmFew, SmMany, BiDisp to LaMany, which corresponds to increasing averaged bubble Reynolds number and/or gas void fraction. This holds for all three directions across all scales from around one channel width to the smallest dissipative scales, and shows that the bubbles modify fluctuations in the flow at scales both larger and smaller than the bubble length scale . Moreover, the increase of the structure functions in the bubble-laden cases is more pronounced at the smaller scales than the larger scales. This behavior is in close agreement with the experimental results in Rensen et al. (2005) which used Taylor's hypothesis to construct the results. We also note that although 22 and 33 are very similar for the Unladen and the SmFew case at the large scales, they are significantly different at the smaller scales, with the bubbles significantly enhancing the small scale fluctuations in the flow.
For the unladen case, departures from 11 = 11 and 22 = 22 are not too strong, with stronger departures at the larger scales, consistent with the Reynolds stress behaviour in the channel centre 1 1 / 2 2 ≈ 1.47. As is decreased the flow becomes more isotropic, although as we shall show later, the small scales do not actually reach an isotropic state. In contrast, for the bubble-laden cases, there are strong departures from isotropy at all scales. For example, for the SmFew case, 11 shows significant deviations from 11 at all scales, and for 22 the deviations from the isotropic form 22 actually become stronger as one goes to scales smaller than the bubble length scale ( ≈ 0.05 for this case). The concept of return to isotropy can therefore be strongly violated for bubbly turbulent flows. This is not surprising however, since the mean velocity of the bubbles is unidirectional, and the fact that means that the bubbles can directly inject fluctuations into the flow at the small scales, in contrast to the unladen case where the fluctuations are injected into the flow at large scales due to the mean shearing of the flow. Figure 5 shows the ratios of the different structure function components in order to see more clearly the anisotropy of the flow. For an isotropic system we would have 11 / 22 = 1 at all scales, and for the unladen case 11 / 22 is quite close to 1, especially at smaller scales, while for the bubble-laden cases 11 / 22 deviates strongly from 1, reaching values (10). Except for the SmFew case, the behaviour of 11 / 22 in the bubble-laden cases monotonically decreases with decreasing , indicating a tendency towards return to isotropy, although the isotropic state is far from achieved. In isotropic turbulence, 11 / 33 → 1 and 22 / 33 → 1 for ( ), while 11 / 33 → 2 and 22 / 33 → 2 for / → 0 (Pope 2000). While deviations from these are not too strong for the unladen case, strong departures are observed for the bubble-laden cases at all scales. The departures are strongest for the ratios involving 11 , which is to be expected since this is the direction of the mean trajectory of the bubbles. Moreover, for each of the ratios plotted in figure 5, the bubble-laden cases reveal a bump at = ( ), indicating the injection of anisotropy into the flow due to the bubbles and their anisotropic motion in the flow due to the buoyancy force acting on them.
While comparisons of the ratios of 11 , 22 , 33 is a standard way to analyse the multiscale anisotropy in the flow (van de Water & Herweĳer 1999;Carter & Coletti 2017), this method provides limited quantitative insight into the degree of anisotropy. For example, while the behaviour of these ratios is well-known for an isotropic flow in the limits / → 0, / (1), their behaviour for intermediate / is not known and cannot be determined simply by the condition of isotropy. It is therefore desirable to provide a simple measure of anisotropy based on that applies at all scales, and also helps to visualize the anisotropic behaviour. In the context of the Reynolds stress tensor, this can be achieved using either the Lumley triangle (Lumley & Newman 1977) or else the barycentric map (Banerjee et al. 2007). As describe in the next subsection, however, these methods cannot be directly applied to quantify multiscale anisotropy associated with , and modification is required.

Quantifying and visualizing scale-dependent anisotropy
Characterizing the scale-dependent anisotropy associated with is considerably more involved than that based on the Reynolds stress . This is because the relationship between componentiality and isotropy breaks down at sub-integral scales. For example, in the case of the Reynolds stress tensor, its components in the isotropic state are ≡ /3, according to which the components in all three Cartesian directions are the same (a "threecomponent flow"). In the case of , its isotropic state is (assuming the flow is incompressible) At the large scales ( / ) = 0, and we have = = /3, and hence the isotropic state corresponds to a three-component flow. On the other hand, in the limit → 0 we have = 2 /15 (Pope 2000) and = 2 − 2 . (4.4) In this case the isotropic state does not correspond to a three-component flow, but rather the components transverse to are twice as large as those parallel to . Hence, in general there is no correspondence between the componentiality of the flow and isotropy. For this reason, measures such as − ( /3), which have previously been used to define scale-dependent anisotropy (e.g. Brugger et al. 2018), do not in fact quantify isotropy, but rather only quantify deviations from the three-component state.
In view of these considerations, anisotropy must be quantified by the deviation of from , rather than from /3. To this end, we define two normalized tensors and with components and ( ) ≡ , (4.6) respectively. Then, following Banerjee et al. (2007), we may express these components in terms of basis matrices = I 1 1 + I 2 2 + I 3 3 , (4.7) = J 1 1 + J 2 2 + J 3 3 , (4.8) where the circumflex (·) denotes the tensor components expressed in the tensor principal axes, with = diag( 1 , 2 , 3 ), where 1 2 3 are the ordered eigenvalues of , and = diag( 1 , 2 , 3 ), where 1 2 3 are the ordered eigenvalues of . The basis matrices 1 , 2 , 3 represent the one-component, two-component, and threecomponent limiting states, and are given by 1 = diag(1, 0, 0), 2 = diag(1/2, 1/2, 0), and 3 = diag(1/3, 1/3, 1/3), respectively. The coefficients in (4.7) and (4.8) then measure how close and are to any of the three limiting states, and can be determined from the eigenvalues of the associated tensors To visualize the componentiality of the flow implied by these coefficients, a barycentric map (Banerjee et al. 2007) can be defined as = I 1 1 + I 2 2 + I 3 3 , (4.11) = I 1 1 + I 2 2 + I 3 3 , (4.12) and similarly for the J components. Here, ( 1 , 1 ) = (1, 0), ( 2 , 2 ) = (0, 0), and ( 3 , 3 ) = (1/2, √ 3/2) are the three corner points corresponding to the limiting states of componentality. They can be chosen arbitrarily, but are commonly set to the corner points of an equilateral triangle, as this aids the interpretation of the visualization. We note that in the traditional Lumley triangle, it is the fact the trace of the Reynolds stress anisotropy tensor, , is zero that allows its properties to be represented in a two-dimensional triangle using two independent eigenvalues of . This restriction does not apply to the Barycentric map approach, which is why we are able to represent the properties of and in a two-dimensional map, even though the trace of these tensors are not zero.
Provided that ( / ) 0 (true for a homogeneous flow), then the projections of in the plane orthogonal to are greater than or equal to . As a result, is then confined to the left side of the triangle, corresponding to a state of axisymmetric contraction. Moreover, as is decreased, moves monotonically away from the three-component state towards the two-component state. On the other hand, may exist anywhere within the triangle.
Following this approach, and can be mapped to a location ( ( ), ( )) in the barycentric map, and the linear distance between the coordinates corresponding to these two tensors then gives a measure of the anisotropy at that scale. In particular, the anisotropy at a given scale is determined by (4.13) where we have used (4.8) and (4.7). For an isotropic flow, = and = 0. At the large scales of an isotropic flow, I 1 = I 2 = 0 and I 3 = 1, so that the coordinates for are ( , ) = (1/2, √ 3/2). Therefore, at the large scales ( ) recovers the property that anisotropy is related to distance from the top of the triangle which corresponds to the threecomponent state, just as for the single-point Reynolds stress tensor (Banerjee et al. 2007). By contrast, for an arbitrary scale in the flow, anisotropy is related to the distance in the barycentric triangle from a point along the left side of the triangle, the location of this point representing . This is illustrated in figure 6. Figure 6: Representation of (4.6) and its isotropic form (4.5) in barycentric map (a) for all cases considered. The transparency of the color reflects the separation from small to large.
(4.13) is represented by the red dashed lines (b) for the case BiDisp, highlighting e.g.
(Δ) for the smallest scale and (1.1 ) for the largest scale.

Application of the new method
We now turn to apply the new method described in § 4.2 to our DNS results. The trajectories for in figure 6(a) all lie on the left side of the triangle, corresponding to the state of axisymmetric contraction, associated with the fact that 11 = 22 33 . The trajectories for at large are near to the three-component upper corner of the triangle, however, for the largest / for which we have data, ( / ) ≠ 0 and so the exact three component state = I/3 is not observed. As is decreased, the trajectory for moves down the left side of the triangle, towards the two-component corner, reflecting the fact that for / → 0, 11 = 22 = 2 33 . The results in figure 6(a) for show that the bubble cases start closer to the one-component corner of the triangle than the unladen case. This reflects the fact that the fluctuations in the gravity direction are much stronger due to the buoyancy forces acting on the bubbles. The trajectories of in the triangle are highly nonlinear, and show a tendency to migrate towards the axisymmetric contraction side of the triangle as is decreased, consistent with an approach to isotropy as is decreased, although never obtaining an isotropic state. The exception to this is the SmFew case for which the trajectory of actually approaches the one-component corner of the triangle as is decreased. This surprising behaviour is consistent with that observed in figure 5.
A quantitative measure of the scale-wise anisotropy is provided by , introduced in § 4.2. This quantity provides a linear measure of anisotropy, and is zero for an isotropic flow. For illustration, in figure 6(b) we show the trajectories of and for the BiDisp case and join the starting and ending points of these trajectories with dashed lines. The length of the line connecting the starting points represents ( = 1.1 ), while the length of the line connecting the ending points represents ( = Δ), where Δ is the grid spacing in the spanwise direction. The plot shows that (1.1 ) > (Δ), such that the flow is more isotropic at the smaller scales in BiDisp case.
In figure 7 we show the results for as a function of / for all of the cases. In general, monotonically decreases with decreasing , consistent with the return towards isotropy prediction. However, for these cases isotropy is never fully recovered, with anisotropy persisting into the dissipative range scales. Such behaviour for unladen turbulent flows has also been observed experimentally, including at much higher Reynolds numbers, as seen in studies by Kurien & Sreenivasan (2000) we note that for 0.02 (most clearly seen in the semi-log plot of figure 7b) the return to isotropy is interrupted for the unladen case, and actually becomes larger as is further decreased, implying increasing anisotropy at these scales. Similar behaviour has also observed in both experimental (Carter & Coletti 2017) and numerical studies (Meneveau 1991;Bos et al. 2007), where it was suggested to be due to anisotropic intermittency at the dissipative scales.
The results in figure 7 show that the departure from isotropy is in general much more pronounced for the bubble-laden cases than for the unladen case. This is again due to the fact that the bubbles introduce a significant source of momentum into the flow in the direction parallel to gravity due to the buoyancy force they experience. For the unladen case, the anisotropy is quite weak in the channel centre because the mean-shear that generates the anisotropy is weak in that region of the flow. We also notice from figure 7(a) that the level of anisotropy decreases with increasing bubble Reynolds number, . In particular, the anisotropy across the scales generally decreases in the sequence SmMany (SmFew), BiDisp to LaMany, which corresponds to increasing . This behaviour is likely due to both the -dependent structure of the wakes produced by the bubbles (see the supplementary materials to Santarelli & Fröhlich 2016) and also the path of the rising bubbles which becomes more chaotic and less uni-directional as increases (Horowitz & Williamson 2010;Ern et al. 2012). It is also interesting to note that cases SmFew and SmMany which have similar , but different , have a similar level of anisotropy. This suggests that plays more of a role than in determining the contribution of the bubbles to the flow anisotropy, at least over some regimes of (of course it cannot be true in general since, for example, when = 0 the bubbles have no effect on the flow). Another observation prompted by figure 7 is that the shape of the curve is remarkably similar for the three cases involving the highest (SmMany, LaMany, and BiDisp). This shape may be approximately divided into three regimes: for > ( ), a slow return-to-isotropy regime occurs in which the flow relaxes towards isotropy as is decreased; a bump (seen more clearly in the semi-log plot) is observed at the scale of the bubble = ( ); and a third regime at < ( ) where the flow again relaxes towards isotropy, but at a much faster rate with decreasing than it does in the first regime. This rapid approach towards isotropy is never fully successful, however, with significant anisotropy persisting at the smallest scales. For the SmFew case, the first two regimes can also be identified, however, we observe scale-independent anisotropy for < ( ). A possible reason for this difference is that due to the low gas void fraction in the SmFew case, scales at < ( ) are influenced more strongly than the other case by the single-phase behaviour, which as discussed, leads to actually increase at the smallest scales. Moreover, we notice that the rate of return-to-isotropy is stronger with increasing , even though the anisotropy becomes weaker with increasing . This is in contrast with the phenomenological notion from single-phase homogeneous anisotropic turbulence that the rate of return is typically faster for more strongly anisotropic flows, and is related to the nature of the slow part of the pressure-strain interaction term (Chung & Kim 1995). The present observation implies, however, that in bubbly flows the additional rapid pressure-strain term arising from the bubble-induced force production (Ma et al. 2020b) may play an important role in the overall return to isotropy in bubbly turbulent flows. More extensive DNS data would be required in order to investigate this in detail.

Energy transfer and third-order structure functions
The third-order structure function is of particular significance since it is related to the mean nonlinear energy transfer among the scales of a turbulent flow (Alexakis & Biferale 2018). In particular, in the Karman-Howarth type equation governing ( 3 , ), the inter-scale energy transfer term is (Hill 2001) ) (here we are writing 3 in the function argument as short-hand for 3 3 ). When F < 0 this corresponds to the nonlinear term supplying energy to the scale, while F > 0 corresponds to the nonlinear term removing energy from the scale. In a single-phase, three-dimensional turbulent flow, there is on average a downscale flux of energy with F < 0 in the inertial and dissipation ranges, representing the transfer of energy injected into the flow at the large-scales down to the smallest scales where it is dissipated due to viscous stresses. However, the presence of bubbles in the flow could change this picture, since bubbles produce kinetic energy at scales ( ) and some of this may then be transferred up towards the larger scales, as well as some being transferred down towards smaller scales. It is of great interest to understand this basic question of how the bubbles modify the energy transfer mechanisms and behaviour in the turbulent flow compared with the single-phase behaviour and mechanisms of strain self-amplification and vortex stretching (Carbone & Bragg 2020; Johnson 2020). However, since our DNS data is for a fixed streamwise and wall-normal location, we are only able to compute the F 3 contribution in (5.1), and thereby cannot fully determine F . Nevertheless, we can still partially explore the question of the effect of the bubbles on the energy transfer by considering F 3 .
In figure 8(a) we plot F 3 for the unladen case. The results show that at larger scales F 3 > 0, while at smaller scales F 3 < 0, corresponding to energy being taken from the large scales and passed down to the small scales via the nonlinear transfer term. At this low Reynolds number, there is no inertial range over which F 3 is constant, indicating that there is no cascade in the strict sense. In figure 8(b) we plot F 3 for all the cases. The first thing that is apparent from these results is that the nonlinear energy transfer is in general much stronger for the bubble cases than the unladen case. For example, for the LaMany case, the peak magnitude of F 3 is three orders of magnitude larger than for the unladen case. The second thing is that the qualitative behaviour of the bubble cases is the same as the unladen case, indicating that energy is passed downscale from the large scales to the small scales in the flow. The study of Pandey et al. (2020) explored a homogeneous bubbly turbulent flow using a Fourier space analysis and also observed that there is a downscale energy transfer, although their flow was considerably different to ours, and as noted in the introduction, there are some issues with their analysis.
Our results indicate then that the introduction of bubbles into the turbulent flow does not lead to an upscale energy transfer, at least for the cases we have considered and for the F 3 contribution. One possible reason for this is that the energy being sent downscale from the large scales (where it  is produced due to the mean-shear production), simply overwhelms any energy being transferred upscale from the bubbles, so that overall the energy transfer is downscale. However, since in most of these bubble cases the dominant source of TKE comes from bubble induced production (Santarelli & Fröhlich 2016), rather than mean-shear production, this does not seem likely. A second possible reason is that since the scales at which the bubbles are producing kinetic energy, ( ), are also scales at which the viscous forces in the flow are strong, then most of the kinetic energy produced by the bubbles may be simply dissipated before it is able to be transferred upscale by nonlinear forces in the flow. A third possible reason relates to the fact that although the bubble length scale is , they may inject energy at scales significantly larger than this. Indeed, the wakes produced by bubbles can have a length that is (10 ). For the Reynolds number of our DNS the scale separation between the bubble diameter and the large scales of the flow is not very large, and therefore the bubbles may actually directly inject TKE into scales on the order of the large scales of the flow, and this energy is then able to be sent downscale via the nonlinear energy transfer term. If this third explanation is true, then an interesting implication is that if BIT flows can develop with ≪ ℓ (where ℓ is the integral length scale of the flow), then the maximum scale at which the bubbles will inject energy into the flow will be much smaller than ℓ . In this case, an upscale energy transfer must occur at scales since otherwise the large scales would have no source of energy.
In figure 9 we plot one of the contributions to F 3 , namely 333 / . In contrast to figure 8, this result shows that for the SmFew and SmMany cases there is an upscale transfer of energy, with energy being extracted from the small scales and passed to the large scales. This reveals a subtle effect of the bubbles; while overall the bubble cases still exhibit a downscale energy transfer, the energy transfer associated with particular components of the velocity field do exhibit an upscale energy transfer. Figure 10 shows 333 / 3/2 33 (the skewness of the longitudinal velocity increment), as well as the unnormalized form 333 . For the present unladen case, in the channel centre we find 333 / 3/2 33 ≈ −0.4 for → 0, which is in close agreement with the value found in isotropic turbulence (Ishihara et al. 2007;Davidson 2015). The Reynolds number of the channel flow is too low for there to be a clear inertial range, and therefore the unladen data for 333 does not show evidence of the "Fourth-Fifths law" 333 = −(4/5) predicted for the inertial range of isotropic turbulence by Kolmogorov (1941a). For the cases involving bubbles, a first observation is that 333 / 3/2 33 is positive for all cases at = ( ), and for the SmFew and SmMany cases it remains positive in the limit → 0, while for others it becomes negative in this limit. The data for 333 indicates that the fluctuations are much larger for the bubble cases than for the unladen case, especially those with higher volume fractions, while the results for the skewness shows that the bubble cases have a skewness that is either larger or smaller in magnitude than the unladen case, depending upon the case. The SmFew case is particularly interesting, since it exhibits a very large value of the skewness for < ( ), even though this case has the smallest volume fraction . However, as expected, the results for 333 for this case reveal that the magnitude of 333 is much smaller for the SmFew case than for the other three bubble cases that have higher values of . This shows how sensitive certain properties of the turbulent flow are to the presence of the bubbles, even for relatively low . It also implies that the bubbles could have a significant effect on tracer particle dispersion in turbulent flows, whose irreversibility is intimately related to the asymmetry in the velocity increment distributions (Jucha et al. 2014;Bragg et al. 2016;Bragg 2017).
The transverse structure functions 111 and 222 were also considered for each case. The associated skewness values were found to be very close to zero for all cases (not shown here).

Fourth and sixth-order structure functions
We now consider the fourth and sixth-order structure functions to explore how the bubbles influence the intermittency of the background turbulence. In figures 11(a,c,e) we consider the normalized fourth-order moments (i.e. the flatness of the velocity increments) in each of the directions, denoted by 1111 / 2 11 , 2222 / 2 22 , and 3333 / 2 33 . If the velocity increments had Gaussian probability distributions, then these quantities would be equal to three at all scales. Departures from three indicate non-Gaussianity in the velocity increment distributions, while spatial dependence of these quantities indicates that the velocity increments are not only non-Gaussian, but also intermittent (Frisch 1995).
For the unladen case, the results in figures 11 (a,c,e) show that at the large scales, the quantities are close to three, but increase as decreases. Furthermore, the results show that the departures from Gaussianity are stronger for the transverse components than they are for the longitudinal component, also observed in isotropic turbulence (Ishihara et al. 2009). The effect of the bubbles on these quantities is non-trivial. At the very smallest scales, the bubble cases show much stronger non-Gaussianity and intermittency in the turbulent flow than for the unladen case. A similar finding was also observed in the experiment of Rensen et al. (2005) for their fully developed turbulent bubbly flow and in Biferale et al. (2012) when comparing the small-scale properties of boiling and non-boiling convective turbulent flows. However, with the exception of the SmFew case, the bubbles tend to suppress intermittency in the flow at scales comparable to their diameter, for the streamwise and wall-normal velocity directions. Roughly speaking, it is only at scales smaller than ( ) that the flatness starts to increase significantly. The flatness results for the SmFew case reveal much stronger non-Gaussianity and intermittency in the flow compared to the other three bubble cases that have much higher gas void fraction (≈ 7 times than SmFew). Some insight into this behaviour can be found by considering the ways in which the bubbles modify the turbulent flow. The most significant qualitative modification the bubbles make to the turbulent flow is the wakes they generate, and the flow properties in these wakes will typically be significantly different from the properties of the background turbulence. When the number of bubbles is not too large (but still large enough to make a statistically significant contribution to the flow properties), there are relatively few regions in the flow where the flow differs due to the wakes from everywhere else (i.e. these regions are "rare"), breaking self-similarity in the flow and enhancing intermittency compared to the case without the bubbles. As the number of bubbles increases, the regions of the flow occupied by the wakes becomes less rare, and hence their contribution to the intermittency reduces. Moreover, as the number of bubbles/size of the bubbles increases, the possibility for wake interaction increases, and this can enhance mixing in the flow, which would reduce the impact of the wakes on the flow intermittency.
The behaviour observed for the flatness is also observed for the sixth-(superflatness, see figures 11b,d,f ) and eighth-order (hyperflatness, not shown here) structure functions. A most surprising result concerns the relationship between and the level of intermittency. For example, as shown in figure 2, the LaMany case has a much larger than the SmFew case, yet the latter exhibits much stronger intermittency in the flow than the former. This is in striking contrast to single-phase turbulence where small-scale intermittency increases with increasing Reynolds number (Frisch 1995). This indicates that the mechanisms generating intermittency in bubbleladen turbulent flows might be significantly different from those in single-phase turbulence. Further data would be required to explore the mechanisms responsible for this behaviour.

Conclusions
In this paper we have analyzed various properties of bubble-laden turbulent flows across the range of scales of the flow, including the anisotropy, energy transfer, and intermittency in the flow.
To explore the anisotropy of the flow at different scales, we developed an extension of the barycentric map approach that was previously developed for analyzing the one-point Reynolds stress tensor (Banerjee et al. 2007). In our approach, the anisotropy at any scale may be quantified and visualized, as well as providing information on the componentiality of the flow at the scale. Using this we were able to explore how the bubbles modify anisotropy in the flow. We found that the bubbles significantly enhance anisotropy in the flow at all scales compared with the unladen case, and that for some bubble cases, very strong anisotropy persists down to the smallest scales of the flow. The strongest anisotropy was observed for the cases involving small bubbles. Deviations from the behaviour that would be expected for an isotropic flow occur not only in the streamwise/gravity direction, but also in the other directions.
Concerning the energy transfer among the scales of the flow, our DNS data did not allow us to thoroughly explore this, but we were able to consider a number of its important aspects. Our results indicate that for the bubble-laden cases, the energy transfer is from large to small scales, just as for the unladen case. However, there is evidence of an upscale transfer when considering the transfer of energy of particular components of the velocity field, rather than the full kinetic energy involving all three components of the flow. We also conjectured, however, that in a bubble-laden channel flow at sufficiently high Reynolds number, and with channel height sufficiently large compared with the bubble size, there may exists an upscale transfer of energy at scales much larger than the bubble diameter. Our DNS does not lie in the parameter range required to observe such behaviour. Another important finding is that although the direction of the energy transfer is the same with and without the bubbles, the energy transfer is much stronger (by several orders of magnitude) for the bubble-laden cases, suggesting that the bubbles play a strong role in enhancing the activity of the nonlinear term in turbulent flows.
We also considered the normalized forms of the fourth and sixth-order structure functions, corresponding to the flatness and superflatness of the fluid velocity increments. These results revealed that the introduction of bubbles into the flow strongly enhances the intermittency of the turbulence in the dissipation range, but suppresses it at scales comparable to the bubble diameter. The SmFew case, however, shows enhancements of intermittency at all scales. We interpreted the effect of the bubbles on the flow intermittency in terms of the contributions of the bubble wakes to the overall properties of the turbulence. This strong enhancement of the dissipation scale intermittency has significant implications for understanding how the bubbles might modify the mixing properties of turbulent flows, and the associated large velocity gradients. These will be the subject of future investigations.