- Email: [email protected]
On the Transfer Matrix of the Modified Power Method
Accepted Manuscript On the Transfer Matrix of the Modified Power Method Peng Zhang, Hyunsuk Lee, Deokjung Lee
PII: DOI: Reference:
S0010-4655(17)30329-6 https://doi.org/10.1016/j.cpc.2017.09.022 COMPHY 6336
To appear in:
Computer Physics Communications
Received date : 20 February 2017 Revised date : 12 September 2017 Accepted date : 29 September 2017 Please cite this article as: P. Zhang, H. Lee, D. Lee, On the Transfer Matrix of the Modified Power Method, Computer Physics Communications (2017), https://doi.org/10.1016/j.cpc.2017.09.022 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
On the Transfer Matrix of the Modified Power Method Peng Zhang1,2, Hyunsuk Lee2, Deokjung Lee2,* 1. School of Power and Mechanical Engineering, Wuhan University, Bayilu 299, Wuchang Dist., Wuhan, Hubei, 430072, P. R. China 2. Ulsan National Institute of Science and Technology, UNIST-gil 50, Eonyang-eup, Ulju-gun, Ulsan, 689-798, Republic of Korea *Corresponding author: [email protected]
Abstract The characteristics of the Transfer Matrix (TM) introduced in the modified power method (MPM) have been studied. Because it can be easily mistaken as the Fission Matrix (FM), the differences between the FM and TM are discussed. Theoretically, it can be concluded that the FM is eigenmode dependent unless a very fine mesh is adopted for the FM tally, whereas the TM is based on the coarse mesh and it can give the correct higher eigenmode solutions if the exact weight cancellation can be done. This is confirmed by comparing the analytical solutions of a one-dimensional monoenergetic homogeneous diffusion problem with the solutions of the 2-by-2 FM and TM. It is further confirmed by the numerical tests that the FM tallied with a coarse mesh cannot give correct higher mode solutions, and the FM tallied with i-th mode neutron weights but on a coarse mesh can only give a correct i-th mode solution. The numerical tests also confirm that the TM of various sizes, when different numbers of modes are considered, can give the first several eigenmode solutions correctly and consistently with the same fine mesh based weight cancellation. The impact of the mesh size on the results of the MPM has also been investigated. In practice, the FM only requires the fundamental mode neutron source, but the TM requires simulating the first several eigenmode fission sources explicitly. The FM and the TM can be used to accelerate the convergence of the fundamental mode. The FM uses its fundamental eigenvector to adjust the neutron weights. The TM is used to calculate the combination coefficients which can then be used to update the neutron sources. All the comparisons clearly prove that the TM is different from the FM and that the TM requires further investigation. Keywords: Transfer Matrix, Fission Matrix, Eigenmodes, Modified Power Method
1. Introduction The Modified Power Method (MPM) has been proposed and studied in recent years [1-18]. It is almost identical to the Fission Matrix (FM) method which has been extensively studied in the past few years [19-26]. The two methods have many applications in common, such as higher eigenmode solutions and acceleration of fission source convergence. The power method is used by many Monte Carlo transport codes to get the fundamental eigensolution. The error associated with the power method is a combination of eigenmodes. Besides, recent studies reveal that the inter-cycle correlation and the real statistical errors of the Monte Carlo tallies are also related to the higher eigenmodes [27-29]. Therefore, the higher eigenmodes are of interest to reactor physicists who want to investigate details of the fission source convergence issues, or to develop methods to accelerate the fission source convergence. They are also of interest to those who want to evaluate the real statistical errors more effectively. It should be noted that the higher eigenmodes are of sufficient interests to the neutron transport community. Some Monte Carlo transport codes contain the FM to get the higher modes. The study on MPM shows that it is an alternative way to get the higher modes. 1
Physically, the elements of the FM are the probabilities of neutrons from one point causing fission in other points. This physical meaning is theoretically exact. In practice, the system is continuous, so it is impossible to get the “exact” FM. Space approximation is required, and the FM is obtained on discrete meshes. This FM is exact for the fundamental mode, and very fine meshes are required for the FM to get the higher modes. Therefore, the FM usually requires big memory for 3D big systems like 3D whole core problems. The Transfer Matrix (TM) was introduced for the general solution strategy of the MPM [15]. It is calculated by the neutron weight integrals before and after the power operation. It can keep the balance of fundamental mode, in which case it is identical to the FM. However, it can also keep the balance of the higher modes, which differentiates it from the FM. The calculation of the TM does not require space approximation, whereas the space discretization is used for the representation of the higher eigenmodes. Thus, the memory requirement can be much lower than that of the FM, especially for 3D big systems like 3D whole core problems. The purpose of this paper was to explain the differences between the TM and FM in detail, to establish that the two matrices are essentially different and that the TM requires further investigation. 2. The Fission Matrix Method 2.1 Review of the theory In this part, we will briefly introduce the theory of the Fission Matrix Method (FMM) based on the work by S.E. Carney, F.B. Brown, J. Colin, etc. [20-26]. The k-eigenvalue form of the neutron transport equation is: 1 E S r , M r, E, k 4
(1)
where M is the net loss operator defined by:
r , E , r, E r, E, M r, E, T ' r , E' E, ' r , E ', ' , dE ' d S
which contains the leakage, collision and scattering loss terms. S r
(2)
is the fission neutron
source defined by: ' r , E ' r , E ', ' . S r dE ' d F
(3)
E is the emission fission neutron energy spectrum, which is generally assumed to be isotropic, and independent of space and the neutron energy that causes the fission. The Green’s function for the neutron transport problem is defined by: 0 r, E, r r0 E E 0 , M G r 0 , E0 , 0
2
(4)
where the subscript “0” represents the initial state point. Based on the principle of linear superposition, the neutron flux can be represented by the Green’s function as: 1 E0 S r 0 G r 0 , E , 0 r, E, . d r r, E, 0 dE0 d 0 0 k 4
Multiplying Eq. (5) by F r , E
(5)
and integrating over all the energy and angular space, it
gives:
1 S r d r0S r0 H r0 r , k
where the kernel H r 0 r
(6)
represents the energy-angle averaged Green’s function:
dE d 0 r , E E0 G r 0 , E , 0 r, E, . H r 0 r dEd 0 0 F 4
(7)
Eq. (7) shows that the kernel H is the weighted integration of Green’s function over the initial and final energies and angles. The weight function is the initial fission neutron energy spectrum and the final fission neutron production. To get the Fission Matrix, the system should be divided into discrete spatial regions. Eqs. (6) and (7) integrated over the initial region J and final region I will give: SI
1 N FI , J S J , k J 1
(8)
where N is the total number of discrete spatial regions, and: S r 0 d r 0 , S I S r d r , r 0 RJ rRI S r0 FI ,J d r d r 0 H r0 r . rRI r 0 RJ SJ
S J
(9)
FI ,J is the element of the Fission Matrix, which is the number of fission neutrons produced in the region I due to one fission neutron starting in the region J. Eq. (8) can be rewritten in the matrix form:
S
1 FS, k
(10)
where S is a vector of length N that represents the fission source, F is a matrix of size N-by-N and is called the Fission Matrix, k is the eigenvalue of the Fission Matrix, and the biggest eigenvalue should be positive and equal to the fundamental eigenvalue of the system. 3
2.2 Considerations on higher modes calculation with FMM For the one-speed integral transport equation, all the higher modes exist, with discrete real eigenvalues and real eigenfunctions. For the multigroup transport equation and continuousenergy transport equation, it is conventional practice to assume that higher modes exist, with real eigenvalues and eigenfunctions, even though that has not been proven. The higher mode eigenfunctions can be solved with
1 Sn r d r 0 Sn r 0 H r 0 r , Kn
where K n denotes the n-th eigenvalue of Eq. (6), S n r
(11)
is the n-th eigenfunction, with
n=0 denoting the fundamental mode and n>0 the higher modes. Since the energy-angle averaged Green’s function H is not a symmetric function of its arguments, the forward eigenfunctions do not form an orthogonal basis set. Identical to Eq. (10), the higher eigenmodes can be determined by: Sn
1 Fn S n , kn
(12)
where the subscript n refers to the mode, with n=0 the fundamental mode. According to Eqs. (10) and (12), if Fi = Fj F (i,j=0,1,...,N-1) , all the N eigenmodes can be calculated with only one matrix F of size N-by-N. This requires: S0 r 0 FI ,J d r0 H r0 r rRI r 0 RJ S0, J Sn r 0 H r 0 r , n=0,1,...,N-1, d r d r 0 rRI r 0 RJ S n, J d r
(13)
which means all the fission source eigenmodes have the same distribution in one discrete mesh. If the spatial mesh used for the FM tally is fine enough, so that:
Sn r 0
1 for r V , n 0,1,, N 1, J
Sn, J VJ
(14)
then Eq. (13) is true and the FM can give all the N eigenmodes correctly. However, if the FM is tallied with a coarse mesh, Eq. (13) cannot be satisfied, and the corresponding FM for every eigenmode is different. In other words, if the FM is obtained with the n-th mode fission source as the weighting function, it can only give the n-th mode solution correctly. 3. The Modified Power Method and the Transfer Matrix 4
The MPM was originally proposed by Booth to obtain the first two eigenmodes at the same time and accumulate the convergence of the fundamental mode by subtracting the first mode from it [1-8]. Zhang implemented the MPM in continuous energy Monte Carlo simulation and extended it to obtain even higher modes with a general solution strategy and the TM [1215]. For completeness of the paper, the MPM will be briefly summarized.
3.1 Review of the theory According to Eqs. (6) and (11), the eigenvalue problem to be solved by the power method is:
ASn r d r 0 Sn r 0 H r 0 r K n Sn r , n=0,1,...,
K n , Sn
where A is the power operator and
(15)
are the eigenpairs of the system. It should be
noted that A is a continuous operator, and it is different from the F matrix of Section 2 which is one possible discretization of A. To solve this problem, an initial distribution is needed, that will finally converge to the fundamental mode if the power operation is applied continuously. To solve the first N eigenmodes with the MPM, N independent initial sources are needed: N 1
j cnj S n , j 0,..., N 1, n 0
(16)
N 1
A j cnj K n Sn , j 0,..., N 1, n 0
where j is the j-th source, and the higher modes of order no less than N are neglected. The idea of MPM is that the linear combination of these N independent sources can produce the first N eigenmodes:
0
x00 x0( N 1) N 1 S0 S N 1 , x( N 1)0 x( N 1)( N 1)
(17)
where xij is the linear combination coefficient, and it’s relation to cij is: 1
x00 x0( N 1) c00 c0(N 1) . x( N 1)0 x( N 1)( N 1) c(N 1)0 c(N 1)(N 1)
(17a)
It is noticed that for any sub-region Rm of the system, there is:
Rm
AS n dr K n S n dr , Rm
5
(18)
where n is the mode index and m is the subregion index. Substituting Eq. (17) into Eq. (18) gives: x00 Rm AS0 AS N 1 dr Rm A0 A N 1 dr x ( N 1)0 K0 K 0 S0 K N 1S N 1 dr S0 S N 1 dr Rm Rm 0
Rm
x0( N 1) x( N 1)( N 1) 0 K N 1
(19)
x00 x0( N 1) K 0 0 N 1 dr . x ( N 1)0 x( N 1)( N 1) 0 K N 1
0
If the whole system is divided into N sub-regions, and denoting Wmn A n dr and Rm
Vmn n dr , it will result in: Rm
W10 W1( N 1) x00 x0( N 1) W N 0 WN ( N 1) x(N 1)0 x(N 1)(N 1) V10 V1( N 1) x00 0 x0( N 1) K 0 , VN 0 VN ( N 1) x(N 1)0 x(N 1)(N 1) 0 K N 1
(20)
WX = VXK. The TM is introduced to help solve the linear combination coefficients. Rewriting Eq. (20) gives:
W = VXKX-1 = VX K VX V. -1
(21)
Denoting Smn Sn dr , and recalling Eq. (17), it gives: Rm
V10 V1( N 1) x00 x0( N 1) S10 S1( N 1) , V N 0 VN ( N 1) x(N 1)0 x(N 1)(N 1) S N 0 S N ( N 1) VX S.
Combining Eqs. (21) and (22) gives: 6
(22)
W SKS 1V.
(23)
P SKS 1 ,
(24)
Defining the TM as:
It can be noted that the eigenvalues of the TM are the eigenvalues of the system, and the eigenvectors of the TM are the integrals of the sources on the N sub-regions. The TM can be solved with the weight integrals: P = WV 1 ,
(25)
and then its eigenvectors can be used to solve the linear combination coefficients: X = V -1S, which will then be substituted into Eq. (17) to update the first N eigenmode solutions. Fig. 1 shows the flow chart of the solution strategy of the MPM.
7
(26)
Figure 1. Flow chart of the MPM.
3.2 The accumulation of the TM Due to the inherent statistical noises in Monte Carlo simulation, a TM based on just onecycle tally contains great statistical error, so it should be avoided. Therefore, the accumulation of TM is adopted to reduce the statistical error. In this work, the TM is accumulated from the 2nd cycle, and it is used to update the neutron weights at the end of every cycle.
3.3 Considerations on higher modes calculation with MPM The MPM requires modeling multiple sources explicitly. Instead of simulating more 8
neutrons with different source distributions, the weight vector is introduced for one neutron, with each element of the vector corresponding to one eigenmode. This will not affect the neutron tracking simulation, but different sources can be modelled at the same time when we look at different elements of the weight vector. This makes the application of the MPM in Monte Carlo more attractive. The neutron weights are all positive for the fundamental eigenmode, whereas both positive and negative weights should be modeled for the higher eigenmodes. During the Monte Carlo neutron transport simulation, neutrons with both positive and negative weights will distribute over the entire system and so the weight cancellation procedure is required. In this study, the fine mesh based weight cancellation is adopted. The fine mesh adopted in MPM is different from that in FMM. In MPM, the fine mesh is used to represent the higher eigenmode distributions. In FMM, the fine mesh is used for the FM tally. If the number of fine mesh is M, and N eigenmodes are modeled with the MPM (N << M), the array size for weight cancellation is (N*M), which is usually much smaller than the array size for tallying the FM, (M*M). To compare intuitively, the array size for FM tally and weight cancellation of the MPM with 8 modes are shown in Fig. 2, with the assumption that the number of meshes are the same for different dimensions.
Figure 2. The size of the array for FM tally or weight cancellation. 4. Analytical Solutions of a 1D1G Homogeneous Diffusion Problem To further compare the differences between the FM and the TM, the analytical solutions of a 1D1G homogeneous diffusion problem will be solved, and the 2-by-2 FM and TM will be discussed in this part. For a 1D1G homogeneous diffusion problem, the following equation should be solved:
9
D
d 2 x 1 a x f x , 2 dx k
(27)
where D is the diffusion coefficient, a is the absorption cross section, f is the fission cross section, is the average number of neutrons produced per fission, k is the multiplication factor and the eigenvalue of the problem, and x is the neutron flux. The left-hand side contains the leakage rate and the absorption rate. The right-hand side is the fission neutron source. The Green’s function for this problem can be solved by:
d 2G x0 , x D a G x0 , x x x0 , dx 2
(28)
and the final solution is: G x0 , x
a 1 K x x0 e , K . D 2 KD
(29)
To obtain the eigenvalue k and the corresponding neutron flux x , the power method is commonly adopted. By assuming an initial eigenvalue k(0) and an initial neutron flux (0) x , the neutron flux can be updated by:
(i 1) x G x0 , x
1 f (i ) x0 dx0 , i 0,1,. (i ) k
(30)
Multiplying both sides of Eq. (30) by f gives:
s (i 1) x
1 dx0 s (i ) x0 H x0 , x , i 0,1, , k (i )
(31)
where s x f x and H x0 , x f G x0 , x . The corresponding eigenvalue can be updated by:
s x dx , s x dx ( i 1)
k
( i 1)
(i )
i 0,1,.
(32)
4.1 The exact eigenvalues According to the diffusion theory, we already know that the first two harmonics of the neutron flux are:
10
x , 1 x c1 sin 2a a
0 x c0 cos
x , a x a,
(33)
where c0 and c1 are the normalization factors, and a is the half-thickness of the slab. According to Eqs. (31) and (32), the fundamental eigenvalue should be: k0
dx dx s x H x , x dx dx x G x , x dxs x dx x 0 0
0
0
0
f
0
f
0
f
a
x a
dx
e
2 4K 2a2 2
a
0
0
K x x0
2 Ka 2 2 2 f 8K a 1 e
f
cos x0 dx0 2a a x a cos 2a x dx
a
x0 a
2 KD
0
(34)
.
The 1st mode eigenvalue should be: k1
dx dx s x H x , x dx dx x G x , x dx s x dx x 0 1
0
0
0
1
f 2 KD
a
x a
dx
a
e
f
0
K x x0
Ka 2 2 2 f 4 K a 1 e
a
0
f 1
sin x0 dx0 x0 a a a x a sin a x dx
f 1
4 K 2a2 2
(35)
2
.
4.2 The analytical solution of the FM Suppose the slab is equally divided into 2 sub-regions, that is, x a, 0 as R1 and x 0, a as R2. According to the symmetry of the problem, the FM should be in the
following form: PFM where
11
,
(36)
fission neutrons produced in one region source neutrons in the same region
R1
dx dx0 s x0 H x0 , x R1
R1
fission neutrons produced in one region source neutrons in the other region
R2
dx dx0 s x0 H x0 , x R1
R1
The eigenvectors of the FM should be 1 1 ' and
dxs x
dxs x
,
(37) .
1 1 ' , as:
1 1 , 1 1
1 1 . 1 1
(38)
Accordingly, the eigenvalues of the FM should be: k0FM , k1FM .
(39)
4.2.1 The FM weighted with the fundamental mode fission source
Substituting the fundamental mode fission source into Eq. (37) gives:
R1
dx dx0 f 0 x0 f G x0 , x R1
R1
dx f 0 x
2 2 Ka f 8K a 2 Ka 1 e
2 4K 2a2 2
a
R2
R1
dx f 0 x
Ka Ka f 2 Ka e 1 e
a
2 4K 2a2 2
f
0
x a
dx
cos x0 dx0 2a 0 x a cos 2a x dx (40)
0
x0 a
2 KD
e
K x x0
,
dx dx0 f 0 x0 f G x0 , x R1
f 2 KD
a
x 0
dx
cos x0 dx0 2a 0 x a cos 2a x dx (41)
0
x0 a
e
K x x0
.
The eigenvalues can be calculated by substituting Eqs. (40) and (41) into Eq. (39):
12
k
FM 0
k1FM
2 2 2 2 Ka f 8K a 1 e
a
2 4K 2a2 2
,
2 2 Ka 2 Ka f 8 K a 4 Ka 1 e 1 e
2 4K 2a2 2
a
(42)
2
.
According to Eqs. (34), (35), and (42), it can be concluded that the fundamental mode eigenvalue of the FM is the same as the exact solution, whereas the 1st mode eigenvalue is not exact. As an example, the 1st mode eigenvalues and the error are shown in Fig. 3 with varying K and a. The cross sections are set as: a=0.20 cm-1, f=0.30 cm-1, D=0.33~0.83 cm. The errors of the eigenvalues are calculated as k1FM k1 105 .
Figure 3. Comparison of the 1st mode eigenvalues. 4.2.2 The FM weighted with the first mode fission source
Substituting the 1st mode fission source into Eq. (37) gives:
R1
dx dx0 f 1 x0 f G x0 , x R1
R1
dx f 1 x
2 2 2 Ka f 2 K a 1 e
2 K 2a2 2
a
R2
R1
dx f 1 x
2 2 Ka f 1 e
a 4 K 2a2 2
x a
dx
sin x0 dx0 a 0 x a sin a x dx (43)
0
x0 a
2 KD
e
K x x0
,
dx dx0 f 1 x0 f G x0 , x R1
f
0
.
13
f 2 KD
a
x 0
dx
sin x0 dx0 a 0 x a sin a x dx (44)
0
x0 a
e
K x x0
Substituting Eqs. (43) and (44) into Eq. (39), the eigenvalues can be calculated: k
FM 0
k1FM
2 2 2 Ka Ka f 4 K a 1 e 3 e
4 K 2a2 2
a
2 2 2 Ka f 4 K a 1 e
a
4 K 2a2 2
, (45)
2
.
According to Eqs. (34), (35), and (45), it can be concluded that the 1st mode eigenvalue of the FM now is exact. On the contrary, the fundamental mode eigenvalue is not exact. As an example, the fundamental mode eigenvalues and the error are shown in Fig. 4 with varying K and a. The cross sections are set as: a=0.20 cm-1, f=0.30 cm-1, D=0.33~0.83 cm. The errors of the eigenvalues are calculated as k0FM k0 105 .
Figure 4. Comparison of the fundamental mode eigenvalues.
4.2.3 Discussion
In the previous sections, the 2-by-2 FM is analyzed. However, in practice, the 2-by-2 coarse mesh is never used. The fine mesh is required for the FM to get the accurate results, even if only the first two eigenmodes are to be obtained. In addition, theoretically, if the first mode can be simulated and is used to calculate the FM, the resulting FM will give exactly the first mode solution. 4.3 The analytical solution of the TM 4.3.1 The determination of the TM
The integrals of the eigenfunctions before and after the power operation can be assumed as: C V 0 C0
C1 k0C0 , W C1 k0C0
According to the MPM, the TM can be calculated as: 14
k1C1 . k1C1
(46)
PTM WV 1 k C 0 0 k0C0
1
k1C1 C0 k1C1 C0
k k 0 1 k0 k1
C1 k0C0 C1 k0C0 2 k0 k1 2 . 2 k0 k1 2
The eigenvectors of the TM will be 1 1 ' and k0 k1 k0 k1 k0 k1 k0 k1
2 2 2 2
k0 k1 k0 k1 k0 k1 k0 k1
k1C1 1 C1 k1C1 2C0C1 C0
C1 C0
(47)
1 1 ' , as:
2 1 1 k0 , 2 1 1 2 1 1 k1 , 2 1 1
(48)
and the corresponding eigenvalues of the TM will be k0 and k1 . It means that the TM can give the exact eigensolutions of the system. 4.3.2 Discussion of the first mode solution The first mode eigenvector has both positive and negative values. To simulate it using Monte Carlo method, neutrons with both positive and negative weights should be adopted. Since the neutron positions are randomly sampled, it is almost impossible for two neutrons with different signs of weights to collide with nuclide in the same position. Therefore, the weight cancellation scheme should be developed, otherwise both the neutrons with positive weight and neutrons with negative weights will be distributed in the space according to the fundamental mode and the information of the higher mode will be lost. In this work, the mesh based weight cancellation is adopted [15]. To reduce the discretization error, fine mesh should be adopted. In this way, the MPM is identical to the FMM since they both require the fine meshes. But as shown in Fig. 2, the memory requirements for the fine mesh are different. Besides, the fine mesh information is not needed during the neutron transport simulation, and the memory requirement can be further reduced, as discussed in [15]. However, it should be made clear that the fine mesh based weight cancellation is not essential. There is exact weight cancellation, according to [8], in which case the mesh discretization error can be eliminated in the first eigensolution. In addition, the MPM can be also applied to deterministic calculations, in which case the weight cancellation is not required and the final converged higher mode solutions would be exact. 5. Numerical Tests In this section, a 1D fuel slab of size 400 cm was modeled, with black boundaries on both sides. The fuel composition is listed in Table 1. Eight eigenmodes were considered for this study. For the simulation of MPM, the slab was equally divided into eight sub-regions, and a fine mesh of size 1.0 cm was adopted for the weight cancellation of higher modes. All the 15
simulations were done with 200 inactive cycles, 600 active cycles, and 1,000,000 histories per cycle. Table 1. Composition of the fuel (units: atoms/barn-cm) Nuclide Density H-1 5.9347E-02 N-14 2.1220E-03 O-16 3.7258E-02 U-235 7.6864E-05 U-238 6.8303E-04 5.1 Comparison of the eigenvalues The eigenvalues of the different methods are shown in Fig. 5. For the simulation of FMM, different number of meshes were used, and the eigenvalues of the corresponding FM are shown and compared. “FM_n” indicates that the FM was tallied with n number of meshes. With the increase of the number of meshes, the eigenvalues of the FM increase gradually. This confirms that the FM can produce accurate higher mode solutions if the mesh is fine enough. In addition, the FM can always give the correct fundamental mode eigenvalue regardless of the number of meshes, because the FM was tallied with the fundamental mode neutron weights. For the MPM, 400 meshes are used for the weight cancellation, but the TM is only an 8-by-8 matrix. For the FMM, at least the 128-by-128 FM is needed to give results with similar precision. This confirms the difference between the TM and the FM.
Figure 5. Eigenvalues of different methods.
5.2 Using eigenvectors of the FM to modify the neutron weights 16
The FM can be used to accelerate the convergence of the fundamental fission source. In practice, the fundamental eigenvector of the FM is used to calculate a correction factor to be multiplied with the neutron weights. From Section 5.1, the FM can always give the correct fundamental mode solution, so it can be used to accelerate the fission source convergence with various meshes. During the simulation of MPM, the N eigenvectors of TM were used to calculate the combination factors according to Eq. (26) and modify the neutron weights accordingly at the end of every cycle. This is different from how the FM is used to correct the neutron weights with only one fundamental eigenvector. To clarify the differences, for the following tests, the first N eigenvectors of the FM were used to calculate the combination factors. If the number of meshes is greater than the number of modes, the linear least square solution of the following equation is solved: VX Φ,
(49)
where V R M N contains the weight integrals at the beginning of the cycle, Φ R M N contains the first N eigenvectors of the FM, X R N N contains the combination coefficients, and M N 8 in this case. The Shannon Entropy results are shown in Fig. 6, where “FM_n1_n2” represents the results of the FM based on n1 meshes while n2 modes are considered. Since the FM can always provide a correct fundamental mode solution, the Shannon Entropy increases quickly at the beginning for all the cases. However, comparing the results of “FM_8_8” and “TM_8_8”, the result with the FM converges very slowly after the first several cycles. This reveals that the coarse mesh based FM cannot provide exact higher mode eigenvectors. With an increase of the number of meshes, the FM can produce the first 8 modes accurately, and so the convergence of the Shannon Entropy becomes faster. In addition to the reduced cycle for the source convergence, there may be concerns about the actual simulation time of the MPM since multiple eigenmodes need to be explicitly simulated. To make it clear, the MPM simulations with different number of modes have been conducted, and the simulation times have been compared with the Original Monte Carlo simulation with the same number of histories and processes. The results are shown in Fig. 7. The simulation time of the MPM is not several times the Original Monte Carlo simulation time, but just ~5% higher per one more mode to be simulated. This is because all the eigenmodes are simulated at the same time by the introduction of the weight vector to one particle, and the time for adjusting the weights are much less than the time of particle tracking information sampling. To better illustrate the differences between the TM and FM, another test has been done, and the results are shown in Fig. 8. “FM_8_1” means the FM is tallied based on 8 meshes, and the fundamental mode of the FM is used to modify the neutron source. This is identical to how the FM is used to accelerate the source convergence, in which a correction factor is calculated for each mesh according to the fundamental eigenvector of the FM. The “FM_8_1” result shows efficient source convergence acceleration at the beginning, but the final steady value deviates from the others. This is because of the eigenvector of the coarse mesh based FM. To confirm that, the “FM_100_1” simulation has been done, in which case 100 uniform meshes are used to tally the FM, and the final steady value is then the same as that of “TM_8_8” and “Original”.
17
Figure 6. Shannon Entropy results with different methods.
Figure 7. The time increase with the number of modes for the MPM.
18
Figure 8. Shannon Entropy results with different methods. Table 2. The eigenvalues by different methods for 1D fuel slab problem. FM tallied with weights of different modes Tally TM w0 w1 w2 w3 w4 w5 w6 k0 k1 k2 k3 k4 k5 k6 k7
1.29268 (1.6) 1.28491 (7.1) 1.27229 (6.8) 1.25505 (7.1) 1.23344 (7.1) 1.20804 (7.1) 1.17915 (7.2) 1.14733 (7.4)
w7
1.29251
1.29268
1.29051
1.28938
1.28690
1.28836
1.28695
1.28747
1.28229
1.28475
1.26374
1.28492
1.28037
1.27667
1.27528
1.27479
1.27458
1.26814
1.27219
1.22604
1.20610
1.27238
1.26528
1.25918
1.25710
1.25644
1.24646
1.25463
1.16553
1.19339
1.14022
1.25501
1.24312
1.23732
1.23388
1.22133
1.23312
1.10715
1.08646
1.12154
1.05271
1.23350
1.22088
1.21152
1.19673
1.20723
1.05009
1.07392
1.10798
1.03867
0.99669
1.20829
1.19254
1.17447
1.17844
1.01183
0.99908
0.97886
1.02747
0.98234
0.94490
1.17943
1.15991
1.14694
0.98650
0.99390
0.96899
1.02104
0.96525
0.92705
0.89508
1.14814
5.3 The FM tallied with weights of different modes Since the neutron weights corresponding to N eigenmodes are tracked explicitly during the simulation of MPM, it is possible to tally the FM with weights of different modes. As illustrated in Section 2, the FM tallied with weights of mode n (FM_n) will give a correct n-th eigenmode solution. The following tests will confirm this. Table 2 shows the eigenvalues of the cycle tallies, the TM, and the FM with weights of different modes. The FM was tallied and accumulated from the beginning of active cycles 19
during the simulation of MPM when all modes are converged. From Table 2 it can be confirmed that the FM tallied with different mode weights can give the correct eigenvalue corresponding to that mode. 5.4 The results of the TM with a different number of modes Several simulations of the MPM with different numbers of modes have been done. The Shannon Entropy results are shown in Fig. 9. All the Shannon Entropy results nearly converge to the same value finally. Besides, when a larger number of modes is considered, generally the convergence of the Shannon Entropy is faster.
Figure 9. Shannon Entropy results of the MPM with different numbers of modes. The Eigenvalue results are illustrated in Fig. 10. When a greater number of modes is considered, more eigenvalues can be produced. The first several eigenvalues are always the same.
20
Figure 10. Eigenvalues of the TM with different numbers of modes. 5.5 The results of the TM with different fine meshes for weight cancellation As discussed previously, the fine mesh is essentially not required for the MPM. In this study, to process weight cancellation with high efficiency, a fine mesh average based weight cancellation procedure is adopted. Therefore, it is necessary to investigate how this fine mesh will affect the results. Several simulations of the MPM with 8 modes but with different fine meshes have been done. The eigenvalue results are shown in Fig. 11. The weight cancellation process with different numbers of meshes will affect the results of the MPM, which is expected. However, the pattern with decreasing the number of meshes is different from that of the FM shown in Fig. 5. There are very small differences between the results of the TM with more than 100 meshes for the weight cancellation, which confirms the reasonability and the accuracy of the 400-mesh based weight cancellation adopted in previous simulations of the MPM.
21
Figure 11. Eigenvalues of the TM with different meshes for the weight cancellation. 6. Conclusions The characteristics of the TM introduced in the MPM have been studied. Due to its similarity to the FM, the TM is usually mistakenly thought of as the same as the FM. In this paper, the FMM and the MPM have been reviewed, and the differences between the TM and FM have been discussed. From the theoretical view, the FM can produce the fundamental mode solution correctly regardless of the mesh size, but for higher mode solutions the fine mesh is required for the tally of the FM. On the other hand, the TM is calculated by solving the linear balance equation system of the N eigenmodes on N predefined sub-regions. The calculation of the TM is based on the coarse mesh. The fine mesh is adopted in this study for the weight cancellation of the higher modes, but it is not essential. The analytical solutions of a one-dimensional monoenergetic homogeneous diffusion problem with the 2-by-2 FM and TM have been calculated. The results confirm that the FM with such coarse mesh can only give one exact eigenmode solution, i.e., the coarse mesh FM is eigenmode dependent. The TM can always give the correct fundamental mode solution, but the higher mode solution depends on the weight cancellation precision. The numerical tests with a one-dimensional fuel slab problem have also been conducted. It is confirmed again that the FM can give more accurate results for the higher modes if a finer mesh is used. This is also confirmed by using the eigenvectors of the FM to update the eigenmodes. In addition, it is also proved that the FM is eigenmode dependent by the numerical results of the FM tallied with weights of different modes. On the other hand, by comparing the TM with different numbers of modes, it is demonstrated that under a certain condition of the weight cancellation, the TM can give all the eigenmode solutions correctly regardless of the number of eigenmodes and coarse meshes. Generally, it can make the convergence of the fundamental mode faster when more modes are considered. The impact of 22
the mesh average based weight cancellation on the results of the MPM has also been investigated, and it is found to be a little different from that of the FMM. Therefore, the TM is essentially different from the FM. Both can give the correct fundamental mode solution regardless of the mesh size. Both can be used to accelerate the convergence of the fission source. Whereas for the FM only the fundamental eigenvector is used, for the TM the first N eigenvectors are used. To get the accurate higher mode solutions, the FM requires very fine mesh and large memory for the tally of the FM, in which case the TM requires an accurate weight cancellation procedure and the memory requirement can be lower. Further investigation will be conducted on the exact weight cancellation and the applications of the MPM. Acknowledgements This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (No. 2017M2B2A9A02049916). References 1. T.E. Booth, Nucl. Sci. Eng. 143 (2003) 291-300. 2. T.E. Booth, Nucl. Sci. Eng. 154 (2006) 48-62. 3. J.E. Gubernatis, T.E. Booth, J. Comp. Phys. 227 (2008) 8508-8522. 4. T.E. Booth, J.E. Gubernatis, Monte Carlo Determination of Multiple Extremal Eigenpairs, LA-UR-07-7672, Los Alamos National Laboratory, 2008. 5. T.E. Booth, J.E. Gubernatis, Multiple Extremal Eigenpairs of Very Large Matrices by Monte Carlo Simulation, LA-UR-08-0043, Los Alamos National Laboratory, 2008. 6. T.E. Booth, J.E. Gubernatis, Improved Criticality Convergence via a Modified Monte Carlo Power Iteration Method, Proc. M&C 2009, Saratoga Springs, New York, May 3-7, 2009. 7. T.E. Booth, J.E. Gubernatis, Phys. Rev. E 80 (2009) 046704. 8. T.E. Booth, J.E. Gubernatis, Nucl. Sci. Eng. 165 (2010) 283-291. 9. T. Yamamoto, Ann. Nucl. Energy 36 (2009) 7-14. 10. B. Shi, B. Petrovic, Ann. Nucl. Energy 38 (2011) 781-787. 11. B. Shi, B. Petrovic, Nucl. Sci. Eng. 172 (2012) 138-150. 12. P. Zhang, H. Lee, D. Lee, Extension of Tom Booth’s Modified Power Method for Higher Eigen Modes, Trans. KNS Spring Meeting, Jeju, Korea, May 7-8, 2015. 13. P. Zhang, H. Lee, D. Lee, The Implementation of Modified Power Iteration Method in Continuous Energy Monte Carlo Simulation, Proc. ANS 2015 Annual Meeting, San Antonio, TX, June 7-11, 2015. 14. P. Zhang, H. Lee, D. Lee, Stabilization Technique of Modified Power Iteration Method for Monte Carlo Simulation of Neutron Transport Eigenvalue Problem, Proc. M&C 2015, Nashville, TN, April 19-23, 2015. 15. P. Zhang, H. Lee, D. Lee, J. Comp. Phys. 305 (2016) 387-402. 16. P. Zhang, H. Lee, D. Lee, J. Comp. Phys. 320 (2016) 17-32. 17. P. Zhang, H. Lee, D. Lee, Nucl. Eng. Tech. 49 (2017) 17-28. 18. P. Zhang, H. Lee, D. Lee, J. Comp. Phys. 344 (2017) 440-450. 19. K Tuttelberg, J. Dufek, Ann. Nucl. Energy 72 (2014) 151-155. 20. S.E. Carney, F.B. Brown, B.C. Kiedrowski, W.R. Martin, Ann. Nucl. Energy 73 (2014) 423-431. 21. F.B. Brown, S.E. Carney, B.C. Kiedrowski, W.R. Martin, Fission Matrix Capability for MCNP Monte Carlo, Proc. SNA + MC 2013, Paris, France, October 27-31, LA-UR-13-26962, Los Alamos National Laboratory, 2013. 22. J. Colin, Veit Max D., Eigenfunction Decomposition of Reactor Perturbations and Transitions Using MCNP Monte Carlo, report LA-UR-13-26449, talk LA-UR-13-26464, Los Alamos National Laboratory, 2013. 23. S.E. Carney, F.B. Brown, B.C. Kiedrowski, W.R. Martin, Higher-Mode Applications of Fission Matrix Capability for MCNP, report LA-UR-13-27078, talk LA-UR-13-26615, Los Alamos National Laboratory, 2013. 24. F.B. Brown, S.E. Carney, B.C. Kiedrowski, W.R. Martin, Fission Matrix Capability for MCNP, Part I – Theory, Proc. M&C 2013, Sun Valley, Idaho, May 5-9, 2013. 25. F.B. Brown, S.E. Carney, B.C. Kiedrowski, W.R. Martin, Fission Matrix Capability for MCNP, Part II – Applications, Proc. M&C 2013, Sun Valley, Idaho, May 5-9, 2013. 23
26. S.E. Carney, F.B. Brown, B.C. Kiedrowski, W.R. Martin, Fission Matrix Capability for MCNP Monte Carlo, LA-UR-12-24533, Los Alamos National Laboratory 2012. 27. J. Miao, B. Forget, K. Smith, Ann. Nucl. Energy 92 (2016) 81-95. 28. M. Nowak, J. Miao, E. Dumonteil, B. Forget, A. Onillon, K. S. Smith, A. Zoia, Ann. Nucl. Energy 94 (2016) 856-868. 29. T. M. Sutton, Nucl. Sci. Eng. 185 (2017) 174-183.
24
PII: DOI: Reference:
S0010-4655(17)30329-6 https://doi.org/10.1016/j.cpc.2017.09.022 COMPHY 6336
To appear in:
Computer Physics Communications
Received date : 20 February 2017 Revised date : 12 September 2017 Accepted date : 29 September 2017 Please cite this article as: P. Zhang, H. Lee, D. Lee, On the Transfer Matrix of the Modified Power Method, Computer Physics Communications (2017), https://doi.org/10.1016/j.cpc.2017.09.022 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
On the Transfer Matrix of the Modified Power Method Peng Zhang1,2, Hyunsuk Lee2, Deokjung Lee2,* 1. School of Power and Mechanical Engineering, Wuhan University, Bayilu 299, Wuchang Dist., Wuhan, Hubei, 430072, P. R. China 2. Ulsan National Institute of Science and Technology, UNIST-gil 50, Eonyang-eup, Ulju-gun, Ulsan, 689-798, Republic of Korea *Corresponding author: [email protected]
Abstract The characteristics of the Transfer Matrix (TM) introduced in the modified power method (MPM) have been studied. Because it can be easily mistaken as the Fission Matrix (FM), the differences between the FM and TM are discussed. Theoretically, it can be concluded that the FM is eigenmode dependent unless a very fine mesh is adopted for the FM tally, whereas the TM is based on the coarse mesh and it can give the correct higher eigenmode solutions if the exact weight cancellation can be done. This is confirmed by comparing the analytical solutions of a one-dimensional monoenergetic homogeneous diffusion problem with the solutions of the 2-by-2 FM and TM. It is further confirmed by the numerical tests that the FM tallied with a coarse mesh cannot give correct higher mode solutions, and the FM tallied with i-th mode neutron weights but on a coarse mesh can only give a correct i-th mode solution. The numerical tests also confirm that the TM of various sizes, when different numbers of modes are considered, can give the first several eigenmode solutions correctly and consistently with the same fine mesh based weight cancellation. The impact of the mesh size on the results of the MPM has also been investigated. In practice, the FM only requires the fundamental mode neutron source, but the TM requires simulating the first several eigenmode fission sources explicitly. The FM and the TM can be used to accelerate the convergence of the fundamental mode. The FM uses its fundamental eigenvector to adjust the neutron weights. The TM is used to calculate the combination coefficients which can then be used to update the neutron sources. All the comparisons clearly prove that the TM is different from the FM and that the TM requires further investigation. Keywords: Transfer Matrix, Fission Matrix, Eigenmodes, Modified Power Method
1. Introduction The Modified Power Method (MPM) has been proposed and studied in recent years [1-18]. It is almost identical to the Fission Matrix (FM) method which has been extensively studied in the past few years [19-26]. The two methods have many applications in common, such as higher eigenmode solutions and acceleration of fission source convergence. The power method is used by many Monte Carlo transport codes to get the fundamental eigensolution. The error associated with the power method is a combination of eigenmodes. Besides, recent studies reveal that the inter-cycle correlation and the real statistical errors of the Monte Carlo tallies are also related to the higher eigenmodes [27-29]. Therefore, the higher eigenmodes are of interest to reactor physicists who want to investigate details of the fission source convergence issues, or to develop methods to accelerate the fission source convergence. They are also of interest to those who want to evaluate the real statistical errors more effectively. It should be noted that the higher eigenmodes are of sufficient interests to the neutron transport community. Some Monte Carlo transport codes contain the FM to get the higher modes. The study on MPM shows that it is an alternative way to get the higher modes. 1
Physically, the elements of the FM are the probabilities of neutrons from one point causing fission in other points. This physical meaning is theoretically exact. In practice, the system is continuous, so it is impossible to get the “exact” FM. Space approximation is required, and the FM is obtained on discrete meshes. This FM is exact for the fundamental mode, and very fine meshes are required for the FM to get the higher modes. Therefore, the FM usually requires big memory for 3D big systems like 3D whole core problems. The Transfer Matrix (TM) was introduced for the general solution strategy of the MPM [15]. It is calculated by the neutron weight integrals before and after the power operation. It can keep the balance of fundamental mode, in which case it is identical to the FM. However, it can also keep the balance of the higher modes, which differentiates it from the FM. The calculation of the TM does not require space approximation, whereas the space discretization is used for the representation of the higher eigenmodes. Thus, the memory requirement can be much lower than that of the FM, especially for 3D big systems like 3D whole core problems. The purpose of this paper was to explain the differences between the TM and FM in detail, to establish that the two matrices are essentially different and that the TM requires further investigation. 2. The Fission Matrix Method 2.1 Review of the theory In this part, we will briefly introduce the theory of the Fission Matrix Method (FMM) based on the work by S.E. Carney, F.B. Brown, J. Colin, etc. [20-26]. The k-eigenvalue form of the neutron transport equation is: 1 E S r , M r, E, k 4
(1)
where M is the net loss operator defined by:
r , E , r, E r, E, M r, E, T ' r , E' E, ' r , E ', ' , dE ' d S
which contains the leakage, collision and scattering loss terms. S r
(2)
is the fission neutron
source defined by: ' r , E ' r , E ', ' . S r dE ' d F
(3)
E is the emission fission neutron energy spectrum, which is generally assumed to be isotropic, and independent of space and the neutron energy that causes the fission. The Green’s function for the neutron transport problem is defined by: 0 r, E, r r0 E E 0 , M G r 0 , E0 , 0
2
(4)
where the subscript “0” represents the initial state point. Based on the principle of linear superposition, the neutron flux can be represented by the Green’s function as: 1 E0 S r 0 G r 0 , E , 0 r, E, . d r r, E, 0 dE0 d 0 0 k 4
Multiplying Eq. (5) by F r , E
(5)
and integrating over all the energy and angular space, it
gives:
1 S r d r0S r0 H r0 r , k
where the kernel H r 0 r
(6)
represents the energy-angle averaged Green’s function:
dE d 0 r , E E0 G r 0 , E , 0 r, E, . H r 0 r dEd 0 0 F 4
(7)
Eq. (7) shows that the kernel H is the weighted integration of Green’s function over the initial and final energies and angles. The weight function is the initial fission neutron energy spectrum and the final fission neutron production. To get the Fission Matrix, the system should be divided into discrete spatial regions. Eqs. (6) and (7) integrated over the initial region J and final region I will give: SI
1 N FI , J S J , k J 1
(8)
where N is the total number of discrete spatial regions, and: S r 0 d r 0 , S I S r d r , r 0 RJ rRI S r0 FI ,J d r d r 0 H r0 r . rRI r 0 RJ SJ
S J
(9)
FI ,J is the element of the Fission Matrix, which is the number of fission neutrons produced in the region I due to one fission neutron starting in the region J. Eq. (8) can be rewritten in the matrix form:
S
1 FS, k
(10)
where S is a vector of length N that represents the fission source, F is a matrix of size N-by-N and is called the Fission Matrix, k is the eigenvalue of the Fission Matrix, and the biggest eigenvalue should be positive and equal to the fundamental eigenvalue of the system. 3
2.2 Considerations on higher modes calculation with FMM For the one-speed integral transport equation, all the higher modes exist, with discrete real eigenvalues and real eigenfunctions. For the multigroup transport equation and continuousenergy transport equation, it is conventional practice to assume that higher modes exist, with real eigenvalues and eigenfunctions, even though that has not been proven. The higher mode eigenfunctions can be solved with
1 Sn r d r 0 Sn r 0 H r 0 r , Kn
where K n denotes the n-th eigenvalue of Eq. (6), S n r
(11)
is the n-th eigenfunction, with
n=0 denoting the fundamental mode and n>0 the higher modes. Since the energy-angle averaged Green’s function H is not a symmetric function of its arguments, the forward eigenfunctions do not form an orthogonal basis set. Identical to Eq. (10), the higher eigenmodes can be determined by: Sn
1 Fn S n , kn
(12)
where the subscript n refers to the mode, with n=0 the fundamental mode. According to Eqs. (10) and (12), if Fi = Fj F (i,j=0,1,...,N-1) , all the N eigenmodes can be calculated with only one matrix F of size N-by-N. This requires: S0 r 0 FI ,J d r0 H r0 r rRI r 0 RJ S0, J Sn r 0 H r 0 r , n=0,1,...,N-1, d r d r 0 rRI r 0 RJ S n, J d r
(13)
which means all the fission source eigenmodes have the same distribution in one discrete mesh. If the spatial mesh used for the FM tally is fine enough, so that:
Sn r 0
1 for r V , n 0,1,, N 1, J
Sn, J VJ
(14)
then Eq. (13) is true and the FM can give all the N eigenmodes correctly. However, if the FM is tallied with a coarse mesh, Eq. (13) cannot be satisfied, and the corresponding FM for every eigenmode is different. In other words, if the FM is obtained with the n-th mode fission source as the weighting function, it can only give the n-th mode solution correctly. 3. The Modified Power Method and the Transfer Matrix 4
The MPM was originally proposed by Booth to obtain the first two eigenmodes at the same time and accumulate the convergence of the fundamental mode by subtracting the first mode from it [1-8]. Zhang implemented the MPM in continuous energy Monte Carlo simulation and extended it to obtain even higher modes with a general solution strategy and the TM [1215]. For completeness of the paper, the MPM will be briefly summarized.
3.1 Review of the theory According to Eqs. (6) and (11), the eigenvalue problem to be solved by the power method is:
ASn r d r 0 Sn r 0 H r 0 r K n Sn r , n=0,1,...,
K n , Sn
where A is the power operator and
(15)
are the eigenpairs of the system. It should be
noted that A is a continuous operator, and it is different from the F matrix of Section 2 which is one possible discretization of A. To solve this problem, an initial distribution is needed, that will finally converge to the fundamental mode if the power operation is applied continuously. To solve the first N eigenmodes with the MPM, N independent initial sources are needed: N 1
j cnj S n , j 0,..., N 1, n 0
(16)
N 1
A j cnj K n Sn , j 0,..., N 1, n 0
where j is the j-th source, and the higher modes of order no less than N are neglected. The idea of MPM is that the linear combination of these N independent sources can produce the first N eigenmodes:
0
x00 x0( N 1) N 1 S0 S N 1 , x( N 1)0 x( N 1)( N 1)
(17)
where xij is the linear combination coefficient, and it’s relation to cij is: 1
x00 x0( N 1) c00 c0(N 1) . x( N 1)0 x( N 1)( N 1) c(N 1)0 c(N 1)(N 1)
(17a)
It is noticed that for any sub-region Rm of the system, there is:
Rm
AS n dr K n S n dr , Rm
5
(18)
where n is the mode index and m is the subregion index. Substituting Eq. (17) into Eq. (18) gives: x00 Rm AS0 AS N 1 dr Rm A0 A N 1 dr x ( N 1)0 K0 K 0 S0 K N 1S N 1 dr S0 S N 1 dr Rm Rm 0
Rm
x0( N 1) x( N 1)( N 1) 0 K N 1
(19)
x00 x0( N 1) K 0 0 N 1 dr . x ( N 1)0 x( N 1)( N 1) 0 K N 1
0
If the whole system is divided into N sub-regions, and denoting Wmn A n dr and Rm
Vmn n dr , it will result in: Rm
W10 W1( N 1) x00 x0( N 1) W N 0 WN ( N 1) x(N 1)0 x(N 1)(N 1) V10 V1( N 1) x00 0 x0( N 1) K 0 , VN 0 VN ( N 1) x(N 1)0 x(N 1)(N 1) 0 K N 1
(20)
WX = VXK. The TM is introduced to help solve the linear combination coefficients. Rewriting Eq. (20) gives:
W = VXKX-1 = VX K VX V. -1
(21)
Denoting Smn Sn dr , and recalling Eq. (17), it gives: Rm
V10 V1( N 1) x00 x0( N 1) S10 S1( N 1) , V N 0 VN ( N 1) x(N 1)0 x(N 1)(N 1) S N 0 S N ( N 1) VX S.
Combining Eqs. (21) and (22) gives: 6
(22)
W SKS 1V.
(23)
P SKS 1 ,
(24)
Defining the TM as:
It can be noted that the eigenvalues of the TM are the eigenvalues of the system, and the eigenvectors of the TM are the integrals of the sources on the N sub-regions. The TM can be solved with the weight integrals: P = WV 1 ,
(25)
and then its eigenvectors can be used to solve the linear combination coefficients: X = V -1S, which will then be substituted into Eq. (17) to update the first N eigenmode solutions. Fig. 1 shows the flow chart of the solution strategy of the MPM.
7
(26)
Figure 1. Flow chart of the MPM.
3.2 The accumulation of the TM Due to the inherent statistical noises in Monte Carlo simulation, a TM based on just onecycle tally contains great statistical error, so it should be avoided. Therefore, the accumulation of TM is adopted to reduce the statistical error. In this work, the TM is accumulated from the 2nd cycle, and it is used to update the neutron weights at the end of every cycle.
3.3 Considerations on higher modes calculation with MPM The MPM requires modeling multiple sources explicitly. Instead of simulating more 8
neutrons with different source distributions, the weight vector is introduced for one neutron, with each element of the vector corresponding to one eigenmode. This will not affect the neutron tracking simulation, but different sources can be modelled at the same time when we look at different elements of the weight vector. This makes the application of the MPM in Monte Carlo more attractive. The neutron weights are all positive for the fundamental eigenmode, whereas both positive and negative weights should be modeled for the higher eigenmodes. During the Monte Carlo neutron transport simulation, neutrons with both positive and negative weights will distribute over the entire system and so the weight cancellation procedure is required. In this study, the fine mesh based weight cancellation is adopted. The fine mesh adopted in MPM is different from that in FMM. In MPM, the fine mesh is used to represent the higher eigenmode distributions. In FMM, the fine mesh is used for the FM tally. If the number of fine mesh is M, and N eigenmodes are modeled with the MPM (N << M), the array size for weight cancellation is (N*M), which is usually much smaller than the array size for tallying the FM, (M*M). To compare intuitively, the array size for FM tally and weight cancellation of the MPM with 8 modes are shown in Fig. 2, with the assumption that the number of meshes are the same for different dimensions.
Figure 2. The size of the array for FM tally or weight cancellation. 4. Analytical Solutions of a 1D1G Homogeneous Diffusion Problem To further compare the differences between the FM and the TM, the analytical solutions of a 1D1G homogeneous diffusion problem will be solved, and the 2-by-2 FM and TM will be discussed in this part. For a 1D1G homogeneous diffusion problem, the following equation should be solved:
9
D
d 2 x 1 a x f x , 2 dx k
(27)
where D is the diffusion coefficient, a is the absorption cross section, f is the fission cross section, is the average number of neutrons produced per fission, k is the multiplication factor and the eigenvalue of the problem, and x is the neutron flux. The left-hand side contains the leakage rate and the absorption rate. The right-hand side is the fission neutron source. The Green’s function for this problem can be solved by:
d 2G x0 , x D a G x0 , x x x0 , dx 2
(28)
and the final solution is: G x0 , x
a 1 K x x0 e , K . D 2 KD
(29)
To obtain the eigenvalue k and the corresponding neutron flux x , the power method is commonly adopted. By assuming an initial eigenvalue k(0) and an initial neutron flux (0) x , the neutron flux can be updated by:
(i 1) x G x0 , x
1 f (i ) x0 dx0 , i 0,1,. (i ) k
(30)
Multiplying both sides of Eq. (30) by f gives:
s (i 1) x
1 dx0 s (i ) x0 H x0 , x , i 0,1, , k (i )
(31)
where s x f x and H x0 , x f G x0 , x . The corresponding eigenvalue can be updated by:
s x dx , s x dx ( i 1)
k
( i 1)
(i )
i 0,1,.
(32)
4.1 The exact eigenvalues According to the diffusion theory, we already know that the first two harmonics of the neutron flux are:
10
x , 1 x c1 sin 2a a
0 x c0 cos
x , a x a,
(33)
where c0 and c1 are the normalization factors, and a is the half-thickness of the slab. According to Eqs. (31) and (32), the fundamental eigenvalue should be: k0
dx dx s x H x , x dx dx x G x , x dxs x dx x 0 0
0
0
0
f
0
f
0
f
a
x a
dx
e
2 4K 2a2 2
a
0
0
K x x0
2 Ka 2 2 2 f 8K a 1 e
f
cos x0 dx0 2a a x a cos 2a x dx
a
x0 a
2 KD
0
(34)
.
The 1st mode eigenvalue should be: k1
dx dx s x H x , x dx dx x G x , x dx s x dx x 0 1
0
0
0
1
f 2 KD
a
x a
dx
a
e
f
0
K x x0
Ka 2 2 2 f 4 K a 1 e
a
0
f 1
sin x0 dx0 x0 a a a x a sin a x dx
f 1
4 K 2a2 2
(35)
2
.
4.2 The analytical solution of the FM Suppose the slab is equally divided into 2 sub-regions, that is, x a, 0 as R1 and x 0, a as R2. According to the symmetry of the problem, the FM should be in the
following form: PFM where
11
,
(36)
fission neutrons produced in one region source neutrons in the same region
R1
dx dx0 s x0 H x0 , x R1
R1
fission neutrons produced in one region source neutrons in the other region
R2
dx dx0 s x0 H x0 , x R1
R1
The eigenvectors of the FM should be 1 1 ' and
dxs x
dxs x
,
(37) .
1 1 ' , as:
1 1 , 1 1
1 1 . 1 1
(38)
Accordingly, the eigenvalues of the FM should be: k0FM , k1FM .
(39)
4.2.1 The FM weighted with the fundamental mode fission source
Substituting the fundamental mode fission source into Eq. (37) gives:
R1
dx dx0 f 0 x0 f G x0 , x R1
R1
dx f 0 x
2 2 Ka f 8K a 2 Ka 1 e
2 4K 2a2 2
a
R2
R1
dx f 0 x
Ka Ka f 2 Ka e 1 e
a
2 4K 2a2 2
f
0
x a
dx
cos x0 dx0 2a 0 x a cos 2a x dx (40)
0
x0 a
2 KD
e
K x x0
,
dx dx0 f 0 x0 f G x0 , x R1
f 2 KD
a
x 0
dx
cos x0 dx0 2a 0 x a cos 2a x dx (41)
0
x0 a
e
K x x0
.
The eigenvalues can be calculated by substituting Eqs. (40) and (41) into Eq. (39):
12
k
FM 0
k1FM
2 2 2 2 Ka f 8K a 1 e
a
2 4K 2a2 2
,
2 2 Ka 2 Ka f 8 K a 4 Ka 1 e 1 e
2 4K 2a2 2
a
(42)
2
.
According to Eqs. (34), (35), and (42), it can be concluded that the fundamental mode eigenvalue of the FM is the same as the exact solution, whereas the 1st mode eigenvalue is not exact. As an example, the 1st mode eigenvalues and the error are shown in Fig. 3 with varying K and a. The cross sections are set as: a=0.20 cm-1, f=0.30 cm-1, D=0.33~0.83 cm. The errors of the eigenvalues are calculated as k1FM k1 105 .
Figure 3. Comparison of the 1st mode eigenvalues. 4.2.2 The FM weighted with the first mode fission source
Substituting the 1st mode fission source into Eq. (37) gives:
R1
dx dx0 f 1 x0 f G x0 , x R1
R1
dx f 1 x
2 2 2 Ka f 2 K a 1 e
2 K 2a2 2
a
R2
R1
dx f 1 x
2 2 Ka f 1 e
a 4 K 2a2 2
x a
dx
sin x0 dx0 a 0 x a sin a x dx (43)
0
x0 a
2 KD
e
K x x0
,
dx dx0 f 1 x0 f G x0 , x R1
f
0
.
13
f 2 KD
a
x 0
dx
sin x0 dx0 a 0 x a sin a x dx (44)
0
x0 a
e
K x x0
Substituting Eqs. (43) and (44) into Eq. (39), the eigenvalues can be calculated: k
FM 0
k1FM
2 2 2 Ka Ka f 4 K a 1 e 3 e
4 K 2a2 2
a
2 2 2 Ka f 4 K a 1 e
a
4 K 2a2 2
, (45)
2
.
According to Eqs. (34), (35), and (45), it can be concluded that the 1st mode eigenvalue of the FM now is exact. On the contrary, the fundamental mode eigenvalue is not exact. As an example, the fundamental mode eigenvalues and the error are shown in Fig. 4 with varying K and a. The cross sections are set as: a=0.20 cm-1, f=0.30 cm-1, D=0.33~0.83 cm. The errors of the eigenvalues are calculated as k0FM k0 105 .
Figure 4. Comparison of the fundamental mode eigenvalues.
4.2.3 Discussion
In the previous sections, the 2-by-2 FM is analyzed. However, in practice, the 2-by-2 coarse mesh is never used. The fine mesh is required for the FM to get the accurate results, even if only the first two eigenmodes are to be obtained. In addition, theoretically, if the first mode can be simulated and is used to calculate the FM, the resulting FM will give exactly the first mode solution. 4.3 The analytical solution of the TM 4.3.1 The determination of the TM
The integrals of the eigenfunctions before and after the power operation can be assumed as: C V 0 C0
C1 k0C0 , W C1 k0C0
According to the MPM, the TM can be calculated as: 14
k1C1 . k1C1
(46)
PTM WV 1 k C 0 0 k0C0
1
k1C1 C0 k1C1 C0
k k 0 1 k0 k1
C1 k0C0 C1 k0C0 2 k0 k1 2 . 2 k0 k1 2
The eigenvectors of the TM will be 1 1 ' and k0 k1 k0 k1 k0 k1 k0 k1
2 2 2 2
k0 k1 k0 k1 k0 k1 k0 k1
k1C1 1 C1 k1C1 2C0C1 C0
C1 C0
(47)
1 1 ' , as:
2 1 1 k0 , 2 1 1 2 1 1 k1 , 2 1 1
(48)
and the corresponding eigenvalues of the TM will be k0 and k1 . It means that the TM can give the exact eigensolutions of the system. 4.3.2 Discussion of the first mode solution The first mode eigenvector has both positive and negative values. To simulate it using Monte Carlo method, neutrons with both positive and negative weights should be adopted. Since the neutron positions are randomly sampled, it is almost impossible for two neutrons with different signs of weights to collide with nuclide in the same position. Therefore, the weight cancellation scheme should be developed, otherwise both the neutrons with positive weight and neutrons with negative weights will be distributed in the space according to the fundamental mode and the information of the higher mode will be lost. In this work, the mesh based weight cancellation is adopted [15]. To reduce the discretization error, fine mesh should be adopted. In this way, the MPM is identical to the FMM since they both require the fine meshes. But as shown in Fig. 2, the memory requirements for the fine mesh are different. Besides, the fine mesh information is not needed during the neutron transport simulation, and the memory requirement can be further reduced, as discussed in [15]. However, it should be made clear that the fine mesh based weight cancellation is not essential. There is exact weight cancellation, according to [8], in which case the mesh discretization error can be eliminated in the first eigensolution. In addition, the MPM can be also applied to deterministic calculations, in which case the weight cancellation is not required and the final converged higher mode solutions would be exact. 5. Numerical Tests In this section, a 1D fuel slab of size 400 cm was modeled, with black boundaries on both sides. The fuel composition is listed in Table 1. Eight eigenmodes were considered for this study. For the simulation of MPM, the slab was equally divided into eight sub-regions, and a fine mesh of size 1.0 cm was adopted for the weight cancellation of higher modes. All the 15
simulations were done with 200 inactive cycles, 600 active cycles, and 1,000,000 histories per cycle. Table 1. Composition of the fuel (units: atoms/barn-cm) Nuclide Density H-1 5.9347E-02 N-14 2.1220E-03 O-16 3.7258E-02 U-235 7.6864E-05 U-238 6.8303E-04 5.1 Comparison of the eigenvalues The eigenvalues of the different methods are shown in Fig. 5. For the simulation of FMM, different number of meshes were used, and the eigenvalues of the corresponding FM are shown and compared. “FM_n” indicates that the FM was tallied with n number of meshes. With the increase of the number of meshes, the eigenvalues of the FM increase gradually. This confirms that the FM can produce accurate higher mode solutions if the mesh is fine enough. In addition, the FM can always give the correct fundamental mode eigenvalue regardless of the number of meshes, because the FM was tallied with the fundamental mode neutron weights. For the MPM, 400 meshes are used for the weight cancellation, but the TM is only an 8-by-8 matrix. For the FMM, at least the 128-by-128 FM is needed to give results with similar precision. This confirms the difference between the TM and the FM.
Figure 5. Eigenvalues of different methods.
5.2 Using eigenvectors of the FM to modify the neutron weights 16
The FM can be used to accelerate the convergence of the fundamental fission source. In practice, the fundamental eigenvector of the FM is used to calculate a correction factor to be multiplied with the neutron weights. From Section 5.1, the FM can always give the correct fundamental mode solution, so it can be used to accelerate the fission source convergence with various meshes. During the simulation of MPM, the N eigenvectors of TM were used to calculate the combination factors according to Eq. (26) and modify the neutron weights accordingly at the end of every cycle. This is different from how the FM is used to correct the neutron weights with only one fundamental eigenvector. To clarify the differences, for the following tests, the first N eigenvectors of the FM were used to calculate the combination factors. If the number of meshes is greater than the number of modes, the linear least square solution of the following equation is solved: VX Φ,
(49)
where V R M N contains the weight integrals at the beginning of the cycle, Φ R M N contains the first N eigenvectors of the FM, X R N N contains the combination coefficients, and M N 8 in this case. The Shannon Entropy results are shown in Fig. 6, where “FM_n1_n2” represents the results of the FM based on n1 meshes while n2 modes are considered. Since the FM can always provide a correct fundamental mode solution, the Shannon Entropy increases quickly at the beginning for all the cases. However, comparing the results of “FM_8_8” and “TM_8_8”, the result with the FM converges very slowly after the first several cycles. This reveals that the coarse mesh based FM cannot provide exact higher mode eigenvectors. With an increase of the number of meshes, the FM can produce the first 8 modes accurately, and so the convergence of the Shannon Entropy becomes faster. In addition to the reduced cycle for the source convergence, there may be concerns about the actual simulation time of the MPM since multiple eigenmodes need to be explicitly simulated. To make it clear, the MPM simulations with different number of modes have been conducted, and the simulation times have been compared with the Original Monte Carlo simulation with the same number of histories and processes. The results are shown in Fig. 7. The simulation time of the MPM is not several times the Original Monte Carlo simulation time, but just ~5% higher per one more mode to be simulated. This is because all the eigenmodes are simulated at the same time by the introduction of the weight vector to one particle, and the time for adjusting the weights are much less than the time of particle tracking information sampling. To better illustrate the differences between the TM and FM, another test has been done, and the results are shown in Fig. 8. “FM_8_1” means the FM is tallied based on 8 meshes, and the fundamental mode of the FM is used to modify the neutron source. This is identical to how the FM is used to accelerate the source convergence, in which a correction factor is calculated for each mesh according to the fundamental eigenvector of the FM. The “FM_8_1” result shows efficient source convergence acceleration at the beginning, but the final steady value deviates from the others. This is because of the eigenvector of the coarse mesh based FM. To confirm that, the “FM_100_1” simulation has been done, in which case 100 uniform meshes are used to tally the FM, and the final steady value is then the same as that of “TM_8_8” and “Original”.
17
Figure 6. Shannon Entropy results with different methods.
Figure 7. The time increase with the number of modes for the MPM.
18
Figure 8. Shannon Entropy results with different methods. Table 2. The eigenvalues by different methods for 1D fuel slab problem. FM tallied with weights of different modes Tally TM w0 w1 w2 w3 w4 w5 w6 k0 k1 k2 k3 k4 k5 k6 k7
1.29268 (1.6) 1.28491 (7.1) 1.27229 (6.8) 1.25505 (7.1) 1.23344 (7.1) 1.20804 (7.1) 1.17915 (7.2) 1.14733 (7.4)
w7
1.29251
1.29268
1.29051
1.28938
1.28690
1.28836
1.28695
1.28747
1.28229
1.28475
1.26374
1.28492
1.28037
1.27667
1.27528
1.27479
1.27458
1.26814
1.27219
1.22604
1.20610
1.27238
1.26528
1.25918
1.25710
1.25644
1.24646
1.25463
1.16553
1.19339
1.14022
1.25501
1.24312
1.23732
1.23388
1.22133
1.23312
1.10715
1.08646
1.12154
1.05271
1.23350
1.22088
1.21152
1.19673
1.20723
1.05009
1.07392
1.10798
1.03867
0.99669
1.20829
1.19254
1.17447
1.17844
1.01183
0.99908
0.97886
1.02747
0.98234
0.94490
1.17943
1.15991
1.14694
0.98650
0.99390
0.96899
1.02104
0.96525
0.92705
0.89508
1.14814
5.3 The FM tallied with weights of different modes Since the neutron weights corresponding to N eigenmodes are tracked explicitly during the simulation of MPM, it is possible to tally the FM with weights of different modes. As illustrated in Section 2, the FM tallied with weights of mode n (FM_n) will give a correct n-th eigenmode solution. The following tests will confirm this. Table 2 shows the eigenvalues of the cycle tallies, the TM, and the FM with weights of different modes. The FM was tallied and accumulated from the beginning of active cycles 19
during the simulation of MPM when all modes are converged. From Table 2 it can be confirmed that the FM tallied with different mode weights can give the correct eigenvalue corresponding to that mode. 5.4 The results of the TM with a different number of modes Several simulations of the MPM with different numbers of modes have been done. The Shannon Entropy results are shown in Fig. 9. All the Shannon Entropy results nearly converge to the same value finally. Besides, when a larger number of modes is considered, generally the convergence of the Shannon Entropy is faster.
Figure 9. Shannon Entropy results of the MPM with different numbers of modes. The Eigenvalue results are illustrated in Fig. 10. When a greater number of modes is considered, more eigenvalues can be produced. The first several eigenvalues are always the same.
20
Figure 10. Eigenvalues of the TM with different numbers of modes. 5.5 The results of the TM with different fine meshes for weight cancellation As discussed previously, the fine mesh is essentially not required for the MPM. In this study, to process weight cancellation with high efficiency, a fine mesh average based weight cancellation procedure is adopted. Therefore, it is necessary to investigate how this fine mesh will affect the results. Several simulations of the MPM with 8 modes but with different fine meshes have been done. The eigenvalue results are shown in Fig. 11. The weight cancellation process with different numbers of meshes will affect the results of the MPM, which is expected. However, the pattern with decreasing the number of meshes is different from that of the FM shown in Fig. 5. There are very small differences between the results of the TM with more than 100 meshes for the weight cancellation, which confirms the reasonability and the accuracy of the 400-mesh based weight cancellation adopted in previous simulations of the MPM.
21
Figure 11. Eigenvalues of the TM with different meshes for the weight cancellation. 6. Conclusions The characteristics of the TM introduced in the MPM have been studied. Due to its similarity to the FM, the TM is usually mistakenly thought of as the same as the FM. In this paper, the FMM and the MPM have been reviewed, and the differences between the TM and FM have been discussed. From the theoretical view, the FM can produce the fundamental mode solution correctly regardless of the mesh size, but for higher mode solutions the fine mesh is required for the tally of the FM. On the other hand, the TM is calculated by solving the linear balance equation system of the N eigenmodes on N predefined sub-regions. The calculation of the TM is based on the coarse mesh. The fine mesh is adopted in this study for the weight cancellation of the higher modes, but it is not essential. The analytical solutions of a one-dimensional monoenergetic homogeneous diffusion problem with the 2-by-2 FM and TM have been calculated. The results confirm that the FM with such coarse mesh can only give one exact eigenmode solution, i.e., the coarse mesh FM is eigenmode dependent. The TM can always give the correct fundamental mode solution, but the higher mode solution depends on the weight cancellation precision. The numerical tests with a one-dimensional fuel slab problem have also been conducted. It is confirmed again that the FM can give more accurate results for the higher modes if a finer mesh is used. This is also confirmed by using the eigenvectors of the FM to update the eigenmodes. In addition, it is also proved that the FM is eigenmode dependent by the numerical results of the FM tallied with weights of different modes. On the other hand, by comparing the TM with different numbers of modes, it is demonstrated that under a certain condition of the weight cancellation, the TM can give all the eigenmode solutions correctly regardless of the number of eigenmodes and coarse meshes. Generally, it can make the convergence of the fundamental mode faster when more modes are considered. The impact of 22
the mesh average based weight cancellation on the results of the MPM has also been investigated, and it is found to be a little different from that of the FMM. Therefore, the TM is essentially different from the FM. Both can give the correct fundamental mode solution regardless of the mesh size. Both can be used to accelerate the convergence of the fission source. Whereas for the FM only the fundamental eigenvector is used, for the TM the first N eigenvectors are used. To get the accurate higher mode solutions, the FM requires very fine mesh and large memory for the tally of the FM, in which case the TM requires an accurate weight cancellation procedure and the memory requirement can be lower. Further investigation will be conducted on the exact weight cancellation and the applications of the MPM. Acknowledgements This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (No. 2017M2B2A9A02049916). References 1. T.E. Booth, Nucl. Sci. Eng. 143 (2003) 291-300. 2. T.E. Booth, Nucl. Sci. Eng. 154 (2006) 48-62. 3. J.E. Gubernatis, T.E. Booth, J. Comp. Phys. 227 (2008) 8508-8522. 4. T.E. Booth, J.E. Gubernatis, Monte Carlo Determination of Multiple Extremal Eigenpairs, LA-UR-07-7672, Los Alamos National Laboratory, 2008. 5. T.E. Booth, J.E. Gubernatis, Multiple Extremal Eigenpairs of Very Large Matrices by Monte Carlo Simulation, LA-UR-08-0043, Los Alamos National Laboratory, 2008. 6. T.E. Booth, J.E. Gubernatis, Improved Criticality Convergence via a Modified Monte Carlo Power Iteration Method, Proc. M&C 2009, Saratoga Springs, New York, May 3-7, 2009. 7. T.E. Booth, J.E. Gubernatis, Phys. Rev. E 80 (2009) 046704. 8. T.E. Booth, J.E. Gubernatis, Nucl. Sci. Eng. 165 (2010) 283-291. 9. T. Yamamoto, Ann. Nucl. Energy 36 (2009) 7-14. 10. B. Shi, B. Petrovic, Ann. Nucl. Energy 38 (2011) 781-787. 11. B. Shi, B. Petrovic, Nucl. Sci. Eng. 172 (2012) 138-150. 12. P. Zhang, H. Lee, D. Lee, Extension of Tom Booth’s Modified Power Method for Higher Eigen Modes, Trans. KNS Spring Meeting, Jeju, Korea, May 7-8, 2015. 13. P. Zhang, H. Lee, D. Lee, The Implementation of Modified Power Iteration Method in Continuous Energy Monte Carlo Simulation, Proc. ANS 2015 Annual Meeting, San Antonio, TX, June 7-11, 2015. 14. P. Zhang, H. Lee, D. Lee, Stabilization Technique of Modified Power Iteration Method for Monte Carlo Simulation of Neutron Transport Eigenvalue Problem, Proc. M&C 2015, Nashville, TN, April 19-23, 2015. 15. P. Zhang, H. Lee, D. Lee, J. Comp. Phys. 305 (2016) 387-402. 16. P. Zhang, H. Lee, D. Lee, J. Comp. Phys. 320 (2016) 17-32. 17. P. Zhang, H. Lee, D. Lee, Nucl. Eng. Tech. 49 (2017) 17-28. 18. P. Zhang, H. Lee, D. Lee, J. Comp. Phys. 344 (2017) 440-450. 19. K Tuttelberg, J. Dufek, Ann. Nucl. Energy 72 (2014) 151-155. 20. S.E. Carney, F.B. Brown, B.C. Kiedrowski, W.R. Martin, Ann. Nucl. Energy 73 (2014) 423-431. 21. F.B. Brown, S.E. Carney, B.C. Kiedrowski, W.R. Martin, Fission Matrix Capability for MCNP Monte Carlo, Proc. SNA + MC 2013, Paris, France, October 27-31, LA-UR-13-26962, Los Alamos National Laboratory, 2013. 22. J. Colin, Veit Max D., Eigenfunction Decomposition of Reactor Perturbations and Transitions Using MCNP Monte Carlo, report LA-UR-13-26449, talk LA-UR-13-26464, Los Alamos National Laboratory, 2013. 23. S.E. Carney, F.B. Brown, B.C. Kiedrowski, W.R. Martin, Higher-Mode Applications of Fission Matrix Capability for MCNP, report LA-UR-13-27078, talk LA-UR-13-26615, Los Alamos National Laboratory, 2013. 24. F.B. Brown, S.E. Carney, B.C. Kiedrowski, W.R. Martin, Fission Matrix Capability for MCNP, Part I – Theory, Proc. M&C 2013, Sun Valley, Idaho, May 5-9, 2013. 25. F.B. Brown, S.E. Carney, B.C. Kiedrowski, W.R. Martin, Fission Matrix Capability for MCNP, Part II – Applications, Proc. M&C 2013, Sun Valley, Idaho, May 5-9, 2013. 23
26. S.E. Carney, F.B. Brown, B.C. Kiedrowski, W.R. Martin, Fission Matrix Capability for MCNP Monte Carlo, LA-UR-12-24533, Los Alamos National Laboratory 2012. 27. J. Miao, B. Forget, K. Smith, Ann. Nucl. Energy 92 (2016) 81-95. 28. M. Nowak, J. Miao, E. Dumonteil, B. Forget, A. Onillon, K. S. Smith, A. Zoia, Ann. Nucl. Energy 94 (2016) 856-868. 29. T. M. Sutton, Nucl. Sci. Eng. 185 (2017) 174-183.
24