Numerical analysis of the neutron multigroup SP _ N equations

The multigroup neutron SP N equations, which are an approximation of the neutron transport equation, are used to model nuclear reactor cores. In their steady state, these equations can be written as a source problem or an eigenvalue problem. We study the resolution of those two problems with an H 1 -conforming ﬁnite element method and a Discontinuous Galerkin method, namely the Symmetric Interior Penalty Galerkin method.


Introduction
The neutron transport equation describes the neutron flux density in a reactor core. It depends on 7 variables: 3 for the space, 2 for the motion direction, 1 for the energy (or the speed), and 1 for the time. The energy variable is discretized using the multigroup theory [9]. In this method, the entire range of neutron energies is divided into G intervals, called energy groups. In each energy group, the neutron flux density is lumped and all parameters are averaged. We denote by I G := {1, · · · , G}, the set of energy group indices. Concerning the motion direction, the P N transport equations are obtained by developing the neutron flux on the spherical harmonics from order 0 to order * erell.jamelot@cea.fr † francois.madiot@cea.fr 1 N . This approach is very time-consuming. The simplified P N (SP N ) transport theory [11] was developed to address this issue. The two fundamental hypotheses to obtain the SP N equations are that locally, the angular flux has a planar symmetry; and that the axis system evolves slowly. The neutron flux and the scattering cross sections are then developed on the Legendre polynomials. From a mathematical point of view, SP N equations correspond to tensorized 1D P N transport equations, so that some couplings are missing. Consequently, the SP N equations do not converge to transport equations. Nevertheless, they are commonly used by physicists since their resolution is cheap in terms of computational cost. The order N is odd, and the number of SP N odd (resp. even) moments is N := N +1 2 . We will denote by I e (resp. I o ) the subset of even (resp. odd) integers of the integer set {0, · · · , N }. Finally, the (motion direction and energy) discretization of the neutron flux is such that there are N × G even and odd moments of the neutron flux. We will denote by φ = ( (φ g m ) m∈Ie ) g∈I G ∈ R N ×G the set of functions containing, for all energy group g, the even moments of the neutron flux.
Likewise, we will denote by p the set of functions containing the odd moments of the neutron flux. Note that while modelling the core of a pressurized water reactor, the number of groups if such that 2 ≤ G 30, physicists usually choose N = 1 or 3, more rarely N = 5.

Setting of the model
The reactor core is modelled by a bounded, connected and open subset R of R d , d = 1, 2, 3, having a Lipschitz boundary which is piecewise regular. The coefficients are piecewise regular, so that we split R into N open disjoint parts For this reason, we will use the following space of piecewise regular functions: For a set of functions ψ = (φ g m ) m,g ∈ R N ×G , we make the following abuse of notation: For a set of vector valued functions q = (q g x,m ) m,g d x=1 , we make the following abuse of notation: Let us use these notations: for E ⊂ R d , L(E) = L 2 (E); L := L 2 (R); V := H 1 0 (R); V := H −1 (R) its dual and Q := H(div , R). For W = L(E), L, V or 2 Q we define the product space W := W N ×G endowed with the following scalar product and associated norm: We also set V := (V ) N ×G , L(E) = (L(E)) d and L p (·) = (L p (·)) N ×G .
. We set q x = (q g x,m ) m,g and we use the notation M q = (Mq x ) d x=1 . Given a source term S f ∈ L, the multigroup SP N equations with zero-flux boundary conditions 1 read as coupled diffusion-like equations set in a mixed formulation: When S f depends on φ, the steady state multigroup SP N equations read as the following generalized eigenproblem: The physical solution to Problem (3) corresponds to the eigenfunction associated to the smallest eigenvalue, which in addition is simple [7]. In neutronics, the multiplication factor k ef f = max λ λ characterizes the physical state of the core reactor: if k ef f = 1: the nuclear chain reaction is self-sustaining; if k ef f > 1: the chain reaction is diverging; if k ef f < 1: the chain reaction vanishes.
The matrices H, T e , T o , M f ∈ R N × N G×G are such that ∀(g, g ) ∈ I G × I G (δ ·,· is the Kronecker symbol): • (T e ) g,g := T g e ∈ R N × N denotes the even removal matrix, such that: T g e = diag t 0 σ g r,0 , t 2 σ g r,2 , ... , denotes the odd removal matrix, such that: where ∀m ∈ I e,o , σ g r,m := σ g t − σ g→g s,m , and ∀m > 0, t m > 0. The coefficient σ g t is the macroscopic total cross section of energy group g, and the coefficients σ g→g s,m denote the P N moments of the macroscopic self scattering cross sections from energy group g to itself.
• For g = g: (T e ) g,g := −S g →g e ∈ R N × N denotes the even scattering matrix, such that: S g →g e = diag t 0 σ g →g s,0 , t 2 σ g →g s,2 , ... , R N × N denotes the odd scattering matrix, such that: where σ g →g s,m are the P N moments of the macroscopic scattering cross sections from energy group g to energy group g. • where the coefficient νσ g f is the product of the number of neutrons emitted per fission times the macroscopic fission cross section; and the coefficient χ g is the fission spectrum of energy group g.
The coefficients of the matrices T e,o , M f are supposed to be such that: |σ g→g s,m | ≤ εσ g r,m a.e. in R.
(4) It happens that the coefficient νσ g f vanishes in some regions. Hypothesis 4−(iii) is valid while modelling the core of a pressurized water reactor: the scattering cross-sections are weaker than the removal cross-sections of an order 0 < ε << 1. Thus, the matrices t T e,o are strictly diagonally dominant matrices: they are invertible. Let us set D = t H T −1 o H. Problem 2 can be written as a set of coupled primal diffusion-like equations with single unknown φ ∈ V : The variational formulation of (5) writes: where: Theorem 1 Suppose that D is positive definite. For a given source term S f ∈ L, it exists a unique φ ∈ V that solves Problem 6. In addition, it holds: φ V S f L .
Proof: The bilinear form c and the linear form are continuous and under the hypothesis on D, the bilinear form c is coercive: we can apply Lax-Milgram theorem to conclude. In the same way, Problem 3 can be written as: The variational formulation of (7) writes: where:

Proof: The bilinear form c is a continuous and under the hypothesis on
V is a subset of L with a compact embedding. We can then apply the work of Babuška and Osborn in [2]. Thus, the couple (φ, λ −1 ) is a solution to Problem 8 iff the couple (φ, λ) is an eigenpair of operator T f . Moreover, Problem 8 admits a countable number of eigenvalues.
We propose first to derive conditions on the macroscopic cross sections so that Problems 5 and 7 are well-posed. Then we obtain a priori error estimates for a discretization performed with some H 1 -conforming FEM and a Discontinuous Galerkin method, namely the Symmetric Interior Penalty Galerkin method (SIPG) [8,Chapter 4]. The outline is as follows: in Section 3, we exhibit some conditions so that the matrix T −1 o and T e are positive definite. Then we study the discretization of the source problem (5) in Section 5, and the discretization of the eigenproblem in Section 6. Finally, we perform in Section 7 a numerical study of convergence on a benchmark representative of a nuclear core.
3 Properties of T e and T −1 o Consider the diagonal matrix containing the even (resp. odd) removal macroscopic cross sections: . We split T e,o so that: is the identity matrix, and: We have then: .
Proof: The Taylor expansion of T −1 o writes: We get that ∀X ∈ R N ×G : Under assumptions of Properties 3 and 4 the matrices T e and T −1 o are positive definite. Moreover, one can show that H ∇ φ L ∇ φ L [12]. We infer that the matrix D is positive definite and that there exists a constant C D > 0 such that for all ξ ∈ R N ×G , From now on, we suppose that this property holds. 6

Discretizations
Let T h be a shape-regular mesh of R, with mesh size h. We denote by K its elements and F its facets. To simplify the presentation, we assume that the meshes are such that in every element, the cross-sections are regular. We define by F i h the set of interior faces of T h , F b h the set of boundary facets and We denote by N ∂ the maximum number of mesh faces composing the boundary of mesh elements We will first consider an H 1 -conforming finite element method (FEM). For k ∈ N * , V k h ⊂ V and V k h ⊂ V are the finite dimension spaces defined by: The discrete variational formulation associated to Problem (6) writes: Similarly, the discrete variational formulation associated to Problem (7) writes: Then, we will consider a non-conforming FEM. We define the broken spaces: For (φ, ψ) ∈ V NC × V NC , and T ∈ R N ×G , we set: , and For F ∈ F i h such that F = ∂K 1 ∩ ∂K 2 , we define the average {D ∇ h ψ} and the jump ψ as: where n i is is the unit outward normal to K i at face F and D i = D| Ki , ψ i = ψ| Ki . For F ∈ F b h such that F ∈ K, we set {D ∇ h ψ}| F = D| K ∇ ψ| K and ψ | F = (ψ K )| F n, where ψ K = ψ| K and n is the unit outward normal to K at face F .
For k ∈ N * , V k h,NC ⊂ H 1 (T h ) and V k h,NC are the finite dimension spaces defined by: with where α is a stabilization parameter. The Symmetric Interior Penalty Galerkin method (SIPG) associated to Problem (6) writes: Similarly, the SIPG method associated to Problem (8) writes:

SIPG discretization
Assumption 5.1 (Regularity of exact solution and space V ) Let us denote by W 2,p (T h ) the broken Sobolev space spanned by those functions v such We assume that d ≥ 2 and that there is 2d/(d + 2) < p ≤ 2 such that, for the exact solution φ ∈ V := V ∩ W 2,p (T h ). This holds for our assumptions on the coefficients, which are piecewise constant with respect to the triangulation [15].
Lemma 6 Suppose that (T h ) h is a shape-and contact-regular mesh sequence. Then, we have for all h > 0: where h K is the diameter of element K.
We aim at asserting the discrete coercivity using the following norm: with the jump semi-norm Under assumption (4), there exists β > 0 we have for all so that Lemma 7 (Discrete coercivity) Let α := C 2 tr N ∂ C D β where • C tr results from the discrete trace inequality (17), • N ∂ is defined in Section 4, • C D is defined in (11).
Thus, it only remains to prove boundedness. To this purpose, we need to define V ,h = V + V k h,NC and the following norm and n K is the unit outward normal to K. Following [8, Section 4.2], we obtain the following results.
Proof: As in the continuous case (Theorem 2), since the discretization is conforming, there exists a unique compact operator . According to Thm. 5, the sequence of the operators (T h ) h is pointwise converging towards T . As T h and T are a compact operators, the sequence of the operators (T h ) h is then converging in L (V ) towards T : T h −T L (V ) → 0. The norm convergence guarantees that there is no spectral pollution (see [17]). Morevover, we can apply Theorem 8.3 in [2] to state the error estimate on the eigenvalue. We remark that (M f φ, φ) L is a norm over V λ := {φ ∈ V | ∀ψ ∈ V , c(φ, ψ) = λ f (φ, ψ)} [12, Section 5.2.2 p. 78].

SIPG discretization
We recall that, in this section, we work under the assumption 5.1.
Theorem 12 Let µ be the regularity of the eigenfunction ϕ associated to λ, and ω = min(µ, k). Let λ h be the discrete eigenvalue associated to Problem (16).
Proof: We apply the theory developed in [1]. The proof is decomposed as follows. We first show that there is no spectral pollution. Then, we derive the error estimate.
Let E : V +V k h,NC → V +V k h,NC be the continuous spectral projector relative to λ defined by where Γ is a circle in the complex plane centred at λ which lies in ρ(T | V +V k h,NC ) and encloses no other points of σ(T | V +V k h,NC ). The absence of spectral pollution relies on two properties. First, using interpolation results [8,Assumption 4.31] we have for all where C is a constant independent of h. Second, we have for all φ h ∈ V k h,NC , where we used Theorem 9 in the second line and regularity results [15] in the third line. Applying [1, Theorem 3.7], we obtain that there is no spectral pollution. Moreover, we apply [1, Theorem 3.14] to state the error estimate on the eigenvalue, is the continuous spectral projector of the adjoint operator T * | V +V k h,NC relative toλ. Using again elliptic regularity results [15] and Theorem 9, we obtain Using elliptic regularity results, we get Applying Theorem 9, we infer that This concludes the proof.

Numerical Results
We consider the test case Model 2, case 1 from the benchmark of Takeda and Ikeda [19]. The geometry of the core is three-dimensional and the domain is {(x, y, z) ∈ R 3 , 0 ≤ x ≤ 140 cm; 0 ≤ y ≤ 140 cm; 0 ≤ z ≤ 150 cm}. This test is defined with 4 energy groups, isotropic scattering and vacuum boundary conditions. Figure 1 represents the cross-sectional geometry on the plane z = 75 cm.
Since the scattering is isotropic, the SP 3 formulation can easily be reformulated as a multigroup diffusion problem with 8 energy groups and an isotropic albedo boundary condition [3]. We then made the computations with the PRIAM solver from the code CRONOS2 [14] for the conforming case and with the MINARET solver [13] from the APOLLO3 R code [18] for the SIPG discretization. In Figure 2, we consider the convergence of the fundamental mode where we used the SP 3 formulation with Q 1 finite elements and a regular cartesian mesh of size h. The approximated order of convergence is 2.22.

Core Radial Blanket Reflector
In Figure 3, we consider the convergence of the fundamental mode for different the SP N formulations with discontinuous P 1 finite elements and a prismatic mesh of size h. The approximated orders of convergence are given in Table 1.  Table 1: Approximated order of convergence associated to Figure 3

Conclusion
We did the numerical analysis of the approximation with an H 1 -conforming finite element method of the neutron multigroup SP N equations. We also studied the numerical analysis of the approximation with the Symmetric Interior Penalty Galerkin method of the neutron multigroup SP N equations. We then illustrated numerically the convergence results on a benchmark representative of a nuclear core. Those results can be extended to a mixed finite element method, see [5] for the diffusion case with an H 1 -conforming finite element method.