- Email: [email protected]
A new nonlinear iterative method for SPN theory
Annals of Nuclear Energy 110 (2017) 920–927
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
A new nonlinear iterative method for SPN theory Pan Qingquan a,⇑, Lu Haoliang b, Li Dongsheng b, Wang Kan a a b
Tsinghua University, Department of Engineering Physics, Beijing 100084, China China Nuclear Power Technology Research Institute, Reactor Engineering and Safety Research Center, Shenzhen, Guangdong 518031, China
a r t i c l e
i n f o
Article history: Received 18 May 2017 Received in revised form 13 July 2017 Accepted 25 July 2017
Keywords: SPN theory Nonlinear iteration Numerical instability
a b s t r a c t The SPN theory shows a great potential for solving the neutron transport equation in the Next Generation Reactor Physics Calculation. However, the numerical instability problem will be encountered in some cases when solving the SPN equation with the conventional nonlinear iterative method. In order to solve the problem, a new coupling corrective relationship which focuses on the angular coupling is put forward, and the new nonlinear iterative method has more margin to adjust the coupling corrective coefficients with that new relationship. As an example, the code named NLSP3 for solving the SP3 equation with the new nonlinear iterative method is developed, and the numerical results demonstrate the feasibility of the new method. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction The conventional nonlinear iterative method (Smith, 1983) is based on the nodal method and nonlinear iterative strategy, which is characterized by the use of the simplified CMFD formula of the neutron diffusion equation containing the nodal coupling corrective coefficients in global calculation. The nodal coupling corrective coefficients determined by solving the two-nodal problem at the interfaces of two adjacent nodes will be updated during the nonlinear iterative process. The method is mainly used to deal with the diffusion equation, however, due to it’s high efficiency, the method is expected to be used in the transport equation, especially in the SPN equation. The SPN theory (Chao, 2015) is the simplified spherical harmonics method for solving the neutron transport equation. It was widely studied (Yu et al., 2014; Yu, 2014) and it shows a great potential in the Next Generation Reactor Physics Calculation because of its better accuracy than diffusion theory. Adopting the conventional nonlinear iterative method to solve the SPN equation, the numerical instability problem is found during the coupling correction of the high-order flux since the high-order coupling corrective coefficients are too large (Yu et al., 2014; Yu, 2014). Even though in other nodal approaches (Berkert and Grundmann, 2008), the numerical instability problem is still existing when CMFD formula is used to speed up the calculation on the SPN equation. So, it is very important to put ⇑ Corresponding author at: Department of Engineering Physics, Tsinghua University, Haidian District, Beijing 100084, China. E-mail addresses: [email protected] (Q. Pan), [email protected] (H. Lu), [email protected] (D. Li), [email protected] (K. Wang). http://dx.doi.org/10.1016/j.anucene.2017.07.030 0306-4549/Ó 2017 Elsevier Ltd. All rights reserved.
forward a new nonlinear iterative method for the SPN theory. In this paper, a new nonlinear iterative method is discussed with the using of a new coupling corrective relationship. The new coupling corrective relationship focuses on the coupling at the angle, and seems a bit more complicated than the conventional one. Except for the calculation at the vacuum boundary, the coupling corrective coefficients have many ways to calculate, which brings in margin to deal with the numerical instability under adequate mathematical derivation or justification. The code named NLSP3 is developed to achieve the new method within the SP3 theory, and the numerical results demonstrate the feasibility of the new method. In this paper, a review summary of Chao’s derivation for the SPN theory is given in Section 2, where the necessary derivations and equations are outlined. And the new coupling corrective relationship for SPN theory is generally introduced. Section 3 presents the special case for SP3 with the new method, and some advices are given after the discussion about the norm of the iterative matrix. Section 4 provides the numerical results and comparisons. The 2D-IAEA transport benchmark and the 3D-TAKEDA benchmark are calculated by the code NLSP3, and the feasibility of the new method is proved. Section 5 makes a conclusion about the work. 2. The new method for SPN theory
2.1. The SPN theory The SPN equation was firstly proposed without rigorous mathematical derivation (Gelbard, 1960). Later, the asymptotic theory
921
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927
approximation and variational method (Brantley and Larsen, 2000) were introduced to deduce a more rigorous equation in SP3. Those year, Professor Chao has done a lot of jobs (Chao, 2015) about SPN theory. In this paper, Chao’s jobs are considered as a reference and a brief introduction about his derivational process is given. Under the case of mono-energetic isotropic scattering, the neutron transport equation shows as follow,
X rwðX; rÞ þ Rt wðX; rÞ ¼
QðrÞ 4p
ð2:1Þ
And the neutron flux are divided into the even and odd parity angular flux,
wðX; rÞ ¼ wE ðX; rÞ þ wO ðX; rÞ
ð2:2Þ
With the expansion in orthogonal spherical harmonics and the assumption of cylindrically symmetric for the azimuth angular dis/n ðrÞ tribution of the nth flux moments, and Xn ¼ jr r/n ðrÞj ; we have,
2n þ 1 Pn ðX Xn Þ/n 4p
R
wE ðX; rÞ ¼
n¼ev en
wo ðX; rÞ ¼
1
Rt
R
n¼ev en
ð2:3Þ
2n þ 1 ½Pn ðX Xn ÞX Xn ½Xn r/n ðrÞ 4p
J n ¼
Z
b Þð n b XÞwE ðX; rÞdX bÞ Pn ðX n Pn ðX n bn X>0 bn X>0 b XÞwo ðX; rÞdX ðn
ð2:11Þ
The final expression for the nth net current as follows.
b Xn2 Þ nðn 1Þ Pn1 ð n b r/n2 ðrÞ n b Xn2 ð2n þ 1Þð2n 1ÞRt n b Xnþ2 Þ ðn þ 1Þðn þ 2Þ Pnþ1 ð n b r/nþ2 ðrÞ n b ð2n þ 1Þð2n þ 3ÞRt n Xnþ2 " # b Xn Þ b Xn Þ 1 ðn þ 1Þ2 Pnþ1 ð n n2 Pn1 ð n b þ n b Xn b Xn ð2n þ 1ÞRt 2n þ 3 2n 1 n n
Jn ¼
r/n ðrÞ ð2:12Þ With the definition of the incoming and outgoing partial currents, such as Eqs. (2.13) and (2.14), the whole SPN theory was deduced.
1 J n ¼ Un J n 2
ð2:4Þ
In order to derive the net current, we integrate the product of X and odd parity angular flux over the solid angle space,
Z
Un ¼
Z J þn þ J n b Þð n b XÞwE ðX; rÞdX Pn ðX n ¼ 2 bn X>0
ð2:13Þ
ð2:14Þ
Z
J¼
2.2. The new coupling corrective relationship
X wo ðX; rÞdX
¼
2k þ 1
1
R Rt k¼ev en 4p
Z
X Pk ðX Xk Þ½X r/k ðrÞdX
ð2:5Þ
With the use of the Legendre polynomial recursion relations, we have
1 ðk 1Þk d0;k þ d0;k2 r/k ðrÞ 2k 1 Rt k¼ev en 3 1 ¼ ½r/0 ðrÞ þ 2r/2 ðrÞ 3Rt
J¼
1
R
ð2:6Þ
where, dm;n is the Dirichlet function. In order to derive the differential equations of nth flux, substitute Eq. (2.3) into (2.1) and integrate it over the angular space,
2k þ 1 k¼ev en 4p
R
Z
1 Pn ðX Xn ÞPk ðX Xk Þ X r X r Rt /k ðrÞ dX
Rt
ð2:7Þ
¼ Qdn0 Making use of Eqs. (2.3) and (2.4), we have,
Xr
1 X r/k ðrÞ ¼ X r ðX Xk ÞðXk r/k ðrÞÞ Rt Rt r/k ðrÞ ¼ ðX Xk ÞX r Xk 1
Rt
ð2:8Þ
kþ1 k kþ1 k J kn;guþ ¼ Dk;FDM n;guþ f n;gu /n;g f n;guþ /n;g kþ1 k kþ1 k Dk;NOD n;guþ f n;gu /n;g þ f n;guþ /n;g
After some variation, SPN equation can be get: 2 2 nðn1Þ 2n þ2n1 ð2nþ1Þð2n1Þ Rt r /n2 ðrÞ ð2nþ3Þð2n1ÞRt r /n ðrÞ 2
2 ðnþ1Þðnþ2Þ ð2nþ1Þð2nþ3Þ Rt r /nþ2 ðrÞ þ Rt /n ðrÞ ¼ Q dn0
ð2:9Þ
The boundary conditions are provided by the relation between the incoming and outgoing partial current. And the nth partial curb at a surface is defined as, rent with the normal vector n b X½wE ðX; rÞ wo ðX; rÞ jbn X>0 ðX; rÞ ¼ n
The nonlinear iterative method is using the results with higherorder precision to correct the Coarse Mesh Finite Difference formula (CMFD) which has a lower-order precision. And the coupling corrective relationship characterizes the difference of the net neutron current calculated between the higher-order nodal method and CMFD. Therefore, the correction process is actually to import a perturbation which is exactly the coupling corrective coefficients to the CMFD formula (Fu, 2001). In that way, the solutions of the CMFD equation will achieve the accuracy of the higher-order node method. When the perturbation is too large, it will make the iterative matrix be a morbid matrix and the nonlinear iteration not converges, which exact is the numerical instability problem. Fig. 1 is the nonlinear iterative calculation strategy, consisting of source iteration and SANM coupling corrective coefficients updating, where the source iteration is divided into internal iteration and outer iteration. In the case of internal iteration, the source term is regarded as a constant, so the iteration is actually solving a linear system of equations. At that situation, the numerical instability problem occurs in the internal iteration, so it is necessary to analyze the iterative matrix of the internal iteration. The conventional coupling corrective relationship is given as Eq. (2.15), it calculates the net current only under the nth flux itself, without the coupling relationship between the different order flux.
ð2:10Þ
Make use of nth even order Legendre polynomial, and integrate Eq. (2.10) over the angular space.
ð2:15Þ
where, ‘‘k” is the number of the node, ‘‘g” is the number of the group. Normally, ‘‘n’ is set to be an odd number. We regard the nth-order flux as the different description about the angle distribution, so, the angular flux is the superposition of all the nth order flux. In order to solve the SPN equation with the current multigroup diffusion codes, The SPN equation can be written as (N + 1)/2 coupled differential equations on the each nth even order flux after some mathematical derivation and substitution, such as the composite form (Zhang, 2013). Adopting the conventional coupling corrective relationship to solve the SPN equation, the numerical
922
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927
2 2 Dk2;g r2 /k2;g ðrÞ þ Rkr2;g /k0;g ðrÞ Rkr0;g /k0;g ðrÞ ¼ Sk0;g 5 5
ð3:1bÞ
where,
Dk0;g ¼
1 3Rkt;g
; Dk2;g ¼
9 5
9 35Rkt;g
; Rkr0;g ¼ Rkt;g Rks;gg0 ;
4 5
Rkr2;g ¼ Rkt;g Rks;gg0
ð3:1cÞ
The Source term expression,
h i G Sk0;g ðrÞ ¼ R Rk0;g0 g /k0;g ðrÞ 2/k2;g ðrÞ g 0¼1 g 0 –g
þ Fig. 1. Nonlinear iterative strategy.
k
kþ1 k Dk;NOD 0;guþ ðf 0;gu /0;g þ f 0;guþ /0;g Þ kþ1 kþ1 Dk;FDM 2;guþ ðf 2;gu /2;g
kþ1 kþ1 Dk;NOD 2;guþ ðf 2;gu /2;g
þ f 2;guþ /k2;g Þ
ð3:2aÞ
J 2 ðrÞ ¼
21 M 1* 3 M / ðrÞ n J M / 2 ðrÞ 80 2 2 80 0
ð3:2bÞ
k
X
kþ1 k Dk;FDM n;guþ ðf n;gu /n;g f n;guþ /n;g Þ
u¼x;y;z
kþ1 k Dk;NOD n;guþ ðf n;gu /n;g þ f n;guþ /n;g Þ
kþ1
k
kþ1
k
1 Duk
ðJ k0;guþ J k0;gu Þ þ Rkr0;g Uk0;g 2Rkr0;g Uk2;g ¼ Sk0;g
1 Duk
ðJ k2;guþ J k2;gu Þ þ Rkr2;g Uk0;g 25 Rkr0;g Uk2;g ¼ 25 Sk0;g
kþ1 k Dk;FDM ðNþ1Þ=2;guþ ðf ðNþ1Þ=2;gu /ðNþ1Þ=2;g f ðNþ1Þ=2;guþ /ðNþ1Þ=2;g Þ
J kn;guþ ¼ Dkn;g
kþ1 k Dk;NOD ðNþ1Þ=2;guþ ðf ðNþ1Þ=2;gu /ðNþ1Þ=2;g þ f ðNþ1Þ=2;guþ /ðNþ1Þ=2;g Þ
n ¼ 0; 2; . . . ; ðN þ 1Þ=2 ð2:16Þ
With the new coupling corrective relationship, the CMFD formula seems a bet more complicated than the conventional one. For simplicity we will take the special case of SP3 as an example, and advices for how to calculate the coupling corrective coefficients will be given through mathematic derivation in next section.
@/kn;g ðuÞ @u
ð3:3bÞ
h i h i G 1 k G Sk0;g ¼ R Rk0;g0 g Uk0;g 2Uk2;g þ vg R mRkf;g0 Uk0;g 2Uk2;g 0 ¼1 g 0 ¼1 g K eff 0
ð3:3cÞ
g –g
Ukn;g ¼
1 Vk
Z V
/kn;g ðrÞdV
ð3:3dÞ
According to Eq. (2.16), we have the specific expression of the new coupling corrective relationship for SP3 equations, kþ1
k
kþ1
k
k;NOD kþ1 kþ1 k k J k0;guþ ¼ Dk;FDM 1;guþ ðf 0;gu /0;g f 0;guþ /0;g Þ D1;guþ ðf 0;gu /0;g þ f 0;guþ /0;g Þ
3. The special case for SP3
kþ1
k
kþ1
k
k;NOD kþ1 kþ1 k k Dk;FDM 3;guþ ðf 2;gu /2;g f 2;guþ /2;g Þ D3;guþ ðf 2;gu /2;g þ f 2;guþ /2;g Þ kþ1
k
kþ1
k
k;NOD kþ1 kþ1 k k J k2;guþ ¼ Dk;FDM 2;guþ ðf 2;gu /2;g f 2;guþ /2;g Þ D2;guþ ðf 2;gu /2;g þ f 2;guþ /2;g Þ
3.1. The SP3 equations
kþ1
k
kþ1
k
k;NOD kþ1 kþ1 k k Dk;FDM 4;guþ ðf 0;gu /0;g f 0;guþ /0;g Þ D4;guþ ðf 0;gu /0;g þ f 0;guþ /0;g Þ
The SP3 equation can be got when N is set to be 3 in the SPN equation. After deformation, a similar form to the neutron diffusion equation can be obtained and the same nodal method can be used to solve the SP3 equation, such as the nonlinear iterative method. Assuming that the neutron scattering and the source term are isotropic, and the high order scattering between the nth groups is zero, then the specific form of the SP3 equation and the boundary conditions show as follow,
Dk0;g r2 /k0;g ðrÞ þ Rkr0;g /k0;g ðrÞ 2Rkr0;g /k2;g ðrÞ ¼ Sk0;g
ð3:3aÞ
Where,
for
1 M 1* 3 M / ðrÞ n J M / 0 ðrÞ 4 0 2 16 2
u¼x;y;z
kþ1
ð3:1dÞ
g ¼1
J 0 ðrÞ ¼
X
k
k
i
The neutron balance equation of the node can be obtained by integrating the SP3 equations in node ‘‘k”, and u as a reference direction as follow,
k f 2;guþ /k2;g Þ
kþ1
h
3.2. The new method in SP3 equations
k
kþ1 k J kn;guþ ¼ Dk;FDM 0;guþ ðf 0;gu /0;g f 0;guþ /0;g Þ kþ1
G
vkg R mRkf;g0 /k0;g0 ðrÞ 2/k2;g0 ðrÞ 0
At the boundary, we have the neutron flux and the current under the Marshak boundary condition as follow, where, the superscript ‘‘M” represents the amounts at the node boundary.
instability problem is found during the coupling correction of the higher-order flux, because the higher-order flux is very small which leads the higher-order coupling corrective coefficients to be too large. So, we redefine the coupling corrective relationship and their coupling at angles is considered. So, the new coupling corrective relationship shows as follow, kþ1
1 K eff
ð3:1aÞ
ð3:4Þ Plugging Eq. (3.4) into Eq. (3.3), then we have the CMFD formula as Eq. (3.5), kþ1 k1 k1 kþ1 k kþ1 k1 k1 k k ak0;gx f kþ1 0;gx /0;gx þ b0;gx f 2;gx /2;gx þ c0;gxþ f 0;gxþ /0;gx þ q0;gx f 2;gxþ /2;gxþ þ kþ1 k1 k1 kþ1 k kþ1 k1 k1 k k ak0;gy f kþ1 0;gy /0;gy þ b0;gy f 2;gy /2;gy þ c0;gyþ f 0;gyþ /0;gy þ q0;gy f 2;gyþ /2;gyþ þ kþ1 k1 k1 kþ1 k kþ1 k1 k1 k k ak0;gz f kþ1 0;gz /0;gz þ b0;gz f 2;gz /2;gz þ c0;gzþ f 0;gzþ /0;gz þ q0;gz f 2;gzþ /2;gzþ þ
C k1;gz /k0;g þ C k3;gz /k2;g ¼ Sk0g ð3:5aÞ
923
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927 kþ1 kþ1 k1 k1 k kþ1 k1 k1 k k ak2;gx f 2;gx /kþ1 2;gx þ b2;gx f 0;gx /0;gx þ c2;gxþ f 2;gxþ /2;gx þ q2;gx f 0;gxþ /0;gxþ þ kþ1 kþ1 k 2;gy f 2;gy /2;gy
a
kþ1 kþ1 k 2;gz f 2;gz /2;gz
a
þ
kþ1 kþ1 bk2;gy f 0;gy /0;gy
þ
kþ1 kþ1 bk2;gz f 0;gz /0;gz
k1 k1 k 2;gyþ f 2;gyþ /2;gy
þc
k1 k1 k 2;gzþ f 2;gzþ /2;gz
þc
k1 k1 k 2;gy f 0;gyþ /0;gyþ þ
þq
k1
k1 þ qk2;gz f 0;gzþ /0;gzþ þ
C k2;gz /k2;g þ C k4;gz /k0;g ¼ 25 Sk0g ð3:5bÞ
"
where,
C k1;gz ¼ Rk0;g þ Rk0;gx þ Rk0;gy þ Rk0;gz
ð3:6aÞ
C k3;gz ¼ 2Rk0;g þ Rk2;gx þ Rk2;gy þ Rk2;gz
ð3:6bÞ
C k2;gz ¼ Rk2;g þ Lk2;gx þ Lk2;gy þ Lk2;gz
ð3:6cÞ
C k4;gz
2 ¼ Rk0;g þ Lk0;gx þ Lk0;gy þ Lk0;gz 5
k 0;gu
1 k;FDM ¼ D þ Dk;NOD 1;guþ Duk 1;guþ
k 2;gu
1 k;FDM ¼ D þ Dk;NOD 2;guþ Duk 2;guþ
a a
A2
2
#
½U ¼ 4
ð3:6eÞ
ð3:6fÞ
ð3:6gÞ
bk2;gu
1 k;FDM ¼ D4;guþ þ Dk;NOD 4;guþ Duk
ð3:6hÞ
k 0;gu
c
1 k;FDM ¼ D1;gu Dk;NOD 1;gu Duk
ck2;gu
1 k;FDM ¼ D2;gu Dk;NOD 2;gu Duk
ð3:6jÞ
qk0;gu ¼
1 k;FDM D3;gu Dk;NOD 3;gu Duk
ð3:6kÞ
qk2;gu ¼
1 k;FDM D4;gu Dk;NOD 4;gu Duk
ð3:6lÞ
3
As þ kF
eff
25 ðAs þ kF Þ
"
5½U ¼
eff
S
# ð3:8aÞ
25 S
where,
U ¼ colðU1 ; U2 ; ; Ug ; ; UG ; Þ Ug ¼ colðU1g ; U2g ; ; Ukg ; ; UNg Þ
ð3:8bÞ
Ukg ¼ colðUk0;g ; Uk2;g Þ 2
1 k;FDM ¼ D þ Dk;NOD 3;guþ Duk 3;guþ
AG1 s
3
0 6 12 6A 6 s 6 13 As ¼ 6 6 As 6 6 4
A21 s
A31 s
0
A32 s
A23 s
0
7 7 AG2 s 7 7 G3 7 As 7 7 7 5
A1G s
A2G s
A3G s
ð3:8cÞ
0
Asg g ¼ diagðR1g g ; 2R1g g ; R2g g ; 2R2g g ; ; Rkg g ; 2Rkg g ; ; RNg g ; 2RNg g Þ 2 ð3:6iÞ
And,
Rk0;gu ¼
1 k;FDM k;FDM k;NOD D1;guþ Dk;NOD 1;guþ þ D1;gu þ D1;gu Duk
ð3:7aÞ
Rk2;gu ¼
1 k;FDM k;FDM k;NOD D3;guþ Dk;NOD 3;guþ þ D3;gu þ D3;gu Duk
ð3:7bÞ
1 k;FDM k;FDM k;NOD ¼ D4;guþ Dk;NOD 4;guþ þ D4;gu þ D4;gu Duk
ð3:7cÞ
1 k;FDM k;FDM k;NOD D2;guþ Dk;NOD 2;guþ þ D2;gu þ D2;gu Duk
ð3:7dÞ
Lk2;gu ¼
A0
ð3:6dÞ
bk0;gu
Lk0;gu
The nonlinear iterative method chooses the source iteration to solve the CMFD equation, in which the numerical instability occurs during the internal iteration. Within the internal iterative process, the calculation of the sources bases on the last results of the outer iteration, that is, the sources are treated as constants. In the case of ‘‘N” nodes and ‘‘G” energy groups, Eq. (3.5) can be rewritten in matrix form, as follow,
In order to solve the CMFD formula, we need the values of the coupling corrective coefficients which should be solved by Eq. (3.4). However, there are two coupling coefficients in each relationship while only one net current can be used. So, except for the calculation at the vacuum boundary, the coupling corrective coefficients have many ways to calculate, which brings in the margin to deal with the numerical instability under adequate mathematical derivation or justification. Then, some advices for how to confirm the exact coefficients are given through numerical derivation.
F 11
F 21
F G1
ð3:8dÞ 3
6 6 F 12 F¼6 6 4
7 F G2 7 7 7 5
F 1G
F 2G
F GG
F 22
ð3:9Þ
F g g ¼ diagðvg mR1fg ; 2vg mR1fg ; vg mR2fg ; 2vg mR2fg ; ; vg mRkfg ; 2vg mRkfg ; ; vg mRNfg ; 2vg mRNfg Þ And,
A0 ¼ M 0 þ N 0
ð3:10aÞ
A2 ¼ M 2 þ N 2
ð3:10bÞ
M0 ¼ diagðM 01 ; M02 ; M 03 ; ; M0g ; ; M 0G Þ
ð3:11aÞ
N0 ¼ diagðN01 ; N02 ; N03 ; ; N0g ; ; N0G Þ
ð3:11bÞ
M2 ¼ diagðM 21 ; M22 ; M 23 ; ; M2g ; ; M 2G Þ
ð3:11cÞ
N2 ¼ diagðN21 ; N22 ; N23 ; ; N2g ; ; N2G Þ
ð3:11dÞ 3
2 1 1 6 C 1;g C 3;g 6 6 C 21;g 2 6 6 C 3;g 6 0 Mg ¼ 6 6 6 6 6 C k1;g C k3;g 4
7 7 7 7 7 7 7 7 7 7 7 5
N N C 1;g C 3;g
ð3:12aÞ
924
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927
3
2
a10;gx
6 b10;gx a b10;gy a 6 2 6 m0;gx q20;gx b10;gz 2 2 2 2 a a 6 b b 0;gx 0;gy 0 0;gx 0;gy Ng ¼ 6 a20;gz k 6 a0;gx bk0;gx ak0;gy bk0;gy ak0;gz bk0;gz 6 mk0;gz qk0;gz 6 k k k k 4 m0;gy q0;gy m0;gz q0;gz mN0;gz qN0;gz mN0;gy qN0;gy mN0;gz qN0;gz 1 0;gy
b20;gz
1 0;gz
3
2 1 1 6 C 4;g C 2;g 6 C 24;g 2 6 6 C 2;g 2 Mg ¼ 6 6 6 6 4 C k4;g C k2;g
7 7 7 7 7 7 7 7 5
ð3:12bÞ
N N C 4;g C 2;g
3
2 b12;gx
6 a22;gz a12;gx b12;gy a12;gy b12;gz a12;gz 6 2 6 q2;gx m22;gx 2 2 2 2 2 6 b2;gx a2;gx b2;gy a2;gy b2;gz N 2g ¼ 6 6 bk2;gx ak2;gx bk2;gy ak2;gy bk2;gz ak2;gz 6 qk2;gz mk2;gz 6 k k 4 q2;gy mk2;gy q2;gz mk2;gz qN2;gz mN2;gz qN2;gy mN2;gy qN2;gz mN2;gz
In the internal iterative process, the Eq. (3.8) can be rewritten as Eq. (3.13),
"
#
"
A0 ½U ¼ A2
M0 þ N 0 2
M þN
#
"
½U ¼
2
#
M0 M
2
"
½U þ
N0
"
#
N2
½U ¼
S
#
25 S ð3:13Þ
The iterative formula for calculating CMFD equation show as follow,
" ½UðKþ1Þ ¼
M0
#1 "
S 25 S
M2
#
"
M0
#1 "
M2
N0 N2
# ½UðKÞ
ð3:14Þ
M¼
M0
#
"
;N ¼
M2
N0 N2
# then E ¼
"
M0
#1 "
M2
N0 N2
# ¼ M 1 N
ð3:15Þ
kEk ¼ kM
Nk < 1
ð3:16Þ
And we know the relationship for all the matrix norm as Eq. (3.17),
kT 1 k kTk P kT 1 Tk ¼ 1
ð3:17Þ
ð3:18Þ
So, We can get,
kMk P
kNk P
1 kN1 k 1 kN1 k
ð3:21Þ
Because the equivalence of matrix norm, we can choose each one from the three kinds of matrix norm. So, there are many ways to conduct the Eq. (3.21). As an example, take the row-norm into consideration, and make a more stringent treatment, then we can easily get a relationship, such as Eq. (3.22) as follow. max fjC nþ1 j þ jC nþ3 jg P max R ðjakn;gu j þ jbkn;gu j þ jmkn;gu j þ jqkn;gu j 16i62NG
16i62NG
u¼x;y;z
ð3:22Þ
3Rk0;g P R ðjak0;gu j þ jbk0;gu j þ jmk0;gu j þ jqk0;gu j þ jRk0;gu j þ jRk2;gu jÞ ð3:23aÞ 2 5
Rk2;g þ Rk0;g P R ðjak2;gu j þ jbk2;gu j þ jmk2;gu j þ jqk2;gu j þ jLk0;gu j þ jLk2;gu jÞ u¼x;y;z
ð3:23bÞ where, n = 0,2. And ‘‘k” represent the number of node, ‘‘g” represent the energy group. Analyzing the inequality Eq. (3.23), the left side of the inequality is a constant while the right side is given by the pseudo-diffusion coefficients Dk;FDM n;guþ and the coupling corrective
According to Eqs. (3.16) and (3.17), we have,
kN1 k kMk P kN1 Mk ¼ kE1 k P 1
" " # # 0 0 M N P 2 M2 N
u¼x;y;z
We know that only if the norm of the iterative matrix is less than one can the calculation reaches convergence. So we should have, 1
7 7 7 7 7 7 7 7 5
The more explicit expression in each nodes of one energy group show as follows,
Let,
"
7 7 7 7 7 7 7 7 5
ð3:19Þ
ð3:20Þ
If kMk P kNk; that is Eq. (3.21), the iteration will reaches into convergence.
coefficients Dk;NOD n;guþ . Because the pseudo-diffusion coefficients are constants in each node, the coupling corrective coefficients should be as small as possible to satisfy the inequality. However, there are two coupling coefficients in each relationship while only one net current can be used, so there is a lot of margin to adjust the coupling corrective coefficients, such as firstly assume that there is functional relationship between the two coupling coefficients, and then combine the functional relationship and coupling corrective relationship together to work out the coefficients. Through the above derivation process on SP3 equation, a brief analysis of the SPN theory is considered. The second-order form such as Eq. (2.9) is useful because it makes the SPN equations look
925
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927
like a set of coupled diffusion equations. It can be deduced into a matrix form (Ciolini et al., 2002) and shows as follow,
A
rt
!
!
!
r2 / þrt / ¼ Q
ð3:24Þ
Z mknu1 ðsinhÞ ¼ 1=N1
where, !
/ ¼ ½/0 ; /2 ; ; /ðNþ1Þ=2 and;
!
Q ¼ ½Q ; 0; ; 0:
ð3:25Þ
where, A is an (N + 1)/2 order square, whose non-vanishing elements are, nðn1Þ An;n1 ¼ ð2nþ1Þð2n1Þ ;
2
2n þ2n1 An;n ¼ ð2nþ3Þð2n1Þ ;
ðnþ1Þðnþ2Þ An;nþ1 ¼ ð2nþ1Þð2nþ3Þ
n ¼ 1; 2; ; ðN þ 1Þ=2
1
Z mknui ðcoshÞ ¼ 1=Ni
1
sinhðakn;gu tÞpn1 ðtÞdt
1
1
coshðakn;gu tÞpni ðtÞdt
Ni ¼ 2=ð2i þ 1Þ i ¼ 0; 1; 2 Solving the simultaneous equations of the 0th-order moment, the 1st-order moment, the 2nd-order moment and the boundary continuity condition, the net neutron currents of the interface are obtained. Dk
J kn;guþ ¼ Dun;g=2 ðakn;gu1 þ 3akn;gu2 þ Hkn;gu akn;gu3 þ Gkn;gu akn;gu4 Þ k
ð3:26Þ
Dk
J kn;gu ¼ Dun;g=2 ðakn;gu1 3akn;gu2 þ Hkn;gu akn;gu3 Gkn;gu akn;gu4 Þ
ð3:31Þ
k
The neutron balance equation of the node can be obtained by integrating the SPN equations in node, ! ! B ! J þrt / ¼ Q Du
ð3:27Þ
And, (n = 0,2)
Hkn;gu ¼
akn;gu coshðakn;gu Þmknu1 ðsinhÞ sinhðakn;gu Þmknu1 ðsinhÞ ak
!
where, J ¼ ½J 0;uþ J0;u ; J2;uþ J2;u ; ; J ðNþ1Þ=2;uþ J ðNþ1Þ=2;u , and matrix B can be simply calculated from the expression of A. Then we substitute the new coupling corrective relationship Eq. (2.16) into the neutron balance equation, the same matrix equations as Eq. (3.8) is get and the subsequent derivation process is the same as the above SP3 equation. Therefore, the main point after derivation is also can be found that, for the numerical stability, to minimize the coupling correction factor as possible and the difference on the net neutron flux should be distributed to each coupling correction factor. 3.3. The Semi-analytical nodal method The code NLSP3 has selected the semi-analytical nodal method (Zimin et al., 1998) to calculate the net currents on the surface. The semi-analytical nodal method has higher precision and efficiency in the calculation of multi-group equations. In this paper, the process of the semi-analytical nodal method in SP3 is briefly introduced. The transverse integral equations of SP3 is marked as Eq. (3.28),
D0 @u@ 2 /0u ðuÞ þ R0 /0u ðuÞ 2R0 /2u ðuÞ ¼ S0u ðuÞ L0u ðuÞ D2 @u@ 2 /2u ðuÞ þ R2 /2u ðuÞ 25 R0 /0u ðuÞ ¼ 25 S0u ðuÞ L2u ðuÞ ð3:28Þ
Gkn;gu ¼ coshðn;gu ak
sinhðakn;gu Þ3mknu2 ðcoshÞ
With the Eqs. (3.4) and (3.31), the coupling corrective coefficients are obtained. After the global calculation of CMFD formula, the angular flux and power distribution are worked out. However, the 2nd-order flux of SP3 is hundreds of times smaller than the 0thorder flux, according to Eq. (3.23) and for numerical stability, it is necessary to avoid the 2nd-order flux being a divisor, otherwise the coefficients will be too large. There are many ways to deal with it. As an example, the coupling corrective relationship shows on k;NOD Dk;NOD 4;guþ for 2nd-order flux, while D2;guþ can be set to be zero, and that is what exactly the code NLSP3 do. Thus, the numerical instability problem of nonlinear iteration is solved.
4. Numerical results 4.1. The 2D-IAEA transport benchmark This benchmark (Stepanek and Auerbach, 1983) simulates a mono-energetic 5-zones light water reactor, surrounded by vacuum boundary conditions. Geometric structure shown in Fig. 2, and regional material parameters in Table 1. In this paper, we use NLSP3 program to calculate the two different meshing of the problem, the results shown in Table 2. The results which is
With the Polynomial expansion form,
/kn;gu ðuÞ ¼ /kn;g þ
4 X 2u akn;gui pni Duk i¼1
ð3:29Þ
The 4th-order polynomial of the semi-analytic nodal method shows as follows:
Pn0 ðuÞ ¼ 1 u Pn1 ðuÞ ¼ Du=2 ¼t
Pn2 ðuÞ ¼ 12 ð3t2 1Þ Pn3 ðuÞ ¼
sinhðakn;gu tÞmknu1 ðsinhÞpn1 ðtÞ sinhðakn;gu Þmknu1 ðsinhÞ
ð3:30Þ
coshðakn;gu tÞmknu0 ðcoshÞpn0 ðtÞmknu2 ðcoshÞpn2 ðtÞ
Pn4 ðuÞ ¼ u 2 a2u ; a2u
aknu
coshðakn;gu Þmknu0 ðcoshÞmknu2 ðcoshÞ
vffiffiffiffiffiffiffiffiffi u k uRn;g Duk ¼t k Dn;g 2
ð3:32Þ
k k n;gu Þmnu0 ðcoshÞmnu2 ðcoshÞ
ðn ¼ 0; 2Þ Fig. 2. Structure of 2D-IAEA Benchmark.
926
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927
Table 1 Material parameters of 2D-IAEA Benchmark. zones
Rt =cm1
Rs =cm1
mRf =cm1
1 2 3 4 5
0.60 0.48 0.70 0.65 0.90
0.53 0.20 0.66 0.50 0.89
0.079 0.0 0.043 0.0 0.0
compared to the results from SURCU (Stepanek and Auerbach, 1983), DNTR (Li et al., 2010) and EFEN (Li et al., 2011) show that NLSP3 has the same accuracy as other transport programs.
4.2. The 3D-TAKEDA benchmark The selected model is the second one in 3D-TAKEDA benchmarks (Takeda and Ikeda, 1991). A small fast breeder reactor is
Table 2 Results of 2D-IAEA Benchmark. Program
mesh
Flux in zone 1
Flux in zone 2
Flux in zone 3
Flux in zone 4
Flux in zone 5
keff
SURCU DNTR(S4) EFEN NLSP3 NLSP3
– 1148D 72 64h 16 16h 96 86h
0.01686 0.01686 0.01686 0.01686 0.01686
0.000125 0.000125 0.000125 0.000124 0.000124
0.000041 0.000035 0.000040 0.000027 0.000036
0.000295 0.000295 0.000297 0.000294 0.000294
0.000791 0.000791 0.000795 0.000715 0.000789
1.0083 1.0085 1.0078 1.0088 1.0089
Note: D Triangular meshing; hGraphical meshing
Fig. 3. Structure of 3D-TAKEDA Benchmark.
Table 3 Results of 3D-TAKEDA Benchmark. Method keff and errors (105)
Exact Monte-Carlo
NSPn 5 5 5(cm)
NLSP3 5 5 5(cm)
0.9732
363
0.9734
20
group
Ref flux
Errors (%)
flux
Errors (%)
fuel zone
1 2 3 4
4.2814E05 2.4081E04 1.6411E04 6.2247E06
0.065 0.014 0.240 0.678
4.2680E05 2.4043E04 1.6388E04 6.2058E06
0.3130 0.1568 0.1378 0.3043
axial propagation zone
1 2 3 4
5.1850E06 4.6912E05 4.6978E05 3.7736E06
2.065 1.272 0.888 0.229
5.2270E06 4.6962E05 4.6944E05 3.7366E06
0.8094 0.1072 0.0715 0.9802
radial propagation zone
1 2 3 4
3.3252E06 3.0893E05 3.2834E05 2.0473E06
2.351 1.110 0.879 0.698
3.3568E06 3.0851E05 3.2765E05 2.0423E06
0.9507 0.1350 0.2095 0.2451
sodium filling zone
1 2 3 4
2.5344E05 1.6658E04 1.2648E04 6.9840E06
3.398 0.361 0.360 0.917
2.6053E05 1.6647E04 1.2658E04 6.9970E06
2.7975 0.0637 0.0809 0.1865
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927
simulated. There are four material zones, namely fuel zone, axial propagation zone, radial propagation zone and sodium filling zone while the control rod has a full lift. The geometrical structure is shown in Fig. 3. The cross-section parameters are described in references and the reference solution is obtained by the exact Monte Carlo method. We used the NLSP3 program to calculate the benchmark under the control rod full-lift condition (5 5 5 cm). The effective multiplicative coefficient and relative power distribution were listed in Table 3. It can be seen from the calculation results that the errors of the effective multiplication coefficient is very small and the calculated results are in good agreement with the reference solution except the flux errors at the sodium-filled zone. The results also show that the NLSP3 program based on nonlinear iteration has a better accuracy than the NSPn (Li et al., 2010) based on the exponential function method in the case of the same meshing for 3D-TAKEDA. 5. Conclusion Due to the conventional coupling corrective relationship’s defect, the numerical instability problem is found during the coupling correction of the higher-order flux when adopting the conventional nonlinear iterative method to solve the SPN equation. Even though in other nodal approaches for SPN equation, the numerical instability problem is still existing when CMFD formula is used to speed up the calculation. To solve the numerical instability problem, a new nonlinear iterative method for SPN theory is put forward by importing a new coupling corrective relationship which focuses on the coupling at the angle. Compared to the conventional nonlinear iterative method, the new method has better application prospect. Except for the calculation at the vacuum boundary, the coupling corrective coefficients have many ways to calculate, which brings in more margin to deal with the numerical instability. In this paper, some inequalities for how to confirm the exact values of the coupling corrective coefficients are given through numerical derivation about the norm of the iterative matrix. The inequality can give the guide and narrow down the range to work out the coupling corrective coefficients. Under the premise of
927
satisfying the inequality, the functional relationships between each coupling corrective coefficients can be selected freely according to the different situations. With that constructor function, the coefficients will be solved by the coupling corrective relationship, and the convergent results are obtained by substituting the coefficients into the CMFD equation. As an example, the code named NLSP3 for solving the SP3 equations with the new nonlinear iterative method was developed, and the numerical results demonstrate the feasibility of the new method. References Berkert, C., Grundmann, U., 2008. Development and verification of a nodal approach for solving the multigroup SP3 equations. Ann. Nucl. Energy 35, 75–86. Brantley, P.S., Larsen, E.W., 2000. The simplified P3 approximation. Nucl. Sci. Eng. 134, 1–21. Chao, Y.A., 2015. A new SPN theory formulation with self-consistent physical assumptions on angular flux. Ann. Nucl. Energy 87, 137–144. Ciolini, R., Coppa, G.G.M., Montagnini, B., Ravetto, P., 2002. Simplified PN and AN methods in neutron transport. Prog. Nucl. Energy 40 (2), 237–264. Fu, X.D., 2001. A numerical stabilization method for nonlinear analytic nodal method. Chin. J. Comput. Phys. 18 (2), 111–114. Gelbard, E.M., 1960. ‘‘Application of spherical harmonics method to reactor problems”. Technical Report WAPD-BT-20, Bettis Atomic Power Laboratory. Li, Y.Z., Wu, H.C., Cao, l.z., 2010. Nodal SP3 method for solving neutron transport equation. Nucl. Power Eng. 31, 73–78. Li, Y.Z., Wu, H.C., Cao, l.z., 2011. ‘‘A polygonal nodal-sp3 method for whole core pinby-pin calculation” // M&C2011. Rio de Janeiro. Brazil. Smith, K.S., 1983. Nodal method storage reduction by nonlinear iteration. Trans. Am. Nucl. Soc 44, 265. Stepanek, J., Auerbach, T., 1983. ‘‘Calculation of four thermal reactor benchmark problems in X-Y geometry”. EPRI NP-2855. Takeda T., Ikeda H., 1991. ‘‘3-D Neutron Transport Benchmarks[R]. OECD/NEA Committee on Reactor Physics” (NEACRP). NEACRP-L-330. Japan:OSAKA University. Yu, L.L., 2014. In-Depth Investigation and Further Development of Homogenization Method and Discontinuity Factor Theory[D]. Shanghai Jiaotong University, Shanghai. Yu, L.L. et al., 2014. The calculation method for SP3 discontinuity factor and its application. Ann. Nucl. Energy 69, 14–24. Zhang, Y.H., 2013. Iterative performance of various formulations of the SPN equations. J. Comput. Phys. 252, 558–572. Zimin, V.G., Ninokata, H., Pogosbekyan, L.R., 1998. Polynomial and semi-analytic nodal method for nonlinear iteration procedure[C]. Proc. Int. Conf. Phys, Nucl. Sci. Tech. (PHYSOR-98), Long Island, NY, October 5–8, pp. 994–1002.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
A new nonlinear iterative method for SPN theory Pan Qingquan a,⇑, Lu Haoliang b, Li Dongsheng b, Wang Kan a a b
Tsinghua University, Department of Engineering Physics, Beijing 100084, China China Nuclear Power Technology Research Institute, Reactor Engineering and Safety Research Center, Shenzhen, Guangdong 518031, China
a r t i c l e
i n f o
Article history: Received 18 May 2017 Received in revised form 13 July 2017 Accepted 25 July 2017
Keywords: SPN theory Nonlinear iteration Numerical instability
a b s t r a c t The SPN theory shows a great potential for solving the neutron transport equation in the Next Generation Reactor Physics Calculation. However, the numerical instability problem will be encountered in some cases when solving the SPN equation with the conventional nonlinear iterative method. In order to solve the problem, a new coupling corrective relationship which focuses on the angular coupling is put forward, and the new nonlinear iterative method has more margin to adjust the coupling corrective coefficients with that new relationship. As an example, the code named NLSP3 for solving the SP3 equation with the new nonlinear iterative method is developed, and the numerical results demonstrate the feasibility of the new method. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction The conventional nonlinear iterative method (Smith, 1983) is based on the nodal method and nonlinear iterative strategy, which is characterized by the use of the simplified CMFD formula of the neutron diffusion equation containing the nodal coupling corrective coefficients in global calculation. The nodal coupling corrective coefficients determined by solving the two-nodal problem at the interfaces of two adjacent nodes will be updated during the nonlinear iterative process. The method is mainly used to deal with the diffusion equation, however, due to it’s high efficiency, the method is expected to be used in the transport equation, especially in the SPN equation. The SPN theory (Chao, 2015) is the simplified spherical harmonics method for solving the neutron transport equation. It was widely studied (Yu et al., 2014; Yu, 2014) and it shows a great potential in the Next Generation Reactor Physics Calculation because of its better accuracy than diffusion theory. Adopting the conventional nonlinear iterative method to solve the SPN equation, the numerical instability problem is found during the coupling correction of the high-order flux since the high-order coupling corrective coefficients are too large (Yu et al., 2014; Yu, 2014). Even though in other nodal approaches (Berkert and Grundmann, 2008), the numerical instability problem is still existing when CMFD formula is used to speed up the calculation on the SPN equation. So, it is very important to put ⇑ Corresponding author at: Department of Engineering Physics, Tsinghua University, Haidian District, Beijing 100084, China. E-mail addresses: [email protected] (Q. Pan), [email protected] (H. Lu), [email protected] (D. Li), [email protected] (K. Wang). http://dx.doi.org/10.1016/j.anucene.2017.07.030 0306-4549/Ó 2017 Elsevier Ltd. All rights reserved.
forward a new nonlinear iterative method for the SPN theory. In this paper, a new nonlinear iterative method is discussed with the using of a new coupling corrective relationship. The new coupling corrective relationship focuses on the coupling at the angle, and seems a bit more complicated than the conventional one. Except for the calculation at the vacuum boundary, the coupling corrective coefficients have many ways to calculate, which brings in margin to deal with the numerical instability under adequate mathematical derivation or justification. The code named NLSP3 is developed to achieve the new method within the SP3 theory, and the numerical results demonstrate the feasibility of the new method. In this paper, a review summary of Chao’s derivation for the SPN theory is given in Section 2, where the necessary derivations and equations are outlined. And the new coupling corrective relationship for SPN theory is generally introduced. Section 3 presents the special case for SP3 with the new method, and some advices are given after the discussion about the norm of the iterative matrix. Section 4 provides the numerical results and comparisons. The 2D-IAEA transport benchmark and the 3D-TAKEDA benchmark are calculated by the code NLSP3, and the feasibility of the new method is proved. Section 5 makes a conclusion about the work. 2. The new method for SPN theory
2.1. The SPN theory The SPN equation was firstly proposed without rigorous mathematical derivation (Gelbard, 1960). Later, the asymptotic theory
921
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927
approximation and variational method (Brantley and Larsen, 2000) were introduced to deduce a more rigorous equation in SP3. Those year, Professor Chao has done a lot of jobs (Chao, 2015) about SPN theory. In this paper, Chao’s jobs are considered as a reference and a brief introduction about his derivational process is given. Under the case of mono-energetic isotropic scattering, the neutron transport equation shows as follow,
X rwðX; rÞ þ Rt wðX; rÞ ¼
QðrÞ 4p
ð2:1Þ
And the neutron flux are divided into the even and odd parity angular flux,
wðX; rÞ ¼ wE ðX; rÞ þ wO ðX; rÞ
ð2:2Þ
With the expansion in orthogonal spherical harmonics and the assumption of cylindrically symmetric for the azimuth angular dis/n ðrÞ tribution of the nth flux moments, and Xn ¼ jr r/n ðrÞj ; we have,
2n þ 1 Pn ðX Xn Þ/n 4p
R
wE ðX; rÞ ¼
n¼ev en
wo ðX; rÞ ¼
1
Rt
R
n¼ev en
ð2:3Þ
2n þ 1 ½Pn ðX Xn ÞX Xn ½Xn r/n ðrÞ 4p
J n ¼
Z
b Þð n b XÞwE ðX; rÞdX bÞ Pn ðX n Pn ðX n bn X>0 bn X>0 b XÞwo ðX; rÞdX ðn
ð2:11Þ
The final expression for the nth net current as follows.
b Xn2 Þ nðn 1Þ Pn1 ð n b r/n2 ðrÞ n b Xn2 ð2n þ 1Þð2n 1ÞRt n b Xnþ2 Þ ðn þ 1Þðn þ 2Þ Pnþ1 ð n b r/nþ2 ðrÞ n b ð2n þ 1Þð2n þ 3ÞRt n Xnþ2 " # b Xn Þ b Xn Þ 1 ðn þ 1Þ2 Pnþ1 ð n n2 Pn1 ð n b þ n b Xn b Xn ð2n þ 1ÞRt 2n þ 3 2n 1 n n
Jn ¼
r/n ðrÞ ð2:12Þ With the definition of the incoming and outgoing partial currents, such as Eqs. (2.13) and (2.14), the whole SPN theory was deduced.
1 J n ¼ Un J n 2
ð2:4Þ
In order to derive the net current, we integrate the product of X and odd parity angular flux over the solid angle space,
Z
Un ¼
Z J þn þ J n b Þð n b XÞwE ðX; rÞdX Pn ðX n ¼ 2 bn X>0
ð2:13Þ
ð2:14Þ
Z
J¼
2.2. The new coupling corrective relationship
X wo ðX; rÞdX
¼
2k þ 1
1
R Rt k¼ev en 4p
Z
X Pk ðX Xk Þ½X r/k ðrÞdX
ð2:5Þ
With the use of the Legendre polynomial recursion relations, we have
1 ðk 1Þk d0;k þ d0;k2 r/k ðrÞ 2k 1 Rt k¼ev en 3 1 ¼ ½r/0 ðrÞ þ 2r/2 ðrÞ 3Rt
J¼
1
R
ð2:6Þ
where, dm;n is the Dirichlet function. In order to derive the differential equations of nth flux, substitute Eq. (2.3) into (2.1) and integrate it over the angular space,
2k þ 1 k¼ev en 4p
R
Z
1 Pn ðX Xn ÞPk ðX Xk Þ X r X r Rt /k ðrÞ dX
Rt
ð2:7Þ
¼ Qdn0 Making use of Eqs. (2.3) and (2.4), we have,
Xr
1 X r/k ðrÞ ¼ X r ðX Xk ÞðXk r/k ðrÞÞ Rt Rt r/k ðrÞ ¼ ðX Xk ÞX r Xk 1
Rt
ð2:8Þ
kþ1 k kþ1 k J kn;guþ ¼ Dk;FDM n;guþ f n;gu /n;g f n;guþ /n;g kþ1 k kþ1 k Dk;NOD n;guþ f n;gu /n;g þ f n;guþ /n;g
After some variation, SPN equation can be get: 2 2 nðn1Þ 2n þ2n1 ð2nþ1Þð2n1Þ Rt r /n2 ðrÞ ð2nþ3Þð2n1ÞRt r /n ðrÞ 2
2 ðnþ1Þðnþ2Þ ð2nþ1Þð2nþ3Þ Rt r /nþ2 ðrÞ þ Rt /n ðrÞ ¼ Q dn0
ð2:9Þ
The boundary conditions are provided by the relation between the incoming and outgoing partial current. And the nth partial curb at a surface is defined as, rent with the normal vector n b X½wE ðX; rÞ wo ðX; rÞ jbn X>0 ðX; rÞ ¼ n
The nonlinear iterative method is using the results with higherorder precision to correct the Coarse Mesh Finite Difference formula (CMFD) which has a lower-order precision. And the coupling corrective relationship characterizes the difference of the net neutron current calculated between the higher-order nodal method and CMFD. Therefore, the correction process is actually to import a perturbation which is exactly the coupling corrective coefficients to the CMFD formula (Fu, 2001). In that way, the solutions of the CMFD equation will achieve the accuracy of the higher-order node method. When the perturbation is too large, it will make the iterative matrix be a morbid matrix and the nonlinear iteration not converges, which exact is the numerical instability problem. Fig. 1 is the nonlinear iterative calculation strategy, consisting of source iteration and SANM coupling corrective coefficients updating, where the source iteration is divided into internal iteration and outer iteration. In the case of internal iteration, the source term is regarded as a constant, so the iteration is actually solving a linear system of equations. At that situation, the numerical instability problem occurs in the internal iteration, so it is necessary to analyze the iterative matrix of the internal iteration. The conventional coupling corrective relationship is given as Eq. (2.15), it calculates the net current only under the nth flux itself, without the coupling relationship between the different order flux.
ð2:10Þ
Make use of nth even order Legendre polynomial, and integrate Eq. (2.10) over the angular space.
ð2:15Þ
where, ‘‘k” is the number of the node, ‘‘g” is the number of the group. Normally, ‘‘n’ is set to be an odd number. We regard the nth-order flux as the different description about the angle distribution, so, the angular flux is the superposition of all the nth order flux. In order to solve the SPN equation with the current multigroup diffusion codes, The SPN equation can be written as (N + 1)/2 coupled differential equations on the each nth even order flux after some mathematical derivation and substitution, such as the composite form (Zhang, 2013). Adopting the conventional coupling corrective relationship to solve the SPN equation, the numerical
922
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927
2 2 Dk2;g r2 /k2;g ðrÞ þ Rkr2;g /k0;g ðrÞ Rkr0;g /k0;g ðrÞ ¼ Sk0;g 5 5
ð3:1bÞ
where,
Dk0;g ¼
1 3Rkt;g
; Dk2;g ¼
9 5
9 35Rkt;g
; Rkr0;g ¼ Rkt;g Rks;gg0 ;
4 5
Rkr2;g ¼ Rkt;g Rks;gg0
ð3:1cÞ
The Source term expression,
h i G Sk0;g ðrÞ ¼ R Rk0;g0 g /k0;g ðrÞ 2/k2;g ðrÞ g 0¼1 g 0 –g
þ Fig. 1. Nonlinear iterative strategy.
k
kþ1 k Dk;NOD 0;guþ ðf 0;gu /0;g þ f 0;guþ /0;g Þ kþ1 kþ1 Dk;FDM 2;guþ ðf 2;gu /2;g
kþ1 kþ1 Dk;NOD 2;guþ ðf 2;gu /2;g
þ f 2;guþ /k2;g Þ
ð3:2aÞ
J 2 ðrÞ ¼
21 M 1* 3 M / ðrÞ n J M / 2 ðrÞ 80 2 2 80 0
ð3:2bÞ
k
X
kþ1 k Dk;FDM n;guþ ðf n;gu /n;g f n;guþ /n;g Þ
u¼x;y;z
kþ1 k Dk;NOD n;guþ ðf n;gu /n;g þ f n;guþ /n;g Þ
kþ1
k
kþ1
k
1 Duk
ðJ k0;guþ J k0;gu Þ þ Rkr0;g Uk0;g 2Rkr0;g Uk2;g ¼ Sk0;g
1 Duk
ðJ k2;guþ J k2;gu Þ þ Rkr2;g Uk0;g 25 Rkr0;g Uk2;g ¼ 25 Sk0;g
kþ1 k Dk;FDM ðNþ1Þ=2;guþ ðf ðNþ1Þ=2;gu /ðNþ1Þ=2;g f ðNþ1Þ=2;guþ /ðNþ1Þ=2;g Þ
J kn;guþ ¼ Dkn;g
kþ1 k Dk;NOD ðNþ1Þ=2;guþ ðf ðNþ1Þ=2;gu /ðNþ1Þ=2;g þ f ðNþ1Þ=2;guþ /ðNþ1Þ=2;g Þ
n ¼ 0; 2; . . . ; ðN þ 1Þ=2 ð2:16Þ
With the new coupling corrective relationship, the CMFD formula seems a bet more complicated than the conventional one. For simplicity we will take the special case of SP3 as an example, and advices for how to calculate the coupling corrective coefficients will be given through mathematic derivation in next section.
@/kn;g ðuÞ @u
ð3:3bÞ
h i h i G 1 k G Sk0;g ¼ R Rk0;g0 g Uk0;g 2Uk2;g þ vg R mRkf;g0 Uk0;g 2Uk2;g 0 ¼1 g 0 ¼1 g K eff 0
ð3:3cÞ
g –g
Ukn;g ¼
1 Vk
Z V
/kn;g ðrÞdV
ð3:3dÞ
According to Eq. (2.16), we have the specific expression of the new coupling corrective relationship for SP3 equations, kþ1
k
kþ1
k
k;NOD kþ1 kþ1 k k J k0;guþ ¼ Dk;FDM 1;guþ ðf 0;gu /0;g f 0;guþ /0;g Þ D1;guþ ðf 0;gu /0;g þ f 0;guþ /0;g Þ
3. The special case for SP3
kþ1
k
kþ1
k
k;NOD kþ1 kþ1 k k Dk;FDM 3;guþ ðf 2;gu /2;g f 2;guþ /2;g Þ D3;guþ ðf 2;gu /2;g þ f 2;guþ /2;g Þ kþ1
k
kþ1
k
k;NOD kþ1 kþ1 k k J k2;guþ ¼ Dk;FDM 2;guþ ðf 2;gu /2;g f 2;guþ /2;g Þ D2;guþ ðf 2;gu /2;g þ f 2;guþ /2;g Þ
3.1. The SP3 equations
kþ1
k
kþ1
k
k;NOD kþ1 kþ1 k k Dk;FDM 4;guþ ðf 0;gu /0;g f 0;guþ /0;g Þ D4;guþ ðf 0;gu /0;g þ f 0;guþ /0;g Þ
The SP3 equation can be got when N is set to be 3 in the SPN equation. After deformation, a similar form to the neutron diffusion equation can be obtained and the same nodal method can be used to solve the SP3 equation, such as the nonlinear iterative method. Assuming that the neutron scattering and the source term are isotropic, and the high order scattering between the nth groups is zero, then the specific form of the SP3 equation and the boundary conditions show as follow,
Dk0;g r2 /k0;g ðrÞ þ Rkr0;g /k0;g ðrÞ 2Rkr0;g /k2;g ðrÞ ¼ Sk0;g
ð3:3aÞ
Where,
for
1 M 1* 3 M / ðrÞ n J M / 0 ðrÞ 4 0 2 16 2
u¼x;y;z
kþ1
ð3:1dÞ
g ¼1
J 0 ðrÞ ¼
X
k
k
i
The neutron balance equation of the node can be obtained by integrating the SP3 equations in node ‘‘k”, and u as a reference direction as follow,
k f 2;guþ /k2;g Þ
kþ1
h
3.2. The new method in SP3 equations
k
kþ1 k J kn;guþ ¼ Dk;FDM 0;guþ ðf 0;gu /0;g f 0;guþ /0;g Þ kþ1
G
vkg R mRkf;g0 /k0;g0 ðrÞ 2/k2;g0 ðrÞ 0
At the boundary, we have the neutron flux and the current under the Marshak boundary condition as follow, where, the superscript ‘‘M” represents the amounts at the node boundary.
instability problem is found during the coupling correction of the higher-order flux, because the higher-order flux is very small which leads the higher-order coupling corrective coefficients to be too large. So, we redefine the coupling corrective relationship and their coupling at angles is considered. So, the new coupling corrective relationship shows as follow, kþ1
1 K eff
ð3:1aÞ
ð3:4Þ Plugging Eq. (3.4) into Eq. (3.3), then we have the CMFD formula as Eq. (3.5), kþ1 k1 k1 kþ1 k kþ1 k1 k1 k k ak0;gx f kþ1 0;gx /0;gx þ b0;gx f 2;gx /2;gx þ c0;gxþ f 0;gxþ /0;gx þ q0;gx f 2;gxþ /2;gxþ þ kþ1 k1 k1 kþ1 k kþ1 k1 k1 k k ak0;gy f kþ1 0;gy /0;gy þ b0;gy f 2;gy /2;gy þ c0;gyþ f 0;gyþ /0;gy þ q0;gy f 2;gyþ /2;gyþ þ kþ1 k1 k1 kþ1 k kþ1 k1 k1 k k ak0;gz f kþ1 0;gz /0;gz þ b0;gz f 2;gz /2;gz þ c0;gzþ f 0;gzþ /0;gz þ q0;gz f 2;gzþ /2;gzþ þ
C k1;gz /k0;g þ C k3;gz /k2;g ¼ Sk0g ð3:5aÞ
923
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927 kþ1 kþ1 k1 k1 k kþ1 k1 k1 k k ak2;gx f 2;gx /kþ1 2;gx þ b2;gx f 0;gx /0;gx þ c2;gxþ f 2;gxþ /2;gx þ q2;gx f 0;gxþ /0;gxþ þ kþ1 kþ1 k 2;gy f 2;gy /2;gy
a
kþ1 kþ1 k 2;gz f 2;gz /2;gz
a
þ
kþ1 kþ1 bk2;gy f 0;gy /0;gy
þ
kþ1 kþ1 bk2;gz f 0;gz /0;gz
k1 k1 k 2;gyþ f 2;gyþ /2;gy
þc
k1 k1 k 2;gzþ f 2;gzþ /2;gz
þc
k1 k1 k 2;gy f 0;gyþ /0;gyþ þ
þq
k1
k1 þ qk2;gz f 0;gzþ /0;gzþ þ
C k2;gz /k2;g þ C k4;gz /k0;g ¼ 25 Sk0g ð3:5bÞ
"
where,
C k1;gz ¼ Rk0;g þ Rk0;gx þ Rk0;gy þ Rk0;gz
ð3:6aÞ
C k3;gz ¼ 2Rk0;g þ Rk2;gx þ Rk2;gy þ Rk2;gz
ð3:6bÞ
C k2;gz ¼ Rk2;g þ Lk2;gx þ Lk2;gy þ Lk2;gz
ð3:6cÞ
C k4;gz
2 ¼ Rk0;g þ Lk0;gx þ Lk0;gy þ Lk0;gz 5
k 0;gu
1 k;FDM ¼ D þ Dk;NOD 1;guþ Duk 1;guþ
k 2;gu
1 k;FDM ¼ D þ Dk;NOD 2;guþ Duk 2;guþ
a a
A2
2
#
½U ¼ 4
ð3:6eÞ
ð3:6fÞ
ð3:6gÞ
bk2;gu
1 k;FDM ¼ D4;guþ þ Dk;NOD 4;guþ Duk
ð3:6hÞ
k 0;gu
c
1 k;FDM ¼ D1;gu Dk;NOD 1;gu Duk
ck2;gu
1 k;FDM ¼ D2;gu Dk;NOD 2;gu Duk
ð3:6jÞ
qk0;gu ¼
1 k;FDM D3;gu Dk;NOD 3;gu Duk
ð3:6kÞ
qk2;gu ¼
1 k;FDM D4;gu Dk;NOD 4;gu Duk
ð3:6lÞ
3
As þ kF
eff
25 ðAs þ kF Þ
"
5½U ¼
eff
S
# ð3:8aÞ
25 S
where,
U ¼ colðU1 ; U2 ; ; Ug ; ; UG ; Þ Ug ¼ colðU1g ; U2g ; ; Ukg ; ; UNg Þ
ð3:8bÞ
Ukg ¼ colðUk0;g ; Uk2;g Þ 2
1 k;FDM ¼ D þ Dk;NOD 3;guþ Duk 3;guþ
AG1 s
3
0 6 12 6A 6 s 6 13 As ¼ 6 6 As 6 6 4
A21 s
A31 s
0
A32 s
A23 s
0
7 7 AG2 s 7 7 G3 7 As 7 7 7 5
A1G s
A2G s
A3G s
ð3:8cÞ
0
Asg g ¼ diagðR1g g ; 2R1g g ; R2g g ; 2R2g g ; ; Rkg g ; 2Rkg g ; ; RNg g ; 2RNg g Þ 2 ð3:6iÞ
And,
Rk0;gu ¼
1 k;FDM k;FDM k;NOD D1;guþ Dk;NOD 1;guþ þ D1;gu þ D1;gu Duk
ð3:7aÞ
Rk2;gu ¼
1 k;FDM k;FDM k;NOD D3;guþ Dk;NOD 3;guþ þ D3;gu þ D3;gu Duk
ð3:7bÞ
1 k;FDM k;FDM k;NOD ¼ D4;guþ Dk;NOD 4;guþ þ D4;gu þ D4;gu Duk
ð3:7cÞ
1 k;FDM k;FDM k;NOD D2;guþ Dk;NOD 2;guþ þ D2;gu þ D2;gu Duk
ð3:7dÞ
Lk2;gu ¼
A0
ð3:6dÞ
bk0;gu
Lk0;gu
The nonlinear iterative method chooses the source iteration to solve the CMFD equation, in which the numerical instability occurs during the internal iteration. Within the internal iterative process, the calculation of the sources bases on the last results of the outer iteration, that is, the sources are treated as constants. In the case of ‘‘N” nodes and ‘‘G” energy groups, Eq. (3.5) can be rewritten in matrix form, as follow,
In order to solve the CMFD formula, we need the values of the coupling corrective coefficients which should be solved by Eq. (3.4). However, there are two coupling coefficients in each relationship while only one net current can be used. So, except for the calculation at the vacuum boundary, the coupling corrective coefficients have many ways to calculate, which brings in the margin to deal with the numerical instability under adequate mathematical derivation or justification. Then, some advices for how to confirm the exact coefficients are given through numerical derivation.
F 11
F 21
F G1
ð3:8dÞ 3
6 6 F 12 F¼6 6 4
7 F G2 7 7 7 5
F 1G
F 2G
F GG
F 22
ð3:9Þ
F g g ¼ diagðvg mR1fg ; 2vg mR1fg ; vg mR2fg ; 2vg mR2fg ; ; vg mRkfg ; 2vg mRkfg ; ; vg mRNfg ; 2vg mRNfg Þ And,
A0 ¼ M 0 þ N 0
ð3:10aÞ
A2 ¼ M 2 þ N 2
ð3:10bÞ
M0 ¼ diagðM 01 ; M02 ; M 03 ; ; M0g ; ; M 0G Þ
ð3:11aÞ
N0 ¼ diagðN01 ; N02 ; N03 ; ; N0g ; ; N0G Þ
ð3:11bÞ
M2 ¼ diagðM 21 ; M22 ; M 23 ; ; M2g ; ; M 2G Þ
ð3:11cÞ
N2 ¼ diagðN21 ; N22 ; N23 ; ; N2g ; ; N2G Þ
ð3:11dÞ 3
2 1 1 6 C 1;g C 3;g 6 6 C 21;g 2 6 6 C 3;g 6 0 Mg ¼ 6 6 6 6 6 C k1;g C k3;g 4
7 7 7 7 7 7 7 7 7 7 7 5
N N C 1;g C 3;g
ð3:12aÞ
924
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927
3
2
a10;gx
6 b10;gx a b10;gy a 6 2 6 m0;gx q20;gx b10;gz 2 2 2 2 a a 6 b b 0;gx 0;gy 0 0;gx 0;gy Ng ¼ 6 a20;gz k 6 a0;gx bk0;gx ak0;gy bk0;gy ak0;gz bk0;gz 6 mk0;gz qk0;gz 6 k k k k 4 m0;gy q0;gy m0;gz q0;gz mN0;gz qN0;gz mN0;gy qN0;gy mN0;gz qN0;gz 1 0;gy
b20;gz
1 0;gz
3
2 1 1 6 C 4;g C 2;g 6 C 24;g 2 6 6 C 2;g 2 Mg ¼ 6 6 6 6 4 C k4;g C k2;g
7 7 7 7 7 7 7 7 5
ð3:12bÞ
N N C 4;g C 2;g
3
2 b12;gx
6 a22;gz a12;gx b12;gy a12;gy b12;gz a12;gz 6 2 6 q2;gx m22;gx 2 2 2 2 2 6 b2;gx a2;gx b2;gy a2;gy b2;gz N 2g ¼ 6 6 bk2;gx ak2;gx bk2;gy ak2;gy bk2;gz ak2;gz 6 qk2;gz mk2;gz 6 k k 4 q2;gy mk2;gy q2;gz mk2;gz qN2;gz mN2;gz qN2;gy mN2;gy qN2;gz mN2;gz
In the internal iterative process, the Eq. (3.8) can be rewritten as Eq. (3.13),
"
#
"
A0 ½U ¼ A2
M0 þ N 0 2
M þN
#
"
½U ¼
2
#
M0 M
2
"
½U þ
N0
"
#
N2
½U ¼
S
#
25 S ð3:13Þ
The iterative formula for calculating CMFD equation show as follow,
" ½UðKþ1Þ ¼
M0
#1 "
S 25 S
M2
#
"
M0
#1 "
M2
N0 N2
# ½UðKÞ
ð3:14Þ
M¼
M0
#
"
;N ¼
M2
N0 N2
# then E ¼
"
M0
#1 "
M2
N0 N2
# ¼ M 1 N
ð3:15Þ
kEk ¼ kM
Nk < 1
ð3:16Þ
And we know the relationship for all the matrix norm as Eq. (3.17),
kT 1 k kTk P kT 1 Tk ¼ 1
ð3:17Þ
ð3:18Þ
So, We can get,
kMk P
kNk P
1 kN1 k 1 kN1 k
ð3:21Þ
Because the equivalence of matrix norm, we can choose each one from the three kinds of matrix norm. So, there are many ways to conduct the Eq. (3.21). As an example, take the row-norm into consideration, and make a more stringent treatment, then we can easily get a relationship, such as Eq. (3.22) as follow. max fjC nþ1 j þ jC nþ3 jg P max R ðjakn;gu j þ jbkn;gu j þ jmkn;gu j þ jqkn;gu j 16i62NG
16i62NG
u¼x;y;z
ð3:22Þ
3Rk0;g P R ðjak0;gu j þ jbk0;gu j þ jmk0;gu j þ jqk0;gu j þ jRk0;gu j þ jRk2;gu jÞ ð3:23aÞ 2 5
Rk2;g þ Rk0;g P R ðjak2;gu j þ jbk2;gu j þ jmk2;gu j þ jqk2;gu j þ jLk0;gu j þ jLk2;gu jÞ u¼x;y;z
ð3:23bÞ where, n = 0,2. And ‘‘k” represent the number of node, ‘‘g” represent the energy group. Analyzing the inequality Eq. (3.23), the left side of the inequality is a constant while the right side is given by the pseudo-diffusion coefficients Dk;FDM n;guþ and the coupling corrective
According to Eqs. (3.16) and (3.17), we have,
kN1 k kMk P kN1 Mk ¼ kE1 k P 1
" " # # 0 0 M N P 2 M2 N
u¼x;y;z
We know that only if the norm of the iterative matrix is less than one can the calculation reaches convergence. So we should have, 1
7 7 7 7 7 7 7 7 5
The more explicit expression in each nodes of one energy group show as follows,
Let,
"
7 7 7 7 7 7 7 7 5
ð3:19Þ
ð3:20Þ
If kMk P kNk; that is Eq. (3.21), the iteration will reaches into convergence.
coefficients Dk;NOD n;guþ . Because the pseudo-diffusion coefficients are constants in each node, the coupling corrective coefficients should be as small as possible to satisfy the inequality. However, there are two coupling coefficients in each relationship while only one net current can be used, so there is a lot of margin to adjust the coupling corrective coefficients, such as firstly assume that there is functional relationship between the two coupling coefficients, and then combine the functional relationship and coupling corrective relationship together to work out the coefficients. Through the above derivation process on SP3 equation, a brief analysis of the SPN theory is considered. The second-order form such as Eq. (2.9) is useful because it makes the SPN equations look
925
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927
like a set of coupled diffusion equations. It can be deduced into a matrix form (Ciolini et al., 2002) and shows as follow,
A
rt
!
!
!
r2 / þrt / ¼ Q
ð3:24Þ
Z mknu1 ðsinhÞ ¼ 1=N1
where, !
/ ¼ ½/0 ; /2 ; ; /ðNþ1Þ=2 and;
!
Q ¼ ½Q ; 0; ; 0:
ð3:25Þ
where, A is an (N + 1)/2 order square, whose non-vanishing elements are, nðn1Þ An;n1 ¼ ð2nþ1Þð2n1Þ ;
2
2n þ2n1 An;n ¼ ð2nþ3Þð2n1Þ ;
ðnþ1Þðnþ2Þ An;nþ1 ¼ ð2nþ1Þð2nþ3Þ
n ¼ 1; 2; ; ðN þ 1Þ=2
1
Z mknui ðcoshÞ ¼ 1=Ni
1
sinhðakn;gu tÞpn1 ðtÞdt
1
1
coshðakn;gu tÞpni ðtÞdt
Ni ¼ 2=ð2i þ 1Þ i ¼ 0; 1; 2 Solving the simultaneous equations of the 0th-order moment, the 1st-order moment, the 2nd-order moment and the boundary continuity condition, the net neutron currents of the interface are obtained. Dk
J kn;guþ ¼ Dun;g=2 ðakn;gu1 þ 3akn;gu2 þ Hkn;gu akn;gu3 þ Gkn;gu akn;gu4 Þ k
ð3:26Þ
Dk
J kn;gu ¼ Dun;g=2 ðakn;gu1 3akn;gu2 þ Hkn;gu akn;gu3 Gkn;gu akn;gu4 Þ
ð3:31Þ
k
The neutron balance equation of the node can be obtained by integrating the SPN equations in node, ! ! B ! J þrt / ¼ Q Du
ð3:27Þ
And, (n = 0,2)
Hkn;gu ¼
akn;gu coshðakn;gu Þmknu1 ðsinhÞ sinhðakn;gu Þmknu1 ðsinhÞ ak
!
where, J ¼ ½J 0;uþ J0;u ; J2;uþ J2;u ; ; J ðNþ1Þ=2;uþ J ðNþ1Þ=2;u , and matrix B can be simply calculated from the expression of A. Then we substitute the new coupling corrective relationship Eq. (2.16) into the neutron balance equation, the same matrix equations as Eq. (3.8) is get and the subsequent derivation process is the same as the above SP3 equation. Therefore, the main point after derivation is also can be found that, for the numerical stability, to minimize the coupling correction factor as possible and the difference on the net neutron flux should be distributed to each coupling correction factor. 3.3. The Semi-analytical nodal method The code NLSP3 has selected the semi-analytical nodal method (Zimin et al., 1998) to calculate the net currents on the surface. The semi-analytical nodal method has higher precision and efficiency in the calculation of multi-group equations. In this paper, the process of the semi-analytical nodal method in SP3 is briefly introduced. The transverse integral equations of SP3 is marked as Eq. (3.28),
D0 @u@ 2 /0u ðuÞ þ R0 /0u ðuÞ 2R0 /2u ðuÞ ¼ S0u ðuÞ L0u ðuÞ D2 @u@ 2 /2u ðuÞ þ R2 /2u ðuÞ 25 R0 /0u ðuÞ ¼ 25 S0u ðuÞ L2u ðuÞ ð3:28Þ
Gkn;gu ¼ coshðn;gu ak
sinhðakn;gu Þ3mknu2 ðcoshÞ
With the Eqs. (3.4) and (3.31), the coupling corrective coefficients are obtained. After the global calculation of CMFD formula, the angular flux and power distribution are worked out. However, the 2nd-order flux of SP3 is hundreds of times smaller than the 0thorder flux, according to Eq. (3.23) and for numerical stability, it is necessary to avoid the 2nd-order flux being a divisor, otherwise the coefficients will be too large. There are many ways to deal with it. As an example, the coupling corrective relationship shows on k;NOD Dk;NOD 4;guþ for 2nd-order flux, while D2;guþ can be set to be zero, and that is what exactly the code NLSP3 do. Thus, the numerical instability problem of nonlinear iteration is solved.
4. Numerical results 4.1. The 2D-IAEA transport benchmark This benchmark (Stepanek and Auerbach, 1983) simulates a mono-energetic 5-zones light water reactor, surrounded by vacuum boundary conditions. Geometric structure shown in Fig. 2, and regional material parameters in Table 1. In this paper, we use NLSP3 program to calculate the two different meshing of the problem, the results shown in Table 2. The results which is
With the Polynomial expansion form,
/kn;gu ðuÞ ¼ /kn;g þ
4 X 2u akn;gui pni Duk i¼1
ð3:29Þ
The 4th-order polynomial of the semi-analytic nodal method shows as follows:
Pn0 ðuÞ ¼ 1 u Pn1 ðuÞ ¼ Du=2 ¼t
Pn2 ðuÞ ¼ 12 ð3t2 1Þ Pn3 ðuÞ ¼
sinhðakn;gu tÞmknu1 ðsinhÞpn1 ðtÞ sinhðakn;gu Þmknu1 ðsinhÞ
ð3:30Þ
coshðakn;gu tÞmknu0 ðcoshÞpn0 ðtÞmknu2 ðcoshÞpn2 ðtÞ
Pn4 ðuÞ ¼ u 2 a2u ; a2u
aknu
coshðakn;gu Þmknu0 ðcoshÞmknu2 ðcoshÞ
vffiffiffiffiffiffiffiffiffi u k uRn;g Duk ¼t k Dn;g 2
ð3:32Þ
k k n;gu Þmnu0 ðcoshÞmnu2 ðcoshÞ
ðn ¼ 0; 2Þ Fig. 2. Structure of 2D-IAEA Benchmark.
926
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927
Table 1 Material parameters of 2D-IAEA Benchmark. zones
Rt =cm1
Rs =cm1
mRf =cm1
1 2 3 4 5
0.60 0.48 0.70 0.65 0.90
0.53 0.20 0.66 0.50 0.89
0.079 0.0 0.043 0.0 0.0
compared to the results from SURCU (Stepanek and Auerbach, 1983), DNTR (Li et al., 2010) and EFEN (Li et al., 2011) show that NLSP3 has the same accuracy as other transport programs.
4.2. The 3D-TAKEDA benchmark The selected model is the second one in 3D-TAKEDA benchmarks (Takeda and Ikeda, 1991). A small fast breeder reactor is
Table 2 Results of 2D-IAEA Benchmark. Program
mesh
Flux in zone 1
Flux in zone 2
Flux in zone 3
Flux in zone 4
Flux in zone 5
keff
SURCU DNTR(S4) EFEN NLSP3 NLSP3
– 1148D 72 64h 16 16h 96 86h
0.01686 0.01686 0.01686 0.01686 0.01686
0.000125 0.000125 0.000125 0.000124 0.000124
0.000041 0.000035 0.000040 0.000027 0.000036
0.000295 0.000295 0.000297 0.000294 0.000294
0.000791 0.000791 0.000795 0.000715 0.000789
1.0083 1.0085 1.0078 1.0088 1.0089
Note: D Triangular meshing; hGraphical meshing
Fig. 3. Structure of 3D-TAKEDA Benchmark.
Table 3 Results of 3D-TAKEDA Benchmark. Method keff and errors (105)
Exact Monte-Carlo
NSPn 5 5 5(cm)
NLSP3 5 5 5(cm)
0.9732
363
0.9734
20
group
Ref flux
Errors (%)
flux
Errors (%)
fuel zone
1 2 3 4
4.2814E05 2.4081E04 1.6411E04 6.2247E06
0.065 0.014 0.240 0.678
4.2680E05 2.4043E04 1.6388E04 6.2058E06
0.3130 0.1568 0.1378 0.3043
axial propagation zone
1 2 3 4
5.1850E06 4.6912E05 4.6978E05 3.7736E06
2.065 1.272 0.888 0.229
5.2270E06 4.6962E05 4.6944E05 3.7366E06
0.8094 0.1072 0.0715 0.9802
radial propagation zone
1 2 3 4
3.3252E06 3.0893E05 3.2834E05 2.0473E06
2.351 1.110 0.879 0.698
3.3568E06 3.0851E05 3.2765E05 2.0423E06
0.9507 0.1350 0.2095 0.2451
sodium filling zone
1 2 3 4
2.5344E05 1.6658E04 1.2648E04 6.9840E06
3.398 0.361 0.360 0.917
2.6053E05 1.6647E04 1.2658E04 6.9970E06
2.7975 0.0637 0.0809 0.1865
Q. Pan et al. / Annals of Nuclear Energy 110 (2017) 920–927
simulated. There are four material zones, namely fuel zone, axial propagation zone, radial propagation zone and sodium filling zone while the control rod has a full lift. The geometrical structure is shown in Fig. 3. The cross-section parameters are described in references and the reference solution is obtained by the exact Monte Carlo method. We used the NLSP3 program to calculate the benchmark under the control rod full-lift condition (5 5 5 cm). The effective multiplicative coefficient and relative power distribution were listed in Table 3. It can be seen from the calculation results that the errors of the effective multiplication coefficient is very small and the calculated results are in good agreement with the reference solution except the flux errors at the sodium-filled zone. The results also show that the NLSP3 program based on nonlinear iteration has a better accuracy than the NSPn (Li et al., 2010) based on the exponential function method in the case of the same meshing for 3D-TAKEDA. 5. Conclusion Due to the conventional coupling corrective relationship’s defect, the numerical instability problem is found during the coupling correction of the higher-order flux when adopting the conventional nonlinear iterative method to solve the SPN equation. Even though in other nodal approaches for SPN equation, the numerical instability problem is still existing when CMFD formula is used to speed up the calculation. To solve the numerical instability problem, a new nonlinear iterative method for SPN theory is put forward by importing a new coupling corrective relationship which focuses on the coupling at the angle. Compared to the conventional nonlinear iterative method, the new method has better application prospect. Except for the calculation at the vacuum boundary, the coupling corrective coefficients have many ways to calculate, which brings in more margin to deal with the numerical instability. In this paper, some inequalities for how to confirm the exact values of the coupling corrective coefficients are given through numerical derivation about the norm of the iterative matrix. The inequality can give the guide and narrow down the range to work out the coupling corrective coefficients. Under the premise of
927
satisfying the inequality, the functional relationships between each coupling corrective coefficients can be selected freely according to the different situations. With that constructor function, the coefficients will be solved by the coupling corrective relationship, and the convergent results are obtained by substituting the coefficients into the CMFD equation. As an example, the code named NLSP3 for solving the SP3 equations with the new nonlinear iterative method was developed, and the numerical results demonstrate the feasibility of the new method. References Berkert, C., Grundmann, U., 2008. Development and verification of a nodal approach for solving the multigroup SP3 equations. Ann. Nucl. Energy 35, 75–86. Brantley, P.S., Larsen, E.W., 2000. The simplified P3 approximation. Nucl. Sci. Eng. 134, 1–21. Chao, Y.A., 2015. A new SPN theory formulation with self-consistent physical assumptions on angular flux. Ann. Nucl. Energy 87, 137–144. Ciolini, R., Coppa, G.G.M., Montagnini, B., Ravetto, P., 2002. Simplified PN and AN methods in neutron transport. Prog. Nucl. Energy 40 (2), 237–264. Fu, X.D., 2001. A numerical stabilization method for nonlinear analytic nodal method. Chin. J. Comput. Phys. 18 (2), 111–114. Gelbard, E.M., 1960. ‘‘Application of spherical harmonics method to reactor problems”. Technical Report WAPD-BT-20, Bettis Atomic Power Laboratory. Li, Y.Z., Wu, H.C., Cao, l.z., 2010. Nodal SP3 method for solving neutron transport equation. Nucl. Power Eng. 31, 73–78. Li, Y.Z., Wu, H.C., Cao, l.z., 2011. ‘‘A polygonal nodal-sp3 method for whole core pinby-pin calculation” // M&C2011. Rio de Janeiro. Brazil. Smith, K.S., 1983. Nodal method storage reduction by nonlinear iteration. Trans. Am. Nucl. Soc 44, 265. Stepanek, J., Auerbach, T., 1983. ‘‘Calculation of four thermal reactor benchmark problems in X-Y geometry”. EPRI NP-2855. Takeda T., Ikeda H., 1991. ‘‘3-D Neutron Transport Benchmarks[R]. OECD/NEA Committee on Reactor Physics” (NEACRP). NEACRP-L-330. Japan:OSAKA University. Yu, L.L., 2014. In-Depth Investigation and Further Development of Homogenization Method and Discontinuity Factor Theory[D]. Shanghai Jiaotong University, Shanghai. Yu, L.L. et al., 2014. The calculation method for SP3 discontinuity factor and its application. Ann. Nucl. Energy 69, 14–24. Zhang, Y.H., 2013. Iterative performance of various formulations of the SPN equations. J. Comput. Phys. 252, 558–572. Zimin, V.G., Ninokata, H., Pogosbekyan, L.R., 1998. Polynomial and semi-analytic nodal method for nonlinear iteration procedure[C]. Proc. Int. Conf. Phys, Nucl. Sci. Tech. (PHYSOR-98), Long Island, NY, October 5–8, pp. 994–1002.