HIGH-ORDER FINITE ELEMENTS FOR THE NEUTRON TRANSPORT EQUATION ON HONEYCOMB MESHES

Presently, APOLLO3 R ©/MINARET solves the transport equation using the multigroup Sn method with discontinuous finite elements on triangular meshes with Lagrange polynomial bases. The goal of this work is to solve the spatial problem on hexagonal geometries in the context of honeycomb lattice reactors, without further refining the computational mesh. The idea here is to construct high-order basis functions on the hexagonal element in order to improve the trade-off between computational cost and accuracy, in particular for multiphysics simulations where, often, thermalhydraulic modelling requires only assembly-average cross-sections to be defined (e.g. severe accident of fast breeder reactors) i.e. the assemblies are assumed homogeneous. One approach to achieve this goal is through the use of generalised barycentric functions such as the Wachspress rational functions. This research endeavour deals with the application of Wachspress rational functions to the neutron transport equation for hexagonal geometries up to order 3. With this method, it is possible to decrease the number of spatial unknowns required for the same accuracy, and thus the computational burden for complex geometries, such as honeycomb lattices is reduced.


INTRODUCTION
The neutron transport equation [1] is solved increasingly in the framework of multiphysics coupling schemes. It has been observed that high-fidelity calculations for such cases are still computationally expensive for the neutronics part, especially in the case of parametric or accidental studies. One possible way of dealing with this situation is through the use of modern computing techniques to benefit from new computer architectures. A second means before optimising the code itself is by working on the discretisation scheme (in the present work, the spatial one).
This work lies in the framework whereby the neutron transport equation is solved with discontinuous Galerkin finite elements, but on meshes that are not traditional triangles or quadrilaterals. Presently, APOLLO3 R /MINARET [2,3] solves the transport equation using the multigroup Sn method with discontinuous finite elements on triangular meshes with Lagrange polynomial bases. The goal of this work is to solve the spatial problem on hexagonal geometries in the context of honeycomb lattice reactors, without further refining the computational mesh. Indeed, despite the lack of regularity of the angular flux, it has been shown in past works [4,5](especially on Adaptive Mesh Refinement techniques) that increasing the discretization order is an effective means to significantly improve the accuracy of the solution.
The idea here is then to construct high-order basis functions on the hexagonal element in order to improve the trade-off between computational cost and accuracy. One approach to achieve this goal is through the use of generalised barycentric functions such as the Wachspress rational functions. They have been employed in the past in neutron transport for quadrilaterals and arbitrary polygonals [6,7]. This research endeavour deals with the application of Wachspress rational functions [8] to the neutron transport equation for hexagonal geometries up to order 3. With this method, it is possible to decrease the spatial unknowns required for the same accuracy, and thus the computational burden for complex geometries, such as honeycomb lattices is reduced [9].

Discontinuous Galerkin form of the neutron transport equation
Considering an open spatial domain D of R 2 with boundary ∂D, meshed in a set M h of N κ hexagonal elements κ, the discrete-ordinates equation for the monoenergetic neutron transport equation in a given direction Ω n is expressed as: Ω n · ∇ψ n (r) + Σ t (r)ψ n (r) = q n (r)∀r ∈ D with Σ t being the macroscopic total cross section and q n (r) being the neutron source in direction Ω n * . Given n(r) the unit outward normal to ∂D at r ∈ ∂D, ∂D − is the inflow boundary defined as The boundary conditions for r ∈ ∂D − are of Dirichlet type for ψ n (r). Taking P (κ) as the space of real-valued functions on κ, the set V h is defined as The unknown angular flux ψ n (r) is thus expanded over the functions v h ∈ V h . These functions form a basis of the approximation space and are non-zero over κ and discontinuous across κ|κ interfaces. Defining ψ + n and ψ − n as the internal and external traces respectively, Eq. (2), the weak formulation for the transport equation, is obtained by multiplying Eq. (1) by v h , integrating by parts over element κ, applying the divergence theorem to the streaming term, and using the upwinding * It contains the external sources and the scattering source from other directions. approximation for the numerical flux at the interfaces. More details can be found in [10]. (2) Usually for meshes with triangular or quadrilateral elements, P is the space of polynomials of degree k. However, in the general case of convex polygons, such as hexagons, it cannot be the case if P were the space of C 0 (κ) functions and such that at order k, the approximation space P k ⊂ P .

Wachspress rational functions on the regular hexagon
J. L. Gout has defined Wachspress rational basis functions on the regular hexagon up to order 3 in [11]. For this particular case, the following geometrical elements are given as shown in Fig. 1 and following the same notations as in [11].  Let us consider a regular "flat-top" hexagon with centre (0,0), a pitch (distance between two parallel sides) of p and thus, a sidelength of p where c i is a constant chosen such that w i (a i ) = 1.
Furthermore, the midpoint of a i and a i+1 is defined as a i,i+1 . Thus, the equation of the straight line d i through a i−1,i and a i,i+1 is given as l i (x) = 0. Using these geometrical constructions, and the previously-defined ones, it is now possible to build the second-order Wachspress shape functions: where c i and c i,i+1 are chosen such that w i (a i ) = 1 and w i,i+1 (a i,i+1 ) = 1.
J. L. Gout proved that these rational functions form a basis for spanning space P , and such that at order k, P k (κ) ⊂ P and P k+1 (κ) ⊂ P [8,11]. Besides, given that these functions are anchored at nodes on the sides only, they can be assimilated to serendipity functions [7].

NUMERICAL RESULTS
In this section, we provide numerical results for a 2D one-group transport benchmark which has been set up in [12].

Benchmark description
The benchmark consists of three mixtures with a radial distribution as given in Fig. 3. The cross sections for the mixtures are as given in Table 1. As stated by A. Hébert, this benchmark is not representative of a real-life problem but has been designed to magnify transport and anisotropy effects (as much as 2500 pcm!) within the domain and on the vacuum boundary condition. In this honeycomb lattice, the length of the side of a hexagonal cell is 19 cm.   [12]. Mixture 0 is here only to complete the hexagonal geometry description, and is equivalent to void.

Results analysis
The one-group benchmark is analysed in two different settings: within the APOLLO3 R /MINARET framework and using the Wachspress shape functions finite elements. Both cases use discrete ordinates with product angular hexagonal quadratures shortened as HQ (n,m) in the results tables -n and m are the polar and azimuthal orders respectively, leading to 6nm directions in 2D geometries.
APOLLO3 R /MINARET [3] solves the neutron transport equations with a Discontinuous Galerkin formulation with the usual Lagrange shape functions on unstructured triangular meshes. Thus, each hexagonal cell within the meshed geometry is at least split in 6 equilateral triangles. In Table 2, splits of s (1, 2 and 3) are given and means that the side of a triangular mesh in the meshed geometry should not exceed sidelength s , thereby leading to an unstructured grid of triangles. The solutions are computed for spatial orders of 1 and 2. The reference case is the best-estimate case computed with HQ (4,5) angular quadrature, order 2, and a split of 3 -number of spatial unknowns N unk amounting to 38100. The best-estimate case with the Wachspress shape functions solver (WSF) is HQ (4,5) angular quadrature and spatial basis order of 3, equivalent to N unk = 2286.
For each case, the effective multiplication factor k eff is given and compared to the best-estimate result. Furthermore, the absorption rates R abs are also analysed † . The maximum and average errors are defined as in [12]. The results are summarised in Table 2 for APOLLO3 R /MINARET and in Table 3 for WSF. The k eff for the reference APOLLO3 R /MINARET and WSF solutions differs by less than 1 pcm whereas the max and are both less than 0.1%.
From Table 2, it can be deduced that APOLLO3 R /MINARET satisfactorily predicts the k eff for low-order quadratures and coarser spatial discretisation, although the errors on the absorption rates are higher. It can also be observed that as the computational mesh is refined, the errors on the absorption rates decreases, yet at the expense of more unknowns in the system, and thus a higher computational cost ‡ . For a given computational mesh, increasing the order of the method is more efficient to decrease the error. Thus, these numerical results hint that increasing the spatial basis order would be a better option for solving the transport problem efficiently.
From Table 3, solutions with first-order Wachspress rational functions are satisfactory considering † The absorption cross section is computed as Σ t − Σ s0 .
‡ Throughout this analysis, the computational cost is not measured as time or arithmetic intensity. The WSF solver has not been implemented optimally with this goal.  the N unk involved in the problem. Increasing the order is cost effective in terms of N unk to the quality of the result. For this problem, an order of 2 for a given angular quadrature is satisfactory with respect to the solution obtained. Considering a purely computational aspect of using Wachspress basis functions on unrefined hexagons, a coarser mesh is sufficient, and thus N unk depends only the order of the method, which increases linearly. Indeed, the matricial problem solved locally grows as N 2 dof s and increasing the order for better accuracy remains competitive so long the computation time remains reasonably low compared to refining meshes. Hence, the ideal framework would be to increase the order locally to avoid increasing the global N unk of the problem.
All these results are consistent with TRIVAC and SNATCH results published in [12], and the new capabilities of DRAGON for hexagonal geometries [13].

CONCLUSIONS
In this work, a high-order discontinuous Galerkin formulation using Wachspress rational basis functions has been applied to the S n neutron transport equation on 2D honeycomb lattices. In the context of precise solutions to the transport equation on coarse meshes as it is the case for multiphysics calculations on hexagonal reactor geometries (e.g. fast-neutron reactor application), it is sufficient to use high-order finite element solutions on unrefined honeycomb meshes. This work contributes to that effort by demonstrating that Wachspress rational basis functions of order 2 or 3 on honeycomb lattices are well adapted. Further work on the use of these functions is to extend them in a polygonal framework and set up p-refinement strategies to reduce the number of unknowns.