Neutron noise simulation by GFEM and unstructured triangle elements

Nuclear Engineering and Design 253 (2012) 238–258

Contents lists available at SciVerse ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

Neutron noise simulation by GFEM and unstructured triangle elements Seyed Abolfazl Hosseini, Naser Vosoughi ∗ Department of Energy Engineering, Sharif University of Technology, Tehran 8639-11365, Iran

h i g h l i g h t s     

We develop a 2-D, 2-group neutron noise simulator based on GFEM. The spatial discretization is performed using unstructured triangle elements. The effects of fluctuations in diffusion coefficient on results are investigated. The developed code can be applied for both rectangular and hexagonal reactor cores. We find that the developed code is reliable tool for neutron noise calculation.

a r t i c l e

i n f o

Article history: Received 28 May 2012 Received in revised form 8 August 2012 Accepted 9 August 2012

a b s t r a c t In the present study, the neutron noise, i.e. the stationary fluctuation of the neutron flux around its mean value, is calculated in 2-group forward and adjoint diffusion theory for both hexagonal and rectangular reactor cores. To this end, the static neutron calculation is performed at the first stage. The spatial discretization of equations is based on linear approximation of Galerkin Finite Element Method (GFEM) using unstructured triangle elements. Using power iteration method, forward and adjoint fluxes with the corresponding eigenvalues are obtained. The results are then benchmarked against the valid results for BIBLIS-2D and IAEA-2D benchmark problems and DONJON computer code. The dynamic calculations are performed in the frequency domain which leads to reducing the dimension of the variable space of the noise equations. The forward/adjoint noise in two energy group is obtained by assuming the neutron noise source as an absorber of variable strength type. The neutron noise induced by a vibrating absorber type of noise source is also obtained from the calculated adjoint Green’s function. Comparison of the calculated neutron noise at zero frequency with the results of static calculation, and utilizing the results of adjoint noise calculations are two different procedures to validate the neutron noise calculations. © 2012 Elsevier B.V. All rights reserved.

1. Introduction The knowledge of the neutron noise, i.e. the difference between the time-dependent neutron flux and its steady state value with the assumptions that the processes are stationary and ergodic in time, is very useful for core monitoring, surveillance and diagnostics (Arzhanov and Pázsit, 2002). Noise analysis method might be applied to determine dynamical parameters of the core, like the Moderator Temperature Coefficient (MTC) in a Pressurized Water Reactor (PWR) (Demazière and Pázsit, 2004) or the Decay Ratio in a Boiling Water Reactor (BWR) (Arzhanov and Pázsit, 2002). It can be applied to identify and localize the noise source such as unseated fuel assemblies, absorber of variable strength, vibrating absorber or vibrations of core internals (Demazière, 2004; Demazière and Pázsit, 2008; Williams, 1974). One of the main advantages of the noise-based methods is that they can be used

∗ Corresponding author. Tel.: +98 21 6616 6130; fax: +98 21 6608 1723. E-mail addresses: [email protected], [email protected] (N. Vosoughi). 0029-5493/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.nucengdes.2012.08.023

on-line without disturbing reactor operation (Sunder et al., 1991). Such monitoring techniques are of interest for the extensive program of power upgrades around the world. The other important advantage of noise analysis is that the calculations of neutron noise are performed in the frequency domain. This leads to reducing the dimension of the variable space of the noise equations. Namely, instead of solving the diffusion equations in both the time and space, after the Fourier transform one needs to solve them only in space whereas the frequency variable acts as a parameter in such a case (i.e. no any derivatives with respect to frequency). As the system of noise equations is not an eigenvalue problem, it is solved easier than the static ones. To calculate the neutron noise, firstly the neutron noise induced by a point-like noise source has to be determined. Although many commercial codes for static calculations are available, there is no dedicated commercial code that is able to calculate the dynamical reactor transfer function. Recently, several researchers have tried to develop convenient methods for calculations of the dynamical reactor transfer function and noise analysis. Some selected researches are presented in the following:

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

Nomenclature FEM g (r) † g (r) keff

Finite Element Method forward flux in the energy group g adjoint flux in the energy group g forward neutron multiplication factor

keff g ı(r, ω) ı† (r, ω) Dg ˙r,g

adjoint neutron multiplication factor neutron spectrum in energy group g forward noise adjoint noise diffusion constant in the energy group g macroscopic removal cross section in the energy group g macroscopic fission cross section in the energy group g macroscopic scattering cross section from energy group g  to g fission neutron yield the nabla operator

†

˙f,g ˙g  →g 

∇

1. Demazière (2004) proposed a 2-D, 2-group neutron noise simulator, which calculates the neutron noise and the corresponding adjoint function in the 2-D heterogeneous systems. The simulator can only be applied to the western LWRs (PWR and BWR cases). Since the used spatial discretization based on finite difference scheme was derived for the rectangular mesh boxes, the VVER-type, some of the Gen-IV reactor cores and totally all the hexagonal-structured reactor cores are not included in its applications. 2. Malmir et al. (2010) reported the development of a 2-D, 2-group neutron noise simulator for hexagonal-structured reactor cores using both the forward and the adjoint methods. The spatial discretization of both the static and dynamical equations was based on a developed box-scheme Finite Difference Method (FDM) for hexagonal mesh boxes. 3. Larsson et al. (2011) calculated neutron noise in 2-group diffusion theory using the Analytical Nodal Method (ANM). The same spatial discretization scheme was applied for the static calculation. 4. Demazière (2011) has developed a multi-purpose neutronic tool, in which both critical and subcritical systems with an external neutron source can be studied, and static and dynamic cases in the frequency domain can be considered. This tool has the ability to determine the different eigenfunctions of any nuclear core. For each situation, the static neutron flux, the different eigenmodes and eigenvalues, the first order neutron noise, and their adjoint functions are estimated. 5. Larsson and Demazière (2012) have developed a coupled neutronics/thermal-hydraulics tool for calculating fluctuations in PWRs. The fluctuations in neutron flux, fuel temperature, moderator density and flow velocity in PWRs are estimated by coupling a dynamic thermal hydraulic module and a dynamic neutron kinetic module. There is also a possibility to directly define the perturbations in the macroscopic cross sections and to supply them to the neutronic part of the model. 6. Tran and Demazière (2012) have developed a neutronic and kinetic solver for hexagonal geometries based on the diffusion theory with multi-energy groups and multi-groups of delayed neutron precursors. The tool is utilized to forward/adjoint problems of static and dynamic states for both thermal and fast system with hexagonal geometries. The spatial discretization of the equations is based on FDM.

239

In the above-mentioned researches, the spatial discretization of noise equations has been performed based on FDM or ANM. Other numerical methods like Finite Element Method (FEM), boundary element or finite volume techniques, Direct Discrete Method (DDM) (Vosoughi et al., 2004; Ayyoubzadeh et al., 2012) may also be used as the spatial discretization scheme. In general, FEM is preferred in most applications to its principal alternative, FDM, due to its flexibility in the treatment of curved or irregular geometries and the high rates of convergence attainable by the use of high order elements. The first application of FEM to the theory of neutron diffusion dates back to 1970s (Kang and Hansen, 1971). The developments in the application of FEM to the neutron diffusion equation have been described in the excellent treatise of Lewis (1981). Recently, several other applications of FEMs including Raviart–Thomas–Schneider, Hybrid, h-adaptivity, Response Matrix, etc. to solve neutron diffusion equation have been introduced (Hebert, 2008; Cavdar and Ozgener, 2004; Wang et al., 2009). In the present study, the Galerkin FEM (GFEM) (Hosseini and Vosoughi, 2012), a weighted residual method, is used to solve the 2-group forward/adjoint diffusion equation in BIBLIS-2D (Maiani and Montagnini, 2004; Nakata and Martin, 1983; Hebert, 1985) and IAEA-2D (Gonzalez et al., 2009; Hebert, 2008) reactor cores. The forward/adjoint equation is solved using a solution algorithm based on a Galerkin-type scheme by considering the linear approximation of a shape function (basic function applied to construct an approximation to the solution of equation). Same spatial discretization scheme is used to solve the forward/adjoint noise equations. Unstructured triangle elements generated by Gambit (FDI, 2000) are used to discretize the equations. Indeed, a key advantage of the unstructured triangular elements is their superiority in mapping the curved boundaries or material interfaces. An outline of the remainder of this contribution is as follows: In Section 2, we briefly introduce the mathematical formulation used to static calculations. In Section 3, the method of dynamical calculations is described. Section 4 presents the main specification of the BIBLIS-2D and IAEA-2D benchmark problems. The results of static and dynamical calculations and discussion on the results are given in Section 5. Section 6 gives a summary and concludes the paper.

2. Static calculations 2.1. Solution of the forward diffusion equation In the absence of an external neutron source, the multigroup forward diffusion equation is as Eq. (1) (Duderstadt and Hamilton, 1991; Lamarsh, 1965):

−Dg ∇ 2 g (r) + ˙r,g g (r) =

+



g keff

G 

˙f,g  g  (r)

g  =1

˙s,g  →g g  (r) g = 1, 2, . . . , G.

(1)

g = / g

where all quantities are defined as usual. Removal cross sections are expressed as ˙r,1 = ˙a,1 + ˙s,1→2 and ˙r,2 = ˙a,2 . Eq. (1) is a linear partial differential equation which may be solved by different numerical methods. All of the methods transform the differential equation into a system of algebraic equations. Here, GFEM is used to discretize the forward diffusion equation. To start the discretization, the whole solution area is divided into the unstructured triangle elements as shown in Fig. 1. The elements have been generated using Gambit mesh generator. In the

240

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

Fig. 1. The domain including the unstructured triangle finite elements.

Fig. 2. Transforming each of triangle elements (˝(e) ) (physical domain) to the calculation domain.

Fig. 3. The reactor core of the BIBLIS-2D.

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

241

Fig. 6. Calculated thermal forward flux distribution using GFEM for BIBLIS-2D. Fig. 4. The reactor core of the IAEA-2D.

and linear approximation of the shape function, the forward flux in each element could be considered as Eq. (2) (Zienkiewicz et al., 2005): (e)

(e)

(e)

(e) (x, y) = Ni (x, y)i + Nj (x, y)j + Nk (x, y)k , (e)

(e)

(2)

(e)

where Ni , Nj and Nk are the components of the shape function. In the linear approximation, these components are equal to the (e) corresponding Ln (x, y) as Eq. (3): (e)

Ln (x, y) =

an + bn x + cn y 2(e)

;

n = i, j, k,

(3)

In which ai = xj yk − yj xk ;

bi = yj − yk ;

ci = xk − xj ,

(4)

aj = xk yi − yk xi ;

bj = yk − yi ;

ci = xi − xk ,

(5)

ak = xi yj − yi xj ;

bk = yi − yj ;

ck = xj − xi ,

(6)

Fig. 5. Calculated fast forward flux distribution using GFEM for BIBLIS-2D.

   1 xi yi      2(e) =  1 xj yj  .   1 x y  k k

(7)

The components of the shape function (N (e) (r)) satisfy the criterion given in Eq. (8) at all points of the domain: (e)

(e)

(e)

Ni (x, y) + Nj (x, y) + Nk (x, y) = 1.

(8)

The GFEM is a weighted residual method in which the purpose is to minimize the residual integral. In the weighted residual methods, the weighting function is considered as W (r) = W T N(r). There are (at least) four sub-methods (Collocation, Sub-domain, Least Squares, Galerkin) according to choices for W T . Since W T is unit in the Galerkin method, the weighting function is considered as

Fig. 7. Calculated fast adjoint flux distribution using GFEM for BIBLIS-2D.

242

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

In the above equation, the differential part may be transformed by applying the Green’s theorem:





d˝W (r)(−Dg ∇ 2 g (r)) = ˝

d˝∇ W (r) · ∇ g (r) ˝



 d˝∇ · (W (r)∇ g (r)) =

−



˝



ds.W (r)∇ g (r) =

−

d˝∇ W (r) · ∇ g (r) ˝

d˝∇ W (r) · ∇ g (r)

∂˝ ˝

 −

dsW (r) ∂˝

∂g (r) . ∂n

(11)

where ∂g (r) = ∇ g (r) · n. ∂n

(12)

Here, n is the normal unit vector on the surface ˝. Two types of boundary conditions (B.C.) are considered. The first B.C. is no incoming neutrons at vacuum boundaries (Marshak B.C.) which is expressed as Eq. (13):

Fig. 8. Calculated thermal adjoint flux distribution using GFEM for BIBLIS-2D.

Eq. (9): W (r) = N(r),

(9)

where N(r) is the global shape function. Multiplying Eq. (1) by the weighting function and integrating the results over the solution space, Eq. (10) is obtained:

 d˝ W (r)(−Dg ∇ 2 g (r) + ˙r,g g (r) − ˝

−



˙s,g  →g h (r)) = 0.

g keff

G 

∂g (r) g (r) . =− 2Dg ∂n

The second B.C. is zero net current or perfect reflective boundary condition which is described by Eq. (14): ∂g (r) = 0. ∂n

˙f,g  g  (r)

g  =1

(10)

(13)

(14)

Substituting Eqs. (9), (11), (13) and (14), and converting the integration on the reactor domain to sum of the integration on finite

g = / g

Fig. 9. Comparison of the calculated and reference power distribution for BIBLIS-2D.

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

243

Fig. 10. Comparison of calculated fast forward flux distribution using DONJON and GFEM for BIBLIS-2D.

elements, the final form of Eq. (10) is obtained as Eq. (15):

  

=

e=1



E

e=1

E 

(e)

(e)

d˝ Dg ∇ N (e) (r) · ∇ N (e)T (r)g + ˙r,g ˝(e)

˝(e)

 (e)

×d˝ N (e) (r)N (e)T (r)g +

d˝ N (e) (r)N (e)T (r) ˝(e)

(e) g



+

G 

g

= / g

⎡ ⎣ g

keff

G 

 (e)

(e)

d˝ N (e) (r)N (e)T (r)g 

˙f,g 

g

˝(e)

 (e) ˙s,g  →g

(e) d˝ N (e) (r)N (e)T (r) g  ˝(e)

2

Fig. 11. Comparison of calculated thermal forward flux distribution using DONJON and GFEM for BIBLIS-2D.

⎤ ⎦.

(15)

244

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

Fig. 12. Comparison of calculated fast adjoint flux distribution using DONJON and GFEM for BIBLIS-2D.

Fig. 13. Comparison of calculated thermal adjoint flux distribution using DONJON and GFEM for BIBLIS-2D.

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

245

Fig. 14. Calculated fast forward flux distribution using GFEM for IAEA-2D.

Fig. 16. Calculated fast adjoint flux distribution using GFEM for IAEA-2D.

To calculate the integrals in Eq. (15), a change of variables is used to transform each of triangle elements (˝(e) ) (physical domain) to the calculation domain, as shown in Fig. 2. In the calculation domain, the variables  (e) and (e) are defined in terms of old variables x(e) and y(e) as Eqs. (16) and (17):

Therefore, the integral of an arbitrary function f (x, y) is obtained from Eq. (21):

 (e) (x, y) =

1

(e)

1

(e) (x, y) =

(e)

(e)

(e)

(e)

(e)

[(y3 − y1 )(x − x1 ) + (x1 − x3 )(y − y1 )], (16)

2(e)

(e)

(e)

(e)

(e)

(e)

(e)

[(y1 − y2 )(x − x1 ) + (x2 − x1 )(y − y1 )]. (17)

2(e)

By some manipulations on the components of the shape function and Eqs. (16) and (17), the shape function components are obtained in the calculation domain as Eqs. (18)–(20): (e) L1 (e) L2 (e)

=1−

(e)

=  (e) =

(e)

−

=

(e) N1 ,

(18)

(e) N2 ,

(19)

(e)

(20)

L3 =  (e) = N3 .



 1 f (x, y)dx dy = ˝(e)

1− (e)

f (x(, ), y(, ))|J|d d , 0

where |J| is the Jacobian determinant according to Eq. (22):

  ∂x(e)   ∂(e) |J| =  (e)  ∂y  ∂(e)

∂x(e) ∂ (e) ∂y(e) ∂ (e)

    .   

(22)

Evaluating the integrals of Eq. (15) for each element, using Eq. (21), and assembling the local matrix into the global matrix, the system of equations which is an eigenvalue problem is obtained. Here,

Fig. 17. Calculated thermal adjoint flux distribution using GFEM for IAEA-2D. Fig. 15. Calculated thermal forward flux distribution using GFEM for IAEA-2D.

(21)

0

246

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

Fig. 18. Comparison of the calculated and reference power distribution for IAEA-2D.

Fig. 19. Comparison of calculated fast forward flux distribution using DONJON and GFEM for IAEA-2D.

Fig. 20. Comparison of calculated thermal forward flux distribution using DONJON and GFEM for IAEA-2D.

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

247

Fig. 21. Comparison of calculated fast adjoint flux distribution using DONJON and GFEM for IAEA-2D.

Fig. 22. Comparison of calculated thermal adjoint flux distribution using DONJON and GFEM for IAEA-2D.

Fig. 24. The arrangement of unstructured triangle elements in the IAEA-2D.

Fig. 23. The arrangement of unstructured triangle elements in the BIBLIS-2D.

248

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

Fig. 25. Fast forward noise distribution at zero frequency in the BIBLIS-2D.

the problem is solved using the power iteration method described as follow: (1) Initialize the multigroup fluxes and the eigenvalue with (0) (0) (0) g , keff ; compute the initial fission source Sf (r) =

G

(0)

g=1

˙f,g (r)g (r) and set the iteration index to n = 1.

Fig. 27. Fast forward noise distribution at zero frequency in the IAEA-2D.

(4) Update

(n−1)

keff

the (n)

˝

d˝Sf (r)/ (n) keff



eigenvalue (n−1)

d˝Sf

˝ (n−1) keff

as:

(n)

keff =

(r). (n)

(n−1)

(5) Compare with and g with g for all energy groups. If the changes are greater than a prescribed tolerance, then set n = n + 1 and repeat the iteration stating at step-(2); otherwise the iteration is over.

(n)

(2) Solve for the multigroup fluxes g using: (n)

(n)

−Dg ∇ 2 g (r) + ˙r,g g (r) = +



(n−1)

˙s,g  →g g 

g (n−1)

keff

Applying the above mentioned algorithm to the system of equations results to the calculation of the forward flux in each energy group and the multiplication factor.

(n−1)

Sf

(r) g = 1, 2, ..., G.

(23)

g = / g

 G

(n)

(3) Update the fission integral as: Sf (r) =

(n)

˙f,g (r)g (r).

g=1

Fig. 26. Thermal forward noise distribution at zero frequency in the BIBLIS-2D.

2.2. Solution of the adjoint diffusion equation For a point detector and a point source, both the forward and the adjoint are equal. If either the source or the detection volume is finite, then only one of the two equations is applicable directly. For the flux in a point with an extended source (i.e. with a point-like

Fig. 28. Thermal forward noise distribution at zero frequency in the IAEA-2D.

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

Fig. 29. Magnitude of thermal neutron noise due to noise source located at X = 195 cm, Y = 195 cm in BIBLIS-2D.

Fig. 30. Phase of thermal neutron noise due to noise source located at X = 195 cm, Y = 195 cm in BIBLIS-2D.

Fig. 31. Magnitude of thermal neutron noise due to noise source located at X = 250 cm, Y = 195 cm in BIBLIS-2D.

249

250

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

Fig. 32. Phase of thermal neutron noise due to noise source located at X = 250 cm, Y = 195 cm in BIBLIS-2D.

Fig. 33. Magnitude of thermal neutron noise due to noise source located at X = 310 cm, Y = 195 cm in BIBLIS-2D.

Fig. 34. Phase of thermal neutron noise due to noise source located at X = 310 cm, Y = 195 cm in BIBLIS-2D.

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

Fig. 35. Magnitude of the neutron noise due to the absorber of variable strength source given by the forward calculation at frequency 1 Hz for BIBLIS-2D.

Fig. 36. Phase of the neutron noise due to the absorber of variable strength source given by the forward calculation at frequency 1 Hz for BIBLIS-2D.

Fig. 37. Magnitude of the neutron noise due to the absorber of variable strength source given by the adjoint calculation at frequency 1 Hz for BIBLIS-2D.

251

252

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

Fig. 38. Phase of the neutron noise due to the absorber of variable strength source given by the adjoint calculation at frequency 1 Hz for BIBLIS-2D.

Fig. 39. Magnitude of the neutron noise due to the absorber of variable strength source given by the forward calculation at frequency 1 Hz for IAEA-2D.

Fig. 40. Phase of the neutron noise due to the absorber of variable strength source given by the forward calculation at frequency 1 Hz for IAEA-2D.

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

253

Fig. 41. Magnitude of the neutron noise due to the absorber of variable strength source given by the adjoint calculation at frequency 1 Hz for IAEA-2D.

or ı function detector) the forward equation is preferred. Likewise, for the point source and extended detector, the adjoint equation is more useful (Pázsit, 2007). To solve the adjoint diffusion equation, it is noted that adjoint operator is the transpose of the forward operator in 2-group diffusion equation (Duderstadt and Hamilton, 1991; Bell and Glasstone, 1970). To this end, Eq. (1) is formed in a matrix notation as Eq. (24): L =

1 F, keff

(24)

1 † keff

F † † ,

3.1. Calculations of forward noise Here, the first-order forward noise approximation of two-group forward diffusion equation is applied to noise calculation. The global form of two-group noise equation which the neutron noise sources is due to the variations of scattering, absorption and fission macroscopic cross sections, is given in Eq. (26) (Demazière, 2004):

where L is called as a loss operator and F as a fission operator. Therefore, the matrix form of the adjoint diffusion equation is written as Eq. (25): L † † =

3. Dynamical calculations

(25)

here L† and F † are the transpose of the L and F, respectively † (Lamarsh, 1965). Also, † refers to adjoint flux, and keff is the adjoint multiplication factor. To solve the adjoint diffusion equation, the same method applied to discretization of forward diffusion equation is used. As the results of calculations, the adjoint flux in each energy group and adjoint multiplication factor are obtained.

[∇ · D(r)∇ + ˙ dyn (r, ω)] ×

+ a (r)

ı˙a,1 (r, ω) ı˙a,2 (r, ω)



ı1 (r, ω) ı2 (r, ω)

+ f (r, ω)

 = s,1→2 (r)ı˙s,1→2 (r, ω)

ı1 ˙f,1 (r, ω)

 ,

(26)

ı2 ˙f,2 (r, ω)

where the above mentioned matrixes and vectors are expressed as Eqs. (27)–(31):

D=

D1 (r)

0

0

D2 (r)

 ,

Fig. 42. Phase of the neutron noise due to the absorber of variable strength source given by the adjoint calculation at frequency 1 Hz for IAEA-2D.

(27)

254

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

⎡ ⎢ ⎣

˙ dyn (r, ω) = ⎢

s,1→2 (r) =

a (r) =



1 (r)

1 (r)

⎤

iωˇeff iω

E   

⎥ ⎥, ⎦

iω +

− ˙a,2 (r) +





(28)



iωˇeff

1−



 −2 (r)

iω +

1−

iωˇeff iω +

iω

v1



1 ˙f,1 (r)

1−

keff

iωˇeff

⎦ . (31)



iω +

.

[∇ .D(r)∇ + ˙ dyn (r, ω)] × ı(r − r  )



(32)

0

g=1

(r, r  , ω)

.

(33)

g=2

S1 (r  , ω)

S(r  , ω) =

+ a (r)



S2 (r  , ω)

= s,1→2 (r  )ı˙s,1→2 (r  , ω)

ı˙a,1 (r  , ω)



+ f (r  ω)

ı˙a,2 (r  , ω)

ı1 ˙f,1 (r  , ω) ı2 ˙f,2 (r  , ω)

 .

(34)

The same applied discretization method for the static calculation is used to the dynamical calculation. By considering the neutron noise source in group 2, the discrete form of Eq. (33) is obtained as Eqs. (35) and (36): E   

(e)

keff



1−

iωˇeff

  (e)

iω +

˝(e) (e)

ds N ∂˝(e)V

⎥ ⎦

= ⎢ N (e) (r) ⎥ , j

(36)

where ∂˝(e)V refers to boundary length with vacuum boundary condition in element e. The analytical solution of the integrals in Eqs. (35) and (36) is given in Table 1. Green’s function components in each energy group are obtained from the calculation of Eqs. (35) and (36). Finally, the fast and thermal forward noises are calculated by an integral over the source and the Green’s function as Eq. (37):



ı1 (r, ω)



ı2 (r, ω)

⎡ ⎢

=⎣

 [G1→1

(r, r  , ω)S

1

(r  , ω)]

+



[G1→2 (r, r  , ω)S1 (r  , ω)] +

⎤ [G2→1 (r, r  , ω)S2 (r  , ω)]dr  [G2→2 (r, r  , ω)S2 (r  , ω)]dr 

⎥. ⎦

In this study, it is assumed that the thermal macroscopic absorption cross section is perturbed. In such a case, there is only thermal neutron noise source (Demazière, 2004) and Eq. (37) reduces to Eq. (38): ı1 (r, ω) ı2 (r, ω)



⎡ ⎢

= ⎣

[G2→1 (r, r  , ω)S2 (r  , ω)]dr  [G2→2 (r, r  , ω)S2 (r  , ω)]dr 

⎤

⎥ ⎦.

(38)

3.2. Calculations of adjoint noise As for the static calculation, the estimation of the neutron noise and the corresponding adjoint function has some practical applications. Hence, it is highly important to calculate the adjoint function of the neutron noise. Considering a point-like neutron noise source located at the position r0 , the adjoint noise equations could be derived by the same process used to obtain the forward noise equations. The matrix form of adjoint noise equations might be expressed as Eq. (39) (Demazière, 2004):



†



∇ · D(r)∇ + ˙ dyn (r, ω) ×

†

ı1 (r, ω) †

ı2 (r, ω)



=

ı(r − r0 ) ı(r − r0 )



(39)

†

matrix ˙ dyn (r, ω). In addition, ı(r − r0 ) refers to Dirac delta function at the position r0 . Like the forward noise calculations, Green’s function technique and GFEM are applied to adjoint noise calculations. The fast and thermal adjoint noises are the results of these calculations.

d˝ N (e) (r)N (e)T (r)G2→1

 +

⎤

†

(e)

˝(e)

2 ˙f,2

⎢ ⎣

(e)

Ni (r)

noises, respectively. ˙ dyn (r, ω) is the transpose of the dynamical

(e)

+

2

⎡



(e)

†

d˝ D1 ∇ N (e) (r)∇ N (e)T (r)G2→1



−˙1

ds N (e) (r)N (e)T (r)

(e) G2→2

where ı1 (r, ω) and ı2 (r, ω) are the fast and thermal adjoint

(e)

˝(e)

e=1

d˝ N (e) (r)N (e)T (r)G2→1

Nk (r)

(r, r  , ω)

in which, Gg→1 and Gg→2 are the Green’s function components in groups 1 and 2 due to neutron noise source in group g, respectively. The neutron noise source may be considered to be in fast or thermal energy group. The general form of neutron noise source is expressed as Eq. (34):

(e)

˝(e)

(37)



ı(r − r  )

v2



Gg→2 (r, r  , ω)

 or

0

Gg→1 (r, r  ,ω)

 

iω

∂˝(e)V

0

−

d˝ N (e) (r)N (e)T (r)G2→1



⎤

In this study, the neutron noise source is considered to be an absorber of variable strength or the vibrating absorber type. To solve the forward noise equation (Eq. (26)), the Green’s function technique is applied. In this method, the Green’s function (neutron noise due to unit point noise source) is calculated from Eq. (33):

=

(e)

+

The coefficient ˙1 (r, ω) applied in Eq. (28) is defined as Eq. (32): ˙1 (r, ω) = ˙r,1 (r) +

(e)

˝(e)

− ˙a,2 +

(30)

0



(29)

2 (r)

−1 (r)

 (e)

,

⎡

d˝ D2 ∇ N (e) (r)∇ N (e)T (r)G2→2

+˙s,1→2



0

0

(e)

˝(e)

e=1

v2

,

−1 (r)

f (r, ω) = ⎣

1−

keff

˙s,1→2 (r)



2 ˙f,2 (r)

−˙1 (r, ω)

(e)

(r)N

(e)T

(r)

G2→1 2

d˝ N (e) (r)N (e)T (r)G2→2



4. Main specification of the benchmark problems = 0,

(35) Calculations are performed for two typical thermal reactors, namely BIBLIS-2D PWR and IAEA-2D benchmark problems in two

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258 Table 1 The analytical solution of the integrals.

˝(e)

d˝ ∇ N (e) (r)∇ N (e)T (r) =

1 4(e)

⎡

(e) (e)

(e) (e)

bi bi + ci ci

⎢ (e) (e) (e) (e) ⎢ bj bi + cj ci ⎣ (e) (e)

(e) (e)

bk bi + ck ci

⎡

1 12

⎢ ⎢ 1 (e) (e)T (e) d˝ N (r)N (r) = 2 ⎢ 24 ˝(e) ⎣ 1 24

1 24

1 24

1 12

1 24

1 24

1 12

(e) (e)

(e) (e)

bi bk + ci ck

(e) (e)

(e) (e)

bj bk + cj ck

(e) (e)

(e) (e)

bk bk + ck ck

bi bj + ci cj bj bj + cj cj bk bj + ck cj

(e) (e)

(e) (e)

(e) (e)

(e) (e)

(e) (e)

(e) (e)

⎤

255

⎤ ⎥ ⎥ ⎦

⎡1⎤

6 ⎥  ⎢ ⎥ ⎥ (e) (e) ⎢ 1 ⎥ ds N (r) = 2 ⎢ ⎥ ⎥ ⎦ ∂˝(e)V ⎣6⎦ 1 6

Table 2 The material cross section of each assembly for BIBLIS-2D reactor. 2-group constants

M1

M2

M3

M4

M5

M6

M7

M8

D1 (cm) D2 (cm) ˙ f,1 (cm−1 ) ˙ f,2 (cm−1 ) ˙ r,1 (cm−1 ) ˙ r,2 (cm−1 )

1.43600 0.36350 0.00587 0.09607 0.02726 0.07506

1.43600 0.36360 0.00619 0.10358 0.02730 0.07844

1.32000 0.27720 0.00000 0.00000 0.02576 0.07160

1.43890 0.36380 0.00745 0.13236 0.02746 0.09141

1.43810 0.36650 0.00619 0.10358 0.02729 0.08483

1.43850 0.36650 0.00643 0.10911 0.02732 0.08731

1.43890 0.36790 0.00619 0.10358 0.02729 0.088024

1.43930 0.36800 0.00643 0.10911 0.02762 0.09051

Table 3 The material cross section of each assembly for IAEA-2D reactor. 2-group constants

M1

M2

M3

M3

D1 (cm) D2 (cm) ˙ f,1 (cm−1 ) ˙ f,2 (cm−1 ) ˙ r,1 (cm−1 ) ˙ r,2 (cm−1 )

1.500 0.400 0.000 0.135 0.030 0.080

1.500 0.400 0.000 0.135 0.030 0.085

1.500 0.400 0.000 0.135 0.030 0.130

1.500 0.400 0.000 0.000 0.040 0.010

energy group approximation. Fig. 3 displays the BIBLIS-2D reactor core. The boundary conditions of the reactor core comprise of no incoming neutrons for the external boundaries. Table 2 represents the material cross section of each assembly for BIBLIS-2D reactor core. The IAEA-2D reactor is an International Atomic Energy Agency (IAEA) benchmark proposed by Argonne Code Center. The IAEA-2D reactor core and the materials distribution are shown in Fig. 4 and Table 3, respectively. 5. Numerical results and discussion 5.1. Results of static calculations The developed program is employed for BIBLIS-2D reactor core simulation. Table 4 shows the calculated forward and adjoint multiplication factors with their Relative Percent Error (RPE) (defined as Eq. (40)) vs. number of the unstructured triangle elements. Table 4 The calculated forward and adjoint multiplication factors with their RPE for BIBLIS2D. Number of elements in 1/8 reactor core

keff GFEM

79 439 892 2018 4887 10,510

† keff

1.02201 1.02450 1.02484 1.02498 1.02504 1.02506

RPE (%)a −0.3034 −0.0605 −0.0273 −0.0137 −0.0078 −0.0059

GFEM

RPE (%)a

1.02201 1.02450 1.02484 1.02498 1.02504 1.02506

−0.3034 −0.0605 −0.0273 −0.0137 −0.0078 −0.0059

a The reference effective multiplication factor is keff = 1.02512 (Maiani and Montagnini, 2004).

Figs. 5 and 6 display the distribution of the calculated fast and thermal forward fluxes in the reactor core, respectively. Figs. 7 and 8 display the distribution of the calculated fast and thermal adjoint fluxes in the reactor core, respectively. The power distribution (corresponding to 10,510 elements in full core modeling) in 1/8th of the reactor core is compared to the reference data (Maiani and Montagnini, 2004; Nakata and Martin, 1983; Hebert, 1985) in Fig. 9. The reference data are the results of valid methods like the Finite Element Response Matrix, Boundary Element Response Matrix and Hermit methods. To more validate the results, the static core calculations are also performed using DONJON computer code (Varin et al., 2005) for BIBLIS-2D reactor core. The calculated forward/adjoint neutron multiplication factor using DONJON is 1.02508. The calculated forward fast and thermal fluxes averaged in each fuel assembly/reflector of BIBLIS-2D using DONJON computer code and GFEM are compared in Figs. 10 and 11, respectively. Also, Figs. 12 and 13 display the calculated adjoint fast and thermal fluxes averaged in each fuel assembly/reflector, respectively. The results of DONJON computer code are assumed as the reference values to calculate the RPEs in Figs. 10–13. As shown in Figs. 10–13, the results of DONJON computer code and GFEM are in a good agreement. RPE (%) =

calculated value − reference value × 100. reference value

(40)

The same results for IAEA-2D reactor are repeated in Table 5 and Figs. 14–22. In Fig. 18, the calculated power distribution (corresponding to 12,073 elements in full core modeling) in 1/12th Table 5 The calculated forward and adjoint multiplication factors with their RPE for IAEA-2D. Number of elements in 1/12 reactor core

1360 3032 5455 6777 8417 12,073 a

†

keff

keff

GFEM

RPE (%)a

GFEM

RPE (%)a

1.00582 1.00565 1.00559 1.00557 1.00556 1.00555

0.0308 0.0139 0.0080 0.0060 0.0050 0.0040

1.00582 1.00565 1.00559 1.00557 1.00556 1.00555

0.0308 0.0139 0.0080 0.0060 0.0050 0.0040

The reference effective multiplication factor is keff = 1.00551 (Hebert, 2008).

256

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

Fig. 43. Real part of the thermal neutron noise due to the vibrating absorber type of noise source given by the adjoint calculation at frequency 1 Hz for BIBLIS-2D.

of the reactor core is compared to the reference data (Gonzalez et al., 2009; Hebert, 2008). The reference data are the results of TRIVAC (Hebert, 2008) and High Order FEM (Gonzalez et al., 2009). TRIVAC is a collection of flux solvers that might be called from DONJON (version 4), which is a code developed at Ecol Polytechnique de Montreal and released under the GNU Lesser General Public License. The calculated forward/adjoint neutron multiplication factor using DONJON computer code is 1.00557 for IAEA-2D reactor core. As shown in Tables 4 and 5, differences between the calculated forward/adjoint multiplication factor and reference value decreases as the number of elements is increased. The calculated † RPEs for keff , keff in the present work are in the range of other similar published works (Maiani and Montagnini, 2004; Nakata and Martin, 1983; Hebert, 1985, 2008; Gonzalez et al., 2009). The difference between the calculated neutron multiplication factor by GFEM (present work) and DONJON computer code is 2 pcm for both BIBLIS-2D and IAEA-2D reactor cores. Comparison of the calculated and reference power distribution gives the maximum RPEs (%) 0.5850 and 0.2661 for BIBLIS-2D and IAEA-2D, respectively. The obtained RPEs for power distribution are the acceptable values in the reactor core calculations (Maiani and Montagnini, 2004; Gonzalez et al., 2009). Also, Figs. 10–13 and 19–22 show a good agreement between the obtained forward/adjoint flux distributions from the GFEM and DONJON computer code for both BIBLIS-2D and IAEA-2D reactor cores.

thermal neutron noise for three different positions of noise source in BIBLIS-2D reactor core. As mentioned in the other similar works (Demazière, 2004; Malmir et al., 2010), if the neutron noise source is assumed to be a point source located in the position rs and adjoint noise source as a Dirac function located in the position r0 , Eq. (41) could be applied to benchmark the forward noise calculation: †

2 ı2 (rs , r0 , ω) = ı1 (r0 , rs , ω) + ı2 (r0 , rs , ω),

(41)

in which, 2 stands for 2,0 (rs ). Therefore, the calculated adjoint noise can be benchmarked against the forward noise by comparing the calculated amplitude and phase of these two approaches in all the elements inside the reactor core. The results of adjoint noise calculations are compared with forward noise ones in Figs. 35–42. Figs. 35 and 36 show the magnitude (absolute value or modulus) and phase of the right hand side of Eq. (41) for BIBLIS-2D reactor core, respectively. Also, Figs. 37 and 38 display the magnitude and phase of the left hand side of Eq. (41) for BIBLIS-2D reactor core, respectively. The same results for IAEA-2D reactor are repeated in Figs. 39–42. The calculations are performed for 1 Hz frequency by considering the fluctuations of the macroscopic absorption cross section as the noise source.

5.2. Results of noise calculations As shown in Figs. 23 and 24, calculations have been performed using unstructured triangle elements. To show the advantages of utilization of the unstructured triangle elements, the size of the elements in the fuel assembly containing the neutron noise source is considered to be smaller than the elements in other fuel assemblies. This leads to achieve more acceptable accuracy with low cost of the calculations. As expected, at zero frequency and without the neutron noise source the fast and thermal forward noises approach corresponding static fluxes. Figs. 25–28 show the obtained fast and thermal forward noises at zero frequency and without the neutron noise source for BIBLIS-2D and IAEA-2D, respectively. The comparison between static calculations and noise calculations in zero frequency is the first benchmarking process of noise calculations. As an another benchmarking procedure, it is more worthy to displace the noise source across the reactor core diameter and compare the calculated neutron noises. Figs. 29–34 show the calculated

Fig. 44. Phase of the thermal neutron noise due to the vibrating absorber type of noise source given by the adjoint calculation at frequency 1 Hz for BIBLIS-2D.

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

257

Fig. 45. Real part of the thermal neutron noise due to the vibrating absorber type of noise source given by the adjoint calculation at frequency 1 Hz for IAEA-2D.

For modeling the vibrating absorber type of noise source, one can use the obtained forward or adjoint noise due to absorber of variable strength type of noise source (Pázsit, 1992; Demazière, 2004; Malmir et al., 2010). In this case, the induced noise at the position r0 can be determined by Eq. (42):

ı1 (r0 , ω) ı2 (r0 , ω)

 †

=  · ε(ω) · ∇ rs [Gı˙ (rs , r0 , ω)], a

(42)

where the vector ε(ω) is the displacement function describing the vibration around the equilibrium position rs at frequency of ω. ∇ rs is the derivative of the adjoint Green’s function with respect to the first variable. If one would like to estimate the response of a single detector, it is more advantageous to use the adjoint approach instead of the forward one. On the other hand, if the determination of the noise in several positions is needed, only the forward approach is applicable (Pázsit and Glockler, 1983; Demazière, 2004). In such a case, the derivative of the forward Green’s function with respect to the

second variable is required and the forward noise induced by the vibrating absorber type of noise source can be calculated as:

ı1 (¯r , ω) ı2 (¯r , ω)

 =  · ε¯ (ω) · ∇ rs [Gı˙a (¯r , r¯ s , ω)].

(43)

Compared to Eq. (42), calculating the derivative of the forward Green’s function with respect to the second variable is more complicated, since the full space dependence of the adjoint Green’s function with respect to the first variable is known while it is not the same for the forward one. Figs. 43–46 show the real part and phase of forward noise due to vibrating absorber type of noise source. Qualitative comparative analysis shows that calculated neutron noises due to both types of noise source have a good agreement with the results of similar published work (Demazière and Andhill, 2005). As expected, in the case of an absorber of variable strength source, all detectors in reactor core present an in-phase behavior. In the case of a vibrating absorber, some detectors exhibit an in-phase

Fig. 46. Phase of the thermal neutron noise due to the vibrating absorber type of noise source given by the adjoint calculation at frequency 1 Hz for IAEA-2D.

258

S.A. Hosseini, N. Vosoughi / Nuclear Engineering and Design 253 (2012) 238–258

behavior, whereas some others exhibit an out-of-phase behavior (Demazière and Andhill, 2005). In summary, results of the static calculations are highly correlated with the results reported in literature as well as those which are derived from DONJON computer code. Moreover, qualifying the results taken from neutron noise calculations confirms that the proposed approach and the previous studies are highly matched. 6. Conclusion In the present study, 2-D, 2-G static and dynamic simulators based on GFEM have been developed using triangle unstructured elements for both rectangular and hexagonal reactor cores. The advantage of unstructured triangle elements to structured ones is the possibility of implementing smaller elements in the boundary layers, regions with high flux gradient or fuel assembly containing neutron noise source. This leads to obtaining the desired accuracy with low cost of calculations. To solve the neutron noise equation, the Green’s function technique was applied. The neutron noise due to the neutron noise source of the type of either absorber of variable strength or vibrating absorber was calculated using the developed simulator. The results of static simulator were benchmarked against the valid reported data for BIBLIS-2D and IAEA-2D reactor cores and DONJON computer code. To validate the dynamical simulation, the fast and thermal forward noises in case of zero frequency were benchmarked against the corresponding forward fluxes. Both, the forward and the adjoint approaches for noise calculations were then benchmarked against each other, which leads to the fact that the dynamical core simulator works satisfactorily. Overall, the main innovations of the present study might be summarized as follows: • Utilizing the unstructured triangle elements which lead to obtain more acceptable accuracy with low cost of calculations. • Utilizing the commercial Gambit computer code as the mesh generator tool. • Performing the neutron noise calculations based on GFEM. • Development of a reliable tool for static and dynamical calculations which is applicable for both rectangular and hexagonal reactor cores. References Arzhanov, V., Pázsit, I., 2002. Detecting Impacting of BWR Instrument Tubes with Wavelet Techniques, Power Plant Surveillances and Diagnostics—Applied Research with Artificial Intelligence, vol. XIV. Springer, Physica Verlag, pp. 157–174. Ayyoubzadeh, S.M., Vosoughi, N., Ayyoubzadeh, S.M., 2012. On an improved Direct Discrete Method and its application in two dimensional multi-group neutron diffusion equation. Ann. Nucl. Energy 44, 1–7. Bell, G.I., Glasstone, S., 1970. Nuclear Reactor Theory. Van Nostrand Reinhold Company, New York. Cavdar, S., Ozgener, H.A., 2004. A finite element/boundary element hybrid method for 2-D neutron diffusion calculation. Ann. Nucl. Energy 31, 1555–1582.

Demazière, C., 2004. Development of a 2-D, 2-group neutron noise simulator. Ann. Nucl. Energy 31, 647–680. Demazière, C., 2011. Core Sim: a multi-purpose neutronic tool for research and education. Ann. Nucl. Energy 38, 2698–2718. Demazière, C., Andhill, G., 2005. Identification and localization of absorbers of variable strength in nuclear reactors. Ann. Nucl. Energy 32, 812–842. Demazière, C., Pázsit, I., 2004. Development of a method for measuring the MTC by noise analysis and its experimental verification in Ringhals-2. Nucl. Sci. Eng. 148, 1–29. Demazière, C., Pázsit, I., 2008. Power Reactor Noise. Lecture Notes. Department of Nuclear Engineering Chalmers University of Technology, Gothenburg, Sweden. Duderstadt, J.J., Hamilton, L.J., 1991. Nuclear Reactor Analysis. John Wiley & Sons, New York. FDI, 2000. Gambit 1.3.2 Users and Reference Manual. Fluid Dynamics International, Fluid Dynamics Analysis Package Revision 1.3.2. Evanston, IL. Gonzalez, S., Ginestar, D., Verdu, G., 2009. High order finite element method for the lambda modes problem on hexagonal geometry. Ann. Nucl. Energy 36, 1450–1462. Hebert, A., 1985. Application of the hermit method for finite element reactor calculations. Nucl. Sci. Eng. 91, 34–58. Hebert, A., 2008. A Raviart–Thomas–Schneider solution of diffusion equation in hexagonal geometry. Ann. Nucl. Energy 35, 363–376. Hosseini, S.A., Vosoughi, N., 2012. Development of two-dimensional, multigroup neutron diffusion computer code based on GFEM with unstructured triangle elements. Ann. Nucl. Energy, http://dx.doi.org/10.1016/j.anucene.2012.07.032. Kang, C.M., Hansen, K.F., 1971. Finite element method for the neutron diffusion equation. Trans. Am. Nucl. Soc. 14, 199–200. Lamarsh, J.R., 1965. Introduction to Nuclear Reactor Theory. Addison-Wesley Publishing Company. Larsson, V., Demazière, C., 2012. A coupled neutronics/thermal-hydraulics tool for calculating fluctuations in Pressurized Water Reactors. Ann. Nucl. Energy 43, 68–76. Larsson, V., Demazière, C., Pázsit, I., Tran, H.N., 2011. Neutron noise calculations using the Analytical Nodal Method and comparisons with analytical solutions. Ann. Nucl. Energy 38, 808–816. Lewis, E.E., 1981. Finite element approximation to even-parity transport equation. Adv. Nucl. Sci. Technol. 13, 155–325. Maiani, M., Montagnini, B., 2004. A Galekin approach to the boundary elementresponse matrix method for the multigroup neutron diffusion equations. Ann. Nucl. Energy 31, 1447–1475. Malmir, H., Vosoughi, N., Zahedinejad, E., 2010. Development of a 2-D 2-group neutron noise simulator for hexagonal geometries. Ann. Nucl. Energy 37, 1089–1100. Nakata, H., Martin, W.R., 1983. The finite element response matrix method. Nucl. Sci. Eng. 85, 289–305. Pázsit, I., 1992. Dynamic transfer function calculations for core diagnostics. Ann. Nucl. Energy 5, 303–312. Pázsit, I., 2007. Transport Theory and Stochastic Process. Lecture note. Department of Nuclear Engineering, Chalmers University of Technology. Pázsit, I., Glockler, O., 1983. On the neutron noise diagnostics of pressurized water reactor control rod vibrations. I: periodic vibrations. Nucl. Sci. Eng. 85, 167–177. Sunder, R., Baleanu, M., Kieninger, K., Kolbasseff, A., Kuhn, W., Rösler, H., 1991. Experiences and Results with COMOS—An On-line Vibration Analysis and Monitoring System, SMORN-VI. Gatlinburg, Tennessee, USA. Tran, H.N., Demazière, C., 2012. Neutron noise calculation in a hexagonal geometry and comparison with analytical solutions. In: Proceedings of PHYSOR 2012, Knoxville, TN, USA, April 15–20. Varin, E., Hebert, A., Roy, R., Koclas, J., 2005. A user guide for DONJON, version 3.01. Vosoughi., N., Salehi, A., Shahriari, M., 2004. Discrete formulation for twodimensional multigroup neutron diffusion equations. Ann. Nucl. Energy 31, 231–253. Wang, Y., Bangerth, W., Ragusa, J., 2009. Three-dimensional h-adaptivity for the multigroup neutron diffusion equation. Prog. Nucl. Energy 51, 543–555. Williams, M.M.R., 1974. Random Processes in Nuclear Reactors. Pergamon Press, Oxford. Zienkiewicz, O.C., Taylor, R.L., Zhu, J.Z., 2005. The Finite Element Method: Its Basis and Fundamentals, 6th ed. CIMNE, The International Center for Numerical Methods in Engineering, Barcelona, Spain, ISBN 0 75066320 0.