Neutronics and Thermohydraulics Coupling Analysis on Novel Organic Cooled Reactor Based on Single-Channel Model

Introduction

Small modular reactors have become a research hotspot in the field of nuclear energy, characterized by low power (less than 300 MWe), modular design and construction, outstanding safety, and economic performance. Among all worldwide designs of small modular reactors, an organic cooled reactor (OCR) is considered one of the most promising concepts (Shirvan et al., 2016).

OCR is a reactor concept that uses an organic fluid as a coolant, allowing high temperature and low-pressure operation, along with negligible corrosion on loop structure materials (Elberg and Fritz, 1963; International Atomic Energy Agency, 1967). OCR is compact, safe, and reliable, which can be used in many fields such as desalination, district heating, and hydrogen production. Since the 1950s, the United States, the Soviet Union, Canada, and other countries have built a number of test OCRs for different purposes (Binstock, 1960; Mccurnin, 1963; Hewson, 1965; Tegart, 1970; Tsykanov et al., 1981). In recent years, the advent of high-temperature fluids, together with advances in hydrocracking and reforming technologies driven by the oil and gas industries, make the OCR concept even more viable. In 2014, the Massachusetts Institute of Technology proposed a core design to be deployed offshore underwater for powering coastal areas, which was moderated by graphite and cooled by an organic fluid (Shirvan and Forrest, 2016). In 2019, the Royal Military College of Canada undertook a conceptual design of an organic cooled and moderated reactor for heating and powering remote areas of northern Canada (Shannon et al., 2020).

The current research on OCR is mostly limited to neutronics analysis (li et al., 2021; Wang et al., 2021). Little research has been carried out on thermohydraulics. In this study, a conceptual design of an OCR core with thermal power of 100 MW(th) is proposed, and the organic fluid HB-40 is used as a coolant. The Monte Carlo particle transport code MCNP is used to develop the physical model of the core for neutronics calculations, and the multi-physics simulation software COMSOL is used to develop the single-channel model for thermohydraulics calculations. By coupling the MCNP code and COMSOL software through mesh mapping and data transfer, a neutronics-thermohydraulics coupling method is established, and the upcoming research on thermohydraulics characteristics is carried out on this basis.

Core Design

Organic Coolant Properties

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 $\rho$ (kg·m⁻³) as a function of temperature T (K) is

$\mathit{\rho }=1.0181\times {10}^3-7.3891\times {10}^{-1}\left(\mathit{T}-273.15\right).(1)$

The viscosity $\mu$ (Pa s) over the temperature range 423–673 K is calculated as

$\mathit{\mu }=\mathit{EXP}\left(-11.6260+2.1856\times {10}^3{\mathit{T}}^{-1}\right).(2)$

The formula for calculating the specific heat $c_p$ (J·kg⁻¹·K⁻¹):

${\mathit{c}}_{\mathit{p}}=1.5049\times {10}^3+3.5476\left(\mathit{T}-273.15\right).(3)$

The formula of the thermal conductivity k (W·m⁻¹·K⁻¹) for temperatures from 398 to 673 K:

$\mathit{k}=1.2677\times {10}^{-1}-1.0864\times {10}^{-4}\left(\mathit{T}-273.15\right).(4)$

Optimal Design of Fuel Element

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.

FIGURE 1

FIGURE 1. Schematic diagram of the optimized fuel element.

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.

Parameters of the Core Design

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.

TABLE 1

TABLE 1. Parameters of the core structure.

TABLE 2

TABLE 2. Parameters at steady-state operation.

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.

FIGURE 2

FIGURE 2. Schematic diagram of the arrangement of 37 FAs. Different types are shown with different colors (dark blue, shutdown control rod assembly bank; blue, regulating control rod assembly bank; yellow, burnable toxic rod assembly bank; white, assemblies that contain neither control rods nor burnable toxic rods).

Methodology of Coupling

Physical Model

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.

FIGURE 3

FIGURE 3. Cross-sectional view of the physical model.

Single-Channel Model

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.

FIGURE 4

FIGURE 4. Schematic diagram of (A) 3 × 3 channels model (B) single-channel model.

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.

Grid Irrelevance Verification

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.

FIGURE 5

FIGURE 5. Axial temperature distribution of the single-channel model with different numbers of grids.

Neutronics-Thermohydraulics Coupling Calculation Method

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.

FIGURE 6

FIGURE 6. Flowchart of the coupling procedure.

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.

FIGURE 7

FIGURE 7. Schematic diagram of the meshes used for data transfer.

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.

Results and Discussions

Power Density Distribution

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%.

FIGURE 8

FIGURE 8. Distributions of neutron fluxes and normalization power in each FA.

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 9

FIGURE 9. Distribution of axial line power density and normalized power density of D-4 assembly.

Temperature Distribution

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.

FIGURE 10

FIGURE 10. Three-dimensional temperature distribution for the single-channel model of D-4 assembly.

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.

TABLE 3

TABLE 3. Steady-state thermohydraulics parameters of D-4 assembly.

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.

FIGURE 11

FIGURE 11. Radial temperature distribution of D-4 assembly (z = 1 m).

Minimum Departure From Nucleate Boiling Ratio

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):

${\mathit{q}}_{\mathit{DNB}}=1000\left(5.1\mathrm{\Delta }{\mathit{T}}_{\mathit{sub}}{\mathit{V}}^{\frac{2}{3}}-409\right).(5)$

In Equation (5), $\Delta T_{sub}=T_{SAT}-T_{BULK}$ is the sub-cooling degree of coolant, where $T_{SAT}$ and $T_{BULK}$ are the coolant saturation temperature and the body temperature, respectively. $T_{SAT}=0.88\times \sqrt{\frac{\ln P-0.9454}{5.568\times {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.

FIGURE 12

FIGURE 12. Axial variation of heat flux on the fuel surface.

Conclusion

In order to obtain accurate neutronics and thermohydraulics parameters of the OCR core, a coupled neutronics-thermohydraulics calculation method is established based on the single-channel model, applying the Monte Carlo code MCNP and the multi-physics software COMSOL. By analyzing the results of coupled calculations, the conclusions are as follows.

1) For the proposed conceptual OCR design, the neutronics-thermohydraulics calculations yielded a total power peaking factor of 1.78, a total neutron flux of 1.15 × 10¹⁴/cm²·s, a maximum channel temperature of 947 K, and a MDNBR value of 2.4, which are within the specification of the neutronics-thermohydraulics design basis limits.

2) After the coupling calculation, the line power density of the assembly, where the hottest channel is located, is reduced. However, the temperature distribution is basically the same because the coolant flow rate of assemblies is distributed according to the non-uniformity power factor.

This research provides a reference for future design and optimization of OCR, but further analyses should consider the effects of cracking products on organic fluid properties to predict the neutronics-thermohydraulics characteristic more accurately.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors without undue reservation.

Author Contributions

XG is responsible for code programming, data processing, and editing. FW is responsible for concept and project management. DG is responsible for editing. WD is responsible for editing.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors acknowledge the support of Science and Technology on Reactor System Design Technology Laboratory (JG2018119). They also gratefully acknowledge the financial support provided by the National Key R&D Program of the People’s Republic of China (2018YFB0104502).

References

Binstock, M. H. (1960). Fuel Element Development for piqua OMR. Canoga Park, CA: Atomics International.

Google Scholar

Elberg, S., and Fritz, G. (1963). “Physical properties of organic nuclear reactor coolants: EUR 400.e,” in 144th National ACS Meeting Division of Industrial and Engineering Chemistry. Symposium on Organic Nuclear Reactor Coolant Technology. (Los Angeles, CA: European Atomic Energy Community).

Google Scholar

Hewson, R. A. (1965). Compilation of piqua Nuclear Power Facility Operating Limits and Controls and post-critical Operational Tests. Canoga Park, CA: Atomics International.

Google Scholar

International Atomic Energy Agency (1967). Organic Liquids as Reactor Coolants and Moderators. Vienna: IAEA.

Google Scholar

li, Q., Xia, B., Li, S., and Lu, D. (2021). Study on the Concept of Organic-Cooled Microreactor. Nucl. Power Eng. 42 (5), 266–270. doi:10.13832/j.jnpe.2021.05.0266

CrossRef Full Text | Google Scholar

Mccurnin, W. R. (1963). Evaluate the Performance of Two OMRE Fourth Core Fuel Elements: NAA-SR-9008. Canoga Park, CA: Atomics International.

Google Scholar

Shannon, C., Chan, P., and Bonin, H. W. (2020). Conceptual Design of an Organic-Cooled Small Nuclear Reactor to Support Energy Demands in Remote Locations in Northern Canada. CNL Nucl. Rev. 9, 39–44. doi:10.12943/CNR.2019.00002

CrossRef Full Text | Google Scholar

Shirvan, K., Ballinger, R., Buongiorno, J., Forsberg, C., Kazimi, M., and Todreas, N. (2016). Technology Selection for Offshore Underwater Small Modular Reactors. Nucl. Eng. Techn. 48 (6), 1303–1314. doi:10.1016/j.net.2016.06.002

CrossRef Full Text | Google Scholar

Shirvan, K., and Forrest, E. (2016). Design of an Organic Simplified Nuclear Reactor. Nucl. Eng. Techn. 48 (4), 893–905. doi:10.1016/j.net.2016.02.019

CrossRef Full Text | Google Scholar

Smee, J. L., Puttagunta, V. R., Robertson, R. F. S., and Hatcher, S. R. (1975). Organic Coolant Summary Report. Pinawa, Manitoba: Atomic Energy of Canada Ltd.

Google Scholar

Snead, L. L., Balden, M., Causey, R. A., and Atsumi, H. (2002). High thermal Conductivity of Graphite Fiber Silicon Carbide Composites for Fusion Reactor Application. J. Nucl. Mater. 307-311, 1200–1204. doi:10.1016/S0022-3115(02)01115-7

CrossRef Full Text | Google Scholar

Tegart, D. R. (1970). “Operation of the WR-1 Organic Cooled Research Reactor,” in American Nuclear Society Conference on Reactor Operating Experience (Pinawa, Manitoba: Atomic Energy of Canada Ltd.).

Google Scholar

Tsykanov, V. A., Chechetkin, Y. V., Kormushkin, Y. P., Polivanov, I. F., Pochechura, V. P., Yakshin, E. K., et al. (1981). Experimental Nuclear Heat Supply Station Based on the Arbus Reactor. At Energy 50 (6), 333–338. doi:10.1007/BF01126338

CrossRef Full Text | Google Scholar

Wang, F., Sun, Y., and Su, X. (2021). Novel Organic Nuclear Reactor Fuel Assembly Physical Calculation and Reactor Core. Mater. Rev. 35 (18), 18195–18199. doi:10.11896/cldb.20070337

CrossRef Full Text | Google Scholar

Weng, M., Liu, S., Liu, Z., Qi, F., Zhou, Y., and Chen, Y. (2021). Development and Application of Monte Carlo and COMSOL Coupling Code for Neutronics/thermohydraulics Coupled Analysis. Ann. Nucl. Energ. 161, 108459. doi:10.1016/j.anucene.2021.108459

CrossRef Full Text | Google Scholar

Yu, P. (2002). Thermal Analysis of Nuclear Reactor. Beijing: Atomic Energy Press.

Google Scholar

Zhang, X., and Huang, Y. (2009). “Single-Channel Model for Steady Thermal-hydraulic Analysis in Nuclear Reactor,” in 2009 Annual Symposium, Chinese Nuclear Society (Beijing: Atomic Energy Press), 254–259.

Google Scholar