A MULTIPHYSICS SIMULATION SUITE FOR SODIUM COOLED FAST REACTORS

A simulation suite has been developed to model reactor power distribution and multiphysics feedback effects in Sodium-cooled Fast Reactors (SFRs). This suite is based on the Finite Element Method (FEM) and employs a general, unstructured mesh to solve the Simplified P3 (SP3) neutron transport equations. In the FEM implementation, two-dimensional triangular elements and three-dimensional wedge elements are selected. Wedge elements are selected for their natural description of hexagonal geometry common to fast reactors. Thermal feedback effects within fast reactors are modeled within the simulation suite. A thermal hydraulic model is developed, modeling both axial heat convection and radial heat conduction within fuel assemblies. A thermal expansion model is included and is demonstrated to significantly affect reactivity. This simulation suite has been employed to model the Advanced Burner Reactor (ABR) benchmark, specifically the MET-1000. It has been demonstrated that these models sufficiently describe the multiphysics feedback phenomena and can be used to estimate multiphysics reactivity feedback coefficients.


INTRODUCTION
Renewed interest in advanced nuclear power reactors, such as the Versatile Test Reactor (VTR) at the Idaho National Laboratory (INL), has encouraged enhanced modeling and simulation of fast nuclear power reactors. In pursuit of this interest, a simulation suite has been developed to simulate fast reactors with an emphasis on Sodium-cooled Fast Reactor (SFR) modeling [1]. Recent additions to this suite include the addition of the Simplified P 3 (SP 3 ) equations.
In addition to calculating the reactor neutron distribution using the SP 3 equations, multiphysics models are included to calculate reactivity feedback effects. Simplified thermal hydraulic and thermal expansion models are developed to simulate these physical effects and calculate how they affect reactivity. All of the multiphysics models are coupled in a power iteration method and use a common geometry mesh.

SIMPLIFIED P 3 NEUTRON TRANSPORT
Originally developed by Gelbard, the SP N transport method is a middle ground between diffusion and transport methods [2]. The foundation of the SP N method is to approximate the multi-dimensional transport equation using relationships from the one-dimensional P N equations. The equations that result from these assumptions are diffusion-like and computationally simple compared to the full transport equation. Initial demonstrations of SP N method were based on numerical results, but the method was later demonstrated to be theoretically consistent with the transport equation using asymptotic derivation [3,4].
The quantity φ 0,g (r) is introduced as a linear combination of the SP 3 moments for simplification. The associated Marshak-like * boundary conditions for zero incoming neutron flux are given as wheren is the outward normal unit vector and ∂Γ is the problem boundary. * Boundary conditions reduce to Marshak boundary conditions in one-dimensional geometries. The SP 3 equations are discretized and solved using an unstructured Finite Element Method (FEM) formulation. For two-dimensional geometries, triangular elements are selected. For threedimensional geometries, wedge elements are selected. Fast reactors typically have hexagonal-z geometry so wedge elements are a natural choice for this coordinate system. Reactor geometries are also typically described in lattices so the wedge element allows for easily "stacking" lattices on top of each other. The linear system of equations associated with the FEM is solved one energy group at-a-time using a Successive Over-Relaxation (SOR) solver with an adaptive over-relaxation parameter. Additional details for the FEM implementation are provided in Reference [1].
Eq. (5) and Eq. (6) represent an eigenvalue problem. In this work, these equations are solved using a power iteration method. The power iteration method can be shown to converge to the fundamental eigenmode [6]. However, in this work, multiphysics models from Section 3 and Section 4 are used to update cross sections periodically during the iteration procedure and therefore convergence is not guaranteed. Cross sections could be recalculated based on thermal hydraulic and thermal expansion updates every power iteration, but this was observed to be unnecessary. It has been observed that updating cross sections every ten iterations provided efficient and convergent behavior. For other work, this number is easily changed in user input. Such a nonlinear update procedure is commonly used in similar reactor simulation applications and no convergence problems have been observed.

THERMAL HYDRAULICS
Using the scalar flux to calculate the reactor power distribution, a simple thermal hydraulics model is developed to estimate material temperatures. These temperatures are used for the calculation of temperature dependent cross sections. A steady-state, one-dimensional, axial heat convection model is developed to calculate coolant temperatures within a flow channel. This model assumes no cross-flow between channels and ideal fluid mixing within the channel. Typical fast reactor assemblies are canned, so one-dimensional flow channels is a reasonable assumption. An example fast reactor hexagonal assembly cross-section view is shown in Fig. 1a. The radial heat conduction model is developed to calculate material temperatures within an assembly and assumes axial heat conduction is negligible at steady-state.
In the axial heat convection model, fluid properties of sodium are functions of temperature [7]. In the radial heat conduction model, sodium and HT9 stainless steel (reactor structural material) are assumed to have constant thermal conductivity as their thermal conductivities vary little over the operating temperature range of SFRs. The fuel is assumed to be metallic uranium with 10% zirconium by weight (i.e. U10Zr). Thermal conductivity within the fuel is modeled as a function of temperature [8]. A plot of the thermal conductivities versus temperature is shown in Fig. 2a.

Axial Convection Model
The axial steady-state coolant enthalpy within each channel is given by an energy balance as where h i,k+1/2 is the enthalpy at the upper elevation of the element, h in is the reactor inlet enthalpy, m i is the mass flow rate through channel i, and q i,k is the total heat generated in channel i at   elevation k. The element-average enthalpy can then be estimated as and the coolant temperature is calculated by an equation of state for sodium [7].

Radial Conduction Model
A one-dimensional radial heat conduction model is used to calculate the average material temperatures in the fuel rods in each flow channel and at each axial elevation. The general methodology is to first calculate the material surface temperatures and then average the temperatures within the materials. A figure of the radial geometry modeled is shown in Fig. 1b. Clad surface temperature is calculated using the coolant temperature, Newton's law of cooling, and the Subbotin-Ushakov convective heat transfer correlation [9,10]. The Subbotin-Ushakov correlation was selected because it has been demonstrated to be useful for sodium cooled reactors at operating conditions. The bond and fuel surface temperatures are simply calculated by using the steady-state heat conduction equation with constant thermal conductivity and no volumetric heat generation.

Cross Section Treatment
After obtaining average material temperatures, cross sections are calculated for each flow channel and axial elevation in the reactor. Temperature dependent cross section tables are generated by a procedure outlined in Section 5.1 and then cross sections are interpolated within the table. In the sodium regions (coolant and bond), the sodium number density is determined from the coolant mass density and the microscopic cross sections are linearly interpolated in temperature. In the cladding material, the macroscopic cross sections are linearly interpolated in temperature. In the fuel, the dominant effect on cross sections due to temperature is Doppler broadening which obeys square-root behavior. Therefore, the macroscopic cross section in the fuel are interpolated in the square-root of fuel temperature. Cross sections are updated periodically during the solution of the SP 3 equations as described in Section 2.

THERMAL EXPANSION
In metal fueled reactors, thermal expansion represents a significant feedback effect. In the Experimental Breeder Reactor II (EBR-II) demonstrations, operators forced the reactor to undergo Unprotected Loss-Of-Flow (ULOF) and Unprotected Loss-Of-Heat-Sink (ULOHS) events. EBR-II shutdown safely solely due to multiphysics feedback effects. These experiments demonstrated the passive safety of fast reactor designs, due in part to the thermal expansion of materials [11].

Material Properties
For the purposes of this simulation, all structural materials in the reactor are thermally expanded as HT9 stainless steel [12]. Fuel material is expanded as U10Zr [13]. The Linear Expansion Factor (LEF) for these materials is plotted in Fig. 2b. The LEF of U10Zr is as much as twice that of HT9 which introduces the possibility of thermal expansion induced contact between fuel and cladding, but this has not been observed in reactor simulations. All sodium in the reactor is assumed to be liquid and effects of thermal expansion within sodium are described by the change in density [7].

Thermal Expansion Model Details
A simplified thermal expansion model is developed. All dimensions within the reactor are expanded from room-temperature according to two user input temperatures, T exp,fuel and T exp,struct corresponding to the average temperatures of fuel and structural material respectively. This is acceptable because the LEF are on the order 10 −2 and the reactivity effect of expanding from room-temperature to operating temperatures is more significant than small, local variations in material temperatures. In the radial (x and y) directions, it is assumed that the reactor expands as HT9 stainless steel. This is because the dominant reactivity effect due to thermal expansion in the radial directions is the expansion of the hexagonal assemblies themselves and the expansion of the grid plate at the base of the reactor, all of which are composed of HT9 stainless steel. In the axial (z) direction, it is assumed that the reactor expands as U10Zr. This is because the dominant reactivity effect due to thermal expansion in the axial direction is the elongation of the metallic fuel [14]. The general thermal expansion scheme is shown in Fig. 3

Figure 3: Thermal Expansion of a General Volume.
The coordinates of the finite elements expand as where x C represents the "cold" coordinate and x H represents the "hot" coordinate. Consistent with this thermal expansion model, the number density, N i , of species i in region j is adjusted as . (15) where a j is the volume/area fraction of region j and is calculated based on material dimensions. Macroscopic cross sections are then updated using number density ratio from Eq. (15) as where Σ H x is the "hot" macroscopic cross section for reaction x.

RESULTS
Detailed verification results of the individual physics models in the simulation suite are provided in Reference [1]. The original simulation suite employed a neutron diffusion model and multiple benchmark problems and mesh convergence rate verifications are provided. The newly implemented SP 3 model has been compared to diffusion and transport results and is undergoing additional verification. Thermal hydraulic and thermal expansion models were verified independently.

Advanced Burner Reactor -MET-1000
To demonstrate the multiphysics effects and coupling of the models, a realistic reactor benchmark is simulated and reactivity coefficients are calculated. The Advanced Burner Reactor (ABR) is a benchmark reactor design proposed by the OECD-NEA [15]. Specifically, the 1, 000 [MWth], metallic-fueled MET-1000 design is selected. Cross sections are not specified in the benchmark and are generated independently for each submission. The benchmark specification provides coolant number densities, but in this work the densities are calculated with the internal thermal hydraulic models. In addition, this work considers thermal expansion which is not included in the original benchmark.
33-group cross sections are generated using MC 2 [16] and tabulated at temperatures listed in Table 1.
Cross section temperatures were chosen to cover the operating range of the reactor at steady-state conditions. The reactor geometry is shown in Fig. 4a. The refined model contains 242,560 elements in 160 axial levels. For visualization † , the 33-group cross sections from MC 2 are collapsed into two-group cross sections and the fast and thermal fluxes from the two-group calculation are shown in Fig. 4b and Fig. 4c. Fast flux is shown to peak in the center of the core in the active fuel region. Thermal flux is shown to peak in core structural material, at the periphery of the active fuel region, as well as in control rod locations where the control rods have been withdrawn.
The multiphysics ABR solution had a runtime of 175 [s] on a single processor with less than 1.0 [GB] of memory required. Comparatively, a diffusion calculation of the same reactor without multiphysics coupling had a runtime of 66 [s] using DIF3D. A full three-dimensional transport calculation using a state-of-the-art method such as Method of Characteristics (MOC) may be expected to run for hundreds of processor-hours and require hundreds of gigabytes of memory. These runtime and resource requirements demonstrate the practical utility of the simulation suite and show that the SP 3 method sits in the middle ground between diffusion and transport solutions.

Reactivity Coefficients
In this work, reactivity is defined as A reactivity coefficient can be defined as a partial derivative with respect to some quantity of interest [17]. Let α x be the reactivity coefficient for quantity x, then where α x (x i ) has been approximated using the finite-difference method.
The multiphysics capabilities are used to calculate the power reactivity coefficient, thermal expansion reactivity coefficient, fuel temperature coefficient (Doppler coefficient), and the Coolant Temperature Coefficient (CTC). Reactivity coefficients are calculated using Eq. (18) and an appropriate choice of Δx. The power reactivity coefficient, α power is calculated by increasing reactor power by ΔQ Rx = 5% of nominal full-power. As thermal expansion is modified by two temperatures, T exp,struct and T exp,fuel , both are increased as a function of reactor power such that the temperature increase corresponds to a 5% increase in reactor power. Doppler reactivity coefficient is calculated by uniformly increasing fuel temperature by 5 [K]. Coolant Temperature Coefficient (CTC) is calculated by uniformly increasing coolant temperature by 5 [K] and this calculation also includes sodium density effects.
Results of the multiphysics reactivity coefficients are presented in Fig. 5. All coefficients are negative, except for the CTC which is slightly positive as expected. However, the net effect is represented by α power as shown in Fig. 5a and α power remains negative, implying the reactor is stable. As expected, the thermal expansion reactivity coefficient shown in Fig. 5b is shown to dominate the power reactivity coefficient.  Multiphysics effects are summarized in Table 2. To separate the contributions of thermal hydraulics and thermal expansion, the reactor was simulated using only the thermal expansion model and then using both thermal expansion and thermal hydraulics models. The data in Table 2 Fig. 5c and Fig. 5d respectively, have similar magnitude and opposite signs. For the ABR design, these coefficients are similar in magnitude but this may not be true for other reactors and thermal hydraulics effects remain important for simulation of an SFR.

CONCLUSIONS AND FUTURE RESEARCH
A multiphysics simulation suite based on the FEM has been developed to simulate SFRs. Physics models are included for SP 3 neutron transport, thermal hydraulics, and thermal expansion calculations. The novel feature of this suite is that all models are inherently coupled and allow for the calculation of multiphysics effects such as those demonstrated in Section 5. Such a suite will be useful for reactor design calculations and future optimization studies.
Future work will include adding capabilities to the simulation suite. Work is ongoing to incorporate a fuel depletion model employing a Chebyshev Rational Approximation Method (CRAM) solver [18]. The incorporation of depletion capabilities will allow for modeling SFRs for full operating cycles and allow for fuel design and loading pattern optimization.
Additionally, the use of higher-order scattering moments is being considered. Typically, isotropic scattering does not adequately describe the neutron distribution in fast reactor materials and higherorder Legendre scattering moments are necessary [19].