Numerical Investigation of Helium Bubble Rising Behavior in Cross-Type Channel

In a certain nuclear reactor, the coolant channel is a cross-type channel. There is little research about bubble rising behavior in this kind of channel, but the neutron detection would be affected by the bubbles. Through comparing with previous experiment, proper numerical model has been determined, and the cross-type channel under different helium flow rates is studied with CFD software Fluent. Research shows that the bubble distribution exhibits some similarities under different helium flow rates, and the bubbles tend to rise near the tube wall; In the axial direction of the flow channel, the bubble coalescence can be observed, which leads the void fraction in the axial direction to increase first and then to be stable; The void fraction increases with increasing helium flow rates.


INTRODUCTION
The cross-section of the coolant channel is cross-shape in a certain nuclear reactor, and in some cases there are hydrogen bubbles in it. Under the effect of various forces such as buoyancy, resistance, pressure, etc., the bubbles are unstable (the trajectory of the rising bubbles is not a line, but more like the "Z," and also the shape of bubbles changes a lot) during the rising process, which leads to significant deformation and coalescence. The existence of these bubbles in the channel will cause errors in neutron detection in the nuclear reactor, furthermore it will affect the safety of the reactor. so it is of great necessity to do some research about the bubble rising behavior in cross-type channel, which can help other researchers know where to put the neutron detector is better (less influence by bubbles).
There are two types of research about submerged bubbles rising: submerged jet upward and submerged jet downward as shown in Figure 1. These researches are mainly carried out in a large water tank or conventional pipe, through submerging a small diameter pipe with gas to generate bubbles to study the laws of bubble formation, detachment, rising behavior, etc. In the submerged jet upward, the effects of surface tension, nozzle diameter, gas flow rates, liquid density, liquid viscosity, rolling conditions, etc. have been extensively studied (Gerlach et al., 2006;Buwa et al., 2007;Guang, 2010;Haijing, 2010;Chen et al., 2019;Khan et al., 2020). There exists bubble coalescence during the rising process of multiple bubbles. Related researches have been performed by Haijing (2010) and Ying (2014); In the submerged jet downward, the influencing factors of bubble formation and detachment are also related to the immersion depth of the nozzle, nozzle shape, the inner and outer diameter, etc. (Tsuge et al., 2005;Wei, 2016;Xuan and Jingjing, 2016;Xuan and Songyang, 2019). In these researches, the focus is on the study of bubbles formation and detachment, but less attention is paid to the distribution of bubbles during the rising process. On the other hand, considering the effect of wall on the rising bubbles, the study about rising bubble in a narrow rectangular channel has been carried out (Jin et al., 2014;Hashida et al., 2019). However, these studies about narrow rectangular channel also can't reflect the situation in the cross-type flow channel. The width of the cross-type flow channel is wider than a narrow rectangular one, and the cross-type flow channel is not a conventional channel which causes bubble motion more complex. On the other hand, the previous researches did not pay much attention to the void fraction in the axial direction of the flow channel, but void fraction can exactly reflect the influence of the bubbles on neutron detection in the reactor. In summary, existing researches, neither bubble rising behavior in large-sized flow channels nor narrow rectangular channels can well reflect the situation in cross-type flow channel.
The most advantage of numerical simulation is cheap, and it can get more detailed parameters than experiments. So the purpose of the present study is to investigate the rising behavior of helium bubbles (considering that the helium is often used instead of hydrogen for research in the laboratory) in cross-type channel through CFD numerical simulations, which can provide a basis for the improvement of neutron detector in the nuclear reactor.

Numerical Model
The VOF model (Jiapeng, 2014) is built on the premise that two or more fluids (or phases) do not mix with each other, which is used to track the moving interface. The basic principle is to use the volume fraction of each phase in each grid cell and the corresponding function to determine the interface of each phase, so as to indirectly determine the changes of the fluids, rather than track the movement of the particles on the   interface. When dealing with the volume fraction of different phases passing through the region, the VOF model assumes that there is a clear interface between the different phases, and there is no interpenetration between each other. This makes the volume fraction of each phase independent of each other, which means when a phase is added, a corresponding volume fraction equation is introduced, combining with each corresponding separate momentum equation, the simulation of the motion of two or more fluids can be realized. In each control unit, the sum of the volume fractions of all phases is 1. When calculating each unit, the properties and volume fractions of each phase are used to obtain a volume average value of the physical parameters in the unit. The obtained average value is shared by every phase in the unit, so that the properties in the grid unit may be a certain phase or a mixture of multiple phases, which depending on the volume fraction. In other words, in the unit, if the volume fraction of the nth phase fluid is α n , there are three kinds of situations: 1. α n = 0: There is no nth phase in the unit; 2. α n = 1: Filled with nth phase in the unit; 3. 0 < α n < 1: There exists interface between nth phase and another phase (or other phases).
The SST (shear stress transport) k-ω model (Ansys Fluent 17.2 Theory Guide, 2016) was developed from the Baseline k-ω model. The SST k-ω model combined the stability of standard k-ω model near the wall and the independence of the k-ε model outside the boundary layer. Also the propagation of turbulent shear stress affected by the turbulent viscosity was considered. The SST kω model has higher accuracy and credibility in the simulation of various complex flow conditions than the standard k-ω and k-ε model.

Validation
For researching the law of bubble formation in the process of gasliquid stirring during the industrial production, an experiment of submerged jet directed upward has been done by Wei (2016), as shown in Figure 2. And the data like bubble size under several flow rates has been achieved. In this experiment, the bubble that has just left the nozzle is selected as the research object, at this time, the bottom of the bubble just rises to be aligned with the bottom end of the nozzle. In the horizontal direction, the size of the bubble is obtained directly by the ruler; in the vertical direction, the distance between the top and bottom of the bubble is measured, which is the vertical size of the bubble. Figure 3 shows the selection criteria.
In order to validate the rationality of the adopted numerical method, a corresponding numerical simulation model according to Wei (2016)'s experiment has been established. Figure 4 shows the schematic diagram.
The similarity in Figure 5 preliminarily illustrates the rationality of the simulation method. Furthermore, the bubble size of experimental results and numerical results have been compared as shown in Figure 6. The maximum error between the numerical results and the experimental results does not exceed 12.2%, and the trend of the increase of the bubble size with the increase of the helium rates is also coincident. Considering that the numerical simulation itself is under the ideal conditions, so the error is acceptable, which verifies the correctness of this numerical method. Figure 7 shows the cross-type channel model, the axial length of the channel is 185.25 l, the outer diameter of the tube is 1.33 l (the inner diameter is l, l = 6.1 mm) in size, and the tube reaches 3.28 l from the bottom of the flow channel. Figure 7 also shows the radial size of flow channel.

Simulation of Cross-Type Channel
Considering that the cross-type channel is a symmetric structure, a quarter symmetry is adopted to improve the calculation efficiency. The calculation domain is shown in of Figure 8 (the region inside the red frame), and there are two symmetric boundaries, the velocity inlet is used at the bottom of the pipe. The top of the cross-type channel is the pressure outlet, and the rest are the walls. We used hexahedral grids, and the wall surface and the inlet part are densified. Because it is a transient problem, the PISO algorithm is adopted. For the convergence condition, the Fluent default values are used.
Aiming to the number of cells of 481181, 845306, and 1622105, the validation of mesh independence has been established. Four different inlet helium flow rates from 0.06 to 1.2 m/s were selected  as the comparison conditions under different number of cells. As shown in Figure 9, the bubble sizes (the bubbles measured at the exit of the nozzle as described above) are almost the same (order of 0.1 mm) with these under the number of cells of 845306 and 1622105, even the gap between 481181 grids and 1622105 grids is small. So in this research a meshing method between 481181 and 845306 grids is adopted.

RESULTS AND DISCUSSION
In this study, a numerical simulation of helium flow rates from 0.01 to 2 L/min was carried out. The study found that some characteristics of bubbles have obvious regularity under different flow rates of helium. Due to space limitations, the following will select some representative working conditions for detailed discussion.
For the sake of analysis, the channel has been divided into three parts, the inlet is z = 0 m, the region of z < 0.2 m is called bottom area, 0.2 m < z < 0.6 m is middle area, z > 0.6 m is top area as shown in Figure 10.

Research on Bubble Distribution Characteristics
Through the contours of phases and animation in the axial direction of the cross-type channel, we found that the bubble size increased with the increase of helium flow rates, manifesting the area occupied by the helium (red one) became larger and larger in the axial section, as shown in Figure 11; Combining the animation of the contours of phases in the axial and radial sections, it can be seen that most bubbles tended to rise near the tube wall (the smaller the flow rates are, the more obvious it will be). There were also a few smaller bubbles distributed far away from the tube wall (some can reach one of the four corners of the cross-type channel), and move downward during the process of bubble rising.
In the radial section, take the helium flow rate of 2 L/min for instance. As shown in Figure 12, at z = 0 m, the distribution of helium was concentrated at the nozzle. Initially, the shape of the bubble was relatively regular, and the helium distribution at the cross-section was approximately circular (quarter symmetry), but then the bubble gradually became irregular, as shown in Figure 12A t = 3.040 s, the helium distribution on the crosssection was biased to one side. The reason is that the rising bubbles detached from the nozzle cause the flow disturbance which makes the bubbles near nozzle under complex stress. Figure 13 shows the radial helium distribution at helium flow rate of 0.5 L/min. It is similar to other working conditions, just the proportion of helium in the cross-section of the flow channel decreasing with reduction of helium flow rates. Overall, at high flow rates (above 1 L/min), the movement of bubbles in the flow channel is similar to the slug flow in a vertical pipe. At low flow rates (below 0.1 L/min), the form is similar to the bubbly flow in a vertical pipe.

Research on the Characteristics of Bubble Coalescence
Through flow animation under different flow rates, it can be found that in the middle area and top area of the flow channel, large bubbles were distributed at a certain interval in the axial direction, and among the large bubbles were small bubbles (not obvious when helium flow rates is lower than 0.1 L/min) as shown in Figure 14.
During the process of bubble rising, the motion of fluid around the bubble will affect the subsequent bubbles. Taking the velocity vector distribution at the flow rate of 2 L/min for instance as shown in Figure 14, the direction of the fluid vector on the side of the rising bubble was opposite. This is why some small bubbles move downward within a certain range. Also, the direction of the fluid vector at the tail of the rising bubble was the same, which accelerated the subsequent bubbles to catch up with the previous one, then bubble coalescence happened. In the axial direction, the bubbles were getting larger and larger due to the bubble coalescence, which leads to an upward trend in void fraction at the bottom area of the flow channel along the axial direction. In the area after the middle area of flow channel, as the bubbles were accelerated by buoyancy, the distance among bubbles became large, which weakened the influence of surrounding fluid on bubbles. Finally, large bubbles were distributed at a certain interval, and some small bubbles were separated among the large bubbles. Figure 15 shows the velocity vector distribution of helium flow rates 0.5 L/min. Compared with the condition of 2 L/min, the rising bubbles had less disturbance on the flow field, but the bubble coalescence still exited, especially at the bottom area of the channel. It was just not as frequent as the helium flow rates of 2 L/min. Figure 16 shows the condition of helium flow rates 0.50 L/min. As the flow rates of helium gas was further decreased, the bubble volume became smaller and smaller, and the flow channel space occupied also became smaller, so that the disturbance caused by the rising bubble became smaller, as a result, the frequency of the bubble coalescence is reduced.

Research on the Influence of Flow Rate on the Void Fraction
Analyzing 10 planes in the axial direction of the flow channel, they are planes of z = 0, z = 0.02, z = 0.04, z = 0.06, z = 0.08, z = 0.1, z = 0.2, z = 0.4, z = 0.6, z = 0.8 m respectively, as shown in Figure 17.
The void fraction of a plane can be obtained through calculating the proportion of the helium share in the selected plane. The graph of the void fraction of each plane with time is shown as followed (The green dotted line is the average of the void fraction). The reason for the above phenomenon is the bubble coalescence, as analyzed in Section "Research on the Characteristics of Bubble Coalescence." Within a certain distance, the latter bubble is accelerated by the disturbance of flow field and will catch up with the previous bubble, which leads to bubble coalescence. As a result, causing an upward trend in void fraction of the bottom area of the flow channel in the axial direction, and in the middle area and top area of the flow channel, due to the accelerated rising bubbles by the buoyancy, the distance among bubbles in the axial direction becomes larger, resulting in the decrease of bubble coalescence frequency. Eventually, large bubbles are distributed at a certain interval, and some small bubbles are separated by large bubbles, this is why the maximum value of the void fraction of a certain cross-section is much larger than the average value.
When the flow rates of helium are 0.5 L/min, the variation law of the void fraction in the axial direction of the flow channel is similar to the situation of 2 L/min helium flow rates: the maximum transient void fraction of a plane of the flow channel is much larger than its average value (maximum transient void fraction can reach 0.14, and maximum average void fraction is only 0.019); in the bottom area of the flow channel, the average void fraction increases rapidly, while in the middle area and top area of the flow channel does not change much and tends to be stable. Figure 19 shows the curve of void fraction vs. time and the average of the void fraction. Figure 20 shows the curve of void fraction vs. time and 0.05 L/min helium flow rates. It can be found that the void fraction of each plane is generally very low under 0.05 L/min helium flow rates. The maximum transient void fraction is only about 0.042, and the maximum average void fraction is even lower than 0.005. So, it is not rather meaningful to further analyze the void fraction.

SUMMARY AND CONCLUSION
In this study, taking the influence of bubbles on neutron detection as the background, a numerical simulation for the rising rules of helium bubbles in the cross-type flow channel has been carried out with CFD software. The major conclusions are as the following: 1. Because of flow disturbance caused by rising bubbles, the bubble coalescence exists widely in the middle area and bottom area of the flow channel (z < 0.6 m). The larger the helium flow rates are, the more frequent the coalescence will be. In the top area of the flow channel, due to the larger distance among the bubbles, the disturbance effect is small, so the coalescence is hard to happen. In the bottom area of the flow channel, the coalescence is most frequent, because of short distance; 2. During the process of bubble rising, the larger bubbles tend to move near to the tube wall, and some smaller bubbles tend to be far away from the tube wall (some can reach the one of the four corners of the cross-type channel), moving downward, due to the disturbance of the larger bubbles on the flow field; 3. When the helium flow rate is greater than 0.1 L/min, the average value of the void fraction increases rapidly along the axial direction, and then tends to stabilize in the middle area and top area of the flow channel (z > 0.2 m), and the greater the helium flow rate, the greater the void fraction.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.

AUTHOR CONTRIBUTIONS
ZM as the tutor guided the whole process of the simulation. RY and DM performed the selection and validation of the numerical method. HJ carried out the validation of mesh independence. CZ planned and carried out the simulation of cross-type channel. XZ finished the data processing. CZ and XZ wrote the manuscript with input from all authors. All authors provided critical feedback and helped shape the research, analysis, and manuscript.