Abstract
COMSOL Multiphysics has been used to conduct thermal-hydraulic analysis in multiple nuclear applications. The aim of this study is to benchmark the prediction accuracy of COMSOL Multiphysics in performing thermal-hydraulic analysis of TRIGA (Training, Research, Isotopes, General Atomics) reactors such as the Geological Survey TRIGA Reactor (GSTR) by comparing its predictions with RELAP5 (a widely used code in nuclear thermal-hydraulic analysis) results and experimental data. The GSTR type is Mark I with a full thermal power of 1 MW, and it resides at the Denver Federal Center (DFC) in Colorado. The numerical investigation of the present work is carried out by developing single-subchannel thermal-hydraulic models of the GSTR utilizing RELAP5 and COMSOL codes. The models estimate the temperatures (fuel, outer clad, and coolant) and water flow patterns in the core as well as fuel element powers at which void starts to form within the coolant subchannels. Then, these models’ predictions are quantitatively evaluated and compared with the measured data.
1. Introduction
The Geological Survey TRIGA Reactor (GSTR) situated in Lakewood, Colorado, is a 1 MW Mark I. It is operated since 1969 by the US Geological Survey (USGS) to provide many services such as fission track experiments, neutron activation analysis, radioisotope production, neutron irradiation, neutron radiography, and delayed neutron analysis [1]. To safely operate in steady-state and transient modes, the GSTR is fueled with uranium-zirconium hydride fuel that has a prompt negative temperature coefficient of reactivity. During pulse mode operations, this reactor can operate with a $3 maximum insertion reactivity [2].
COMSOL Multiphysics has been utilized in many nuclear reactors investigations to perform multiphysics simulations such as thermal hydraulics. The aim of the present work is to study the prediction accuracy of COMSOL Multiphysics in reactor cores with natural convection cooling such as TRIGA reactors by benchmarking its predictions with the experimental data as well as RELAP5 (extensively used for thermal-hydraulic evaluations of TRIGA reactors) results. Previously, a detailed neutronic model of the GSTR core is developed using MCNP5 in [3, 4], which calculates the neutronic parameters of the GSTR core. Then, the models in the present study use these neutronic parameters calculated by the MCNP5 model to predict the temperatures and void formation for the subchannel.
Thus, the goal of this numerical and experimental study is to develop single-subchannel thermal-hydraulic models using a selection of applicable thermal-hydraulic codes such as RELAP5 and COMSOL, calculate the temperatures (fuel, outer clad, and coolant) and fluid flow direction with different heat rates, ranging from 5 to 30 kW, conduct an experiment to measure the water temperature in the flux monitor hole situated between B- and C-rings, assess the models’ predictions and compare them with the measured data, and benchmark the COMSOL Multiphysics single-subchannel thermal-hydraulic evaluation of a TRIGA Reactor with RELAP5 results and experimental data.
2. Geological Survey TRIGA Reactor
The GSTR situated in Lakewood, Colorado, is a 1 MW Mark I TRIGA (Figure 1). This reactor can operate with thermal powers up to 1 MW maximum in steady-state mode operations and up to $3 maximum reactivity insertions during pulse mode operations [2]. The GSTR core, like other Mark I TRIGA reactors, is composed of seven rings (labeled A through G) and situated at the bottom of a tank (2.34 m inner diameter and 7.50 m high). The water inside this tank yields effective radiation shielding [5].
(a)
(b)
This reactor is fueled with uranium-zirconium hydride (U-ZrH) that has a prompt negative temperature coefficient of reactivity and <20% concentration of [6]. The cladding material of the GSTR fuel elements is made of stainless steel (located in the B-, C-, D-, E-, and G-rings) or aluminum (located in the F- and G-rings) (Figure 2) [3]. To hold the fuel elements in their specified positions and in vertical orientation, the core top and bottom grid plates are used. A detailed description of the GSTR including core, control rods, fuel elements, graphite reflector, and core support grids can be retrieved from [7].
The GSTR core is safely cooled by natural convection with steady-state thermal powers up to 1 MW. While operating in steady-state mode, the NRC requires a continuous monitoring of the pool temperature. The maximum allowable pool temperature is [8]; thus, tank cooling is needed to run the reactor for long time at powers above 100 kW. As Mark I TRIGA reactors, primary cooling loop pump pulls out the water from the core to cool it down by removing the heat using a heat exchanger. Then, water in the secondary loop is cooled by a cooling tower, which consists of sprayers and fans to cool the water circulating through the secondary cooling system [9]. A schematic representation of the GSTR coolant circuit is illustrated in Figure 3.
3. Numerical Methods
This section presents a detailed overview of the single-subchannel thermal-hydraulic models developed in this research. To correctly compare between the results of these models, the same boundary and initial conditions described in detail in the following paragraphs are incorporated within each model.
The numerical simulations are carried out using a one- and three-dimensional domains utilizing RELAP5 and COMSOL, respectively. All thermodynamic properties for the water and fuel are prescribed as a function of temperature. The gravity term is acting in the z direction as the flow is moving vertically upward in this direction. The inlet and outlet are specified as pressure and temperature boundary conditions, and the mass flow and fluid directions of the coolant depend on the natural convection caused by the heat flux generated within the fuel element. The remaining boundaries are as follows: the symmetry boundary at the connected adjacent subchannels and nonslip wall boundary with equivalent heat flux (see Figure 4 for the axial and radial power profiles of the fuel element). The generated power in the fuel elements depends on the fuel cladding material, location, and reactor-operating power. Figure 5 presents the power generated in each fuel element at reactor full power [3]. Simulation are considered converged when residuals fall below 10−6 threshold value. Moreover, an additional check is also made to ensure that mass and energy balances are globally satisfied.
(a)
(b)
This research also calculates the average cross-sectional area of the coolant for each GSTR ring to anticipate the temperatures and void development within coolant subchannels at each ring (see Table 1 and Figure 6). The cross-sectional areas and GSTR fuel elements pitches presented in Table 1 are computed utilizing the following equations:
3.1. RELAP5/MOD3.3 Model
The single-subchannel thermal-hydraulic model is developed utilizing RELAP5/MOD3.3 Patch 04, as this code has been widely used in the literature in [10–16] to do TRIGA Reactor thermal-hydraulic analysis. The RELAP5 model developed in the present research is consistent with the modeling approach used by these previous efforts. The RELAP5 model is composed of a cold leg, hot subchannel, horizontal connector, coolant source, and coolant sink (Figure 7). The cold leg is included to impose the pressure differential between the cold coolant at the inlet and the coolant passing through the heated subchannel. The cold leg has the same height and volume as the hot subchannel to communicate correctly to the hot subchannel. This modeling approach is important to differentiate between the gravitational forces imposed on the cold leg versus that imposed on the hot subchannel permitting for the method of cooling to be only natural convection. The hot subchannel is the volume that contains the heat structure (fuel element). The horizontal connector is a nonphysical connector to permit communications between cold leg and hot subchannel volumes during simulation. The coolant source and sink (time dependent volumes) are incorporated in the model to impose the temperature and temporal pressure boundary conditions during the computational simulations.
The hot subchannel volume represents the coolant flow due to the natural convection driven by the heat generated in the fuel element. This volume is modeled to have 24 axial connected volumes. Thus, the flow in the pipe components is driven solely by the density gradients resulting from variations in temperature resulting from the input power. This approach has been used in multiple previous models of TRIGA-type reactors [3, 10, 11, 13–15]. The heat structure, a two-dimensional component, is represented with 24 radial mesh cells and 24 axial mesh cells to correctly model the fuel element geometry and power profiles predicted by the MCNP5 model. Figure 8 displays the two-dimensional representation of the heat structure in the model.
3.2. COMSOL Multiphysics Model
This research also develops a three-dimensional single-subchannel thermal-hydraulic model (shown in Figure 9) using COMSOL Multiphysics code, as this code has the capacity to model fluid flow on heated volumes (fuel) in one-, two-, and three-dimensional domains. Three different meshes are generated for the COMSOL mesh convergence tests. The total number of the quadrilateral control volumes for these meshes are 4,306,584, 2,725,000, and 1,741,360 for fine, medium, and coarse meshes, respectively. The temperature predictions results from different meshes show a satisfactory monotonic convergence, and the profiles obtained from these runs are overlapping each other, which means even the coarse mesh is adequately refined for the present numerical simulation and no further improvement is to be anticipated by using the medium or fine meshes. Hence, the results obtained from the three different meshes should give enough confidence that mesh convergence is satisfied and numerical predictions are reliable. Therefore, as no significant difference are noticed between the predictions of fine, medium, or coarse meshes, the coarse mesh is used to reduce the associated computational cost. The gap and clad materials in the three-dimensional COMSOL model are substituted with a single effective material to reduce the solution time and required number of mesh cells (see Figure 10). The effective thermal conductivity of the composite (gap and clad) is calculated as a function of the gap and clad thermal resistance as shown in the following equation [17]:
(a)
(b)
The heat flux rate going through the gap and clad layers equals to the heat flux going through the homogenous material with an effective thermal conductivity [17]:
Equating equation (6) with equation (7) leads to [17]
Solving for effective thermal conductivity in (equation (8)) results in [17]
Mass of the homogenous material with effective density equals to the mass of the gap and clad regions; therefore, the total mass of the composite equals [17]
The total mass of the homogenous material with effective density is equal to
Setting equation (16) equal to equation (11) leads to
The heat capacity of the homogenous material with effective heat capacity equals to the heat capacity of the gap and clad regions; therefore, the total heat capacity of the composite equals [17]
The total heat capacity of the homogenous material with effective heat capacity is equal to [17]
Setting equation (13) equal to equation (14) and solving for effective heat capacity leads to [17]
4. Model Results and Discussion
The single-subchannel thermal-hydraulic models that are developed in this research predict the temperatures of the fuel, outer clad, and coolant to benchmark the COMSOL predictions with RELAP5 results. These temperatures are predicted at normal inlet temperature () and maximum allowable inlet temperature (). To study the effects of changing inlet temperature boundary condition on the temperatures’ predictions, Figure 11(a) displays the temperatures predicted by the RELAP5 model versus the power generated by a fuel element located in the B-ring. The resulting temperatures consistently increase by , as the inlet temperature increases from to , as the impact of changing the inlet temperature is substantially linear (see Figure 11(a)). As changing inlet temperature results in largely linear change, the temperatures can be predicted even for different inlet temperatures without needing to rerun the models for the new inlet temperature. In order to anticipate the temperatures and void formations inside the core at the maximum allowable coolant temperature, the remaining results in this section are presented with an inlet temperature of .
(a)
(b)
This research also studies the effects of changing the flow cross-sectional areas (Table 1) in these single-subchannel models on the temperatures’ predictions. Figure 11(b) presents these temperatures predicted by the RELAP5 model versus the fuel element power for each GSTR ring at inlet temperatures of . As demonstrated in Table 1, the flow cross-sectional areas from B- through F-rings differ only marginally; thus, the calculated temperatures in these rings are roughly identical for the same fuel element powers. Differently, the G-ring flow cross-sectional area is higher and therefore results in different temperatures in this ring. As changing the flow cross-sectional areas of coolant in the models to predict temperatures in the B- through F-rings results in extremely minor temperatures’ differences, the temperatures in these rings versus the fuel element power can be anticipated without needing to rerun the models more than one time.
4.1. Benchmarking with RELAP5/MOD3.3 Results
4.1.1. Temperature Predictions
Figure 12 presents the maximum fuel temperatures calculated by all of the models versus the fuel element power. The predicted temperatures by these models are approximately the same with minor differences. This present research also runs the models to produce a graph that displays the maximum outer clad temperatures versus the fuel element power (Figure 13). The results of these models are approximately similar. The outer clad temperatures calculated by the COMSOL model are approximately the same as the temperatures predicted by the RELAP5 with slight differences. The COMSOL model correctly calculates the outer clad temperature, as the outer clad temperatures are affected by the heat transfer from the outer clad to coolant. The small differences in fuel and outer clad temperatures predicted by RELAP5 and COMSOL models in Figures 12 and 13 are due to the variations in the analysis of flow behavior including field equation formulations, numerical scheme order, and geometry domain (one-dimensional for RELAP5 and three-dimensional for COMSOL).
(a)
(b)
(a)
(b)
Due to the hydrostatic pressure discussed in Section 3, the coolant boiling temperature inside the GSTR core increases to be in the range of (top of the core) and (bottom of the core) due to the water column above the core. Formation of void within the core starts at outer clad temperatures higher than the saturation temperature (between and ) as a result of the onset of incipient boiling. Based on this void-initiation temperature, the COMSOL model anticipate that void development happens in the GSTR core when the power generated by the fuel element is above 17.90 kW as a result of outer clad temperatures that are higher than the saturation temperature. The COMSOL models’ predictions of the outer clad temperatures are conservative, as these models do not include two-phase physics. Lacking two-phase physics does not completely invalidate the COMSOL models’ predictions, as convective heat transfer initially increases with the onset of incipient boiling. This is considered in detail in [18]. To elaborate on void formation within the GSTR core and to compare with the COMSOL model that lack two-phase physics, this paper runs the RELAP5 model to calculate the local maximum void fraction within the thermal-hydraulic subchannels versus the fuel element power. The RELAP5 model predicts that void starts to form within the thermal-hydraulic subchannels at fuel element powers more than 18.90 kW.
Figure 14(a) shows the outlet temperatures calculated by the RELAP5 and COMSOL models versus the fuel element power. The average coolant temperature’s predictions by the RELAP5 and COMSOL models differ slightly with maximum variations. There is a maximum of differences between the COMSOL and RELAP5 results as shown in Figure 14(a). As explained in Section 3.1, the subchannel flow volume in RELAP5 code is discretized by specifying the subchannel cross-sectional flow area and subchannel length. The subchannel length is discretized into 24 nodes. The outlet temperature in RELAP5 model is the temperature of the last node volume at the outlet of the subchannel. On the contrary, the outlet temperature in the three-dimensional COMSOL model is calculated by selecting the outlet surface of the subchannel geometry and computing the average temperature of this surface using postprocessing tools within the COMSOL model. Different schemes of calculating outlet temperatures in RELAP5 and COMSOL models lead to slightly different predicted outlet temperatures.
(a)
(b)
Figure 14(b) presents the average coolant temperatures for each of the rings at reactor full power for fuel elements with stainless steel cladding in the B-, C-, D-, E-, and G-rings and with aluminum cladding in the F-ring. The highest outlet temperatures happen in the B- and C-rings and this is anticipated, as the fuel elements in these rings generate the highest power in the GSTR core. As displayed in Figure 14(b), the average outlet temperature predicted by the RELAP5 and COMSOL has a big difference for the G-ring as the subchannel flow area and fuel pitch in the G-ring is higher than all other rings (see Table 1) due to space between the graphite reflector that is surrounding the core and the fuel elements in the G-ring. When incorporating this larger flow area in the three-dimensional COMSOL model, it predicts higher temperatures than the temperatures predicted by the one-dimensional RELAP5 model. Increasing the flow cross-sectional area in the 24 axial nodes of the subchannel in the RELAP5 model reduces the outlet temperature, but the drop in outlet temperature is higher than the drop in the COMSOL model. The RELAP5 model predicts higher mass flow rate in the G-ring when increasing the cross-sectional flow area as shown in Figure 15, which leads to lower predicted outlet temperature in the RELAP5 model than outlet temperature predicted by the COMSOL model.
This work also generates three figures to predict temperatures of the GSTR including fuel (Figure 16(a)), outer clad (Figure 16(b), and average outlet (Figure 16(c)) temperatures at each GSTR ring versus the reactor power. To prepare these Figures, graphs of fuel, outer clad, and outlet temperatures versus the fuel element power (similar to Figures 12(a), 13(a), and 14(a)) are generated for stainless steel-clad (located in the B-, C-, D-, E-, and G-rings) and aluminum-clad (located in F- and G-rings) fuel elements. As shown in Figure 16, the hottest fuel elements in the GSTR core are located in the B- and C- rings, as these elements generate the highest power [7].
(a)
(b)
(c)
4.1.2. Coolant Flow
The coolant flow and average mass flow rate calculated by the single-subchannel models are presented in this subsection. In these single-subchannel models, water (coolant) tends to flow axially from the bottom of the subchannel to the top of the subchannel due to natural convection driven by the heat generated within the fuel element, with no interaction with the neighboring subchannels. Figure 17 presents the two-dimensional coolant flow direction predicted by the models.
In the single-subchannel models, the coolant subchannel is not connected to the neighboring thermal-hydraulic subchannels. In real sense, the subchannels in the core of the reactor are interlinked, thus the coolant could flow between neighboring subchannels, based on variations of the coolant density in the core. This is particularly true within nonuniform power density cores, for instance, the GSTR. Single-subchannel models restrict the flow directions, as the heat produced from neighboring fuel elements does not contribute to the buoyancy forces within the subchannel (coolant does not flow to or from the adjacent subchannels), which may limit the coolant flow rate and flow directions in these models. For instance, driving forces in B-, C-, and D-rings, which have the highest power density, do not contribute to the flow in the other subchannels. Figure 15 displays the calculated average mass flow by the models in each ring of the GSTR. This research verifies that the calculated power from equation (16) using the predicted average mass flow rate (Figure 15) and the average inlet and outlet temperatures (Figure 14) is the same as the inputted heat flux in the models:
4.2. Benchmarking with Experimental Data
This present work conducts an experiment to measure the coolant temperature inside the core at five vertical positions (see Figure 18) with reactor steady-state full power and an inlet temperature of (inlet temperature at the day of experiment). This experiment was done by inserting the thermocouple (with 1.8 m sheath to permit temperature measuring within the core) inside the flux monitor hole. As the flux monitor hole is situated between the B- and C-rings (see Figure 19), average of the predicted coolant temperatures at B- and C-rings with an inlet temperature of (same as experiment) by the RELAP5 and COMSOL models is calculated and then compared it with the measured coolant temperatures.
Figure 20 compares the measured coolant temperature at five vertical locations inside flux monitor hole with the predicted coolant axial temperatures by the RELAP5 and COMSOL models. This experiment shows that the measured temperature increases versus the core height. RELAP5 and COMSOL models also predict that coolant temperature rises along the core height with maximum temperature at the subchannel outlet. The measured outlet temperature in the flux monitor hole equals to 41.7 , while the predicted outlet temperatures by the RELAP5 and COMSOL models are and , respectively. These differences are due to the flow limitation of single-subchannel thermal-hydraulic models. Within the single-subchannel models, the flow of the coolant occurs directly from the subchannel’s bottom part to the subchannel’s top part due to natural convection caused by the heat generated within the fuel element. In real sense, the subchannels in the core of the reactor are interlinked; thus, the coolant could flow between neighboring subchannels, based on variations of the coolant density in the core. This is particularly true within nonuniform power density cores, for instance, the GSTR. To study the effects of flow patterns inside the core and improve the temperatures predictions, future research will develop a detailed computational fluid dynamics (CFD) multichannel thermal-hydraulic model.
5. Conclusions and Future Work
The paper presents a numerical and experimental study to develop single-subchannel thermal-hydraulic models using a selection of applicable thermal-hydraulic codes such as RELAP5 and COMSOL, calculate the fuel, outer clad, and coolant temperatures and fluid flow direction with different heat rates, ranging from 5 to 30 kW, conduct an experiment to measure the coolant temperature in flux monitor hole situated between B- and C-rings of the GSTR core, assess the models’ predictions and compare them with the measured data, and benchmark the COMSOL Multiphysics single-subchannel thermal-hydraulic evaluation of a TRIGA Reactor with RELAP5 results and experimental data.
The temperatures and reactor powers at which void starts to form predicted by RELAP5 and COMSOL models are approximately similar with minor differences. The experiment to measure the coolant temperature inside the GSTR core at five vertical positions in the flux monitor hole (located between B- and C-rings) shows that the measured temperature increases as a function of core height. RELAP5 and COMSOL models, in addition, predict that coolant temperature rises along the core height with maximum temperature at the outlet. The measured outlet temperature in the flux monitor hole equals to 41.7 , while the predicted outlet temperatures by the RELAP5 and COMSOL models are and , respectively. These differences are due to the flow limitation of single-subchannel thermal-hydraulic models without taking into account the coolant cross flow between neighboring subchannels, based on variations of the coolant density in the core.
The next step of this research is to improve the temperature and void formation predictions by developing multichannel thermal-hydraulic models, as these models will incorporate the coolant cross flow between adjacent subchannels, depending on coolant density variations within the GSTR core. This will improve the models’ accuracy, as the GSTR has a nonuniform power density core with multiple coolant flow paths. In contrast to the single-subchannel thermal-hydraulic models, multichannel thermal-hydraulic models will consider the complex heat transfer and fluid flow patterns within the core to more precisely predict the temperatures and fluid flow paths. A detailed experimental verification of the temperature and flow conditions in the GSTR core will also improve the validation of the different models. The results of both of these efforts will be presented in future papers.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.