The Eulerian two-fluid model is employed to simulate the liquid and vapor phases within the studied rod bundle channel. The phases are treated as interpenetrating continua coexisting everywhere in the flow domain. Governing equations for mass, momentum, and energy for each phase are solved and formulations describing phase interactions at interfaces are introduced to close the equations.

Continuity equation: ∂ ∂ t α k ρ k + ∇ ⋅ α k ρ k u → k = m ˙ i k − m ˙ k i where α k , ρ k and u → k are the volume fraction, density and velocity of phase k, respectively. m ˙ i k denotes the mass transfer from phase i to phase k, whereas m ˙ k i stands for the mass transfer from phase k to phase i.

Momentum equation: ∂ ∂ t α k ρ k u → k + ∇ ⋅ α k ρ k u → k u → k = − α k ∇ p + α k ρ k g → + ∇ ⋅ α k τ k + τ k t + m ˙ i k u → i − m ˙ k i u → k + M → k where p is the pressure, which is assumed to be equal in all phases. g → denotes the gravity vector. τ k and τ k t are the molecular and turbulent stresses, respectively. M → k is the interphase momentum transfer.

Energy equation: ∂ ∂ t α k ρ k h k + ∇ ⋅ α k ρ k u → k h k = ∇ ⋅ α k λ k ∇ T k + μ k Pr t ∇ h k + ∂ ∂ t α k p + m ˙ i k h i − m ˙ k i h k + Q k where h k , λ k and T k are the enthalpy, thermal conductivity and temperature of phase k, respectively. μ k and Pr t are the turbulent viscosity and the turbulent Prandtl number. Q k is the heat transfer to phase k.

In the momentum equation, the interphase momentum transfer represents the sum of the forces that the phases exert on one another, which is composed of the drag force and the non-drag forces. The non-drag forces include the lift force, the turbulent dispersion force, the wall lubrication force and the virtual mass force. Many researchers convinced that the influences of the wall lubrication force and the virtual mass force could be ignored in the calculations (Li et al., 2011; Rabiee et al., 2017). Therefore, in this study, only the drag force ( F → D ), the lift force ( F → L ) and the turbulent dispersion force ( F → T D ) are considered: M → k = F → D + F → L + F → T D

The drag force is determined by the relative velocity between the liquid and vapor, and it is calculated as a function of the drag coefficient: F → D = C D 1 2 ρ l u → r A i 4 u → r where ρ l is the density of the liquid phase. u → r = u → l − u → g is the relative velocity between the continuous and dispersed phases. A i stands for the interfacial area. The drag coefficient C D is calculated by Tomiyama formulation (Tomiyama et al., 1998) in this study.

The lift force is calculated as follow: F → L = C L ρ l α g u → r × ∇ × u → l where u → l is the velocity of the liquid phase and C L is the lift coefficient. α g is the volume fraction of the vapor phase. The classical lift forc derived analytically C L = 0.5 using a simple shear model. However, this calculation was based on a weakly sheared inviscid flow and wake effects were not taken into account. It was reported by many studies that the lift coefficient became smaller for viscous flow. Morage et al. (1999) pointed out the necessity for a negative lift coefficient to fit certain experimental data. According to the theory proposed by Tomiyama et al. (2002), a positive lift coefficient indicates that the bubble will tend to migrate towards the wall while a negative lift coefficient indicates that the bubble will tend to migrate away from the wall. They also proposed a correlation of lift coefficient based on experiments on vertical tubes, which can reach a minimum value of −0.29. Hence, in order to reflect the non-uniformity effect introduced by mixing vanes of spacer grids, the lift coefficient is supposed to have a smaller value and it is set to C L = − 0.5 in this study.

The turbulent dispersion force is calculated as follow (burns et al., 2004): F → T D = C T D ν l t Pr t ∇ α l α l − ∇ α g α g = − C T D ν l t Pr t 1 α l + 1 α g ∇ α g = − C T D ν l t α l α g ⁡ Pr t ∇ α g where α l is the volume fraction of the liquid phase. ν l t is the turbulent kinematic viscosity of the liquid phase. C T D is a constant coefficient and set to 1 by default.

The classical numerical method for calculating flow boiling is the Rensselaer Polytechnic Institute (RPI) wall boiling model developed by Kurul and Podowski (1990). The RPI model partitions the total heat flux ( q t ″ ) into three components: the liquid phase convection heat flux ( q l ″ ), the evaporation heat flux ( q e ″ ) and the quenching heat flux ( q q ″ ). According to the Weisman-Pei bubble crowding model (Weisman and Pei, 1983), when the vapor volume along heated walls becomes high enough, the capability of the liquid to remove heat from the wall is restricted. The vapor phase convection heat flux ( q g ″ ) is taken into account and the original RPI wall boiling model is extended as follow: q t ″ = q l ″ + q e ″ + q q ″ 1 − K d r y + q g ″ K d r y where K d r y is the wall contact area fraction for the vapor that is estimated as: K d r y = 0 , α δ ≤ α d r y β 2 3 − 2 β , α δ > α d r y where β = α δ − α d r y 1 − α d r y where α δ is the vapor volume fraction averaged over the bubbly layer thickness. α d r y is the critical vapor volume fraction.

The four components of total heat flux are calculated as follow: q l ″ = ρ l c p l u l * T l + T W − T l q e ″ = N w f π d w 3 6 ρ g h lg q q ″ = 2 F A π d w 2 4 N w f ρ l c p l λ l t w π T W − T q q g ″ = ρ g c p g u g * T g + T W − T g where c p l and c p g are the specific heat capacities of the liquid and vapor phases, respectively. u l * and u g * are the near-wall velocities of the liquid and vapor, respectively. T l + and T g + are the dimensionless temperatures of the liquid and vapor, respectively. T W , T l and T g are the wall temperature as well as the liquid temperature and the vapor temperature close to the wall, respectively. T q is the quenching temperature, that is, the temperature at which liquid is brought to the wall by the departure of a bubble. N w , f and d w are the nucleation site density, the bubble departure frequency and the bubble departure diameter, respectively, and these parameters are derived by additional auxiliary models. h lg is the latent heat. λ l is the liquid conductivity. t w is the time elapsed between bubble departure and the nucleation of the next bubble. F A is an area scaling coefficient and is set to 2.0 according to Bartolomei and Chanturiya (1967).

Turbulence in the continuous phase is modelled using realizable two-layer k-epsilon model. For the dispersed phase, turbulence can be modelled with two alternative methods: either with full turbulence models or with response models. Many researchers recommended the Issa turbulence response model in their studies (Liu et al., 2021; Povolny and Cuhra, 2014; Yang P. et al., 2021), which used a semi-empirical correlation to link the turbulent fluctuations of the continuous and dispersed phases. Thus, aiming to improve computing efficiency, the Issa turbulence response model is employed in this study to model turbulence in the dispersed phase.

The detail settings of all the principal and auxiliary models for the CHF prediction are summarized in Table 1.

The geometric configuration utilized in this study was a 5 × 5 PWR fuel rod bundle assembled in a square channel as shown in Figure 1. The rod bundle consisted sixteen cold rods, eight hot rods and one guide thimble at the center. The cladding outer diameter was 9.5 mm and the cladding thickness was 0.6 mm. The thimble outer diameter was 12.5 mm. The pitch of the rod bundle was 12.6 mm.

A total of ten spacer grids were positioned along the 3.7 m heated length, including seven Mixing Grids (MG) with a height of 33 mm and three Mid Span Mixing Grids (MSMG) with a height of 18 mm, as shown in Figure 2. MG and MSMG were primarily composed of inner strips, outer strips, steel dimples, springs and mixing vanes. The strips arranged in a crisscross pattern formed a basic framework of the spacer grid. The steel dimples and springs maintained the fuel rods and the guide thimble in their proper radial position, preventing vibration under the impact of the fluid. Furthermore, the mixing vanes located at the end of the spacer grid were designed with a deflection angle with respect to the flow direction, so as to generate lateral flow and enhance fluid mixing. The inlet section was extended by an additional 0.3 m to ensure that the fluid was fully developed before the heating length. A schematic of the arrangement of spacer grids in the rod bundle is illustrated in Figure 3.

Due to the complexity of the spacer grids geometry, the flow domain where spacer grids were positioned was discretized with unstructured polyhedral meshes while a stretch mesh generator was employed to obtain extruded meshes for the fluid domain of rod bundles. In order to capture the flow details close to the wall, two layers of wall prism meshes were applied. For independence verification, five sets of meshes were tested under the same working condition. Figure 4 shows the variation of pressure drop with the number of cells. As can be seen, the pressure drop decreased gradually with the increase of the cell number. When the cell number exceeded 21.6 million, the influence on the calculation result of pressure drop could be neglected. Therefore, the mesh configuration with 21.6 million cells was selected for CHF prediction in this study. Figure 5 shows the final generated meshes of one of the regions of MG and MSMG, respectively.

Four experimental cases with significantly different working parameters were selected for the validation of the CFD analysis. Their working conditions are presented in Table 2, where p, G, and T _sub represent pressure, mass flux, and inlet subcooling, respectively.

As shown in Figure 1, fuel rods in the rod bundle were divided into two types: the sixteen fuel rods on the periphery were low-power rods, also known as cold rods, while the other eight fuel rods were high-power rods, also known as hot rods. The power ratio between a hot rod and a cold rod was 1.28. Having the heat flux shape identical to the experiment, the power distribution of fuel rods followed a cosine distribution in the axial direction. In addition, the guide thimble had no heat flux boundary and was set to adiabatic condition in CFD calculation. The heat flux distribution for all fuel rods was defined through the user-defined field function in STAR-CCM+.

The axial heat flux distribution applied to the inner surface of the cladding of cold rods was as follow: q c o l d ″ z = q r e f ″ × f z where z is the axial elevation and f is a normalized fitted cosine distribution function. q r e f ″ is named reference heat flux, which was progressively increased in CFD calculation.

The axial heat flux distribution applied to the inner surface of the cladding of hot rods was as follow: q h o t ″ z = 1.28 × q r e f ″ × f z

To calculate the total heating power, the axial heat flux distributions of all cold rods and hot rods were integrated and summed together as follow: Q = 16 × ∫ q c o l d ″ z π D i n d z + 8 × ∫ q h o t ″ z π D i n d z = 16 + 8 × 1.28 × q r e f ″ π D i n ∫ f z d z where D i n is the cladding inner diameter.

Then the average heat flux could be derived by dividing the total heating power by the total heating area as follow: q a v g ″ = Q 16 + 8 × π D o u t H h e a t where D o u t is the cladding outer diameter and H h e a t is the heating length.

In the CFD analysis, a method similar to that in the experiment was adopted. The reference heat flux q r e f ″ was set to a specific value and a steady calculation was carried out. Then important parameters, such as the wall superheat and the maximum volume fraction of vapor, were monitored. When the monitored parameters were stable and the calculation was convergent, the reference heat flux was increased to the next level and the next steady calculation was carried out. Such a process was repeated until the determination of the CHF. Between two steady calculations, the reference heat flux was increased in steps of 100 kW/m².

In the studies of CHF prediction using CFD, different methods have been used to determine whether DNB occurred. A common method is estimating the maximum void fraction of heating wall. DNB is considered to have occurred when the void fraction exceeds the critical value. However, for different researchers, the critical void fraction may take different values, such as 0.74 proposed by Podowski and Podowski (2009), 0.8 proposed by Vyskocil and Macek (2010b) and 0.82 proposed by Weisman and Pei (1983). Another common method is based on the wall temperature excursion. When a large temperature increase is monitored, the corresponding heat flux is considered as CHF (Yan et al., 2013; Karoutas et al., 2015). Considering that temperature excursions in numerical simulations could occur at low-mesh-quality cells due to numerical variations, Xu et al. (2019) proposed a new method by defining a variable called DNB indicator which was related to multiple thermal parameters. DNB occurred when the DNB indicator was larger than the indicator threshold, that need to be determined empirically with experimental cases. The studies of Li et al. (2018) and Liu et al. (2021) reported that the transition of the boiling curve could be used as an effective criterion of DNB occurrence. When plotting the curve of the wall superheat against the heat flux, in the subcooled boiling area, the wall superheat increases linearly with the increasing heat flux. As the heat flux continues to increase and reaches a certain value, a turn in the boiling curve can be noticed. This is because when the CHF is reached, the wall temperature will increase rapidly with the increasing heat flux, resulting in a rapid decrease in the slope of the boiling curve. In this study, the final method was adopted to determine the CHF values of the experimental cases.

The velocity distribution of liquid in the region between the inlet section ( z = 0   m ) and the beginning of heating length ( z = 0.3   m ) was assessed to verify the validity of the additional 0.3 m extension. Taking case C1 as an example, the contours of the velocity of liquid at different sections are shown in Figure 6. As can be seen, the liquid velocity distribution of the beginning of heating length had a high similarity with the section at 0.2 m elevation, which indicated that the fluid was fully developed.

A series of steady-state CFD calculations were conducted for each of the four experimental cases. Based on the calculation results, the curves of the wall superheat against the reference heat flux q r e f ″ of the experimental cases were plotted, as illustrated in Figure 7. As can be seen, for each case, the points of the boiling curve were initially located nearly in a straight line. Then the boiling curve turned when the wall superheat increased to a higher level. The transition of the boiling curve was taken as an indication of the occurrence of DNB, and the first point after the transition was considered the predicted CHF point. Therefore, the predicted values of the reference heat flux corresponding to the CHF points for cases C1 to C4 are 1,200 kW/m², 1700 kW/m², 1,600 kW/m² and 1900 kW/m², respectively.

For each case, based on the reference heat flux corresponding to the predicted CHF point, the corresponding average heat flux was derived, which represented the predicted CHF value. Predicted CHF values of all calculated cases were compared to the experimental CHF values, as shown in Figure 8. As can be seen, the predicted results agreed well with experimental ones. Deviations between the predicted and the experimental CHF values for the four cases were within 15%.

It is well known that the enhancement of the CHF by mixing vanes is attributed to the generation of swirl and cross flows and the increase in the turbulent intensity. Taking the CHF point of case C1 as an example (when reference heat flux q r e f ″ was set to 1,200 kW/m²), the influence of the mixing vanes was illustrated. Four observation sections, naming Section A to Section D, were defined as shown in Figure 9, which were located in the upstream and downstream areas of the 7^th and 8^th spacer grids, respectively. The contours of secondary flow distribution ( v x 2 + v y 2 ) in these sections were shown in Figure 10. As can be seen, the maximum secondary flow velocity in the cross-section at the inlet of the 7^th spacer grid was 0.70 m/s. When the fluid passed through the 7^th spacer grid, due to the effect of the mixing vanes, the secondary flow velocity at Section B increased significantly to 2.23 m/s. From the downstream area of the 7^th spacer grid to the upstream area of the 8^th spacer grid, the secondary flow velocity decreased to 0.77 m/s because the grid effect started disappearing. After flowing through the 8^th spacer grid, the secondary flow velocity at Section D increased again and reached 2.39 m/s. The presence of mixing vanes enhanced the lateral flow within the channel and in general, the location of CHF was most likely to be in the upstream area of the spacer grid.

Calculation results of case C1 were used as an example to verify the accuracy of the predicted CHF locations. The temperature values of all cells of hot rod surfaces along with the vertical location were extracted and plotted. Figure 11 shows the extracted temperature distribution of hot rod surfaces with the reference heat flux setting to 1,100 kW/m², which was the working condition before the CHF was reached. Each point in Figure 11 represented the value of one cell on the surface of the hot rods. The positions of all spacer grids were indicated in the figure. Data in different spans were shown in different colors for suitable distinguishment. It can be seen that due to the cosine shape axial power profile, the highest temperature did not occur at the outlet of the rod bundle channel, but in the region close to the axial power peak. It can also be noticed that after flowing through a spacer grid, the temperature would drop slightly.

Figure 12 shows the temperature of hot rod surfaces with the reference heat flux setting to 1,200 kW/m², where the CHF was reached. The shape of the temperature distribution was similar to that in Figure 11. However, a large temperature increase was observed in the upstream area of the 8^th spacer grid, indicating that the predicted CHF occurred at this location. Besides, temperature excursion could also be observed in the upstream area of the final spacer grid and in the end region of the rod bundle channel, but to a lesser extent, which indicated the possible location of the next CHF.

To visually assess the thermal-hydraulic parameter distributions around the CHF location, the contours of the temperature and the vapor volume fraction on the surface of the rod with the highest temperature in the channel close to the 8^th spacer grid were shown in Figures 13, 14, respectively. As can be seen, at the CHF location, the temperature was the highest and the vapor volume fraction was the larger. The larger vapor volume fraction was 0.9947.

The comparison between CFD simulation results and experimental data was summarized in Table 3. The deviation between predicted CHF values and experimental CHF values of the four cases was within 15%. All the predicted CHF locations were in the upstream area of a spacer grid, which was consistent with the experiment. The deviation of predicted and experimental CHF locations was within one grid-to-grid span length.