Abstract

As an inevitable trend for the sustainable development of the global economy, saving energy and reducing emissions are key goals for the entire world. The use of air bubbles to reduce viscous friction is one of the most effective approaches to achieve this goal, as it may significantly reduce the frictional drag of ships. However, the injection of air bubbles will change flow characteristics near the wall due to the significant differences in density and viscosity between air and water. In addition, parameters such as bubble size, bubble surface tension, bubble number and bubble position also affect the flow near the wall, resulting in significant diversity and instability in two-phase flow. To clarify the mechanism of these effects, a two-dimensional channel flow with air bubbles is studied using Computational Fluid Dynamics (CFD). The interactions between bubbles and water and between bubbles and wall are studied, and the detailed characteristics of bubbles moving in fully developed flow are considered. This study shows that the velocity gradient is the main factor influencing wall shear stress, and the presence of bubbles has a marked impact on the local velocity gradient distribution of the nearby flow. It is also found that shorter distance between a bubble and the wall enhances the flow interaction and leads to more significant perturbations of wall shear stress.

1. Introduction

The use of air bubbles is one of the most effective approaches to reduce drag for ships. The interactions during air bubbles, water flow and ship walls play an important role in the application of such bubble drag-reduction technology. In this paper, the interaction mechanism of bubbles, water and wall surface in two-dimensional pipeline flow is studied, and its mechanism is analyzed by numerical simulation. Bubbles can reduce frictional drag of a moving body by changing the properties of the near-wall fluid flow [13]. The maritime shipping industry is thus interested in employing this drag-reduction technology to reduce energy use and greenhouse gas emissions [4, 5]. The use of bubbles to reduce skin frictional drag of ships’ external turbulent flow is a long-term engineering goal [6, 7]. Although the industrial community mainly cares about drag-reduction effects for full-scale ships, understanding the mechanisms of complex water-bubble interactions and bubble-wall interactions is one of the key components for this goal [8]. It is generally acknowledged that drag reduction using bubbles depends on the properties of bubbly flows, such as bubbles’ size, void fraction and their interactions with water [9].

Experimental study is one of the most important methods for the investigation of bubbly flow, due to its reliability and robustness. A series of studies have been carried out on drag reduction using bubbles, and the influence of injected bubble was analyzed [1013]. Verschoof et al. [14] found that bubble deformability was a fundamental factor of bubble drag reduction in turbulent flow. Rawat et al. [15] discussed the influence of micro-bubbles on turbulent flow, and they discovered that the bubble effect plays a major role in the improvement of flow structures. Shen et al. [16] experimentally revealed that drag reduction with bubbles was obviously related to the volumetric flow rate of injected gas and the static pressure of boundary layer.

However, the limited physical space and overlap of bubbles in mixed-regime flows make it difficult to observe the complicated but important bubble-water-wall interactions in experiments. Complementary to experimental study, numerical simulation is an important option, since it can provide valuable flow details that are challenging to measure experimentally. In addition, numerical simulation has been applied to various fields, including heat transfer and nanofluids [1722]. Hence, a great deal of research has been carried out in the recent years through numerical simulations [2326]. Xu et al. [27] discovered that a sustained level of drag reduction can be achieved by seeding small bubbles in a turbulent channel flow through direct numerical simulation. Pang et al. [28] numerically simulated the drag reduction through the use of small bubbles, and found rules for drag reduction according to the bubbles’ size and the liquid-phase Reynolds number. Lu et al. [29, 30] used direct numerical simulations to model bubble deformation and drag reduction with bubbles, and found that bubble deformability usually resulted in different void fraction distributions. Lyu et al. [31] assessed the roles of several bubble parameters in the drag reduction for a SUBOFF submarine model, and found that the air injection velocity ratio and air volume fraction are critical for drag performance of the submarine. Marinho et al. [32] numerically studied the flow in curved ducts with Taylor bubbles using a commercial software package CFX, and discussed the transient effects of the air concentration on the volume fraction and viscosity of bubbly flows. Fox [33] reviewed the development of large-eddy-simulation tools for dispersed multiphase flows including use of bubbly flow for drag reduction. It is found that CFD modeling of bubbly flow still faces some significant challenges.

The prior studies in the area of drag reduction with air bubbles usually focused on bubbly flow with a mass of bubbles [3436]. Bubble drag reduction is cumulative with the action of discrete bubbles, so that study focusing on the change in drag of a plate due to discrete air bubbles, that is, a single bubble or two bubbles, is essential. However, the experimental measurement and data analysis for application of discrete small bubbles are seldom implemented, due to the insufficient precision and great difficulty of controlling such small discrete bubbles. Furthermore, few investigations have been conducted to characterize discrete bubbles in the boundary layer using CFD methods as well. Therefore, a simple two-dimension CFD model coupled with the volume of fluid (VOF) method is applied to study the bubbles’ effects on the flow characteristics and the wall shear stress in this paper. Several parameters including bubble size, bubble number, surface tension, and initial bubble position are taken into account to demonstrate the interesting bubbly flow.

2. Numerical Configuration

In this study, the ANSYS Fluent commercial software package is applied for two-dimension numerical simulation to study bubbly flow characteristics in a two-dimensional channel. A finite volume method is used to discretize computational domain and the VOF method is applied for bubble interface capture. The standard k-ɛ turbulence model is used to demonstrate turbulent flow near the channel’s boundaries at the walls.

The equation for conservation of mass is given as

The momentum equation is

In these equations, u is the velocity, the subscripts i and j are 1 or 2 and represent x or y components, respectively, and μ is the viscosity of water.

The standard k-ɛ model is based on equations describing the turbulence kinetic energy k and turbulence dissipation rate ɛ. The transport equation for the turbulent kinetic energy k is [37]

The transport equation for the turbulence dissipation rate ε iswhere is the generation of turbulence kinetic energy due to buoyancy, is the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, and , , and are constants. Parameters and are the turbulent Prandtl numbers for k and ɛ, respectively. Terms and are user-defined source terms, and the term is the generation of turbulence kinetic energy due to the mean velocity gradients and can be defined as

The turbulent viscosity, , is determined by combining k and ε as follows:where is a constant.

The VOF method implemented in Fluent uses a multidimensional universal limiter with an explicit solution algorithm. Together with an interface compression algorithm, this method ensures a sharp interface and limits the volume fraction of the phase to values between 0 and 1. The momentum equation of the phase takes the following formwhere is the mass transfer from phase to phase and is the mass transfer from phase to phase . The continuum surface force model has been implemented so that the addition of surface tension to the VOF calculation results in a source term in the momentum equation. In this model, the surface tension is constant along the surface and where only the forces normal to the interface are considered. The pressure drop across the surface depends upon the surface tension coefficient, and the surface curvature is measured by two radii in orthogonal directions.

The ANSYS ICEM commercial software package is used for geometry modeling and mesh generation. The region of flow simulation, which is set to have a 10 mm height and 100 mm length, is discretized with uniform quadrilateral elements as shown in Figure 1. Boundary conditions including velocity inlet, wall, and outflow are used at the boundaries of the flow region, as shown in Figure 1. Mesh dependence and time-step dependence are studied, and a mesh number of 0.4 million (Figure 2) and time step of 10−5 s are verified to be adequate for this two-phase flow simulation. Transient numerical simulation is first conducted for the channel with 100 mm length until a fully developed two-dimension turbulent flow is achieved. The water flow simulation is conducted with parameters as follows: the water density is ρ = 998.2 kg/m3, the dynamic viscosity is μ = 0.001003 N·s/m2, the flow velocity is u = 0.5 m/s, and the Reynolds number Re = ρuh/μ = 4976 is based on the channel height . Typical velocity profiles at the end of every simulation period are presented in Figure 3(a), and the velocity profiles between the 3rd and 4th differs a little (Figure 3(b)). A more clear illustration of the velocity profiles for the 4th period is also shown in Figure 4 to further illustrate that the flow is fully developed at that time.

Experimental data for two parallel plates by Whan and Rothfus [38] are introduced to validate the present CFD framework, with comparison of velocity profiles at Reynolds numbers around 5000, as shown in Figure 5, where is the maximum of velocity along the x direction. With consideration of the differences of Reynolds numbers and the imperfection of present two-dimensional model, it is found that the results from CFD and experiment agree with each other in an acceptable degree. As a consequence, the CFD framework is applied in the subsequent simulations of bubbly flows.

3. Results and Discussion

To investigate the influence of air bubbles on two-phase flow characteristics, parameters including bubble size, surface tension, bubble number, and initial bubble position are taken into account during the simulation. Effects of those parameters on the wall shear stress at a certain monitoring point and frictional drag at the wall are monitored, as shown in Figure 6. The simulation cases and corresponding drag at the wall are listed in Table 1, and the condition descriptions for these cases are listed in Table 2.

It is well-known that friction at the wall is decided by wall shear stress, and the wall shear stress is described by the Law of Newton inner friction , in which is dynamic viscosity of the fluid, is the local velocity along the tangential axis of the wall and y is the normal axis of the wall. To analyze the influence of the dynamic viscosity and velocity gradient on the wall shear stress, three typical time points in Case 9 are selected, as shown in Figure 7. The wall shear stress decreases rapidly after reaching a peak point (t > 0.16 s) and then increases quickly to another peak point. That is because the distance from the bubble to the wall dh is 0.5 mm in Case 9 and the bubble stays close to the wall, as the bubble approaches the monitoring point, the effects of bubbles on the flow distribution gradually increase, and the velocity gradients at time points a1 and a3 are bigger than those at other points. This makes the wall shear stress at time points a1 and a3 bigger than that without bubbles. As shown in Figure 8, the velocity gradients at time points a1, a2, and a3 have little difference, while the wall shear stress at time point a2 is much less than that of the other two time points. The reason for this phenomenon may be the fact that the bubble is right under the monitoring point; as a result the fluid viscosity is mainly decided by the air viscosity. In other words, the major factor affecting wall shear stress is the fluid viscosity when the bubble is close to the wall; the effect of the velocity gradient on the wall shear stress is very weak.

3.1. Effects of Bubble Size

Bubbles with different diameters, that is, 0.5 mm and 1 mm, are modeled to study the effects of bubble size on the wall shear stress at the monitoring point and drag at the wall, corresponding to Cases 1 and 2 as listed in Table 1. The initial distance from the bubble to the wall is set to 5 mm. It is found that larger bubble makes a little more perturbation on drag than the smaller bubble does, and both bubbles lead to a drag increase instead of reduction, as shown in Figures 9 and 10. At the earlier stage, the drag at the wall clearly starts to increase at time point b1 for Case 1, and at time point b2 for Case 2, respectively. Furthermore, an interesting phenomenon occurs that the drag curves intersect at time point b3, and the drag at the wall for Case 2 is larger than that for Case 1 after time point b3 in Figure 9. When the bubble goes through the monitoring point at time point b4, the drag at the wall shows little change. One possible reason for this phenomenon is that when the bubbles are introduced into the fully developed water flow with zero speed, it takes time for them to join the flow. It is not difficult to understand that a larger bubble takes more time and consequently results in different time delays for Cases 1 and 2. The velocity distribution and profiles near the monitoring point are shown in Figure 11. It is clear that velocity gradient of Case 2 is slightly bigger than that of Case 1. In addition, both the velocity contours and velocity profiles around the bubble at three different positions are shown in Figure 12. As the bubble moves along the channel, the velocity contour and velocity profiles at x = 50 mm change significantly compared to those at x = 20 mm, while they change little as the bubble moves from x = 50 mm to x = 80 mm as shown in Figure 12. This further proves that the two-phase flow is fully developed before it reaches the monitoring point (x = 50 mm).

3.2. Effects of Surface Tension

Surface tension usually plays a certain role in multiphase interface. Oishi and Murai [33] found that bubbles have rotating effects that engage with the two neighboring vortices to reduce drag at high Weber number, which is set as We =ρU2D/σ, where σ is surface tension. To clarify such effect, some cases with different surface tensions, that is, 0.01, 0.072, and 0.15 (corresponding to Weber number varying from 1.67 to 24.96), are taken into account, as listed in Table 1. The bubble deformation ratio is defined as , where is the major diameter of the bubble and is the minor diameter. It is found that the bubble achieves a significant deformation under a surface tension of σ = 0.01 N/m, while it deforms only a little under two other conditions of σ = 0.0728 N/m and 0.15 N/m, as shown in Figure 13. Furthermore, it is observed that surface tension has an important influence on wall shear stress at the monitoring point in Figure 14, and the wall shear stresses of Cases 2 and 4 are bigger than that of Case 3. The variation trend of drag at the wall accords with that of the wall shear stress, as shown in Table 1. In order to analyze the mechanism for this phenomenon, three typical time points (c1, c2, and c3) are selected for further investigation, as shown in Figure 14. The velocity gradient at time point c1 is bigger than that at other time points and the velocity at time point c3 is the lowest one, as shown in Figure 15. The effects of bubbles on flow structure at the above three time points are shown in Figure 16. When the bubble just passes the monitoring point, the wall shear stress reaches its maximum value. Furthermore, the disturbance to the flow field is more significant under the surface tension of 0.15 than that of the other two cases. This effect is slightest for Case 3, which also confirms the variations of wall shear stresses at different time points. The reason for this phenomenon could be that the interaction of a bubble and the fluid is neutralized by the bubble deformation for Case 3. To reveal more details, a variation of turbulent kinetic energy due to the change in surface tension is shown in Figure 17. The variation trend is similar to the wall shear stress at the monitoring point. With the comparison between cases of presence and absence of bubbles, it is found that the change of turbulent kinetic energy around the bubble is obvious for the surface tension of 0.15 but there is little change when the surface tension is 0.01.

3.3. Effects of Bubbles Number

The effects of bubbles number on the wall shear stress at the monitoring point are depicted in Figure 18, corresponding to Cases 2 and 10 in Table 1. It is found that wall shear stress at the monitoring point increases with the increase of the number of bubbles. Due to the enhanced action with multiple bubbles, the peak value of the wall shear stress with two bubbles is more than twice the peak value with a single bubble, as shown in Figure 18. The wall shear stress achieves a peak value at about 0.005 s after bubbles pass the monitoring point. The wall shear stress variation is mainly influenced by the presence or absence of bubbles: they change the local flow structures. The velocity gradients near the monitoring point for Cases 2 and 10 have a similar tendency towards variation with the wall shear stress, as shown in Figure 19. As a result, drag at the wall for Case 10 is bigger than that for Case 2, as shown in Table 1.

3.4. Effects of the Initial Position of a Single Bubble

To study the influence of initial bubble position on the change in drag at the wall within a bubbly flow, simulations on a single bubble and two bubbles moving with water in the channel are conducted, corresponding to the case list in Table 1. For Case 2 and Cases 5 to 8 for a single bubble, as the distance of the bubble to the wall dh is lowered, and the interaction between the bubble and the wall is enhanced and leads to more significant perturbations of wall shear stress, as demonstrated in Figure 20. The wall shear stress at the monitoring point increases rapidly after reaching a low point and then decreases quickly to initial levels for Cases 5 to 8. Furthermore, the amplitude of the variation gradually increases when the distance of the bubble to the wall decreases. To uncover the interesting trend in this variation, Case 8 is chosen for further study due to its remarkable variation of wall shear stress.

To analyze the degree of influence of the velocity gradient at different bubble positions, nine typical time points from Case 8 are selected for discussion, as shown in Figure 21. The velocity distributions for different time points are shown in Figure 22. It is obvious that the velocity gradient has a maximum value at time point d9 when the wall shear stress reaches its peak value. Also, the gradient has a minimum value at time point d7 when the shear wall stress reaches the valley. The reason for this situation is the fact that the bubble is located just below the monitoring point at time point d7 as shown in Figure 23, and the flow structure is disturbed by the bubble so the velocity gradient is affected near the monitoring point. Briefly, the trends of the velocity gradients and wall shear stresses are similar. Furthermore, as the distance between the bubble and the wall decreases, drag at the wall increases, as shown in Table 1. In addition, it is clear that the turbulent kinetic energy around the bubble is bigger than those at other positions for Case 2 and Cases 5 to 8 as shown in Figure 24. The change of turbulent kinetic energy supports the trends of wall shear stress as well.

3.5. Effects of Bubbles’ Initial Arrangement

To investigate the effect of the bubbles’ initial arrangement on the flow characteristics, numerical simulation are conducted with two bubbles arranged with the same distance of 2.5 mm vertically (Case 11) and horizontally (Case 10) respectively. The wall shear stresses at the monitoring point show obvious fluctuations, and several time points are selected to analyze the cause of this interesting phenomenon as shown in Figure 25. For Case 11, the wall shear stress at the monitoring point has two peak values and reaches a minimum value at timepoint e3 (Figure 25). The reason for this phenomenon is the fact that the velocity profile in the y direction is non-uniform due to the existence of the wall boundary layer, resulting in a time difference when the two bubbles pass the monitoring point. As a result, the relative distance between the two bubbles in the x direction grows gradually (Figure 26). Compared with those at time points e1 and e3, the velocity gradients at the other three time points are much bigger, as shown in Figure 27. Furthermore, the velocity gradient at time point e4 is slightly bigger than that at time point e5, and both of them are distinctly bigger than that at time point e2. This also explains the variation of wall shear stress and drag at the wall (Table 1) for Cases 10 and 11.

Flow characteristics near the monitoring point are applied to detail the above phenomena. The positions of bubbles at different time points and the velocity distribution near the monitoring point are shown in Figure 26. It is clear that bubble 1 has just passed the monitoring point at time point e2, and at that time bubble 2 has not arrived yet. It takes time for the flow to respond to bubbles, and the velocity gradient fluctuation at the monitoring point at that time point is mainly affected by bubble 1. At time point e3, the influence of the two bubbles is coincidentally neutralized, leading to little fluctuation of the wall shear stress. As shown in Figure 26, at time point e4, both bubbles have passed the monitoring point and bubble 2 dominates the influence on wall shear stress at the monitoring point. Since bubble 2 is closer to the wall, the wall shear stress at time point e4 is much bigger than that at time point e2. This phenomenon has been verified in Section 3.4 as well. For Case 10, the superposition of the influence of the two bubbles makes the peak shear stress value appear at time point e5.

4. Conclusions

In this paper, the bubbly flow in a two-dimension channel is numerically investigated using ANSYS Fluent. The bubble-liquid interfaces are tracked with the VOF method, providing a way to clarify the bubble’s effects on the flow structures and wall shear stresses under various bubble parameters such as bubble size, surface tension, bubbles number, initial positions of a single bubble and initial arrangements of multiple bubbles. Through the present investigation with CFD, the following conclusions can be drawn:(1)The effect of larger bubbles on the drag at the wall is more significant than that of smaller ones, and both bubbles for the two cases result in a drag increase, in accordance with the studies summarized by Murai [7] as well.(2)Higher bubble-liquid interfacial surface tension causes more drastic changes in wall shear stress and turbulent kinetic energy around the bubble, while lower surface tension causes a more pronounced deformation of the bubble and has little effect on the flow structures. Also, study of the effects of the bubble numbers reveals that more bubbles leads to greater influence on the wall shear stress.(3)The change of turbulent kinetic energy reflects the trend of wall shear stresses. It is concluded that shorter distances between bubbles and the wall can enhance the flow interaction and lead to more significant flow perturbations due to effects from the bubbles from the wall shear stress and viscous drag.

Nomenclature

u:Velocity
μ:Viscosity of water
k:Turbulence kinetic energy
ɛ:Turbulence dissipation rate
:Generation of turbulence kinetic energy due to buoyancy
:Contribution of the fluctuating dilatation to the overall dissipation rate
:Turbulent Prandtl numbers for k
:Turbulent Prandtl numbers for ɛ
:Generation of turbulence kinetic energy due to the mean velocity gradients
:Turbulent viscosity
:Volume fraction
:Mass transfer
ρ:Water density
Re:Reynolds number
:Channel height
:The maximum of velocity along the x direction
τ:Wall shear stress.

Data Availability

All other data in this paper were generated at Fishery Machinery and Instrument Research Institute, Key Laboratory of Ocean Fishing Vessel and Equipment, Joint Research Laboratory for Open Sea Fishery Engineering, Pilot National Laboratory for Marine Science and Technology, Nanjing University of Science and Technology, and State Key Laboratory of Ocean Engineering. Derived data supporting the findings of this study are available from the corresponding author at [email protected] upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was financially supported by the National Key R&D Program of China under Project no. 2019YFD0901504, the National Natural Science Foundation of China (NSFC) under Project no. 51609115, the Fundamental Research Funds for the Central Universities under Project no. 30918012201, and the Open Fund of State Key Laboratory of Ocean Engineering under Project no. 1818.