Thermal-hydraulic modelling of the Cold Neutron Source thermosiphon system

Annals of Nuclear Energy 90 (2016) 135–147

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Thermal-hydraulic modelling of the Cold Neutron Source thermosiphon system Tom Pavlou a, Mark Ho b, Guan Heng Yeoh a,b,⇑, Weijian Lu b a b

School of Mechanical and Manufacturing Engineering, University of New South Wales, Sydney, NSW 2052, Australia Australian Nuclear Science and Technology Organisation (ANSTO), Locked Bag 2001, Kirrawee DC, NSW 2232, Australia

a r t i c l e

i n f o

Article history: Received 17 March 2015 Received in revised form 2 September 2015 Accepted 21 November 2015 Available online 23 December 2015 Keywords: Cold Neutron Source CFD Cryogenic Thermosiphon system

a b s t r a c t In the OPAL research reactor at ANSTO, the Cold Neutron Source (CNS) generates cold neutrons through the interaction of the thermal neutrons and the liquide deuterium inside a moderator cell for neutron scattering experiments. In this paper, a computational fluid dynamics (CFD) model is developed to provide a better understanding of the fluid flow and heat transfer characteristics of the CNS thermosiphon system. A complete model for the system is constructed which includes the closed circulation loop of the flow of deuterium driven solely by buoyancy through the moderator cell, heat exchanger and two interconnecting pipes. The loop for the helium coolant which is required to remove the heat deposited in the CNS thermosiphon system is also included in the model. Based on assessment of the full-scale model against experimental data, good agreement is achieved between measurements and predictions. The model predictions are also consistent with theoretically calculated values for the mass flow rates of the circulation flow of deuterium and the total heat removal of the entire CNS thermosiphon system. This validated CFD model will allow for the assessment of the behaviour of system under different operating conditions and transients due to inaccessibility by process instrumentation within the reactor environment. It also provides us with a design tool for future modifications without the need for costly prototype testing. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The Australian Nuclear Science and Technology Organisation’s (ANSTO) Cold Neutron Source (CNS) (Thiering and Lu, 2006), which was commissioned in 2006, is located in the heavy water reflector of the 20 MW Open-Pool Australian Light-water (OPAL) research reactor as shown in Fig. 1. This CNS generates cold neutrons through the interaction of the thermal neutrons from the reactor core and the cryogenic deuterium inside the moderator cell for neutron scattering experiments. Neutron scattering is a highly useful tool for studying all types of condensed matter and the CNS in OPAL supports numerous instruments (Kennedy, 2006) including: small angle neutron scattering, cold quasi-Laue diffractometer and a polarisation analysis spectrometer (see Fig. 1). Much advancements made in the field of neutron scattering have been beneficial to the aerospace, automotive, healthcare, electronics and even food industries (Gil, 2011; Lopez-Rubio and Gilbert, 2009). ⇑ Corresponding author at: School of Mechanical and Manufacturing Engineering, University of New South Wales, Sydney, NSW 2052, Australia. Tel.: +61 2 9385 4099; fax: +61 2 9663 1222. E-mail address: [email protected] (G.H. Yeoh). http://dx.doi.org/10.1016/j.anucene.2015.11.034 0306-4549/Ó 2015 Elsevier Ltd. All rights reserved.

Besides the CNS in OPAL, there are a number of other CNSs at various reactors around the world such as in Germany (Axmann et al., 1997), France (Ageron, 1980), Korea (Kim et al., 2008), Hungary (Grósz et al., 2000) and the United States of America (Alsmiller and Lillie, 1992). Among all these, the CNSs at the FRM-II in Germany and ILL in France use liquide deuterium as the moderator, similar to the OPAL CNS. For the FRM-II CNS, Gaubatz and Gobrecht (2000) indicates that the total heat load is around 4 kW which is comparable to the heat load of the OPAL CNS. ILL in Grenoble is home to the world’s only other operational liquide deuterium CNS. It has both a Vertical Cold Source (VCS) and a Horizontal Cold Source (HCS) (Ageron, 1980). Having been constructed in the 1980s, the CNS at ILL was seen as a benchmark for both CNSs in FRM-II and OPAL. From a thermo-hydraulic standpoint, the CNS thermo-siphon system in OPAL consists of one closed deuterium loop, which operates under natural convection, two helium coolant flows, which are supplied by a Brayton-cycle helium refrigerator, a stainless steel heat exchanger and an aluminium moderator. The moderator cell is where the neutrons from the reactor core interact with the cryogenic deuterium moderator. Part of the neutron population would come into thermal equilibrium with the moderator. It is a design criterion of the OPAL CNS that

136

T. Pavlou et al. / Annals of Nuclear Energy 90 (2016) 135–147

Nomenclature Cp G H k p PrT Q q00 T Tref U

specific heat of constant pressure generation of kinetic energy due to buoyancy specific enthalpy turbulent kinetic energy pressure turbulent Prandtl number volumetric heating source Reynolds flux term temperature reference temperature velocity vector

lT s00

q qref e

k

turbulent or eddy viscosity Reynolds stress flux term density reference density dissipation of turbulent kinetic energy thermal conductivity

Subscripts g gas l liquid s solid

Greek symbols l dynamic viscosity

the deuterium remains in liquid phase that there can be no boiling, whether bulk or localised, as this would decrease the ability of the deuterium to moderate the neutrons. There have been a number of publications in literature outlining the numerical simulations of different components of the CNS thermosiphon system in isolation but not the modelling of the entire system. Based on a CFD model, Ozden and Tari (2010) presented findings of the effects of baffle spacing and orientation on the turbulent shell-side fluid flow of a heat exchanger. They stated that the effects of bypass and leakage flows on the overall performance of a heat exchanger are significant but for a small heat exchanger, such as the one modelled, these effects are minimal and can be neglected in comparison to the main flow stream. As their analysis was only concerned with the shell-side configuration, the tubeside heat exchanger was not modelled. Three different turbulence models were tested: standard k
e, realisable k
e and the SpalartAllmaras. The CFD results gave reasonable agreement with those predicted by the analytical Bell-Delaware method. The total heat transfer rate varied from the analytical calculations by less than 10% whereas the pressure drop varied by up to 34% for the configuration using six baffles. Note that the CFD results were in closer agreement with the analytical predictions from the BellDelaware model when more baffles were used. Mohammadi

Fig. 1. Plan view of OPAL’s reflector vessel: (1) reactor core, (2) CNS thermosiphon system, and (3) hot neutron source flange and labelled neutron beam tubes (TG1, TG3, TG4, CG1, CG3, CG4, HB1 and HB2).

et al. (2009) presented the CFD modelling of a shell and tube heat exchanger with the effects of leakage and bypass flows incorporated into a model of 76 tubes. RNG k
e turbulence model due to its improved accuracy for rapidly strained flows was utilised. The model was validated against the Verein Deutscher Ingenieure (VDI) method showing close agreement in terms of the shell-side heat transfer coefficient as well as the pressure drop using a geometrical arrangement with horizontally orientated baffles. They concluded that the bypass flows and leakage effects played an important role in the performance of such a heat exchanger. Wang et al. (2009) proposed a new type of heat exchanger that has a continuous helical baffle. Although the geometry of the heat exchanger being modelled in this case was different to the heat exchanger of the CNS in OPAL, the numerical approach and methods used to calculate the overall heat transfer coefficient and pressure drop are comparable to the studies utilising the CFD approach. Realisable k
e turbulence model was adopted and they have ascertained that their novel design with a helical baffle outperformed the conventional baffle design in terms of increasing the rate of heat transfer and decreasing the pressure drop across the shellside configuration. There are several key points to observe on CFD modelling of heat exchangers. Although leakage and bypass flows play a role in the characteristics of a heat exchanger, they can be ignored without introducing significant errors. CFD can indeed be employed to model a heat exchanger to the accuracy of obtaining a usable engineering solution and that the k
e turbulence model (in one of its forms) appears to present the best agreement with existing analytical methods for calculating the pressure drop and rate of heat transfer within a heat exchanger for the shell-side configuration. For the moderator cell, Buscaglia et al. (2014) presented a Finite Element Analysis (FEA) of the moderator cell of a CNS. This paper is of particular interest as it signifies the only prior numerical modelling of any kind for the moderator cell being modelled in this current study. A transient analysis was performed to capture the turbulent flow and identify areas of flow recirculation which could lead to localised heating and possible boiling of the deuterium. The numerical simulations showed the incoming jet of deuterium hitting the bottom of the cell formed a fast moving thin sheet that followed the contour of the bottom of the cell. Although no areas of recirculation were observed, there was significant mixing of the deuterium in the moderator cell. The CFD results showed that the average temperature throughout the cell did not exceed the inlet temperature by more than 4.5 K. It was also observed that the maximum bulk fluid temperature, even near the walls, was approximately 27 K (7.5 K greater than the inlet temperature) and hence

T. Pavlou et al. / Annals of Nuclear Energy 90 (2016) 135–147

well within the required operating limits. This model was effective in proving that under the conditions specified there should not be any localised boiling of the deuterium and that the helium flow is effective in removing the required heat from the deuterium. However, the biggest assumption made was the inlet temperature and velocity of the deuterium. This analysis was only performed with one inlet velocity and it is possible that the behaviour of the flow and local temperatures could be very different given different inlet conditions. The model being developed in this project aims to address this uncertainty by removing the assumption of an inlet velocity and temperature and modelling the entire deuterium loop. Also, the model presents a number of significant simplifications: the helium coolant was not modelled, instead an empirical equation was used to model the heat transfer between the deuterium and the wall was developed. An inlet velocity for the deuterium was specified, whereas in reality, the closed deuterium loop operates under natural convection and hence there is no specified inlet condition. The nuclear heat load on both the deuterium fluid and aluminium moderator cell were functions of distance from the reactor core derived from neutronics calculations. Sun et al. (2012) also performed an FEA study of the aluminium structure of a moderator cell. Unlike the model considered in Axmann et al. (1997), both the hydrogen moderator and helium coolant flows were modelled using CFD. The moderator cell being modelled in this case is located in the China Advanced Research Reactor (CARR) and is geometrically very similar to the configuration in OPAL. In fact the entire in-pile structure is similar consisting of a closed loop (with hydrogen rather than deuterium) operating under natural convection and two helium coolant flows: one through the heat exchanger and one through the moderator cell itself. The turbulent flow within the CNS was resolved using the standard k
e model. The results showed that the maximum temperature in the hydrogen moderator was 21.9 K which meant that boiling did not occur. It should be noted that the paper mainly focussed on the mechanical stresses developed in the aluminium structure rather than the flow characteristics and heat transfer which are of more concern in the current study. To the best of our knowledge, there does not appear to be any literature outlining thermo-hydraulic testing or CFD modelling of an entire CNS thermosiphon system. There is certainly considerably value to have a better understanding of the flow characteristics of the deuterium circulation loop, to be able to assess its behaviour under different operating conditions due to the low cryogenic temperatures (approximately 20 K) and inaccessibility of the system of which useful instrumentation could be placed within to provide specific information on the flow rate, temperature or any other flow characteristics. The OPAL CNS has experienced several serious process failures which have been very costly in terms of repairs and down-time (Lu and Thiering, 2012). Having a full thermo-hydraulic CFD model of the CNS thermosiphon system would certainly assist in better understanding its behaviour and improve its reliability and reduce down-time and repair costs. In addition to this, having such a model of the current CNS thermosiphon system would be a useful starting point for the design of the next generation CNS in the coming years at ANSTO. 2. Description of system physics The CNS thermosiphon system is a complex thermo-hydraulic system that comprises a number of subsystems as illustrated in Fig. 2. The three main subsystems are the moderator cell, the heat exchanger and the pipes connecting the moderator chamber and heat exchanger. Each of these subsystems is complicated in its own right and exhibits intricate flow structures and heat transfer.

137

From a neutronics perspective, the moderator cell represents the most important component of the in-pile system. The main purpose of a CNS is to produce cold neutrons, which are produced in the moderator cell as thermal neutrons from the heavy water region surrounding the CNS interact with the cryogenic deuterium. During these interactions, the neutrons lose energy and slow down into the cold spectrum. After leaving the moderator cell and travelling down guide tubes they are then used in neutron scattering experiments. The moderator cell is approximately 400 mm in height and 300 mm in diameter and is made from aluminium alloy. It is cylindrical in shape with toro-spherical shell caps on the top and bottom. The liquide deuterium enters the moderator cell downwards through the inlet tube and exits upwards through the outlet which is located on the vertical axis about which the moderator cell is rotationally symmetric. The helium jacket is very similar in shape to the moderator cell but is slightly larger and is also made from aluminium alloy. There exists a gap of approximately 2 mm between the helium jacket and the moderator cell. The helium jacket contains the cryogenic helium that is used to cool the CNS. The helium enters the moderator chamber through an annular region at the top of the inlet pipe and travels downwards through the displacer before flowing upwards through the 2 mm annular gap between the jacket and the moderator cell. It exits the moderator chamber through an annular region surrounding the outlet of the deuterium. The displacer is a design feature, commonly known as a reentrant cavity that is aimed to optimise the yield of cold neutrons. The helium coolant, which flows through the displacer, has very low scattering and absorption cross sections meaning that it is practically ‘‘invisible” to the neutrons. From a thermo-hydraulic standpoint, the displacer does not serve any particular purpose; however, it is critical to the fluid flow and heat transfer within the moderator chamber and is thus included in the CFD model. The function of the stainless steel cryogenic heat exchanger is to cool the natural convection deuterium loop. It is a counter-flow, shell and tube heat exchanger, with the warm deuterium entering at the top and flowing downwards and the helium coolant entering at the bottom and flowing upwards. It should be noted that the helium flowing through the heat exchanger is separate to the helium flowing through the moderator chamber. The total heat removal of the heat exchanger under normal operation is estimated to be 1700 W. There are two main pipes that connect the moderator cell to the heat exchanger: one for the downwards arm (cold leg) and one for the upwards arm (hot leg) as seen in Fig. 1. These pipes are also made from the aluminium alloy and are essentially parallel-flow tube-in-tube heat exchangers. The deuterium fluid flows through the circular pipe and the helium coolant flows through the annular cross section surrounding the deuterium. The flow of the liquide deuterium within the CNS thermosiphon system takes place due to natural convection, driven by density differential in the fluid due to temperature gradients. The natural convection loop is established as the liquide deuterium warms up under nuclear heat load in the moderator chamber at the bottom of the loop and cools down in the heat exchanger at the top of loop, and the hot and cold legs provide the flow paths. The most important consideration in the CFD model was applying the correct heat load to the liquide deuterium which dictates the characteristics of the natural convection loop. Both of the helium coolant flows operate under forced convection. They are supplied by a 5 kW helium refrigerator operating on Brayton Cycle (Lu and Thiering, 2012). The amount of heat removed by the helium flows is dependent on the flow characteristics of the helium. It is for this reason that every effort is made during the geometry and mesh creation stages to ensure that the model was as accurate as possible and that it truly represents the physical system in the CNS.

138

T. Pavlou et al. / Annals of Nuclear Energy 90 (2016) 135–147

Fig. 2. Schematic diagram of the CNS thermosiphon system.

While there is significant conductive heat transfer within the CNS thermosiphon system, there is very little conduction from the outside to the CNS. The thermosiphon is thermally insulated in a vacuum containment during cryogenic operation. Within the structure heat is conducted through the aluminium moderator cell and pipes as well as the stainless steel heat exchanger. The material properties of deuterium are perhaps the most important modelling consideration since at cryogenic temperatures, a small change in the temperature can result in a big change in density or specific heat capacity at constant pressure. Appropriate data for the thermo-physical properties of deuterium and helium can be found in Lemmon et al. (2005), Richardson and Leachman (2012), Marquardt et al. (2002), Misiorek et al. (1981). Material properties of the aluminium alloy moderator cell and jacket and pipes, and the stainless steel heat exchanger are extracted from NIST (Misiorek et al., 1981). It should be noted that the thermophysical properties of helium, aluminium alloy and stainless steel are rather well known and thus do not require further quantification on the use of these properties during the numerical simulation. For deuterium, it was however ascertained that the average

deviations of calculated values from experimental data for the saturation density, heat conductivity and viscosity are estimated as ±0.2%, ±1.5% and ±1% respectively.

3. Mathematical formulation In order to resolve the fluid flow and heat transfer of the cryogenic deuterium and helium coolant within the CNS thermosiphon system, the Favre-Averaged Navier-Stokes (FANS) equations governing the conservation of mass, momentum and energy for helium gas (g) and liquide deuterium (l) are solved:

   g;l @q e g;l ¼ 0  g;l U þr q @t   e g;l U @ q @t

      T  e g;l ¼
rp e g;l  U e g;l þ r U e g;l U g;l þ r  lg;l r U þr q    2 e g;l I
r  s  g;l
qref 00g;l þ q
lg;l r  U g;l g 3

ð1Þ

ð2Þ

T. Pavlou et al. / Annals of Nuclear Energy 90 (2016) 135–147



e g;l  g;l H @ q



@t

    e g;l H e g;l ¼ r  kg;l r Te g;l
r  q  00g;l þ Q g;l  g;l U þr q ð3Þ

e is the velocity  and qref the mean and reference densities, U where q  is the mean pressure, g is the gravitational vector, l is the vector, p dynamic viscosity, k is the thermal conductivity and Q is the volumetric heating source term. The Reynolds stress flux and Reynolds  00 in Eqs. (2) and (3) that need to be resolved 00 and q flux terms s via the consideration of a two-equation turbulence model can be written as

  T  2 e g;l I
2 q e g;l þ r U e g;l  g;l kg;l I 00g;l ¼ lTg;l r U
lTg;l r  U
s 3 3  00 ¼
q

lTg;l Pr Tg;l

e g;l rH

ð4Þ

ð5Þ

In Eq. (5), the turbulent Prandtl numbers for helium gas and liqeg uide deuterium are set to a value of 0.7. The specific enthalpy H e l for liquide deuterium are defined as for helium gas and H

eg ¼ H

Z eT g T ref g

el ¼ C p;g dT or H

Z eT l T ref l

C p;l dT

ð6Þ

where C p is the mean specific heat of constant pressure dependent e and T ref is the reference temperature. on the mean temperature T Conjugate heat transfer across the solid (s) regions of the CNS thermosiphon system is accounted for of which the transport equation of the temperature can be written as:

   C p;s Te s @ q @t

  ¼ r  ks r Te s þ Q s

2

kg;l

eg;l

ð8Þ

where Cl is a constant with a value set at 0.09. Also, the effect of buoyancy is accounted for the recirculating deuterium flow within the closed loop through the consideration of a source term in the turbulence models via

G¼


lTg;l  g;l g  rq q g;l

meshed as part of a larger mesh, and (2) it enables certain sections of the mesh to be improved and refined without having to re-mesh the entire geometry. This has proven to be useful as initial CFD testing have shown areas of large gradients, particularly in the pipes leading into and out of the displacer of the moderator cell. Individual sections are re-meshed to better capture the flow effects in these regions in the computational geometries before producing the final results. Fig. 3 shows the cross-sectional view of the structured mesh of hexahedral elements generated for the aluminium moderator solid. Owing to the fact that many of the solids being modelled are of thin regions (1–10 mm), care has been taken to ensure that multiple mesh elements are present across these thicknesses. This ensures that the temperature is adequately predicted from one side of the surface to the other. However, the fluid regions within the moderator cell have been filled with unstructured mesh of tetrahedral elements (see Fig. 3), which have been created using Hypermesh 12.0, CFD Tetramesh generator. It should be reiterated that the heat exchanger in the CNS thermosiphon system has 12 baffles and 133 tubes. A novel solution is developed to create the effect of baffles in the helium coolant flow without actually meshing the solid domain. This involves meshing cylindrical slices corresponding to the height of the baffle spacing of the heat exchanger. These slices are stacked with appropriate boundary conditions – no fluid is allowed to flow across the baffles but conjugate heat transfer is present through these thin solid domains. An example of a slice can be seen in Fig. 4 where the yellow1 region represents a baffle. A final assembly of all meshes generated via Hypermesh 12.0 is created to position the meshes of the different slices correctly relative to one another. Generated meshes of the top section of the heat exchanger with 133 tubes for the flow of deuterium are also shown in Fig. 4.

ð7Þ

Two-equation turbulence models are adopted to handle the turbulent fluid flow within the CNS thermosiphon system. In this study, the standard k
e (Launder and Spalding, 1974) and Shear Stress Transport (SST) turbulence models are assessed (Menter, 1993). In both models, the turbulent or eddy viscosity is determined according to

lTg;l ¼ C l q g;l

139

ð9Þ

4. Computational details 4.1. Mesh generation Computational geometries of the subsystems of the CNS thermosiphon system have been created using SolidWorks 2013. Following the creation of these geometries, Hypermesh 12.0 is employed to generate the computational meshes. It should be stressed that the generation of the computational meshes constitutes the most arduous process in the development of the computational geometries. Owing to the size of the entire CNS thermosiphon system, the meshing process is broken down into smaller, more manageable components. This serves two purposes: (1) it enables components to be meshed that could not have been

4.2. Assembly of subsystems and boundary conditions Owing to the fact that the mesh has been created in smaller manageable components, a significant amount of patching is required to bring together the different subsystems by the logical connection of interfaces. This allows mass and heat transfer for the fluid regions and heat transfer to be communicated across the liquid and solid regions. The computational code of ANSYS CFX V14.5 is adopted to assemble the three fluid domains and one solid domain of the CNS thermosiphon system. Subdomains are also created to allow for the specification of certain boundary conditions that are specific to certain regions of the model such as the nuclear heat load on the deuterium fluid in the moderator cell. A grid sensitivity analysis is performed to assess three computational meshes for the entire CNS thermosiphon system – approximately 4 million elements (coarse mesh), 11 million elements (medium mesh) and 23 million elements (fine mesh). The difference for the simulation results between the coarse and medium meshes is found to be 7.7% and between the medium and fine meshes to be 3.6% for the predicted maximum temperatures. As a compromise between computational accuracy and efficiency, the medium mesh (11 million elements) is chosen to determine the rest of the numerical results in this work. Except for the flow of recirculating deuterium within the closed loop that is driven solely by buoyancy, appropriate boundary conditions are applied to the helium loop which includes the specifications of the mass flow rates, inlet temperatures and the outlet relative pressures. 1 For interpretation of colour in Figs. 4–6 and 11, the reader is referred to the web version of this article.

140

T. Pavlou et al. / Annals of Nuclear Energy 90 (2016) 135–147

Deuterium Deuterium Inlet Outlet Helium Inlet

simulation of the CNS thermosiphon system, the heat balance is monitored in addition to the residuals so that the power input equals to the power output at pseudo steady-state operations. With the computational model solving at approximately 4 min/time step on 64 parallel processors on the Trentino High Performance Cluster at University New South Wales, a typical solution takes about 3 weeks for this balance to be reached.

Helium Outlet

5. Results and discussion 5.1. Validation against PNPI tests

Displacer Boom connecon pipe

Fig. 3. Cross-sectional view of moderator cell, showing the displacer, alongside with the generated mesh comprising hexahedral and tetrahedral elements.

4.3. Numerical methods Finite volume method is employed to discretise the transport equations governing the fluid flow and heat transfer. 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 are approximated by the high resolution differencing scheme. For all regions, the diffusion terms are approximated via the second order central differencing scheme. Linkage between the pressure and velocity is achieved through a coupled solver strategy to enhance the convergence of the pressure. The time derivative terms in all of the transport equations are approximated via a second order differencing scheme. An implicit procedure is adopted to obtain the transient solution of the fluid flow and heat transfer with a time step of 0.025 s. For the iterative procedure at each time step, target normalised residuals are set to 10
5 for mass, momentum and turbulence transport equations while a target normalised residual is set to 10
6 for the enthalpy and temperature transport equations in the fluid and solid regions respectively. During the

Petersburg Nuclear Physics Institute (PNPI) conducted prototype testing of a full-scale CNS thermosiphon system in Russia during 2003 (Thiering and Lu, 2006). The main purpose of the testing was to demonstrate the capability of the thermosiphon system in removing the heat from the moderator cell subject to varying coolant flow rates and varying heat loads and proving that there deuterium would remain in single phase. For the purpose of validation, numerical predictions are presented for the main test which is the heat removal capacity of the cryogenic coolant system under estimated operational heat loads. The schematic diagram of the experimental rig can be seen in Fig. 5. In addition to the CNS thermosiphon system itself, this figure shows the helium and deuterium buffer tanks which are coloured brown and green respectively. In the current validation study, these are not modelled during the CFD simulation. The relevant parameters for the experimental rig seen in Fig. 5 as well as the measurement uncertainties are summarised in Table 1. Since a closed circulation loop is considered for the flow of deuterium, no boundary conditions are required to be imposed. For the flow of helium, mass flow rates of 80 g/s are specified and at the inlet locations of HeF1 and HeF2 and temperatures of 19 K at the inlet locations of HeT1 and HeT2 of the thermosiphon system. During the tests, the power was increased in steps up to a total value of 5000 W. At each incremental power step, the power level was held steady for a period of 15 min. For computational efficiency and not over burdening the computational resources, pseudo steady-state calculations are performed for the validation of the CFD model. Two cases are simulated: one with a nominal total power of approximately 2000 W and the other with a nominal total power of approximately 4000 W. Relevant parameters measured during the tests at these two power levels are compared against the numerical predictions of the CFD model. The comparison between relevant parameters of PNPI’s experimental testing and the CFD model simulations at nominal powers of 2000 W and 4000 W for the k
e turbulence model can be seen in Tables 2 and 3 and at nominal power of 4000 W for the SST turbulence model in Table 4, respectively. It can be seen that there is very close agreement for all of the parameters at both power levels. The largest percentage difference is 14.795 and this is recorded for the total heat removal at 2000 W for the k
e turbulence model. It has been reported that there is a discrepancy of approximately 500 W between the heat input and the heat removed. There is a likelihood that this difference could be due to the other heat loads that are indirectly applied to the experimental rig as opposed to the heat load that is applied through the heating coils. These extra heat loads are likely a radiative heat load due to the large temperature difference between the CNS rig and its surroundings and a conductive heat load at the fixtures supporting the CNS thermosiphon system. At each successive power load there existed this discrepancy in the experimental results. It is therefore this discrepancy that have led to the largest percentage difference between the experimental and CFD results. In the current study, this difference is accounted by increasing the heat load on the helium by 10% so as

T. Pavlou et al. / Annals of Nuclear Energy 90 (2016) 135–147

141

Fig. 4. Cross-sectional view of heat exchanger alongside with the generated meshes of cylindrical slices as well as the top section with 133 tubes for the deuterium flow.

to best estimate the extra non-nuclear heat load. Nevertheless, most of the other differences are less than 5%. In Table 3, the deuterium average temperature from the k
e turbulence model is approximately 5% higher than the experimental results. However, the SST turbulence model demonstrates much closer agreement with the experimental results with a difference of only 2.67% (see Table 4). One possible reason for this discrepancy may be the SST model’s treatment of the near-wall regions particularly in predicting of the flow in the heat exchanger. The deuterium is in the range of low Reynolds number flow as it travels through the tubes of the heat exchanger. The SST turbulence model performs well in boundary layer flows Menter (1993). It uses the k
x turbulence model in the viscous sub-layer and gradually introduces the k
e turbulence model in the bulk flow away from the wall. The SST turbulence model shows that more heat being removed via the heat exchanger and hence a lower average temperature in the deuterium. The predicted results strongly suggest that the CFD model developed in this project is indeed capable of simulating the PNPI prototype testing conditions very accurately and is thus successfully validated. The SST turbulence model produces results that are in closer agreement with the experimental results than the k
e turbulence model. Nevertheless, the k
e turbulence model predictions are still within 5% of the experimental results.

5.2. Natural convection flow of deuterium in the circulation loop of the CNS thermosiphon system in OPAL Based on a validated CFD model, one of the most important aspect of this model is the ability to predict the direction of the natural convection deuterium flow in the thermosiphon. As the flow is driven by only density differential of the fluid within the closed circulation loop, the fluid flow is governed entirely by the fundamental physics of the buoyancy-driven problem. As can be seen in Fig. 6, the fluid flow rises along the ‘‘hot leg” and then falling through the heat exchanger in the ‘‘cold leg” of the CNS thermosiphon system. This effect is best demonstrated where the hot leg is coloured in red and the cold leg in green. Another important result from the model simulations is the predictions of the mass flow rate of the deuterium natural convection closed loop. Table 5 illustrates the comparison between the predictions of the k
e and SST turbulence models and theoretical calculated value. Not only did the flow of deuterium travels in the ‘‘correct” direction, the simulated mass flow rate is within 1% of the theoretically calculated value. This result is not surprising as the correct mass flow rate is essential to obtain a correct heat balance throughout the CNS thermosiphon system and the material properties adopted in the CFD simulation are based on a thorough review of the literature.

142

T. Pavlou et al. / Annals of Nuclear Energy 90 (2016) 135–147 Table 3 Comparison between measurements and predictions by the k
e turbulence model for nominal power of 4000 W.

DT1, (K) DT2, (K) DT3, (K) HeT1, (K) HeT2, (K) HeT3, (K) Tav, (K) HeQTot (W)

PNPI

CFD

% Difference

25.10 22.00 25.20 19.25 19.00 24.10 23.60 4300

26.00 22.64 26.23 19.01 19.00 23.98 25.21 4208

3.58 2.90 4.10
1.27 0.00
0.50 6.83
2.15

Table 4 Comparison between measurements and predictions by the SST turbulence model for nominal power of 4000 W.

DT1, (K) DT2, (K) DT3, (K) HeT1, (K) HeT2, (K) HeT3, (K) Tav, (K) HeQTot (W)

PNPI

CFD

% Difference

25.10 22.00 25.20 19.25 19.00 24.10 23.60 4300

25.08 22.48 25.15 19.01 19.00 23.96 24.23 4205


0.06 2.20
0.20
1.27 0.00
0.57 2.67
2.21

Fig. 5. Schematic diagram of the full-scale experimental setup of the CNS thermosiphon system at PNPI.

Hot leg Table 1 Experimental parameters for the PNPI tests.

DT1, (K) DT2, (K) DT3, (K) HeT1, (K) HeT2, (K) HeT3, (K) HeF1, (g/s) HeF2, (g/s) HeF3, (g/s) Tav (K) HeQ, (W) DQ, (W) HeQTot

Description

Measurement uncertainties

D2 temperature before HX D2 temperature after HX D2 temperature after moderator He temperature to moderator He temperature to HX He temperature after hot leg He mass flow rate to moderator He mass flow rate to HX He mass flow rate after hot leg D2 average temperature Heat load on He Heat load on D2 Total He heat removal

±0.2 ±0.2 ±0.2 ±0.2 ±0.2 ±0.2 ±5.0 ±5.0 ±5.0 ±0.2 ±50 ±50 ±50

Cold leg

Note: HX – heat exchanger, D2 – deuterium, He – helium.

Table 2 Comparison between measurements and predictions by the k
e turbulence model for nominal power of 2000 W.

DT1, (K) DT2, (K) DT3, (K) HeT1, (K) HeT2, (K) HeT3, (K) Tav, (K) HeQTot (W)

PNPI

CFD

% Difference

22.60 20.50 22.60 19.25 19.00 21.75 21.30 2500

22.83 21.39 23.03 19.01 19.00 21.52 22.36 2130

1.02 4.33 1.89
1.27 0.00
1.05 4.97
14.79

Fig. 6. Natural convection flow of deuterium through the ‘‘hot leg” and ‘‘cold leg” within the closed circulation loop of the CNS thermosiphon system represented by the predicted temperature.

143

T. Pavlou et al. / Annals of Nuclear Energy 90 (2016) 135–147 Table 5 Comparison of mass flow rate of deuterium between predictions by the k
e and SST turbulence models and theoretical calculated value. Mass flow rate (g/s) k
e model SST model Theoretical prediction

104.7 102.5 103.8

Table 6 Heat load predictions by the k
e and SST turbulence models by the helium flow.

k
e model SST model

Lower He (W)

Upper He (W)

Total (W)

2669 2461

1538 1744

4208 4205

5.3. Fluid flow and heat transfer within the moderator cell of the CNS thermosiphon system in OPAL A total heat load of about 4200 W is applied to the CNS thermosiphon system. Under steady-state operation the entirety of this load must be removed by the two helium coolant flows in order to achieve an energy balance throughout the system. The split in total heat removal between the upper helium coolant and the lower helium coolant for different turbulence models can be seen in Table 6. It can be seen that the SST turbulence model removes around 8% more heat through the upper helium flow i.e. the heat exchanger, than the k
e turbulence model. Fig. 7 shows the streamlines coloured by the helium velocity and temperature within the moderator cell. The helium coolant enters the moderator vertically through an annular region, passes through the displacer and then flows through the small gap between the moderator cell and the cell jacket. As can be seen in Fig. 7(a), the helium accelerates very quickly as it passes through the 90° bend before entering the displacer, resulting in velocities over 40 m/s. The helium swirls around inside the displacer with a clear region of recirculation near the top of the displacer. From streamlines coloured by velocity, there appears to be axial symmetry of the flow as it passes through the thin gap between the cell and the jacket. This is a very important result as it means that both sides should receive similar cooling. This effect is further elucidated in Fig. 7(b). The helium heats up inside the displacer as it swirls around and gradually heats up as it flows upwards through the small annular gap. The helium experiences an increase in temperature of approximately 6 K through the moderator cell. Fig. 8 depicts the collective temperature contour plots of the flow of deuterium where three of these planes are equal-spaced in the XY plane with the central plane in the geometric centre of the moderator and the fourth plane is located in the geometrical centre of the moderator on the YZ plane. These planes are highlighted separately in Fig. 9 accompany with the corresponding velocity contour plots. The deuterium enters the moderator vertically through a pipe at a temperature of approximately 22 K and a mass flow rate of approximately 100 g/s. It can be seen in Fig. 9 that the deuterium enters the chamber at a uniform velocity of approximately 0.7 m/s and continues downwards until hitting the curved bottom of the moderator cell. Fig. 9(b), (d) and (f) shows that the majority of the flow then accelerates tangentially along the aluminium cell towards the bottom of the cell and a small portion accelerates back up the cell wall. These three figures also show that there is a strongly swirling region located in the geometric centre of the cell. Fig. 9(a), (c) and (e) shows that there is an increase in temperature from the bottom to the top of the cell. This is expected as the cold fluid enters the cell and blasts its way to the bottom before slowing down after the initial acceleration along the walls. It then slowly swirls around in the main body of the moderator, giving it time to heat up. Fig. 9(g) and (h) shows that the nearwall regions of the displacer are cooler than the bulk flow. This is due to the heat transfer through the aluminium from the warmer deuterium to the cooler helium. Fig. 10 shows the temperature distribution of the aluminium solid region of the moderator cell. The lowest temperatures are

(a)

(b) Fig. 7. Streamlines of helium flow within the moderator cell coloured by: (a) helium velocity and (b) helium temperature.

in the pipe through which the helium coolant enters the chamber and the highest temperatures are around the deuterium outlet at the top of the cell. The maximum wall temperature of the aluminium is 26.6 K which is below the boiling point of the liquide deuterium. This indicates that under these steady-state conditions localised boiling is highly unlikely to occur in any of the near-wall regions.

144

T. Pavlou et al. / Annals of Nuclear Energy 90 (2016) 135–147

is at its fastest, namely in the pipes entering and leaving the displacer. There is also a region of high heat transfer on the back side of the displacer where the cold helium flow hits the wall and accelerates towards the bottom of the displacer. 5.4. Fluid flow and heat transfer within the heat exchanger of the CNS thermosiphon system in OPAL

Fig. 8. Orientation and different cross-sectional planes of temperature contours cutting the moderator cell.

The contour plots of the wall heat transfer coefficient, as seen in Fig. 11, are an indicator of where the heat is being transferred from the aluminium solid. The areas with the highest wall heat transfer coefficient on the deuterium side can be seen in red in Fig. 11(a). There are two main regions: one at the bottom of the moderator where the incoming cold stream hits the aluminium wall and accelerates along the surface and another at the outlet. Both of these regions are where the flow is the fastest. Areas of high wall heat transfer coefficient are also seen in Fig. 11(b) where the flow

The heat exchanger is where the upper helium coolant flow removes the heat from the hot deuterium natural convection flow. Temperature contour plots of both the shell-side and the tube-side fluid can be seen in Fig. 12. Not surprisingly, the hot deuterium fluid enters at the top of the heat exchanger and cools as it flow downwards through the tubes, while the helium fluid does the opposite in the shell. Both fluids experience a change in temperature of approximately 3 K across the heat exchanger. Although not shown in a figure, the average velocity of the deuterium remains constant at approximately 0.09 m/s. It takes the deuterium of approximately 11 s to pass through the heat exchanger. The velocity fields of the helium coolant, as seen in Fig. 13, indicate that the heat exchanger is operating near its optimum. The velocity streamlines in Fig. 13(a) pass smoothly from baffle to baffle without significant recirculation. The best performance of a heat exchanger is obtained when there are no recirculation regions (Mohammadi et al., 2009). The cyclic behaviour of the shell-side flow can be seen across the whole heat exchanger in Fig. 13(b). The helium velocity peaks to a value of approximately 8 m/s as it flow around each baffle with the bulk of the flow moving at approximately 4 m/s. 5.5. Phenomenological discussion on the subcooled boiling characteristics of liquide deuterium In the event that the supply of helium to the CNS thermosiphon system is reduced, there is a likelihood that the insufficient

(a)

(c)

(e)

(g)

(b)

(d)

(f)

(h)

Fig. 9. Temperature contours: (a) front plane, (c) central plane, (e) back plane and (g) YZ plane and velocity vectors: (b) front plane, (d) central plane, (f) back plane and (h) YZ plane of the moderator cell.

145

T. Pavlou et al. / Annals of Nuclear Energy 90 (2016) 135–147

(a)

(a)

(b)

(b)

Fig. 10. Distribution of temperature within the aluminium solid material of the moderator cell.

Fig. 11. Distribution of wall heat transfer coefficient within the moderator cell.

removal of heat from the moderator cell may lead to the occurrence of subcooled boiling as the temperature of the cell wall raises above the saturation temperature of the liquide deuterium. To aptly simulate the subcooled boiling phenomena, considerations for modelling requires the physics to be predicted not only at the heated cell wall but also away from the cell heated wall. In the context of CFD, modelling subcooled nucleate boiling can be broadly classified into two categories: (1) Heat transfer mechanistic modelling at the heating wall and (2) Two-phase flow and bubble behaviours in the bulk liquid away from the heating wall. For the first category, suitable models are developed to predict the heat transfer at the heated wall. Amongst the many different models, a mechanistically approach of partitioning the wall heat flux such as the model proposed by Kurul and Podowski (1990) and Podowski (2012) is commonly adopted which entails partitioning the wall heat flux into three heat flux components of single-phase convection, transient conduction and evaporation. Yeoh et al. (2008) extended the model to account for the transient conduction due to sliding bubbles along the heated wall. The

closure of these models, which require accurate predictions of the number of active nucleation sites at the cell heated wall, bubble sizes and their frequencies leaving the cell heated wall, remains the greatest challenge due to the lack of measurement data especially for liquids at cryogenic temperatures. For the second category, the prediction of local bubble sizes in the bulk liquid flow is expected to be strongly influenced by the complex bubble behaviours near the heated wall. In essence, complex mechanistic behaviours of bubble coalescence and bubble collapse due to condensation in subcooled nucleate boiling exist during the boiling of cryogenic liquids (Li et al., 2006, 2009). Physical insights strongly indicate that the numerical solutions of two-phase cryogenic systems can be significantly improved through consideration of the bubble mechanistic behaviours in determining the local bubble size such as the Sauter mean bubble diameter, which appears prominently in the interphase exchange terms of mass, momentum and energy. Because of the importance of accurately determining the Sauter mean bubble diameter, which strongly governs the inter-phase heat and mass transfer through the

146

T. Pavlou et al. / Annals of Nuclear Energy 90 (2016) 135–147

(a)

(b)

Fig. 12. Distribution of temperature: (a) deuterium (tube-side) and (b) helium (shell-side) within the heat exchanger.

(a)

(b) Fig. 13. Cross-sectional planes of: (a) streamlines of helium flow within a cylindrical slice and (b) velocity vectors of the entire heat exchanger.

interfacial area as well as the momentum drag, the potential to extend the modelling of boiling of cryogenic liquids in order to better predict the non-uniform bubble size distribution via the

population balance approach (Yeoh and Tu, 2004, 2005) is of enormous significance. The presence of subcooled boiling within the moderator cell as well as other components of the CNS

T. Pavlou et al. / Annals of Nuclear Energy 90 (2016) 135–147

thermosiphon system will be investigated and analysed in our future work with respect to the consideration of the mechanistic wall heat flux partitioning and population balance models. 6. Conclusions Investigating the relevant fluid flow and heat transfer processes within the CNS thermosiphon system has given us an understanding of the fundamental physics of the buoyant-driven flow of liquide deuterium and the force convection flow of helium coupled with the conduction heat transfer across different materials. For the first time, a sophisticated CFD model is developed to solve a closed-circulation loop of the flow of liquide deuterium within the moderator cell, heat exchanger and two inter-connecting pipes representing the hot leg and cold leg of CNS thermosiphon system. A novel technique for modelling the shell and tube heat exchangers is adopted that allows both the shell-sides and tube-sides to be meshed without needing to generate a mesh for the intricate solid region. This permits the baffles to be included within the heat exchanger for the CFD model. The full-scale CNS thermosiphon system is validated against experimental data obtained from PNPI prototype tests. The SST turbulence model demonstrates a much closer agreement than the k
e turbulence model with the experimental data with a difference of 2.67%. However, the k
e turbulence model predictions are still within 5% of the experimental data of the predicted temperatures. For the simulation of the CNS thermosiphon system under nominal heat load, the CFD model is able to predict the ‘‘correct” direction of the flow of deuterium within the CNS thermosiphon system by rising along the hot leg after exiting the moderator cell and falling along the cold leg after it is cooled by the heat exchanger before entering the moderator cell. The CFD predicted mass flow rate is within 1% of the theoretically calculated value as the CFD model has managed to obtain a correct heat balance throughout the CNS thermosiphon system. The CFD model also demonstrates the correct split in total heat removal between the upper and lower helium flow of the CNS thermosiphon system with the majority of the heat being removed in the moderator cell. Future work in the modelling of the CNS thermosiphon system will include the prediction of voids due to subcooled boiling within the different components of the system through the consideration of the mechanistic wall heat flux partitioning and population balance models, as well as characterising the system’s behaviour during operational transients. References Ageron, P., 1980. Cold neutron sources at ILL. Nucl. Instrum. Methods Phys. Res., Sect. A 284 (1), 197–199. Alsmiller Jr, R.G., Lillie, R.A., 1992. Design calculations for the ANS cold neutron source Part II. Heating rates. Nucl. Instrum. Methods Phys. Res., Sect. A 321 (1), 265–270. Axmann, A., Böning, K., Rottmann, M., 1997. FRM-II: the new German research reactor. Nucl. Eng. Des. 178 (1), 127–133.

147

Buscaglia, G.C., Dari, E.A., Martín, J.E., Arnica, D.L., Bonetto, F.J., 2014. Finite element modeling of liquid deuterium flow and heat transfer in a cold-neutron source. Int. J. Comput. Fluid Dyn. 18 (5), 355–365. Gaubatz, W., Gobrecht, K., 2000. The FRM-II cold neutron source. Phys. B: Condens. Matter 276, 104–105. Gil, A., 2011. New sources and instrumentation for neutron science. J. Phys.: Conf. Ser. 289, paper 012034. Grósz, T., Hargitai, T., Mityukhlyaev, V.A., Rosta, L., Serebrov, A.P., Zaharov, A.A., 2000. Thermo-hydraulic test of the moderator cell of LH2 cold neutron source at BNC. Phys. B: Condens. Matter 276, 214–215. Kennedy, S.J., 2006. Construction of the neutron beam facility at Australia’s OPAL research reactor. Phys. B: Condens. Matter 385, 949–954. Kim, Y.K., Lee, K.H., Kim, H.R., 2008. Cold neutron source at KAERI, Korea. Nucl. Eng. Des. 238 (7), 1664–1669. Kurul, N., Podowski, M.Z., 1990. Multi-dimensional effects in forced convection subcooled boiling. Proceedings of the 9th Heat Transfer Conference, vol. 2. Hemisphere Publishing Corporation, Jerusalem, Israel, pp. 21–26. Launder, B.E., Spalding, D.B., 1974. The numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng. 3, 269–289. Lemmon, E.W., McLinden, M.O., Friend, D.G., 2005. Thermophysical properties of fluid systems, NIST chemistry webbook, NIST standard reference database, 69. Li, X.D., Wang, R.S., Huang, R.G., Shi, Y., 2006. Numerical investigation of boiling flow of nitrogen in a vertical tube using the two-fluid model. Appl. Therm. Eng. 26, 2425–2432. Li, X.D., Wei, W., Huang, R.G., Shi, Y., 2009. Numerical and experimental investigation of heat transfer on heating surface during subcooled boiling flow of liquid nitrogen. Int. J. Heat Mass Transfer 52, 1510–1516. Lopez-Rubio, A., Gilbert, E.P., 2009. Neutron scattering: a natural tool for food science and technology research. Trends Food Sci. Technol. 20 (11), 576–586. Lu, W., Thiering, J.R., 2012. Using a multi-parameter monitoring methodology to predict failures in the cryogenic plant of the cold neutron source at Australia’s OPAL reactor. Adv. Cryog. Eng.: Trans. Cryog. Eng. Conf.-CEC 57, 1537–1542. Marquardt, E.D., Le, J.P., Radebaugh, R., 2002. Cryogenic material properties database. Cryocoolers 11, 681–687, Springer. Menter, F.R., 1993. Zonal two equation k–x turbulence models for aerodynamics flows, AIAA paper 93-2906. Misiorek, H., Zakrzewski, T., Rafalowicz, J., 1981. The influence of neutron irradiation on the thermal conductivity of aluminum in the range 5–50 K. Int. J. Thermophys. 2 (4), 341–353. Mohammadi, K., Heidemann, W., Müller-Steinhagen, H., 2009. Numerical investigation of the effect of baffle orientation on heat transfer and pressure drop in a shell and tube heat exchanger with leakage flows. Heat Transfer Eng. 30 (14), 1123–1135. Ozden, E., Tari, I., 2010. Shell side CFD analysis of a small shell-and-tube heat exchanger. Energy Convers. Manage. 51 (5), 1004–1014. Podowski, M.Z., 2012. Towards mechanistic modelling of boiling heat transfer. Nucl. Eng. Technol. 44, 889–896. Richardson, I.A., Leachman, J., 2012. Thermodynamic properties status of deuterium and tritium. Adv. Cryog. Eng.: Trans. Cryog. Eng. Conf.-CEC 57, 1841–1848. Sun, S., Liu, Y.T., Wang, H.L., Chen, D.F., Feng, Q.K., Yu, X.L., 2012. Thermal and mechanical characteristics of the moderator cell in a cold neutron source. Ann. Nucl. Energy 50, 1–7. Thiering, J.R., Lu, W., Ullah, R., 2006. Commissioning of the OPAL reactor cold neutron source. In: Pacific Basin Nuclear Conference, Australian Nuclear Association. Wang, Q.W., Chen, Q.Y., Chen, G.D., Zeng, M., 2009. Numerical investigation on combined multiple shell-pass shell-and-tube heat exchanger with continuous helical baffles. Int. J. Heat Mass Transfer 52 (5), 1214–1222. Yeoh, G.H., Tu, J.Y., 2004. Population balance modeling for bubbly flows with heat and mass transfer. Chem. Eng. Sci. 59, 3125–3139. Yeoh, G.H., Tu, J.Y., 2005. Thermal hydraulic modeling of bubbly flows with heat and mass transfer. AIChE J. 51, 8–27. Yeoh, G.H., Cheung, S.C.P., Tu, J.Y., Ho, M.K.M., 2008. Fundamental consideration of wall heat partition of vertical subcooled boiling flows. Int. J. Heat Mass Transfer 51, 3840–3853.