A 3-D FINITE ELEMENT FEW-GROUP DIFFUSION CODE AND ITS AP- PLICATION TO GENERATION IV REACTOR CONCEPTS

In this paper, a recently developed 3-d few-group finite element-based diffusion code is described. Its geometrical flexibility allows future modeling of complex and irregular geometries of (very) small and medium size reactor concepts – (v)SMRs – being in the spotlight for energy provision in remote residential and industrial regions or for space applications, and also liquid metal cooled Generation IV reactors where thermally induced core deformation results in localized assembly lattice distortions which cannot be treated by traditional 3-d neutron kinetics codes devoted to the regular lattices of LWR and Generation IV systems. The description of the implemented FEM solution method is followed by first applications to a prismatic (or block type) high-temperature reactor MHTGR-350MW within an OECD/NEA benchmark activity and to the sodium cooled fast reactor concept ASTRID within the past EU project ESNII+. Finally, an outlook to planned further code development activities is given.


Galerkin Representation of the Steady-state Few-group Diffusion Equation
Finite elements methods have been used to solve the diffusion and transport equation in rectangular geometry for simulation of PWR and BWR [2,3] and in hexagonal geometries [4,5]. Starting point is the within-group form of the 3-d steady state diffusion equation Using the theorem of Gauß, the second order derivative is transformed into first order derivatives (partial integration) and the volume integral over the first term is transformed into a surface integral, leading to the following Galerkin formulation: Therein,  denotes the surface of the problem domain, and ≔ ⃗ ( )  ⃗ ( ) is the Neumann boundary condition representing the gradient of the neutron flux along the surface normal on . Using a spatial discretization of the problem domain into N finite elements (more details on the elements follow in the next section), the spatial flux distribution ( ) is expanded in terms of the N test functions  ( ) with expansion coefficients fj as follows (the energy group index g is suppressed for better readability): Thus, the test functions  ( ) are also used as basis or shape functions. Inserting this expansion into (2), a linear equation system for the unknown expansion coefficients fj is obtained:

Finite-Element Basis Functions
The finite element representation of the diffusion equation is attributed to the calculation of the elements aij of the coefficient matrix and the elements bi of the vector of the right hand side of Eq. (4). This can be done analytically, provided that the basis functions  ( ) are suitably defined and the neutronic properties (e.g. cross sections, diffusion coefficients, sources) are constant within a finite element. In the following, straight triangular prisms are considered as finite elements. In order to be able to analytically evaluate the integrals occurring in the coefficients, each prism element is mapped to an isoparametric (standard) prism as shown in Figure 1. The Jacobian J for this coordinate transformation is In the simplest case of linear basis functions, there is one basis function for every vertex of the prism, i.e. in total six basis functions per finite element. For the coefficient matrix to be sparse in the linear equation system, the basis functions must, in addition to the property required above, be different from zero only within the respective finite element. One possible choice of linear basis functions, is given by the barycentric coordinates introduced above: For the expansion coefficients fj to represent the neutron flux of energy group g at the nodes j = 1, …, N, the shape functions  ( ) also obey the following condition (k = 1, …, N): The analytic evaluation of the integrals in the coefficients of the equation system is done by the isoparametric coordinate transformation over the standard prism. As an example, the integral of the type (5) is considered in the following. Since there are six basis functions, the result is a 66 matrix of coefficients cij. With the Jacobi-matrix J, the volume element transforms according to d 3 r = det J dr ds dt, so that for c11 the following is obtained: Using block matrix notation, the result for all 36 matrix elements can be written as: with the abbreviation ∶= (   2 1 1  1 2 1  1 1 2 ).
In this way, all volume and surface integrals in the coefficients of the system matrix and the vector of the right-hand side in (4)

Meshing of the Problem Domain and Numerical Solution Scheme
For spatial meshing using straight triangular prisms, FEM-Diff-3d includes a simple meshing module for regular Cartesian and hexagonal lattices. In radial direction, this is effectively a triangulation task. For Cartesian core lattices, each fuel assembly is represented by at least to prisms, whereas in hexagonal core configurations, each fuel assembly is formed by at least six prisms (Figure 2, left). This allows for asymmetric material distribution within one hexagon, which can be useful for, e.g., radially off-centered placed control rods. Spatial mesh refinement in radial and axial direction is possible by subdividing each prism into an arbitrary number of prismatic sub elements as illustrated in Figure 2 (right). The meshing module generates a list of nodes with Dirichlet boundary conditions and a list of elements connectivities which also contains cross section library assignments and Neumann boundary condition information. A standard inner-outer iteration scheme is applied for solving the FEM-discretized diffusion equations. A preconditioned conjugate gradient scheme with incomplete LU preconditioning of the system matrix is used for the numerical solution of the within-group equations. The ILU decomposition is performed by the SPARSKIT2 numerical library [6]. For the solution of the eigenvalue problem, the power iteration method with optional Wielandt shift in connection with a simple global overrelaxation scheme is applied.

The Prismatic High-Temperature Gas-Cooled Reactor MHTGR-350MW
The OECD/NEA Prismatic Coupled Neutronics Thermal Fluids benchmark of the MHTGR-350 MW Core design [7] provides a means for assessing state-of-the-art methods and codes for the assessment of prismatic gas-cooled high temperature reactors (HTR). The radial and axial core layout of the MHTGR-350 which has a nominal thermal power of 350 MW and an average power density of 5.9 W/cm 3 , is shown in Figure 3. The core consists of an annular active zone around a central graphite reflector, followed by a replaceable reflector with hexagonal graphite blocks which is surrounded by a permanent reflector. The active core is formed by fuel columns of 7.9 m height, each being a stack of ten fuel blocks. There are six start-up control rod locations in the central reflector and 24 operating control rod positions in the outer replaceable reflector. The fuel (15.5% enriched UC0.5O1.5) is comprised in TRISO particles bonded in a cylindrical graphite matrix forming a compact. The compacts are placed in hexagonal graphite blocks with a flat-to-flat pitch of 36 cm which also contain coolant holes and burnable poisons. The neutronic specification of Phase I, Exercise 1, is based on a fixed cross section set (no thermalhydraulic parameterization) in 26 energy groups for the EOEC state and is provided to the benchmark participants. It comprises block-wise homogenized cross sections for 12 axial layers with 91 hexagonal blocks per layer in one sixth radial core segment. In total, there are 234 materials where 220 materials represent the active core region. There is one control rod group inserted in the EOEC state with an insertion depth of one fuel block height from the top reflector (axial position 1105.5 cm). The controlled block can be either represented by a full hexagon or by a decomposition of the hexagon into six triangular prisms (denoted by "hexagonal" and "triangular CR model"); appropriate cross section sets are provided. Solutions for Phase I, Exercise 1, have been obtained with FEM-Diff-3d and -for comparison -also with PARCS [8] and TRIVAC [9]; they have been published in [10] and are reproduced in Table 1. Both multiplication factor and average power densities for the hexagonal control rod model (first row in Table 1) are very close to the respective PARCS (using the Triangle-based Polynomial Expansion Nodal (TPEN) method) and TRIVAC (using the Thomas-Raviart-Schneider finite element discretization with parabolic polynomials) solutions. Small differences are observed for the control rod worth, the axial power offset and the maximum power density. The origin for the variations in the maximum power density values may be attributed to the degree of spatial flux resolution within a hexagonal node. Since in PARCS distributions are reported nodewise, the maximum power density is close to the node average and so gives the lowest value.  As seen from Figure 4 (left), the radial power profile compares very well with the solution obtained by PARCS. The higher control rod worth and the smaller axial power offset obtained for the triangular control rod model originate from the more detailed radial resolution of the control rod and the associated different cross section sets. This can be also seen in the axial power profile differences (Figure 4, right).

The Sodium-Cooled Fast Reactor Concept ASTRID
The ASTRID core [11] features an axially highly heterogeneous design (see Figure 5) with four axial layers, a lower axial blanket followed by an inner fissile, inner blanket and upper fissile axial zone for the 177 subassemblies of the inner radial core zone. The 114 subassemblies of the outer radial core zone consist of only one outer fissile axial zone above the lower axial blanket. At operating condition, the active core heights are 110.863 cm and 120.912 cm for the inner and outer core, respectively. Each subassembly consists of a hexagonal array of 217 pins. The pitch of the subassembly lattice equals 17.611 cm. Using the 11 energy groups structure in ref. [13], macroscopic cross sections have been generated for ASTRID with the 2-d lattice and depletion code HELIOS [14] using the 112 energy group data library. For the discretization of the ASTRID core, 6 primatic elements per fuel assembly and 61 axial levels have been used, which resulted in a total of 170810 nodes and 323178 prismatic elements. The calculation time is about 3 minutes on one Intel Xeon CPU E5-2680 v2 at 2.8 GHz with 97 standard (unaccelerated) power iterations. As an application, the voiding scenarios listed in Table 2 have been calculated by FEM-Diff-3d and the sodium void reactivity effects (SVRE) have been compared with a reference solution that was derived from averaging the solutions (Monte Carlo as well as deterministic) obtained by all institutions in the frame of the ESNII+ project. In Table 2, the reference solution for each scenario is given in column 3, the corresponding results obtained with FEM-Diff-3d are provided in column 4, and the last column contains the deviation from the reference solution. Obviously, the FEM-Diff-3d solutions are consistent with the respective references, with differences of up to 354 pcm found only for scenarios 1 and 7 which represent voiding of the upper plenum. The larger differences might be attributed to the limitations of the diffusion approach to describe voided zones (high leakage, neutron streaming effects).

CONCLUSION
The few-group 3-d finite element-based diffusion code described in this paper provides the geometric flexibility required for the safety assessment of not only Generation-IV but also innovative micro reactors like (v)SMR. Whilst the selected HTR and SFR benchmarks have proven its applicability, implementation of an SP3 like transport solver, convergence acceleration methods, application of meshing tools for irregular geometries, extensions for time-dependent simulations and coupling with thermal-hydraulic codes are foreseen and already partially completed, e.g. acceleration of the outer iteration loop by implementing a Wielandt shifted power method, coupling with the GRS thermal hydraulic system code ATHLET [15] and basic implementation of transient simulation capability including treatment of delayed neutrons. The respective work provides the basis for application of FEM-Diff-3d to (v)SMRs with irregular geometry.