- Email: [email protected]
Development of deterministic code based on the discrete ordinates method for the third-order neutron correlation technique
Available online at www.sciencedirect.com annals of
NUCLEAR ENERGY Annals of Nuclear Energy 35 (2008) 927–936 www.elsevier.com/locate/anucene
Development of deterministic code based on the discrete ordinates method for the third-order neutron correlation technique Tomohiro Endo, Akio Yamamoto *, Yoshihiro Yamane Department of Materials, Physics and Energy Engineering, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan Received 6 October 2006; received in revised form 21 August 2007; accepted 1 September 2007 Available online 10 January 2008
Abstract The third-order neutron correlation technique is one of the subcriticality measurement techniques. The significant feature of this technique is that an absolute value of subcriticality can be evaluated from the measured neutron correlation factors, namely Y1 and X1. In order to extend the applicability of the third-order neutron correlation technique to actual complicated subcritical systems, we previously derived the spatial and neutron energy dependent theoretical formulas of Y1 and X1 by using the first-, second-, and third-order importance functions. These theoretical formulas depend on the locations of both an external neutron source and a neutron detector. The numerical solutions of these formulas are very useful to design the arrangements of the external neutron source and the neutron detector and to investigate the contamination of the spatial effect in the measured neutron correlation factors, Y1 and X1. Therefore, we developed a novel deterministic code that can calculate both Y1 and X1 based on the SN method. In this paper, we focus on the description of the numerical procedure of this code and its verification by comparing it with the results of Monte Carlo simulations. 2008 Elsevier Ltd. All rights reserved.
1. Introduction The third-order neutron correlation technique is one of the subcriticality measurement techniques based on the reactor noise analysis (Furuhashi and Izumi, 1968; Dragt, 1967). In this technique, first, we measure the time series data of neutron counts is measured. Second, we estimate the second- and third-order neutron correlation factors Y and X by analyzing the measured time series data. Finally, we can estimate the absolute subcriticality from Y1 and X1, which are saturation values of Y and X when the counting gate width T tends to infinity. In order to extend the applicability of the third-order neutron correlation technique to actual complicated subcritical systems, we previously derived the theoretical formulas that considered the spatial and neutron energy effects (Endo et al., 2006a,b). Through these derivations,
*
Corresponding author. Tel.: +81 52 789 5121; fax: +81 52 789 3608. E-mail address: [email protected] (A. Yamamoto).
0306-4549/$ - see front matter 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2007.09.014
we clarified that Y1 and X1 can be expressed by the first-, second-, and third-order importance functions related to the neutron detections. Here, the first-order importance function I 1 ð~ r; E; ~ XÞ is a conventional one, and the second- and third-order importance functions, namely, ~ and I 3 ð~ ~ are newly introduced in order to I 2 ð~ r; E; XÞ r; E; XÞ, derive the theoretical formulas of Y1 and X1. These importance functions satisfy the adjoint neutron transport equations and depend on the arrangement of the neutron detector. By integrating these importance functions with respect to the spatial distribution and the energy spectrum of an external neutron source, the theoretical formulas of Y1 and X1 can be expressed. If we obtain the numerical solutions of Y1 and X1 for the actual complicated systems, it will be very useful to design arrangements of the external neutron source and the neutron detector and also to investigate the contamination of the spatial effects in the measured values of Y1 and X1. Although some deterministic and Monte Carlo codes can calculate Y1 (Humbert, 2003; Valentine, 1997), none of them can estimate both Y1 and X1. Therefore, in order to apply the third-order
928
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
neutron correlation technique to the actual subcriticality measurement, a numerical simulation code that can estimate both Y1 and X1 is necessary. In the present study, we focus on the development of the numerical simulation code that can estimate both Y1 and X1 and confirm the correctness of the results. We consider the fact that the importance functions satisfy the adjoint neutron transport equations. Then, we develop a new deterministic code to calculate Y1 and X1 by SN method using the importance function. In Section 2, we derive the theoretical formulas of Y1 and X1 by using the importance functions. In Section 3, we describe the computational procedures of this newly developed code. In Section 4, we show the numerical results of Y1 and X1 obtained using this code. The validity of this code is confirmed by a comparison with the results of Monte Carlo simulations. 2. Derivation of theoretical formulas of Y‘ and X‘ by using the importance functions In this section, we derive the spatial and neutron energydependent theoretical formulas of Y1 and X1 by using the importance functions (Endo et al., 2006b). In the third-order neutron correlation technique, the second- and third-order neutron correlation factors, namely, Y(T) and X(T), respectively, are defined as follows: l ðT Þ 1; Y ðT Þ ¼ 2 lðT Þ l ðT Þ l ðT Þ X ðT Þ ¼ 3 3 2 þ 2; lðT Þ lðT Þ
ð1Þ ð2Þ
where l(T), l2(T), and l3(T) represent the mean, variance, and third-order central moment of neutron counts C(T) detected during the counting gate width T, respectively. Their definitions are as follows: lðT Þ ¼ hCðT Þi;
ð3Þ 2
ð4Þ
3
ð5Þ
l2 ðT Þ ¼ hðCðT Þ lðT ÞÞ i; l3 ðT Þ ¼ hðCðT Þ lðT ÞÞ i;
where Æxæ implies the expected value of x. In a subcritical system, as T increases, Y(T) and X(T) become saturated. By measuring the saturation values of Y(T) and X(T) that are denoted by Y1 and X1, respectively, we can estimate the absolute value of subcriticality. The details of the measurement are reported in our previous paper (Endo et al., 2006a). In our previous study, we have derived the theoretical formulas of Y1 and X1 by using the first-, second-, and third-order importance func~ and I 3 ð~ ~ respections, i.e., I 1 ð~ r; E; ~ XÞ, I 2 ð~ r; E; XÞ r; E; XÞ, tively. Here, these importance functions satisfy the following adjoint transport equations: ~ ¼ Rd ð~ B y I 1 ð~ r; E; XÞ r; EÞ; 2 ~ ¼ Rf ð~ B y I 2 ð~ r; E; XÞ r; EÞhmðm 1ÞiðI 1;f ð~ rÞÞ ;
ð6Þ ð7Þ
~ ¼ Rf ð~ r; E; XÞ r; EÞf3hmðm 1ÞiI 1;f ð~ rÞI 2;f ð~ rÞ By I 3 ð~ 3
þ hmðm 1Þðm 2ÞiðI 1;f ð~ rÞÞ g;
ð8Þ
where B is the adjoint Boltzmann operator defined as Z Z 1 ~ þ Rt ð~ ~!X ~0 Þ By ¼ Xr r; EÞ dE0 dX0 Rs ð~ r; E ! E0 ; X 0 4p Z Z 1 vf ð~ r; E0 Þ ; ð9Þ mRf ð~ r; EÞ dE0 dX0 4p 0 4p Rd ð~ r; EÞ is the macroscopic detection cross-section of the neutron detector; I 1;f ð~ rÞ and I 2;f ð~ rÞ are the fission spectrum-weighted average importance functions: Z 1 Z vf ð~ r; EÞ ~ I 1;f ð~ I 1 ð~ rÞ ¼ dE dX r; E; XÞ; ð10Þ 4p 4p 0 Z 1 Z vf ð~ r; EÞ ~ I 2;f ð~ I 2 ð~ rÞ ¼ dE dX r; E; XÞ; ð11Þ 4p 4p 0 and Æm(m 1)æ and Æm(m 1)(m 2)æ the second- and thirdorder factorial moments of fission neutrons, respectively. The other notations are the commonly used ones. Here, the above-mentioned importance functions satisfy the following boundary conditions at a free surface: ~ ¼ 0; I 1 ð~ rboundary ; E; XÞ ð12Þ ~ ¼ 0; rboundary ; E; XÞ I 2 ð~ ~ ¼ 0; I 3 ð~ rboundary ; E; XÞ
ð13Þ ~ > 0; for ~ nX
ð14Þ
where ~ n is a unit vector in the direction of the outward normal at a boundary ~ rboundary on the surface. These boundary conditions are the same as those of the adjoint neutron flux. By using these importance functions, the theoretical formulas of Y1 and X1 can be expressed as follows: Z 1 Y1 ¼ dV Sð~ rÞfhqiI 2;s ð~ rÞ þ hqðq 1Þi CR V
X1
ðI 1;s ð~ rÞÞ2 g; ð15Þ Z 1 ¼ dV Sð~ rÞfhqiI 3;s ð~ rÞ þ 3hqðq 1ÞiI 1;s ð~ rÞI 2;s ð~ rÞ CR V 3
rÞÞ g; þ hqðq 1Þðq 2ÞiðI 1;s ð~ where CR is the neutron count rate defined as Z 1 X dV Sð~ rÞ hqiI 1;s ð~ rÞ: CR ¼ V
ð16Þ
ð17Þ
q¼0
Here, Sð~ rÞ is the spatial distribution of the source strength, i.e., the expected number of source events occurring per unit time at position ~ r; q is the number of neutrons emitted per external source event; Æqæ, Æq(q 1)æ, and Æq(q 1)(q 2)æ are the factorial moments of the external source neutrons per source event; I 1;s ð~ rÞ, I 2;s ð~ rÞ, and I 3;s ð~ rÞ are the external source spectrum-weighted importance functions, which are given as follows: Z 1 Z v ð~ r; EÞ ~ I 1 ð~ I 1;s ð~ rÞ ¼ dE dX s r; E; XÞ; ð18Þ 4p 4p 0 Z 1 Z v ð~ r; EÞ ~ I 2 ð~ rÞ ¼ dE dX s r; E; XÞ; ð19Þ I 2;s ð~ 4p 4p 0
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
I 3;s ð~ rÞ ¼
Z
Z
1
dE
dX 4p
0
vs ð~ r; EÞ ~ I 3 ð~ r; E; XÞ; 4p
ð20Þ
X1 ¼
1 CR
Z
Z dV
dERf ð~ r; EÞ/ð~ r; EÞf3hmðm 1Þi
CR ¼
Z
Rs;g!g0 ðx; y; zÞ
M X
ð27Þ
wm I 2;g0 ;m ðx; y; zÞ
m
þ mRf ;g ðx; y; zÞI 2;f ðx; y; zÞ þ Rd;g ðx; y; zÞ hmðm 1ÞiðI 1;f ðx; y; zÞÞ2 ; G X
vf ;g0
M X
ð28Þ
wm I 2;g0 ;m ðx; y; zÞ;
ð29Þ
m
þ Rt;g ðx; y; zÞI 3;g;m ðx; y; zÞ ¼ Q3;g ðx; y; zÞ; ð22Þ
G X
Rs;g!g0 ðx; y; zÞ
g0 ¼1
ð23Þ
M X
ð30Þ
wm I 3;g0 ;m ðx; y; zÞ
m
þ mRf ;g ðx; y; zÞI 3;f ðx; y; zÞ þ Rf ;g ðx; y; zÞ
3. Computational procedures We have developed a deterministic code to calculate Y1 and X1 based on the theoretical formulas shown in Section 2. The computational procedures of this code are described as follows. ~ I 2 ð~ ~ and As shown in Eqs. (6)–(8), I 1 ð~ r; E; XÞ, r; E; XÞ, ~ I 3 ð~ r; E; XÞ, satisfy the adjoint neutron transport equations; hence, we can calculate these importance functions by the deterministic transport method, e.g., the SN or PN method (Bell and Glasstone, 1970). In our code, we used the SN method to calculate the importance functions. In this case, Eqs. (6)–(8) are expressed by multi-group adjoint transport ~m ¼ equations with respect to the discrete direction X ðlm ; gm ; nm Þ as follows: o o o lm þ g m þ n m I 1;g;m ðx; y; zÞ ox oy oz þ Rt;g ðx; y; zÞI 1;g;m ðx; y; zÞ ¼ Q1;g ðx; y; zÞ;
g0 ¼1
G X g0 ¼1
0
G X
ð26Þ
o o o lm þ gm þ nm I 3;g;m ðx; y; zÞ ox oy oz
It should be noted that the delayed neutron is not explicitly taken into account in the above derivation. However, we are interested in the saturated value of Y1 and X1 that are obtained for the infinite gate width T. For a large gate width, most of the fission chains will die away. When this physical consideration is taken into account, it is not necessary to independently treat prompt and delayed neutrons, i.e. they can treat all together.
Q1;g ðx; y; zÞ ¼
Q2;g ðx; y; zÞ ¼
Q3;g ðx; y; zÞ ¼ dERd ð~ r; EÞ/ð~ r; EÞ:
wm I 1;g0 ;m ðx; y; zÞ;
m
g0 ¼1
1
dV V
M X
þ Rt;g ðx; y; zÞI 2;g;m ðx; y; zÞ ¼ Q2;g ðx; y; zÞ;
I 2;f ðx; y; zÞ ¼
3
I 1;f ð~ rÞI 2;f ð~ rÞ þ hmðm 1Þðm 2ÞiðI 1;f ð~ rÞÞ g Z þ dV Sð~ rÞf3hqðq 1ÞiI 1;s ð~ rÞI 2;s ð~ rÞ V 3 þhqðq 1Þðq 2ÞiðI 1;s ð~ rÞÞ g ; Z
vf ;g0
o o o lm þ gm þ nm I 2;g;m ðx; y; zÞ ox oy oz
1
0
V
G X g0 ¼1
r; EÞ is the normalized energy spectrum of the where vs ð~ external source. Furthermore, the theoretical formulas of Y1 and X1 can also be expressed by using the total neutron flux /ð~ r; EÞ: Z Z 1 1 dV dERf ð~ r; EÞ/ð~ r; EÞhmðm 1ÞiðI 1;f ð~ rÞÞ2 Y1 ¼ CR V 0 Z 2 þ dV Sð~ rÞhqðq 1ÞiðI 1;s ð~ rÞÞ ð21Þ V
I 1;f ðx; y; zÞ ¼
929
Rs;g!g0 ðx; y; zÞ
M X
ð24Þ
wm I 1;g0 ;m ðx; y; zÞ
m
þ mRf ;g ðx; y; zÞI 1;f ðx; y; zÞ þ Rd;g ðx; y; zÞ;
ð25Þ
f3hmðm 1ÞiI 1;f ðx; y; zÞI 2;f ðx; y; zÞ 3
þ hmðm 1Þðm 2ÞiðI 1;f ðx; y; zÞÞ g; I 3;f ðx; y; zÞ ¼
G X g0 ¼1
vf ;g0
M X
wm I 3;g0 ;m ðx; y; zÞ;
ð31Þ ð32Þ
m
where common notations are used. It should be noted that wm is a quadrature weight of the solid angle and it is normalized as follows: M X
wm ¼ 1:
ð33Þ
m
I1,g,m(x, y, z), I2,g,m(x, y, z) and I3,g,m(x, y, z) can be solved in the same algorithm as the usual angular neutron flux of the fixed source problem. We solve these importance functions by the diamond differencing scheme with a negative flux fix-up (Lewis and Miller, 1984). For example, we explained the algorithm for calculating the first-order importance function I1,g,m(x, y, z). By volume-averaging Eq. (24) over a certain cell, we obtained the balance equation of the importance function I1,g,m: g lm xin yout m I 1;g;m I xout I yin 1;g;m I 1;g;m 1;g;m þ Dx Dy nm zin I 1;g;m I zout þ ð34Þ 1;g;m þ Rt;g I 1;g;m ¼ Q1;g ; Dz where Dx, Dy, and Dz are the widths of the cell in the x-, y-, yin zin and z-directions, respectively; I xin 1;g;m , I 1;g;m , and I 1;g;m are the surface-averaged incoming importance functions across the yout zout yz-, zx-, and xy-planes, respectively; I xout 1;g;m , I 1;g;m , and I 1;g;m are the surface-averaged outgoing importance functions across the yz-, zx-, and xy-planes, respectively; I 1;g;m is
930
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
the cell-averaged importance function; and Q1;g , is the cellaveraged source term. In the diamond differencing scheme, we assume the following relationship between the incoming and outgoing importance functions in each direction: xout I xin 1;g;m þ I 1;g;m ¼ I 1;g;m ; ð35Þ 2 yout I yin 1;g;m þ I 1;g;m ¼ I 1;g;m ; ð36Þ 2 zout I zin 1;g;m þ I 1;g;m ¼ I 1;g;m : ð37Þ 2 By substituting Eqs. (35)–(37) into Eq. (34), we obtain the cell-averaged importance function as follows: gm yout nm zout 2 lDxm I xout þ I þ I þ Qt;g 1;g;m Dy 1;g;m Dz 1;g;m I 1;g;m ¼ ð38Þ 2 lDxm þ gDym þ nDzm þ Rt;g
The term Q1;g contains I 1;g;m , as shown in Eq. (25). Therefore, iterative calculations are necessary to solve Eq. (38) in the actual implementation. This equation shows that we can calculate I 1;g;m from the outgoing importance functions. Further, the incoming importance functions are calculated using the relationship shown in Eqs. (35)–(37), i.e., xout I xin 1;g;m ¼ 2I 1;g;m I 1;g;m ;
ð39Þ
I yin 1;g;m
¼ 2I 1;g;m
I yout 1;g;m ;
ð40Þ
¼ 2I 1;g;m
I zout 1;g;m :
ð41Þ
I zin 1;g;m
According to the above equations, when twice the value of I 1;g;m is less than those of the outgoing importance functions, the incoming importance function becomes negative. If a negative incoming importance function is obtained, the importance function is fixed-up to be zero, and then, I 1;g;m is recalculated by substituting the fixed incoming importance functions into Eq. (34). The incoming importance functions obtained at a particular cell are used as the outgoing importance functions at the adjacent cells. The first outgoing importance function is given by the boundary condition as shown in Eq. (12). Thus, the importance functions can be solved from the outgoing to the incoming ones. The same algorithms are also applied to the case of the higher-order importance functions, namely I2,g,m(x, y, z) and I3,g,m(x, y, z). These importance functions can be solved from a low order to high order by using Eqs. (24)–(32). First, I1,g,m(x, y, z) is obtained by solving the adjoint fixed source problem in which the fixed source is provided by the macroscopic detection cross-section Rd;g ðx; y; zÞ. Second, I2,g,m(x, y, z) is obtained by solving the adjoint fixed source problem in which the fixed source Rf ;g ðx; y; zÞ hmðm 1Þi ðI 1;f ðx; y; zÞÞ2 is obtained by I1,g,m(x, y, z), which is calculated beforehand. Finally, I3,g,m(x, y, z) is obtained by the adjoint fixed source problem with the fixed source problem with the fixed source Rf ;g ðx; y; zÞf3hmðm 1ÞiI 1;f ðx; y; zÞ 3 I 2;f ðx; y; zÞ þ hmðm 1Þðm 2ÞiðI 1;f ðx; y; zÞÞ g.
In order to accelerate the transport sweep of Eqs. (24), (27) and (30), we use the coarse-mesh finite-difference (CMFD) acceleration (Smith and Rhodes, 2000). In our code, the neutron current equation in the diffusion sweep is corrected to reproduce the net neutron current obtained by a transport sweep. For simplicity, we explain this acceleration method for the forward neutron flux in a slab system. The neutron current in the diffusion sweep is treated as follows: J iþ1=2 ¼
2aiþ1=2 Di Diþ1 e iþ1=2 ð/iþ1 þ /i Þ; ð/ /i Þ þ D Di Diþ1 þ Diþ1 Di iþ1 ð42Þ
where i is the mesh number; /i, the total neutron flux of the ith mesh; Ji+1/2, the net neutron current at the (i + 1/2)th plane between the ith and (i + 1)th mesh; Di, the width of the ith mesh; Di, the conventional diffusion coefficient; e iþ1=2 , the correction factor to reproduce the net current D obtained by the transport sweep; and ai+1/2, a stabilization factor (Yamamoto et al., 2004; Yamamoto, 2005). Assum~ n and the net neutron curing that the total neutron flux / i n rent e J iþ1=2 are obtained by the nth transport sweep, the e iþ1=2 is determined as follows: correction factor D 2aiþ1=2 Di Diþ1 ~ n n n ~n þ / ~ ~n e e / /iþ1 /i D iþ1=2 ¼ J iþ1=2 þ iþ1 i Di Diþ1 þ Diþ1 Di ð43Þ It is noted that the stabilization factor ai+1/2 is set on the basis of our experience. If the optical thickness of a cell, which is denoted by s, is smaller than 1, we may set ai+1/2=1. If s > 1, we must set ai+1/2 > s to stably accelerate the transport sweep. After I1,g,m(x, y, z), I2,g,m(x, y, z) and I3,g,m(x, y, z) are calculated, Y1 and X1 can be obtained by using Eqs. (15)– (17). The external source-weighted importance functions in Eqs. (15)–(17) are calculated as follows: I 1;s ðx; y; zÞ ¼
G X
vs;g0
g0 ¼1
I 2;s ðx; y; zÞ ¼
G X
G X g0 ¼1
wm I 1;g0 ;m ðx; y; zÞ;
ð44Þ
wm I 2;g0 ;m ðx; y; zÞ;
ð45Þ
wm I 3;g0 ;m ðx; y; zÞ:
ð46Þ
m
vs;g0
g0 ¼1
I 3;s ðx; y; zÞ ¼
X X m
vs;g0
X m
In the case of Eqs. (15)–(17), if we calculate the importance functions corresponding to the specific neutron detector, then we can estimate Y1 and X1 corresponding to any neutron source by changing only Sð~ rÞ, Æqæ, Æq(q 1)æ, Æq(q 1)(q 2)æ, and vs ð~ r; EÞ, without the recalculation of the importance functions. Alternatively, we can also calculate Y1 and X1 by using the total neutron flux / as shown in Eqs. (21)–(23). An advantage of this method is that Y1 can be calculated without I 2;g0 ;m ðx; y; zÞ and X1 can be calculated without I 3;g0 ;m ðx; y; zÞ. However, it is observed that the recalculation
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
of total neutron flux is required when we calculate Y1 and X1 for different neutron sources. 4. Numerical analysis In this section, we show the numerical results of Y1 and X1 calculated by our code. Then, the validity of this code is confirmed by a comparison with the results of Monte Carlo simulations. Since the numerical simulation code that can estimate both Y1 and X1 is newly developed, verification calculations that confirm the correctness of the derived formulae are performed in a simple geometry. 4.1. Calculation model The system analyzed in this study is a bare homogeneous rectangular parallelepiped system (30 cm · 40 cm · 60 cm) shown in Fig. 1. We assume that this system is a thermal reactor system, and we calculate it with two energy group constants, as listed in Table 1. These group constants are obtained from the Takeda benchmark model 1, a small LWR core (Takeda and Ikeda, 1991). Here, we adjust the group constants of fission and capture under the constraint such that their sum is constant in order to simulate the three subcritical systems whose neutron multiplication factor, keff, values are approximately 0.99, 0.95, and 0.80. We assume that the fissile material is only 235U. The probability distribution of the fission neutrons and the factorial moments, namely Æmæ, Æm(m 1)æ, and Æm(m 1)(m 2)æ, are listed in Table 2 (Gwin et al., 1984). We assume that the external neutron source is 252Cf. This nuclide is a spontaneous fission source, i.e., multiple neutrons are emitted per fission. The probability distribution of the source neutrons and the factorial moments, namely, Æqæ, Æq(q 1)æ and Æq(q 1)(q 2)æ, are also listed
931
in Table 2 (Gwin et al., 1984). The external source spectrum is identical to the fission R spectrum given in Table 1. rÞ, which is the total The total source strength V dV Sð~ spontaneous fission rate in this calculation, is set to be 10,000 fissions/s. Neutron detection is achieved by thermal neutron capture reaction occurring in detector regions, i.e., Rc;g if g ¼ 2 and ðx; y; zÞ 2 detector region; Rd;g ¼ 0 otherwise: ð47Þ Eq. (47) indicates that no physical detector is assumed in the present calculation, i.e., no additional regions or cross-sections are used to represent the detector. The capture reaction is used as an ‘‘ideal’’ detector instead of a physical detector. In order to examine the spatial effects of Y1 and X1 due to the arrangement of the external neutron source and neutron detector, we calculate Y1 and X1 for two types of arrangement: Case I: The spatial distribution of the source strength is uniform over the entire system. The (15 6 x 6 15, 20 6 y 6 20, 25 6 z 6 25). The neutron detector region spans over the entire area of the system (15 6 x 6 15, 20 6 y 6 20, 25 6 z 6 25). In this case, it is considered that the spatial effect is weak. Case II: A point-like source distributed in a volume 1 cm · 1 cm · 1 cm is located at the corner of the system (x = 14.5, y = 19.5, z = 24.5). The detector region spans one-eighth of the area of the system (0 6 x 6 15, 0 6 y 6 20, 0 6 z 6 25). It is considered that the spatial effect in this case is strong. 4.2. SN method
Detector
In this calculation, we divide the geometry into the 120 · 160 · 200 spatial meshes with 0.25-cm mesh widths. We used an S12 quadrature set, i.e., the total solid angle 4p is divided into 168 directions. The convergence criteria are 107 with respect to the total adjoint fission rate over the entire system and the solid angle-averaged importance function at each mesh point. We calculate Y1 and X1 by using both procedures based on Eqs. (15)–(17) and (21)–(23), in order to confirm the equivalence of these procedures.
50 cm
4.3. Monte Carlo method
Source
40 cm
30 cm Fig. 1. Bare homogeneous rectangular parallelepiped system.
The existing Monte Carlo code for the neutron transport calculation does not have the function to calculate both Y1 and X1. In this study, we program a two-group Monte Carlo simulation code that is applicable only to the calculation geometry, as shown in Fig. 1; however, it
932
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
Table 1 Two energy group constants Group
Rt (1/cm)
Ra (1/cm)
Rs;g!1 (1/cm)
Rs;g!2 (1/cm)
Fission spectrum
Neutron velocity (cm/s)
1 2
2.23775E01 1.03864E+00
8.52709E03 1.58196E01
1.92423E01 0.00000E+00
2.28253E02 8.80439E01
1.00000 0.00000
1.00000E+07 2.20000E+05
Group
Subcritical system keff = 0.99
1 2
keff = 0.95
keff = 0.80
Rc (1/cm)
Rf (1/cm)
mRf (1/cm)
Rc (1/cm)
Rf (1/cm)
mRf (1/cm)
Rc (1/cm)
Rf (1/cm)
mRf (1/cm)
4.54992E03 3.12760E02
3.97717E03 1.26920E01
9.69157E03 3.09279E01
4.71061E03 3.64041E02
3.81648E03 1.21792E01
9.29999E03 2.96783E01
5.31321E03 5.56344E02
3.21388E03 1.02562E01
7.83157E03 2.49922E01
Table 2 Neutron emission probability and factorial moments P(v) v 235
U
0
1
2
3
4
5
6
7
8
9
0.0291
0.1660
0.3361
0.3073
0.1333
0.0259
0.0021
0.0002
0.0000
0.0000
0
1
2
3
4
5
6
7
8
9
0.0021
0.0249
0.1225
0.2722
0.3067
0.1876
0.0680
0.0141
0.0018
0.0001
P(q) q 252
Cf
235
U
252
Cf
Æmæ
Æm(m 1)æ
Æm(m 1)(m 2)æ
2.437
4.705
6.891
Æqæ
Æq(q 1)æ
Æq(q 1)(q 2)æ
3.773
12.051
32.026
has a time-dependent tally function for the third-order neutron correlation technique. In this code, we use the Mersenne Twister as the random number generator (Matsumoto and Nishimura, 1998). We perform the analog Monte Carlo simulation with a fixed source. We simulate a measurement of the third-order neutron technique by taking the time dependent tally of the capture reactions that occurred in the detector regions. It should be noted that we skip the tally during time tskip to detect neutrons after
the calculation system reaches a stationary state. We suppose 100,000 successive time gates with the fundamental counting gate width T0. By using the bunching method (Misawa et al., 1990), we calculate Y(kT0) and X(kT0) corresponding to the gate width (kT0) that is proportional to T0 (k = 1, 2, . . .). We repeated these procedures N times until the statistical errors of Y(kT0) and X(kT0) reach adequately small values. Table 3 shows the parameters tskip, T0 and N for three subcritical systems.
Table 3 Parameters of Monte Carlo simulations tskip, T0, and N for three subcritical systems
4.4. Numerical results
Subcritical system
tskip (s)
T0 (s)
N
keff = 0.99 keff = 0.95 keff = 0.80
4.0E02 8.0E03 2.0E03
1.0E04 2.0E05 5.0E06
10,000 50,000 200,000
First, we show the equivalence of the two SN procedures based on Eqs. (15)–(17) and (21)–(23) in Tables 4 and 5. As shown in these tables, the differences between the two procedures are zero in most cases. The maximum relative difference is 5 · 104% for Y1 and 8 · 104%
Table 4 Numerical results of Y1 Solution
Subcritical system keff = 0.99
SN
Eqs. (15)–(17) Eqs. (21)–(23) Monte Carlo a
keff = 0.95
keff = 0.80
Case I
Case II
Case I
Case II
Case I
Case II
1980.08 1980.07 (5 · 104%)a 1990.30 (0.52%)
258.750 258.749 (4 · 104%) 260.470 (0.66%)
87.3862 87.3862 (0%) 87.4970 (0.13%)
13.1105 13.1105 (0%) 13.0220 (0.68%)
8.19280 8.19280 (0%) 8.19150 (0.02%)
1.49227 1.49226 (7 · 104%) 1.48920 (0.21%)
The relative difference with respect to the result of SN method based on Eqs. (15)–(17).
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
933
Table 5 Numerical results of X1 Solution
Subcritical system keff = 0.99
SN
Eqs. (15)–(17) Eqs. (21)–(23) Monte Carlo a
keff = 0.95
keff = 0.80
Case I
Case II
Case I
Case II
Case I
Case II
11,724,270 11,724,270 (0%)a 11,891,000 (1.42%)
191,866 191,866 (0%) 196,570 (2.45%)
22397.9 22397.9 (0%) 22441.0 (0.19%)
436.069 436.070 (2 · 104%) 433.900 (0.50%)
172.630 172.630 (0%) 172.220 (0.24%)
5.27259 5.27255 (8 · 104%) 5.25720 (0.29%)
The relative difference with respect to the result of SN method based on Eqs. (15)–(17).
300 2000 250 200 Y [-]
Y [-]
1500
1000
150 100
500
Monte Carlo
0 1.0E-04
1.0E-03
1.0E-02
1.0E-01
Monte Carlo
50
Yinf, SN
Yinf, SN 0 1.0E-04
1.0E+00
1.0E-03
T [s]
(keff = 0.9, Case I)
1.0E-02 T [s]
1.0E-01
1.0E+00
(keff = 0.99, Case II) 14
100
12 80 10
Y [-]
Y [-]
60
8 6
40
4 Monte Carlo
Monte Carlo
20
2
Yinf, SN 0 1.0E-05
1.0E-04
1.0E-03
1.0E-02
1.0E-01
0 1.0E-05
1.0E+00
Yinf, SN
1.0E-04
1.0E-03
1.0E-02
1.0E-01
1.0E+00
T [s]
T [s]
(keff = 0.95, Case I)
(keff = 0.95, Case II) 1.6
8 1.2
Y [-]
Y [-]
6
4
2
Monte Carlo
0.8
0.4
Monte Carlo
Yinf, SN 0 1.0E-06
1.0E-05
1.0E-04
1.0E-03
1.0E-02
Yinf, SN
1.0E-01
0.0 1.0E-06
1.0E-05
1.0E-04
1.0E-03
1.0E-02
1.0E-01
T [s]
T [s]
(keff = 0.80, Case I)
(keff = 0.80, Case II) Fig. 2. Numerical results of Y.
for X1 (keff = 0.80, Case II). It is considered that these minute differences arise from the diamond differencing scheme with the negative flux fix-up. In the procedure
based on Eqs. (15)–(17), the negative incoming importance functions I1,g,m, I2,g,m and I3,g,m, are fixed-up to be zero. On the other hand, in the procedure based on
934
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
2.0E+05 1.2E+07
X [-]
X [-]
1.5E+05 8.0E+06
1.0E+05
4.0E+06 5.0E+04
Monte Carlo Xinf, SN 0.0E+00 1.0E-04
1.0E-03
1.0E-02
1.0E-01
Monte Carlo Xinf, SN
0.0E+00 1.0E-04
1.0E+00
1.0E-03
1.0E-02
T [s]
1.0E-01
1.0E+00
T [s]
(keff = 0.99, Case II)
(keff = 0.99, Case I) 500
20000
400
15000
300
X [-]
X [-]
25000
10000
200
5000
0 1.0E-05
100
Monte Carlo Xinf, SN
1.0E-04
1.0E-03 1.0E-02
1.0E-01 1.0E+00
0 1.0E-05
Monte Carlo Xinf, SN
1.0E-04
T [s]
1.0E-03
1.0E-02
1.0E-01
1.0E+00
T [s]
(keff = 0.95, Case II)
(keff = 0.95, Case I) 6
200
5
160
4 X [-]
X [-]
120 3
80 2 Monte Carlo
40
0 1.0E-06
1.0E-05
1.0E-04
1.0E-03
1.0E-02
Monte Carlo
1
Xinf, SN
Xinf, SN
1.0E-01
0 1.0E-06
1.0E-05
1.0E-04
1.0E-03
1.0E-02
1.0E-01
T [s]
T [s]
(keff = 0.80, Case II)
(keff = 0.80, Case I) Fig. 3. Numerical results of X.
Eqs. (21)–(23), the negative incoming importance functions, I1,g,m and I2,g,m and the negative outgoing flux are fixed-up to be zero. The fixed-up condition depends on the outgoing importance function or the incoming flux and the cell-averaged importance function or flux at each mesh point in the transport sweep. Therefore, the numerical results obtained by the two procedures may be different. Fortunately, these differences are not serious as shown in our numerical results. Second, we compare the numerical results of Y1 and X1 calculated by the SN code with Y(T) and X(T) obtained by Monte Carlo simulations. Figs. 2 and 3
show the numerical results of Y and X, respectively, for each subcritical system and case. In these figures, the solid lines show Y1 or X1 calculated by the SN code. As shown in Figs. 2 and 3, the results of the Monte Carlo simulations effectively converge on the solid lines as T increases, regardless of the subcriticality and the arrangement of neutron source and detector. In order to compare the absolute values of Y1 and X1, we estimate Y1 and X1 by Monte Carlo simulations using a fitting procedure. The following analytical formulas of Y(T) and X(T) are used to estimate Y1 and X1 by least square fitting:
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
1 eaT Y ðT Þ ¼ Y 1 1 ; aT 1 eaT aT 2 X ðT Þ ¼ X 1 d 1 þ e þ ð1 dÞ aT 3 4eaT e2aT 1 ; 2aT
ð48Þ
ð49Þ
where a and d are additional fitting coefficients. The estimated values of Y1 and X1 are also listed in Tables 4 and 5, respectively. As shown in these tables, the estimated Y1 and X1 of Monte Carlo simulations agree well with the results of the SN code in most cases. The maximum relative difference between the results of the SN code and Monte Carlo simulations is 0.7% for Y1 (keff = 0.95, Case II) and 2.5% for X1 (keff = 0.99, Case II). It is considered that these differences are due to the following factors: (1) Statistical errors of the Monte Carlo simulation. As shown in Fig. 2, the statistical errors of Y(T) are very small. However, as shown in Fig. 3, the statistical errors of X(T) are larger than those of Y(T), particularly when keff = 0.99. These statistical errors degrade the values of Y1 and X1 that are estimated by least square fitting. (2) Systematic errors of the SN code due to the space and solid angle discretization. In the SN code, the numbers of spatial meshes and solid angle directions slightly influence keff of the calculation system. The absolute values of the total flux and importance functions are sensitive to keff, i.e., these values are inversely proportional to the subcriticality (1 keff)/keff. Therefore, the effects of discretization errors become significant as the subcriticality of the system approaches to zero. Although there are some differences between the results of the SN code and Monte Carlo simulations, it is fair to say that these differences are in the range of statistical errors of Monte Carlo simulations, as shown in Figs. 2 and 3. Therefore, we can confirm the validity of our SN code. 5. Conclusion In order to utilize the numerical solutions of Y1 and X1 for the investigation of the spatial and neutron energy effects in the measured neutron correlation factors, Y1 and X1, we have developed a new deterministic code to calculate Y1 and X1 based on the SN method. The calculation procedures of this code are based on the theoretical formulas of Y1 and X1 derived by using the first-, second-, and third-order importance functions. Through the numerical analysis, we confirm the equivalence of the two SN procedures. One is based on Eqs. (15)–(17), which are expressed by only the importance functions, and the other
935
is based on Eqs. (21)–(23), which contain the total neutron flux. Furthermore, we confirm the validity of our SN code by a comparison with the results of Monte Carlo simulations. In general, Monte Carlo simulation requires a long calculation time. In particular, a considerably long calculation time is necessary to estimate Y1 and X1 because we must estimate not only the mean but also the second- and thirdorder factorial moments of neutron counts with low statistical uncertainty. Furthermore, we need to carry out the Analog Monte Carlo simulation in order to estimate Y1 and X1 rigorously so that a larger calculation time is required and this is directly proportional to the source strength and the inversely proportional to the subcriticality. Moreover, a smaller macroscopic detection cross-section causes larger statistical errors in Y1 and X1. On the other hand, our developed code is based on the SN method, which is one of the deterministic transport methods. In the case of the SN method, the calculation time of Y1 and X1 does not strongly depend on the subcriticality, source strength, and detector cross-section, although it depends on the number of spatial meshes, solid angle directions, and neutron energy group. We emphasize that the characteristics of importance functions, namely, Eqs. (15)–(17), decrease the calculation time. Once we calculate the importance functions corresponding to a specific neutron detector, we can calculate Y1 and X1 corresponding to any neutron source without recalculating the importance functions. Therefore, our developed SN code is an efficient tool to investigate the spatial and energy effects of Y1 and X1 for actual complicated subcritical systems. Since the verification calculations are performed in a simple configuration, the applicability of our code to more complicated problems will be confirmed in a future study. References Bell, G.I., Glasstone, S., 1970. Nuclear Reactor Theory. Van Nostrand Reinhold Company, New York. Dragt, J.B., 1967. Threefold correlations and third order moments in reactor noise. Nukleonik 10, 7–13. Endo, T., Yamane, Y., Yamamoto, A., 2006a. Space and energy dependent theoretical formula for the third order neutron correlation technique. Ann. Nucl. Energy 33, 521–537. Endo, T., Yamane, Y., Yamamoto, A., 2006b. Derivation of theoretical formula for the third order neutron correlation technique by using importance function. Ann. Nucl. Energy 33, 857–868. Furuhashi, A., Izumi, A., 1968. Third moment of the number of neutrons detected in short time intervals. J. Nucl. Sci. Technol. 5, 48–59. Gwin, R., Spencer, R.R., Ingle, R.W., 1984. Measurements of the energy dependence of prompt neutron emission from 233U, 235U, 239Pu, and 241 Pu for En = 0.005–10 eV relative to emission from spontaneous fission of 252Cf. Nucl. Sci. Eng. 87, 381–404. Humbert, P., 2003. Neutron noise computation using PANDA deterministic code. In: Proceedings of the International Conference on Supercomputing in Nuclear Applications SNA2003, Paris, France, September 22–24, 2003 (CD-ROM). Lewis, E.E., Miller Jr., W.F., 1984. Computational Methods of Neutron Transport. Wiley–Interscience, New York.
936
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
Matsumoto, M., Nishimura, T., 1998. Mersenne Twister: a 623-dimensionally equidistributed uniform pseudorandom number generator. ACM Trans. Model. Comput. Simul. 8, 3–30. Misawa, T., Shiroya, S., Kanda, K., 1990. Measurement of prompt neutron decay constant and large subcriticality by the Feynmanmethod. Nucl. Sci. Eng. 104, 53–65. Smith, K.S., Rhodes, J.D., 2000. CASMO-4 characteristics method for two-dimensional PWR and BWR core calculations. Trans. Am. Nucl. Soc. 83, 292.
Takeda, T., Ikeda, H., 1991. 3D neutron transport benchmarks. J. Nucl. Sci. Technol. 28, 656–669. Valentine, T.E., 1997. MCNP-DSP Users Manual, ORNL/TM-13334, Oak Ridge National Laboratory. Yamamoto, A., Kitamura, Y., Ushio, T., Sugimura, N., 2004. Convergence improvement of coarse mesh rebalance method for neutron transport calculations. J. Nucl. Sci. Technol. 41, 781–789. Yamamoto, A., 2005. Generalized coarse mesh rebalance method for acceleration of neutron transport calculations. Nucl. Sci. Eng. 151, 274–282.
NUCLEAR ENERGY Annals of Nuclear Energy 35 (2008) 927–936 www.elsevier.com/locate/anucene
Development of deterministic code based on the discrete ordinates method for the third-order neutron correlation technique Tomohiro Endo, Akio Yamamoto *, Yoshihiro Yamane Department of Materials, Physics and Energy Engineering, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan Received 6 October 2006; received in revised form 21 August 2007; accepted 1 September 2007 Available online 10 January 2008
Abstract The third-order neutron correlation technique is one of the subcriticality measurement techniques. The significant feature of this technique is that an absolute value of subcriticality can be evaluated from the measured neutron correlation factors, namely Y1 and X1. In order to extend the applicability of the third-order neutron correlation technique to actual complicated subcritical systems, we previously derived the spatial and neutron energy dependent theoretical formulas of Y1 and X1 by using the first-, second-, and third-order importance functions. These theoretical formulas depend on the locations of both an external neutron source and a neutron detector. The numerical solutions of these formulas are very useful to design the arrangements of the external neutron source and the neutron detector and to investigate the contamination of the spatial effect in the measured neutron correlation factors, Y1 and X1. Therefore, we developed a novel deterministic code that can calculate both Y1 and X1 based on the SN method. In this paper, we focus on the description of the numerical procedure of this code and its verification by comparing it with the results of Monte Carlo simulations. 2008 Elsevier Ltd. All rights reserved.
1. Introduction The third-order neutron correlation technique is one of the subcriticality measurement techniques based on the reactor noise analysis (Furuhashi and Izumi, 1968; Dragt, 1967). In this technique, first, we measure the time series data of neutron counts is measured. Second, we estimate the second- and third-order neutron correlation factors Y and X by analyzing the measured time series data. Finally, we can estimate the absolute subcriticality from Y1 and X1, which are saturation values of Y and X when the counting gate width T tends to infinity. In order to extend the applicability of the third-order neutron correlation technique to actual complicated subcritical systems, we previously derived the theoretical formulas that considered the spatial and neutron energy effects (Endo et al., 2006a,b). Through these derivations,
*
Corresponding author. Tel.: +81 52 789 5121; fax: +81 52 789 3608. E-mail address: [email protected] (A. Yamamoto).
0306-4549/$ - see front matter 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2007.09.014
we clarified that Y1 and X1 can be expressed by the first-, second-, and third-order importance functions related to the neutron detections. Here, the first-order importance function I 1 ð~ r; E; ~ XÞ is a conventional one, and the second- and third-order importance functions, namely, ~ and I 3 ð~ ~ are newly introduced in order to I 2 ð~ r; E; XÞ r; E; XÞ, derive the theoretical formulas of Y1 and X1. These importance functions satisfy the adjoint neutron transport equations and depend on the arrangement of the neutron detector. By integrating these importance functions with respect to the spatial distribution and the energy spectrum of an external neutron source, the theoretical formulas of Y1 and X1 can be expressed. If we obtain the numerical solutions of Y1 and X1 for the actual complicated systems, it will be very useful to design arrangements of the external neutron source and the neutron detector and also to investigate the contamination of the spatial effects in the measured values of Y1 and X1. Although some deterministic and Monte Carlo codes can calculate Y1 (Humbert, 2003; Valentine, 1997), none of them can estimate both Y1 and X1. Therefore, in order to apply the third-order
928
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
neutron correlation technique to the actual subcriticality measurement, a numerical simulation code that can estimate both Y1 and X1 is necessary. In the present study, we focus on the development of the numerical simulation code that can estimate both Y1 and X1 and confirm the correctness of the results. We consider the fact that the importance functions satisfy the adjoint neutron transport equations. Then, we develop a new deterministic code to calculate Y1 and X1 by SN method using the importance function. In Section 2, we derive the theoretical formulas of Y1 and X1 by using the importance functions. In Section 3, we describe the computational procedures of this newly developed code. In Section 4, we show the numerical results of Y1 and X1 obtained using this code. The validity of this code is confirmed by a comparison with the results of Monte Carlo simulations. 2. Derivation of theoretical formulas of Y‘ and X‘ by using the importance functions In this section, we derive the spatial and neutron energydependent theoretical formulas of Y1 and X1 by using the importance functions (Endo et al., 2006b). In the third-order neutron correlation technique, the second- and third-order neutron correlation factors, namely, Y(T) and X(T), respectively, are defined as follows: l ðT Þ 1; Y ðT Þ ¼ 2 lðT Þ l ðT Þ l ðT Þ X ðT Þ ¼ 3 3 2 þ 2; lðT Þ lðT Þ
ð1Þ ð2Þ
where l(T), l2(T), and l3(T) represent the mean, variance, and third-order central moment of neutron counts C(T) detected during the counting gate width T, respectively. Their definitions are as follows: lðT Þ ¼ hCðT Þi;
ð3Þ 2
ð4Þ
3
ð5Þ
l2 ðT Þ ¼ hðCðT Þ lðT ÞÞ i; l3 ðT Þ ¼ hðCðT Þ lðT ÞÞ i;
where Æxæ implies the expected value of x. In a subcritical system, as T increases, Y(T) and X(T) become saturated. By measuring the saturation values of Y(T) and X(T) that are denoted by Y1 and X1, respectively, we can estimate the absolute value of subcriticality. The details of the measurement are reported in our previous paper (Endo et al., 2006a). In our previous study, we have derived the theoretical formulas of Y1 and X1 by using the first-, second-, and third-order importance func~ and I 3 ð~ ~ respections, i.e., I 1 ð~ r; E; ~ XÞ, I 2 ð~ r; E; XÞ r; E; XÞ, tively. Here, these importance functions satisfy the following adjoint transport equations: ~ ¼ Rd ð~ B y I 1 ð~ r; E; XÞ r; EÞ; 2 ~ ¼ Rf ð~ B y I 2 ð~ r; E; XÞ r; EÞhmðm 1ÞiðI 1;f ð~ rÞÞ ;
ð6Þ ð7Þ
~ ¼ Rf ð~ r; E; XÞ r; EÞf3hmðm 1ÞiI 1;f ð~ rÞI 2;f ð~ rÞ By I 3 ð~ 3
þ hmðm 1Þðm 2ÞiðI 1;f ð~ rÞÞ g;
ð8Þ
where B is the adjoint Boltzmann operator defined as Z Z 1 ~ þ Rt ð~ ~!X ~0 Þ By ¼ Xr r; EÞ dE0 dX0 Rs ð~ r; E ! E0 ; X 0 4p Z Z 1 vf ð~ r; E0 Þ ; ð9Þ mRf ð~ r; EÞ dE0 dX0 4p 0 4p Rd ð~ r; EÞ is the macroscopic detection cross-section of the neutron detector; I 1;f ð~ rÞ and I 2;f ð~ rÞ are the fission spectrum-weighted average importance functions: Z 1 Z vf ð~ r; EÞ ~ I 1;f ð~ I 1 ð~ rÞ ¼ dE dX r; E; XÞ; ð10Þ 4p 4p 0 Z 1 Z vf ð~ r; EÞ ~ I 2;f ð~ I 2 ð~ rÞ ¼ dE dX r; E; XÞ; ð11Þ 4p 4p 0 and Æm(m 1)æ and Æm(m 1)(m 2)æ the second- and thirdorder factorial moments of fission neutrons, respectively. The other notations are the commonly used ones. Here, the above-mentioned importance functions satisfy the following boundary conditions at a free surface: ~ ¼ 0; I 1 ð~ rboundary ; E; XÞ ð12Þ ~ ¼ 0; rboundary ; E; XÞ I 2 ð~ ~ ¼ 0; I 3 ð~ rboundary ; E; XÞ
ð13Þ ~ > 0; for ~ nX
ð14Þ
where ~ n is a unit vector in the direction of the outward normal at a boundary ~ rboundary on the surface. These boundary conditions are the same as those of the adjoint neutron flux. By using these importance functions, the theoretical formulas of Y1 and X1 can be expressed as follows: Z 1 Y1 ¼ dV Sð~ rÞfhqiI 2;s ð~ rÞ þ hqðq 1Þi CR V
X1
ðI 1;s ð~ rÞÞ2 g; ð15Þ Z 1 ¼ dV Sð~ rÞfhqiI 3;s ð~ rÞ þ 3hqðq 1ÞiI 1;s ð~ rÞI 2;s ð~ rÞ CR V 3
rÞÞ g; þ hqðq 1Þðq 2ÞiðI 1;s ð~ where CR is the neutron count rate defined as Z 1 X dV Sð~ rÞ hqiI 1;s ð~ rÞ: CR ¼ V
ð16Þ
ð17Þ
q¼0
Here, Sð~ rÞ is the spatial distribution of the source strength, i.e., the expected number of source events occurring per unit time at position ~ r; q is the number of neutrons emitted per external source event; Æqæ, Æq(q 1)æ, and Æq(q 1)(q 2)æ are the factorial moments of the external source neutrons per source event; I 1;s ð~ rÞ, I 2;s ð~ rÞ, and I 3;s ð~ rÞ are the external source spectrum-weighted importance functions, which are given as follows: Z 1 Z v ð~ r; EÞ ~ I 1 ð~ I 1;s ð~ rÞ ¼ dE dX s r; E; XÞ; ð18Þ 4p 4p 0 Z 1 Z v ð~ r; EÞ ~ I 2 ð~ rÞ ¼ dE dX s r; E; XÞ; ð19Þ I 2;s ð~ 4p 4p 0
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
I 3;s ð~ rÞ ¼
Z
Z
1
dE
dX 4p
0
vs ð~ r; EÞ ~ I 3 ð~ r; E; XÞ; 4p
ð20Þ
X1 ¼
1 CR
Z
Z dV
dERf ð~ r; EÞ/ð~ r; EÞf3hmðm 1Þi
CR ¼
Z
Rs;g!g0 ðx; y; zÞ
M X
ð27Þ
wm I 2;g0 ;m ðx; y; zÞ
m
þ mRf ;g ðx; y; zÞI 2;f ðx; y; zÞ þ Rd;g ðx; y; zÞ hmðm 1ÞiðI 1;f ðx; y; zÞÞ2 ; G X
vf ;g0
M X
ð28Þ
wm I 2;g0 ;m ðx; y; zÞ;
ð29Þ
m
þ Rt;g ðx; y; zÞI 3;g;m ðx; y; zÞ ¼ Q3;g ðx; y; zÞ; ð22Þ
G X
Rs;g!g0 ðx; y; zÞ
g0 ¼1
ð23Þ
M X
ð30Þ
wm I 3;g0 ;m ðx; y; zÞ
m
þ mRf ;g ðx; y; zÞI 3;f ðx; y; zÞ þ Rf ;g ðx; y; zÞ
3. Computational procedures We have developed a deterministic code to calculate Y1 and X1 based on the theoretical formulas shown in Section 2. The computational procedures of this code are described as follows. ~ I 2 ð~ ~ and As shown in Eqs. (6)–(8), I 1 ð~ r; E; XÞ, r; E; XÞ, ~ I 3 ð~ r; E; XÞ, satisfy the adjoint neutron transport equations; hence, we can calculate these importance functions by the deterministic transport method, e.g., the SN or PN method (Bell and Glasstone, 1970). In our code, we used the SN method to calculate the importance functions. In this case, Eqs. (6)–(8) are expressed by multi-group adjoint transport ~m ¼ equations with respect to the discrete direction X ðlm ; gm ; nm Þ as follows: o o o lm þ g m þ n m I 1;g;m ðx; y; zÞ ox oy oz þ Rt;g ðx; y; zÞI 1;g;m ðx; y; zÞ ¼ Q1;g ðx; y; zÞ;
g0 ¼1
G X g0 ¼1
0
G X
ð26Þ
o o o lm þ gm þ nm I 3;g;m ðx; y; zÞ ox oy oz
It should be noted that the delayed neutron is not explicitly taken into account in the above derivation. However, we are interested in the saturated value of Y1 and X1 that are obtained for the infinite gate width T. For a large gate width, most of the fission chains will die away. When this physical consideration is taken into account, it is not necessary to independently treat prompt and delayed neutrons, i.e. they can treat all together.
Q1;g ðx; y; zÞ ¼
Q2;g ðx; y; zÞ ¼
Q3;g ðx; y; zÞ ¼ dERd ð~ r; EÞ/ð~ r; EÞ:
wm I 1;g0 ;m ðx; y; zÞ;
m
g0 ¼1
1
dV V
M X
þ Rt;g ðx; y; zÞI 2;g;m ðx; y; zÞ ¼ Q2;g ðx; y; zÞ;
I 2;f ðx; y; zÞ ¼
3
I 1;f ð~ rÞI 2;f ð~ rÞ þ hmðm 1Þðm 2ÞiðI 1;f ð~ rÞÞ g Z þ dV Sð~ rÞf3hqðq 1ÞiI 1;s ð~ rÞI 2;s ð~ rÞ V 3 þhqðq 1Þðq 2ÞiðI 1;s ð~ rÞÞ g ; Z
vf ;g0
o o o lm þ gm þ nm I 2;g;m ðx; y; zÞ ox oy oz
1
0
V
G X g0 ¼1
r; EÞ is the normalized energy spectrum of the where vs ð~ external source. Furthermore, the theoretical formulas of Y1 and X1 can also be expressed by using the total neutron flux /ð~ r; EÞ: Z Z 1 1 dV dERf ð~ r; EÞ/ð~ r; EÞhmðm 1ÞiðI 1;f ð~ rÞÞ2 Y1 ¼ CR V 0 Z 2 þ dV Sð~ rÞhqðq 1ÞiðI 1;s ð~ rÞÞ ð21Þ V
I 1;f ðx; y; zÞ ¼
929
Rs;g!g0 ðx; y; zÞ
M X
ð24Þ
wm I 1;g0 ;m ðx; y; zÞ
m
þ mRf ;g ðx; y; zÞI 1;f ðx; y; zÞ þ Rd;g ðx; y; zÞ;
ð25Þ
f3hmðm 1ÞiI 1;f ðx; y; zÞI 2;f ðx; y; zÞ 3
þ hmðm 1Þðm 2ÞiðI 1;f ðx; y; zÞÞ g; I 3;f ðx; y; zÞ ¼
G X g0 ¼1
vf ;g0
M X
wm I 3;g0 ;m ðx; y; zÞ;
ð31Þ ð32Þ
m
where common notations are used. It should be noted that wm is a quadrature weight of the solid angle and it is normalized as follows: M X
wm ¼ 1:
ð33Þ
m
I1,g,m(x, y, z), I2,g,m(x, y, z) and I3,g,m(x, y, z) can be solved in the same algorithm as the usual angular neutron flux of the fixed source problem. We solve these importance functions by the diamond differencing scheme with a negative flux fix-up (Lewis and Miller, 1984). For example, we explained the algorithm for calculating the first-order importance function I1,g,m(x, y, z). By volume-averaging Eq. (24) over a certain cell, we obtained the balance equation of the importance function I1,g,m: g lm xin yout m I 1;g;m I xout I yin 1;g;m I 1;g;m 1;g;m þ Dx Dy nm zin I 1;g;m I zout þ ð34Þ 1;g;m þ Rt;g I 1;g;m ¼ Q1;g ; Dz where Dx, Dy, and Dz are the widths of the cell in the x-, y-, yin zin and z-directions, respectively; I xin 1;g;m , I 1;g;m , and I 1;g;m are the surface-averaged incoming importance functions across the yout zout yz-, zx-, and xy-planes, respectively; I xout 1;g;m , I 1;g;m , and I 1;g;m are the surface-averaged outgoing importance functions across the yz-, zx-, and xy-planes, respectively; I 1;g;m is
930
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
the cell-averaged importance function; and Q1;g , is the cellaveraged source term. In the diamond differencing scheme, we assume the following relationship between the incoming and outgoing importance functions in each direction: xout I xin 1;g;m þ I 1;g;m ¼ I 1;g;m ; ð35Þ 2 yout I yin 1;g;m þ I 1;g;m ¼ I 1;g;m ; ð36Þ 2 zout I zin 1;g;m þ I 1;g;m ¼ I 1;g;m : ð37Þ 2 By substituting Eqs. (35)–(37) into Eq. (34), we obtain the cell-averaged importance function as follows: gm yout nm zout 2 lDxm I xout þ I þ I þ Qt;g 1;g;m Dy 1;g;m Dz 1;g;m I 1;g;m ¼ ð38Þ 2 lDxm þ gDym þ nDzm þ Rt;g
The term Q1;g contains I 1;g;m , as shown in Eq. (25). Therefore, iterative calculations are necessary to solve Eq. (38) in the actual implementation. This equation shows that we can calculate I 1;g;m from the outgoing importance functions. Further, the incoming importance functions are calculated using the relationship shown in Eqs. (35)–(37), i.e., xout I xin 1;g;m ¼ 2I 1;g;m I 1;g;m ;
ð39Þ
I yin 1;g;m
¼ 2I 1;g;m
I yout 1;g;m ;
ð40Þ
¼ 2I 1;g;m
I zout 1;g;m :
ð41Þ
I zin 1;g;m
According to the above equations, when twice the value of I 1;g;m is less than those of the outgoing importance functions, the incoming importance function becomes negative. If a negative incoming importance function is obtained, the importance function is fixed-up to be zero, and then, I 1;g;m is recalculated by substituting the fixed incoming importance functions into Eq. (34). The incoming importance functions obtained at a particular cell are used as the outgoing importance functions at the adjacent cells. The first outgoing importance function is given by the boundary condition as shown in Eq. (12). Thus, the importance functions can be solved from the outgoing to the incoming ones. The same algorithms are also applied to the case of the higher-order importance functions, namely I2,g,m(x, y, z) and I3,g,m(x, y, z). These importance functions can be solved from a low order to high order by using Eqs. (24)–(32). First, I1,g,m(x, y, z) is obtained by solving the adjoint fixed source problem in which the fixed source is provided by the macroscopic detection cross-section Rd;g ðx; y; zÞ. Second, I2,g,m(x, y, z) is obtained by solving the adjoint fixed source problem in which the fixed source Rf ;g ðx; y; zÞ hmðm 1Þi ðI 1;f ðx; y; zÞÞ2 is obtained by I1,g,m(x, y, z), which is calculated beforehand. Finally, I3,g,m(x, y, z) is obtained by the adjoint fixed source problem with the fixed source problem with the fixed source Rf ;g ðx; y; zÞf3hmðm 1ÞiI 1;f ðx; y; zÞ 3 I 2;f ðx; y; zÞ þ hmðm 1Þðm 2ÞiðI 1;f ðx; y; zÞÞ g.
In order to accelerate the transport sweep of Eqs. (24), (27) and (30), we use the coarse-mesh finite-difference (CMFD) acceleration (Smith and Rhodes, 2000). In our code, the neutron current equation in the diffusion sweep is corrected to reproduce the net neutron current obtained by a transport sweep. For simplicity, we explain this acceleration method for the forward neutron flux in a slab system. The neutron current in the diffusion sweep is treated as follows: J iþ1=2 ¼
2aiþ1=2 Di Diþ1 e iþ1=2 ð/iþ1 þ /i Þ; ð/ /i Þ þ D Di Diþ1 þ Diþ1 Di iþ1 ð42Þ
where i is the mesh number; /i, the total neutron flux of the ith mesh; Ji+1/2, the net neutron current at the (i + 1/2)th plane between the ith and (i + 1)th mesh; Di, the width of the ith mesh; Di, the conventional diffusion coefficient; e iþ1=2 , the correction factor to reproduce the net current D obtained by the transport sweep; and ai+1/2, a stabilization factor (Yamamoto et al., 2004; Yamamoto, 2005). Assum~ n and the net neutron curing that the total neutron flux / i n rent e J iþ1=2 are obtained by the nth transport sweep, the e iþ1=2 is determined as follows: correction factor D 2aiþ1=2 Di Diþ1 ~ n n n ~n þ / ~ ~n e e / /iþ1 /i D iþ1=2 ¼ J iþ1=2 þ iþ1 i Di Diþ1 þ Diþ1 Di ð43Þ It is noted that the stabilization factor ai+1/2 is set on the basis of our experience. If the optical thickness of a cell, which is denoted by s, is smaller than 1, we may set ai+1/2=1. If s > 1, we must set ai+1/2 > s to stably accelerate the transport sweep. After I1,g,m(x, y, z), I2,g,m(x, y, z) and I3,g,m(x, y, z) are calculated, Y1 and X1 can be obtained by using Eqs. (15)– (17). The external source-weighted importance functions in Eqs. (15)–(17) are calculated as follows: I 1;s ðx; y; zÞ ¼
G X
vs;g0
g0 ¼1
I 2;s ðx; y; zÞ ¼
G X
G X g0 ¼1
wm I 1;g0 ;m ðx; y; zÞ;
ð44Þ
wm I 2;g0 ;m ðx; y; zÞ;
ð45Þ
wm I 3;g0 ;m ðx; y; zÞ:
ð46Þ
m
vs;g0
g0 ¼1
I 3;s ðx; y; zÞ ¼
X X m
vs;g0
X m
In the case of Eqs. (15)–(17), if we calculate the importance functions corresponding to the specific neutron detector, then we can estimate Y1 and X1 corresponding to any neutron source by changing only Sð~ rÞ, Æqæ, Æq(q 1)æ, Æq(q 1)(q 2)æ, and vs ð~ r; EÞ, without the recalculation of the importance functions. Alternatively, we can also calculate Y1 and X1 by using the total neutron flux / as shown in Eqs. (21)–(23). An advantage of this method is that Y1 can be calculated without I 2;g0 ;m ðx; y; zÞ and X1 can be calculated without I 3;g0 ;m ðx; y; zÞ. However, it is observed that the recalculation
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
of total neutron flux is required when we calculate Y1 and X1 for different neutron sources. 4. Numerical analysis In this section, we show the numerical results of Y1 and X1 calculated by our code. Then, the validity of this code is confirmed by a comparison with the results of Monte Carlo simulations. Since the numerical simulation code that can estimate both Y1 and X1 is newly developed, verification calculations that confirm the correctness of the derived formulae are performed in a simple geometry. 4.1. Calculation model The system analyzed in this study is a bare homogeneous rectangular parallelepiped system (30 cm · 40 cm · 60 cm) shown in Fig. 1. We assume that this system is a thermal reactor system, and we calculate it with two energy group constants, as listed in Table 1. These group constants are obtained from the Takeda benchmark model 1, a small LWR core (Takeda and Ikeda, 1991). Here, we adjust the group constants of fission and capture under the constraint such that their sum is constant in order to simulate the three subcritical systems whose neutron multiplication factor, keff, values are approximately 0.99, 0.95, and 0.80. We assume that the fissile material is only 235U. The probability distribution of the fission neutrons and the factorial moments, namely Æmæ, Æm(m 1)æ, and Æm(m 1)(m 2)æ, are listed in Table 2 (Gwin et al., 1984). We assume that the external neutron source is 252Cf. This nuclide is a spontaneous fission source, i.e., multiple neutrons are emitted per fission. The probability distribution of the source neutrons and the factorial moments, namely, Æqæ, Æq(q 1)æ and Æq(q 1)(q 2)æ, are also listed
931
in Table 2 (Gwin et al., 1984). The external source spectrum is identical to the fission R spectrum given in Table 1. rÞ, which is the total The total source strength V dV Sð~ spontaneous fission rate in this calculation, is set to be 10,000 fissions/s. Neutron detection is achieved by thermal neutron capture reaction occurring in detector regions, i.e., Rc;g if g ¼ 2 and ðx; y; zÞ 2 detector region; Rd;g ¼ 0 otherwise: ð47Þ Eq. (47) indicates that no physical detector is assumed in the present calculation, i.e., no additional regions or cross-sections are used to represent the detector. The capture reaction is used as an ‘‘ideal’’ detector instead of a physical detector. In order to examine the spatial effects of Y1 and X1 due to the arrangement of the external neutron source and neutron detector, we calculate Y1 and X1 for two types of arrangement: Case I: The spatial distribution of the source strength is uniform over the entire system. The (15 6 x 6 15, 20 6 y 6 20, 25 6 z 6 25). The neutron detector region spans over the entire area of the system (15 6 x 6 15, 20 6 y 6 20, 25 6 z 6 25). In this case, it is considered that the spatial effect is weak. Case II: A point-like source distributed in a volume 1 cm · 1 cm · 1 cm is located at the corner of the system (x = 14.5, y = 19.5, z = 24.5). The detector region spans one-eighth of the area of the system (0 6 x 6 15, 0 6 y 6 20, 0 6 z 6 25). It is considered that the spatial effect in this case is strong. 4.2. SN method
Detector
In this calculation, we divide the geometry into the 120 · 160 · 200 spatial meshes with 0.25-cm mesh widths. We used an S12 quadrature set, i.e., the total solid angle 4p is divided into 168 directions. The convergence criteria are 107 with respect to the total adjoint fission rate over the entire system and the solid angle-averaged importance function at each mesh point. We calculate Y1 and X1 by using both procedures based on Eqs. (15)–(17) and (21)–(23), in order to confirm the equivalence of these procedures.
50 cm
4.3. Monte Carlo method
Source
40 cm
30 cm Fig. 1. Bare homogeneous rectangular parallelepiped system.
The existing Monte Carlo code for the neutron transport calculation does not have the function to calculate both Y1 and X1. In this study, we program a two-group Monte Carlo simulation code that is applicable only to the calculation geometry, as shown in Fig. 1; however, it
932
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
Table 1 Two energy group constants Group
Rt (1/cm)
Ra (1/cm)
Rs;g!1 (1/cm)
Rs;g!2 (1/cm)
Fission spectrum
Neutron velocity (cm/s)
1 2
2.23775E01 1.03864E+00
8.52709E03 1.58196E01
1.92423E01 0.00000E+00
2.28253E02 8.80439E01
1.00000 0.00000
1.00000E+07 2.20000E+05
Group
Subcritical system keff = 0.99
1 2
keff = 0.95
keff = 0.80
Rc (1/cm)
Rf (1/cm)
mRf (1/cm)
Rc (1/cm)
Rf (1/cm)
mRf (1/cm)
Rc (1/cm)
Rf (1/cm)
mRf (1/cm)
4.54992E03 3.12760E02
3.97717E03 1.26920E01
9.69157E03 3.09279E01
4.71061E03 3.64041E02
3.81648E03 1.21792E01
9.29999E03 2.96783E01
5.31321E03 5.56344E02
3.21388E03 1.02562E01
7.83157E03 2.49922E01
Table 2 Neutron emission probability and factorial moments P(v) v 235
U
0
1
2
3
4
5
6
7
8
9
0.0291
0.1660
0.3361
0.3073
0.1333
0.0259
0.0021
0.0002
0.0000
0.0000
0
1
2
3
4
5
6
7
8
9
0.0021
0.0249
0.1225
0.2722
0.3067
0.1876
0.0680
0.0141
0.0018
0.0001
P(q) q 252
Cf
235
U
252
Cf
Æmæ
Æm(m 1)æ
Æm(m 1)(m 2)æ
2.437
4.705
6.891
Æqæ
Æq(q 1)æ
Æq(q 1)(q 2)æ
3.773
12.051
32.026
has a time-dependent tally function for the third-order neutron correlation technique. In this code, we use the Mersenne Twister as the random number generator (Matsumoto and Nishimura, 1998). We perform the analog Monte Carlo simulation with a fixed source. We simulate a measurement of the third-order neutron technique by taking the time dependent tally of the capture reactions that occurred in the detector regions. It should be noted that we skip the tally during time tskip to detect neutrons after
the calculation system reaches a stationary state. We suppose 100,000 successive time gates with the fundamental counting gate width T0. By using the bunching method (Misawa et al., 1990), we calculate Y(kT0) and X(kT0) corresponding to the gate width (kT0) that is proportional to T0 (k = 1, 2, . . .). We repeated these procedures N times until the statistical errors of Y(kT0) and X(kT0) reach adequately small values. Table 3 shows the parameters tskip, T0 and N for three subcritical systems.
Table 3 Parameters of Monte Carlo simulations tskip, T0, and N for three subcritical systems
4.4. Numerical results
Subcritical system
tskip (s)
T0 (s)
N
keff = 0.99 keff = 0.95 keff = 0.80
4.0E02 8.0E03 2.0E03
1.0E04 2.0E05 5.0E06
10,000 50,000 200,000
First, we show the equivalence of the two SN procedures based on Eqs. (15)–(17) and (21)–(23) in Tables 4 and 5. As shown in these tables, the differences between the two procedures are zero in most cases. The maximum relative difference is 5 · 104% for Y1 and 8 · 104%
Table 4 Numerical results of Y1 Solution
Subcritical system keff = 0.99
SN
Eqs. (15)–(17) Eqs. (21)–(23) Monte Carlo a
keff = 0.95
keff = 0.80
Case I
Case II
Case I
Case II
Case I
Case II
1980.08 1980.07 (5 · 104%)a 1990.30 (0.52%)
258.750 258.749 (4 · 104%) 260.470 (0.66%)
87.3862 87.3862 (0%) 87.4970 (0.13%)
13.1105 13.1105 (0%) 13.0220 (0.68%)
8.19280 8.19280 (0%) 8.19150 (0.02%)
1.49227 1.49226 (7 · 104%) 1.48920 (0.21%)
The relative difference with respect to the result of SN method based on Eqs. (15)–(17).
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
933
Table 5 Numerical results of X1 Solution
Subcritical system keff = 0.99
SN
Eqs. (15)–(17) Eqs. (21)–(23) Monte Carlo a
keff = 0.95
keff = 0.80
Case I
Case II
Case I
Case II
Case I
Case II
11,724,270 11,724,270 (0%)a 11,891,000 (1.42%)
191,866 191,866 (0%) 196,570 (2.45%)
22397.9 22397.9 (0%) 22441.0 (0.19%)
436.069 436.070 (2 · 104%) 433.900 (0.50%)
172.630 172.630 (0%) 172.220 (0.24%)
5.27259 5.27255 (8 · 104%) 5.25720 (0.29%)
The relative difference with respect to the result of SN method based on Eqs. (15)–(17).
300 2000 250 200 Y [-]
Y [-]
1500
1000
150 100
500
Monte Carlo
0 1.0E-04
1.0E-03
1.0E-02
1.0E-01
Monte Carlo
50
Yinf, SN
Yinf, SN 0 1.0E-04
1.0E+00
1.0E-03
T [s]
(keff = 0.9, Case I)
1.0E-02 T [s]
1.0E-01
1.0E+00
(keff = 0.99, Case II) 14
100
12 80 10
Y [-]
Y [-]
60
8 6
40
4 Monte Carlo
Monte Carlo
20
2
Yinf, SN 0 1.0E-05
1.0E-04
1.0E-03
1.0E-02
1.0E-01
0 1.0E-05
1.0E+00
Yinf, SN
1.0E-04
1.0E-03
1.0E-02
1.0E-01
1.0E+00
T [s]
T [s]
(keff = 0.95, Case I)
(keff = 0.95, Case II) 1.6
8 1.2
Y [-]
Y [-]
6
4
2
Monte Carlo
0.8
0.4
Monte Carlo
Yinf, SN 0 1.0E-06
1.0E-05
1.0E-04
1.0E-03
1.0E-02
Yinf, SN
1.0E-01
0.0 1.0E-06
1.0E-05
1.0E-04
1.0E-03
1.0E-02
1.0E-01
T [s]
T [s]
(keff = 0.80, Case I)
(keff = 0.80, Case II) Fig. 2. Numerical results of Y.
for X1 (keff = 0.80, Case II). It is considered that these minute differences arise from the diamond differencing scheme with the negative flux fix-up. In the procedure
based on Eqs. (15)–(17), the negative incoming importance functions I1,g,m, I2,g,m and I3,g,m, are fixed-up to be zero. On the other hand, in the procedure based on
934
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
2.0E+05 1.2E+07
X [-]
X [-]
1.5E+05 8.0E+06
1.0E+05
4.0E+06 5.0E+04
Monte Carlo Xinf, SN 0.0E+00 1.0E-04
1.0E-03
1.0E-02
1.0E-01
Monte Carlo Xinf, SN
0.0E+00 1.0E-04
1.0E+00
1.0E-03
1.0E-02
T [s]
1.0E-01
1.0E+00
T [s]
(keff = 0.99, Case II)
(keff = 0.99, Case I) 500
20000
400
15000
300
X [-]
X [-]
25000
10000
200
5000
0 1.0E-05
100
Monte Carlo Xinf, SN
1.0E-04
1.0E-03 1.0E-02
1.0E-01 1.0E+00
0 1.0E-05
Monte Carlo Xinf, SN
1.0E-04
T [s]
1.0E-03
1.0E-02
1.0E-01
1.0E+00
T [s]
(keff = 0.95, Case II)
(keff = 0.95, Case I) 6
200
5
160
4 X [-]
X [-]
120 3
80 2 Monte Carlo
40
0 1.0E-06
1.0E-05
1.0E-04
1.0E-03
1.0E-02
Monte Carlo
1
Xinf, SN
Xinf, SN
1.0E-01
0 1.0E-06
1.0E-05
1.0E-04
1.0E-03
1.0E-02
1.0E-01
T [s]
T [s]
(keff = 0.80, Case II)
(keff = 0.80, Case I) Fig. 3. Numerical results of X.
Eqs. (21)–(23), the negative incoming importance functions, I1,g,m and I2,g,m and the negative outgoing flux are fixed-up to be zero. The fixed-up condition depends on the outgoing importance function or the incoming flux and the cell-averaged importance function or flux at each mesh point in the transport sweep. Therefore, the numerical results obtained by the two procedures may be different. Fortunately, these differences are not serious as shown in our numerical results. Second, we compare the numerical results of Y1 and X1 calculated by the SN code with Y(T) and X(T) obtained by Monte Carlo simulations. Figs. 2 and 3
show the numerical results of Y and X, respectively, for each subcritical system and case. In these figures, the solid lines show Y1 or X1 calculated by the SN code. As shown in Figs. 2 and 3, the results of the Monte Carlo simulations effectively converge on the solid lines as T increases, regardless of the subcriticality and the arrangement of neutron source and detector. In order to compare the absolute values of Y1 and X1, we estimate Y1 and X1 by Monte Carlo simulations using a fitting procedure. The following analytical formulas of Y(T) and X(T) are used to estimate Y1 and X1 by least square fitting:
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
1 eaT Y ðT Þ ¼ Y 1 1 ; aT 1 eaT aT 2 X ðT Þ ¼ X 1 d 1 þ e þ ð1 dÞ aT 3 4eaT e2aT 1 ; 2aT
ð48Þ
ð49Þ
where a and d are additional fitting coefficients. The estimated values of Y1 and X1 are also listed in Tables 4 and 5, respectively. As shown in these tables, the estimated Y1 and X1 of Monte Carlo simulations agree well with the results of the SN code in most cases. The maximum relative difference between the results of the SN code and Monte Carlo simulations is 0.7% for Y1 (keff = 0.95, Case II) and 2.5% for X1 (keff = 0.99, Case II). It is considered that these differences are due to the following factors: (1) Statistical errors of the Monte Carlo simulation. As shown in Fig. 2, the statistical errors of Y(T) are very small. However, as shown in Fig. 3, the statistical errors of X(T) are larger than those of Y(T), particularly when keff = 0.99. These statistical errors degrade the values of Y1 and X1 that are estimated by least square fitting. (2) Systematic errors of the SN code due to the space and solid angle discretization. In the SN code, the numbers of spatial meshes and solid angle directions slightly influence keff of the calculation system. The absolute values of the total flux and importance functions are sensitive to keff, i.e., these values are inversely proportional to the subcriticality (1 keff)/keff. Therefore, the effects of discretization errors become significant as the subcriticality of the system approaches to zero. Although there are some differences between the results of the SN code and Monte Carlo simulations, it is fair to say that these differences are in the range of statistical errors of Monte Carlo simulations, as shown in Figs. 2 and 3. Therefore, we can confirm the validity of our SN code. 5. Conclusion In order to utilize the numerical solutions of Y1 and X1 for the investigation of the spatial and neutron energy effects in the measured neutron correlation factors, Y1 and X1, we have developed a new deterministic code to calculate Y1 and X1 based on the SN method. The calculation procedures of this code are based on the theoretical formulas of Y1 and X1 derived by using the first-, second-, and third-order importance functions. Through the numerical analysis, we confirm the equivalence of the two SN procedures. One is based on Eqs. (15)–(17), which are expressed by only the importance functions, and the other
935
is based on Eqs. (21)–(23), which contain the total neutron flux. Furthermore, we confirm the validity of our SN code by a comparison with the results of Monte Carlo simulations. In general, Monte Carlo simulation requires a long calculation time. In particular, a considerably long calculation time is necessary to estimate Y1 and X1 because we must estimate not only the mean but also the second- and thirdorder factorial moments of neutron counts with low statistical uncertainty. Furthermore, we need to carry out the Analog Monte Carlo simulation in order to estimate Y1 and X1 rigorously so that a larger calculation time is required and this is directly proportional to the source strength and the inversely proportional to the subcriticality. Moreover, a smaller macroscopic detection cross-section causes larger statistical errors in Y1 and X1. On the other hand, our developed code is based on the SN method, which is one of the deterministic transport methods. In the case of the SN method, the calculation time of Y1 and X1 does not strongly depend on the subcriticality, source strength, and detector cross-section, although it depends on the number of spatial meshes, solid angle directions, and neutron energy group. We emphasize that the characteristics of importance functions, namely, Eqs. (15)–(17), decrease the calculation time. Once we calculate the importance functions corresponding to a specific neutron detector, we can calculate Y1 and X1 corresponding to any neutron source without recalculating the importance functions. Therefore, our developed SN code is an efficient tool to investigate the spatial and energy effects of Y1 and X1 for actual complicated subcritical systems. Since the verification calculations are performed in a simple configuration, the applicability of our code to more complicated problems will be confirmed in a future study. References Bell, G.I., Glasstone, S., 1970. Nuclear Reactor Theory. Van Nostrand Reinhold Company, New York. Dragt, J.B., 1967. Threefold correlations and third order moments in reactor noise. Nukleonik 10, 7–13. Endo, T., Yamane, Y., Yamamoto, A., 2006a. Space and energy dependent theoretical formula for the third order neutron correlation technique. Ann. Nucl. Energy 33, 521–537. Endo, T., Yamane, Y., Yamamoto, A., 2006b. Derivation of theoretical formula for the third order neutron correlation technique by using importance function. Ann. Nucl. Energy 33, 857–868. Furuhashi, A., Izumi, A., 1968. Third moment of the number of neutrons detected in short time intervals. J. Nucl. Sci. Technol. 5, 48–59. Gwin, R., Spencer, R.R., Ingle, R.W., 1984. Measurements of the energy dependence of prompt neutron emission from 233U, 235U, 239Pu, and 241 Pu for En = 0.005–10 eV relative to emission from spontaneous fission of 252Cf. Nucl. Sci. Eng. 87, 381–404. Humbert, P., 2003. Neutron noise computation using PANDA deterministic code. In: Proceedings of the International Conference on Supercomputing in Nuclear Applications SNA2003, Paris, France, September 22–24, 2003 (CD-ROM). Lewis, E.E., Miller Jr., W.F., 1984. Computational Methods of Neutron Transport. Wiley–Interscience, New York.
936
T. Endo et al. / Annals of Nuclear Energy 35 (2008) 927–936
Matsumoto, M., Nishimura, T., 1998. Mersenne Twister: a 623-dimensionally equidistributed uniform pseudorandom number generator. ACM Trans. Model. Comput. Simul. 8, 3–30. Misawa, T., Shiroya, S., Kanda, K., 1990. Measurement of prompt neutron decay constant and large subcriticality by the Feynmanmethod. Nucl. Sci. Eng. 104, 53–65. Smith, K.S., Rhodes, J.D., 2000. CASMO-4 characteristics method for two-dimensional PWR and BWR core calculations. Trans. Am. Nucl. Soc. 83, 292.
Takeda, T., Ikeda, H., 1991. 3D neutron transport benchmarks. J. Nucl. Sci. Technol. 28, 656–669. Valentine, T.E., 1997. MCNP-DSP Users Manual, ORNL/TM-13334, Oak Ridge National Laboratory. Yamamoto, A., Kitamura, Y., Ushio, T., Sugimura, N., 2004. Convergence improvement of coarse mesh rebalance method for neutron transport calculations. J. Nucl. Sci. Technol. 41, 781–789. Yamamoto, A., 2005. Generalized coarse mesh rebalance method for acceleration of neutron transport calculations. Nucl. Sci. Eng. 151, 274–282.