On the Influences of Air Bubbles on Water Flow in a Two-Dimensional Channel

Fishery Machinery and Instrument Research Institute, Chinese Academy of Fishery Science, Shanghai 200092, China Key Laboratory of Ocean Fishing Vessel and Equipment, Ministry of Agriculture, Shanghai 200092, China Joint Laboratory for Open Sea Fishery Engineering, Pilot National Laboratory for Marine Science and Technology (Qingdao), Shandong, Qingdao 266237, China School of Energy and Power Engineering, Nanjing University of Science and Technology, Nanjing 210094, China State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200240, China


Introduction
e use of air bubbles is one of the most effective approaches to reduce drag for ships. e 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 [1][2][3]. e maritime shipping industry is thus interested in employing this dragreduction technology to reduce energy use and greenhouse gas emissions [4,5]. e 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 [10][11][12][13]. 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 [17][18][19][20][21][22]. Hence, a great deal of research has been carried out in the recent years through numerical simulations [23][24][25][26]. 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. e prior studies in the area of drag reduction with air bubbles usually focused on bubbly flow with a mass of bubbles [34][35][36]. 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.
erefore, 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.

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. e standard k-ε turbulence model is used to demonstrate turbulent flow near the channel's boundaries at the walls. e equation for conservation of mass is given as e 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. e standard k-ε model is based on equations describing the turbulence kinetic energy k and turbulence dissipation rate ε. e transport equation for the turbulent kinetic energy k is [37] e transport equation for the turbulence dissipation rate ε is where G b is the generation of turbulence kinetic energy due to buoyancy, Y M is the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, and C 1ε , C 2ε , and C 3ε are constants. Parameters σ k and σ ε are the turbulent Prandtl numbers for k and ε, respectively. Terms S k and S ε are user-defined source terms, and the term G k is the generation of turbulence kinetic energy due to the mean velocity gradients and can be defined as e turbulent viscosity, μ t , is determined by combining k and ε as follows: where C μ is a constant. e 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 α q of the q th phase to values between 0 and 1. e momentum equation of the q th phase takes the following form where _ m qp is the mass transfer from phase q to phase p and _ m pq is the mass transfer from phase p to phase q. e 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. e pressure drop across the surface depends upon the surface tension coefficient, and the surface curvature is measured by two radii in orthogonal directions. e ANSYS ICEM commercial software package is used for geometry modeling and mesh generation. e 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 twodimension turbulent flow is achieved. e water flow simulation is conducted with parameters as follows: the water density is ρ � 998.2 kg/m 3 , the dynamic viscosity is μ � 0.001003 N·s/m 2 , the flow velocity is u � 0.5 m/s, and the Reynolds number Re � ρuh/μ � 4976 is based on the channel height h. Typical velocity profiles at the end of every simulation period are presented in Figure 3(a), and the velocity profiles between the 3 rd and 4 th differs a little (Figure 3(b)). A more clear illustration of the velocity profiles for the 4 th 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 U max 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.

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. e 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 τ � μ(dU x /dy), in which μ is dynamic viscosity of the fluid, U x 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. e wall shear stress decreases rapidly after reaching a peak point (t > 0.16 s) and then increases quickly to another peak point.
at 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 a 1 and a 3 are bigger than those at other points. is makes the wall shear stress at time points a 1 and a 3 bigger than that without bubbles. As shown in Figure 8, the velocity gradients at time points a 1 , a 2 , and a 3 have little difference, while the wall shear stress at time point a 2 is much less than that of the other two time points. e 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.

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. e 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 b 1 for Case 1, and at time point b 2 for Case 2, respectively. Furthermore, an interesting phenomenon occurs that the drag curves   intersect at time point b 3 , and the drag at the wall for Case 2 is larger than that for Case 1 after time point b 3 in Figure 9.

Mathematical Problems in Engineering
When the bubble goes through the monitoring point at time point b 4 , 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. e 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. is further proves that the two-phase flow is fully developed before it reaches the monitoring point (x = 50 mm).  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. e 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 (c 1 , c 2 , and c 3 ) are selected for further investigation, as shown in Figure 14. e velocity gradient at time point c 1 is bigger than that at other time points and the velocity at time point c 3 is the lowest one, as shown in Figure 15. e 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. is effect is slightest for Case 3, which also confirms the variations of wall shear stresses at different time points. e 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. e 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.

Effects of Bubbles Number.
e 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. e wall shear stress achieves a peak value at about 0.005 s after bubbles pass the monitoring point. e wall shear stress variation is mainly influenced by the presence or absence of bubbles: they change the local flow structures. e 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.

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. e 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. e 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 d 9 when the wall shear   stress reaches its peak value. Also, the gradient has a minimum value at time point d 7 when the shear wall stress reaches the valley. e reason for this situation is the fact that the bubble is located just below the monitoring point at time point d 7 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. e change of turbulent kinetic energy supports the trends of wall shear stress as well.

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. e 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 e 3 (Figure 25). e 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 e 1 and e 3 , the velocity gradients at the other three time points are much bigger, as shown in Figure 27. Furthermore, the velocity gradient at time point e 4 is slightly bigger than that at time point e 5 , and both of them are distinctly bigger than that at time point e 2 . is also explains the variation of wall shear stress and drag at the wall (Table 1) for Cases 10 and 11.        Mathematical Problems in Engineering 13 in Figure 26, at time point e 4 , 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 e 4 is much bigger than that at time point e 2 . is 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 e 5 .

Conclusions
In this paper, the bubbly flow in a two-dimension channel is numerically investigated using ANSYS Fluent. e bubbleliquid 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. rough the present investigation with CFD, the following conclusions can be drawn: (1) e 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) e 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 G b : Generation of turbulence kinetic energy due to buoyancy Y M : Contribution of the fluctuating dilatation to the overall dissipation rate σ k : Turbulent Prandtl numbers for k σ ε : Turbulent Prandtl numbers for ε G k : Generation of turbulence kinetic energy due to the mean velocity gradients μ t : Turbulent viscosity α q : Volume fraction _ m: Mass transfer ρ: Water density Re: Reynolds number h: Channel height U max : e maximum of velocity along the x direction τ: Wall shear stress.

Data Availability
All other data in this paper were generated at Fishery

Conflicts of Interest
e authors declare that they have no conflicts of interest.