The OCR core designed in this study uses HB-40 (18wt% o-terphenyl, 82wt% hydrogenated terphenyl) as a coolant, an organic mixture that remains liquid at room temperature. HB-40 was used in the WR-1 Reactor in Canada between 1965 and 1985, where it would undergo a cracking reaction when exposed to high temperature along with radiation and would generate volatile gases as well as high boiling substances (Tegart, 1970). It is important to emphasize that HB-40 is assumed to remain stable in this research. Its physical properties only change with temperature, ignoring the cracking reaction or the effect of cracking products on the flow and heat transfer performance. The main properties of HB-40 are summarized as follows (Smee et al., 1975).

The density ρ (kg·m⁻³) as a function of temperature T (K) is ρ = 1.0181 × 10 3 − 7.3891 × 10 − 1 ( T − 273.15 ) . (1)

The viscosity μ (Pa s) over the temperature range 423–673 K is calculated as μ = EXP ( − 11.6260 + 2.1856 × 10 3 T − 1 ) . (2)

The formula for calculating the specific heat c p (J·kg⁻¹·K⁻¹): c p = 1.5049 × 10 3 + 3.5476 ( T − 273.15 ) . (3)

The formula of the thermal conductivity k (W·m⁻¹·K⁻¹) for temperatures from 398 to 673 K: k = 1.2677 × 10 − 1 − 1.0864 × 10 − 4 ( T − 273.15 ) . (4)

Because HB-40 triggers cracking when exposed to radiation, in order to decrease the fluid damage caused by fast neutrons and photons, the structure of the fuel element is optimized by adding a graphite layer to the original Zr-4 cladding as moderator, as shown in Figure 1. It should be emphasized here that the optimized fuel element is designed to be held in the under-moderated region, although the calculation details are not presented to avoid a lengthy article.

Considering that the thermal conductivity of graphite will be reduced by lattice defects due to neutron irradiation, an average effective thermal conductivity of 10 W/(m K) is assumed (Snead et al., 2002). In addition, because of the fine material compatibility of HB-40, which eliminates the possibility of fuel-coolant interaction in case of cladding rupture, UC is chosen as the fuel material for its greater thermal conductivity compared to UO₂. The former has a thermal conductivity of 23 W/(m K) at 1,000°C, while the latter is only about 2.6 W/(m K) (Yu, 2002). The use of UC can improve the heat transfer characteristics of fuel pellets to compensate for the bad thermal conductivity of HB-40.

The OCR core is rated at 100 MW(th) and consists of 37 fuel assemblies (FAs). Each assembly includes 264 fuel rods, 24 guide tubes for control rod movement, and 1 measuring tube, arranged in square arrays of 17 × 17. The active height of FAs is 200 cm, and the equivalent diameter of the core is 147.5 cm. The detailed dimensions of the core are shown in Table 1, and the steady-state operation parameters are shown in Table 2.

In order to flatten the power distribution, the FAs are loaded in three zones according to ²³⁵U enrichment. 13 FAs with enrichment of 1.95 w/o are arranged in the inner core, 16 FAs with enrichment of 4.25 w/o are arranged in the outer core, the remaining 8 FAs with enrichment of 2.75 wo/o are arranged in mid-place of the core.

To meet the demand of reactivity control, 21 of the 37 FAs contain control rod bundles with Ag-In-Cd as absorber material and 8 FAs contain burnable toxic rod bundles with borosilicate glass as absorber material. The arrangement of FAs is centrosymmetric, as shown in Figure 2. It must be emphasized here that all calculations in this study are performed at hot full-power operation, and all control rods are assumed to be withdrawn, also with a boric acid concentration of 0 ppm.

The Monte Carlo code MCNP5 is developed by Los Alamos National Laboratory, United States, which can treat an arbitrary three-dimensional configuration of materials in geometric cells. In this research, MCNP5 is used to build the physical model for the proposed core. Figure 3 shows the cross-sectional view of the core model. Besides the FAs, the baffle, barrel, and pressure vessel are also considered for the accuracy of neutronics calculations. F4 tally card is applied for neutron flux calculations, and F7 tally card is applied to obtain the radial and axial power distribution. In order to obtain axial power distribution, all FAs in the model are divided into 20 segments along the axial direction, and the power density of each segment in different FAs is calculated. For each calculation in MCNP, there are 330 generations, including 30 inactive generations, with an initial neutron number of 100,000 per generation, and a total of 30 million particle behaviors are recorded in every calculation.

The single-channel model is commonly applied for the steady-state thermohydraulics analysis of the core, which assumes that each coolant channel is isolated and closed, and there is no exchange of energy and momentum between the channels (Zhang and Huang, 2009). Figure 4 shows the schematic diagram of the 3 × 3 channels model and the single-channel model. Compared with the 3 × 3 channels model, only one channel needs to be taken for calculation in the single-channel model, which reduces the calculation cost. Based on the single-channel model, the temperature distributions of fuel elements and corresponding coolant channels in each assembly are calculated in turn.

Thermohydraulics boundary conditions are defined for the solid and fluid domains separately in COMSOL (Weng et al., 2021). The temperature boundary of 553 K and “outflow boundary” are applied to the coolant inlet and outlet, respectively. The four surrounding faces are “symmetric boundaries” and the upper and lower surfaces for the solid part are set as adiabatic. The volumetric heat source is added to the fuel pellet, which is the calculated axial power density distribution. The heat transfer field of both solid and fluid domains is coupled with the turbulent field of the fluid domain, and the contact boundary between the fluid and the solid is set as a “wall” by default. The inlet boundary condition of the turbulent field is set as the inlet flow rate, and the outlet boundary condition is set as a constant pressure of P ₀ = 1.15 MPa with backflow suppression.

Due to the irradiated swelling problem of UC, which will force the fuel pellet to come into contact with the inner wall of the cladding during reactor operation, a contact heat conduction model is used for the gap heat conduction between the fuel pellet and the Zr-4 cladding. The equivalent heat transfer coefficient at the interface of irradiated UO₂ fuel pellet and Zr-4 cladding is generally taken to be 5678 W/(m²·K) in the current fuel element design for pressurized water reactor power plants (Yu, 2002). Because the fuel element structures are similar and both UC and UO₂ suffer from irradiated swelling, the empirical gap equivalent heat transfer coefficient h _g = 5678 W/(m²·K) is adopted in this research to deal with the heat conduction of the gap.

The neutronics properties of each fuel rod in the same assembly are similar, so the coolant temperature and density distribution of each channel in the same assembly are considered to be approximately equal in the thermohydraulics calculations. The assembly power obtained by neutronics calculation can be converted into single-channel power in thermohydraulics calculation. And when the results of thermohydraulics calculation are fed back to neutronics calculation, the temperature and density parameters of single-channel will correspond to the parameters of assemblies in neutronics calculation. In addition, the coolant inlet flow rate in the single-channel model is determined by the coolant flow rate of the corresponding assembly, which is assigned according to the assembly power non-uniformity factor.

Three grid quantities are provided for the grid irrelevance analysis of the single-channel model. The average temperatures of the cross sections at different axial positions are calculated under different grid quantities, as shown in Figure 5. It is found that the simulated temperature values are basically the same under the same calculation conditions, which meet the requirements of grid irrelevance. In this study, the calculation is completed using the model with a grid number of 116,000.

In the coupled neutronics-thermohydraulics calculation, an initial temperature field is first assumed. We regard the average temperature of the coolant inlet and outlet as the initial coolant temperature, which is 573 K. The axial power density distributions of each assembly are obtained from the neutronics calculation, and the temperature fields of fuel elements and coolant channels in different assemblies are solved in turn by the thermohydraulics calculation. As the coolant density is determined by its temperature, the coolant density parameters of each assembly in MCNP input file can be updated after thermohydraulics calculation and the neutron transport process can be simulated again to obtain the new axial power density distribution of each assembly. The iterative calculation converges when the difference in the calculated power distribution is less than 5% compared with the previous step. Otherwise, the updated power density distribution will be imported into COMSOL for further calculation. The neutronics-thermohydraulics coupling calculation process is shown in Figure 6.

The coupling strategy between MCNP and COMSOL included mesh mapping and data transfer. In the core model built by MCNP, the fuel pellets and coolant channels in each assembly are divided into 20 segments along the axial direction. Each segment of the fuel pellet and coolant channel is considered a cell. Meanwhile, in the single-channel model built by COMSOL, the solid and fluid domains are divided into 20 corresponding segments along the axial direction to produce 20 fuel pellet cells and coolant channel cells, as shown in Figure 7. Make sure that the fuel pellet cells and coolant channel cells of different axial heights in the MCNP model can find overlapping cells in the COMSOL model. In the subsequent coupling calculation, the data of temperature, density, and power distribution would be transferred manually based on these overlapping cells.

Regarding the data transfer process, firstly, the power density of fuel pellet cells at different axial heights of each assembly is calculated by MCNP. Secondly, the data are imported to COMSOL manually as the heat source. Then, the temperature fields of fuel elements and coolant channels of each assembly are obtained by calculating COMSOL, and the average temperature and density of the coolant cells at different axial heights are derived. Finally, the obtained coolant density values at different axial heights of each assembly are used to update the original coolant density values of the corresponding assembly’s cells in the MCNP input file so that the axial power density distribution of each assembly can be calculated again.

The neutron flux for each assembly after the coupling calculation is shown in Figure 8 (only 1/4 of the core is shown due to symmetry). We can see that the closer to the core center, the higher the neutron flux. It is also worth noting that the fast neutron flux of each assembly is always higher than resonant neutrons and thermal neutrons. The total neutron flux of the core is calculated as 1.15 × 10¹⁴/cm²·s, of which the thermal neutron flux is 3.14 × 10¹³/cm²·s, accounting for 27.3%; the resonant neutron flux is 3.96 × 10¹³/cm²·s, accounting for 34.4%; and the fast neutron flux is 4.42 × 10¹³/cm²·s, accounting for 38.3%.

Finally, the normalized power distribution of the assemblies is given in the fifth row. As presented in Figure 8, the radial power peaking factor of the core is 1.239, and the assembly with the highest power density is D-4, located in the center of the core. Coolant channel with the highest power output is called the hottest channel, so the hottest channel is considered to exist in assembly D-4. The axial line power density distribution of D-4 before and after the coupling calculation is shown in Figure 9. After the coupling calculation, the axial power distribution of D-4 tends to be more like a cosine shape with less fluctuation, and the maximum line power density decreases to 9.10 kW/m. The axial power peaking factor of D-4 is 1.434. Therefore, the total power peaking factor of the core is 1.78.

Figure 10 shows the three-dimensional temperature for the single-channel model of D-4 assembly after the coupling calculation. Due to the uneven axial power distribution, the temperature in the center is significantly higher, with a maximum temperature of 946.8 K.

Regarding the steady-state thermohydraulics parameters of D-4, a comparison of the results before and after coupling is shown in Table 3. After the coupling calculation, the line power density of D-4 is reduced. Because the coolant flow rate of assemblies is distributed according to the power non-uniformity factor, the coolant flow rate of D-4 is consequently lower. Therefore, the calculated temperature distribution after coupling is basically the same. For example, the maximum temperature of the fuel pellet is only 7.6 K lower than that before coupling.

Figure 11 shows the temperature variation along the radial direction for the single-channel model of D-4 at axial height z = 1 m after coupling calculation. The curve shows three stages, whose structure in the order of temperature from the highest to the lowest are the fuel pellet, cladding (Zr-4 cladding and graphite cladding), and the coolant. There is a large temperature gradient between every two stages, and their formation is due to the existence of a gas gap and flow boundary layer, respectively.

For reactor design, the maximum heat flux at the surface of the fuel element is required to be less than the critical heat flux. Furthermore, the Departure from Nucleate Boiling Ratio (DNBR) concept is introduced to express the ratio of the critical heat flux q _DNB (W/m²) to the actual heat flux q (W/m²). The Atomic Energy of Canada Limited recommends the following correlation to calculate the critical heat flux of fresh HB-40 (Smee et al., 1975): q DNB = 1000 ( 5.1 Δ T sub V 2 3 − 409 ) . (5)

In Equation (5), Δ T s u b = T S A T − T B U L K is the sub-cooling degree of coolant, where T S A T and T B U L K are the coolant saturation temperature and the body temperature, respectively. T S A T = 0.88 × ln P − 0.9454 5.568 × 10 − 6 + 273.15 , where P is the coolant channel pressure in terms of bar and V is the coolant flow rate in terms of m/s.

The power distribution analysis indicates that the hottest channel is located in the D-4 assembly. The variation of the critical heat flux along the axial height of the fuel element in D-4 is calculated with the correlation. The variation of actual heat flux is exported from the single-channel model in COMSOL. The results are compared as shown in Figure 12. Influenced by the power distribution, the actual heat flux is high in the middle and low at both ends. Because the body temperature of coolant keeps increasing from inlet to outlet, the critical heat flux decreases along the axial height. Based on the ratio of the critical heat flux and actual heat flux, the variation of DNBR at rated power and steady state is obtained. The minimum departure from the nucleate boiling ratio is calculated as 2.41 at axial height z = 1.15 m, behind the maximum heat flux point, which is at axial height z = 0.95 m.