# Numerical Simulation of Bubble Free Rise after Sudden Contraction Using the Front-Tracking Method

Authors: Ying Zhang; Min Lu; Wenqiang Shang; Zhen Xia; Liang Zeng; Peisheng Li
Journal: International Journal of Photoenergy
 (2017)

---

## Abstract

Based on the front-tracking method (FTM), the movement of a single bubble that rose freely in a transverse ridged tube was simulated to analyze the influence of a contractive channel on the movement of bubbles. The influence of a symmetric contractive channel on the shape, speed, and trajectory of the bubbles was analyzed by contrasting the movement with bubbles in a noncontractive channel. As the research indicates, the bubbles became more flat when they move close to the contractive section of the channel, and the bubbles become less flat when passing through the contractive section. This effect becomes more obvious with an increase in the contractive degree of the channel. The symmetric contractive channel can make the bubbles first decelerate and later accelerate, and this effect is deeply affected by Reynolds number (Re) and Eötvös number (Eo).


---

## Body

## 1. Introduction

The corrugated tube has been the focus of scholars because the internal structure of the pipe can disturb the movement of the fluid and increase the surface area of the pipe; thereby, it can increase the heat transfer of the pipe. So it is widely used in various fields such as the solar hot water system because there are a number of solar collector tube and heat transfer tubes in this system. In order to improve the heat absorption and heat transfer efficiency, different types of corrugated tubes have been studied. Vicente et al. [1] studied heat transfer of the threaded pipe under different sizes by experimentation, he shows the guidelines to choose which roughness geometry offers the best performance for specific flow conditions. García et al. [2] analyzed the thermal-hydraulic behavior of three types of enhancement technique based on artificial roughness, and the results show that shape of the artificial roughness exerts a greater influence on the pressure drop characteristics than on the heat transfer augmentation. Saha [3] has studied the pressure drop characteristics of pipe with internal transverse rib turbulators on two opposite surfaces and with wire-coil inserts. Pethkool et al. [4] experimentally studied the heat transfer performance of helically corrugated tube and got the empirical formula. Kareem et al. [5] studied the three-start spirally corrugated tube by experiment and simulation, he said that this geometry with a creative spiral corrugation profile can improve the heat transfer significantly with reasonable increase in friction factor. It is obvious that scholars have made a lot of research on corrugated tubes in different structural forms and got some good results. However, most of these studies have been studied for the effects of different pipe shapes on heat transfer and pressure drop with a single fluid, while the study of multiphase flow in the tube has been rarely reported.In practical cases, such as solar collectors and heat exchangers, the temperature of water in the tube is generally higher than the natural temperature. The gas in the water will be released and the bubbles will form. The interaction of the bubble with the contraction structure can also cause negligible disturbance to the flow in the pipe and affect its heat transfer. Therefore, we believe that the study of bubble flow in the corrugated tube has a certain reference value for improving the heat transfer of the heat exchanger, reducing the pressure drop of the helically corrugated tube and obtaining the empirical formula of thermal performance.The method of numerical simulation of bubble motion has been well developed [6–11]. In order to describe the movement characteristics of the bubbles in the corrugated tube intuitively, it is convenient to study it with the simulation method. FTM is a classic method used in numerical simulations that has been widely used in multiphase analyses [12–14]. It can accurately capture the complex topological changes of the moving interface, and it is very important in tracking the phase interface during the research process of multiphase flow.Therefore, in this paper, the FTM was used to simulate the free rise of a single bubble in the transverse ridged tube. However, the cross-sectional shape of the transverse ridged tube is not a regular rectangle like an ordinary pipe, so we set the contraction channel by modifying the boundary position. And we analyzed the coupling effect between the bubble and the solid wall of the channel and studied the influences of sudden contraction on the shape, velocity, and trajectory of the bubble. It is hoped that the analysis of the flow field under coupling effect between the bubble and the solid wall of the channel will provide a theoretical help for the design of the heat transfer pipe in a solar hot water system, so that the internal flow of the transverse ridged tube can be more reasonable and achieve better heat transfer effect.


## 2. Mathematical Model, Numerical Methods, and Physical Model

### 2.1. Tracking and Processing the N-S Equations at the Phase Interface

The surface tension is concentrated on the interface when the control equation is applied to the whole calculation area. We consider introducing a source term and is applied by multiplying the volume force and theδ function. The δ function is nonzero only on the interface. Using the amendments above, the two-dimensional incompressible flow momentum equation can be expressed as
(1)∂∂tρu+∇⋅ρuu=−∇p+ρg+∇⋅μ∇u+∇uT+∫sσκnδx−xfds,where δx−xf is the two-dimensional Dirac function, xf is the position of the surface, the integral is along the interface S, n is the interface normal, σ is the interface tension coefficient, κ is the secondary interface average curvature, and ρ and μ are noncontinuous density field and viscosity field.To solve the equation above, we split (4) and remove the pressure introducing the temporary speed u∗. Thus, an equation containing only the A convection term and the D diffusion term is obtained:
(2)u∗−unΔt=−An+g+1ρnDn.By solving the equation, we obtain the temporary velocity field without pressure, and then add the pressure term to obtain the general velocity discrete equation:(3)un+1−u∗Δt=−∇hpρn.In the formula,∇h is a discrete form of the Hamiltonian operator whose step size is h and because of the speed it needs to meet the discrete form of mass conservation:
(4)∇hun+1=0.Then the pressure Poisson equation is obtained:(5)∇h⋅1ρn∇hp=1Δt∇h⋅u∗.The CSF (continuous surface force method) method given by Brackbill [15] was used to calculate the surface tension.
(6)Fσ=∫sσTkn+∂σT∂stx−xfds.Take the surface tension on the unit line to study, for the two-dimensional flow is
(7)κn=∂t∂s.Get the surface tension on the unit interface element:
(8)δFσl=∫△sσ∂t∂s+∂σT∂stds=σtl+1/2−σtl−1/2.Here, t is the tangent vector on the interface.


### 2.2. Reconstruction of Density and Viscosity

For immiscible and incompressible fluids, the fluid maintains its characteristics at both sides of the interface, so there is a step in the physical properties at the interface. When the interface moves, the density and viscosity distribution of the interface also change, so the Heaviside function is introduced as an indicator to characterize the change:(9)Hφ=1,φ>α121+φα+1ππφα0,φ<−α,,φ≤αwhere φ is the distance from the given point to the interface, and 2α is the absolute value of the thickness of the transition zone between the two fluids. So the density and viscosity distribution is
(10)ψφ=Hφψ1+1−Hφψ2.Here, ψ1,ψ2 are the density and viscosity of the two fluids.


### 2.3. Interface Movement

The moving interface is typically combined with a fixed mesh, and the information exchange between the interface and the grid is achieved by the area weight function. With the bilinear interpolation method, the abrupt change density and the surface tension of the interface are changed from the interface to the fixed grid. Finally, we use a simple first-order time integral to obtain the interface position of the next time.(11)xfn+1=xfn+Δtufn.Here, xfn is the position of surface and u is velocity.


### 2.4. Physical Model

In the actual production industry, the geometric conditions of the contractive channel are varied, and to facilitate the research, only the unit shown in Figure1 is extracted as the object of this paper and the upper and lower borders are set to periodic boundary conditions to satisfy the structure of the transverse ridged tube. Since the three-dimensional model will spend much time, so we use the two-dimensional model. The contractive ratio is defined as E = L/M and is used to describe the contractive degree. D is the diameter of the bubble.Figure 1
Physical model.Re is defined as(12)Re=ρluDυl.Eo is defined as(13)Eo=ρlgD2σ.Mo is defined as(14)Mo=ρl−ρggυl4σ3ρl2.Here,ρl is the density of liquid phase, υl is the viscosity of liquid phase, σ is the coefficient of the surface tension, ρg is the density of the gas phase, and g is the acceleration due to gravity.


### 2.5. Initial Conditions

GivenD/L = 0.2, the center coordinates of the bubble are (0.5 L, 0.5 L). In this paper, the ratio of the two-phase density is 5 : 1, and the wetting condition (contact angle effect) is not considered. The calculating area is the 12 rectangular region, and the grid resolution is 256 × 512 (including the contraction region, which is not involved in the calculation). The X direction was defined using a no-slip boundary condition, and the Y direction was defined using periodic boundary conditions. All the results in this paper are simulated in two-dimensional space.


## 2.1. Tracking and Processing the N-S Equations at the Phase Interface

The surface tension is concentrated on the interface when the control equation is applied to the whole calculation area. We consider introducing a source term and is applied by multiplying the volume force and theδ function. The δ function is nonzero only on the interface. Using the amendments above, the two-dimensional incompressible flow momentum equation can be expressed as
(1)∂∂tρu+∇⋅ρuu=−∇p+ρg+∇⋅μ∇u+∇uT+∫sσκnδx−xfds,where δx−xf is the two-dimensional Dirac function, xf is the position of the surface, the integral is along the interface S, n is the interface normal, σ is the interface tension coefficient, κ is the secondary interface average curvature, and ρ and μ are noncontinuous density field and viscosity field.To solve the equation above, we split (4) and remove the pressure introducing the temporary speed u∗. Thus, an equation containing only the A convection term and the D diffusion term is obtained:
(2)u∗−unΔt=−An+g+1ρnDn.By solving the equation, we obtain the temporary velocity field without pressure, and then add the pressure term to obtain the general velocity discrete equation:(3)un+1−u∗Δt=−∇hpρn.In the formula,∇h is a discrete form of the Hamiltonian operator whose step size is h and because of the speed it needs to meet the discrete form of mass conservation:
(4)∇hun+1=0.Then the pressure Poisson equation is obtained:(5)∇h⋅1ρn∇hp=1Δt∇h⋅u∗.The CSF (continuous surface force method) method given by Brackbill [15] was used to calculate the surface tension.
(6)Fσ=∫sσTkn+∂σT∂stx−xfds.Take the surface tension on the unit line to study, for the two-dimensional flow is
(7)κn=∂t∂s.Get the surface tension on the unit interface element:
(8)δFσl=∫△sσ∂t∂s+∂σT∂stds=σtl+1/2−σtl−1/2.Here, t is the tangent vector on the interface.


## 2.2. Reconstruction of Density and Viscosity

For immiscible and incompressible fluids, the fluid maintains its characteristics at both sides of the interface, so there is a step in the physical properties at the interface. When the interface moves, the density and viscosity distribution of the interface also change, so the Heaviside function is introduced as an indicator to characterize the change:(9)Hφ=1,φ>α121+φα+1ππφα0,φ<−α,,φ≤αwhere φ is the distance from the given point to the interface, and 2α is the absolute value of the thickness of the transition zone between the two fluids. So the density and viscosity distribution is
(10)ψφ=Hφψ1+1−Hφψ2.Here, ψ1,ψ2 are the density and viscosity of the two fluids.


## 2.3. Interface Movement

The moving interface is typically combined with a fixed mesh, and the information exchange between the interface and the grid is achieved by the area weight function. With the bilinear interpolation method, the abrupt change density and the surface tension of the interface are changed from the interface to the fixed grid. Finally, we use a simple first-order time integral to obtain the interface position of the next time.(11)xfn+1=xfn+Δtufn.Here, xfn is the position of surface and u is velocity.


## 2.4. Physical Model

In the actual production industry, the geometric conditions of the contractive channel are varied, and to facilitate the research, only the unit shown in Figure1 is extracted as the object of this paper and the upper and lower borders are set to periodic boundary conditions to satisfy the structure of the transverse ridged tube. Since the three-dimensional model will spend much time, so we use the two-dimensional model. The contractive ratio is defined as E = L/M and is used to describe the contractive degree. D is the diameter of the bubble.Figure 1
Physical model.Re is defined as(12)Re=ρluDυl.Eo is defined as(13)Eo=ρlgD2σ.Mo is defined as(14)Mo=ρl−ρggυl4σ3ρl2.Here,ρl is the density of liquid phase, υl is the viscosity of liquid phase, σ is the coefficient of the surface tension, ρg is the density of the gas phase, and g is the acceleration due to gravity.


## 2.5. Initial Conditions

GivenD/L = 0.2, the center coordinates of the bubble are (0.5 L, 0.5 L). In this paper, the ratio of the two-phase density is 5 : 1, and the wetting condition (contact angle effect) is not considered. The calculating area is the 12 rectangular region, and the grid resolution is 256 × 512 (including the contraction region, which is not involved in the calculation). The X direction was defined using a no-slip boundary condition, and the Y direction was defined using periodic boundary conditions. All the results in this paper are simulated in two-dimensional space.


## 3. Calculation Results and Discussion

Hua and Lou simulated several typical bubble morphology in Bhaga and Weber [16] with a two-dimensional model, and they compared the Re number between the experiment and simulation and results from simulations agree with those of experiments very well within 10% difference. So, we simulated one of the case (Eo = 116, Re = 13.95) from Hua and Lou [12] and the result is shown in Figure 2. We can find that the results of our model are very close to the result reported in [12].Figure 2
The evolution of a rising bubble with Eo = 116 and Re = 13.95.
(a)
Present work
(b)
Reported in [12]In the process of numerical simulation, the grid density has a significant effect on the calculation results. In order to ensure that there is no correlation between the grid density and the calculated results, we have carried out the four cases of the grid density of 64 × 128, 128 × 256, 256 × 512, 512 × 1024 to simulate the bubble rises in the contractive channel and we extract the position of bubble centroid at different time. The results are shown in Figure3. It can be seen that the difference between the numerical results is not large at different grid densities. The calculation results of the other three cases are very close except that the grid density is 64 × 128. To meet the accuracy of the calculation results and to save the calculation time as much as possible, we select the 256 × 512 grid density to calculate the object of our study.Figure 3
Grid convergence test.


## 4. Effect of a Symmetric Contractive Channel on Bubble Deformation

In this part, we primarily analyze the effect of the symmetric contractive channel on the bubble’s deformation. Figure4 shows the movements of bubbles with a contractive ratio of 1, 1.18, 1.45, 1.88, 2.2, 2.67, 4.57, and 16 under the same situation (Re = 28.2, Eo = 16, and Mo = 0.01). In the noncontractive channel (E = 1), the bubble changes from circular into a hat, and then, its shape becomes stable.Figure 4
Free rise of the bubbles in different contractive channels (Re = 28.2, Eo = 16, and Mo = 0.01).Because the contractive degree is small (E = 1.18, 1.45), the deformation of the bubble in this channel does not differ significantly from a bubble in the noncontractive channel. When the contractive ratio is 1.88, 2.2, 2.67, and 4.57 and before entering into the contractive portion of the channel, the shape of the bubble is more flat compared with the shape of bubbles in the noncontractive channel. Even when the contractive ratio is greater (E = 2.67, 4.57), the bubble becomes wider than that section of the contraction. When the bubble just enters the contractive portion of the channel, the middle portion of the bubble is raised up. When the bubbles enter into the contractive portion of the channel, the distance of the vertical wall on both sides of the bubble becomes smaller, and the shape of the bubble is not as flat as before. When the contractive ratio is 16, the bubble cannot pass through the contraction, and its shape becomes more flat.To analyze the deformation of the bubble more directly, the deformation of the bubble was measured using the aspect ratio. Figure5 shows the variation of aspect ratio with the centroid height of the bubble during the bubble motion given the condition in Figure 4. In a sudden contractive channel, the deformation of the bubble is primarily affected by the horizontal wall of the contraction above the bubble and the vertical wall in the contractive channel. When the bubble is close to the contractive part of the channel, the bubble is primarily affected by the horizontal wall at the top of the bubble. When the bubble enters into the contractive portion, it is primarily affected by the vertical wall on both sides of the bubble.Figure 5
Aspect ratio of the bubbles with different contractive ratios (Re = 28.2, Eo = 16, and Mo = 0.01).When the contractive ratio is 1.18 and 1.45 and the bubble is near the contractive portion, because the horizontal wall above the bubble is small, there is a slight effect on the shape of the bubble, and there is only a slight increase in the aspect ratio. When the bubbles remain in the contractive part, the aspect ratio of the bubble is smaller because of the effect of the vertical wall.When the contractive ratio is 1.88, 2.2, 2.67, and 4.57 and the bubbles are close to the contractive part, the bubbles are primarily affected by the horizontal wall above the bubble. Therefore, when the aspect ratio of the bubble is compared with the motion in the noncontractive channel, the aspect ratio of the bubble is greater, and with an increase in the contractive degree, a larger horizontal wall area results in a stronger effect. For example, when the contractive ratio is 4.57, the aspect ratio achieved the maximum value of 5.2, which is 1.85 times greater than the motion of the bubbles in the noncontractive channel at the same time. When the bubbles enter the contractive portion, the shape of the bubble will be affected by the vertical wall on both sides of the channel, and the aspect ratio of the bubble becomes smaller. When the bubbles come out from the contractive portion, the aspect ratio of the bubble will become larger.When the contractive ratio is 16, because the contractive degree of the channel is too large, the channel is too narrow and the bubbles cannot enter the contractive portion. The bubble is primarily affected by the horizontal wall on the top of the bubble. The aspect ratio of the bubble increases constantly, and the maximum aspect ratio is close to 8 and would then decrease, but the bubble still cannot pass through the channel.


## 5. Effect of the Symmetric Contractive Channel on the Velocity of the Bubbles

When the bubbles move upward from the bottom under the action of gravity, the channel which contracted suddenly will affect the speed of the bubbles, and this effect is primarily influenced by Eo and Re. In this paper, we only analyze the situation in which the bubbles can pass through the contractive channel, and we do not analyze the situation in which the bubble cannot pass through the channel when the contractive ratio is too large, as shown in Figure4 (E = 16). In this section, we study the case which the contractive ratio is 2.2, because in this case the effect of the wall becomes obvious and the bubble can go through the channel smoothly. We compared the motion of the bubble with that of the noncontractive channel. To react to the speed and the deformation of the bubble, this paper analyzed the bottom and the top velocity of the bubble.The influence of the contractive channel on the velocity of the bubble was primarily related to Re and Eo. Figure6 is the change chart for the bottom and the top velocity of the bubble in the case of Re = 70 and Eo = 8, and Figure 7 is in the case of Re = 70 and Eo = 3.2.Figure 6
The free rising velocity of bubbles in the contractive and noncontractive channels (Re = 70, Mo = 8 × 10−5, and Eo = 8).Figure 7
The velocity of the bubbles in the contractive and noncontractive channels (Re = 70, Mo = 8 × 10−5, and Eo = 3.2).As Figure6 shows, during the early process of the bubble rising, the bubble was far away from the section of the contraction. The velocity of the bubble was less affected, so the speed of the bubble in the contraction channel was the same as the speed in the noncontractive channel. When the bubble was close to the contractive channel, it was affected by the horizontal wall. Therefore, the aspect ratio of the bubble becomes greater, the resistance to the bubble rising becomes greater, and the velocity of the bubble becomes smaller. When the bubble enters into the contraction, the aspect ratio of the bubble decreases, the resistance to the bubble rising was also reduced, and the velocity of the bubble increased. When the bubble was coming out of the contractive portion, the velocity of the bubble decreases.When the bubble is located near the contractive part, the contractive channel decelerates the bubble, and when the bubble enters the contraction, the contraction accelerates the bubble. Figure8 shows the height-time diagrams of the bubble from which we can observe the influence of the contraction clearly under the combined influence of the “deceleration” and the “acceleration.” The time it takes the bubble to go through the contractive and noncontractive channels is approximately the same.Figure 8
The centroid height of the bubble in the contractive and noncontractive channels (Re = 70, Mo = 8 × 10−5, and Eo = 8).It can be observed from Figure7, in the case of a larger Eo value, the contractive channel also has the effect of “decelerating” and “accelerating” the bubbles compared with the front situation. The contractive channel had a weak influence on the acceleration of the bubbles. Therefore, as shown in Figure 9, it takes the bubble a longer time to pass through the contractive channel compared with the bubble going through the noncontractive channel.Figure 9
The centroid height of the bubble in the contractive and noncontractive channels (Re = 70, Mo = 1 × 10−5, and Eo = 3.2).The first case shows the influence of the contractive channel on the velocity of the bubble for the two different values of Eo. This influence on the velocity is strongly affected by Eo and Re. According to the analysis above and due to the influences of the contractive channel, the velocity of the bubbles decreased when the bubble nears the section of the contractive channel, and then when the bubble enters the contraction, the velocity of the bubble increased. To measure the comprehensive influence of the contraction channel on the velocity of the bubble, we defined the average velocity ratio as, in the same case, the average velocity of a bubble that moves from the initial position to a position that is completely through the channel; the contraction channel has a constant cross-sectional area. Figure10 shows the distribution of the average velocity ratio of the bubbles when Re and Eo are different.Figure 10
Relationship between the ratio of average velocity and Re and Eo.As seen in Figure10, with an increase in Re, the ratio of the average velocity increased given the same value of Eo. According to the research from Clift [1], with an increase in Re, the shape of the bubbles has the tendency to become more flat. A flat bubble shape results in a greater resistance to rising. However, when the bubbles pass through the contraction channel, the bubbles are squeezed due to the effect of the wall making the aspect ratio of the bubble decrease, and the resistance of the bubble to rise has a decreasing trend. Therefore, the ratio of the average velocity increases with an increase in Re. Figure 10 shows that a decrease in Eo results in the average velocity ratio of the bubble being reduced. A smaller Eo results in a more difficult deformation of the bubble, so the effect that the vertical wall makes the aspect ratio of the bubble become smaller is not obvious. That is, the resistance of the bubble to rise is not obvious. A smaller Eo results in the bubble having a smaller average velocity ratio.


## 6. Conclusion

(1)
The bubble rises freely in a symmetric contractive channel when the distance between the bubbles and the contractive part of the channel is considerable, and the deformation of the bubble is similar to the situation in which the bubbles move in a noncontractive channel. When the bubble is close to the contractive channel, due to the effect of the wall, the shape of the bubbles will be more flat; that is, the aspect ratio becomes greater. When the bubbles enter the contractive part of the channel, the shape of bubbles will not be as flat as before, and the aspect ratio becomes small. A greater contractive degree of the channel results in the effect becoming more obvious.
(2)
Compared with the noncontractive channel, when the bubble rises freely in a symmetric contractive channel, the contractive channel causes the velocity of the bubble to decrease at first and then increase combining both “deceleration” and “acceleration.” The average velocity of the bubble in the contractive channel may be greater than in the noncontractive channel or it may be smaller. The main influential factors are Re and Eo, and when Re or Eo is greater, the bubble will have an accelerating trend.



---

*Source: 5128345-2017-10-24.xml*