Abstract
Although bubbly flows have been extensively applied to many industrial fields, efficiency improvements of bubbly equipments are greatly limited due to the absence of accurate understandings of liquid turbulence modulation by bubbles. In order to obtain an accurate mechanism explanation, an Euler-Lagrange two-way method was used to investigate and analyze deeply the liquid turbulence modulation by microbubbles in this paper. The present investigation shows that for the bubbly flow laden with microbubbles the liquid turbulence modulation is mainly caused by the drag force, and the turbulence modulation level is related to the microbubble diameter. A detailed mechanism explanation on the physical phenomenon was presented here.
1. Introduction
Due to good physical properties, bubbly flows are frequently encountered in many industrial fields such as smelting of metals, two-phase heat exchanger, and drag-reduction engineering by microbubbles [1–3]. In order to understand hydrodynamic characteristics of bubbly flows, many numerical and experimental investigations have been carried out in the past 30 years [1–21]. As pointed out by some researchers [4–6], of many hydrodynamic properties of the bubbly flow, it is very important for the system design of the drag reduction by microbubbles to deeply understand the liquid-phase turbulence modulation by microbubbles. In order to understand the turbulence modulation mechanism by bubbles, many experimental and numerical investigations had also been performed for different kinds of turbulence. A detailed review on the liquid turbulence modulation by bubbles can be referred to [7]. A review on physical phenomena of bubbly flows had been presented in [8]. Ceccio [9] presented a detailed review on the drag reduction mechanism by microbubbles. Thus, the related content cannot be repeated again here.
As a matter of fact, the mechanism on the liquid turbulence modulation by bubbles cannot be fully understood yet and the existing mechanisms cannot be generally accepted so far. From the above review, one can see that the present literature only showed a modulation trend on the liquid-phase turbulence by bubbles and did not analyze why the liquid turbulence is enhanced or damped by bubbles. That means that the modulation mechanism of the liquid-phase turbulence by bubbles is not essentially clear. In order to bring an effective control of bubbly flows into effect in industrial processes, it is very important to deeply understand interactions between the bubble and liquid phases. In view of the previous reason, a mechanism of the liquid-phase turbulence modulation by microbubbles was presented in this paper. In order to understand the liquid turbulence modulation by microbubbles, we had investigated influences of interphase force closures on the phase distribution and the liquid-phase turbulence so as to find out the most important force influencing both of them in the previous work [10]. Different from the previous job, we want to investigate how the liquid turbulence is modulated by microbubbles with different diameters under the action of the key force in this paper.
2. Methodology
The numerical method can independently achieve the full control of various effects on the physical phenomena and can obtain very detailed data compared with the experimental one. Especially, the Euler-Lagrange model can provide the most detailed insight into single-bubble dynamics as well as into bubble-bubble and bubble-fluid interactions. It can be used to evaluate influences of basic physical and geometric parameters like inertia, viscosity, surface tension, bubble size, and gas volume fraction on the evolution of bubble swarms, interactions between bubbles, and the induced flow structure of liquid [11]. And thus, in this paper, we investigated liquid turbulence modulation by microbubbles with a developed Euler-Langrage two-way method in our previous job [8]. The flow field was simulated with direct numerical simulations (DNS) in the Euler frame of reference, while the bubble dynamics were fully analyzed by Newtonian motion equations. The computational domain (see Figure 1) is a vertical upward channel with respect to the gravity direction, in which x, y, and z correspond to the streamwise, wall-normal, and spanwise direction, respectively. In order to understand the turbulence modulation mechanism induced by microbubbles, to learn about the main reason why the liquid-phase turbulence is modulated by microbubbles, and to understand the relationship between the microbubble size and the liquid turbulence modulation, different computational cases were designed, as listed in Table 1. In Table 1, FB, FD, FM, and FL denote the buoyant force, the drag force, the added mass, and the lift force; h is half the channel height, d is the bubble diameter; We is the Weber number, We = ρ f uτd/2σ, where σ is the surface tension coefficient, uτ is the frictional velocity, ρ f is the liquid density; ν f is the kinematic viscosity of liquid. In the computational process, interactions between microbubbles were neglected with respect to the fact that the global void fraction (α0 = 1.34 × 10−4) is less than 10−3 according to [12]. Additionally, microbubbles can be thought as spherical in shape; that is, the bubble deformation was also neglected due to too small in size (i.e., We ≪ 1) [13, 14]. In order to simulate the bubbly flow, the single-phase turbulent channel flow is firstly prepared, and the friction Reynolds number is Reτ = 150, which is based on the wall friction velocity (uτ) and half the channel width (h). When the steady turbulence is statistically reached, a set of bubbles, of which initial velocities are set to be zero, are injected into the flow field. Note that the shear velocity is still defined with the wall shear stress as usually done for the single-phase problem so as to investigate how the liquid turbulence will be modulated by interactions between microbubbles and liquid turbulence. Namely, for the present definition of the shear velocity, the turbulence modulation is only dependent on interactions between microbubbles and the liquid turbulence if that happens. Readers can refer to the literature [15] to learn about the other definition on the shear velocity. Additionally, according to the investigation of Molin et al. [15], the wall lift force has a small influence on the total lift force, and it cannot cause the reversal direction of the total lift force for microbubbles. Therefore, the wall lift force was neglected for the present investigation too.
Computational cases.

Flow geometry and coordinate system.
For the bubbly flow, the liquid phase is assumed to satisfy a modified Navier-Stokes equation in which the feedback force on liquid is regarded as an effective body force. The nondimensional governing equations for the gas and liquid phase can be written as follows:
where Reτ is the frictional Reynolds number (
The lift coefficient, CLF, is calculated by
where
Here Sr
b
and Re
b
are nondimensional shear rate and bubble Reynolds number. They can be defined as
Presently, a finite difference scheme was used to discretize the governing equations. A second-order central difference scheme is used for the spatial discretization. For the time integration, the Adams-Bashforth scheme was used for all the terms except that the implicit method was used for the pressure term. The velocity-pressure coupling was made with the MAC method. A staggered grid system was used to prevent a checkboard pressure field, which causes the fact that velocity components are located at cell faces and other variables are located at cell centers. To calculate interphase forces related to the relative velocity between bubbles and local liquid, three-dimensional 8-node combined with a two-dimensional 4-node interpolation polynomials are used to compute liquid velocities at the same positions of bubbles (near the wall, the interpolation scheme switches to one-sided). The pressure Poisson equation was solved by using the multigrid method. Uniform grids were used in the streamwise and spanwise directions. Nonuniform grids were used in the wall-normal direction. The grid system for the present investigation was absolutely similar to our previous jobs [8, 10]. Bubble velocities and bubble displacements were calculated by the acceleration integral equation with the second order Crank-Nicholson method. A detailed description can be referred to the literature [8].
3. Validation of Grid Resolution
For the present simulation, the computational domain in wall units is 1500 × 300 × 750. A grid system of 64 × 64 × 64 is adopted. The grid spacing of Δy+ varies from about 0.45 near the wall to 9 at the center of channel in wall-normal direction, and the grid spacings are Δx+ = 23.4 and Δz+ = 11.7 in the x and z directions, respectively. The Kolmogorov length scale is calculated by the following:
For the present simulation, for the single-phase liquid, ε+ varies from 0.03 at the center of channel to 0.34 near the wall as seen from Figure 5. The Kolmogorov length scale in wall units has the smallest value of about 1.3 near the wall and the largest value of 2.4 at the center of the channel. The above analyses show that the grid size in the wall-normal direction is smaller than the Kolmogorov length scale at the near-wall region and comparable to the Komogorov length scale in the centre region. In the streamwise and spanwise direction, grid sizes are larger than the Kolmogorov length scale. The Kolmogorov time scale in wall units is
For the current investigation, it has the smallest value of 1.43 ν/uτ2 near the wall and the largest value of 3.20 ν/uτ2 at the channel center, respectively. However, the present time step is 0.015 ν/uτ2, which is much smaller than that of the Kolmogorov time scale. The present computational criteria of DNS are in good agreement with [10]. Additionally, in order to testify accuracy of the present DNS for the single-phase flow, several turbulent statistics from our results were compared with those of Kim et al. [21] in [10]. Although the frictional Reynolds number, Reτ = 180 in [21], is larger than that in the present Reτ = 150, profile trends for all statistics are very similar. In view of a prolix depiction, comparisons are not presented again.
4. Results and Discussion
Since the microbubble modulation on the streamwise component of turbulence intensity is uneven along the wall-normal direction and the modulation on the spanwise one is similar to that of the wall-normal one, therefore the wall-normal component profile as the function of the wall coordinate (y+ = yuτ/ν f ) is only presented in Figure 2. Figure 2(a) exhibits the influence of interphase force closure on the wall-normal component of liquid turbulence intensity, and Figure 2(b) displays its dependence on the microbubble size. As shown in Figure 2, the liquid-phase turbulence intensity component is suppressed by microbubbles in the wall-normal direction. That is in good agreement with results of spherical bubbles in the literature [4]. One can see from Figure 2(a) that the suppression of microbubbles on the liquid-phase turbulence is mainly caused by the drag force, and the lift force seems to reduce this suppression. However, the added mass force has a very weak influence on the liquid-phase turbulence intensity, which is in good agreement with conclusions of [17, 18]. When the buoyant force (FB) is neglected in the computation, the mean slip velocity between microbubbles and liquid is absent, which leads to the absence of the drag and lift force. Under the circumstance (i.e., for Case E), microbubbles have almost no influence on the liquid-phase turbulence intensity (namely, the curves overlap for Cases A and E in Figure 2(a)) as pointed out by Chahed et al. [19]. The previous phenomenon further shows that the added mass force, if only being considered, cannot cause any changes in the liquid turbulence intensity.

Liquid wall-normal turbulence intensity profiles.
Figure 2(b) shows that the suppression of microbubbles on the liquid-phase turbulence intensity component in the wall-normal direction is tightly related to the microbubble size for the same global void fraction. The suppression becomes very evident with increasing the microbubble diameter initially. However, the suppression decreases again when the bubble diameter increases to a certain threshold value. For the present investigation, the threshold value seems to be d b + ∼ 2.25. The reason on the wall-normal turbulence intensity suppression by microbubbles is showed subsequently.
Figure 3 shows the liquid-phase instantaneous velocity vector field at the z-y plane, and Figure 4 displays the liquid-phase streamwise vorticity isosurfaces corresponding to −0.1 and 0.1 values. One can see from Figures 3 and 4 that modulations of microbubbles on the bursting event and the streamwise vorticity of liquid are similar, and both of them are tightly related to the microbubble diameter. It can be seen from Figures 3 and 4 that both the turbulence bursting event and the streamwise vortex are suppressed due to the microbubble injection. With increasing the microbubble diameter, the suppression degrees increase at first and then decrease, which corresponds to the modulation of microbubbles on liquid turbulence intensity.

Liquid instantaneous velocity vector field distribution at z-y plane.

Isosurfaces of liquid streamwise vorticity corresponding to −0.1 and 0.1 values.

Liquid production and dissipation.
As is well known to all of us, different from the laminar flow, the turbulence flow contains different-scale eddies. As a matter of fact, the liquid velocity fluctuations correspond to different scale eddies in the turbulence flow. If the different scale eddies are suppressed, the turbulence intensity will be reduced. It is well known that the large-scale eddies are responsible for the turbulence production extract energy from the mean flow and the Kolmogorov microscale eddies denote the turbulence dissipation translate part of energy into heat due to friction. If different scale eddies are suppressed, it will lead to reductions of turbulence production and dissipation as shown in Figure 5.
In order to analyze the influence of microbubbles on the turbulence eddy, a single bubble modulation on an eddy is analyzed in detail based on the mechanical theory. Figure 6 shows the mechanism of a microbubble influencing a single liquid eddy. When the diameter of microbubbles is very small, they can follow liquid well. They can penetrate into small eddies and move together with them. In the moving process, the surrounding liquid will exert a so-called drag force on the microbubble. According to the law of action and reaction, when liquid exerts a force on the microbubble, it will also bear an equal and opposite force from the microbubble. For Case B, when microbubbles move with liquid eddies, there are two kinds of drag forces exerting on them. One is the normal drag force (F N ) that is caused by the slip velocity under the action of the buoyant force as shown in Figure 6. This component of the drag force is the most important one. It is present only if the buoyant force is not neglected. Its reaction force will drag the local liquid away from the eddy and make it move together with the microbubble. That means that the eddy will be broken down in advance under the action of the reaction force before it continues to turn into the smaller scale eddy and even is transferred to heat in the end. The other is the tangential drag force as shown in Figure 6. Unequal liquid velocities on both sides of the microbubble (u l A >u l B ) will cause the fact that the tangential force on both sides of the microbubble (FTA = −FTB) shows different behaviors in the same direction. The reaction force of FTA will reduce the rotation speed of liquid located at point A. However, the reaction force of FTB will enhance the rotation speed of liquid located at point B. That shows that the tangential reactions will make the rotation speed of liquid become even from the vortex core to the outer in a liquid eddy. If a series of microbubbles are randomly scattered in the same eddy, the liquid rotation speed at each longitudinal position should be the same under the influence of the tangential reaction. That will also reduce the probability that the eddy is evolving into the smaller scale eddies or even into heat. And thus, the reaction forces decay the eddy in advance and prevent it from evolving and developing further. The decay and disappearance of eddies mean the turbulence intensity reduction. If the added mass force is taken into account in the computation, according to its definition, microbubbles will drag more liquid surrounding them to move, which will strengthen the eddy decay. That corresponds to the fact that the suppression of the wall-normal turbulence intensity component by microbubbles is slightly enhanced as shown in Figure 2(a). Additionally, microbubbles will spin close to the wall under the influence of the lift force. In this process, microbubbles will make liquid surrounding microbubbles rotate together with them, which will offset the influence of the drag force on eddies. That means that the lift force reduces the decay of the drag force on eddies. Besides, most microbubbles will accumulate near the wall due to the lift force as shown in Figure 7. In Figure 7, the local void fraction is defined as the ratio of the total volume of local microbubbles to one of local meshes in the wall-normal direction. One can see from Figures 2(a) and 7(a) that the more the microbubbles near the channel wall are, the weaker the liquid turbulence modulation by microbubbles is. One reason may be that the eddy scale near the wall is too small to be modulated by the present scale microbubble, and the other may be that the turbulence modulation does not happen in the viscous sublayer but in the other layer. Ortiz-Villafuerte and Hassan [20] pointed out that the microbubble layer near the wall has no major role in the turbulence modulation but the most important aspect is the accumulation of microbubbles in the buffer layer. It can also be seen from Figures 2(a) and 7(a) that Cases B and C for the strong turbulence modulation have relatively high void fractions in the buffer layer. In short, the lift force modulation on the turbulence modulation might be related to two facts: one is that it offsets the modulation of the drag force, and the other is that it can drive most of bubbles close to the wall, which reduce interactions between bubbles and the turbulence in the buffer layer.

Influence of microbubbles on an eddy.

Local void fraction profile.
Additionally, another turbulence modulation might happen when microbubbles are caught in two vortexes bordering on each other with the contrary rotation direction, as shown in Figure 8. The magnitude and direction of forces exerting on both sides of microbubbles are also different due to different instantaneous velocity fields on both sides of microbubbles. And thus, the two vortexes with different rotation directions will influence and entangle each other, and even decay under the connection of microbubbles until microbubbles escape from the bondage of one of them. It must be pointed out that from the transient point of view, the liquid turbulence modulation by microbubbles can happen in the arbitrary direction, but from the average point of view, the modulation may only occur in the streamwise direction because the average slip velocity has the definite value in the streamwise direction due to buoyancy.

Influence of microbubbles on two adjacent vortexes with the contrary direction.
As is shown in Figure 2(b), the suppression of microbubbles on the liquid wall-normal turbulence intensity component is directly related to the microbubble diameter. With increasing the microbubble diameter, the inhibition degree increases firstly and then decreases. According to the literature [18], the slip velocity increases with increasing the microbubbles diameter. The large slip velocity corresponds to the great normal drag force, which leads to comparatively large decay of eddies; thus, microbubbles with the big diameter have the great suppression on the turbulence intensity. When the microbubble diameter increases to a critical value, the following behavior of microbubbles to liquid becomes bad. Microbubbles cannot well penetrate into eddies and cannot effectively decay eddies, so the suppression of them on the turbulence intensity also become weak. As shown in Figures 2(b) and 7(b), microbubbles for Case D-1 are fewer than ones for Case D in the buffer layer but the turbulence suppression degree for Case D-1 is more evident than that for Case D due to the fact that the drag force exerting on microbubbles for Case D-1 is larger than one for Case D. And the microbubbles diameter for Case D-2 is the largest but the turbulence suppression degree is the weakest, which may be caused by the following facts: one is that the bubble diameter might be too large to penetrate into the smaller eddies and the other is that microbubbles are so few in the buffer layer that interactions between microbubbles and the turbulence are reduced.
If the microbubbles diameter continues to increase, the wake flow will become more and more serious and even the shedding vortex appears behind microbubbles. It will lead to the fact that the turbulence is not suppressed but is enhanced due to the addition of microbubbles. Based on the previous analysis, Figure 9 displays the relationship between the microbubble diameter and the suppression degree for the dilute bubbly flow. The vertical coordinate variable denotes the integral value of the suppression degree of the wall-normal turbulence intensity component along the height of the channel, and it is defined as
where vrms+,sp denotes the wall-normal turbulence intensity component of the single-phase liquid and vrms+,sp denotes one of the liquid phases. For the present investigation, the critical one (db,cr+) seems to be larger than 3. It was investigated that there is also a critical diameter of microbubbles for the drag reduction occurrence by microbubbles. Kanai and Miyata [22] and Shen et al. [23] pointed out that the drag reduction can be attained only when the bubble diameter is less than about 1 mm; and the drag reduction rate is generally higher when the bubble diameter is smaller.

The degree of liquid turbulence modulation as a function of the bubble diameter (only limited to the dilute bubbly flow).
5. Conclusions
Through the numerical investigation, it is discovered that microbubbles suppress the liquid-phase turbulence, the degree of which is tightly associated with the microbubble diameter. For the present investigation, the liquid turbulence suppression is mainly caused by the drag force exerted on microbubbles. The reaction of the drag force prevents the turbulence eddies from developing and evolving further into Kolmogorov scale ones, which directly leads to reduction of the turbulence fluctuation.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Footnotes
Acknowledgment
The authors gratefully acknowledge the financial support from the NSFC Fund (nos. 51376026 and 51225601).
