Thermal hydraulics modeling of On-line refueling for the Advanced High Temperature Reactor (AHTR)

Nuclear Engineering and Design 358 (2020) 110440

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

Thermal hydraulics modeling of On-line refueling for the Advanced High Temperature Reactor (AHTR)

T

P. Avignia, , B. Petrovica ⁎

a

Nuclear and Radiological Engineering, Georgia Institute of Technology, 770 State St., Atlanta, GA 30332-0745, United States

ARTICLE INFO

ABSTRACT

Keywords: AHTR FHR TRISO fuel Plank fuel Online refueling ANSYS Fluent CFD

As part of the effort to move forward the Fluoride salt-cooled High-temperature Reactor (FHR) technology, the Advanced High Temperature Reactor (AHTR) is being developed at Oak Ridge National Laboratory (ORNL) and several other institutions. Due to its plank fuel design and inherently low heavy metal loading, which challenges fuel utilization, online refueling is considered as an option for improving the economic viability of this reactor. This work presents the thermal hydraulic modeling of the online refueling for the AHTR. Analyses at the single channel and fuel assembly level have been performed in order to develop a model of the reactor. The reactor modeling has been integrated with CFD studies for characterizing the steady state conditions and the refueling transient. The simulations of the transient have been complemented by the analysis of the stability of the refueling operation. This work demonstrates the viability of online refueling from a thermal hydraulic standpoint and develops a modeling approach for this type of operational transient.

1. Introduction and background information

1.1. AHTR reactor design features

The Advanced High Temperature Reactor (AHTR) is a Fluoridecooled High-temperature Reactor (FHR) moderated by graphite and cooled by FLiBe (Holcomb, 2011; Varma, et al., 2012). The large amount of graphite required to thermalize neutrons, in relation to the amount of fuel loaded in the core, results in a reduced heavy metal loading, which in the batch reloading mode tends to negatively impact fuel utilization (Petrovic and Maldonado, 2016). Online refueling has been proposed as a potential solution for achieving better fuel utilization (Avigni and Petrovic, 2018) and improving the economic viability of this reactor design. Online refueling offers several advantages in terms of fuel utilization but it also presents several challenges from the design, operation and safety points of view. Nonproliferation aspects could also represent a potential challenge, in any case limited by the large size of the fuel assembly, which is 6 m tall and cannot be easily diverted. In order to investigate the viability of this solution and assure that the reactor can be operated safely, this work analyzes the thermal hydraulic aspects of the online refueling.

The Advanced High Temperature Reactor (AHTR) is a Fluoride-Saltcooled High-temperature Reactor (FHR) developed by ORNL (Holcomb, 2011; Varma, et al., 2012). Table 1 presents a list of main design parameters for the reactor; the thermal power of the core is 3400 MWth and the nominal coolant flow rate is roughly 28400 kg/s. The core is constituted by 252 fuel assemblies 6 m tall and hexagonally shaped (Fig. 1, left). Each assembly contains 18 fuel planks (Fig. 1, middle and right); the moderator material is graphite and the fuel is based on TRISO particles. The baseline fuel enrichment was initially 19.75%, but lower values have been considered since in order to minimize fuel and enrichment costs and 9% has been selected as the new reference enrichment (Varma, et al., 2012), hence it is used in our studies. The coolant flows through a 7 mm wide gap that separates each pair of fuel plates.

⁎

1.2. Online refueling as solution for fuel utilization issues An intrinsic characteristic of the AHTR, deriving from its fuel design based on TRISO particles and graphite as moderator, is the tendency to low fuel utilization compared to other reactor designs. The core and

Corresponding author at: CERN, Route de Meyrin, 1211 Genève, Switzerland. E-mail address: [email protected] (P. Avigni).

https://doi.org/10.1016/j.nucengdes.2019.110440 Received 7 August 2019; Received in revised form 6 November 2019; Accepted 11 November 2019 0029-5493/ © 2019 Elsevier B.V. All rights reserved.

Nuclear Engineering and Design 358 (2020) 110440

P. Avigni and B. Petrovic

Nomenclature AHTR CAD CFD FHR

FLiBe NEUP ORNL SAS TRISO UDF

Advanced High Temperature Reactor Computer-Aided Design Computational Fluid Dynamics Fluoride-cooled High-temperature Reactor

• A similar procedure, but reversed, is used to position the new fuel

fuel design of the AHTR require a relatively high moderator-to-heavymetal volumetric ratio, resulting in a limited amount of volume available for the fuel in the core; combined with the constraints on enrichment (Holcomb, 2011; Varma, et al., 2012), the amount of fissile material loaded in the core is limited, resulting in a short cycle length and reduced fuel discharge burnup. Several improvements of the fuel design have been developed in order to address the fuel utilization issue (Gentry et al., 2017; Kingsbury and Petrovic, 2016; Huang and Petrovic, 2018); the online refueling represents a valuable option for addressing the issue, by eliminating the constraints relative to the low amount of fuel loaded in the core. A comprehensive analysis of the online refueling approach for the AHTR has been developed and presented in (Avigni and Petrovic, 2018). The analysis demonstrates that online refueling eliminates the issues related to cycle length, refueling outages and relative cost, providing better fuel moderation and utilization; these contributions result in a reduction of the overall fuel cycle cost (Avigni and Petrovic, 2018). Moreover, this approach is compatible with the need of performing the refueling at high temperature. The analysis presented in (Avigni and Petrovic, 2018) addresses the online refueling for the AHTR on several levels, including general thermal hydraulics, neutronics, and economics. This work presents the detailed thermal–hydraulic analysis and the analysis of the stability of the fuel assembly during the refueling operation.

assembly in the core.

A thermal hydraulic analysis of the online refueling has been performed in order to show that this procedure can be performed safely and the flow distribution and the temperature of the fuel assemblies can be maintained within appropriate operational conditions. In particular, the analysis involves the estimation of the thermal margin with respect to fuel failure temperature and the verification of flow stability during the operation. The analysis required the development of a novel methodology involving CFD analysis, system codes and reduced-order modeling of hydraulic aspects. The CFD evaluations have been performed using the ANSYS Fluent code (“ANSYS Fluent User’s Guide”, 2013) and the Fluent Meshing mode as well as ANSYS Meshing have been used for the model preprocessing and mesh generation. System modeling has been performed using the RELAP5-3D code (RELAP5-3D code development team, 2012), whereas 1D reduced-order modeling and coupling of simulations has been developed via MATLAB scripts. The ANSYS Mechanical package has also been used for transient thermal evaluations. The flow variations that characterize the online refueling transient affect the neutronics behavior and the flow and temperature distribution in the lower and upper plenum, thus requiring a modeling approach that spans from local scale to full system scale. The flow distribution also has an impact on the stability of the fuel assembly during its extraction from the core, mainly affected by the fact that the flow through the channel of the removed assembly strongly increases with respect to its nominal value, thus producing a non-negligible lift force that needs to be released via a dedicated support system. The work has been organized starting from the single channel level steady analysis, followed by the modeling of the assembly in steady conditions; this model has been used for developing the characterization of the time dependence of the coolant flow rate through the channel of the removed assembly; the results of this analysis have been integrated with the CFD modeling of the lower and upper plenum.

1.3. Methodology for the thermal hydraulic analysis of the online refueling The online refueling procedure is based on the replacement of a used fuel assembly from the core with unused fresh one as the reactor is operating. The procedure works as follows:

• The assembly lifting mechanism is positioned over the assembly to be replaced; • The control rod is inserted in the assembly to be replaced and the power is adjusted consequently; • The assembly lifting mechanism grasps the assembly and lifts it to the upper plenum, keeping it covered with coolant; • The assembly lifting mechanism moves the assembly horizontally to the refueling lobe, a cavity in the upper portion of the reactor vessel; • The assembly is extracted from the coolant and lifted to the argon atmosphere; • The assembly is lowered into a spent fuel pool;

2. Thermal-hydraulic analysis of the coolant channel and fuel plate A steady state thermal–hydraulic analysis has been performed at the coolant channel and fuel plate level in order to provide a basic characterization of the flow and temperature distributions. The analysis is based on the following assumptions:

Table 1 AHTR parameters list. Power

3400

MWth

Efficiency

45%

Cold leg temperature Hot leg temperature

650 700

°C °C

3 12.9

MW/m3

Mass flow rate

28,400

kg/s

52

W/cm3

Flow velocity Core diameter

2 7.81

m/s m

2.55 7

cm mm

Core height Number of fuel assemblies

6 252

m

Primary loops Volumetric core power density Power density in fueled region Fuel plate thickness Coolant channel thickness Fuel stripe thickness TRISO particle diameter

7 847

mm μm

Lithium Fluoride and Beryllium Fluoride molten salt Nuclear Energy University Program Oak Ridge National Laboratory Scale-Adaptive Simulation Tristructural-Isotopic particle User-Defined Function

• Axial conduction in the vertical direction of the fuel plate can be neglected; • The coolant (FLiBe) is incompressible with constant properties; • The channel is approximated as an infinite parallel plates configuration; • Entrance length effects are neglected. The effects of axial heat conduction have been evaluated via finite volume modeling, by comparison of two identical models, one including the axial conduction term and another not including it. The results show that the influence of axial conduction (in steady state conditions) affects the temperature distribution and its maximum values by less than 1%. This difference is acceptable, considering the 2

Nuclear Engineering and Design 358 (2020) 110440

P. Avigni and B. Petrovic

Fig. 1. Cross-sectional view of the AHTR core (left), fuel assembly (middle) and fuel plate (right).

uncertainties in the properties of the materials (particularly thermal conductivity). The properties of the salt are relatively well known for the AHTR operational conditions and reference formulas for the interpolation of experimental data are provided in (Richard et al., 2014). This fluid can be approximated as incompressible and most of the properties are not strongly dependent on the temperature. Table 2 provides the average value and uncertainty for specific properties in the operating temperature range of the AHTR (650 °C-700 °C). Friction factor and heat transfer coefficient correlations have been developed based on CFD (Fluent) simulations of a single coolant channel. The simulations also show that a 2D rectangular channel approximates the full 3D coolant channel with less than 1% error in both flow and temperature distribution. The friction and heat transfer correlations have been used to simulate the thermal hydraulic behavior of the coolant channel and fuel plate in nominal operating conditions. Fig. 2 shows the temperature and heat transfer coefficient profiles of the coolant along the active portion of the average coolant channel. The average heat transfer coefficient is about 7500 W/m2K, derived from a Reynolds number of about 7300, which corresponds to about 1.9 m/s mean flow velocity through the coolant channel. The temperature profile in the direction transversal to the plate is parabolic with maximum temperature at the interface between the graphite meat and each fuel stripe. The temperature profile is strongly affected by the thermal conductivity of graphite in irradiated conditions. The maximum temperature depends on several factors, mainly the fission density distribution and the burnup distribution. However, the thermal limit for TRISO fuel is notably high, maintaining the hottest fuel assembly far from unsafe conditions. The effect of irradiation of graphite on the fuel plate has been analyzed; the graphite considered for the analysis is H-451 grade graphite, developed for the HTGR program, having similar properties to the one planned to be used for the AHTR (Ingersoll et al., 2004). The conductivity of graphite strongly decreases with irradiation but it remains above 20 W/(m*K) even for high neutron fluence (Fig. 3). The maximum fuel temperature increases with the irradiation by less than 100 °C from non-irradiated to irradiated conditions, but the integrity of the fuel is not challenged, due to a large safety margin inherent to this type of fuel. A 2D heat conduction model has been developed to evaluate the temperature distribution in the average fuel plate of the assembly (Fig. 4). This distribution has been obtained using 675 °C for the coolant temperature of both the internal channels and the bypass channel, and 7500 W/(m2*K) for the heat transfer coefficient. The temperature

distribution in the central portion of the plate practically does not depend on the x location, so the maximum temperature can be evaluated, assuming no heat conduction in the x direction, by solving the 1D heat conduction equation along the y direction. 3. Thermal hydraulic characteristics of the fuel assembly The finite-volumes model developed for the fuel plate has been extended to the calculation of flow and temperature distribution of the fuel assembly for the simulation of the heat transfer and the localization of hot spots. Further analyses at the assembly level include the CFD study of the flow distribution at the outlet of the assembly and a simplified 2D transient calculation of the thermal behavior of the assembly during the removal. 3.1. Thermal hydraulic modeling of the fuel assembly A three-dimensional calculation at the fuel assembly level has been developed in MATLAB, using friction and heat transfer correlations for the characterization of the flow along the vertical direction and 2D heat conduction calculations for the characterization of the temperature distribution at each axial level. This model is based on the preliminary modeling presented in (Avigni and Petrovic, 2014a; Avigni et al., 2013; Avigni and Petrovic, 2014b). The main characteristics of the temperature distribution at the assembly level are shown in Fig. 5. Since the purpose of this section is to present the general characteristics of the average fuel assembly, the average radial power density is considered for this calculation. This assumption is conservative compared to considering the used fuel assembly; the removal of the assembly will be analyzed later. The assembly in nominal steady state conditions has a large margin to the fuel thermal limit; hot spots are localized in the first and last plate of each third of assembly; the maximum is not coincident with the mid Table 2 Summary of FLiBe coolant properties in the temperature range 650 °C-700 °C.

3

Property

Mean Value

Uncertainty

Deviation from mean

Mass Density [kg/m3] Dynamic viscosity [Pa*s] Thermal conductivity [W/(m*K)] Heat capacity [J/(kg/K)] Volumetric heat capacity [J/(m3*K)]

1950.0 0.006 1.1

0.05% 20% 15%

0.36% 6.08% 0.66%

2416.0 4.71e6

2% –

0% 0.36%

Nuclear Engineering and Design 358 (2020) 110440

Temperature [°C]

700 690 680 670 660 650

0

1

2 3 Elevation [m]

4

5

Heat transfer coeff. [W/(m²*K)]

P. Avigni and B. Petrovic

7800 7650 7500 7350 7200

0

1

2

3 Elevation [m]

4

5

Fig. 2. Temperature (left) and heat transfer coefficient (right) as function of elevation in the active portion of the coolant channel.

plane but it is not far from it. Also, the two external plates, the one next to the support and the one next to the box, appear to be 20 °C hotter than the remaining 4 plates. This temperature peaking effect would be even more important if a corrected, non-uniform power distribution on the horizontal plane (normal to the flow direction), derived from neutronic considerations, was implemented. With respect to the axial coolant temperature distribution, the assumed 0.5% bypass flow results in very similar profiles for internal and bypass channels; the two smaller channels are 20 to 25 °C hotter than the other channels. The 3D analysis of the fuel assembly suggests that the fuel is substantially far from the thermal limit and the design allows adequate cooling of the system. Minor geometrical adaptions, such as the increase in the dimension of the two smaller channels, may help to reduce the peaking factors.

the turbulent spectrum in unstable flow conditions by adjusting to the already resolved scales in a dynamic way. The boundary layer is treated via the enhanced wall treatment model (EWT-omega) available in Fluent (“ANSYS Fluent Theory Guide”, 2013). The simulation has been run in transient mode, using time steps of 0.02 s. Fig. 7 shows the vorticity contours for six consecutive time steps spaced 0.2 s apart. The solution is not steady and shows strong oscillatory interaction of the multiple jets exiting the top of the assembly. The non-uniformities in flow distribution caused by the jets are important in the initial part of the outlet volume, but gradually homogenize towards the upper boundary; the resulting flow at a distance of 30 cm from the assembly outlets (corresponding to the inlet of the upper plenum) is uniform, justifying the fact that the assembly channels do not need to be explicitly modeled in the upper plenum model, presented in section 5.

3.2. Flow at the outlet of a fuel assembly

3.3. Control rod insertion, assembly extraction and transfer to spent fuel pool

The flow at the outlet of the assembly shows more complicated characteristics from a turbulence standpoint, since each channel strongly interacts with the nearby parallel channels. The behavior of the flow at this location has been investigated using a 2D model. Fig. 6 shows the 2D model and the mesh used for the simulation. The model includes the outlets of the assembly channels and an outlet volume representing the perforation through the upper support plate (30 cm think), through which the coolant flows before reaching the upper plenum. The turbulence modeling approach is based on the Scale-Adaptive Simulation formulation by Menter, available in Fluent (“ANSYS Fluent Theory Guide”, 2013); this approach is an improved unsteady RANS formulation, based on k-omega SST model, allowing the resolution of

The assembly removal and transfer to spent fuel pool has been simulated using ANSYS Mechanical; a simplified 2D approximation of the assembly is modeled. The simplified fuel assembly cross section model neglects the presence of the cladding and assumes that 100% of the power is deposited in the fuel stripes. With respect to the internal channels, the reference temperature and heat transfer coefficient selected for the steady state conditions are 675 °C and 7500 W/m2K respectively. The power decreases exponentially. The calculation assumes that the assembly is removed immediately after the control rod is inserted; an alternative option, not considered in this work, would be to allow for some delay between control rod insertion and extraction, in

Fig. 3. Graphite thermal conductivity as a function of temperature for various irradiations. 4

Nuclear Engineering and Design 358 (2020) 110440

P. Avigni and B. Petrovic

Fig. 4. Temperature distribution of the average fuel plate.

order to reduce the residual decay power and thus the maximum assembly temperature during the removal procedure. The amount of coolant in the upper plenum (equivalent to roughly 7 m height) is sufficient to keep the fuel assembly covered in the coolant, as long as the assembly is kept in the upper plenum. Following the extraction to the upper plenum, the assembly is moved horizontally to the refueling lobe, where a lifting mechanism allows to move the assembly to the spent fuel pool; it is at this point that the assembly, from being fully immersed in the coolant, emerges to the argon atmosphere. During the removal transient the heat transfer coefficient decreases in time to simulate the flow reduction of the internal channels. The flow reduction, represented by the reduction of the heat transfer coefficient, is assumed to happen instantaneously, when the assembly starts to be lifted; this assumption is conservative in terms of peak fuel temperature: the actual evolution of the heat transfer coefficient, gradually decreasing as the assembly is lifted, would result in lower average temperature levels. The transient simulates the insertion of the control rod inside the assembly (60 s for a reduction of power density to roughly 60% of nominal), the extraction to the upper plenum and eventually the extraction to the argon gas buffer above the plenum for transfer to the spent fuel pool. Similarly to the assembly extraction from

the core, the reduction of heat transfer coefficient due to the transition to the Ar atmosphere is conservatively modeled as an instantaneous decrease. The values of heat transfer coefficient for the stay of the assembly in the upper plenum and in argon atmosphere have been estimated based on correlations for free convection on vertical wall. The heat transfer coefficient for the stay in the upper plenum has been conservatively set at 100 W/m2*K, and the heat transfer coefficient for the stay in argon atmosphere has been set to 10 W/m2*K. The analysis demonstrates that the thermal margin is sufficient to maintain the fuel temperatures below a critical level, as long as the assembly is kept covered in the coolant. The fuel temperature presents a peak roughly 10 min after the insertion of control rods, but the peak is lower than under the nominal operating conditions. Instead, if the fuel assembly is extracted to the upper plenum (argon gas atmosphere), the temperature of the fuel keeps rising. Fig. 8 shows the minimum and maximum fuel assembly temperature as function of time. The plot shows the initial temperature decrease due to the insertion of the control rod during the first 60 s, followed by the 1 min extraction of the assembly to the argon atmosphere in the upper plenum, highlighted by the increase in minimum temperature between 1 and 2 min. For times greater than 2 min, the temperatures steadily

Fig. 5. 3D temperature [°C] distribution for 1/3 assembly. 5

Nuclear Engineering and Design 358 (2020) 110440

P. Avigni and B. Petrovic

Fig. 6. Mesh details for the 2D model of the assembly outlet.

Fig. 7. 2D vorticity contours for three time steps, 0.2 s spaced.

increase, and the maximum temperature reaches 1000 °C in roughly an hour. Despite the large thermal margin, a limited time is available for transfer of the assembly from the core to the pool.

been conducted. The analysis is coupled with 1D calculations to correctly represent the boundary conditions of the simulated fluid volume. These 1D calculations are performed via user-defined functions that simulate the pressure drop along the internal and bypass channels of the assembly using basic pressure drop correlations, allowing to compute the flow splitting between the two types of channels (Fig. 10). A generic assembly in the core is considered for this analysis, assuming that the pressure drop characteristics for a channel from which the fuel assembly is being extracted are substantially independent from the

3.4. Flow in the channel of the removed assembly In order to analyze the transient removal of the assembly, a CFD analysis of the flow distribution in the volume between the lower support plate and the lower face of the removed assembly (Fig. 9) has 6

Nuclear Engineering and Design 358 (2020) 110440

P. Avigni and B. Petrovic

Fig. 8. Time evolution of maximum (orange) and minimum (blue) temperatures for an average assembly: insertion of control rod (1 min), extraction to upper plenum (1 min) and stay in argon atmosphere.

position in the core. The purpose of the analysis is the characterization of the coolant flow rate through the channel of the removed assembly as a function of the assembly elevation in the core, which is useful for evaluating the behavior of the core and the plena during the removal transient, as described in sections 4 and 5. It is assumed in the modeling approach that the removal of the assembly from the core can be treated as a sequence of steady state flow solutions, due to the relatively large time scale of the transient. As for the upper plenum simulations, moving and deforming mesh approaches have been investigated, but considered too computationally demanding. Assembly extraction simulations have been conducted using Fluent CFD simulation software. The model includes the entire zone between the upper surface of the lower support plate and the lower face of the extracted assembly (Fig. 9), as well as partial modeling of the outlets from the lower plate, the inlets of the assembly and the inlets of the bypass channels. Several simulations have been run using different heights to simulate the extraction of the assembly. The final mesh is polyhedral, and it is generated in Fluent solver mode starting from a tetrahedral mesh with boundary layers. On the solver side, the k-ε realizable model has been used for turbulence,

integrated with the enhanced wall treatment for the boundary layer. Y + is between 5 and 25 in the upper portion of the model, and approaches values lower than 1 in the lower portion, thus requiring a flexible modeling for the boundary layer; the enhanced wall treatment, used for this simulation, is capable of combining the two-layer model for fine near-wall mesh and enhanced wall function for coarse near-wall mesh. UDFs have been used to compute pressure and velocity values for the boundary conditions. The boundary conditions are velocity-inlet for the lower support plate channels and pressure-outlet for the assembly internal channels and the bypass channels. Fig. 11 shows the typical flow configuration in the channel of the extracted assembly at an elevation of 20 cm. The figure shows the pathlines colored by velocity magnitude. The pattern of the inlet channels is clearly visible, and it can be noticed that a large part of fluid tends to flow through the bypass channel for this specific elevation of the assembly. Parametric simulations have been run for a set of values of flow rate and elevation in order to correlate these input parameters to the pressure drop via an interpolated correlation. This correlation has two input parameters, the elevation of the assembly and the total mass flow rate through the channel, and two output parameters, the pressure drop across the channel of the removed assembly and the fraction of the total channel flow that flows through

Fig. 9. CAD model configuration for various time steps. As the time increases, the height changes to simulate the extraction of fuel assembly from the core, i.e., the figure shows the portion above the core. 7

Nuclear Engineering and Design 358 (2020) 110440

P. Avigni and B. Petrovic

drops across the channel of the removed assembly and the rest of the core are the same. The functional dependence between the channel mass flow rate and the elevation of the assembly in the core is shown in Fig. 12. The flow rate strongly increases from 110 to 270 kg/s during the initial 20 cm lifting, then the slope of the curve decreases, but the mass flow keeps increasing up to 360 kg/s. The cause of this strong increase in flow rate is the fact that, following the initial lifting of the assembly, the coolant is allowed to flow also through the bypass channels surrounding the assembly box, which present a friction coefficient that is lower than for the assembly channels. A strong increase also happens when the assembly is extracted from the channel, bringing the flow rate through the channel to about 620 kg/s. This is caused by the fact that, when the assembly is fully extracted to the upper plenum, it does not constitute an obstacle for the flow through the channel anymore.

Internal flow

Assembly

Assembly

Bypass flow

Assembly

Assembly

Assembly

Bypass flow

Lower Support Plate

4. Analysis of flow distribution in the lower plenum The calculated flow through the channel of the removed assembly has been used as an input boundary condition for the simulation of the behavior of the flow in the lower plenum. The objective of the simulations is to demonstrate that the assembly removal can be performed online in such a way that the flow distribution is not strongly altered during the transient. Due to the computational size of the problem, the simulation approach has initially been developed for a reduced-scale version of the problem, and it has been subsequently applied to the full-size version. The simulation approach is structured in two parts, a steady-state simulation and a transient simulation; the steady version of the problem is targeted at the understanding of the reference operational characteristics of the flow in the lower plenum, whereas the transient simulation aims at showing the changes that the removal of the assembly produces on these nominal steady conditions. A 3D model of the lower plenum has been prepared to perform the simulations; Fig. 13 shows a meshed version of it. The flow inlets for the model are three downcomer sections distributed at a 90° angle span around the circumferential perimeter of the vessel; the remaining fourth quadrant is not modeled because it is occupied by the maintenance cooling system. The outlets of the model consist of 253 hexagonal faces corresponding to the inlet boundary of each assembly. The mesh has been created via Fluent meshing based on the polyhedral meshing approach; Fig. 13 provides an illustration of the generated mesh for the full-size problem and Table 3 provides the characteristics of such mesh. The final mesh used for the simulation is on the order of 5 million cells and 11 GB memory requirement. The simulations are run on 40 parallel processes, with about 130 k cells per process.

Fig. 10. Schematic of flow configuration during assembly removal.

the bypass channels (complementary of the flow through the internal channels). The correlation has been developed by selecting specific values for the elevation, ranging from 2 mm (slightly lifted assembly) to 5.8 m (almost fully extracted assembly), and total channel flow rate, ranging from 110 kg/s (nominal flow) to 500 kg/s (expected upper bound on mass flow rate through the channel with assembly not completely extracted). At each elevation, the dependence of the pressure drop on the mass flow rate has been interpolated (represented) using the following power function:

P = A (E )

mflow B (E )

(1)

where ΔP is the channel pressure drop, mflow is the channel mass flow rate, A is the interpolation coefficient and B is the interpolation exponent and E is the assembly elevation. Both A and B are function of the assembly elevation E; the correlations for A(E) and B(E) have been derived from interpolation on the simulation data, keeping a constant flow rate. The correlation has been coupled to a simplified model of the rest of the core; the coupling assumes that the total mass flow rate through the core is fixed by design at 28400 kg/s and that the flow splits between the channel and the rest of the core in such a way that the pressure

Fig. 11. Pathlines colored by velocity magnitude [m/s] 8

Nuclear Engineering and Design 358 (2020) 110440

P. Avigni and B. Petrovic

Fig. 12. Channel mass flow rate as function of assembly extraction elevation.

In terms of turbulence modeling, in order to gain some flexibility in terms of y+, the k-ε realizable model has been used together with the two-layer wall treatment approach, which is a blend between the wallfunction and the full resolution of the viscous layer, making it less sensitive to y+ value. The y+ values for this model range between 15 and 82 in the most important regions of the model, thus resulting in the application of the enhanced wall function approach; the remaining parts of the model, for which y+ is lower, use a two-layer modeling approach. In nominal conditions the flow rate through the reactor is 28400 kg/ s at 650 °C inlet temperature, resulting in 2.7 m/s inlet velocity through

Table 3 Meshing and simulation parameters. Parameter

Value

Unit

Minimum mesh size Maximum mesh size Number of inflation layers Number of cells Mesh memory usage Number of cells per core Average runtime

5 40 3 5,085,520 10,731 129,000 30

mm mm – – MB #/core min

Fig. 13. 3D model (left) and mesh (right) for the flow simulation in the lower plenum. 9

Nuclear Engineering and Design 358 (2020) 110440

P. Avigni and B. Petrovic

the downcomer sections. A pressure outlet boundary condition is used on the outlet zones, with target mass flow rates calculated based on 1D correlations. Fig. 14 shows the pathlines colored by velocity magnitude in the lower plenum in nominal conditions. The peak velocity is located in the downcomers. The flow distribution is compartmentalized into three zones, corresponding to the three downcomer sections. The three zones are not symmetric, due to the asymmetric distribution of the downcomers along the vessel circumference. The consequence of this asymmetry is the fact that the reactor is divided into three main sections which behave somewhat independently in terms of flow and temperature; if for example the temperature increases in one of the mentioned sectors, only a subsection of the fuel assemblies is affected by this, with different consequences on the behavior of the core. The transient simulation has been run starting from the steady conditions. The extraction of the assembly is modeled by setting a time dependent target mass flow rate for the outlet boundary conditions of each assembly channel. The outlet mass flow rate is calculated via UDF as function of time based on the following steps:

Polyhedral mesh has been generated based on the parameters listed in Table 4. The mesh size is constituted of 25 million cells and the simulations have been run on 14 nodes and 112 parallel processes, resulting in 233,000 cells per node. The realizable k-epsilon model with enhanced wall treatment has been used for turbulence modeling; a steady simulation has been run to characterize the nominal conditions, followed by a transient simulation for the extraction of the assembly. The steady-state solution has been calculated assuming uniform mass flow rate distribution among the assemblies, corresponding to roughly 110 kg/s per assembly and a total of 28400 kg/s; uniform static pressure boundary condition has been modeled at the outlets. A parabolic profile in the radial direction has been assumed at the lower boundary for the incoming temperature distribution, with a reference outer radius of 8.45 m and min/max temperatures of 650 °C and 750 °C. Fig. 17 shows the pathlines colored by the temperature and the temperature distribution at three elevations in the core. The pathlines are directed vertically and bend by 90° towards the three outlets in the upper part of the plenum; the flow is distinctly compartmentalized in three sectors corresponding to the three cooling circuits. The presence of control rod guide tubes in the upper plenum does not seem to enhance mixing of fluid. Due to the flow compartmentalization and reduced mixing, it is expected that the sector (and associated cooling loop) that is associated with the assembly being removed will see a lower temperature, whereas the remaining two sectors will behave in a complementary way, thus receiving coolant at higher temperature. This consideration is based on the assumption that the power of the core is constant during the transient and the amount of flow circulated by each cooling loop remains unchanged. In terms of the temperature distribution in the upper plenum in steady condition, recirculation zones at low coolant temperature are located at the side boundaries of the plenum, and the planar temperature distribution remains relatively unchanged with elevation. The transient of the assembly removal has been simulated starting from steady conditions. An explicit modeling of the assembly removal requires the use of moving and deforming mesh approach, which has been considered too demanding for the available computational resources. For this reason, the assembly removal has been approximated by implementing a transient inlet flow-rate for the location of the removed assembly according to the results presented in section 3.4; moreover, the insertion of control rod is modeled by implementing a transient reduction of inlet temperature for the coolant at the inlet boundary. The advantage of this approximated approach is the fact that the mesh does not need to be recomputed at each time-step, which could be extremely demanding for such a large problem. The insertion of the control rod does not substantially alter the temperature distribution in the upper plenum. Fig. 18 shows the velocity magnitude and the temperature in the upper plenum at the end of

• The lifting velocity of the removed assembly is calculated as a

• • •

function of time. The time-profile of the velocity is linearly increasing from 0 to nominal and subsequently it is assumed to be constant until the end of the transient. The initial gradual increase of lifting velocity aims at avoiding sudden changes in flow conditions and improving stability. The nominal extraction velocity selected for the transient is 0.5 m/s; The elevation of the removed assembly is calculated according to the time, time step and current extraction velocity; the elevation increases quadratically at the beginning of the transient, and then linearly; The mass flow rate through the channel of the removed assembly is calculated based on the correlation presented in section 3.4; The mass flow rate through the remaining assemblies is calculated based on the total flow and the flow through the channel of the removed assembly.

The mass flow variation during the transient is substantial for the channel of the removed assembly, but the variation is very small for the remaining 252 assemblies. Fig. 15 shows the velocity magnitude contours at three elevations in the lower plenum, at the beginning and at the end of the transient. The plots show that the effects of the flow rate distribution at the assembly channels inlets during the assembly extraction are limited to the volume in proximity of the inlet of the extracted assembly. The levels at 1 m and 2 m below the assembly inlets are not strongly affected by the flow change. 5. Analysis of flow distribution in the upper plenum A CFD model of the upper plenum has been developed for the simulation of the flow and temperature distribution following the removal of a fuel assembly during operation, in order to demonstrate that the assembly extraction can be performed without strongly altering the operational conditions of the reactor or compromising its safety. The model, based on the geometry shown in Fig. 16, uses the conclusions from the analysis of the flow in the assembly channel as an inlet boundary condition. The model shown in Fig. 16 includes the following components:

• 253 hexagonal velocity-inlet boundaries, representing the outlets of the core assemblies; • 3 circular outlet faces, representing the 3 cooling loops of the reactor; • Control rod guide tubes, one for each assembly.

Fig. 14. Lower plenum velocity streamlines colored by velocity magnitude (front view). 10

Nuclear Engineering and Design 358 (2020) 110440

P. Avigni and B. Petrovic

Fig. 15. Horizontal velocity magnitude contours at different elevations before and after the assembly extraction.

Table 4 Meshing parameters for the upper plenum model. Parameter

Value

Units

Minimum mesh size Maximum mesh size Growth rate Number of boundary layers Aspect ratio of first layer Mesh size Minimum orthogonal quality Maximum aspect ratio Mesh memory usage

1.0 7.0 1.08 3 10 25,104,557 0.164 24.99 42,112

cm cm – – – #cells – – MB

6. Lift force on the assembly during extraction and viability of online refueling 6.1. Challenges to online refueling Thermal and hydraulic aspect of online refueling has been addressed in previous sections. However, another critical aspect for the evaluation of feasibility of online refueling is the evaluation of the structural stability handling viability of the assembly removal operation. The handling stability is substantially dependent on the overall forces acting on the fuel assembly as it is extracted that are caused by flow, buoyancy and handling equipment. The viability of the online refueling also requires a redesign of the upper support plate and the control rod guide tubes layout.

Fig. 16. 3D model of upper plenum, showing 253 hexagonal perforated faces (constituting the inlet boundaries of the model) and 3 circular outlets.

the transient; at this stage the alteration of the flow distribution is more consistent. The increased flow at the outlet of the channel of the removed assembly results in higher flow velocity and colder coolant temperature in the vertical direction along the plenum; when the coolant reaches the upper part of the plenum, flow mixing results in a homogenization of velocity and temperature.

6.2. Components of the lift force A basic 1D analysis based on the hydraulic analysis of the coolant channel (section 2) has been performed to evaluate the overall force acting on the assembly during removal. The analysis accounts for the fact that nuclear-grade graphite, which is the main component of the fuel assembly, is on average less dense than FLiBe, potentially causing the assembly to float if not properly counterbalanced. A second aspect 11

Nuclear Engineering and Design 358 (2020) 110440

P. Avigni and B. Petrovic

Fig. 17. Pathlines colored by temperature (left) and contours of temperature (right) [K] for the steady-state simulation.

Fig. 18. Upper plenum velocity magnitude [m/s] (left) and temperature [K] (right) after assembly extraction.

to be accounted for is the fact that the assembly is planned to be removed without stopping the flow, which can potentially provide a nonnegligible vertical push. The lift force is caused by three contributions:

two contributions, the drag due to internal channels and due to bypass channels; as shown in Fig. 19, it tends to zero as the assembly approached full extraction. Fig. 19 shows the components of the lift force as function of the assembly elevation in the channel; the peak force is at the beginning of the transient. The overall force acting on static assembly completely inserted in channel is equal to roughly 10000 N and directed upward. In absence of mechanical constraints, the assembly would be ejected from the core; this explains the need for an upper support plate.

1. Weight. The weight of the assembly is directed downwards and it is mainly based on the density of graphite (1900 kg/m3), combined with the weight of the control rod. The average density of the control rod is roughly 2050 kg/m3, corresponding to a total mass of 1530 kg; 2. Buoyancy. The buoyancy force is determined by the pressure difference between the inlet and the outlet of the assembly. The pressure difference arises from the contribution of the static column of coolant (constant with elevation) plus the pressure difference due to the friction pressure drop along the channel (dependent on elevation). As shown in Fig. 19, the effect of buoyancy is on the order of 2000 kg, depending on the elevation. 3. Drag. It is calculated based on friction correlations and it is made of

6.3. Means of reducing the lift force To ensure stable online extraction and replacement of a fuel assembly, lift force (at its maximum value during the operation) should be reduced to a null-force or a negative (downward oriented) force. A reduction of the lift force can be achieved via the following means: 12

Nuclear Engineering and Design 358 (2020) 110440

P. Avigni and B. Petrovic

refueling for the AHTR. The analysis provides insights into the thermal–hydraulic characteristics of the reactor at the fuel plank level and the assembly level. Furthermore, a methodology was developed to integrate the analysis at the assembly level with the full core behavior in order to prove that the online refueling operations do not compromise the safety of the reactor. A series of recommendations is provided to improve the design of the fuel assembly in such a way that the extraction of the fuel assembly from the core can be performed reliably in a stable manner. These recommendations include an increase of the assembly weight combined with an increased width of the coolant channels; this updated assembly design needs to be coupled with a reduction of mass flow rate and potentially with a slight reduction of core power during refueling, in order to maintain the reactor in the vicinity of its nominal operational parameters. The models developed in this work provide a methodology by which further, more detailed studies can be performed. The study indicates that the online refueling is a valid option for improving the fuel utilization and ultimately the economics of the AHTR. Acknowledgements

Fig. 19. Total force and its components acting on the assembly as a function of extraction elevation. All forces are directed upward, except for the weight component.

This work has been partially funded by the following US Department of Energy under the Nuclear Energy University Program (NEUP) projects: “Fuel and Core Design Options to Overcome the Heavy Metal Loading Limit and Improve Performance and Safety of Liquid Salt Cooled Reactors” and “Integrated Approach to Fluoride High Temperature Reactor (FHR) Technology and Licensing Challenges”.

1. Reduction of mass flow rate to roughly 25% of its nominal value. In this condition, the maximum overall force acting on the assembly is equal to 0. 2. Increase of the assembly weight. This could be achieved by oversizing the control rod, providing as well and increase of the amount of hafnium for each control rod, which is insufficient in the current design. An alternative method for adding weight to the assembly could potentially be the addition of beryllium oxide to the assembly structure, which would also improve the moderation capabilities. Addition of weight to the assembly alone is not sufficient for reaching autonomous stability of the assembly during refueling. 3. Reduction of drag by increasing the flow area. Increase of flow area has two counterindications: the reduction of heat transfer coefficient and the reduction of moderation capacity. The reduction of heat transfer coefficient is compensated by the large thermal margin intrinsic to fuel design, whereas the reduction of moderation is acceptable because of substantially higher fuel utilization provided by online refueling. In order to achieve a null lift force, the thickness of the fuel plate needs to be reduced from 25.5 mm to 16.5 mm, with constant fuel assembly size. Alternatively, coolant channels may be made wider, which would however increase the fuel assembly and core size. 4. Reduction of drag by temporarily reducing the mass flow rate through the core and allowing for a temporarily higher core outlet temperature to allow continued operation at full power.

References Holcomb, D.E., at al., 2011. Core and Refueling Design Studies for the Advanced High Temperature Reactor, ORNL/TM-2011/365, ORNL, September 2011. Varma, V.K., et al., 2012. AHTR Mechanical, Structural, and Neutronic Preconceptual Design, ORNL/TM-2012/320, ORNL, September 2012. Gentry, C.A., Maldonado, G.I., Chvala, O., Petrovic, B., 2017. Neutronic Evaluation of a Liquid Salt Cooled Reactor Assembly. Nuclear Science Engineering 187, 166–184. Kingsbury, C.W., Petrovic, B., 2016. A model for Assessing FHR Fuel Fabrication and Fuel Cycle Cost. In: Proc. ICAPP 2016, San Francisco, CA, April 17-20, 2016. Huang, M., Petrovic, B., 2018. Development of Methodology for Efficient Fuel Design Evaluation of the Advanced High Temperature Reactor. Annals of Nuclear Energy 121, 646–660. RELAP5-3D code development team, 2012. RELAP5-3D Code Manual Volume I: Code Structure, System Models and Solution Methods, INEEL-EXT-98-00834, INL, June 2012. “ANSYS Fluent User’s Guide”, Release 15.0, ANSYS, Inc., Canonsburg, PA, US, November 2013. “ANSYS Fluent Theory Guide”, Release 15.0, pp. 122-128, ANSYS, Inc., Canonsburg, PA, US, November 2013. Avigni, P., Petrovic, B., 2014a. Fuel element and full core thermal-hydraulic analysis of the AHTR for the evaluation of the LOFC transient. Ann. Nucl. Energy 64, 499–510. Avigni, P., Petrovic, B., 2014b. Preliminary Evaluation of the Fuel Assembly At-Power Removal Transient of the Advanced High Temperature Reactor. In: Proc. 2014 Intl. Congress on Advances in Nuclear Power Plants (ICAPP 2014), Charlotte, NC, April 69, 2014, Paper 14311, pp. 78-85, 2014. Avigni, P., Petrovic, B., Ricotti, M. E., 2013. Preliminary Thermal-Hydraulic Analysis of the AHTR Fuel Element. In: Proc. 15th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-15), Pisa, Italy, May 12-15, 2013. Petrovic, B., Maldonado, G.I., 2016. Final Report: Fuel and Core Design Options to Overcome the Heavy Metal Loading Limit and Improve Performance and Safety of Liquid Salt Cooled Reactors, SRC#00128483, Georgia Institute of Technology & University of Tennessee, Atlanta, Georgia, April 2016. Richard, J., Wang, D., Yoder, G., Carbajo, J., Williams, D., Forget, B., Forsberg, C., 2014. Implementation of Liquid Salt Working Fluids into TRACE. In: Proceedings of ICAPP 2014, paper 14214. Charlotte, USA, April 6-9, 2014. Avigni, P., Petrovic, B., 2018. On-line refueling for the Advanced High Temperature Reactor (AHTR). Nuclear Engineering and Design 340, 166–182. Ingersoll, D.T., et al., 2004. Status of Preconceptual Design of the Advanced HighTemperature Reactor (AHTR), ORNL/TM-204/104, ORNL, May 2004.

While none of the measures by itself cannot achieve the target without leading to other problems, a combination of the aforementioned solutions is capable of fully counterbalancing the lift force on the assembly during removal, i.e., making it a null-force or negative (downward oriented). A potential solution could be the addition of 200 kg weight to the assembly coupled with a doubled coolant channel thickness, which would not require any reduction of power or flow rate during the online refueling. 7. Conclusions This work has analyzed the thermal hydraulic aspects of the online

13