Abstract
Two-phase flow instability may occur in nuclear reactor systems, which is often accompanied by periodic fluctuation in fluid flow rate. In this study, bubble rising and coalescence characteristics under inlet flow pulsation condition are analyzed based on the MPS-MAFL method. To begin with, the single bubble rising behavior under flow pulsation condition was simulated. The simulation results show that the bubble shape and rising velocity fluctuate periodically as same as the inlet flow rate. Additionally, the bubble pairs’ coalescence behavior under flow pulsation condition was simulated and compared with static condition results. It is found that the coalescence time of bubble pairs slightly increased under the pulsation condition, and then the bubbles will continue to pulsate with almost the same period as the inlet flow rate after coalescence. In view of these facts, this study could offer theory support and method basis to a better understanding of the two-phase flow configuration under flow pulsation condition.
1. Introduction
Two-phase flow phenomena could occur in nuclear reactor systems (steam generators, boiling water reactor cores, condensers, etc.), conventional power plants, refrigeration, and chemical and aerospace. In a two-phase system, when mass flow density, void fraction, and pressure drop are coupled, the static flow instability (the flow is drifted by a slight disturbance) and the kinetic instability (constant/variable amplitude flow oscillations at a particular frequency) are collectively referred to as flow instability [1]. These flow instability phenomena are usually accompanied by periodic fluctuation in flow rate. As is well known, two-phase flow instability can degrade the performance of the equipment, which will cause many safety problems; for example, the mechanical forces caused by the flow oscillations will make the components undergo harmful forced mechanical vibrations, and continuous mechanical vibration could lead to fatigue failure of components. In addition, bubble dynamics is an indispensable part of the two-phase flow research and the intrinsic driving force of gas-liquid two-phase flow evolution. Therefore, it is necessary to study the bubble behavior under flow pulsation condition.
In experimental research, a great deal of achievements on bubble dynamics have been made until now. In the bubble coalescence visualization experiment of Sanada et al. [2], the high-speed camera was used to capture the rising and coalescence process of bubbles in the stationary liquid phase, and the microscopic characteristics such as the trajectory, velocity, and interface configuration evolution of the bubble coalescence process were obtained. Liu and Zheng [3] used the method of particle image velocimetry to study the motion behavior of bubbles in a stationary liquid with a rectangular column, and the diversity of bubble rising shapes in different viscosity liquids was finally observed. Sathe et al. [4] used the discrete wavelet transform algorithm to reconstruct the PIV flow field distribution around the bubble, which reduced the interference signal during the PIV flow field acquisition process and further obtained the flow field characteristics like phase drift velocity. Basarova et al. [5] experimentally investigated the ascending motion behavior of spherical bubbles in different concentrations of water ethanol and water propanol and obtained characteristic values such as bubble rising velocity and bubble size, which proved the effect of the concentration of the two solutions on the liquid properties and the gas-liquid interface. Maldonado et al. [6] experimentally studied a single bubble of 2.5 mm diameter rises in different liquids. They established a method for calculating the aspect ratio and the rising velocity of bubble in experiment. Meanwhile, they found that the bubble aspect ratio has a unique relationship with the bubble rise velocity under different conditions. Liu et al. [7] performed experiments to study the rise of bubbles in stationary water and glycerol aqueous solution. Specifically, the experiment detected that the bubble shape has a strong connection and interaction with the terminal velocity. At the same time, the inertial force and the surface tension mainly affect the rising shape of the bubble in water, but in the aqueous glycerin solution, the influence of the viscous force needs to be considered. Raymond and Rosant [8] performed experiments on the movement of moderate deformation bubbles in a viscous liquid mixed with glycerin and water, where the dimensionless Reynolds number and the Weber number range are [,100] and [,5], respectively. The bubble aspect ratio and the bubble terminal velocity, which describe the bubble motion characteristics, were obtained by these experiments and the experimental results were also compared with numerical simulation. Sugrue et al. [9] investigated the influences of channel orientation on bubble dynamics. Besides, the bubble dynamics experiments under different fluid conditions [10–13] have been further conducted.
Due to the uncertainty of the experiment and the limitations of the experimental conditions, numerical simulation and mechanism studies play an increasingly important role in study of bubble dynamics. Tomiyama et al. [14] investigated the bubble terminal velocity in an unlimited large static liquid dominated by surface tension and proposed a theoretical model applicable to predict the terminal velocity of a twisted spherical bubble with a high nondimensional Reynolds number. Based on the balance of all forces acting on the bubble, Wang et al. [15] analyzed the related forces of bubble and illustrated the bubble interface and the shape of the bubble. Moreover, by theoretical analysis, a conclusion was drawn that the mass and heat transfer in the gas-liquid phase interface provides the possibility of continuous motion of the bubble. Gumulya et al. [16] performed a numerical analysis on the rising behavior of bubbles in viscous liquids. On the basis of the three-dimensional VOF method, the correlation between bubble aspect ratio and several dimensionless numbers in different bubble shapes was emphatically studied. Additionally, they also paid attention to the formation and expansion features of bubble wakes. Yu et al. [17] used the new kind of lattice Boltzmann method to simulate a pair of bubbles with different initial positions and initial velocities, and the interactions between bubbles have been successfully obtained, such as attraction and repulsion. Sanada et al. [18] principally studied the effects of viscosity on coalescence of a bubble upon impact with a free surface; it was concluded that both the liquid film and the bubble wake have a significant impact on the bubble coalescence or bounce.
In the 1990s, a meshless method called the Moving Particle Semi-implicit (MPS) method was proposed by Koshizuka and Oka [19]. It is applicable to the simulation of incompressible fluid, especially those with free surface. In order to solve the problem of particle tracking in the presence of inlet and outlet flows in fluid mechanics, Yoon et al. [20] came up with a hybrid method called MPS-MAFL. At this stage, two-phase simulation based on the MPS-MAFL method has been widely used [21–25] because of its outstanding advantages in capturing phase interface.
Under the condition of inlet flow pulsation, the gas-liquid phase interface changes rapidly and dramatically due to the fluctuation of flow velocity. When using the traditional grid method to simulate the bubble motion characteristics under the condition of velocity pulsation, the mesh quality and the stability of the calculation near the phase interface are very challenging. Nowadays, most of the researches on bubble dynamics were carried out under the static condition, and the study of bubble behavior under flow pulsation condition is insufficient. Therefore, it is reasonable and feasible to study the bubble rising and bubble pairs coalescence characteristics under flow pulsation condition by using the predominance of MPS-MAFL in tracking phase interface. Furthermore, this research will also provide a strong theoretical support for understanding the two-phase flow configuration.
2. Numerical Method and Program Validation
2.1. Governing Equations
In MPS-MAFL method, the continuity equation and the Navier-Stokes equation for the incompressible viscous fluid are as follows:where u, , , n, , and are velocity, surface tension coefficient, steam liquid interface curvature, unite vector of phase interface, viscosity, and the motion of a computing point that is adaptively configured during the calculation, respectively. Due to the permission of arbitrary calculation between fully Lagrangian () and Eulerian () calculations, the sharp fluid front is calculated accurately by moving the computing points in Lagrangian coordinates, while the computing point is fixed at the boundaries that are described with Eulerian coordinates [20].
2.2. Particle Interaction Models
In MPS method, the governing equations are solved by particle interaction models such as gradient model, Laplacian model, and divergence model, where a particle interacts with nearby particles covered with a weight function , where is the distance between two particles and is the radius of interaction area. In this study, the kernel function is adopted as weight function, and its expression is as follows:
The gradient and Laplacian of physical quantity on particle are expressed by the summation of physical quantities over its neighboring particles according to the weight function.
The gradient model and Laplacian model are as follows, respectively:where represents spatial dimension and is diffusion coefficient and can be determined as
In the MPS-MAFL method, the divergence model is shown in (7). The velocity divergence between two particles is defined by and the velocity divergence at one particle is calculated from the weighted average of the individual velocity divergences.
2.3. Computational Domain and Boundary Condition
In this study, the calculation domain of the single bubble rising process under flow pulsation condition is a 0.08m wide and 0.47m high two-dimensional region, as shown in Figure 1(a). The initial radius of the bubble is 0.01m and the bubble is located in the center of the channel. The distance between the center of the initial bubble and the inlet of the flow channel is 0.037m. This geometric setting ensures that the influence of the wall can be neglected during bubble motion. The vapor phase property is given based on ideal gas state equation and the liquid area is described by using discretized particles with nonuniform scheme. And the particle number density decreases as the distance from the phase interface increases in the liquid area.
(a) Single bubble
(b) Bubble pairs
As for the boundary condition setting, the left and the right sides of the flow path are set as nonslip boundary conditions, and the upper side of the flow path adopts a free surface boundary condition. The bubble interface is set as free surface, and the pressure values at various points on the interface are obtained by solving the ideal gas state equation. Moreover, at the inlet boundary of the channel, a flow pulsation condition is introduced. The flow at the inlet changes in a sinusoidal cycle over time as shown in (8). The equation contains three variables: the initial flow rate of the liquid phase , the relative amplitude of the flow pulsation , and the flow pulsation period .
The calculation domain of bubble pairs coalescence under flow pulsation condition is set in a similar way, as shown in Figure 1(b). The calculation domain is 0.01m wide and 0.008m high, and the two bubbles are symmetrically distributed. The boundary conditions and particle placement are the same as single bubble conditions.
In the process of bubble dynamics simulation using MPS method, the particle size (the distance between neighboring particles) could affect the simulation efficiency and convergence. Specifically, the optimal particle size could be figured out by sensitivity analysis. In this study, the simulation of single bubble rising in a static fluid was performed with different particle size. As shown in Figure 2, the interfacial particle size was set to 1/8, 1/9, 1/10, and 1/11 of the bubble radius, respectively. It can be seen that the bubble centroid velocity calculation results have good convergence under these particle sizes. At the same time, the shape change during the bubble rising progress was obtained. It is found that the terminal shapes are nearly identical in each case. Therefore, based on the sensitivity analysis results, the interfacial particle size was set to 1/10 of the initial bubble radius to ensure the efficiency of the program simulation and the convergence of the calculation.
2.4. Program Validation
The single bubble rising in static water was simulated by MPS-MAFL, and the simulation results were compared with the experiment of Liu et al. [26]. The physical property parameters of the relevant cases are listed in Table 1. In MPS-MAFL, the calculation conditions and parameters were set as same as the experiment, and the terminal rising velocity of the bubble under different initial bubble diameter conditions was obtained, as shown in Figure 3, where scatter points, respectively, represent experimental data points for two temperatures.
As shown in Figure 3, the calculated bubble terminal velocity value is slightly lower than the experimental one with the error about 15%. One reason is that the lateral movement of the bubbles existing in the experiment had been simplified in the two-dimensional simulations of this study. Within the range of bubble diameter larger than 6 mm, the bubble terminal rising velocity slightly increases with the bubble diameter, which is consistent with the experimental conclusion of Liu et al. [26]. Thus, the overall trend of the simulation results is in reasonable agreement with the experiment.
In the process of single bubble rising in a static liquid, the stable shape of the bubble is another important characteristic that should be concerned. On the basis of previous experimental results, Grace [27, 28] proposed a graphical correlation that predicts the final stable shape of the rising bubble in an infinite static liquid. It can be seen from this diagram that the terminal shape of the bubble can be roughly divided into four regions, namely, Sphere region, Ellipsoid region, Dimples region, and Spherical-cap region. In order to verify the accuracy of the simulation, the following four cases are, respectively, set: Re = 8.854, Eo = 0.392; Re = 88.54, Eo = 3.92; Re = 10, Eo = 100; and Re = 1000, Eo = 300. Moreover, the stable shape in the bubble rising process in each case is obtained, and the stable shape of the bubble is marked at the corresponding position according to the Reynolds number and the Eotvos number (also called Bond number), as shown in Figure 4. In view of these facts, it is completely consistent with the results predicted by Grace.
This study aims to discuss the bubble rising and coalescence characteristics under inlet flow pulsation. Thereby, it is necessary to verify the accuracy of the added inlet flow pulsation condition before the simulation. First, a flow channel only with liquid is set, and its geometric parameters are the same as those for studying single bubble rising under flow pulsation condition, as shown in Figure 5. Then a flow rate that changes periodically with time was added to the inlet of the flow path. The initial flow rate of the liquid phase, the flow pulsation period, and the relative amplitude of the flow pulsation were 0.05 m/s, 0.2s, and 0.7, respectively. Figure 5 shows the velocity values of the particles at the inlet and outlet of the flow channel when the flow stability is reached. As shown in Figure 5, the velocity values of the outlet particles are the same as those of the inlet particles at any of the same moments and they all satisfy the analytical formula of the added flow pulsation condition.
3. Results and Discussions
The Reynolds number (Re) and the Bond number (Bo) are commonly adopted to study the bubble dynamics. The expressions are, respectively, as follows:where is the density of liquid, is liquid viscosity, and is the bubble diameter. represents the surface tension. Through these two dimensionless expressions, the Reynolds number represents the ratio of the inertial force to the viscous force of the bubble, and the Bond number represents the ratio of the gravity of the bubble to its surface tension. During the movement of the bubble, these four forces have a significant effect on the bubble. Besides, the density ratio and viscosity ratio of the gas-liquid two phases also have an impact. In this study, the density ratio and viscosity ratio of liquid phase to gas phase are set as 1000:1 and 100:1 in all bubble simulations, respectively.
3.1. Single Bubble Rising under Flow Pulsation Condition
In this study, the influence of the pulsation condition on the rising motion characteristics of the bubble is investigated. Due to the addition of the pulsating flow in the form of (8), three independent variables are added. Together with the two dimensionless numbers mentioned above, there are actually five independent variables that may affect the motion behavior of the bubble under the inlet flow pulsation condition. Therefore, the control variable method can be used to discuss the bubble deformation process and the surrounding flow field changes under the flow pulsation condition. The specific analysis is as follows.
3.1.1. Single Bubble Deformation Process
The cases of studying the single bubble rising process under flow pulsation condition are listed in Table 2. Initially, a case under static condition is taken as a reference example, while its dimensionless number and other conditions were as same as case in Table 2. And the bubble rising process under static condition was shown in Figure 6(a) with red-dashed line. To sum up, in this case, the bubble will only gradually stabilize during the ascent, and the terminal shape is consistent with the prediction of Grace [27].
(a) Re = 5, Bo = 10, = 0.05m/s, = 0.7, τ = 0.2s/Re = 5, and Bo = 10 (static condition: red-dashed line)
(b) Re = 10, Bo = 10, = 0.05m/s, = 0.7, and τ = 0.2s
By simulating the various cases set in Table 2, it is possible to observe the change of the shape and the surrounding flow field during the bubble rising in the liquid under the periodic variation of the inlet flow rate. Figure 6 shows the change process of the rising bubble under cases and in two cycles (t = 0.4s), and the flow field velocity vector information and pressure distribution in the red frame (around the bubble) are given by Figure 7. The bubble has a wake during the ascent, and there are two symmetric vortices around the bubble. It can also be observed in Figure 6 that there is a significant shape and volume fluctuation in each cycle of the bubble rising process due to the introduction of the inlet periodic change flow rate. This is one of the differences between the bubble deformation processes under the inlet flow pulsation condition and the static condition.
At the same time, in order to numerically verify the fluctuation behavior of the bubble under pulsation condition, the bubble equivalent diameter at each moment is calculated. The trend of the bubble equivalent radius over time under cases and is given by Figure 8. It can be seen from the curve change in Figure 8 that, during the whole process of bubble rising, the bubble equivalent radius almost fluctuates in the same cycle ( = 0.2s) as the flow pulsation condition. It appears that the periodic fluctuation of the equivalent radius also corresponds to the volatility of the bubble shape and volume in one cycle. In addition, the findings indicate that the pressure at the location of the bubble will decrease as it gets closer to the upper free surface during the rising process. Based on the ideal gas state equation, the bubble volume and equivalent radius tend to increase throughout the process. Moreover, since the bubble rises faster in case b, the overall trend of the bubble equivalent radius curve in case is higher than that in case a, as shown in Figure 8.
3.1.2. Single Bubble Rising Velocity
The bubble rising velocity is another important characteristic that should be concerned to reflect the bubble motion state. Under static condition, setting the same dimensionless number and other conditions as case b-d in Table 1, the change of the bubble centroid velocity with time is obtained, as shown in Figure 9. It appears that the single bubble rising velocity gradually increases to a certain stable value and then fluctuates in a small range for different dimensionless number of cases. However, when the inlet flow rate changes periodically, it will have a significant impact on bubble rising velocity. According to the simulation results of each case in Table 1, the influence of five independent variables on the bubble rising velocity under flow pulsation condition is found. In general, the bubble centroid velocity during the bubble rising process varies sinusoidally in the same cycle as the inlet flow rate under pulsation condition, as shown in Figure 10. This is another difference from the bubble rising process under static condition.
(a) The effect of Re on the velocity (case a-d)
(b) The effect of Bo on the velocity (case a, e, f)
(c) The effect of on the velocity (case a, , h)
(d) The effect of on the velocity (case k-o)
(e) The effect of τ on the velocity (case a, i, j)
The effects of each variable on bubble rising velocity under flow pulsation are discussed below. Based on the Reynolds number expression of (9), the Reynolds number increases as the viscosity of the corresponding liquid decreases, while other conditions are equal. When the viscosity of the corresponding liquid becomes smaller, the bubble rises faster in the liquid. Therefore, the results revealed that there is a positive correlation between the Reynolds number and the pulsation amplitude of the bubble velocity under flow pulsation condition, as shown in Figure 10(a). Another dimensionless Bond number characterizes the ratio of gravity to surface tension. As can be seen from Figure 10(b), the Bond number has slight effect on the bubble centroid velocity.
Figure 10(c) shows that the larger the initial flow velocity of the liquid phase is, the greater amplitude of periodic change in bubble centroid velocity becomes. Based on the result of Figure 10(d), it reveals that the fluctuation amplitude of the bubble centroid velocity becomes greater as the amplitude of the inlet flow pulsation increases. Figure 10(e) illustrates more information about the flow pulsation period; when the flow pulsation period is different, the bubble rising velocity fluctuates as almost the same pulsation period as the respective inlet flow rate. The velocity variation curves representing the respective pulsation periods are interlaced in the figures.
3.1.3. Single Bubble Aspect Ratio
In addition to calculating the bubble equivalent radius at each moment as aforementioned to numerically verify the fluctuating behavior of the bubble under flow pulsation condition, the variation of the bubble aspect ratio (the ratio of width to height of the bubble) with time under different variables is acquired, as shown in Figure 11. From all the results, the aspect ratio of the bubble also changes periodically under flow pulsation condition.
(a) The effect of Re on the aspect ratio (cases a-d)
(b) The effect of Bo on the aspect ratio (cases a, e, and f)
(c) The effect of on the aspect ratio (cases a, , and h)
(d) The effect of on the aspect ratio (cases k-o)
(e) The effect of τ on the aspect ratio (cases a, i, and j)
Figure 11(a) shows that the baseline of periodic fluctuation of the bubble aspect ratio increases with the Reynolds number when other conditions are the same. This is because the degree of lateral elongation of the bubble increases with the Reynolds number. Actually, this conclusion could also be derived from the study of Grace [27] and Hua [29]. When it comes to Bond numbers, it has slight effect on bubble aspect ratio, as in Figure 11(b). Figure 11(c) shows the influence of the initial flow rate of the liquid phase on the aspect ratio of the bubble. Specifically, the larger the initial flow rate of the liquid phase is, the greater periodic change effects on the bubble becomes. This will give rise to a more dramatic change to the bubble aspect ratio. The effect of the flow pulsation amplitude on the aspect ratio of the bubble is similar to the initial flow rate of the liquid phase, as shown in Figure 11(d). Similar to the aforementioned results of the flow pulsation period affecting the bubble rising velocity, Figure 11(e) shows that the bubble aspect ratio fluctuates as almost the same pulsation period as the respective inlet flow rate.
3.2. Bubble Pairs Coalescence under Flow Pulsation Condition
In this study, the coalescence process of bubble pairs under the condition of inlet flow pulsation is also simulated. The geometric and boundary conditions are shown in Figure 1(b). Other related parameter values are set as Re = 400, Bo = 0.4, = 0.03m/s, = 1, and τ = 5 ms. Firstly, the bubble coalescence process under flow pulsation condition is shown in Figure 12. Furthermore, in order to compare the differences of bubble coalescence process under flow pulsation condition and static condition, a bubble pair coalescence case under static condition is set, in which the geometric conditions and the Re/Bo numbers are the same as those of pulsation condition. This bubble coalescence process under static condition is depicted in Figure 12 with red-dashed line.
As shown in Figure 12, numerical simulations in both static and flow pulsations cases illustrate the appearance of gas bridges between the bubbles during coalescence. The large surface tension of the gas bridge drives the coalescence bubble to expand rapidly in the vertical direction and to shrink in the horizontal direction. Then the bubble will shrink longitudinally, and the lower surface will accelerate to the upper surface. However, the bubble change progress is not completely consistent under the two conditions. Under static condition, the bubble will gradually reach a steady state and continue to rise. As to pulsation condition, due to the existence of a strong pressure field caused by the inlet flow pulsation condition, the volume of both two bubbles shrinks during the initial rising phase as shown in Figure 12 (t = 2.50ms, 5.34ms). Then the two bubbles get closer and coalesce together. This makes the bubble coalescence slightly delayed. Similar to the static condition, the shape of the bubble will also undergo a series of changes due to the surface tension after coalescence. However, under the influence of the inlet pulsation condition, cyclical fluctuation in shape and volume will still occur during the final stage of bubble rising as shown in Figure 12 (t = 17.74ms ~ 29.30ms), which is another difference in the bubble coalescence processes under the two conditions.
Moreover, the bubble centroid velocity during bubble coalescence under two conditions is obtained, as shown in Figure 13. It can be seen that the two bubbles have a low velocity before coalescence under the static condition, but in the initial rising phase of the pulsation condition, the bubble centroid velocity fluctuates periodically. After the two bubbles coalesce and experience a short-term fluctuation, the bubble velocity is finally stable under static condition (after 0.023ms). However, the bubble centroid velocity still fluctuates under the flow pulsation condition. At the same time, the effect of the flow pulsation condition about the bubble coalescence time mentioned above can also be verified from Figure 13.
4. Conclusions
In the present study, the single bubble rising and coalescence of bubble pairs under flow pulsation condition are numerically investigated by using the MPS-MAFL method. The key conclusions of the present work are as follows.
Different from the single bubble rising process under static condition, the bubble shape and rising velocity will have periodic fluctuation under the flow pulsation condition, and the fluctuation period is almost the same as the flow pulsation period. Correspondingly, the shape and velocity of bubble will only gradually stabilize during the rising process under the static condition. For the coalescence process of bubble pairs, due to the introduction of flow pulsation condition, the bubble pairs begin to coalesce slightly later than the static condition, and after coalescence, the bubbles will continue to fluctuate with the same cycle as the inlet flow pulsation.
The calculation and analysis results of this study could help to reveal the bubble behavior under flow pulsation condition, and they establish the basis for further application of MPS method to the two-phase flow analysis.
Nomenclatures
| : | Bond number |
| : | Bubble diameter, m |
| : | Dimension number |
| : | Eotvos number |
| g: | Gravity acceleration, m/s2 |
| n: | Unit vector of phase interface |
| r: | Particle location vector |
| : | Particle distance, m |
| : | Reynolds number |
| : | Time, s |
| u: | Particle velocity vector |
| : | Weight function. |
| : | Diffusion coefficient |
| : | Local curvature |
| : | Viscosity, Pa·s |
| : | Density, kg/m3 |
| : | Surface tension coefficient, N/m. |
| : | Particle number. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The present study is supported by the National Natural Science Foundation of China (nos. 11875217, 11505177, and 11505134), the National Natural Science Foundation of China outstanding youth fund (11622541), and the Young Elite Scientists Sponsorship Program by CAST (2018QNRC001).