720p Movies | Kodi/XBMC and other plugins | Anti Vírus

Neutronics modeling of the High Flux Isotope Reactor using COMSOL

Neutronics modeling of the High Flux Isotope Reactor using COMSOL

Annals of Nuclear Energy 38 (2011) 2594–2605 Contents lists available at ScienceDirect Annals of Nuclear Energy journal homepage: www.elsevier.com/l...

3MB Sizes 1 Downloads 189 Views

Annals of Nuclear Energy 38 (2011) 2594–2605

Contents lists available at ScienceDirect

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

Neutronics modeling of the High Flux Isotope Reactor using COMSOL David Chandler a,⇑, G. Ivan Maldonado a, R.T. Primm III b,1, J.D. Freels b a b

University of Tennessee, Department of Nuclear Engineering, 311 Pasqua Engineering, Knoxville, TN 37996-2300, USA Oak Ridge National Laboratory, Research Reactors Division, P.O. Box 2008, Oak Ridge, TN 37831-6399, USA

a r t i c l e

i n f o

Article history: Received 21 February 2011 Accepted 3 June 2011 Available online 2 August 2011 Keywords: HFIR COMSOL Multiphysics Diffusion theory Neutronics SCALE NEWT

a b s t r a c t The High Flux Isotope Reactor located at the Oak Ridge National Laboratory is a versatile 85 MWth research reactor with cold and thermal neutron scattering, materials irradiation, isotope production, and neutron activation analysis capabilities. HFIR staff members are currently in the process of updating the thermal hydraulic and reactor transient modeling methodologies. COMSOL Multiphysics has been adopted for the thermal hydraulic analyses and has proven to be a powerful finite-element-based simulation tool for solving multiple physics-based systems of partial and ordinary differential equations. Modeling reactor transients is a challenging task because of the coupling of neutronics, heat transfer, and hydrodynamics. This paper presents a preliminary COMSOL-based neutronics study performed by creating a two-dimensional, two-group, diffusion neutronics model of HFIR to study the spatially-dependent, beginning-of-cycle fast and thermal neutron fluxes. The 238-group ENDF/B-VII neutron cross section library and NEWT, a twodimensional, discrete-ordinates neutron transport code within the SCALE 6 code package, were used to calculate the two-group neutron cross sections required to solve the diffusion equations. The two-group diffusion equations were implemented in the COMSOL coefficient form PDE application mode and were solved via eigenvalue analysis using a direct (PARDISO) linear system solver. A COMSOL-provided adaptive mesh refinement algorithm was used to increase the number of elements in areas of largest numerical error to increase the accuracy of the solution. The flux distributions calculated by means of COMSOL/SCALE compare well with those calculated with benchmarked three-dimensional MCNP and KENO models, a necessary first step along the path to implementing two- and three-dimensional models of HFIR in COMSOL for the purpose of studying the spatial dependence of transient-induced behavior in the reactor core. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Neutron diffusion theory is one of the simplest and most widely used methods to determine the neutron distribution within a reactor and can be used to characterize as many neutron energy groups as the user desires (Stacey, 2001). The two-group, spatially-dependent neutron diffusion equations were implemented in COMSOL Multiphysics v3.5a (COMSOL, 2008) to simulate neutron transport in a two-dimensional model of the High Flux Isotope Reactor (HFIR) in order to determine the thermal and fast neutron distributions within the reactor at steady-state beginning-of-cycle (BOC) conditions. This task was the first step to accomplishing a time-dependent neutronics solution in COMSOL using a twodimensional and then a three-dimensional model of HFIR. The purpose of the model development is to study the spatial dependence of transient-induced behavior in the reactor core. A COMSOL-based ⇑ Corresponding author. Tel.: +1 865 974 7562; fax: +1 865 974 0668. E-mail addresses: [email protected], [email protected] (D. Chandler), [email protected] (G.I. Maldonado), [email protected] (R.T. Primm), [email protected] (J.D. Freels). 1 Present address: Primm Consulting, LLC, USA. 0306-4549/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2011.06.002

thermal hydraulic and structural analysis model of the HFIR is under development in an independent but parallel project (Primm, 2011). That model is expected to eventually be merged with the current work to form a comprehensive multiphysics solution for transient-induced behavior. The cross sections needed to solve the diffusion equations were calculated by means of NEWT (DeHart, 2009), a two-dimensional neutron transport code in the SCALE 6 package (SCALE, 2009). The same geometry used in NEWT was created in COMSOL and the cross sections calculated by NEWT were inserted into COMSOL. The diffusion equations and associated boundary conditions were coded into COMSOL by means of the partial differential equation (PDE) coefficient application mode and the flux profiles in the HFIR core were solved via eigenvalue analysis and a direct (PARDISO) linear system solver. A COMSOL-provided adaptive mesh refinement algorithm was used to solve the diffusion equations using a sequence of refined meshes by increasing the number of elements in areas where the previous calculation (same PDE problem, but different mesh) yielded the largest numerical errors. Similar diffusion analyses have been performed for a molten salt breeder reactor (MSBR) core channel (Memoli et al., 2009) and a CANDU lattice (Gomes, 2008).

D. Chandler et al. / Annals of Nuclear Energy 38 (2011) 2594–2605

2595

1.1. High Flux Isotope Reactor description HFIR is a versatile research reactor located at the Oak Ridge National Laboratory (ORNL). HFIR was constructed in the mid-1960s for the purpose of producing heavy (transuranic) isotopes like 252 Cf. Today, the steady-state neutron fluxes produced in this 85 MWth reactor are utilized for cold and thermal neutron scattering, materials irradiation, isotope production, and neutron activation analysis. HFIR is a pressurized, light water-cooled and -moderated reactor that was designed with an over-moderated flux-trap (FT) in the center of the core and a large beryllium reflector on the outside of the core in order to produce a large thermal flux to power ratio for the purpose of transuranic isotope production. The central FT is surrounded by two concentric fuel annuli containing highly enriched uranium (HEU) in aluminum clad and water coolant channels. On the outside of the fuel elements (FE) are two concentric poison bearing control elements (CE), a large beryllium reflector, and light water. A mockup of HFIR is presented in Fig. 1. The FT consists of 37 target rod locations that accommodate materials to be activated. The FE consists of an inner fuel element (IFE) and an outer fuel element (OFE), each constructed of aluminum involute plates (171 in IFE and 369 in OFE) containing uranium enriched to approximately 93 weight percent in 235U/U in the form of U3O8 in an aluminum matrix. The total loading of a fresh HFIR core is about 9.4 kg of 235U and a typical fuel cycle length ranges between 22 and 26 days depending on the experiment loading. The CEs are located between the fuel and the reflector and each consist of three sections: a black region (Eu2O3–Al), a grey region (Ta–Al), and a white region (Al), and are named based on their neutron absorbing capability. 1.2. Computer code description NEWT is a multi-group discrete-ordinates code that is run within SCALE, a code package developed and maintained at ORNL. NEWT performs two-dimensional (2-D) neutron transport calculations and utilizes the Extended Step Characteristic (ESC) approach for spatial discretization on an arbitrary mesh structure. The primary function of NEWT is to calculate the spatial flux distributions within a nuclear system and collapse the cross sections into multiple (or single) energy groups as specified by the user. These collapsed cross sections can be supplied to ORIGEN-S for depletion

Fig. 1. Mockup of the High Flux Isotope Reactor.

Fig. 2. Schematic representation of NEWT and COMSOL models.

calculations, or in the case of this study, can be extracted from the output for use in another application (DeHart, 2009). COMSOL Multiphysics is a software package that uses the finite element method for spatial discretization to solve physics-based systems of PDEs and/or ODEs. Steady-state and time-dependent multiphysics simulations can be set up using the predefined (i.e. built-in) physics/engineering modules or by specifying a system of user-specific PDEs. Additional built-in modules can be modified to the user’s needs through equation based modeling capabilities and can be coupled together with other modules, with user defined PDEs, or with external coding through a MATLAB interface, all of which makes COMSOL a versatile simulation tool.

2. NEWT model development A two-dimensional NEWT model of HFIR was developed by modifying an existing NEWT model of HFIR that was created for low enriched uranium (LEU) conversion studies (Primm, 2009). The major modifications to the LEU model included changing the fuel to HEU, modeling the CEs in the control region, and by modeling the bottom half of HFIR rather than using symmetry across the core horizontal midplane (y = 0). The geometry utilized in the LEU model is a two-dimensional quarter configuration of HFIR such that symmetry was utilized at the core horizontal midplane and the core vertical centerline (x = 0, cartesian geometry). The HEU model developed for these studies utilized the same radial and axial boundaries and atomic densities as used in benchmarked threedimensional HEU TRITON (Chandler and Primm, 2009) and MCNP (Primm, 2005) models. The HEU input models a half configuration of HFIR such that symmetry was only utilized across the core centerline (x = 0, cartesian geometry). The flux trap is modeled as multiple homogenized regions in order to incorporate aluminum cladding, targets, and structure, water coolant, and curium target rods. The IFE and OFE are modeled as 8 and 9 radial regions, respectively, in order to incorporate the non-uniform distribution of HEU along the arc of the involute fuel plates. The non-fueled upper and lower regions of the FEs are modeled by homogenizing the water channels and the aluminum plates while the plate that separates the FEs is a mixture of aluminum and water coolant.

2596

D. Chandler et al. / Annals of Nuclear Energy 38 (2011) 2594–2605

Table 1 Geometric and material descriptions of model regions. Region description

Material description

Inner radius (cm)

Outer Radius (cm)

Axial bottom (cm)

Axial top (cm)

Target structure 1 Target fuel Target structure 2 Target structure 3 IFE 1 IFE 2 IFE 3 IFE 4 IFE 5 IFE 6 IFE 7 IFE 8 OFE 1 OFE 2 OFE 3 OFE 4 OFE 5 OFE 6 OFE 7 OFE 8 OFE 9 FE side plates FE extensions FE middle plate White CE Grey CE Black CE RBE PBE Water

Al, water Al, water, Cm target Al, water Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water HEU, Al, water Al, water Al, water Al, water Al, water Ta, Al, water Eu2O3, Al, water Be, Al, water Be, Al, water Water

0.000 0.887 3.547 0.000 7.14 7.50 8.00 8.50 9.50 10.50 11.50 12.00 15.13 15.50 16.00 16.50 17.50 18.50 19.50 20.00 20.50 6.40 (20.98) 7.14 (15.13) 12.60 22.02434 (22.987) 22.02434 (22.987) 22.02434 (22.987) 23.8125 33.3375 0.00

0.887 3.547 6.400 3.547 7.50 8.00 8.50 9.50 10.50 11.50 12.00 12.60 15.50 16.00 16.50 17.50 18.50 19.50 20.00 20.50 20.98 7.14 (21.7475) 12.6 (20.98) 15.13 22.65934 (23.622) 22.65934 (23.622) 22.65934 (23.622) 33.3375 54.61 84.61

25.40 25.40 25.40 30.48 (25.40) 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 30.48 30.48 (25.40) 30.48 0.00 (50.48) 12.70 (0.00) 50.48 (12.70) 30.48 30.48 50.48

25.40 25.40 25.40 25.40 (30.48) 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 25.40 30.48 25.40 (30.48) 30.48 50.48 (0.00) 0.00 (12.70) 12.70 (50.48) 30.48 30.48 50.48

The white, grey, and black regions of the CEs are modeled by homogenizing them with the clad and water gap regions. The CEs are inserted such that the faces of the gray regions are posi-

tioned at the core horizontal midplane as shown in Fig. 2. The white regions of the CEs are positioned above and below the ICE and OCE grey regions, respectively, and the black regions of the CEs are positioned below and above the ICE and OCE grey regions, respectively. The beryllium reflector is modeled as two separate regions, a removable beryllium reflector (RB) and a permanent beryllium reflector (PB), each composed of beryllium, water, and aluminum. The radial and axial boundaries of the regions defined in the NEWT geometry are defined in Table 1 with brief material descriptions of each region. The 238-group ENDF/B-VII neutron cross section library and NEWT were used to generate the two-group macroscopic cross sections needed for the COMSOL input. The 238-group cross sections and fluxes were collapsed into two-group form: fast flux (group 1) consisting of neutrons with energies between 3 eV and 20 MeV and thermal flux (group 2) consisting of neutrons with energies between 105 eV and 3 eV. A cutoff energy of 3 eV was

Fig. 3. Grid structure and material placement in NEWT model.

Fig. 4. Diffusion theory estimate of the extrapolation distance.

D. Chandler et al. / Annals of Nuclear Energy 38 (2011) 2594–2605

Fig. 5. COMSOL geometry drawing of HFIR.

2597

Fig. 6. NEWT half core thermal neutron flux distribution.

chosen such that upscattering (scattering from the thermal energy group to the fast energy group) could be neglected. The two-group macroscopic cross sections required to solve the diffusion theory equations include the transport cross section (Rtr), the absorption cross section (Ra), the average number of neutrons emitted per fission event times the fission cross section (vRf), and the downscatter (fast ? thermal) cross section (Rs1 ? 2). Each mixture’s two-group macroscopic cross sections were printed in the output through the utilization of the homogenization block in NEWT. After the collapsed energy structure defined in the previous paragraph was developed, the homogenized cross sections were created by combining the flux weighted collapsed cross sections with the number densities and added such that reaction rates in the homogenized materials were conserved (DeHart, 2009). The spatial and eigenvalue convergence criteria were both set to 104 and coarse-mesh finite-difference acceleration was activated to speed convergence since there is a significant amount of scattering in the system (DeHart, 2009). Two options are available in NEWT to set the regions within which convergence testing is applied: (1) force converged scalar fluxes in every computational cell and (2) relax convergence such that averaged scalar fluxes within a mixture are converged (DeHart, 2009). The second option is useful for mixtures in which fluxes become very small (large reflectors or near a vacuum boundary condition) and since this model utilizes three vacuum boundary conditions and includes a large beryllium reflector and water surrounding the core, the second option was utilized. The NEWT geometry, material number assignments, and fine rectangular mesh are depicted in Fig. 3. 3. Derivation of diffusion theory Fig. 7. NEWT half core fast neutron flux distribution.

The derivation of diffusion theory is based on Fick’s law, which governs that there will be a net flow of neutrons in a reactor from a region of greater neutron density into a region of lower neutron

density (Lamarsh and Baratta, 2001), and the equation of continuity, which governs that the net number of neutrons in a nuclear

2598

D. Chandler et al. / Annals of Nuclear Energy 38 (2011) 2594–2605

Fig. 8. COMSOL half core thermal neutron flux distribution.

system must be conserved. The expression for Fick’s law is shown in Eq. (1) and the expression for neutron continuity is shown in Eq. (2). For a complete derivation of these relationships refer to Stacey (2001) and Lamarsh and Baratta (2001).

J ¼ Dr/

ð1Þ

where J is the current density vector, D is  neutron   the diffusion  d d d coefficient ¼ 3R1tr , r is the gradient operator ¼ dx þ dy þ dz in rectangular coordinates, / is the neutron flux.

Z

dn dV ¼ dt

Z

SdV 

Z

Ra /dV 

Z

r  JdV

ð2Þ

The left-hand side (LHS) of Eq. (2) represents the time rate of change of the number of neutrons in volume V. The production rate in V, absorption rate in V, and the net leakage from the surfaces of V are shown from left to right on the right-hand side (RHS) of Eq. (2). The diffusion equation is developed by substituting Fick’s law into the equation of neutron continuity and is shown in Eq. (3)

dn ¼ S  Ra /  Dr2 / dt

ð3Þ

The focus of this study is to obtain the steady-state BOC fluxes and thus the time-dependence term can be neglected. Also, no external sources are present in HFIR and therefore fission is the only contributor to the production rate, S. Thus, the steady-state one-group diffusion equation with no external sources is written in the form of

v Rf /  Ra / þ D r 2 / ¼ 0

ð4Þ

Since two energy groups are being studied in this analysis, scattering from one energy group to another must be included into the diffusion equation. The two-group neutron diffusion equations for fast (group 1: noted with subscript 1) and thermal (group 2: noted with subscript 2) fluxes assuming all fission neutrons are born as fast neutrons are shown in Eqs. (5) and (6), respectively. The equations were rearranged such that the LHS of the equations describe the neutron loss mechanism and the RHS of the equations describe the neutron production mechanism. The effective multiplication factor, keff, was also inserted to balance the equations and describes how the population of neutrons varies from one generation to another.

D1 r2 /1 þ ðRa1 þ R1!2 Þ/1 ¼ s

1 ðv Rf1 /1 þ v Rf2 /2 Þ þ R2!1 /2 s keff ð5Þ

D. Chandler et al. / Annals of Nuclear Energy 38 (2011) 2594–2605

2599

Fig. 9. COMSOL half core fast neutron flux distribution.

D2 r2 /2 þ ðRa2 þ R2!1 Þ/2 ¼ R1!2 /1 s s

ð6Þ

2!1 s

By neglecting upscatter (R ¼ 0) and simplifying Eqs. (5) and (6) based on the problem definition of steady-state, BOC two-group analysis of HFIR, the following equations are applicable for the fast flux in the multiplying regions (Eq. (7)), the fast flux in the nonmultiplying regions (Eq. (8)), and the thermal flux in all the regions (Eq. (9)). 2

D1 r /1 þ ðRa1 þ R

1!2 Þ/1 s

1 ¼ ðv Rf1 /1 þ v Rf2 /2 Þ keff

surfaces (the top, bottom, and outer edge boundaries) where it is assumed that the neutrons that pass through these boundaries will not reenter the system. The vacuum BC is described pictorially in Fig. 4 and mathematically in Eq. (11). Interface BCs are used to show that the flux is continuous across the boundary interface between two different media, A and B, and is used at all of the interface boundaries (Eq. (12)).

Dr2 /jboundary ¼ 0

ð10Þ

/boundary ¼ ð2:1312ÞðDÞjr/jboundary

ð11Þ

/A jinterface ¼ /B jinterface

ð12Þ

ð7Þ

Þ/1 ¼ 0 D1 r2 /1 þ ðRa1 þ R1!2 s

ð8Þ

D2 r2 /2 þ Ra2/2 ¼ R1!2 /1 s

ð9Þ

Symmetry boundary conditions (BC), vacuum BCs, and continuity BCs are to be applied to the COMSOL model. A symmetrical BC is governed by Eq. (10) and shows that the divergence of the gradient of the neutron flux is equal to zero at the boundary. This BC is used at the core centerline since only half of the reactor is modeled. A vacuum BC uses the flux slope at the boundary to extrapolate the flux outside of the physical boundary a distance, d, where the flux vanishes to zero and assumes that no neutrons are reflected back through the boundary. The vacuum BC is used at the three pool

4. COMSOL model development The first step in developing a COMSOL model is to select the application mode that is best for the problem being solved. Since COMSOL does not have a built-in neutronics application-mode module, the PDE coefficient application mode and eigenvalue analysis is chosen so that the diffusion equations described in the previous paragraphs can be implemented. The COMSOL model was developed by using the graphical user interface (GUI) of the COMSOL client to:

2600

D. Chandler et al. / Annals of Nuclear Energy 38 (2011) 2594–2605

r  ðcru  au þ cÞ þ au þ b  ru ¼ da ðk  k0 Þu  ea ðk  k0 Þ2 u þ f

ð13Þ

where c is the diffusion coefficient, u is the dependent variable, a is the conservative flux convection coefficient, c is the conservative flux source term, a is the absorption coefficient, b is the convection coefficient, da is the damping/mass coefficient, k is the eigenvalue, k0 is the linearization point for the eigenvalue, ea is the mass coefficient, and f is the source term.

r  ðD1 r/1 Þ þ ðRa1 þ R1!2 Þ/1 s ¼ k/1 þ

1 ðv Rf1 /1 þ v Rf2 /2 Þ keff

ð14Þ

r  ðD1 r/1 Þ þ ðRa1 þ R1!2 Þ/1 ¼ k/1 s

ð15Þ

r  ðD2 r/2 Þ þ Ra2 /2 ¼ k/2 þ R1!2 /1 s

ð16Þ

Two PDE coefficient modes are coupled together and dependent upon each other since the thermal and fast neutron fluxes are being solved. The dependent variable for the first mode is the fast flux and the dependent variable for the second mode is the thermal flux. Thus, Eqs. (14) and (15) are input for the multiplying and nonmultiplying regions, respectively, in the fast flux PDE coefficient mode and Eq. (16) is input for the multiplying and non-multiplying regions in the thermal flux PDE coefficient mode. Since there are 31 unique materials in the model there are 31 distinct equations specified for the fast flux and 31 distinct equations for the thermal flux (62 distinct sets of cross sections and diffusion coefficients). The boundary conditions were defined in the boundary settings GUI. Boundary conditions were defined for both the fast and thermal flux physics modes. Symmetry was used at the core centerline, vacuum BCs were applied at the three outer pool boundaries, and continuity BCs were used at all of the inner surfaces. The governing equations for the boundary conditions are shown in Eqs. (17)–(22). The first equation listed for each BC is written in ‘‘generic’’ terms and the second equation listed for each BC is the equation applied to the COMSOL model. Symmetry – Neumann boundary condition: Fig. 10. Thermal (middle) and fast (right) flux in the flux trap (h = 60.96 cm, r = 6.4 cm).

(a) create the identical geometry that was utilized in the NEWT model, (b) import the macroscopic cross sections previously calculated by NEWT, (c) code the diffusion theory equations into the subdomain settings, (d) define the appropriate boundary conditions in the boundary settings, (e) create an appropriately fine mesh, (f) set up the solver, and finally, (g) perform the actual calculation. The two-dimensional HFIR COMSOL model as it appears in the draw mode of the GUI is shown in Fig. 5. The same dimensions defined in the NEWT model (Table 1) were utilized in the COMSOL model. The PDE coefficient mode allows the user to enter a system of PDEs into the software in the form expressed in Eq. (13). The equations entered in the subdomain settings GUI for the fast flux in the multiplying regions, the fast flux in the non-multiplying regions, and the thermal flux in all the regions are shown in Eqs. (14)–(16), respectively. The subscripts 1 and 2 indicate fast and thermal energy groups, respectively.

n  ðcru þ au  cÞ þ qu ¼ g

ð17Þ

n  ðDr/Þ ¼ 0

ð18Þ

Vacuum – Dirichlet boundary condition:

hu ¼ r

ð19Þ

/ ¼ ð2:1312ÞðDÞjr/jboundary

ð20Þ

Continuity – Neumann boundary condition:

n  ððcru þ au  cÞ1  ðcru þ au  cÞ2 Þ þ qu ¼ g

ð21Þ

n  ððDr/Þ1  ðDr/Þ2 Þ ¼ 0

ð22Þ

COMSOL provides an adaptive mesh refinement algorithm that provides an iterative solution scheme to update the mesh based on the results from the solution. In this manner, machine accuracy may be obtained in the solution, thus yielding only round-off error left in the solution. The automatic mesh refinement algorithm was used to solve the diffusion equations using a sequence of refined meshes. The predefined ‘‘extremely course’’ triangular mesh was used as the initial mesh and was refined five times by the ‘‘longest edge’’ refinement method. When using the ‘‘longest edge’’ refinement method, the longest edge of the triangles that are determined to have the largest errors are bisected in order to increase the number of elements in areas of largest numerical error. The initial mesh was used to solve the system of PDEs and was improved by

D. Chandler et al. / Annals of Nuclear Energy 38 (2011) 2594–2605

2601

Fig. 11. Thermal (middle) and fast (right) flux in the fuel regions (active fuel h = 50.8 cm).

computing local mesh element error indicators, (f(i, j)h(j)b(i))aVol(j), and refining the mesh where the errors are largest. The local error indicators depend on the equation number (i), the mesh element number (j), the mesh element size (h), and the mesh element volume (Vol). The global error is estimated by taking the sum of the local error estimates. The mesh refinement loop continues until the maximum number of refinements is reached or the maximum number of elements is obtained, both of which are user-defined values (COMSOL, 2008). Mesh quality, as described in Knupp (2007), concerns the characteristics of a mesh that permit a particular numerical PDE simulation to be efficiently performed, with fidelity to the underlying physics, and with the accuracy required for the problem. The mesh quality influences the convergence, accuracy, and efficiency of finite-element-based simulations. Since the mesh quality is based on the geometry and the aspect ratio of the element, and not the actual physics being solved, the mesh quality shown is not necessarily an indicator of solution ability. However, the COMSOL documentation (COMSOL, 2008) states that if the quality is at least greater than 0.3, then the mesh quality is sufficient to obtain an accurate solution. A direct (PARDISO) linear system solver was chosen in the solver parameters dialog box because it is more efficient and uses less memory than other solvers. The PARDISO solver is a direct solver for sparse symmetric and unsymmetric linear systems (Ax = b) of equations. In diffusion theory, the matrix A represents the destruction matrix [leakage, absorption, and downscatter (for the fast flux group only)], x is the neutron flux, and b is the production vector [fission and downscatter (for the thermal group only)]. Referring to Eqs (13)–(16), static eigenvalue problems in COMSOL are setup similar to time-dependent problems by linking the time derivative, @//@t, to the eigenvalue, k. Thus, for a critical system (keff = unity), the eigenvalue, k, is equal to zero according

to the set of PDEs being solved and therefore a desired eigenvalue of 0 was defined and the eigenvalue linearization point was set to 0. 5. Results The NEWT code in the SCALE 6 system was used to model HFIR and generate the two-group macroscopic cross sections needed to solve the diffusion equations by calculating zone-averaged neutron fluxes and then using them to collapse the 238-group ENDF/B-VII cross section library to a two-group structure. The effective multiplication factor for the defined configuration was calculated to be 1.0033. The thermal (105 < E < 3 eV) and fast (3 eV < E < 20 MeV) neutron flux distributions as calculated with NEWT are depicted in Figs. 6 and 7, respectively. In these figures, the neutron flux is viewed by the color spectrum scale whereby the red color represents the largest flux and the blue represents the smallest flux. Due to the very fine mesh applied around the CEs, it is difficult to see the flux profile in and around them. A very fine mesh is needed for the discrete representation of the model of HFIR since it is a very compact HEU loaded core and thus there are sharp gradients in the fluxes. The steepest gradients occur near the edges of the fuel elements where the neutrons are moving from fuel regions to moderating [water (FT and flow channels) and beryllium] and absorbing (control elements) regions. The thermal and fast neutron flux distributions as calculated with COMSOL are shown in Figs. 8 and 9, respectively, and are in very good agreement with the NEWT results shown in Figs. 6 and 7. Figs. 8 and 9 are surface plots that show the continuous distribution of the fluxes and contour lines overlay the plots to capture the discrete curves of the solution field. Again, the neutron flux is viewed by the color spectrum scale whereby red represents the largest flux and blue represents the smallest flux. The square root of the thermal flux is plotted along with the thermal flux surface plot in Fig. 8 for the

2602

D. Chandler et al. / Annals of Nuclear Energy 38 (2011) 2594–2605

Fig. 12. Thermal (left) and fast (right) flux in the control elements (h = 102 cm, w = 0.635 cm).

sole purpose of showing more variability in the color spectrum. It is important to note that the square root of the flux has no physical meaning. As mentioned in the previous paragraph, steep flux gradients are unique to the compact HFIR core, which is emphasized in Fig. 8 where the FT is red and all other regions are blue. Region specific thermal and fast neutron flux surface plots for the FT, FE, CEs, and the beryllium reflector regions are illustrated in Figs. 10–13, respectively. The FT, FE, and beryllium reflector plots are bounded by the y = 30.48 cm and y = 30.48 cm planes (active fuel length is only 50.8 cm in length) and the CE surface plot shows the entire length of the elements as modeled. The width-toheight ratio of the CEs was increased for better visibility of the plot and CE drawings were placed next to the plots such that the three regions (white, grey, and black) could be easily identified.

Fast neutrons are born in the fuel regions and leak out into the FT and beryllium reflector regions where they are moderated to lower energies. The fast flux decreases with increasing penetration into the FT and the beryllium reflector regions because they are being thermalized. The thermal flux increases with increasing penetration into the FT and is greatest at center of the core. The thermal flux also increases with increasing penetration into the beryllium reflector and is greatest (in the reflector) at a distance of approximately 4 cm into the reflector (at the horizontal midplane) and then exponentially decreases with distance out of the reflector and into the pool. The fast flux at the horizontal midplane is greatest at the outer edge of the IFE and the inner edge of the OFE and dips slightly in the region between the FEs since fast neutrons are produced in the fuel regions and slow down in non-fuel regions due to scattering and moderating. The fast flux decreases exponentially at the horizontal midplane as a function of distance out of the OFE and into the reflector and out of the IFE and into the FT and this is again due to these fast neutrons being moderated and scattered in the hydrogenous and beryllium regions. The effect of the black CE regions is very apparent in the thermal flux plot where the color darkens around the upper and lower sections of the CEs where the black regions are located. The thermal flux is much larger at the inner edge of the beryllium reflector at the core horizontal midplane than it is at the upper and lower sections of the reflector’s inner edge because the grey regions (moderate neutron absorbers) of the CEs are located in the center 25.4 cm region and the black regions (strong thermal neutron absorbers) of the CEs are located above and below the grey regions for the OCE and ICE, respectively. This effect isn’t as apparent in the fast flux profile because europium has a much larger absorption cross section in the thermal group in comparison to the fast group. The thermal and fast neutron fluxes at the horizontal midplane are compared to benchmarked MCNP (Primm, 2005) and KENO (Chandler and Primm, 2009) axially averaged fluxes in Fig. 14. The fluxes shown in Fig. 14 are normalized such that /fast,max = /thermal,max = 1. It is important to note that the MCNP and KENO fluxes are axially averaged because averaging impacts the flux profile. The two-group MCNP and KENO fluxes were calculated for specific regions since they are transport calculations whereas COMSOL calculated the fluxes at mesh intervals inside regions. The MCNP and KENO fluxes in the FT were averaged over the entire length of a few targets (50.8 cm in length), the fluxes in the FEs were averaged over their active length (50.8 cm), and the fluxes in the reflector were averaged over their length of 60.96 cm. Also, the MCNP input is specific to cycle 400 where no transuranic targets were loaded into the FT and the KENO input was set up for depleting the reflector for numerous cycles and therefore utilized smeared poisons in the CE channel rather than explicitly modeling the CEs in order to produce cycle-averaged fluxes in the beryllium reflector. Although the three models were created for three unique analyses, the flux profiles are in good agreement with each other. The thermal flux profiles for all three models are consistent with each other, but there are small discrepancies in the fast flux profiles. These discrepancies can be attributed to diffusion theory approximations and because the COMSOL-generated fluxes are values for a plane along the axial centerline of the core but are being compared to the axially averaged fluxes generated in MCNP and KENO (see Figs. 6–13). Pertinent mesh data from the COMSOL simulations for the solution based on the initial ‘‘extremely course’’ mesh and each of the five iterative solutions are listed in Table 2 and include the number of DOF, the number of mesh points, the number of elements, the minimum element quality, the global error, the memory usage, and the solution time. During the global adaptive-mesh outer

D. Chandler et al. / Annals of Nuclear Energy 38 (2011) 2594–2605

2603

Fig. 13. Thermal (left) and fast (right) flux in the beryllium reflector (h = 60.96 cm, w = 30.80 cm).

iterations, the number of DOF, mesh points, and elements increased while the global error decreased, which shows that the accuracy of the solution is increasing as the mesh is adapting and the memory usage and solution time are increasing. The mesh quality associated with the initial mesh and the mesh used on the fifth (final) refinement cases are shown in Fig. 15. The mesh quality is viewed by the color spectrum scale whereby the red color represents the highest quality and the blue represents the lowest quality. In the six cases studied here, the minimum element quality ranges from 0.438 to 0.529, which is in the green color range. Compute nodes from an ORNL Research Reactors Division cluster (named Betty) operating on a Linux platform were used for the calculations described in previous paragraphs. The compute nodes used for these calculations each have dual AMD Opteron 2350 (2.0 GHz) quad-core 64-bit processors (total of 8 processors per node) and contain 64 GB of ram on each node. Through the utilization of only a single compute node of the cluster, the solution time for the COMSOL problem with five mesh updates was 17.6 min. The solution time required to run NEWT to generate the cross sections was approximately 9 h and this calculation was also executed on one of the compute nodes in serial mode. The detailed, benchmarked, cycle 400 HFIR MCNP5 model (50 million neutron histories) requires approximately 4 h when running in parallel and distributed over 14 processors (2.5 days in serial). In comparison, it takes approximately 4 h of run time for the KENO V.a/CSAS5 and KENO-VI/CSAS6 models of HFIR to complete when running in serial and simulating 50 million and 1 million neutron histories, respectively. The SCALE 6.0 code system, including the NEWT and KENO codes, only run in serial mode, but future releases of SCALE will allow for parallel processing. 6. Conclusions A two-dimensional, two-group, diffusion neutronics model of the High Flux Isotope Reactor was constructed with COMSOL

Multiphysics. NEWT, a two-dimensional, discrete-ordinates neutron transport code in the SCALE 6 code package, was used to calculate the thermal (105 eV < E < 3 eV) and fast (3 eV < E < 20 MeV) group cross sections. The multi-group cross sections calculated by NEWT were then used in COMSOL to build a diffusion model of HFIR. The PDE coefficient form application mode and eigenvalue analysis were used to implement and solve the diffusion equations. A COMSOL-provided adaptive mesh refinement algorithm was used to increase the number of elements in areas of largest numerical error to increase the accuracy of the solution. The COMSOL simulation of steady-state, beginning-of-cycle HFIR conditions proved that COMSOL is capable of performing neutronic analyses for the compact HFIR core. Fast neutrons are born in the fuel regions due to fission reactions in the highly enriched uranium and are moderated and scattered to thermal energies as they leak from the core into hydrogenous and beryllium regions. The greatest thermal neutron flux is located at the center of the core in the over-moderated flux trap since fast neutrons leak from the fuel elements into the flux trap where they become thermalized. The black regions of the control elements proved to be very absorbing, especially for thermal neutrons. The thermal neutrons were unable to penetrate through the black regions, but were able to penetrate through the grey and white regions and into the beryllium reflector. This model was primarily developed to establish the basis of using COMSOL with neutronics data computed by NEWT to perform HFIR core physics analyses. Since COMSOL proved to be a powerful FEA simulation tool, it will be adopted for more complex and computationally intense neutronics studies in the future. COMSOL is also currently being used at HFIR to update thermal hydraulic and structural methods. The next step in this study is to develop a time-dependent, two-dimensional, three-group, neutron kinetics model of HFIR for the purpose of studying the spatial dependence of transient-induced behavior in the reactor core. The COMSOL models will eventually be merged together to form a comprehensive multiphysics solution for transient-induced behavior.

2604

D. Chandler et al. / Annals of Nuclear Energy 38 (2011) 2594–2605

1.1 FT

COMSOL (thermal at horizontal midplane) MCNP (thermal axially averaged) KENO VI (thermal axially averaged) COMSOL (fast at horizontal midplane) MCNP (fast axially averaged) KENO VI (fast axially averaged)

OFE

IFE

1.0 0.9

Note: flux trap in MCNP model has no transuranic targets and MCNP/KENO fluxes are axially averaged whereas COMSOL fluxes are at the core horizontal midplane.

0.7

Beryllium Reflector

0.6 0.5

Water Reflector

Normalized Flux

0.8

0.4 0.3 0.2 0.1 0.0 0

5

10

15

20

25

30

35

40

45

50

Distance from Core Centerline Fig. 14. Normalized two-group flux profiles.

Fig. 15. COMSOL mesh quality (mesh refinement 0 on left, mesh refinement 5 on right).

55

60

2605

D. Chandler et al. / Annals of Nuclear Energy 38 (2011) 2594–2605 Table 2 Automatic mesh refinement statistics/parameters. Mesh refinement

DOF

Mesh points

Elements

Min. element quality

Global error

Memory (GB)

Cumulative solution time (s)

Solution time (s)

0 1 2 3 4 5

5.016E+04 1.615E+05 3.510E+05 7.010E+05 1.335E+06 2.501E+06

6.293E+03 2.024E+04 4.398E+04 8.778E+04 1.671E+05 3.129E+05

1.250E+04 4.025E+04 8.757E+04 1.749E+05 3.331E+05 6.244E+05

0.529 0.438 0.438 0.438 0.464 0.438

3.024E01 5.147E03 1.591E03 6.924E-04 3.229E04 1.646E04

4.7 5.0 5.4 6.4 8.2 11.0

5.2 31.4 82.6 185.6 392.0 1054.0

5.2 26.2 51.1 103.0 206.5 662.0

Acknowledgements This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC05-00OR22725 with the US Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. References Chandler, D., Primm III, R.T., Maldonado, G.I., 2009. Reactivity Accountability Attributed to Beryllium Reflector Poisons in the High Flux Isotope Reactor, ORNL/TM-2009/188, December 2009. COMSOL, 2008. COMSOL, Inc., COMSOL Multiphysics User’s Guide, Version 3.5a, Burlington, MA. DeHart, M.D., 2009. NEWT: A New Transport Algorithm for Two-Dimensional Discrete Ordinates Analysis in Non-Orthogonal Geometries, ORNL/TM-2005/39, Version 6, vol. II.

Gomes, G., 2008. Comparison between COMSOL and RFSP-IST for a 2-D Benchmark Problem, Atomic Energy of Canada Limited, Excerpt from the Proceedings of the COMSOL Conference 2008 Hannover, Germany, November 2008. Knupp, P.M., 2007. Remarks on Mesh Quality, Sandia National Laboratory, Paper from the 45th American Institute of Aeronautics and Astronautics Aerospace Sciences Meeting and Exhibit 2007, Reno, NV. Lamarsh, J.R., Baratta, A.J., 2001. Introduction to Nuclear Engineering, third ed. Prentice Hall, New Jersey. Memoli, V., et al. 2009. A Preliminary Approach to the Neutronics of the Molten Salt Reactor by means of COMSOL Multiphysics, Politecnico di Milano, Excerpt from the Proceedings of the COMSOL Conference 2009 Milan. Primm III, R.T., Xoubi, N., 2005. Modeling of the High Flux Isotope Reactor Cycle 400, ORNL/TM-2004/251, Oak Ridge National Laboratory, August 2005. Primm III, R.T., et al., 2009. Design Study for a Low Enriched Uranium Core for the High Flux Isotope Reactor, ORNL/TM-2009/87, Oak Ridge National Laboratory, March 2009. Primm III, R.T., et al., 2011. Design Study for a Low-Enriched Uranium Core for the High Flux Isotope Reactor, Annual Report for FY 2010, ORNL/TM-2011/06, Oak Ridge National Laboratory. SCALE: A Modular Code system for Performing Standardized Computer Analyses for Licensing Evaluations, ORNL/TM-2005/39, Version 6, vols. I–III, January 2009. Available from Radiation Safety Information Computational Center at Oak Ridge National Laboratory. Stacey, W.M., 2001. Nuclear Reactor Physics. John Wiley and Sons, Inc., New York.