# Computational Methods for Multidimensional Neutron Diffusion Problems

Authors: Song Han; Sandra Dulla; Piero Ravetto
Journal: Science and Technology of Nuclear Installations
 (2009)

---

## Abstract

A neutronic module for the solution of two-dimensional steady-state
multigroup diffusion problems in nuclear reactor cores is developed. The module can produce
both direct fluxes as well as adjoints, that is, neutron importances. Different numerical schemes
are employed. A standard finite-difference approach is firstly implemented, mainly to serve
as a reference for less computationally challenging schemes, such as nodal methods and
boundary element methods, which are considered in the second part of the work. The validation of the methods proposed is carried out by comparisons of results for reference structures. In particular a critical problem for a homogeneous reactor for which an analytical solution exists is considered as a benchmark. The computational module is then applied to a fast spectrum system, having physical characteristics similar to the proposed lead-cooled ELSY project. The results  show
the effectiveness of the numerical techniques presented. The flexibility and the possibility to
obtain neutron importances allow the use of the module for parametric studies, design assessments,
and integral parameter evaluations as well as for future sensitivity and perturbation
analyses and as a shape solver for time-dependent procedures.


---

## Body

## 1. Introduction and Objective of the Work

The European ELSY Project focuses on the assessment of the fast spectrum lead-cooled reactor technology [1]. Various sophisticated numerical tools are today available for the detailed neutronic design of such a system. However, a simple multidimensional computational tool may in any case be rather useful for its flexibility and the possibility to retrieve information on the physical features of the system and to carry out extensive parametric studies. Furthermore, it can be appropriate for an easy and efficient coupling with a thermal-hydraulic module to allow investigations on nonlinear effects. For the same reasons, while it may be very complicate to couple high-performance computational codes to a kinetic module within a quasistatic framework, the tool developed by the authors in this work can be easily used for the evaluation of neutron shapes for system dynamic computations.The neutronic module is based on a coarse-mesh approach for the calculation of the neutron distribution of the full core. While several nodal methods are available for neutronic calculations, it is interesting to compare and benchmark their performance and to assess their convergence trends. The present work aims at the evaluation of the performance of different computational schemes, starting from a standard finite volume approach in two-dimensional multigroup diffusion on a fine-mesh structure and then proceeding to a nodal formulation of the same problem on a coarse-mesh discretization [2], with the inclusion of a computational option based on the boundary element method (BEM) [3, 4] on the same coarse mesh. The interest of such methods is due to the possibility to develop a computational algorithm based on a response matrix formulation, thus allowing an efficient treatment of a full core geometry including the possibility to easily improve the physical model described (i.e., passing from diffusion to transport), retaining the same structure of the numerical code. The comparison of the performances of the two algorithms (nodal and BEM) for the evaluation of the flux distribution in a full core, both in terms of accuracy of prediction and efficiency of calculation, constitutes an important step in the assessment of coarse-mesh methods for the application to fast reactors. The introduction of the BEM approach in the code developed can also be of interest for further developments that may lead to the possibility of solving the problem with a higher-order neutronic model with respect to diffusion, such as a second-order form of the transport equation [5, 6].As pointed out, the availability of a steady-state flux solver is rather important also for the development of a time-dependent module which is needed for safety evaluations and assessments, within a quasi-static treatment [7]. For that purpose an adjoint solution is also needed. The adjoint can also be used for sensitivity analyses which are important for the preliminary design of a new nuclear system.In a first stage a fine-mesh finite-difference module in 2D cartesian geometry is developed, as a reference validation tool for coarse mesh and boundary element algorithms. On its turn, such module is validated against analytical results for homogeneous configurations.In a second stage the implementation of a nodal method is carried out. A polynomial representation of the neutron flux within a node is used up to second and third order. Eigenvalue and flux results are reported for different number of meshes. The results show the good performance of the coarse-mesh technique associated to a reduced computational effort.The introduction of a boundary element formulation is then performed and results are presented and compared with the previous techniques. The flexibility of the scheme could be of great importance for future applications to different configurations of lead-cooled fast reactors.


## 2. The Physical Model

For the neutronic design of a nuclear reactor core the solution of the system of multigroup diffusion equations is usually sufficient. The physical characteristics of a lead-cooled reactor call for a multidimensional treatment. In this paper the diffusion model is adopted and the equations are solved for a two-dimensional configuration. The neutron steady-state multigroup diffusion equations without and with an external independent source, respectively, are the following:(1)-∇·Dg∇ϕg(r)+Σrgϕg(r)=∑g′=1,g′≠gGΣsg′→gϕg′(r)+1kχg∑g′=1GνΣfg′ϕg′(r),g=1,2,…,G,(2)-∇Dg∇ϕg(r)+Σrgϕg(r)=∑g′=1,g′≠gGΣsg′→gϕg′(r)+χg∑g′=1GνΣfg′ϕg′(r)+sg(r),g=1,2,…,G,
where Dg is the diffusion coefficient of group g; Σrg≡Σag+Σsg-Σsg→g is the removal cross-section of group g; Σsg→g′ is the scattering cross-section from group g to group g′ ; νΣfg′ is  the mean number of neutrons per fission times the fission cross-section of group g′;χg is the neutron spectrum of group g; ϕg is  the neutron scalar flux for group g; sg is the neutron source for group g; k is the effective multiplication constant.For the homogeneous equation (1) the standard effective multiplication constant is introduced as an eigenvalue to guarantee solvability.


## 3. Finite-Difference and Nodal Approach

The multigroup diffusion equations without independent source in two-dimensional Cartesian coordinate systems are considered. The solution of the source-driven problem is derived as one iteration step in the power method used to obtain the eigenvalue. The equations can be written explicitly in the following form:(3)-∂∂xDg∂∂xϕg(x,y)-∂∂yDg∂∂yϕg(x,y)+Σrgϕg(x,y)=∑g′=1,g′≠gGΣsg′→gϕg′(x,y)+1kχg∑g′=1GνΣfg′ϕg′(x,y),g=1,2,…,G.A numerical solution can be readily obtained by the standard finite-difference method, once the domain withx-dimension a and y-dimension b is subdivided into I×J meshes that may be chosen as Δx=a/I and Δy=b/J. The finite-difference solution is used here as a reference for the assessment of the coarse-mesh schemes. The multiplication eigenvalue is determined by the classical power method [8].Nodal methods have proved to be very effective for the accurate solution of reactor core physics problems in nuclear engineering applications. The general principle of nodal diffusion theory methods require to subdivide the reactor core into a relatively small number of regions, or nodes, which present dimensions larger than the physical characteristic length associated to the diffusion process. The detailed flux distribution within each node is approximated by the superposition of a suitable set of spatial functions. The global flux distribution is then determined by the application of a weighted residual technique associated to proper conservation principles and coupling conditions for the nodes. In the following the derivation of the nodal scheme is presented for a general 3D system. The node is defined by its sidesΔx, Δy, and Δz. For each spatial coordinate a one-dimensional equation is obtained by integration over the transverse coordinates. For instance, by integration over y and z, the following equation is obtained for node n, centered at point (xn,yn,zn), and group g:(4)ddxJ¯gxn(x)+1ΔyLnyg(x)+1ΔzLnzg(x)+Σtngϕ¯gxn(x)=∑g′=1GΣng′→gϕ¯g′xn(x)+1kχg∑g′=1GνΣfng′ϕ¯g′xn(x),g=1,…,G,
where
(5)ϕ¯gxn(x)≡1ΔyΔz∫yn-Δy/2yn+Δy/2dy∫zn-Δz/2zn+Δz/2dzϕgn(x,y,z),J¯gxn(x)≡1ΔyΔz∫yn-Δy/2yn+Δy/2dy∫zn-Δz/2zn+Δz/2dzJgn(x,y,z),Lnyg(x)=1Δz∫zn-Δz/2zn+Δz/2dzny·[Jg(x,y2,z)-Jg(x,-y2,z)],Lnzg(x)=1Δy∫yn-Δy/2yn+Δy/2dynz·[Jg(x,y,z2)-Jg(x,y,-z2)].
Following a well-assessed nodal procedure [9, 10] the unknown transverse-averaged flux ϕ¯gxn is expressed as a superposition of suitable polynomial functions as(6)ϕ¯gxn(x)≅ϕ¯gnf0(x)+∑i=13agxinfi(x),xn-Δx2≤x≤xn+Δx2.
If the polynomials herewith appearing are chosen according to the following formulae:(7)f0(x)=1,f1(x)=xΔx≡ξ,f2(x)=3ξ2-14,f3(x)=2ξ(ξ-12)(ξ+12),
the coefficient of the zeroth-order terms turns out to play the role of node-averaged flux. Furthermore, since f3(xn±Δx/2)=0, the values of the fluxes at the node boundary can be related to the expansion coefficients by(8)agx1n=ϕgx+n-ϕgx-n,agx2n=ϕgx+n+ϕgx-n-2ϕ¯gn.
Outgoing surface currents can also be related to the expansion coefficients and to the incoming currents:(9)(Jgx+n)out=-DngΔx(agx1n+3agx2n+12agx3n)+(Jgx+n)in,(Jgx-n)out=DngΔx(agx1n-3agx2n+12agx3n)+(Jgx-n)in.The weighted residual method is now applied for the determination of a relationship involving the unknown coefficientsagx3n. This can be achieved by spatial integration after the application of the weighting function ωpg=ξ2p+1 for any group, with p being any positive integer. The following expression is readily obtained:(10)-12DngΔx2agx3n+Σtng[agx1n+14(2p+32p+5-1)agx3n]=(∑g′=1GΣng′→g+χgk∑g′=1GνΣfng′)×[ag′x1n+14(2p+32p+5-1)ag′x3n].
It is worth taking the limit for p→∞, yielding(11)-12DngΔx2agx3n+Σtngagx1n=(∑g′=1GΣng′→g+χgk∑g′=1GνΣfng′)ag′x1n.
The physical meaning of this process has been discussed in the previous reference works [9, 10], where the present nodal approach is described in full detail.Combining (11) with (8)-(9), a relationship connecting incoming and outgoing partial currents with the average flux is obtained(12)(Jgx+n)out(1+Δx12Σtng+8ΔxDng)+(Jgx-n)out(-Δx12Σtng+4ΔxDng)-6ΔxDngϕ¯gn=(Jgx+n)in(1-Δx12Σtng-8ΔxDng)+(Jgx-n)in(Δx12Σtng-4ΔxDng)+Δx24(∑g′=1GΣng′→g+χgk∑g′=1GνΣfng′)ag′x1n,(13)(Jgx-n)out(1+Δx12Σtng+8ΔxDng)+(Jgx+n)out(-Δx12Σtng+4ΔxDng)-6ΔxDngϕ¯gn=(Jgx-n)in(1-Δx12Σtng-8ΔxDng)+(Jgx+n)in(Δx12Σtng-4ΔxDng)-Δx24(∑g′=1GΣng′→g+χgk∑g′=1GνΣfng′)ag′x1n.
Alternatively the above formulae can be given the following form:(14)s11(Jgx-n)out+s12(Jgx+n)out+s13ϕ¯gn=s14(Jgx-n)in+s15(Jgx+n)in+s16,s12(Jgx-n)out+s11(Jgx+n)out+s13ϕ¯gn=s15(Jgx-n)in+s14(Jgx+n)in-s16,
with appropriate definition of the coefficients. The formula above clearly shows the connection between boundary partial currents and node average flux. The algebraic problem is then solved by standard techniques and outer iterations are carried out for the determination of the multiplication eigenvalue, as discussed before.


## 4. Boundary Elements Approach

The Boundary Elements Method (BEM) [11] constitutes an innovative approach for the effective solution of problems in applied sciences. In the previous works the method has been exploited for applications to the solution of multigroup diffusion problems [3, 4]. Also applications to transport problems in space second-order forms have been recently proposed. Results have proven the effectiveness of the technique and its excellent properties [5, 6].In the following, the basics of a formulation of the boundary element method for the diffusion problem are outlined. The neutron diffusion equation in a homogeneous regionV with a closed smooth boundary surface A is written as(15)D∇r2ϕ(r)-Σaϕ(r)+q(r)=0.
The fission production and the group-to-group transfer in the multigroup model may be included in the source term. The fundamental solution, that is, the Green function, is found by solving the problem:(16)D∇r2ϕ̃(r,r′)-Σaϕ̃(r,r′)+δ(r-r′)=0.
Using the properties of the Green function and integrating over the volume equation (15), the following integral formulation is obtained:(17)c(r)ϕ(r)+∫A[J̃n′(r,rA′)ϕ(rA′)+ϕ̃(r,rA′)Jn′(rA′)]dA′=∫Vϕ̃(r,r′)q(r′)dV′,
where(18)Jn′(rA′)=-D∂ϕ(rA′)∂nA′,J̃n′(r,rA′)=D∂ϕ̃(r,rA′)∂nA′,
and c(r) is equal to 0, 1/2, or 1, depending on whether r is outside V, on the boundary of V or inside V, respectively. Equation (17) allows to obtain ϕ(r) and Jn(r) at any point inside V, provided that the flux and the normal current at the points of the boundary surface are known. If we take r=rA, (17) becomes an integral equation(19)c(rA)ϕ(rA)+∫A[J̃n′(rA,rA′)ϕ(rA′)+ϕ̃(rA,rA′)Jn′(rA′)]dA′=∫Vϕ̃(rA,r′)q(r′)dV′,
in terms of the boundary values of ϕ and Jn. Introducing the partial currents(20)Jn±(rA)=14ϕ(rA)±12Jn(rA),J̃n±(rA)=14ϕ̃(rA)∓12J̃n(rA)
into (19), one has(21)c(rA)2Jn+(rA)+∫AJ̃n′+(rA,rA′)Jn′+(rA′)dA′=-c(rA)2Jn-(rA)+∫AJ̃n′-(rA,rA′)Jn′-(rA′)dA′+14∫Vϕ̃(rA,r′)q(r′)dV′.


### 4.1. Multigroup Boundary Integral Equations

The multigroup equations can be cast into matrix form as(22)(∇r2I+Q)ϕ(r)+q(r)=0,
where ϕ and q are G-dimensional vectors. In the case of a subcritical system driven by an external source the following definition holds:(23)q(r)=s(r)+Fϕ(r),
where the fission operator takes the form(24)F=((χ1D1)ν1Σf1(χ1D1)ν2Σf2…(χ1D1)νGΣfG(χ2D2)ν1Σf1(χ2D2)ν2Σf2…(χ2D2)νGΣfG⋮⋮⋱⋮(χGDG)ν1Σf1(χGDG)ν2Σf2…(χGDG)νGΣfG).
On the other hand, for the homogeneous problem,(25)q(r)=1kFϕ(r).The eigenequation(26)Qσ=λσis assumed to haveG real distinct eigenvalues λh(h=1,2,…,G). S is defined as the matrix which has the normalized eigenvectors σh as its columns. Thus, the following equality is verified:(27)S-1QS=Λ,whereΛ=diag[λ1,…,λG].Now setting:(28)ϕ(r)=Sψ(r),q(r)=Sη(r),
after left multiplication by S-1, (22) reads
(29)∇r2ψ(r)+Λψ(r)+η(r)=0.
By redefining λh=-γh2, the following set of equations is obtained for each component of vector ψ:
(30)∇r2ψh(r)-γh2ψh(r)+ηh(r)=0,h=1,2,…,G.
The corresponding Green fundamental solutions of the following equations:
(31)∇r,r′2ψ̃h(r,r′)-γh2ψ̃h(r,r′)+δ(r-r′)=0,h=1,2,…,G,
are easily found out as
(32)ψ̃h(r,r′)=12πK0(γh|r-r′|),h=1,2,…,G,
together with their derivatives(33)∂ψ̃h(rΓ′)∂n′=γh2πK1(γh|r-r′|)(r-r′)·n′|r-r′|,h=1,2,…,G,
where K0 and K1 are the modified Bessel functions of orders 0 and 1, respectively.Boundary integral equations can be generated with the same procedure leading to (17), obtaining(34)c(r)ψh(r)+∫Γ[ψ̃h(r,rΓ′)∂nΓ′ψh(rΓ′)-ψ̃h(r,rΓ′)∂ψh(rΓ′)∂nΓ′]dΓ′=∫Ωψ̃h(r,r′)ηh(r′)dΩ′,h=1,2,…,G.
By left multiplying the set of (34) written in matrix form by S, one gets back to relations involving ϕ, namely(35)c(r)ϕg(r)+∑g′=1G∫Γ[Ψ̃gg′(r,rΓ′)∂nΓ′ϕg′(rΓ′)-Ψ̃gg′(r,rΓ′)∂ϕg′(rΓ′)∂nΓ′]dΓ′=∑g′=1G∫ΩΨ̃gg′(r,r′)qg′(r′)dΩ′,g=1,2,…,G,
where(36)Ψ̃gg′(r,rΓ′)=∑h=1Gσghψ̃h(r,rΓ′)σhg′*.Letr be taken on the boundary, r=rΓ. Then, finally, the multigroup boundary integrated equations can be explicitly written down:(37)c(rΓ)2Jg+(rΓ)+∑g′=1G∫ΓJ̃gg′+(rΓ,rΓ′)Jg′+(rΓ′)dΓ′=-c(rΓ)2Jg-(rΓ)+∑g′=1G∫ΓJ̃gg′-(rΓ,rΓ′)Jg′-(rΓ′)dΓ′+14∑g′=1G∫ΩΨ̃gg′(rΓ,r′)qg′(r′)dΩ′,g=1,2,…,G,
where(38)Jg±(rΓ)=14ϕ(rΓ)∓12Dg∂ϕ(rΓ)∂nΓ,J̃gg′±(rΓ,rΓ′)=14Dg′Ψ̃(rΓ,rΓ′)±12∂Ψ̃(rΓ,rΓ′)∂nΓ′.


### 4.2. Reduction of the Domain Integrals of the Source Term

The domain integrals involving the independent and flux dependent fission sources in (37) can be reduced into a boundary integrated form as(39)∫ΩΨ̃gg′(rΓ,r′)qg′(r′)dΩ′=∑h=1G-σghσhg′Bh2×{c(rΓ)qg′(rΓ)+∫Ωψ̃h(rΓ,r′)∇2qg′(r′)dΩ′+∫Γ[∂ψ̃h(rΓ,rΓ′)∂nΓ′qg′(rΓ′)-ψ̃h(rΓ,rΓ′)∂qg′(rΓ′)∂nΓ′]dΓ′}=∑h=1G-σghσhg′Bh2{c(rΓ)qg′(rΓ)ψ̃h∂nΓ′+∫Γ[∂ψ̃h(rΓ,rΓ′)∂nΓ′qg′(rΓ′)-ψ̃h(rΓ,rΓ′)∂qg′(rΓ′)∂nΓ′]dΓ′}h=1,2,…,G.
The source distributions qg(r) are now assumed to be either spatially constant or first-order polynomials. In the latter case for 2D geometry the following expression is assumed:(40)q(x,y)=l1x+l2y+q̅,
where, obviously, q̅ is the average source on the domain x∈[0,a] and y∈[0,b]. The following relationship holds:(41)qg′(x,y)=l1x+l2y+q̅g′-(al12+bl22).
As q̅ in the presence of fission emissions is evaluated by the flux of the previous generation on the boundary, the coefficients of expansion (40) can be evaluated by using a least-square approach. The problem can be nicely written in matrix form(42)Ml=N,
having set the following definitions:(43)M=(x1-Δx2y1-Δy2x2-Δx2y2-Δy2⋮⋮xN-Δx2yN-Δy2),l=(l1l2),N=(q1-q̅q2-q̅⋮qN-q̅),
where (xi,yi) are points on the boundary at which the source takes the value qi. By left multiplication of (42) by MT, the normal equations are obtained for the problem
(44)MTMl=MTN,
or, explicitly(45)(∑n=1N(xn-a2)2∑n=1N(xn-a2)(yn-b2)∑n=1N(xn-a2)(yn-b2)∑n=1N(yn-b2)2)(l1l2)=(∑n=1N(xn-a2)(qg,n-q̅g)∑n=1N(yn-b2)(qg,n-q̅g)).
The above system of equations yields the value of the two-dimensional vector l.


### 4.3. Discretization of the Boundary Integral Equations

In this section the numerical solution of (37) is outlined. The boundary Γ is partitioned into I elements of length di. After denoting the two extremes and the middle point by rΓ,i-1/2, rΓ,i+1/2, and rΓ,i, respectively, a local dimensionless coordinate τ is introduced to denote each point r along Γi as(46)r=rΓ,i+τ(rΓ,i+1/2-rΓ,i-1/2).
Thus, it is assumed that(47)Jg±(rΓ)=Jg,i±.
Then, assuming a linear behaviour on each boundary mesh of the integrand functions, the boundary integrals are approximated by finite sums as(48)∫ΓJ̃gg′±(rΓ,i,rΓ′)Jg′±(rΓ′)dΓ′≅∑j=1I∫-11J̃gg′±(rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2))Jg′,j±dj2dτ.
By indicating(49)qg(rΓ)=qg,i,∂qg(rΓ)∂nΓ=qg,i′,
one has(50)∫Γ[∂ψ̃h(rΓ,rΓ′)∂nΓ′qg′(rΓ′)-ψ̃h(rΓ,rΓ′)∂qg′(rΓ′)∂nΓ′]dΓ′≅∑j=1I∫-11{∂ψ̃h∂nΓ′[rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2)]qg′,j-ψ̃h[rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2)]qg′,j′∂ψ̃h∂nΓ′}dj2dτ.
The integrals herewith appearing are obtained by application of a numerical integration formula, namely(51)∫-11f(τ)dτ≃∑k=1Kωkf(τk),
where ωk and τk are the quadrature weights and abscissae, respectively. In the present paper a standard Gauss-Legendre formula is used.At last the problem is cast into the following algebraic response-matrix form:(52)∑g′=1G∑j=1IMgg′,ijJg′j+=∑g′=1G∑j=1INgg′,ijJg′j-+∑g′=1G∑j=1IKgg′,ijqg′j+∑g′=1G∑j=1ILgg′,ijqg′j′,i=1,2,…,I,g=1,2,…,G,
where(53)Mgg′,ij=c(r)2δi,jδg,g′+dj2∫-11J̃gg′-[rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2)]dτ,Ngg′,ij=-c(r)2δi,jδg,g′+dj2∫-11J̃gg′+[rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2)]dτ,Kgg′,ij=-∑h=1Gσgh14Bh2σhg′*×{c(r)δi,j+dj2∫-11∂ψ̃h∂nΓ′[rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2)]dτ},Lgg′,ij=∑h=1Gσgh14Bh2σhg′*dj2∫-11ψ̃h[rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2)]dτ.
The problem for the whole system where nodes are coupled is then solved by numerical iteration schemes.


## 4.1. Multigroup Boundary Integral Equations

The multigroup equations can be cast into matrix form as(22)(∇r2I+Q)ϕ(r)+q(r)=0,
where ϕ and q are G-dimensional vectors. In the case of a subcritical system driven by an external source the following definition holds:(23)q(r)=s(r)+Fϕ(r),
where the fission operator takes the form(24)F=((χ1D1)ν1Σf1(χ1D1)ν2Σf2…(χ1D1)νGΣfG(χ2D2)ν1Σf1(χ2D2)ν2Σf2…(χ2D2)νGΣfG⋮⋮⋱⋮(χGDG)ν1Σf1(χGDG)ν2Σf2…(χGDG)νGΣfG).
On the other hand, for the homogeneous problem,(25)q(r)=1kFϕ(r).The eigenequation(26)Qσ=λσis assumed to haveG real distinct eigenvalues λh(h=1,2,…,G). S is defined as the matrix which has the normalized eigenvectors σh as its columns. Thus, the following equality is verified:(27)S-1QS=Λ,whereΛ=diag[λ1,…,λG].Now setting:(28)ϕ(r)=Sψ(r),q(r)=Sη(r),
after left multiplication by S-1, (22) reads
(29)∇r2ψ(r)+Λψ(r)+η(r)=0.
By redefining λh=-γh2, the following set of equations is obtained for each component of vector ψ:
(30)∇r2ψh(r)-γh2ψh(r)+ηh(r)=0,h=1,2,…,G.
The corresponding Green fundamental solutions of the following equations:
(31)∇r,r′2ψ̃h(r,r′)-γh2ψ̃h(r,r′)+δ(r-r′)=0,h=1,2,…,G,
are easily found out as
(32)ψ̃h(r,r′)=12πK0(γh|r-r′|),h=1,2,…,G,
together with their derivatives(33)∂ψ̃h(rΓ′)∂n′=γh2πK1(γh|r-r′|)(r-r′)·n′|r-r′|,h=1,2,…,G,
where K0 and K1 are the modified Bessel functions of orders 0 and 1, respectively.Boundary integral equations can be generated with the same procedure leading to (17), obtaining(34)c(r)ψh(r)+∫Γ[ψ̃h(r,rΓ′)∂nΓ′ψh(rΓ′)-ψ̃h(r,rΓ′)∂ψh(rΓ′)∂nΓ′]dΓ′=∫Ωψ̃h(r,r′)ηh(r′)dΩ′,h=1,2,…,G.
By left multiplying the set of (34) written in matrix form by S, one gets back to relations involving ϕ, namely(35)c(r)ϕg(r)+∑g′=1G∫Γ[Ψ̃gg′(r,rΓ′)∂nΓ′ϕg′(rΓ′)-Ψ̃gg′(r,rΓ′)∂ϕg′(rΓ′)∂nΓ′]dΓ′=∑g′=1G∫ΩΨ̃gg′(r,r′)qg′(r′)dΩ′,g=1,2,…,G,
where(36)Ψ̃gg′(r,rΓ′)=∑h=1Gσghψ̃h(r,rΓ′)σhg′*.Letr be taken on the boundary, r=rΓ. Then, finally, the multigroup boundary integrated equations can be explicitly written down:(37)c(rΓ)2Jg+(rΓ)+∑g′=1G∫ΓJ̃gg′+(rΓ,rΓ′)Jg′+(rΓ′)dΓ′=-c(rΓ)2Jg-(rΓ)+∑g′=1G∫ΓJ̃gg′-(rΓ,rΓ′)Jg′-(rΓ′)dΓ′+14∑g′=1G∫ΩΨ̃gg′(rΓ,r′)qg′(r′)dΩ′,g=1,2,…,G,
where(38)Jg±(rΓ)=14ϕ(rΓ)∓12Dg∂ϕ(rΓ)∂nΓ,J̃gg′±(rΓ,rΓ′)=14Dg′Ψ̃(rΓ,rΓ′)±12∂Ψ̃(rΓ,rΓ′)∂nΓ′.


## 4.2. Reduction of the Domain Integrals of the Source Term

The domain integrals involving the independent and flux dependent fission sources in (37) can be reduced into a boundary integrated form as(39)∫ΩΨ̃gg′(rΓ,r′)qg′(r′)dΩ′=∑h=1G-σghσhg′Bh2×{c(rΓ)qg′(rΓ)+∫Ωψ̃h(rΓ,r′)∇2qg′(r′)dΩ′+∫Γ[∂ψ̃h(rΓ,rΓ′)∂nΓ′qg′(rΓ′)-ψ̃h(rΓ,rΓ′)∂qg′(rΓ′)∂nΓ′]dΓ′}=∑h=1G-σghσhg′Bh2{c(rΓ)qg′(rΓ)ψ̃h∂nΓ′+∫Γ[∂ψ̃h(rΓ,rΓ′)∂nΓ′qg′(rΓ′)-ψ̃h(rΓ,rΓ′)∂qg′(rΓ′)∂nΓ′]dΓ′}h=1,2,…,G.
The source distributions qg(r) are now assumed to be either spatially constant or first-order polynomials. In the latter case for 2D geometry the following expression is assumed:(40)q(x,y)=l1x+l2y+q̅,
where, obviously, q̅ is the average source on the domain x∈[0,a] and y∈[0,b]. The following relationship holds:(41)qg′(x,y)=l1x+l2y+q̅g′-(al12+bl22).
As q̅ in the presence of fission emissions is evaluated by the flux of the previous generation on the boundary, the coefficients of expansion (40) can be evaluated by using a least-square approach. The problem can be nicely written in matrix form(42)Ml=N,
having set the following definitions:(43)M=(x1-Δx2y1-Δy2x2-Δx2y2-Δy2⋮⋮xN-Δx2yN-Δy2),l=(l1l2),N=(q1-q̅q2-q̅⋮qN-q̅),
where (xi,yi) are points on the boundary at which the source takes the value qi. By left multiplication of (42) by MT, the normal equations are obtained for the problem
(44)MTMl=MTN,
or, explicitly(45)(∑n=1N(xn-a2)2∑n=1N(xn-a2)(yn-b2)∑n=1N(xn-a2)(yn-b2)∑n=1N(yn-b2)2)(l1l2)=(∑n=1N(xn-a2)(qg,n-q̅g)∑n=1N(yn-b2)(qg,n-q̅g)).
The above system of equations yields the value of the two-dimensional vector l.


## 4.3. Discretization of the Boundary Integral Equations

In this section the numerical solution of (37) is outlined. The boundary Γ is partitioned into I elements of length di. After denoting the two extremes and the middle point by rΓ,i-1/2, rΓ,i+1/2, and rΓ,i, respectively, a local dimensionless coordinate τ is introduced to denote each point r along Γi as(46)r=rΓ,i+τ(rΓ,i+1/2-rΓ,i-1/2).
Thus, it is assumed that(47)Jg±(rΓ)=Jg,i±.
Then, assuming a linear behaviour on each boundary mesh of the integrand functions, the boundary integrals are approximated by finite sums as(48)∫ΓJ̃gg′±(rΓ,i,rΓ′)Jg′±(rΓ′)dΓ′≅∑j=1I∫-11J̃gg′±(rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2))Jg′,j±dj2dτ.
By indicating(49)qg(rΓ)=qg,i,∂qg(rΓ)∂nΓ=qg,i′,
one has(50)∫Γ[∂ψ̃h(rΓ,rΓ′)∂nΓ′qg′(rΓ′)-ψ̃h(rΓ,rΓ′)∂qg′(rΓ′)∂nΓ′]dΓ′≅∑j=1I∫-11{∂ψ̃h∂nΓ′[rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2)]qg′,j-ψ̃h[rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2)]qg′,j′∂ψ̃h∂nΓ′}dj2dτ.
The integrals herewith appearing are obtained by application of a numerical integration formula, namely(51)∫-11f(τ)dτ≃∑k=1Kωkf(τk),
where ωk and τk are the quadrature weights and abscissae, respectively. In the present paper a standard Gauss-Legendre formula is used.At last the problem is cast into the following algebraic response-matrix form:(52)∑g′=1G∑j=1IMgg′,ijJg′j+=∑g′=1G∑j=1INgg′,ijJg′j-+∑g′=1G∑j=1IKgg′,ijqg′j+∑g′=1G∑j=1ILgg′,ijqg′j′,i=1,2,…,I,g=1,2,…,G,
where(53)Mgg′,ij=c(r)2δi,jδg,g′+dj2∫-11J̃gg′-[rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2)]dτ,Ngg′,ij=-c(r)2δi,jδg,g′+dj2∫-11J̃gg′+[rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2)]dτ,Kgg′,ij=-∑h=1Gσgh14Bh2σhg′*×{c(r)δi,j+dj2∫-11∂ψ̃h∂nΓ′[rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2)]dτ},Lgg′,ij=∑h=1Gσgh14Bh2σhg′*dj2∫-11ψ̃h[rΓ,i,rΓ,j′+τ(rΓ,j+1/2-rΓ,j-1/2)]dτ.
The problem for the whole system where nodes are coupled is then solved by numerical iteration schemes.


## 5. Verification of Numerical Consistency and Convergence of Methods

It is worth going through a verification of the performance of various numerical approaches by comparison of results with highly accurate benchmark results that may be obtained directly by an analytical formula. Therefore a homogeneous rectangular system is now considered, whose dimensions are 160 cm and 140 cm. The material data for a three-group model are reported in Table1. A critical calculation is carried out and errors with respect to the analytical value of the effective multiplication constant are reported in Table 2. The various models are identified as follows: FD indicates the finite-difference scheme, N is nodal, BEM αβγ denotes a distribution of order α for the source (0, constant, or 1, linear, in the present work), β the number of points used for the discretization of each edge, and γ the order of the quadrature formula for integration. The number reported in the node column indicates the number of nodes assumed for each coordinate; hence the number of nodes covering the 2D domain increases as the square of the indicated number. The performance of the nodal method is quite remarkable, while good results are obtained by BEM only with a first-order approximation of the source distribution. Also, it is important to call the attention on the important role associated to the number of points used for the discretization of the boundary.Table 1
Material data for the homogeneous core.
gΣrgνΣfgDgΣsg→g+1Σsg→g+2χg(cm-1)(cm-1)(cm)(cm-1)(cm-1)—10.3046E-010.1076E-012.110.2578E-010.5164E-030.773720.7060E-020.3309E-021.290.4705E-02—0.219330.6316E-020.4366E-020.86——0.0070Table 2
Errors in the calculation of the effective multiplication factor for the homogeneous benchmark; the analytical value isk= 0.90145. Relative errors are expressed in pcm for the various models.
NodesFDNBEMBEMBEMBEM0 10 241 10 241 10 2401 20 2404929189070494495498823512654212633121030128393516570787149020370554161022425042617112281803491811232140299181223611026418122408023918123Furthermore, it is interesting to observe the convergence pattern for the methods studied, see Figure1. While FD shows a monotonic behaviour, Nodal and some BEM methods are characterized by nonmonotonic trends and reach a saturation. Such properties for the convergence of nodal methods are well-known to exist [12] and would deserve deeper investigation. They are connected to the assumption of the spatial distribution within the node, which, in general, may not be consistent with the physical model as described by the differential equations to be solved.Figure 1
Convergence pattern for the calculation of the homogeneous system. FD (○); N (□); BEM 0 10 24 (
); BEM 1 10 24 (*); BEM 1 10 240 (•); BEM 1 20 240 (+).Some final considerations concerning the computational effort are in order. Considering Figure2, it is clear that to obtain relative errors below 10 pcm nodal and BEM are by far more effective.Figure 2
Relative errors versus computational times for the homogeneous reactor calculation. Identification of curves as in the previous figure.A nonhomogeneous system is also considered, to test the performance of the schemes in the presence of a spatial interface. The symmetrical core is 180 cm wide and characterized by the data reported in Table1. The data for the 90 cm thick reflector are specified in Table 3. The system is reflected only in the x-direction. The number of nodes herewith indicated refers to the discretization in the x-direction only, since along y only two nodes are used with zero-current on the boundary condition, to simulate an infinite medium. Results concerning k values are reported in Table 4. A very fine finite-difference solution (320 nodes), to be assumed as reference, yields the value k=1.04907. The enhanced performance of nodal and boundary element methods is clearly visible.Table 3
Cross-sections for the reflector in the heterogeneous case.
gΣrgDgΣsg→g+1Σsg→g+2(cm-1)(cm)(cm-1)(cm-1)10.9667E-022.520.9506E-020.1270E-0320.1508E-021.560.1385E-02—30.1407E-031.00——Table 4
Results for the calculation ofk for the reflected system.
NodesFDNBEMBEMBEM1 10 241 10 2401 20 24041.040261.049231.049051.049141.0491681.043731.049241.049061.049211.04922121.045321.049241.049011.049211.04922161.046201.049241.048951.049211.04922201.046761.049241.048891.049201.04922241.047151.049241.048831.049191.04922281.047431.049241.048771.049191.04921321.047641.049241.048711.049181.04921361.047811.049241.048651.049181.04921401.047951.049241.048591.049171.04920


## 6. Results for a Fast-Spectrum Configuration

A realistic fast-spectrum system is considered, with physical characteristics of the same type as the ones characterizing the ELSY project. The core is constituted by three regions, with increasing multiplicativities passing from the inner to the outermost regions. Lead is assumed for the coolant and the reflector. Cross-sections are provided by spectrum and homogenization calculations carried out by the ENEA group participating to the ELSY project. The sketch of the system and the nuclear data for the core regions and for the reflector are reported in Figure3 and Table 5. Although only three energy groups are used, the test is useful to verify the performance of the computational schemes in a highly heterogeneous system where spectral effects are important.Table 5
Material data for the ELSY-like reactor.
gΣrgνΣfgDgΣsg→g+1Σsg→g+2χg(cm-1)(cm-1)(cm)(cm-1)(cm-1)—CORE110.3303E-010.1076E-012.110.2578E-010.5164E-030.773720.8638E-020.3309E-021.290.47052E-02—0.219330.7369E-020.4366E-020.86——0.0070CORE210.3306E-010.1146E-012.110.2562E-010.5135E-030.773720.8651E-020.3698E-021.300.4576E-02—0.219430.7460E-020.4766E-020.86——0.0069CORE310.3314E-010.1267E-012.110.2531E-010.5071E-030.773520.8711E-020.4551E-021.310.4324E-02—0.219530.7675E-020.5629E-020.87——0.0070REFLECTOR10.1274E-01—2.520.9506E-020.1270E-03—20.3413E-02—1.560.1385E-02——30.1358E-02—1.00———Figure 3
Geometrical configuration of the ELSY-like system. The transverse lines along which graphs are plotted in the following figures are indicated.As already anticipated, adjoints may well serve in many interesting and useful applications in nuclear core analysis. Once the physical problem is cast into an algebraic form, the adjoint problem is easily constructed by taking the mathematical adjoint of the matrix operators appearing in the model. It is then possible to compute neutron importances that may be used for sensitivity and perturbation analyses and for the construction of integral parameters for time-dependent calculations. The results of some such calculations are now presented in Figures4, 5, 6, 7, 8, and 9. For brevity, only results for the nodal method are presented, having previously proved its consistency with BEM.Figure 4
Graphs of the first-group (high energy) fluxes. Symbols on curves identify fluxes plotted along the transverses bearing the same symbol indicated in Figure3.Figure 5
Graphs of the second group (intermediate energy) fluxes.Figure 6
Graphs of the third group (low energy) fluxes.Figure 7
Graphs of the neutron importance for the first group.Figure 8
Graphs of the neutron importance for the second group.Figure 9
Graphs of the neutron importance for the third group.


## 7. Conclusions and Future Developments

A neutronic module for multidimensional multigroup diffusion calculations is developed. Different discretization algorithms are implemented and tested. Some benchmark calculations allow to assess the consistency and applicability of various techniques. Nodal methods are shown to perform adequately for the simulation of fast spectrum systems, such as lead-cooled reactors. Calculations are carried out for a system characterized by properties similar to the system adopted as the basis of the ELSY project.As a future development, it is foreseen to couple the neutronic module developed in the frame of this work with a thermohydraulic code. This would allow to have a flexible and efficient computational tool for the assessment of the ELSY project and for parametric calculations aimed at the neutronic-thermal-hydraulic optimization. Furthermore, the treatment of full three-dimensional Cartesian geometries is envisaged.A further future development will consider the extension of the computational tool to allow time-dependent evaluations within a quasi-static approach. To that end, the module can perform already adjoint calculations. The extension is to be done following the same philosophy as the one already employed for other recent applications [13].



---

*Source: 973605-2009-09-13.xml*