Abstract
The cold neutron source (CNS) system of the Open Pool Australian Light-Water (OPAL) reactor is a 20 L cryogenically cooled liquid deuterium thermosiphon system. The CNS is cooled by forced convective helium which is held at room temperature during stand-by (SO) mode and at ~20 K during normal operation (NO) mode. When helium cooling stops, the reactor is shut down to prevent the moderator cell from overheating. This computational fluid dynamics (CFD) study aims to determine whether the combined effects of conduction and natural convection would provide sufficient cooling for the moderator cell under the influence of reactor decay heat after reactor shutdown. To achieve this, two commercial CFD software packages using an identical CFD mesh were first assessed against an experimental heat transfer test of the CNS. It was found that both numerical models were valid within the bounds of experimental uncertainty. After this, one CFD model was used to simulate the thermosiphon transient condition after the reactor is shut down. Two independent shutdown conditions of different decay-heat power profiles were simulated. It was found that the natural convection and conduction cooling in the thermosiphon were sufficient for dissipating both decay-heat profiles, with the moderator cell remaining below the safe temperature of 200°C.
1. Introduction
Using neutrons for the study of materials has advanced greatly over the last 60 years and the introduction of cold neutron sources (CNS) have allowed the investigation of long-range atomic structures found in polymers, ceramics, and other functional materials. Ever since Butterworth et al. [1] built the first CNS at Harwell, UK in 1957, CNS moderator media have evolved from ice sources (composed of CH4 or H2O), to hydrogen sources, to optimised liquid deuterium sources found at FRM-II [2] in Germany, Institut Laue-Langevin (ILL) [3] in France and the OPAL [4] reactor at the Australian Nuclear Science and Technology Organisation (ANSTO).
The principle of a reactor driven CNS is relatively simple. Neutrons streaming out of the reactor enter a low-temperature moderating volume positioned immediately adjacent to the core. The neutrons are thermalised to the Maxwellian energy distribution of the cold moderator which is usually cooled to ~20 K. Some of these neutrons travel down the neutron guide that connects the CNS with the neutron instruments. To ensure the rapid cooling of these neutrons, the moderator must operate efficiently, preferably with maximal neutron scattering and minimal parasitic neutron capture. Deuterium is an excellent moderator with a minimal neutron capture cross-section compared with hydrogen. Thus, liquid deuterium has become the favoured CNS moderating medium.
Many of the initial studies into the use of liquid deuterium for cold neutron science were pioneered by the French in the late 60s. Harig [5] published in French experimental results comparing hydrogen and deuterium as the moderator media for a cold neutron source installed in the SILOETTE reactor. The purpose of these experiments was to determine the optimal CNS design that was to be installed in the French-German High Flux Reactor (HFR) at ILL at the time. Ageron et al. [6] surmised this work in an English publication which concluded that the liquid deuterium source gave the best overall performance with the highest flux and gain factor at wavelengths greater than 8 Å. This finding was repeated by Egelstaff [7] who examined many CNS benchmarking studies and concluded that deuterium, in a sufficiently large volume, was the best moderator for cold neutrons. Despite the findings of Ageron and Egelstaff, hydrogen moderated cold neutron sources outnumber deuterium source on the basis of cost and compactness.
Grósz et al. [8] made a more recent study of a hydrogen-based cold neutron source for a small (0.4–0.5 L), double-walled, helium cooled moderator cell fitted at the end of the beam tube of the Budapest Research Reactor (BRR). Grósz et al. studied two shapes: a disc and a cylinder for the moderator cell and determined the cylinder to be the more favourable configuration for heat removal. The study was followed up with an experimental investigation of the moderator cell’s heat transfer characteristic [9] and a final report on the operation of the cold neutron source at BRR published in 2002 [10].
The South Korean HANARO reactor [11] also has a hydrogen-based CNS which has a total H2 inventory of 355 grams. The Hanaro CNS, similar to the CNS at NIST [12], operates with a two-phase flow hydrogen moderator. The use of hydrogen is popular as it is relatively inexpensive compared to deuterium. Hydrogen has both a larger scattering cross-section and a higher capture cross-section than deuterium which places an upper limit to the volume of H2 installed lest too many neutrons be lost. Thus, H2 systems usually have a smaller moderator inventory compared to deuterium systems and require less cryogenic cooling. Thermal analysis of two-phase H2-based CNS systems are predominantly analytical as CFD simulations of H2/D2 two-phase flows with heat and mass transfer (boiling) are somewhat challenging. For example, Li et al. [13] used a simplified analytical model to analyse a proposed two-phase H2 CNS for the Chinese Advanced Research Reactor (CARR) as did Han et al. [14] for the Hanaro CNS.
More recently a 12 L deuterium CNS was installed at FRM-II Germany and OPAL was commissioned with a 20 L deuterium CNS; both designs use a thermosiphon for moderator cooling. In the near future, NIST plans to upgrade its hydrogen CNS with a deuterium source. Between FRM-II and OPAL, the FRM-II CNS differs by operating with a two-phase flow regime whereas the OPAL CNS operates in a completely liquid phase. A novel design feature of the FRM-II CNS is the inclusion of a small percentage (<10%) of hydrogen in the deuterium to enhance neutron moderation.
Heat and mass transfer modelling of the cold neutron source is of great interest to CNS designers who aim to optimise the heat transfer characteristics of the CNS system. By using computational fluid dynamics (CFD) as a design tool, moderator cooling can be maximised in a system that is both simple and cost effective. The use of CFD for CNS design often begins with subsystem analyses focusing on the CNS moderator cell containing the moderating medium. For example, Sun et al. [15] used CFD to calculate the steady state temperature field of the moderator cell in the China Advanced Research Reactor (CARR) reactor during stand-by operation. The temperature field was then used as the input parameter set for a Finite Element Analysis (FEA) calculation of moderator cell stress distribution to confirm that local stress concentrations did not exceed material design limits.
The Open Pool Australian Light-Water (OPAL) reactor, located at the Lucas Heights Research Laboratory in Sydney, has a vertically installed, 20 L, liquid deuterium cold neutron source. The liquid deuterium circulates in a closed thermosiphon loop cooled by gaseous cryogenic helium. The moderator cell, heat exchanger, and connecting tubes are all surrounded by a helium jacket. Both the helium leg that circulates around the moderator cell and the helium leg cooling the heat exchanger are supplied by a common Brayton cycle cryogenic system. The moderator cell is fabricated in aluminium to maximise neutron transparency, while the heat exchanger, located above and away from the reactor core, is made of stainless steel.
The transient moderator temperature of the OPAL CNS was first modelled by Buscaglia et al. [16] in 2004. The aim of the calculation was to check that there were no “flow dead spots” which could attain a temperature higher than the deuterium boiling point. Buscaglia et al.’s initial simulation encompassed the moderator cell only. Subsequent simulations of the OPAL CNS by Pavlou et al. [4] using ANSYS-CFX included the complete thermosiphon (Figure 1), with both helium and deuterium flows. This removed the limiting assumptions of inlet velocity and inlet temperature for the moderator-only simulation.
In this paper, the same CFD mesh model (Figure 2) of the complete thermosiphon system was recalculated using the commercial CFD code STARCCM+ V11.02. These results were evaluated against the previous ANSYS CFX V14.5 results attained by Pavlou et al. The first objective of this study was to cross-verify two CFD computer codes commonly used in the nuclear industry for the analysis of the CNS system. The second objective was to use the validated STARCCM+ model to examine two separate shutdown conditions to determine the transient peak temperatures that would be experienced within the entire thermosiphon system. The results will demonstrate the capability of CFD for detailed nuclear engineering and design by showing the heat transfer dynamics of the thermosiphon and its internal flow structure.
2. Mathematical Equations for the Thermosiphon System
The recirculation of deuterium in the thermosiphon loop is driven by the action of natural convection established under the influence of nuclear heating. This nuclear heat creates a thermal gradient across the vertical axis of the moderator cell and the hot pipe (return leg of the thermosiphon; see Figure 1). This results in density stratification in the deuterium which creates an upward flow under the force of buoyancy. The deuterium rises in the hot pipe and cools in the heat exchanger at the top of the thermosiphon loop before draining into the cold leg and returning to the moderator cell. During normal operating condition, the two helium legs operate under forced convection with one leg of helium circulating through the heat exchanger and the other through the moderator cell jacket. If helium flow stops, the reactor automatically shuts down. In such a state, reactor decay heat still drives a significant amount of conductive and natural convective heat transfer in the thermosiphon system, whether the deuterium is in liquid or gaseous state. Thus, based on these characteristics, the effective conservation equations governing mass, momentum, and energy in Favre-Averaged Navier-Stokes (FANS) form must be solved for assessing the fluid flow and heat transfer processes. They can be written aswhere is the mean density, is the velocity vector, is the mean pressure, is the gravitational vector, is the turbulent kinetic energy, is the specific enthalpy, is the scalar solid temperature, is the dynamic viscosity of the fluid, is the thermal conductivity of the fluid, is the mean specific heat of constant pressure of the solid material, is the thermal conductivity of the solid material, and and are the volumetric heating source of the fluid and solid material, respectively. The turbulent or eddy viscosity was calculated according towhere is the dissipation of turbulent kinetic energy of the fluid. The turbulent Prandtl number was set to a value of 0.9. Appropriate thermophysical properties such as dynamic viscosity, thermal conductivity, and mean specific heat of constant pressure of deuterium and helium were considered [17–20]. Material properties of the aluminium alloy moderator cell and jacket and pipes and the stainless heat exchanger were extracted from NIST data [20].
In this study, the Shear Stress Transport (SST) turbulence model by Menter [21] was adopted. The SST turbulence model entails a combination of favoured features of the model with the model so that the inner region of the boundary layer is adequately resolved by the latter while the former is employed to obtain numerical solutions in the outer part of the boundary layer. From previous study by Pavlou et al. (2016), assessment of different turbulence models has revealed that the SST turbulence model performed better than the standard turbulence model [22] with logarithmic wall functions to bridge the bulk fluid flow with the boundary layer flow near the walls of the thermosiphon system. The favourable feature of the model in resolving the wall-bounded flow allows the CFD wall-model to handle nonequilibrium boundary layer regions.
3. Numerical Details
A one-to-one assessment between the CFD models of ANSYS CFX and STARCCM+ was performed. This entailed exporting the mesh structure used in the ANSYS CFX V14.5 simulation into the STARCCM+ V11.02 simulation. Mesh convergence for the common CFD model used was established by testing three meshes with nominal steady-state operating parameters. A coarse mesh of 4 million elements, a medium mesh of 11 million elements, and a fine mesh of 23 million elements were tested using the CFX software. In the results, a difference of 7.7% in maximum temperature was found between the coarse mesh and the medium mesh and a difference of 3.6% was found between the medium mesh and fine mesh. To strike a balance between computing solution time and accuracy, the medium resolution mesh was adopted for both the CFX and StarCCM+ simulations. Further details of the CFD model are included in Pavlou et al. [4]. In the CFD model, the finite thickness of the baffle structure within the heat exchanger was not considered; appropriate boundary conditions were alternatively imposed whereby no fluid was allowed to flow across the baffles but conjugate heat transfer was captured in these thin solid domains by way of 1D conduction modelling. The mesh of the entire thermosiphon system was approximately 11 million elements. Except for the recirculation of deuterium in the thermosiphon loop which was driven solely by buoyancy, the same boundary conditions were applied to the helium loop which included the specifications of the mass flow rates, inlet temperatures, and the outlet relative pressures. Because of the identical mesh adopted in both CFD models, the assessment was focused on the differences in numerical discretisation and models implemented within the two computer codes.
Finite volume method was employed to discretise the effective transport equations governing the fluid flow and heat transfer processes in ANSYS CFX and STARCCM+. For the fluid region, the advection terms in the equations governing the conservation of mass, momentum, and energy and turbulence transport of turbulent kinetic energy and dissipation of turbulent kinetic energy were approximated by the high resolution differencing scheme in ANSYS CFX while they were approximated by the second-order upwinding differencing scheme in STARCCM+. These difference schemes were in keeping with the recommended settings of each code. The diffusion terms in both computer codes were approximated via the second-order central differencing scheme. Linkage between the pressure and velocity was achieved through a coupled solver strategy to enhance the convergence of the pressure.
The time derivative terms in all of the transport equations were approximated via a second-order differencing scheme. An implicit procedure was adopted to obtain the transient solution of the fluid flow and heat transfer. For the verification and validation of CFD codes against prior experimental tests, the solutions were marched to steady state with a target normalised residual level of for all variables. During the simulation of the thermosiphon system, the heat balance was monitored in addition to the residuals to ensure the power input equals the power output in the pseudo steady-state simulation. For reactor shutdown conditions, the solution was marched with a varying time stepping ranging from 0.06 to 0.6 seconds.
4. Results and Discussion
4.1. Verification and Validation of CFD Codes against PNPI Tests
A full-scale prototype testing of the CNS thermosiphon system was performed at Petersburg Nuclear Physics Institute (PNPI) [23]. The capability of the thermosiphon system in removing the heat from the moderator cell subject to varying coolant flow rates and varying heat loads was demonstrated. For the purpose of verification and validation of the computational results, numerical predictions were assessed against the measured heat removal capacity of the cryogenic coolant system under estimated operational heat loads.
Figure 3 shows the schematic diagram of the experimental rig. This figure also shows the helium and deuterium buffer tanks (coloured in brown and green) besides the CNS thermosiphon system. In this assessment, the buffer tanks have not been modelled in the CFD simulations. Relevant parameters for the experimental rig depicted in Figure 3 are summarised in Table 1.
For the closed circulation loop of the deuterium, no boundary conditions are required or imposed. For the flow of helium, experimentally measured mass flow rates of 80 g/s were prescribed at the inlet locations of HeF1 and HeF2 and a temperature of 19 K was set at both inlet locations at HeT1 and HeT2 of the thermosiphon system. The power was raised in incremental steps, up to a total value of 5000 W during the experiments. At each incremental power step, the power level was held steady for a period of 15 minutes. For this study, pseudo steady-state calculations were performed for the validation of two independent CFD codes against the PNPI experimental data which operated at a nominal total power of approximately 4000 W. Comparison between relevant parameters of PNPI’s experimental testing and predictions of STARCCM+ and CFX at a nominal power of 4000 W for the SST turbulence model is illustrated in Table 2. Between the computational predictions and experimental measurements, it can be seen that there was a remarkably close agreement for all assessment criteria. The largest percentage difference for both CFD models was recorded to be around 3%. In the evaluation of the thermosiphon system by Pavlou et al. [4], the SST turbulence model of CFX has been shown to perform better than the turbulence model.
The use of SST turbulence model for STARCCM+ has also been shown to yield comparable results. One possible explanation for the close agreement with experimental measurements is the accuracy of the SST model for capturing the flow dynamics of near-wall regions which are a dominant flow feature around the tube bundles of heat exchangers. Since the deuterium flow is in the range of low Reynolds number as it travels through the tubes of the heat exchanger, the SST turbulence model, which uses the turbulence model in the viscous sublayer, performs well when boundary layer flows are involved [9]. In both CFD models, the simulations using the SST turbulence model show more heat is removed via the heat exchanger and hence a lower average temperature is predicted for the deuterium flow.
Within the CNS thermosiphon system, the closed circulation loop is governed entirely by the physics of buoyancy forces as determined by the density difference in the deuterium. As can be seen in Figure 4, both CFD models predict the fluid rising along the “hot leg” before passing through the top of the heat exchanger and draining through the “cold leg” of the system. A closer examination on the predicted temperatures reveals that STARCCM+ yields marginally lower temperatures than CFX. This is primarily due to the increased heat removal within the moderator chamber by a higher helium flow velocity as predicted by STARCCM+ compared to the CFX result as seen in Figure 5(a). The discrepancy can be explained by the discretisation schemes being adopted in the respective computer codes. STARCCM+ utilises a second-order upwind differencing scheme while CFX uses a blending of first order upwind and central differencing schemes for the advection term. Thus, the fluid flow could be anticipated to be more diffused through the use of the first order upwind differencing scheme in CFX; and hence a lower velocity would be predicted for the helium flow. As a result of the higher helium flow velocity running through the moderator of the STARCCM+ simulation, the temperature of the deuterium entering the hot leg and subsequently in the heat exchanger was comparatively lower (see Figure 6(a)). In Figure 6(b), both CFD codes were able to indicate the heat exchanger works well and without any excessive recirculating regions appearing between each baffle. The velocity streamlines pass smoothly from baffle to baffle and the cyclic behaviour of the shell-side flow can be seen across the whole heat exchanger. In the heat exchanger, the helium velocity peaks at a value of approximately 7.25 m/s for the STARCCM+ simulation and around 8 m/s for the CFX simulation. However, the STARCCM+ simulations show a wider streamline distribution which reduces the regions of recirculation and predicts better flow mixing characteristics in the shell-and-tube structure.
(a)
(b)
(a)
(b)
4.2. Determining the Peak Temperatures of the Thermosiphon System during Reactor Shutdown Conditions
To demonstrate the versatility of CFD tools for nuclear subsystem design, the CFD model in STARCCM+ was used to determine the transient heat transfer characteristics after the reactor and cryogenic system is shutdown. The OPAL reactor was designed to operate at full power () even when the helium cryogenic system is not operating (i.e., in Standby Operation or SO mode). In SO mode, the nuclear heat is removed by forced convective helium circulating at a room temperature. If the forced convective helium flow is stopped, the reactor will trip to ensure the moderator chamber does not overheat. After the reactor is shutdown, a low level of decay heat persists as a result of fission product decay gamma emitted from the core into the moderator chamber. The action of this decay heat plus the absence of forced convective helium cooling means that the natural circulation of gaseous deuterium assume the task of transferring heat away from the moderator cell and into the thermosiphon heat exchanger. It should be noted that the CFD model used in this transient simulation had been validated for single liquid phase flows during normal operation as shown in Section 4.1 of this paper. For this transient simulation, an initial temperature of 60°C is adopted for the thermosiphon system. This initial temperature of 60°C applies to the heat exchanger, the moderator cell, the connecting pipes, and all gaseous domains, namely, the helium and deuterium volumes. Under the assumption of a completely gaseous fluid domain, the transient simulation model remains validated, given the same single phase conservation equations that were solved previously for a liquid system are solved for the gaseous system. The only requirement was in ensuring that the appropriate physical properties for gaseous deuterium are accounted for.
The OPAL reactor has two independent shutdown systems to provide redundancy and thus two trip scenarios with different power transients are analysed. Table 3 shows the transient decay heat inputs for the moderator chamber after the initiation of the first shutdown system (Trip 1) and second shutdown system (Trip 2) of the OPAL reactor.
Imposing these transient heat sources results in a maximum aluminium moderator cell temperature of 120°C in the case of Trip 1 and 132°C in the case of Trip 2 (see Figure 7). Both Trip 1 and Trip 2 maximum temperatures occur roughly 300 seconds after the cessation of helium cooling and reactor shutdown. As time progresses, the temperature drops, indicating heat is transferred from the moderator cell to the rest of the thermosiphon system. Most of the heat flows to the heat exchanger which serves as the local heat sink. The heat exchanger is connected by numerous piping to the ultimate heat sink which is the process system(s) outside the reactor block. Figure 8 shows the temperature contours of the thermosiphon system predicted for Trips 1 and 2 at 10 min.
(a)
(b)
(a)
(b)
The importance of natural convection heat transfer is shown in Figure 9 which shows a plot of three simulation scenarios:(a)A scenario with a total absence of heat transfer processes (i.e., excluding all conduction, convection, and radiation) resulting in a continual rise in moderator-cell temperature. The system is adiabatic in that no heat is allowed to escape the thermosiphon system.(b)A conduction-only modelling scenario whereby heat is allowed to conduct within the thermosiphon system itself but is not allowed to escape from the thermosiphon.(c)A conduction and natural convection modelling scenario whereby heat is allowed to migrate from its concentrated source in the moderator cell to the heat exchanger acting as the heat sink: again, heat is not allowed to migrate out of the thermosiphon system.
All three scenarios are conservatively bound by the notion that despite the presence of a heat source in the moderator cell, heat is not allowed to migrate out of the system. Of the three simulations shown in Figure 9, the first scenario with an absence-of-heat-transfer gives an unphysically conservative result, resulting in the maximum temperature increase of 145 K at 2000 seconds. Given this trend, a temperature of 200°C could be attained, the point at which aluminium begins to soften and reduces in strength. The second conduction-only scenario fairs slightly better, resulting in a maximum temperature increase of 105 K at 2000 seconds. Overall, the conduction-only model increases with the same temperature trend as the first scenario with no heat transfer. The third scenario with the most realistic representation including both conduction and convection predicts only a maximum temperature rise of 60 K at 250 s and subsequently plateaus to a rise of 50 K. Given an initial model temperature of 60°C the temperature of the thermosiphon system will not exceed 120°C, thus remaining within the safe operating limits of the CNS material and structural design. As shown in Figure 9, it is readily apparent that the heat transfer process is markedly improved with the inclusion of natural convection cooling in the numerical modelling and therefore was found to be an indispensable component of the heat transfer modelling.
5. Conclusions
Behind the day to day operations of the OPAL cold neutron source much work is done to ensure its continued safe operation and operational readiness. Here the question of the OPAL CNS temperature was posited in the event of a halt in helium coolant flow. To answer this question a thorough CFD study was carried out. Two commercial CFD codes using an identical volumetric mesh were used to simulate the pseudo steady state operation of the OPAL CNS thermosiphon cooling. The results were validated against experimental results borne from the CNS prototype test. Both CFD simulations, using a combination of expansionary and unstructured mesh and by using the Shear Stress Transport (SST) turbulence model, attained good results that were within a difference of 3% for helium and deuterium temperatures compared against the experiments. Small difference in temperatures attained by ANSYS-CFX and Star-CCM+ could be attributed to a different selection in differencing scheme for solving the set of conservation and turbulence transport equations. However, these small differences turned out to be small and within the bounds of experimental uncertainty.
After the STARCCM+ steady state simulation was cross validated against experimental data, the CFD model was used to conduct transient studies of thermosiphon temperature increases under the influence of reactor decay heat and the absence of helium convective cooling. For the analysis of transient conditions with an initial temperature of 60°C, the maximum temperature attained by the thermosiphon in a Trip 1 scenario was 120°C and 132°C for a Trip 2 scenario. These high resolution results in time and space show the CNS thermosiphon system will remain safe after the reactor is shut down and without active helium cooling. This study also demonstrates the versatility of CFD, as a preventative maintenance tool for investigating the demanding heat transfer characteristics of a nuclear engineering system.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.