- Email: [email protected]

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 ﬁnite-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 ﬂuxes. 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 coefﬁcient form PDE application mode and were solved via eigenvalue analysis using a direct (PARDISO) linear system solver. A COMSOL-provided adaptive mesh reﬁnement algorithm was used to increase the number of elements in areas of largest numerical error to increase the accuracy of the solution. The ﬂux distributions calculated by means of COMSOL/SCALE compare well with those calculated with benchmarked three-dimensional MCNP and KENO models, a necessary ﬁrst 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 ﬁrst 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) coefﬁcient application mode and the ﬂux proﬁles in the HFIR core were solved via eigenvalue analysis and a direct (PARDISO) linear system solver. A COMSOL-provided adaptive mesh reﬁnement algorithm was used to solve the diffusion equations using a sequence of reﬁned 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 ﬂuxes 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 ﬂux-trap (FT) in the center of the core and a large beryllium reﬂector on the outside of the core in order to produce a large thermal ﬂux 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 reﬂector, 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 reﬂector 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 ﬂux distributions within a nuclear system and collapse the cross sections into multiple (or single) energy groups as speciﬁed 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 ﬁnite 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 predeﬁned (i.e. built-in) physics/engineering modules or by specifying a system of user-speciﬁc PDEs. Additional built-in modules can be modiﬁed to the user’s needs through equation based modeling capabilities and can be coupled together with other modules, with user deﬁned 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 modiﬁcations 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 conﬁguration 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 conﬁguration of HFIR such that symmetry was only utilized across the core centerline (x = 0, cartesian geometry). The ﬂux 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 reﬂector is modeled as two separate regions, a removable beryllium reﬂector (RB) and a permanent beryllium reﬂector (PB), each composed of beryllium, water, and aluminum. The radial and axial boundaries of the regions deﬁned in the NEWT geometry are deﬁned 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 ﬂuxes were collapsed into two-group form: fast ﬂux (group 1) consisting of neutrons with energies between 3 eV and 20 MeV and thermal ﬂux (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 ﬂux 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 ﬁssion event times the ﬁssion 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 deﬁned in the previous paragraph was developed, the homogenized cross sections were created by combining the ﬂux 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 ﬁnite-difference acceleration was activated to speed convergence since there is a signiﬁcant 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 ﬂuxes in every computational cell and (2) relax convergence such that averaged scalar ﬂuxes within a mixture are converged (DeHart, 2009). The second option is useful for mixtures in which ﬂuxes become very small (large reﬂectors or near a vacuum boundary condition) and since this model utilizes three vacuum boundary conditions and includes a large beryllium reﬂector and water surrounding the core, the second option was utilized. The NEWT geometry, material number assignments, and ﬁne rectangular mesh are depicted in Fig. 3. 3. Derivation of diffusion theory Fig. 7. NEWT half core fast neutron ﬂux distribution.

The derivation of diffusion theory is based on Fick’s law, which governs that there will be a net ﬂow 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 ﬂux 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 coefﬁcient ¼ 3R1tr , r is the gradient operator ¼ dx þ dy þ dz in rectangular coordinates, / is the neutron ﬂux.

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 ﬂuxes and thus the time-dependence term can be neglected. Also, no external sources are present in HFIR and therefore ﬁssion 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) ﬂuxes assuming all ﬁssion 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 ﬂux 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 deﬁnition of steady-state, BOC two-group analysis of HFIR, the following equations are applicable for the fast ﬂux in the multiplying regions (Eq. (7)), the fast ﬂux in the nonmultiplying regions (Eq. (8)), and the thermal ﬂux 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 ﬂux 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 ﬂux 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 ﬂux slope at the boundary to extrapolate the ﬂux outside of the physical boundary a distance, d, where the ﬂux vanishes to zero and assumes that no neutrons are reﬂected back through the boundary. The vacuum BC is used at the three pool

4. COMSOL model development The ﬁrst 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 coefﬁcient 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 coefﬁcient, u is the dependent variable, a is the conservative ﬂux convection coefﬁcient, c is the conservative ﬂux source term, a is the absorption coefﬁcient, b is the convection coefﬁcient, da is the damping/mass coefﬁcient, k is the eigenvalue, k0 is the linearization point for the eigenvalue, ea is the mass coefﬁcient, 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 coefﬁcient modes are coupled together and dependent upon each other since the thermal and fast neutron ﬂuxes are being solved. The dependent variable for the ﬁrst mode is the fast ﬂux and the dependent variable for the second mode is the thermal ﬂux. Thus, Eqs. (14) and (15) are input for the multiplying and nonmultiplying regions, respectively, in the fast ﬂux PDE coefﬁcient mode and Eq. (16) is input for the multiplying and non-multiplying regions in the thermal ﬂux PDE coefﬁcient mode. Since there are 31 unique materials in the model there are 31 distinct equations speciﬁed for the fast ﬂux and 31 distinct equations for the thermal ﬂux (62 distinct sets of cross sections and diffusion coefﬁcients). The boundary conditions were deﬁned in the boundary settings GUI. Boundary conditions were deﬁned for both the fast and thermal ﬂux 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 ﬁrst 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) ﬂux in the ﬂux 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) deﬁne the appropriate boundary conditions in the boundary settings, (e) create an appropriately ﬁne mesh, (f) set up the solver, and ﬁnally, (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 deﬁned in the NEWT model (Table 1) were utilized in the COMSOL model. The PDE coefﬁcient 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 ﬂux in the multiplying regions, the fast ﬂux in the non-multiplying regions, and the thermal ﬂux 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 reﬁnement 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 reﬁnement algorithm was used to solve the diffusion equations using a sequence of reﬁned meshes. The predeﬁned ‘‘extremely course’’ triangular mesh was used as the initial mesh and was reﬁned ﬁve times by the ‘‘longest edge’’ reﬁnement method. When using the ‘‘longest edge’’ reﬁnement 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) ﬂux 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 reﬁning 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 reﬁnement loop continues until the maximum number of reﬁnements is reached or the maximum number of elements is obtained, both of which are user-deﬁned 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 efﬁciently performed, with ﬁdelity to the underlying physics, and with the accuracy required for the problem. The mesh quality inﬂuences the convergence, accuracy, and efﬁciency of ﬁnite-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 sufﬁcient to obtain an accurate solution. A direct (PARDISO) linear system solver was chosen in the solver parameters dialog box because it is more efﬁcient 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 ﬂux group only)], x is the neutron ﬂux, and b is the production vector [ﬁssion 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 deﬁned 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 ﬂuxes 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 deﬁned conﬁguration was calculated to be 1.0033. The thermal (105 < E < 3 eV) and fast (3 eV < E < 20 MeV) neutron ﬂux distributions as calculated with NEWT are depicted in Figs. 6 and 7, respectively. In these ﬁgures, the neutron ﬂux is viewed by the color spectrum scale whereby the red color represents the largest ﬂux and the blue represents the smallest ﬂux. Due to the very ﬁne mesh applied around the CEs, it is difﬁcult to see the ﬂux proﬁle in and around them. A very ﬁne 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 ﬂuxes. The steepest gradients occur near the edges of the fuel elements where the neutrons are moving from fuel regions to moderating [water (FT and ﬂow channels) and beryllium] and absorbing (control elements) regions. The thermal and fast neutron ﬂux 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 ﬂuxes and contour lines overlay the plots to capture the discrete curves of the solution ﬁeld. Again, the neutron ﬂux is viewed by the color spectrum scale whereby red represents the largest ﬂux and blue represents the smallest ﬂux. The square root of the thermal ﬂux is plotted along with the thermal ﬂux 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) ﬂux 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 ﬂux has no physical meaning. As mentioned in the previous paragraph, steep ﬂux 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 speciﬁc thermal and fast neutron ﬂux surface plots for the FT, FE, CEs, and the beryllium reﬂector regions are illustrated in Figs. 10–13, respectively. The FT, FE, and beryllium reﬂector 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 identiﬁed.

Fast neutrons are born in the fuel regions and leak out into the FT and beryllium reﬂector regions where they are moderated to lower energies. The fast ﬂux decreases with increasing penetration into the FT and the beryllium reﬂector regions because they are being thermalized. The thermal ﬂux increases with increasing penetration into the FT and is greatest at center of the core. The thermal ﬂux also increases with increasing penetration into the beryllium reﬂector and is greatest (in the reﬂector) at a distance of approximately 4 cm into the reﬂector (at the horizontal midplane) and then exponentially decreases with distance out of the reﬂector and into the pool. The fast ﬂux 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 ﬂux decreases exponentially at the horizontal midplane as a function of distance out of the OFE and into the reﬂector 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 ﬂux plot where the color darkens around the upper and lower sections of the CEs where the black regions are located. The thermal ﬂux is much larger at the inner edge of the beryllium reﬂector at the core horizontal midplane than it is at the upper and lower sections of the reﬂector’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 ﬂux proﬁle because europium has a much larger absorption cross section in the thermal group in comparison to the fast group. The thermal and fast neutron ﬂuxes at the horizontal midplane are compared to benchmarked MCNP (Primm, 2005) and KENO (Chandler and Primm, 2009) axially averaged ﬂuxes in Fig. 14. The ﬂuxes shown in Fig. 14 are normalized such that /fast,max = /thermal,max = 1. It is important to note that the MCNP and KENO ﬂuxes are axially averaged because averaging impacts the ﬂux proﬁle. The two-group MCNP and KENO ﬂuxes were calculated for speciﬁc regions since they are transport calculations whereas COMSOL calculated the ﬂuxes at mesh intervals inside regions. The MCNP and KENO ﬂuxes in the FT were averaged over the entire length of a few targets (50.8 cm in length), the ﬂuxes in the FEs were averaged over their active length (50.8 cm), and the ﬂuxes in the reﬂector were averaged over their length of 60.96 cm. Also, the MCNP input is speciﬁc to cycle 400 where no transuranic targets were loaded into the FT and the KENO input was set up for depleting the reﬂector for numerous cycles and therefore utilized smeared poisons in the CE channel rather than explicitly modeling the CEs in order to produce cycle-averaged ﬂuxes in the beryllium reﬂector. Although the three models were created for three unique analyses, the ﬂux proﬁles are in good agreement with each other. The thermal ﬂux proﬁles for all three models are consistent with each other, but there are small discrepancies in the fast ﬂux proﬁles. These discrepancies can be attributed to diffusion theory approximations and because the COMSOL-generated ﬂuxes are values for a plane along the axial centerline of the core but are being compared to the axially averaged ﬂuxes 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 ﬁve 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) ﬂux in the beryllium reﬂector (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 ﬁfth (ﬁnal) reﬁnement 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 ﬁve 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 coefﬁcient form application mode and eigenvalue analysis were used to implement and solve the diffusion equations. A COMSOL-provided adaptive mesh reﬁnement 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 ﬁssion 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 ﬂux is located at the center of the core in the over-moderated ﬂux trap since fast neutrons leak from the fuel elements into the ﬂux 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 reﬂector. 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 ﬂux proﬁles.

Fig. 15. COMSOL mesh quality (mesh reﬁnement 0 on left, mesh reﬁnement 5 on right).

55

60

2605

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

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 Reﬂector 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.

Copyright © 2018 KUNDOC.COM. All rights reserved.