Virtual element methods for the spatial discretisation of the multigroup neutron diffusion equation on polygonal meshes with applications to nuclear reactor physics

Highlights

• First application of the virtual element method to the multigroup neutron diffusion equation.

• Application of the virtual element method to arbitrary numbers of coupled PDEs.

• Solution of multigroup neutron diffusion equation on general polygonal grids.

• Solution of the 2D IAEA and 2D C5G7 nuclear reactor physics benchmarks using the virtual element method.

Abstract

The Continuous Galerkin Virtual Element Method (CG-VEM) is a recent innovation in spatial discretisation methods that can solve partial differential equations (PDEs) using polygonal (2D) and polyhedral (3D) meshes. This paper presents the first application of VEM to the field of nuclear reactor physics, specifically to the steady-state, multigroup, neutron diffusion equation (NDE). In this paper the theoretical convergence rates of the CG-VEM are verified using the Method of Manufactured Solutions (MMS) for a reaction-diffusion problem in the presence of both highly distorted and non-convex elements and also in the presence of discontinuous material data. Finally, numerical results for the 2D IAEA and the 2D C5G7 industrial nuclear reactor physics benchmarks are presented using both block-Cartesian and general polygonal meshes.

Introduction

The modern design and analysis of nuclear reactor cores, within contemporary nuclear power plants (NPPs), requires the solution of a complex transport equation which models the migration and interaction of neutrons with the materials within the nuclear reactor core (Lewis and Miller, 1993). This neutron transport equation (NTE) is a partial-integro-differential equation (PIDE) which is a function of a seven-dimensional phase of solution variables which are: position $(x,y,z)$, energy (E), angle $(\theta ,\chi )$ and time (t) (Ornstein and Uhlenbeck, 1937). A mathematical approximation to the NTE is often used in whole core nuclear reactor physics simulations and is called the neutron diffusion equation (NDE). The NDE is derived from the NTE by assuming that the angular variation of neutrons is linearly-anisotropic. Using this approximation to the angular variation of neutrons means that the NDE is only a function of position $(x,y,z)$, Energy (E) and time (t) (Hèbert, 2016). This means that the NDE is much easier to solve than the NTE for whole core nuclear reactor physics simulations. The NDE is a reasonable approximation to the NTE for neutron migration in weakly absorbing media several mean free paths away from any isolated sources, boundaries of the domain or at the interfaces between materials with significantly different material properties (Hèbert, 2016). The focus of this paper will be the solution of the NDE for various nuclear reactor physics verification test cases using the virtual element method (VEM) which is a relatively recent innovation in spatial discretisation techniques (Beirão da Veiga et al., 2013).

A typical industrial nuclear reactor design process involves multiple stages, beginning with the solution of the multi-group NTE on a single nuclear reactor fuel assembly for a large number of energy groups, this is commonly referred to as the nuclear reactor lattice calculation (Enrico Sartori, 2010). The solution from the lattice calculation is then homogenised over space and condensed over energy to reduce the number of energy groups to between 2 and 4 energy groups (Sanchez, 2009, Smith, 1986). The space-homogenised and energy-condensed macroscopic neutron cross-sections are then used to perform whole-core nuclear reactor physics simulations; which typically involves solving the multi-group NDE. Spatial discretisation of the space-homogenised multi-group NDE has typically been accomplished using nodal methods (Enrico Sartori, 2010, Lawrence, 1986). The nodal method is a locally conservative, coarse Cartesian mesh, discretisation method that calculates integral quantities over volumes and surfaces such as reaction rates in a region and net currents in or out of a “control volume” or “node”. Usually the “control volume” or “node” is an entire nuclear fuel assembly with homogenised material properties. The homogenisation of nuclear fuel assembly material properties necessitates a loss of spatial information of the solution field. If a high-fidelity, geometry conforming, spatial discretisation method is required to solve the multi-group NTE the most widespread approach, adopted by industry, is the method of characteristics (MoC) (Brough et al., 1980). The main alternative to the MoC is the finite element (FE) method which saw its first uses in nuclear reactor physics in the early 1970s (Semenza et al., 1972, Kang and Hansen, 1973).

In recent years the FE method has been extended to include adaptive mesh refinement (AMR) algorithms with the adaptivity driven by dual weighted residual (DWR) error estimators (Wang and Ragusa, 2009, Wang et al., 2009). The benefits of the FE method over other spatial discretisation methods such as finite differences (FD) and nodal methods is the ability to discretise differential operators using unstructured meshes. Therefore, the mesh can more closely approximate the complex geometry of the problem when compared to the purely Cartesian meshes employed in FD and nodal methods. However, the FE method is limited to meshes comprised of triangles or quadrilaterals in 2D and tetrahedra, prisms, pyramids or hexahedra in 3D. In addition to this limitation the convergence properties degrade in the presence of highly distorted or stretched elements (Braess, 2007). Nuclear reactor design often generates very complex geometries in which geometric details at small length scales relative to the size of the full nuclear reactor core have to be represented by the mesh. Representing these details with a FE mesh can result in highly distorted elements that can cause degradation of the CG-FEM solution (Vartziotis et al., 2013). Determining when this happens and performing operations to improve the quality of the FE mesh can be challenging; albeit AMR approaches can perform some of these operations (Frey, 2008).

It has been estimated that geometry decomposition, meshing and mesh manipulation can take up to 77% of the total analysis time (Cottrell et al., 2009). Isogeometric analysis (IGA) is a spatial discretisation method that attempts to streamline the computer aided design (CAD) to computer aided engineering (CAE) analysis process by removing the mesh generation step altogether. It uses the same basis functions that are used by CAD software to represent the geometry of the physical domain (Cottrell et al., 2009). Therefore, in theory, an engineer could perform CAE analysis of the physical problem directly on the geometric description of their engineering design without any processing by another software package to perform ancillary mesh generation. Research has been conducted on applying IGA spatial discretisation methods to nuclear reactor physics problems (Owens et al., 2016, Hall et al., 2012, Owens et al., 2017, Welch et al., 2017, Welch et al., 2017, Latimer et al., 2020, Latimer et al., 2020). Other attempts to mitigate or ameliorate the challenge of mesh generation have involved relaxing the strict requirements on mesh topology by constructing polygonal and polyhedral spatial discretisations. These methods include the polygonal finite element method (PFEM) (Wachspress, 1975, Sukumar and Tabarraei, 2004, Sukumar and Malsch, 2006, Tabarraei and Sukumar, 2008), the piecewise linear finite element discretisation (PLFE) (Bailey et al., 2008, Ragusa, 2015) and the mimetic finite difference method (MFD) (Beirão da Veiga et al., 2011, Beirão da Veiga and Manzini, 2008, Lipnikov et al., 2014, Kuznetsov et al., 2004, Hyman et al., 2002). Each of these methods, while solving the problem of restrictive mesh topology, have their own limitations. Within both the PFEM and PLFE methods higher order schemes are impossible while in high-order MFD the method is complex compared to Galerkin schemes. The virtual element method (VEM) represents a viable alternative to all of these methods.

The VEM is a relatively recent innovation in spatial discretisation methods which has been used to solve a wide range of partial differential equations (PDEs) (Beirão da Veiga et al., 2013). Its mathematical formulation is largely the same as the FEM; however it utilises many of the concepts developed in the MFD literature. In particular the VEM has inherited the notion of a discrete field from MFD whereby degrees of freedom (DoFs) are distributed to topological entities such as vertices, edges, faces and cells which are then used to define discrete bilinear forms. The MFD formulation that most closely resembles the VEM presented in this paper may be found in Beirão da Veiga et al. (2011), which was the extension of MFD to arbitrary order for elliptic systems. However, the VEM is sufficiently different to MFD to merit the creation of a new name with a number of recent publications within the literature. The central idea in VEM is to take what are effectively MFD DoFs but to map shape functions to them explicitly (Beirão da Veiga et al., 2013). The form of these shape functions is not known on the interior of the element except that they are composed of both polynomial and non-polynomial parts (Beirão da Veiga et al., 2013). The method avoids performing integrations on the shape functions directly by using projections onto the space of polynomials that may be computed using the degrees of freedom alone (Beirão da Veiga et al., 2013). Because the shape functions need not be known explicitly on the interior of an element the VEM has far greater flexibility in terms of element geometry compared to the finite element method.

Approximation spaces of arbitrary order and arbitrary regularity (Beirão da Veiga and Manzini, 2013) may be constructed on very general element geometries and numerous numerical experiments demonstrate the robustness of the VEM to element degeneracy. Thus far various aspects of VEM theory have been developed including the choice of internal DoFs (Russo, 2016), the development of projection operators (Ahmad et al., 2013), VEM with arbitrary regularity (Beirão da Veiga and Manzini, 2013), higher order VEM on polyhedral meshes (Beirão da Veiga et al., 2017), non-conforming VEM (Cangiani et al., 2016, de Dios et al., 2016, Liu et al., 2017, Berrone et al., 2018) and SUPG stabilised VEM (Berrone et al., 2018, Benedetto et al., 2016). Additionally H(div) and H(curl) conforming VEM were presented in Beirão da Veiga et al. (2015). More recently some effort has been expended in developing orthogonal polynomial bases and alternative stabilisation terms for higher order VEM (Dassi and Mascotto, 2018, Mascotto, 2018, Berrone and Borio, 2017). The main reason for this research is the mitigation of ill-conditioning for higher order VEM schemes. The VEM has been applied to numerous problems in science and engineering including general second order elliptic problems in both primal and mixed form (Beirão da Veiga et al., 2016, Beirão da Veiga et al., 2016), hyperbolic problems (Vacca, 2017), parabolic problems (Vacca and Beirão da Veiga, 2015), polyharmonic problems (Antonietti et al., In Press) and biharmonic problems (Antonietti et al., 2017). Additionally some applied problems have been presented such as linear elasticity and plate bending problems (Gain et al., 2014, Brezzi and Marini, 2013), friction and contact problems (Wang and Wei, 2018, Wriggers et al., 2016). This list is far from exhaustive but provides a sufficient overview of the current range of applicability of the VEM in science and engineering.

The next section will discuss the continuous Galerkin formulation of the VEM (CG-VEM). Section 3 will discuss the CG-VEM discretisation procedure for the multi-group NDE, including a few details on how the final block matrix system of equations is assembled. Section 4 presents two numerical verification tests: the 2D IAEA benchmark and the 2D OECD-NEA C5G7 benchmark. These problems were solved with various types of meshes in order to demonstrate the versatility of the VEM with respect to mesh topology. The centroidal Voronoi tessellation algorithm (Aurenhammer et al., 2013, Okabe et al., 2000), is a computationally efficient and numerically robust mesh generation algorithm which is used to generate the polygonal meshes presented in this paper. These types of mesh generation algorithms can tessellate very complex 2D and 3D geometrical domains. Therefore, this approach can reduce and ameliorate the challenge of mesh generation which was described earlier in this paper, some recent advances may be found in (Landier, 2017, Lee, 2015, Contreras and Hitschfeld-Kahler, 2014, Liu et al., 2017, Kim et al., 2019). A method of manufactured solutions (MMS) is presented in the appendix to verify the order of convergence of the CG-VEM in the presence of both element distortion and concavity. Additionally, an MMS involving a reaction–diffusion equation with piecewise-constant coefficients and piecewise-smooth forcing function was performed. Section 5 will provide a discussion of results and some concluding thoughts.

Note, that we presented an early version of the research presented in this paper at the international conference on nuclear engineering (ICONE) in 2018 which was hosted in London, United Kingdom (UK). The conference proceedings may be found in Ferguson et al. (2018).