Numerical Study of Bubble Coalescence and Breakup in the Reactor Fuel Channel with a Vaned Grid

: The characteristics of bubbles of different sizes in fuel assembly are vital to two-phase ﬂow resistance and heat transfer capacity. However, due to the swirl ﬂow caused by the mixing vane, bubbles can crowd at the heated surface, which may anticipate the occurrence of departure from nucleation boiling. In the current work, the adiabatic two-phase ﬂow in a simpliﬁed fuel assembly was analyzed by using the Eulerian two-ﬂuid model and the MUSIG (MUltiple SIze Group) model. This computational domain consists of two coolant channels and two sets of vaned spacers, with three sets of periodic boundary conditions at the side faces of the domain. The distributions of vapor phase and bubble diameters were obtained, based on which the effects of mixing vanes on the bubble characteristics were analyzed. Vapor phase crowded at the rod surface in the higher inlet vapor fraction case, but crowded in the channel center in the lower inlet vapor fraction cases. This work can be used as a reference for the design of mixing vanes to avoid the anticipation of departure of nucleation boiling that may be caused by unreasonable design.


Introduction
Two-phase boiling flow has been widely used in power systems because of its high heat transfer performance compared with single-phase convective flow. In a newly designed pressurized water reactor (PWR), such as AP1000 (Westinghouse Co., Pittsburg, PA, USA), subcooled boiling and even saturated boiling is allowed in the core due to increasing confidence in the design and analysis technologies of PWR and the increasing demand on the economy of PWR. Under highly pressurized conditions, slight boiling in the reactor can significantly increase the heat transfer capability with no deterioration on the flow stability performance. However, the enhancement of the heat transfer ability is limited by the boiling crisis, where the wall heat flux reaches the critical heat flux and the local heat transfer coefficient decreases by one or two orders of magnitude. The wall temperature will increase dramatically by a hundred degrees under the condition where wall heat flux is specified, i.e., the cladding outer surface. Thus, the cladding and even the fuel pin may be burnt or melt, releasing vast radioactivity into the coolant.
Spacer grids and mixing vanes were placed in the fuel assembly to support the fuel rods and enhance the heat transfer capacity by generating the vortexes. However, as pointed by Krepper et al. [1] and Zhang et al. [2], bubbles will aggregate at parts of the rod surface due to the swirling caused by the mixing vane and that will lead to a vapor fraction peak at the rod surface, which may prevent the heat from being transferred to the liquid phase and result in the anticipation of the burnout. Thus, it is essential to investigate the vapor phase distributions for the two-phase flow in fuel assembly.
Two-phase flow has drawn lots of attention in the past few decades due to its essential role in the safety analysis of reactor systems. It has been studied by experimental and numerical methods

Eulerian Two-Fluid Model
There is no mass and heat transfer between the liquid and vapor phases for the adiabatic flow. The mass and momentum equations for liquid phase are given as: (α l ρ l ) + ∇ · (α l ρ l U l ) = 0 (1) ∂ ∂t (α l ρ l U l ) + ∇ · (α l ρ l U l U l ) = −α l p + ∇ · τ l + α l ρ l g + F lg (2) where α, ρ, U and τ denote the volume fraction, density, velocity and stress tensor for a specific phase, respectively; the subscript l means the liquid phase; g is the gravity; p is the pressure for both liquid and vapor phases, since these two phases share an equal pressure; F lg is the interphase forces caused by the non-equilibrium between two phases, including the drag force, lift force, wall lubrication force, turbulent dispersion force and virtual mass force, i.e., where F D , F L , F WL , F TD and F VM are the drag force, lift force, wall lubrication force, turbulent dispersion force and virtual mass force, respectively, which are defined in Section 2.3. All these forces are calculated based on the local parameters in each computational cell, such as velocity, density and vapor fraction of each phase and the bubble diameter. These phasic parameters, except the bubble diameter, can be obtained by solving the two-phase governing equations. The bubble diameter is one of the key parameters in estimation of F lg since all the interphase forces are calculated based on bubble diameter. However, the bubble diameter is not easy to determine due to the polydisperse characteristics of bubbles. The bubble size in a two-phase flow system may vary by one or two orders of magnitude, especially in a complicated geometry such as fuel assembly with mixing vanes. Many works were proposed to simplify the specification of bubble diameters in the two-phase CFD simulation. Thomas et al. [8] set the bubble diameters as a constant in the whole computation domain and tuned this value for different test conditions based on experience. Zhang et al. [3] and Krepper et al. [1] defined the bubble diameter as a function of localized fluid subcooling in each computational cell, which was better than Thomas et al.'s [8] work, but it still cannot present the polydisperse characteristics of bubble size. In the recent work on bubbly flow simulation, the bubble size distribution has been considered by using the PBM or the MUSIG models [9,10]. In the PBM and MUSIG models, bubble size fraction equations are solved along with a set of momentum equation for all of the gas phase. The bubble mean diameter, such as the Sauter mean diameter, is calculated based on the bubble size fraction and then used to solve the interphase transfer terms. These two models can deal with poly-dispersed bubbly flow with a limited range of bubble diameters. However, the direction of the interphase forces may change when bubble diameter exceeds a certain value, such as the lift force estimated by the Tomiyama model [11]. The method that averages bubble sizes to a single value cannot present different directions of interphase forces on bubbles with different sizes. The inhomogeneous MUSIG (iMUSIG) model was proposed to solve these problems with wide range of bubble diameters. Two sets of momentum equations were solved for small-diameter bubbles and large-diameter bubbles, respectively. The bubble diameters used in each momentum equation for the interphase actions were averaged from a set of bubble size fractions respectively. The diameter where lift force changes its sign is used as the critical value to divide the sets of bubble diameters. Thus, the governing equations for the poly-dispersed bubbles are given as follows: where the subscript i means the bubble group. S i and M i are the mass source and momentum source caused by the bubble breakup and coalescence. Due to the momentum conservation between liquid and vapor phases, the net interphase forces should be zero, i.e., As mentioned above, the iMUSIG model was proposed to deal with the two-phase flow with bubble size varying over a very large range exceeding the critical bubble size where lift force changes its direction. The homogeneous MUSIG model is qualified to model the bubble size distribution and its effects on the momentum transfer between two phases if most of the bubble sizes are lower than this critical size. It was found that size of most bubbles in the bubbly flow in fuel assembly with vaned spacers is lower than the critical value estimated by the Tomiyama model [11] after a series of simulation on the two-phase flow in fuel channels. Thus, in this work, the homogeneous MUSIG model was employed, and Equations (4)- (6) can be deduced into the following forms, i.e.,

MUSIG Model
The MUSIG model compiled in ANSYS CFX code is a framework based on the population balance model that can be used to handle the polydisperse bubbly flow with large size variation. Bubbles with different sizes are divided into several groups. The continuous bubble sizes are divided into several bins, each with a characteristic bubble size, i.e., d g,i . The size fraction function (f i ) is used to present the fraction of bubbles located in each bin. Effects of bubble breakup and coalescence can be modeled by the fraction transfer among different bins. Thus, the governing equations for the size fraction function can be given as following equation: where S f,i is the source term of size fraction function for each group of bubbles, which is caused by the bubble breakup and coalescence and modeled by following equation: where B B,i and B C,i are the birth rates due to larger bubble breakup and smaller bubble coalescence into group i; B B,i and D C,i are the death rates caused by the breakup into smaller bubbles and form into larger bubbles due to coalescence with other bubbles. These birth and death terms can be expressed as: where m i is the mass fraction of bubble size i. g(m i :m j ) denotes the frequency of bubbles of size i breakup into bubbles of size j and the left part of the bubble. Q(m i :m j ) represents the frequency of bubble coalescence between two bubbles from group i and j. X jki is the coalescence mass matrix, defined by: The bubble breakup frequency g(m i :m j ) is given by the Luo and Svendsen model [12], i.e., where ε l and ν l are the liquid phase eddy dissipation rate and kinematic viscosity, respectively. σ is the surface tension. The bubble coalescence frequency Q(m i :m j ) is modeled by the Prince and Blanch model [13]; that is, where

Phasic Interaction Models
Momentum transfer can happen between the two phases due to the interphase non-equilibrium caused by the velocity difference and the gradients of velocity and vapor fraction. In this work, the drag force, lift force, wall lubrication force, turbulent dispersion force and virtual mass force are taken into consideration to model the adiabatic two-phase flow in the fuel assembly. The models are detailed in this section.

Drag Force
The drag force presents the resistance introduced by the relative velocity between bubbles and the surrounding liquid, acting in the opposite direction of the fluid flow and having a significant effect on the two-phase velocity fields. It is calculated by [7]: where C D is the drag force coefficient, estimated by the Ishii-Zuber correlation [14], considering of the bubbly flow with different bubbly sharps, i.e., where C D (sphere), C D (ellipse) and C D (cap) are the drag force coefficients for sphere bubble, ellipse bubble and cap bubble, respectively, defined by: where the bubble Reynolds number (Re g ) and Eötvös Number (E o ) are defined by: and where σ is the surface tension. dg means the Sauter mean diameter of bubble estimated by the bubble sizes (d i ) and bubble number density (n i ), which can be calculated from bubble size fraction. The Sauter mean diameter is calculated by:

Lift Force
The lift force represents the interphase actions on the bubbles caused by the liquid shear velocity field. This force is always orthogonal to the flow direction; that is, pointing to the wall or center of the channel. That means the lift force has a significant impact on the radial distribution of the vapor phase volume fraction. Lift force is calculated by [7]: where C L is the lift force coefficient and calculated by the Tomiyama correlation modified by Frank et al. [15].
where E o is a modified Eötvös number, defined by:

Wall Lubrication Force
In the bubbly flow, small bubbles should be pushed to the wall since the lift force on small bubbles pointing towards the wall from fluid. However, it is observed that the bubbles are concentrated in the near wall region, but not immediately adjacent to the wall. It is because of the force acting on the dispersed bubble phase that push bubble away from the wall, named wall lubrication force, which is calculated by [7]: where n W is the unit normal vector pointing away from the wall and C WL is the wall lubrication force coefficient calculated by the Antal correlation [16]; that is, where y W is the distance to the nearest wall from the cell center. The coefficients C W1 and C W2 are set as −0.01 and 0.05 respectively. It should be noted that the lubrication force only acts in a very thin layer next to the wall, with a cut-off distance of five time of Sauter diameter, beyond which the lubrication force will be zero.

Turbulent Dispersion Force
The distribution of bubbles dispersed in the liquid phase is affected by the turbulent eddies of liquid phase. The turbulent dispersion force is such a force to account the impacts of turbulent fluctuations of liquid velocity on the gas bubbles. It has a significant role in the radial vapor fraction and can be estimated by [7]: where C TD and σ lg are coefficients, defined as 1.0 and 0.9, respectively [17]. µ t,l is the turbulent viscosity of liquid phase; C l,g is the momentum transfer coefficient caused by the interphase drag force, i.e.,

Virtual Mass Force
The virtual mass force will act on the bubble phase when bubble is accelerating relative to the continuous primary phase. This force is defined by [7]:

Model Validation
The experiments carried out by Vial et al. [18] were used as the benchmark to validate the models for two-phase flow. In their work, Laser Doppler technology was employed to measure the velocity and gas distribution in a bubble column. Gas was injected into a tube with 2 m length and 0.1 m diameter from a single orifice nozzle with 5 mm diameter at atmosphere pressure with gas injection velocity of 3.7 m/s. The profiles for velocity and gas volume fraction were measured in the downstream of the nozzle along the radial direction with uncertainties of ±15% and ±0.01 m/s, respectively. These data were used to validate current models. Only one-fourth of the fluid domain inside the tube was modeled given the symmetry of the geometry. Twenty groups of bubbles were set in the MUSIG model with initial bubble diameter of about 5 mm. The calculated vapor fraction and axial profiles were compared with experimental data in Figure 1. As can be seen, the predicted axial velocity agrees well with the experiment, generally. The gradient for predicted profile of axial velocity is slightly lower than that for experiment velocity profile, which might be because of the interfacial non-drag force, which is slightly related to the bubble shapes. Besides, it seems there is a slight increase for the axial velocity at the location of tube inner wall. Actually, it is a velocity decrease rather than a velocity increase. Reversed flow in the near wall region happens due to the concentrated vapor distribution in the tube center, which accelerates the two-phase mixture in the center of the channel and causes reverse flow in the near wall region. Thus, the axial velocity in the near wall region is negative. However, on the tube inner surface, the axial velocity is zero due to the non-slip boundary condition. Thus, it results in an increase in the profile of axial velocity near the wall. Besides, as can be noted, the predicted profile for vapor fraction locates in the uncertainties range of experiment data; while the predicted velocity profile doesn't. As for the two-phase modelling, to predict the velocity profile is always more challenge than to predict the vapor fraction profile. It is because there are still some flaws in the models to mimic the interactions between the two phases, such the turbulent dispersion force and the lift force, which effects the velocity profile in the radial direction. in the profile of axial velocity near the wall. Besides, as can be noted, the predicted profile for vapor fraction locates in the uncertainties range of experiment data; while the predicted velocity profile doesn't. As for the two-phase modelling, to predict the velocity profile is always more challenge than to predict the vapor fraction profile. It is because there are still some flaws in the models to mimic the interactions between the two phases, such the turbulent dispersion force and the lift force, which effects the velocity profile in the radial direction.

Geometry, Grid, Boundary Conditions and Numerical Scheme
In the typical reactor core, the reactor assembly comprises 17 by 17 fuel rods and about ten spacer grids. It requires too much computational cost to model the whole fuel assembly with 17 by 17 fuel rods, especially for two-phase in fuel assembly. In the past, researchers chose 5 by 5 fuel rods as a simplified geometry to study the single-phase and two-phase thermohydraulics by experiments and simulations. However, it still took up too many computational resources for the two-phase simulation in the assembly with 5 by 5 fuel rods with complicated physical models, such as the two-phase flow with MUSIG model. Fortunately, the arrangement of spacer grid and mixing vane shows some regularity and periodicity. In the work of Li and Gao [19], a two-channel model with three sets of periodic boundary conditions was proposed and compared with the results from a 5 by 5 rod geometry to prove the rationality of this simplification. This two-channel model was also employed in the current work for the purposes of reducing computational cost and obtaining reasonable results. The computational domain is shown in Figure 2. It consists of two channels, two sets of spacer grids and mixing vane. The cross-section connected to the adjacent channels is set to be periodic boundaries to allow the transverse flow to cross these boundaries. Besides, the weld bead as well as components to hold the bundles was omitted to simplify the modeling, since this work is focused on the effects of spacer grids and mixing vanes on thermohydraulics.
The computational domain was meshed by a mixing grid strategy, i.e., being divided into two kinds of subdomains, including the regions containing upstream and downstream of mixing vanes and the regions of mixing vanes. The upstream and downstream of mixing vanes were meshed by the swept hexahedral grids, while the subdomains of mixing vane vicinity were meshed by tetrahedral grids. These five subdomains are connected by interfaces with conformal grids to ensure the conservation at the interfaces.
As noted by Mahaffy et al.
[20], mesh independence can be checked in two ways, i.e., increasing the grid quantity or increasing the solver precision. In the current research, two steps were taken to check the mesh independence. First, three sets of mesh strategies containing 0.92 M, 1.66 M, and 2.25 M cell elements were employed for simulation with single precision solver. The total pressure drop and outlet size fractions obtained from these three sets of meshes were compared and the results form mesh with 1.66 M elements shown almost no difference (less than 0.1%) with results from mesh with 2.25 M cells. After that, a double precision solver was employed to compare with single precision results based on the mesh with 1.66 M elements. No difference has been found between these two sets of results. Thus, the mesh strategy with 1.66 M cell elements was chosen as the independent one and used for further analysis.
Energies 2018, 11, 256 9 of 16 d outlet size fractions obtained from these three sets of meshes were compared and the results form sh with 1.66 M elements shown almost no difference (less than 0.1%) with results from mesh with 5 M cells. After that, a double precision solver was employed to compare with single precision ults based on the mesh with 1.66 M elements. No difference has been found between these two s of results. Thus, the mesh strategy with 1.66 M cell elements was chosen as the independent one d used for further analysis. Liquid-vapor two-phase mixture with a given vapor volume fraction flows into the in boundary with a uniform velocity profile of 5.0 m/s. Twenty sets of bubble sizes varying from 0.0 mm to 5 mm were employed to capture the size distribution of vapor phase. The bubble sizes a discretized by the geometric mass discretization method, i.e., where dmin and dmax are minimum and maximum bubble diameter limitations; N is the total numb of bubble group. All the governing equations are discretized by the high-resolution format and solved by CF 17.1 code. Besides, the volume fraction was solved coupled with the velocity-pressure equations f a better convergence performance. The convergence criteria were set to be

Cases Setup
In the two-phase flow with dense bubbles, two bubbles may coalesce if they collide with ea other in the fluid domain under the actions of turbulence, eddy capture, velocity gradient or bo force. Besides, large bubbles may break up due to forces from the surrounding liquid, such as t shear force caused by the liquid velocity gradient, the turbulence fluctuation and the interfacial sl All these mechanisms on the breakup and coalescence are related to the density and the size of bubb The baseline case was set with the inlet vapor fraction of 0.3 and all the bubbles in group 2; that the diameter of all the bubbles is 0.071 mm. Besides, three more cases were carried out as the contr group to analyze the effects of inlet VOF and the inlet bubble diameter on the three-dimension thermohydraulics parameters (Table 1). The reference inlet VOF is 0.1 and the reference bubb diameter is 0.451 mm (all bubbles in group 10).  Liquid-vapor two-phase mixture with a given vapor volume fraction flows into the inlet boundary with a uniform velocity profile of 5.0 m/s. Twenty sets of bubble sizes varying from 0.005 mm to 5 mm were employed to capture the size distribution of vapor phase. The bubble sizes are discretized by the geometric mass discretization method, i.e., where d min and d max are minimum and maximum bubble diameter limitations; N is the total number of bubble group. All the governing equations are discretized by the high-resolution format and solved by CFX 17.1 code. Besides, the volume fraction was solved coupled with the velocity-pressure equations for a better convergence performance. The convergence criteria were set to be 1.0 × 10 −4 for all the variables.

Cases Setup
In the two-phase flow with dense bubbles, two bubbles may coalesce if they collide with each other in the fluid domain under the actions of turbulence, eddy capture, velocity gradient or body force. Besides, large bubbles may break up due to forces from the surrounding liquid, such as the shear force caused by the liquid velocity gradient, the turbulence fluctuation and the interfacial slip. All these mechanisms on the breakup and coalescence are related to the density and the size of bubble. The baseline case was set with the inlet vapor fraction of 0.3 and all the bubbles in group 2; that is, the diameter of all the bubbles is 0.071 mm. Besides, three more cases were carried out as the control group to analyze the effects of inlet VOF and the inlet bubble diameter on the three-dimensional thermohydraulics parameters (Table 1). The reference inlet VOF is 0.1 and the reference bubble diameter is 0.451 mm (all bubbles in group 10).

Results and Discussions
The results can be analyzed after reaching the convergence results for each case. All the cases were carried out with a transient solver strategy for better convergence performance. It will take about 50 h with 48 cores on a HPC for each case. First of all, the pressure drop characteristics along the flow channel are given in Figure 3. As can be seen, the initial bubble diameter makes no difference on the flow resistance for two-phase flow, while the vapor fraction makes a significant difference. Besides, most of the pressure drop in the channel was presented from the down surface of the spacer to the down surface of the mixing vanes. That means the main part of the pressure drop was a consequence of the flow area reduction at the low surface of spacer, rather than of the mixing vane. The results can be analyzed after reaching the convergence results for each case. All the cases were carried out with a transient solver strategy for better convergence performance. It will take about 50 hours with 48 cores on a HPC for each case. First of all, the pressure drop characteristics along the flow channel are given in Figure 3. As can be seen, the initial bubble diameter makes no difference on the flow resistance for two-phase flow, while the vapor fraction makes a significant difference. Besides, most of the pressure drop in the channel was presented from the down surface of the spacer to the down surface of the mixing vanes. That means the main part of the pressure drop was a consequence of the flow area reduction at the low surface of spacer, rather than of the mixing vane. The distribution of vapor fraction, which is crucial for heat transfer capacity at the rod surface, is shown in Figures 4-6. Figures 4 and 5 present the vapor fraction profiles at the planes downstream from the spacer with different distances for cases 1 and 4, respectively. As can be seen, there are significant differences between the results from these two cases. As for the first case with inlet vapor fraction of 0.3 (high vapor fraction case), there was obvious vapor crowding in the wake region of the transverse flow around the rod surface, especially at the plane 0.08 m from the spacer. Besides, a large amount of bubbles were aggregated at the center of the channel with vapor fraction lower than the vapor fraction near the rod surface. That means distinct bubble aggregation can be observed in the high vapor fraction case, which may lead to the departure from nucleation boiling (DNB) on the rod surface. When it comes to the low vapor fraction case, i.e., case 4 with inlet vapor fraction of 0.1, most of the bubbles aggregated at the center of channel, instead of the rod surface. Thus, there is no risk for the DNB in the low vapor fraction case. The vapor fraction distributions on the rod surface are also illustrated in Figure 6. In case 1, the vapor fraction at rod surface in the upstream of mixing vane was about 0.5. However, the vapor fractions on more of the rod surface in the downstream were much larger than 0.5, and can be as high as 0.9. In this case, the mixing vane would increase the vapor fraction at the wall, rather than decrease it, which is one of the purposes of the mixing vanes. On the other hand, the vapor fraction in the downstream region can be significantly reduced by the mixing vanes in the case of low inlet vapor fraction.
In addition to the vapor fraction, bubble diameters are also critical in the theory of bubble crowding for DNB. Figure 7 presented the diameter of the bubble at the plane of 0.08 m, which is the plane with maximum vapor fraction on the wall. As can be seen, for case 1, the bubble diameters near the rod surface and in the channel centers were about 0.35 mm. However, the bubble diameters in the channel centers for the low vapor fraction case were as high as 0.45 mm, which was because most of the bubbles were entrained into the center of the channels and coalesced into larger bubbles. Besides, the Sauter averaged diameter along the flow direction for each case was analyzed and is presented in Figure 8. It can be observed from this figure that the vapor fraction had obvious effects on the bubble diameter for the whole flow path, while inlet bubble diameter only made a big difference in the upstream of the first vaned spacer. After passing the first mixing vane, the Sauter diameters were almost the same for the cases with same inlet vapor fraction but different bubble size. That means the bubble breakup caused The distribution of vapor fraction, which is crucial for heat transfer capacity at the rod surface, is shown in Figures 4-6. Figures 4 and 5 present the vapor fraction profiles at the planes downstream from the spacer with different distances for cases 1 and 4, respectively. As can be seen, there are significant differences between the results from these two cases. As for the first case with inlet vapor fraction of 0.3 (high vapor fraction case), there was obvious vapor crowding in the wake region of the transverse flow around the rod surface, especially at the plane 0.08 m from the spacer. Besides, a large amount of bubbles were aggregated at the center of the channel with vapor fraction lower than the vapor fraction near the rod surface. That means distinct bubble aggregation can be observed in the high vapor fraction case, which may lead to the departure from nucleation boiling (DNB) on the rod surface. When it comes to the low vapor fraction case, i.e., case 4 with inlet vapor fraction of 0.1, most of the bubbles aggregated at the center of channel, instead of the rod surface. Thus, there is no risk for the DNB in the low vapor fraction case. The vapor fraction distributions on the rod surface are also illustrated in Figure 6. In case 1, the vapor fraction at rod surface in the upstream of mixing vane was about 0.5. However, the vapor fractions on more of the rod surface in the downstream were much larger than 0.5, and can be as high as 0.9. In this case, the mixing vane would increase the vapor fraction at the wall, rather than decrease it, which is one of the purposes of the mixing vanes. On the other hand, the vapor fraction in the downstream region can be significantly reduced by the mixing vanes in the case of low inlet vapor fraction.
In addition to the vapor fraction, bubble diameters are also critical in the theory of bubble crowding for DNB. Figure 7 presented the diameter of the bubble at the plane of 0.08 m, which is the plane with maximum vapor fraction on the wall. As can be seen, for case 1, the bubble diameters near the rod surface and in the channel centers were about 0.35 mm. However, the bubble diameters in the channel centers for the low vapor fraction case were as high as 0.45 mm, which was because most of the bubbles were entrained into the center of the channels and coalesced into larger bubbles. Besides, the Sauter averaged diameter along the flow direction for each case was analyzed and is presented  Figure 8. It can be observed from this figure that the vapor fraction had obvious effects on the bubble diameter for the whole flow path, while inlet bubble diameter only made a big difference in the upstream of the first vaned spacer. After passing the first mixing vane, the Sauter diameters were almost the same for the cases with same inlet vapor fraction but different bubble size. That means the bubble breakup caused by the mixing vane will erase the effects of inlet bubble size. Furthermore, the strategy to investigate the bubble diameter characteristics by using only two sets of spacer grids was reasonable, since the parameters show periodic tendency after flowing pass the mixing vanes.
Energies 2018, 11,256 12 of 16 by the mixing vane will erase the effects of inlet bubble size. Furthermore, the strategy to investigate the bubble diameter characteristics by using only two sets of spacer grids was reasonable, since the parameters show periodic tendency after flowing pass the mixing vanes.   by the mixing vane will erase the effects of inlet bubble size. Furthermore, the strategy to investigate the bubble diameter characteristics by using only two sets of spacer grids was reasonable, since the parameters show periodic tendency after flowing pass the mixing vanes.      As mentioned in Section 2.3.2, the lift force was estimated by the Tomiyama model. However, as can be noted from Equation (29), the directions of lift forces for large and small bubble were not the same. That is, the bubbles can be pushed by the lift force towards the bulk or the wall, depending on the bubble diameter. Thus, there should be a critical value to divide these two sets of bubbles. For the current case with pressure of 15.5 MPa, this critical diameter was 2.1 mm. The bubbles larger or smaller than this critical diameter should be modeled by two sets of momentum equations, since they had different transverse velocity. However, in the current work, we use the homogenous MUSIG model with just one set of momentum equation, rather than the in-homogenous MUSIG model. This is because, according to our pre-analysis, most of the bubbles were smaller than 2.1 mm. Thus, we assumed that all the bubbles were in the smaller group of bubbles. Here, to validate the assumption, we calculated the fraction of each set of bubbles in Figure 9. In the upstream of the first spacer, about one-third of the vapor fraction located in the large group. However, the volume fractions of the large bubbles in the downstream of the mixing vane were rather small when compared with the fractions of small bubbles. Given the purpose of the current work (to investigate bubble information in the fuel assembly under the effects of vaned spacers), the MUSIG model was qualified for the current simulation. In the downstream of the mixing vane, the bubbles break up due to the strong shear and turbulence caused by the mixing effects of the vanes. Then, the small bubbles may coalesce due to the collision of two adjacent bubbles under the actions of turbulence, eddy capture, velocity gradient, body force and wake. The process of bubble breakup and coalescence was a dynamic process, and even in that region a balance between breakup and coalescence was obtained. The dynamic evolution for each group of bubbles was investigated by analyzing the size fraction of bubbles. Figure 10   As mentioned in Section 2.3.2, the lift force was estimated by the Tomiyama model. However, as can be noted from Equation (29), the directions of lift forces for large and small bubble were not the same. That is, the bubbles can be pushed by the lift force towards the bulk or the wall, depending on the bubble diameter. Thus, there should be a critical value to divide these two sets of bubbles. For the current case with pressure of 15.5 MPa, this critical diameter was 2.1 mm. The bubbles larger or smaller than this critical diameter should be modeled by two sets of momentum equations, since they had different transverse velocity. However, in the current work, we use the homogenous MUSIG model with just one set of momentum equation, rather than the in-homogenous MUSIG model. This is because, according to our pre-analysis, most of the bubbles were smaller than 2.1 mm. Thus, we assumed that all the bubbles were in the smaller group of bubbles. Here, to validate the assumption, we calculated the fraction of each set of bubbles in Figure 9. In the upstream of the first spacer, about one-third of the vapor fraction located in the large group. However, the volume fractions of the large bubbles in the downstream of the mixing vane were rather small when compared with the fractions of small bubbles. Given the purpose of the current work (to investigate bubble information in the fuel assembly under the effects of vaned spacers), the MUSIG model was qualified for the current simulation. As mentioned in Section 2.3.2, the lift force was estimated by the Tomiyama model. However, as can be noted from Equation (29), the directions of lift forces for large and small bubble were not the same. That is, the bubbles can be pushed by the lift force towards the bulk or the wall, depending on the bubble diameter. Thus, there should be a critical value to divide these two sets of bubbles. For the current case with pressure of 15.5 MPa, this critical diameter was 2.1 mm. The bubbles larger or smaller than this critical diameter should be modeled by two sets of momentum equations, since they had different transverse velocity. However, in the current work, we use the homogenous MUSIG model with just one set of momentum equation, rather than the in-homogenous MUSIG model. This is because, according to our pre-analysis, most of the bubbles were smaller than 2.1 mm. Thus, we assumed that all the bubbles were in the smaller group of bubbles. Here, to validate the assumption, we calculated the fraction of each set of bubbles in Figure 9. In the upstream of the first spacer, about one-third of the vapor fraction located in the large group. However, the volume fractions of the large bubbles in the downstream of the mixing vane were rather small when compared with the fractions of small bubbles. Given the purpose of the current work (to investigate bubble information in the fuel assembly under the effects of vaned spacers), the MUSIG model was qualified for the current simulation. In the downstream of the mixing vane, the bubbles break up due to the strong shear and turbulence caused by the mixing effects of the vanes. Then, the small bubbles may coalesce due to the collision of two adjacent bubbles under the actions of turbulence, eddy capture, velocity gradient, body force and wake. The process of bubble breakup and coalescence was a dynamic process, and even in that region a balance between breakup and coalescence was obtained. The dynamic evolution for each group of bubbles was investigated by analyzing the size fraction of bubbles. Figure 10 shows the size fractions of the bubbles in groups 2, 6, 8, 11 and 14 for cases 1 and 4, respectively. At the inlet of computational domain, all the bubbles were in group 2. Prior to the first vaned spacer, the fraction of group 2 decreased sharply to a low value of about zero. Meanwhile, a fraction of groups 6 and 8 In the downstream of the mixing vane, the bubbles break up due to the strong shear and turbulence caused by the mixing effects of the vanes. Then, the small bubbles may coalesce due to the collision of two adjacent bubbles under the actions of turbulence, eddy capture, velocity gradient, body force and wake. The process of bubble breakup and coalescence was a dynamic process, and even in that region a balance between breakup and coalescence was obtained. The dynamic evolution for each group of bubbles was investigated by analyzing the size fraction of bubbles. Figure 10 shows the size fractions of the bubbles in groups 2, 6, 8, 11 and 14 for cases 1 and 4, respectively. At the inlet of computational domain, all the bubbles were in group 2. Prior to the first vaned spacer, the fraction of group 2 decreased sharply to a low value of about zero. Meanwhile, a fraction of groups 6 and 8 increased right after flow into the computational domain. Then, group 11 became the dominated group with the largest size fraction for both of the two cases. In this period, the bubble coalescence was the dominating mechanism for the bubble evolution. After flowing past the vaned spacer, the fractions of larger bubbles (groups 11 and 14) dropped dramatically while fractions of moderately large bubbles (groups 6 and 8) increased immediately for case 1, due to the swirl flow caused by mixing vane. In this stage, the bubble breakup dominated. However, for case 4, only a fraction of group 6 increased sharply after the mixing vane, which means the bubbles after the breakup in the case of the larger inlet vapor fraction were larger than bubbles in the case of the lower inlet vapor fraction. In the region of about 0.2 m downstream of the mixing vane, bubbles coalesced again into larger bubbles in case 1, while this distance becomes about 0.3 m in case 4 with a lower inlet vapor fraction. The same phenomena can be observed in the region near the second spacer due to the periodicity of the geometry and the flow characteristics. group with the largest size fraction for both of the two cases. In this period, the bubble coalescence was the dominating mechanism for the bubble evolution. After flowing past the vaned spacer, the fractions of larger bubbles (groups 11 and 14) dropped dramatically while fractions of moderately large bubbles (groups 6 and 8) increased immediately for case 1, due to the swirl flow caused by mixing vane. In this stage, the bubble breakup dominated. However, for case 4, only a fraction of group 6 increased sharply after the mixing vane, which means the bubbles after the breakup in the case of the larger inlet vapor fraction were larger than bubbles in the case of the lower inlet vapor fraction. In the region of about 0.2 m downstream of the mixing vane, bubbles coalesced again into larger bubbles in case 1, while this distance becomes about 0.3 m in case 4 with a lower inlet vapor fraction. The same phenomena can be observed in the region near the second spacer due to the periodicity of the geometry and the flow characteristics.

Conclusions
In the current work, the distribution of vapor phase in a simplified fuel assembly, which is related to the anticipation of critical heat flux, were investigated by a two-phase CFD code CFX with consideration of the characteristics of bubble coalescence and breakup. To model the two-phase flow with bubble size evolution, the Eulerian two-fluid model was employed with the interphase action models for momentum and turbulence transfer and the MUSIG model for bubble size evaluation. After validation, the models were used to simulate the two-phase flow in the assembly channel with periodic boundary conditions and two sets of vaned spacers. The following conclusions can be drawn from the results analysis.
(1) Vaned spacers caused sharp pressure drop, which is associated with the vapor volume fraction, rather than the inlet bubble size. (2) Vapor phase crowded at the rod surface for the higher inlet vapor fraction case, but crowded in the channel center for the lower inlet vapor fraction cases, i.e., increasing the inlet fraction can increase the risk of bubble aggregation at the rod surface, which might anticipate the critical heat flux under the condition of boiling two-phase flow. (3) Mixing vanes can reduce the average bubble size by breakup. However, the bubble size distribution downstream of the mixing vane was related to the inlet vapor fraction, rather than the inlet bubble size.