Numerical Investigation of Conjugated Heat Transfer of the Plate-Type Fuel Assembly in the Research Reactor

In nuclear reactors, the research of conjugated heat transfer between the fuel and coolant in the fuel assembly is fundamental for improving the safety, reliability and economy. The numerical approach based on Computational Fluid Dynamics (CFD) can be used to realize the rapid analysis of the conjugated heat transfer. Besides, the numerical simulation can provide detailed physical fields that are useful for the designing and optimizing of the fuel assembly. The plate-type fuels are generally used to enhance heat transfer in research reactors with high power density. In this study, a standard plate-type fuel assembly in the research reactor was taken into consideration. The solid-fluid conjugated heat transfer of the fuel assembly and coolant was numerically investigated. In the fluid region, the subcooled flow boiling simulation model was established by implementing the Rensselaer Polytechnic Institute model into the Euler multi-phase flow method. The results show that the conjugated heat transfer of the fuel assembly and coolant can be simulated using the model established in this work. The influence of fluid velocity, power density and the width of the flow channel on the temperature distribution and the conjugated heat transfer was investigated and discussed.


INTRODUCTION
Research reactors are widely used to produce high neutron flux for research, training, education, and irradiation test (Gong et al., 2015;Guo et al., 2018). The fuel assembly used in the research reactors is generally the plate-type fuel assembly, including several fuel plates, two side plates and narrow rectangular channels formed between the fuel plates. The research of heat transfer of the fuel assembly is fundamental for improving the safety, reliability and economy. However, the coolant flow and heat transfer characteristics in the narrow rectangular channels are significantly different from those of conventional channels and tubes (Sun et al., 2020). Therefore, investigation of coolant flow and heat transfer characteristics in the plate-type fuel assembly is important for research reactors designing.
In the last few years, 1-dimensional codes, 1.5-dimensional codes and 2-dimensional codes, such as the RELAP5/MOD3 (Lu et al., 2009;Son et al., 2015), FRAPCON (Phillippe et al., 2012), FEMAX (Yoo et al., 2006) and FROBA (Deng et al., 2016) were used for the thermo-mechanical coupling analysis and safety analyses of the fuels. In order to improve the accuracy of the calculation, 3dimensional simulations based on MOOSE (He et al., 2018) and COMSOL (Liu and Zhou, 2017) were used. However, the fluid was considered as 1-dimensional flow in the calculations and the detailed flow and heat transfer characteristics cannot be obtained.
The numerical investigation based on Computational Fluid Dynamics (CFD) can provide 3-dimensional detailed physical fields that are useful for the designing and optimizing of the fuel assembly. Recently, with the help of CFD, some researchers investigated the thermal-hydraulic characteristics of the fuel assembly.  numerically investigated the flow fields and phase distributions in the subchannel of a pressurized water reactor fuel assembly based on the CFD model of subcooled flow boiling. Kumar et al. (2018) and Zhang and Liu (2020) carried out CFD simulations to predict the void and temperature distribution inside the rod bundle. The influence of the spacer grid on the flow and heat transfer characteristics of subcooled flow boiling was studied by Zhang and Liu (2020). Shirvan (2016) numerically investigated the boiling crisis for helical cruciform-shaped rods using commercial software Star CCM+.
As for the plate-type fuel assembly, some researchers focused on the single-phase flow in the channels (Gong et al., 2015;González Mantecón and Mattar Neto, 2018;Liao et al., 2020). Gong et al. (2015) investigated the heat transfer characteristics of the plate-type fuel assembly. The single-phase turbulence flow was calculated using the commercial software FLUENT. The temperature distribution of the fuel assembly was obtained with different coolant inlet velocity and the assembly with a hot spot was specially studied. Liao et al. (2020) conducted the fluid-solid coupling simulation with the consideration of irradiation effects for plate-type fuel assembly. The thermo-hydraulic model in the fluid domain was established based on FLUENT. González Mantecón and Mattar Neto (2018) numerically investigated fluid-structure interaction of the plate-type fuel assemblies using the commercial software ANSYS CFX for modeling fluid flow and ANSYS Mechanical for modeling the plates. The maximum deflection of the plates was detected at the leading edge. An extra deflection peak was observed near the trailing edge of the plates for fluid velocities greater than the Miller's velocity.
The research of subcooled flow boiling heat transfer in the fuel assembly is fundamental for improving the safety, reliability and economy. However, fewer thermal-hydraulic investigations of the plate-type fuel assembly have taken the subcooled flow boiling into consideration. Guo et al. (2018) conducted the thermal-hydraulic analysis of flow blockage in the fuel assembly in the JRR-3M research reactor. The 3dimensional model of the fuel assembly was built and the flow and heat transfer characteristics were simulated using FLUENT. Guo et al. (2018) found when the blockage ratio of the fuel assembly is 70%, the departure from nucleate boiling will occur. Park et al. (2021) investigated the thermalhydraulic behaviors of flow channels of the plate-type fuel assembly by the CFD simulation with a multi-phase flow model. The wall boiling model proposed by Kurul and Podowski (1991) was used and the influence of the flow channel blockage on the flow instability and heat transfer was studied.
In this work, a standard plate-type fuel assembly in a research reactor was taken into consideration. The solid-fluid conjugated heat transfer of the fuel assembly and the coolant was numerically investigated. In the simulation, the Euler multi-phase flow method was used with the Rensselaer Polytechnic Institute (RPI) wall boiling model implemented. The temperature distribution and the flow field of the fuel assembly were obtained. The influence of power density, fluid inlet velocity and subchannel width on the conjugated heat transfer was investigated.

MATHEMATIC MODEL AND NUMERICAL METHOD
The simulation of conjugated heat transfer considering the phase change process is complicated. In this work, the conjugated heat transfer model was established. The simulation domain was divided into the solid region and the fluid region. The heat transfer of each region was calculated using different models and the heat transfer between the two regions was calculated.

Solid Region
In the solid region, heat is transferred by conduction. Since the thermal conductivity of the solid was assumed as constant, the heat transfer of the solid region was calculated by where ρ s is the density of the solid, c s is the specific heat of the solid, T is the temperature, t is time, k s is thermal conductivity of the solid, _ q s is the volume heat source.

Fluid Region
In the fluid region, the subcooled flow boiling was numerically simulated based on the Euler multi-phase flow method with liquid water being the primary phase and water vapor being the second phase. The volume fraction was used to represent the content of each phase and only two phases were set in the fluid region.
where α l is the liquid phase volume fraction, α v is the vapor phase volume fraction i.e., the void fraction. The governing equations included the mass conservation equations, the momentum conservation equations and the energy conservation equations for each phase. In order to simplified the description of the governing equations, the subscripts p and q were used denote the two phases.

1) Mass conservation equations
The mass conservation equation for phase q in the Euler multi-phase flow method was where ρ q and u q are the density and the velocity of phase q respectively, _ m pq and _ m qp are the mass transfer from phase p to phase q and the mass transfer from phase q to phase p.

2) Momentum conservation equations
The momentum conservation equation for phase q in the Euler multi-phase flow method was where τ q is the stress-strain tensor for phase q where μ q and λ q are the shear and bulk viscosity of phase q, p is pressure, g is gravity, _ m pq u p − _ m qp u q is the momentum transfer caused by phase change mass transfer and M q is the momentum transfer due to the interfacial forces.
The standard k-ε turbulence model was chosen to calculate the effect of turbulent flow and the two-layer model with all y+ wall function was employed. Generally, the interfacial momentum transfer between liquid and vapor should be the sum of the drag force, turbulent dispersion force, lift force, wall lubrication force and virtual mass force.
The sensitivity analysis of the forces on the simulation result of subcooled flow boiling was conducted. The simulation object and  conditions were the same as those in the experiment conducted by Zeitoun and Shoukri (1996). The channel geometry was an annulus with a 306 mm length. The outer wall diameter was 25.4 mm and the heated inner wall diameter was 12.7 mm. The heat flux of the inner wall of the channel was 508 kW/m 2 , the mass flux of the coolant was 264.3 kg/(m 2 s), the coolant inlet temperature was 94.6°C and the system pressure was 150 kPa. The result shows that the influence of the lift force, wall lubrication force and virtual mass force on the calculated wall temperature was less than 1 K. Similar to the results of the sensitivity analysis we finished in our previous work (Li et al., 2018), the effects of the three forces are to move the vapor bubbles radially away from or towards the heated surface. While the forces have little influence on the void fraction distribution along the coolant flow direction.
Since the lift force, wall lubrication force and virtual mass force have little effect on the simulation results and the inclusion of the three forces takes more time to obtain convergence, only the drag force and turbulent dispersion force were taken into consideration in this work.
where F D,k is the drag force and F TD,k is the turbulent dispersion force.The drag force can by calculated using where u r u l -u v is the relative velocity between the two phases, d B is the bubble diameter and C D is the drag force coefficient. In this work, the bubble diameter was calculated using the Anglart's method (Anglart and Nylund, 1996) and the drag force coefficient was calculated using the Tomiyama's method (Tomiyama, 2004). The turbulent dispersion force plays an important role in taking the vapor bubbles from the near wall region to the bulk liquid region. The turbulent dispersion force was calculated using the method proposed by Lopez de Bertodano (1991).
where C TD is the turbulent dispersion coefficient and k l is the turbulent kinetic energy of the liquid phase.

3) Energy conservation equations
The energy conservation equation for phase q in the Euler multi-phase flow method was where h q is the specific enthalpy of phase q, q q is the heat flux, h qp and h pq are the interphase enthalpy.

Solid-Fluid Boundary
At the solid-fluid boundary, heat transferred from the solid region to the fluid region.
where _ q w fluid is the heat flux of the fluid region and _ q w solid is the heat flux of the solid region at the solid-fluid boundary.
Wall Boiling Model Figure 1 shows the schematic of wall boiling heat flux partitioning. The Rensselaer Polytechnic Institute (RPI) wall boiling model (Kurul and Podowski, 1991) divided the total heat flux from the hot wall to the fluid region _ q w to the liquid convection heat flux _ q c , the quenching heat flux _ q q and the evaporation heat flux _ q e . Lavieville et al. (2005) modified the RPI model by considering the heat flux directly from the hot wall to the vapor _ q v . Thus, the total heat flux can be calculated by where K dry is the wall contact area fraction for vapor. K dry can be calculated by where α δ is the average void fraction in the mesh near the hot wall, α dry is the critical void fraction.The four parts heat flux can be calculated by where u + l and u + v are the frictional velocities of liquid and vapor near the hot wall, T + l and T + v are the non-dimensional parameters for liquid and vapor, K quench is the bubble influence wall area fraction and was calculated using the model proposed by Kurul and Podowski (1991), t w is the bubble waiting time, D d and f d are the bubble departure diameter and bubble departure frequency, N w is the active nucleation site density. In this work, the widely used active nucleation site density model proposed by Lemmert and Chawla (1977), the bubble departure diameter model proposed by Tolubinsky and Kostanchuk (1970) and the bubble departure frequency model proposed by Cole (1960) were used. In Tolubinsky's model, D 0 was 0.6 mm.

Model Validation
Due to the lack of experimental results of subcooled flow boiling in the plate-type fuel assembly, the conjugated heat transfer model was validated using the experimental results of subcooled flow boiling in a tube conducted by Bartolemei and Chanturiya (1967) and the experimental results of subcooled flow boiling in an annulus channel conducted by Zeitoun and Shoukri (1996). In Bartolemei's   work, the tube diameter was 15.4 mm, the heat flux of the tube wall was 0.57 MW/m 2 , the mass flux of the coolant was 900 kg/(m 2 s), the subcooled temperature of the coolant was 60 K and the system pressure was 4.5 MPa. The coolant temperature, wall temperature and void fraction were measured by Bartolemei and Chanturiya (1967). Since the system pressure in Bartolemei's work was larger than that in our simulation for plate-type fuel assembly, the experimental results at lower system pressure obtained by Zeitoun and Shoukri (1996) were also used to validate the model. Figures 2A,B show the comparison of the numerical results of the coolant temperature, wall temperature and void fraction along the coolant flow direction with the experimental results obtained by Bartolemei and Chanturiya (1967). Figure 2C shows the comparison of the numerical results of void fraction along the coolant flow direction with the experimental results obtained by Zeitoun and Shoukri (1996). It can be seen that the calculated results agree well with the measured data and the calculated void fraction is slightly larger than the measured data at lower system pressure.

PHYSICAL MODEL
In this work, a standard plate-type fuel assembly of the research reactor was chosen as the simulation object. The structure, material and thermal-hydraulic parameters of the fuel assembly were the same as those used in Gong's work (Gong et al., 2015). Table 1 shows the structure, material and thermal-hydraulic parameters of the standard plate-type fuel assembly.
A three-dimensional model of the standard plate-type fuel assembly was built on the Star CCM+ platform. Figure 3 shows the simulation domain of the conjugated heat transfer of the plate-type fuel assembly. The simulation domain included the inlet section, the test section and the outlet section. In the test section, there were 20 fuel plates with 1.52 mm thickness and 19 subchannels with 2.28 mm width and two by-pass subchannels with 1.24 mm width. The fuel, cladding and coolant parts were considered in the simulation domain. As shown in Figure 3B, half of the fuel assembly was included in the simulation domain due to the symmetry of the geometry and the fluid field.
The hexahedral mesh was used for the physical model. In the mesh independence analysis, a single calculation unit of the assembly was used. A flow channel with two half of the fuel plates was included in the unit. The power distribution of the fuel was assumed to be uniform with 2.35 × 10 9 W/m 3 . The coolant inlet temperature, inlet velocity and the system pressure were 35°C, 4.5 m/s and 152 kPa, respectively. Figure 4 shows the calculation results of coolant outlet temperature and the coolant temperature at the center of the calculation region with different number of grid cells. The calculation results of the coolant outlet temperature change a little with the number of grid cells due to the energy conservation. The coolant temperature at the center of the calculation region changes a little when the number of grid cells exceeds 180,000. Thus, the mesh with a 189,000 grid cells was used in this work and the same meshing method was used for the whole fuel assembly. The cell size was chosen as Δx ≈ 0.23 mm, Δy ≈ 1.1 mm, Δz ≈ 2 mm in the fluid domain and Δx ≈ 0.19 mm, Δy ≈ 2.4 mm, Δz ≈ 3 mm in the solid domain.

Conjugated Heat Transfer of the Assembly
The conjugated heat transfer of the plate-type fuel assembly was carried out. The power distribution of each plate along height direction was assumed to be a cosine distribution as where q max is the maximum power density and was set as 3.654 × 10 9 W/m 3 , H is the height of the fuel.  Frontiers in Energy Research | www.frontiersin.org September 2021 | Volume 9 | Article 663533 Figure 5 shows the temperature distribution of the fuel plates and Figure 6 shows the temperature distribution of the coolant.
The maximum temperature appears near the center of the fuel plate due to the non-uniform power density distribution of the fuel. Since the power density is zero, the temperature near the ends of the fuel plate equals to the coolant temperature approximately.

Conjugated Heat Transfer of a Single Unit
In order to reduce the time consuming of the simulation, the conjugated heat transfer of a single unit was investigated in this section. The influence of power density, inlet velocity and channel width on the conjugated heat transfer was investigated. Figure 7 shows the coolant outlet temperature, maximum coolant temperature, maximum cladding temperature and maximum fuel temperature with different inlet velocities. The uniform power density of the fuel was 2.35 × 10 9 W/m 3 and the coolant inlet temperature was 35°C. As shown in Figure 7, the temperatures of the coolant, cladding and fuel decrease with the increase of the inlet velocity due to the increase of the coolant mass flux. Within the inlet velocity range of 1.5-7.5 m/s, the maximum coolant temperature was lower than the saturation temperature of the coolant. The fuel plate can keep safe since no local nucleate boiling occurs. The difference between the coolant temperature and the cladding temperature is the smallest with 7.5 m/s inlet velocity due to the enhancement of single phase convective heat transfer at the highest inlet velocity.

Influence of Inlet Velocity
With 0.5 m/s inlet velocity, the maximum coolant temperature exceeds the saturation temperature of the coolant, which means the boiling phase change occurs in the flow channel. As shown in Figure 7, the difference between the coolant temperature and the cladding temperature is small with 0.5 m/s inlet velocity. This is because the nucleate boiling occurs and the heat transfer between the coolant and the cladding was enhanced by nucleate boiling. Figure 8 shows the void fraction distribution in the flow channel with 0.5 m/s inlet velocity. As shown in Figure 8, vapor appears near the hot fuel plate. The vapor fraction increases along the flow direction due to the increase of the coolant temperature. The maximum void fraction appears at the leeward of the plate due to the small velocity. Figure 9 shows the coolant outlet temperature, maximum coolant temperature, maximum cladding temperature and maximum fuel temperature with the different uniform power density of the fuel. The inlet velocity was 4.5 m/s and the coolant inlet temperature was 35°C. The coolant temperature, the cladding temperature and the fuel temperature increase with the power density as shown in Figure 9. Since the maximum coolant temperature reaches above the saturation temperature with 9.40 × 10 9 and 11.75 × 10 9 W/m 3 power density, the boiling phase change occurs. As shown in Figure 10, due to the high power density, the void fraction near the fuel plate at the rear of the flow channel was very large, which  weaken the heat transfer and lead to higher maximum fuel temperature and higher maximum cladding temperature. The influence of power density distribution on the conjugated heat transfer was numerically investigated. Three kinds of power density distributions along the height direction were used with the total power of the fuel plate maintained, as shown in Figure 11. Figure 12 shows the simulation results of the temperature distribution of the calculation region with q avg 2.35 × 10 9 W/ m 3 . As shown in Figure 12, the coolant outlet temperatures were equal for these three cases due to the conservation of energy. The maximum temperature of the cladding and the fuel is the largest for q 2 due to the concentrate distribution of the power density. The maximum cladding temperature and the maximum fuel temperature were lower for the uniform distribution of the power density. With the uniform power density distribution, the maximum temperature of the fuel plate appears at the rear due to the increase of the coolant temperature along the flow direction.

Influence of Channel Width
The conjugated heat transfer with different channel widths was simulated. In the simulations, the coolant inlet mass flux was maintained. Figure 13 shows the coolant outlet temperature, maximum coolant temperature, maximum cladding temperature and maximum fuel temperature with different channel widths, 1.75, 2.28 and 2.75 mm. Since the coolant inlet mass flux was maintained, the inlet velocities were 5.23, 4.50, and 4.00 m/s for these three cases, respectively. As shown in Figure 13, the coolant outlet temperature does not change with the channel width due to the conservation of energy. The maximum cladding temperature and maximum fuel temperature increase with the channel width. This is because, with the channel width increase, the coolant velocity decreases and the convective heat transfer weakens.

CONCLUSION
The conjugated heat transfer model between the plate-type fuel assembly and the coolant was established in this paper. With the RPI wall boiling model added into the Euler multi-phase flow method, the heat transfer and boiling phase change process in the whole fuel assembly and a single unit of the assembly of the research reactor were simulated. In the model validation, the simulated wall temperature variation, coolant temperature variation and void fraction variation were in good agreement with the experimental results in the literature.
The influence of the coolant inlet velocity, power density and flow channel width on the conjugated heat transfer was numerically investigated. According to the simulation results, the coolant outlet temperature, maximum coolant temperature, maximum cladding temperature and maximum fuel temperature increase with the decrease of the inlet velocity and increase with the increase of the fuel power density. Three kinds of power density along the axial direction were considered in this work. The maximum cladding temperature and the maximum fuel temperature were lower for the uniform distribution of the power density. As for the investigation of the influence of channel width on the conjugated heat transfer, the coolant inlet mass flux was maintained. The maximum cladding temperature and maximum fuel temperature increase with the increase of channel width.
Among the simulation conditions in this work, the boiling phase change occurs with 0.5 m/s inlet velocity under 2.35 × 10 9 W/m 3 power density, 4.5 m/s inlet velocity under 9.4 × 10 9 W/m 3 power density and 4.5 m/s inlet velocity under 11.75 × 10 9 W/m 3 power density. The maximum void fraction appears at the leeward of the plate and the rear of the flow channel near the plate.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.