Polynomial axial expansion in the Method of Characteristics for neutron transport in 3D extruded geometries

. Recently a solver based on the Method Of Characteristics (MOC) for 3D extruded geometries has been developed in the APOLLO3 R (cid:2) project [1]. In this method the domain is divided in regions and the Step Characteristics (SC) approximation is used in each of them to represent the solution of the transport problem. Since the biggest degree of heterogeneities is found along the radial direction, the idea proposed in this paper is to keep the SC approximation to compute the solution over the radial plane and to implement a polynomial expansion of the ﬂux along the vertical direction. In fact the results of the previous works [1] show that the ﬂux gradient along the z axis are likely to be represented by a classical polynomial basis. In this paper we weigh up the beneﬁts of using a less reﬁned mesh.


Introduction
In the recent years a solver based on the Method of Characteristics (MOC) allowing the treatment of 3D extruded geometries has been developed inside the TDT module of APOLLO3 R [1][2] [3]. The standard Step Characteristics (SC) approximation is used and results show an excellent agreement with Monte-Carlo simulations. However a fine mesh refinement is needed to converge, due to strong flux gradients. In this paper an improvement of this method is proposed: the results of the previous work show that the flux gradients are likely to be represented by a polynomial basis along the vertical direction. Since most of the geometrical and physical heterogeneities are radially located, the SC approach is preserved to represent the solution over the radial plane. As a matter of fact the strong irregularities in the geometrical meshes prevent from an efficient use of a polynomial expansion. On the contrary along the axial direction the computational meshes assume a Cartesian shape, well suited for a polynomial representation of sources and fluxes, as shown in Figure 1. A suitable polynomial development in this direction allows us to approximate the strong flux gradients without the help of a large number of axial meshes. The seminal ideas of this method are already given in [3]. They are here developed in a different form.

Equation derivation 2.1 Polynomial basis definition
The starting point is the mono-group, steady state neutron transport equation is a geometrical domain D and bounde-mail: laurent.graziano@cea.fr e-mail: simone.santandrea@cea.fr e-mail: daniele.sciannandrone@cea.fr Figure 1: Axial distribution of the flux for a thermal (red) and a fast (blue) group inside the fuel pin. The figure also shows the reference Monte Carlo solution with associated error bars [3].
A local polynomial basis can be used to express sources and flux moments: where Δz r is the height of the given region, z r ∈ [−Δz r /2, Δz r /2] is the local coordinate for the r th region, z r is the value of the axial coordinate at the region center and N p is the chosen order for the polynomial expansion. Thanks to this definitionz r ∈ [−1, 1] for every region, independently on the region height. In this way passing from a region to the next one through an horizontal border the polynomial value only change is sign, switching from +1 to −1. This property allows to avoid a certain number of operations during the trajectory sweep, as it will be shown later.
From the local polynomial basis P(z) is possible to retrieve an orthonormalized basis ξ(z) through the Gram-Schmidt orthonormalization algorithm. The coefficients of this algorithm can be used to form a lower triangular matrix that allows to switch from the orthonormal to the polynomial basis: where: The terms < P p , ξ j > are computed using the geometrical matrixP: which elements read: f or p + p even 0 f or p + p odd and V r denotes the volume of the r region.

Source expansion
The source term can now be expanded over the orthonormalized basis: where θ r is the characteristic function for the r th region and q r are the orthonormal coefficients of the source moments: Σ n r is the macroscopic in-group transfer cross section, A n the real spherical harmonics, φ n r the coefficients of the flux moments expansion and S n r is the external source, accounting for fissions and transfer from other groups. The sum over n stands for: Where K is the anisotropy order.
Thanks to the orthonormality property, by projecting q( r, Ω) over the polynomial basis ξ(z) we get: This important property is very useful, as shown later, and justify the choice of a orthonormalized basis.
To switch from the orthonormal source coefficients to the polynomial one the following relation is used, retrieved from the definition (3): From which we get: q r,pol = tM −1 · q r ( Ω) and q r = tM · q r,pol ( Ω) (7) where q r,pol ( Ω) is the coefficient vector referring to the polynomial basis P(z).

Transmission equation
The SC numerical MOC transmission equation obtained from the integral transport equation reads: where Ψ r (0, Ω) is the entering flux at the trajectory beginning, Ψ r (t, Ω) is the flux in a generic point t along the trajectory and: are called respectively the transmission coefficient and the optical length. The source expansion of Eq.(5) is now used. Eq.(3) allows us to switch to the standard polynomial basis and a local system of coordinates along the trajectory is used to express thez component: z in r is the z coordinate of the entering point of the trajectory in the r region, θ is the polar component of the Ω direction. The polynomial transmission equation obtained in this way reads: where: and P k (z in r ) is defined in (2). For numerical reasons it is useful to rearrange the transmission equation as follows: where: This final form of the transmission equation allows to minimize the number of floating point operations to be executed on-the-fly during the sweep, for each chord. The term T , on the other hand, can be pre-computed during the classification phase. The classification process is widely explained in [3]. In a few words, the chord classification is a procedure that exploits the regularity of extruded geometries to recognize chords with the same length. Starting for a given 2D chord, a given polar angle and moving vertically, every chord that lies on that basic 2D chord and intersect two successive vertical surface will have the same length. Fig. 2 gives a visual interpretation of the concept. These chords intersecting two successive vertical surfaces are referred to as the V-chords in [3] . In the same way Hchords are defined as chords intersecting two successive horizontal surfaces, and M-chords, as chords intersecting a horizontal and then a vertical surface (or vice versa). The advantage of the classification process is that when the same chords is found several times, its coefficient can be computed just once. Since the T term is computationally expensive, the possibility to strongly reduce the number of times that it must be computed represent a big advantage.

Computational cost
The transmission sweep is one of the most expensive part of the iterative scheme in computational terms. The transmission equation adapted to the polynomial case derived in this section is clearly more expensive than the classical SC transmission, which asks for one floating point operation per chord (for classified chords). On the other hand the polynomial equation needs 1+N p operations (N p being the polynomial degree). Moreover, it requires to know at each point during the trajectory sweep the value of the axial coordinate and all its powers. If this information is computed on-the-fly it requires 1+N p additional operations. In addition, even the coefficients that ca be precomputed during the classification phase are much more expensive. For these reason for an equivalent tracking (so the same amount and type of chords) the polynomial solution will be more expensive in computational terms.

Escape coefficient computation
The escape coefficients E p−k (τ) can be computed integrating Eq. (9). However, by integrating by parts this equation, a useful recursive relation can be obtained. This allows to compute only the 0 th order coefficient and to retrieve the others: The generic escape coefficient E b can be expressed as a function of the first one through the relation: Unfortunately the recursive relation Eq. (11) is illconditioned for small values of the total cross section. Another possible approach is to start with the highest order coefficient and retrieve the others with the corresponding backward recursive relation. This is a more stable approach for small τ values. In this case the highest order coefficient must be computed. The more efficient way to do this is through the use of tabulated values. Again, small values of the optical length cause big instabilities not only in the use of the recursive relation, but also in the initial computation of the highest order coefficient.
To avoid this problem the strategy adopted is to tabulate a different function. The escape coefficient Eq. (9) is computed as: The function E T b (τ) happens to be well-conditioned even for small τ values. Integrating (12) by parts a forward and a backward relation can be retrieved: Also in this case the backward relation is the more stable for small values of the optical length. A table must be precomputed for the highest order escape coefficient. Instead of integrating the relation (12) a more efficient approach can be the following: the forward recursive relation is used starting from the 0 th order value for large values of the optical length and a Maclaurin series is used to expand the function e −τ for small τ values. For large τ: For small τ:

Balance equation
A balance equation can be obtained by projecting the transport equation (1) over ξ(z): where the integral is performed over the r region.
The following definition is given for the orthonormal coefficients of the angular flux: The equation is rearranged, the basis switch is done again with Eq. (3) and the following useful property of the polynomial basis P(z) is used: The final balance equation obtained is: where the current term is: Here, t Ω j t∩r means that the sum is performed for all the trajectories intersecting the r region parallel to the direction Ω j . TheC matrix reads: ·S ·M and the shift matrixS is:

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
Since the profile of theC matrix is lower triangular with zeros on the main diagonal the components of Ψ( r, Ω) can be computed treating the system as a standard lower triangular system of equations, solving starting from the 0 order polynomial component.

Rebalancing
The tracking strategies used in the Method of Characteristics for three-dimensional geometries are widely explained in [1]. The same strategies apply to the polynomial version, since the same type of geometries are treated. As already said, the use of a polynomial development leads to a lower numbers of axial meshes. This means regions with bigger volumes that could be integrated, with similar relative numerical errors, using a bigger axial integration parameter. Fig.2 represents the radial and axial integration parameters, respectively Δr and Δs. As supposed, to a bigger Δs value corresponded an equivalently precise volume integration. Unfortunately, the same thing can not be said about the precision of the solution. Releasing the Δs value caused big instability issues, that lead in most of the cases to a divergent solution. The problem has been identified in a difference between the balance used to compute the flux moments and the one coming from the transmission equation. A rebalancing technique has been implemented to force the two balances to coincide. The equations for this technique are retrieved starting by integrating the transmission equation along a chord of length L: The only terms depending on t are: So we obtain: Integrating now the equation over the volume of each region we obtain an average angular flux. Multiplying by the real spherical harmonic and integrating over the angular variable we obtain a balance equation for the angular flux moments, only for the constant polynomial value: where the volume integral has been decomposed coherently with the trajectory-based spatial discretization, and computed numerically as: Eq.(6) and Eq.(7) are now used to expand the source term and to switch to the orthonormal basis, obtaining (written as the numerical equivalent form): The last term can be expressed in a tensorial form, defining the operator Γ with dimension N p ← N p * N mom : where ":" refers to the double tensorial reduction. This new balance equation can be compared to the standard balance equation, that reads: where the Spurious matrixS p reads: A rebalancing equation can be obtained forcing the two equations to give the same result. In particular, since only the last term is different in the two equations, a correction to be applied to the 0 th polynomial order transport source can be obtained, in order to let this two terms to coincide: V rS p · q 0,r =Γ : q r + δ q 0,r from which we obtain: δ q 0,r =Γ −1 0 · V rS p · q 0,r −Γ :q r whereΓ 0 is: where theM −1 and his related sum over p disappeared since the matrix is lower triangular.

Results
The polynomial axial MOC method has been validated using the same calculation scheme and the same case of a previous work [3]. This case is a heterogeneous threedimensional assembly of the innovative Sodium-cooled Fast Breeder Reactor ASTRID, developed at the CEA in France. The accuracy of the solution is analyzed comparing the results with a reference Monte Carlo simulation obtained with the TRIPOLI4 R code, also developed at the CEA.
The two-dimensional section of the ASTRID assembly used in the calculation is shown in Fig. 3a, while Fig. 3b shows the actual computational domain. Axially, the geometry is composed by a fertile part superposed to a fissile section.
The multi-group cross sections have been computed with the Subgroup method of the AUTOP module of APOLLO3 R . The self-shielding procedure and the 2D-3D equivalence applied is explained in [2] [3]. A 1968 groups discretization of the energy is used for the multi-group calculation.
The presence of two parts with different materials and composition leads to the strong flux gradient, shown in Fig. 1. For this reason the Step approximation needs a big number of axial meshes, in particular in the interface section, and justifies the idea of a polynomial development. Fig. 4 shows the amount of meshes needed both in the Step and in the polynomial case to properly describe the flux gradients. The use of a polynomial development allowed to decrease the meshes number, from 30 to 2.
Although the number of axial meshes decreased of about one order of magnitude, the number of chords did not decreased proportionally (for the same values of the integration parameters described in section 2.6). In fact, the most abundant chords even with 30 planes in this geometry (and probably in most of them) are the V-chord, connecting two vertical surfaces. The reason lies in the length of the average 2D chord, smaller than the average plane height. Reducing the number of meshes the number of M-chords decreased, while the H-chords basically disappeared, leading to a percentage of V-chords > 95%. Nevertheless the total number of chords slightly decreased. The V-chords became almost the only type. The bigger percentage of V-chords in comparison with the step case represent an advantage, since it is possible to pre-compute the transmission coefficients for classified chords. However, since the amount of V-chords was already high in the step case and since the number of floating operation per chord in the polynomial case is much higher than in the SC approximation, as explained is section 2.3, the hoped advantage in terms of computational time has not been achieved. Calculations have been performed on a Xeon E5-2680@2.8 GHz with two CPUs sharing 256GB of main memory. Each CPU is composed of 10 cores, with 12GB of memory each.  Figure 3: On the left the two dimensional section of the ASTRID assembly with 2π/12 symmetry. The two different colors of the fuel pins represents the two self-shielding areas. On the right the actual computational domain, that corresponds to one twelve of the full assembly.

Conclusions and perspectives
The polynomial axial development has been implemented for the 3D Method of Characteristics. Results show a good agreement with the Monte Carlo simulation and with the results obtained in the previous work, where the 3D SC method has been developed[2] [3].
The tracking size for the Polynomial case, thanks to a lower amount of meshes, is smaller and the percentage of classified chords is bigger. Nevertheless the computational cost is still higher for the polynomial case. Two factors may come in use to increase the computational performances of the Polynomial method, that we will try to develop. First, the DP N acceleration technique that will be developed should take advantage of the lower amount of meshes both for the important memory occupation that this method requires, and for the computational cost that strongly depends on the number of unknowns. Secondly, thanks to the bigger axial size of the meshes, an equivalent precise solution should be obtained using a bigger axial integration parameter (Δs) and a significant lower number of chords. This analysis is already in progress, that is why the rebalancing technique explained in Sec. 2.6 has been implemented. Tests carried out with bigger Δs seem to show equally precise values in term of reactivity. However we cannot consider these results conclusive, since the reactiv-ity and the integrated reactions rates proved to be well calculated, even if the axial profile of the flux was not. This phenomena observed while using the Step method is not yet well understood, however it may cause errors for example in an evolution calculation. For this reason we leave the discussion of the possibility to increment the value of the axial integration parameters for the future, when conclusive results will be available.