A new procedure for solving neutron transfer problems

Annals of Nuclear Energy 138 (2020) 107141

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

A new procedure for solving neutron transfer problems Jiannan Tang a,b, Mei Huang a,b,⇑, Yuanyuan Zhao b, Xiaoping Ouyang b,c, Chihiro Morita d a

State Key Laboratory of Disaster Prevention and Reduction for Power Grid Transmission and Distribution Equipment, State Grid Hunan Electric Power Co., Hunan 410000, China School of Nuclear Science and Engineering, North China Electric Power University, Beijing 102206, China c Northwest Institute of Nuclear Technology, Xi’an 710024, China d Dept. of Civil and Environmental Eng., University of Miyazaki, Japan b

a r t i c l e

i n f o

Article history: Received 26 May 2019 Received in revised form 30 August 2019 Accepted 16 October 2019

Keywords: Half-boundary method Neutron transfer problems Boltzmann equation

a b s t r a c t In this paper, we propose an accurate and efficient procedure named half boundary method (HBM) to solve one-dimensional neutron transfer problems. In the process of solving the Boltzmann transfer equation, the relationships between every two neighboring nodes are solved first. Then, starting from a vacuum boundary, the neighboring relationships are iteratively applied to solve the neutron distribution. In two reflective boundaries problems, the reflective boundaries are applied to solve the neutron flux at the boundaries first, then the neutron distribution is solved similar to the vacuum boundary problems. Different from the traditional method, no matter how many segments are dispersed, at most two-order matrixes are involved in HBM calculations. So, HBM shows great potential in saving computational memory, particularly for problems with huge girds. In this paper, we introduce the fundamental theory and investigate the applicability, accuracy, and efficiency of HBM by solving five examples with different situations and boundaries. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction The particle transfer problem is one of the most critical problems in nuclear engineering. Particularly, the neutron distribution affects the heat generation distribution, furthermore influences the temperature distribution of the nuclear core and the radiation range, and finally impacts the safety and lifetime of the whole fission nuclear plant (Duderstadt et al., 1976). In methods to simulate the neutron transfer problems, deterministic and nondeterministic methods are the important and popular methods to solve neutron transfer problems. As a nondeterministic method, Monte Carlo methods (Dunn and Shultis, 2012) are widely applied in nuclear engineering such as the detector field, mediation field, and nuclear reactor, etc. However, in order to obtain accurate enough results, Monte Carlo methods need to cost tons of time because massive samples need to be calculated, particularly in deep penetration problems (Zheng et al., 2018). Compared with the nondeterministic methods, the deterministic method could solve these problems without samples, thus massive calculation time could be saved. However, along with the higher requirement of precision and more ⇑ Corresponding author at: State Key Laboratory of Disaster Prevention and Reduction for Power Grid Transmission and Distribution Equipment, State Grid Hunan Electric Power Co., Hunan 410000, China. E-mail address: [email protected] (M. Huang). https://doi.org/10.1016/j.anucene.2019.107141 0306-4549/Ó 2019 Elsevier Ltd. All rights reserved.

complex geometric model, deterministic methods cast more and more computational resource. In order to simplify the problems, the diffusion approximation is used (Hébert et al., 2016). Nevertheless, the diffusion approximation is valid only for regions that are much more scattering than absorbing (Capilla et al., 2005). Thus, in problems of the nonscattering region the diffusion approximation fails. Therefore, the Boltzmann transfer equation needs to be solved to widen applications and obtain results that are more accurate. In order to solve these problems, the angular neutron flux needs to be solved first, thus some methods such as the discrete ordinates method (SN) and spherical harmonics (PN) are developed, and they are two of the most popular methods. In SN, there are a finite set of angular directions in the unit sphere, and their corresponding weights define the appropriate quadrature. This method was applied in several codes such as PENTRAN (Sjoden and Haghighat, 1996), TORT (Rhoades and Childs, 1991) and DRAGON (Marleau et al., 2016). Besides, in PN, it is assumed that the angular dependence of the neutron flux and the cross-sections can be taken into account by means of a finite expansion in terms of Legendre polynomials (Stacey, 2001). Then, SN and PN need to be coupled with some spatial discrete methods like the finite difference method (Zhu et al., 2016; Vagheian et al., 2016), nodal collocation method (Capilla et al., 2016; Jamalipour et al., 2018), finite element method (Aydin et al., 2005; Yousefi et al., 2017), finite volume method (Hu

2

J. Tang et al. / Annals of Nuclear Energy 138 (2020) 107141

et al., 2017; Ge et al., 2015), etc., the Boltzmann transfer equation finally could be solved. With the increase in precision, the traditional methods such as the FDM, FEM, and FVM need to solve matrices with higher-order, and massive computational memory will be occupied. Therefore, HBM is proposed to solve this problem. In previous works, the HBM (Huang et al., 2013, 2017; Tang et al., 2017, 2018a,b) was developed for the heat conduction equations, this method was implemented into some computer codes which could solve the thermal problems of the double vessels and the nuclear fuel rod of the nuclear reactor. In this work, this method is applied to study the neutron transfer problems in different geometries and various boundary conditions. Besides, neutron transport and temperature are coupled to obtain results that are more accurate and realistic. 2. Half boundary method The time-independent Boltzmann transfer equation can be written as follows:

!

!

X r/ð! r ; X ; EÞ þ

X! !! ð r ; EÞ/ð r ; X ; EÞ

Fig.2. Schematic of the discrete direction.

t

!! !! !! ¼ Q s ð r ; X ; EÞ þ Q f ð r ; X ; EÞ þ Sð r ; X ; EÞ

ð1Þ

! !! where Uð r ; X ; EÞ is the neutron angular flux in location r ¼ ðx; y; zÞ P ! (in Cartesian coordinates) and direction X in energy E status, t is the total macroscopic cross-section, Qs, Qf, and S are the scattering source term, the source of neutrons by fission term and the fixed source term, respectively. In this paper, what we paid attention is the one-dimensional neutron transfer problems with steady-state, single energy group, and fixed neutron source, thus Eq. (1) could be simplified to:

l

@/ðx; lÞ X þ /ðx; lÞ ¼ Sðx; lÞ @x t

ð2Þ

where, x is the position in Cartesian coordinates, l ¼ cosh and h is the angle between the direction of the neutron velocity and the xaxis positive direction. In order to solve the neutron problems, first we disperse the one-dimensional geometry model in both space and direction of the neutron velocity, the schematics are shown in Figs. 1 and 2, respectively. Fig. 1 shows us that the one-dimensional model is dispersed into (n-1) parts, there are n nodes, besides the distance between arbitrary node i and (i + 1) is Dxi ¼ xðiþ1Þ  xi . Similar to the strategy of space-dispersion method, due to the geometrical symmetry on x-axis, the four quadrants are dispersed into M parts, an arbitrary part Xm is shown in Fig. 2, hm is the angle between the ! direction of the X m and the x-axis (0 < h < 2p), and the range of ! ! X m is DXm . Then Eq. (2) in the direction of X m could be written as:

lm

@/ðx; lm Þ X þ /ðx; lm Þ ¼ Sðx; lm Þ @x t

Z

xðiþ1Þ

xi

@/ðx; lm Þ dx þ @x

Z

xðiþ1Þ

xi

X

/ðx; lm Þdx ¼

Z

t

xðiþ1Þ xi

Sðx; lm Þdx ð4Þ

The first term in Eq. (4) is

R xðiþ1Þ xi

@/ðx;lm Þ dx @x

¼ /ðiþ1Þ;m  /i;m ,

where /i;m indicates the neutron angular flux of node i in lm Rx P direction. The second term in Eq. (4) is xiðiþ1Þ t /ðx; lm Þdx ¼ P P 1 t;i 2 ð/ðiþ1Þ;m þ /i;m ÞDxi , where the t;i is the total macroscopic cross-section between node i and node (i + 1). The last term in   Rx Eq. (4) is xiðiþ1Þ Sðx; lm Þdx ¼ Si;m Dxi , where Si;m is the average fixed ! neutron source between i and node (i + 1) in the X m direction. Thus the Eq. (4) could be simplified as follows:

lm ð/ðiþ1Þ;m  /i;m Þ þ

X1  ð/ðiþ1Þ;m þ /i;m ÞDxi ¼ Si;m Dxi 2 t;i

ð5Þ

By rearranging the variables at node (i + 1) to the left side of Eq. (5), Eq. (5) becomes:

X X  1 1 ðlm þ Dxi Þ/ðiþ1Þ;m ¼ ðlm  Dxi Þ/i;m þ Si;m Dxi 2 2 t;i t;i

ð6Þ

Eq. (6) is the relationship between arbitrary node i and node (i + 1). Thus, if we repeatedly use Eq. (6), then we can obtain the relationship between node 1 and node n (n > 2) in the model shown in Fig. 1, and the relationship could be expressed as below:

/n;m ¼

ð3Þ þ

Base on the dispersed geometrical model, we integrate Eq. (3) from node xi to node x(i+1), the following expression is obtained:

lm

2 P

n1 Q

lm 0:5Dxk

k¼1

lm þ0:5Dxk

½ð

l Q

l¼ðn1Þ h¼ðn1Þ

Fig.1. Schematic of the discrete space.

P Pt;k /1;m t;k

lm 0:5Dxh lm þ0:5Dxh

P Pt;h Þ t;h





S l;m Dxðl1Þ

lm þ0:5Dxðl1Þ

P t;ðl1Þ

þl

S ðn1Þ;m Dxðn1Þ

m þ0:5Dxðn1Þ

P

t;ðn1Þ

ð7Þ

3

J. Tang et al. / Annals of Nuclear Energy 138 (2020) 107141  P If the t;i , Dxi and Si;m are location independent and constant, the Eq. (7) would be simplified as:

 /n;m ¼

am bm

n1 /1;m þ

1

 n1 am bm

bm  am



Si;m Dxi

ð8Þ

P P where, am ¼ lm  0:5Dxi t;i , bm ¼ lm þ 0:5Dxi t;i . It is obvious, compare with Eq. (8), Eq. (7) could be applied to more problems that have a varied total cross-section, neutron source and space step, however, Eq. (8) could save more computational resource and short the calculation time for homogeneous problems. Therefore, we can choose an appropriate equation to solve the particular neutron transfer problem. We can also couple Eq. (7) and Eq. (8) together to solve one issue. After we obtain the relationship between node 1 and node n, we still need the boundaries which are located in node 1 and node n to solve the equation. There are two kinds of common boundaries which are vacuum and reflective. In the vacuum boundary, the neutron angular flux with a direction which toward the inside of the model is zero. Besides, in the reflective boundary, the neutron angular flux has the velocity toward the inside of the model is on x axisymmetric. A least one of the two boundaries in node 1 or n is a vacuum boundary, above Eqs. (7) or (8) is solvable. However, what if all two boundaries are reflective, in most existing works, such as FEM, a fake value of neutron angular flux need to be set first, and then stable results are obtained by applying iterative computations. In this paper, we solve the two-reflective boundaries problems directly, no iteration involved, thus massive computational time could be saved. Fig. 3 shows us the angular schematic of node ! ! ! ! i in the x-y plan, the X m1 , X m2 , X m3 and X m4 are symmetric, because of the reflective boundary in node i, thus /i;m1 ¼ /i;m2 or /i;m4 ¼ /i;m3 . Then, we apply Eq. (7) (or Eq. (8)) to obtain the rela! ! ! tionship between node 1 and node n in X m1 and X m2 ( X m3 and ! X m4 ) respectively, and we will get two equations. By rearranging the variables to the left side of the above two equations into a matrix, and further simplify them into a matrix form:



2 3 /n;m1 6   S1n;m1 1 cn;m1 0 0 6 /1;m1 7 7 6 7¼ 0 0 1 cn;m2 4 /n;m2 5 S1n;m2 /1;m2

Fig.3. Schematic of the discrete direction.

P P Q lm1 0:5Dxk t;k Qn1 lm2 0:5Dxk t;k P P , S1n;m1 ¼ where, cn;m1 ¼ n1 , c ¼ n;m2 k¼1 l þ0:5Dx k¼1 l þ0:5Dx k k m1 m2 t;k t;k P      lm1 0:5Dxh t;h Ql P2 S l;m1 Dxðl1Þ S ðn1Þ;m1 Dxðn1Þ P P þ l þ0:5Dx P l¼ðn1Þ h¼ðn1Þ l þ0:5Dx lm1 þ0:5Dxðl1Þ t;ðl1Þ ðn1Þ h m1 m1 t;ðn1Þ P   t;h   lm2 0:5Dxh t;h P2 Ql S l;m2 Dxðl1Þ P P and S1n;m2 ¼ l¼ðn1Þ þ h¼ðn1Þ l þ0:5Dx l þ0:5Dx m2



S ðn1Þ;m2 Dxðn1Þ

lm2 þ0:5Dxðn1Þ

P

.

Because

of

the

h

t;h

reflective

m2

ðl1Þ

t;ðl1Þ

boundary,

thus

t;ðn1Þ

/n;m1 ¼ /n;m2 and /1;m1 ¼ /1;m2 . Further, simplify Eq. (10), and we can obtain the following equation:



1 cn;m1 1 cn;m2

"

/n;m1 /1;m1

#

 ¼

S1n;m1



S1n;m2

ð11Þ

There are two equations and two variables in Eq. (11), since it is solvable, and the iterative calculations are avoided. In order to illustrate the computational processes, the flowchart is shown in Fig. 4. First, the geometry, material properties, numerical model and boundaries are put in the calculation. Then we need to calculate the angular neutron flux in the boundaries by HBM. After that, we get the results (i is from 1 to n) as requirements by HBM. Finally, we solve the angular neutron flux at all angles (m is from 1 to M) and plot the results we want. The HBM is very flexible, different from the traditional methods, we can obtain the results we want without solving all unknown variables. Need to be noted that the angular neutron flux at the direction inside the model is known (it is zero) in vacuum boundary, the boundary calculation procedure could be skipped. In this section, we have elaborated the derivation process of the HBM for steady neutron transfer problems with constant energy and fixed neutron source. In addition, the neutron transfer problems which include multiple energy groups and neutron sources from the scattering and fission could be easily solved base

ð10Þ

Fig. 4. Flow chart for solving neutron transfer problems with HBM.

4

J. Tang et al. / Annals of Nuclear Energy 138 (2020) 107141

on this fundamental research, and we will develop the application of HBM for solving more complicated neutron problems in future works.

/ðx; lÞ ¼

P

8 > > S

> :

1

e l

t

! ðx2LþsgnðlÞ2LÞ ;

if

l–0

ð14Þ

t

S P ;

if

l¼ 0

t

where, sgnðlÞ ¼

3. Results and discussions After illustrating the process of the HBM in neutron transfer problems, in this part, two numerical examples are tested to show the applicability of HBM with different boundaries as compared with the analytical solution. Besides, an example with a more complex geometric model and an example which couples thermal distribution is simulated. All numerical calculations of onedimensional problems are carried out in Matlab 2016a software with double precision in which the digit’s environment variable is assigned to be 64, and all the calculations are performed on a PC with a 3.60 GHz CPU and 8.0 GB RAM. Example 3.1. As the first example, the one-dimensional steadystate neutron transfer problem with a constant total cross-section P ( t ¼ 1 cm1) and a uniform fixed neutron source (S = 1 n/cm2 s) is considered to assess the veracity and efficiency of HBM. The geometric sketch is shown in Fig. 5. The thickness of the infinite plate is L = 1 cm. The boundary of face A and face B are all vacuum boundaries. In order to validate the accuracy of the HBM in neutron transfer problems, we need to obtain the exact results first as the reference. The mathematical description of this example 3.1 is as follow:

8 > <

P

lÞ l @/ðx; þ t /ðx; lÞ ¼ Sðx; lÞ @x /ð0; lÞ ¼ 0; for 0 < l  þ1 > : /ðL; lÞ ¼ 0; for  1  l < 0

If

ð12Þ

l–0 is kept fixed, the general solution of Eq. (12) is:

/ðx; lÞ ¼ CðlÞe

P l

t

S ðx2LÞ þ P

ð13Þ

t

where, CðlÞ is a function only depends on l. If we put two boundaries into Eq. (13), we will get CðlÞ, and then the solution of Eq. (12) is given by:

Fig. 5. Geometric sketch for example 3.1.

/ðxÞ ¼

1 2

Z

þ1

1

l . jlj

Then, the neutron distribution is given by:

S 1 /ðx; lÞdl ¼ P ½1  2 t

Z

1 0

P e

t l x

P þ e

t l ðLxÞ

! dl ð15Þ

After obtaining the exact results, we simply illustrate the process of how to apply the HBM into this example with the instructions from the flow chart in Fig. 4. In Example 3.1, it is a onedimensional neutron transfer problems with uniform macroscopic cross-section, constant fixed neutron source, and two vacuum boundaries. Then we need to choose the n and M. Because the number of M will affect the accuracy of HBM significantly. In order to verify the sensibility of M, we will calculate example 3.1 with different values of M and variable angular dispersion schemes. Now, we first just set M = 4. Because of the vacuum boundaries, thus we have obtained the /ð0; lm Þ when lm > 0 and /ðL; lm Þ when lm < 0 now we actually have finished the left part in Fig. 4. After that, if we want to obtain the angular neutron flux results in every 0.1 cm in this problem, we set n = 11, by using Eq. (6) iteratively, we can obtain the angular neutron flux in direction lm . Finally, we just need to solve the angular neutron flux in other angles, respectively. If we want the neutron distribution in x, we just need to sum the angular neutron flux of every node together. We have illustrated the whole process of how to apply HBM in example 3.1. Now, the results got by HBM with different M as compared with the analytical solution are shown in Fig. 6 (the n = 21), and the results are put in Table 1. As it is shown in Fig. 6 and Table 1, when M = 4, the results are unacceptable, the biggest error is more than 24.8% on where x = 0.3 cm. As the value of M increases, the average error decrease, particular when M increases from 4 to 8, the average error reduces more than 16%, it is a significant improvement. However, when the M is more then 16, if we continue to increase the value of M, the accuracy only raise a little. In addition, from Fig. 6, we can see that when M = 16, there are 8 results don’t coincide the exact curve very well. We further test the results with M = 20, M = 40 and M = 60, and there is a small difference between the results from M = 40 and M = 60, so we choose 40 as the value of M in this paper. It is obvious that the value of M is an important parameter to improve the accuracy of neutron transfer problem with HBM, larger M will bring more accurate results, but this influence decreases when M is already big enough. Besides, the number of angle deviation M is irrelevant to the order of matrixes during the calculating process, only one-order matrixes are used in this problem. As we can see from Fig. 6, the parameter M is critical in the process of solving neutron transfer problems, larger M could lead to better results. To invest the accuracy of the HBM, the results from HBM and FDM are compared and shown in Table 1, we can see that the average error decrease as the M increases and the HBM is more accurate than FDM in this problem. However, it simultaneously costs more computational resources as the value of M increasing. Therefore, we want to find a better disperse scheme in angle to reduce the computational cost. In the above works, we disperse the angle in the same DXm (scheme A, start from 0.5DXm , in other words, h1 ¼ 0:5DX). Now we have other two schemes, one scheme disperses the angle in same surface area of a fake sphere (scheme B) which is shown in Fig. 2, another scheme (scheme C) is the inverse of scheme B. In order to show the difference obviously, we set M = 16, the sketch of both schemes A and scheme B is shown in

5

J. Tang et al. / Annals of Nuclear Energy 138 (2020) 107141

Fig. 6. Neutron flux distribution with different M.

Table 1 Neutron flux distribution with M = 4, M = 8, M = 16, and M = 40. M

Location (cm)

Exact results

FDM

HBM

Result

Error (%)

Average Error (%)

Result

Error (%)

Average Error (%)

4

0.0 0.1 0.2 0.3 0.4 0.5

0.42575 0.55253 0.61247 0.64797 0.66722 0.67336

0.34381 0.36184 0.37619 0.39146 0.40617 0.41721

19.24596 34.51219 38.57812 39.58734 39.12487 38.04129

34.8483

0.37851 0.42601 0.46201 0.48721 0.50214 0.50708

11.09571 22.89830 24.56610 24.80979 24.74146 24.69407

22.13417

8

0.0 0.1 0.2 0.3 0.4 0.5

0.42575 0.55253 0.61247 0.64797 0.66722 0.67336

0.42284 0.46482 0.52338 0.56826 0.59979 0.60614

0.68423 15.87453 14.54637 12.30175 10.10579 9.98243

10.58252

0.42458 0.50762 0.56741 0.60767 0.63087 0.63845

0.27481 8.12807 7.35710 6.21942 5.44798 5.18445

5.43530

16

0.0 0.1 0.2 0.3 0.4 0.5

0.42575 0.55253 0.61247 0.64797 0.66722 0.67336

0.42562 0.52506 0.60038 0.64204 0.66605 0.67105

0.03041 4.97158 1.97485 0.91468 0.17592 0.34308

1.401753

0.42584 0.53395 0.60263 0.64467 0.66740 0.67459

0.02114 3.36271 1.60661 0.50928 0.02698 0.18267

0.95157

40

0.0 0.1 0.2 0.3 0.4 0.5

0.42575 0.55253 0.61247 0.64797 0.66722 0.67336

0.42568 0.55102 0.61129 0.64714 0.66614 0.67310

0.01550 0.27301 0.19302 0.12756 0.16129 0.03842

0.13480

0.42570 0.55089 0.61354 0.64880 0.66759 0.67356

0.01174 0.29682 0.17470 0.12809 0.05545 0.02970

0.11608

Fig. 7, and the results with different schemes are shown in Fig. 8. In Fig. 7, the uniform angles which are dispersed by the solid lines are the scheme A, scheme A includes four parts which are DXA1 , DXA2 , DXA3 and DXA4 in the first quadrant, and hA1 ¼ p=16, hA2 ¼ 3p=16, hA3 ¼ 5p=16 and hA4 ¼ 7p=16. At the same time, green dash lines in Fig. 7 draw the scheme B, they have the same surface area, scheme B includes four parts which are DXB1 , DXB2 , DXB3 and DXB4 in the first quadrant, and hB1 ¼ 0:5arccosð3=4Þ, hB2 ¼0:5ðarccosð3=4Þþarccosð1=2ÞÞ, hB3 ¼0:5ðarccosð1=2Þþarccosð1=4ÞÞ

and hB4 ¼ 0:5ðarccosð1=4Þ þ p=2Þ. Scheme C is the inverse of scheme B, thus hC1 ¼ 0:5arcsinð1=4Þ, hC2 ¼ 0:5ðarcsinð1=4Þþ arcsinð1=2ÞÞ, hC3 ¼ 0:5ðarcsinð1=2Þ þ arcsinð3=4ÞÞ and hC4 ¼ 0:5ðarcsinð3=4Þ þ p=2Þ. The results are shown in Fig. 8, and the data are put in Table 2. From Fig. 8, we can see that the results from scheme A are similar to that from scheme B, in the location between 0.3 cm and 0.7 cm, the scheme A has better results, but in other location, the results from scheme B are more accurate. Scheme A has a slightly better

6

J. Tang et al. / Annals of Nuclear Energy 138 (2020) 107141 Table 2 Neutron flux distribution with scheme A, scheme B, and scheme C. Scheme

Location (cm)

Exact results

HBM Result

Error (%)

Average Error

A

0.0 0.1 0.2 0.3 0.4 0.5

0.42575 0.55253 0.61247 0.64797 0.66722 0.67336

0.42584 0.53395 0.60263 0.64467 0.66740 0.67459

0.02114 3.36271 1.60661 0.50928 0.02698 0.18267

0.95157%

B

0.0 0.1 0.2 0.3 0.4 0.5

0.42575 0.55253 0.61247 0.64797 0.66722 0.67336

0.42244 0.53953 0.60557 0.64322 0.66287 0.66900

0.77745 2.35281 1.12659 0.73306 0.65196 0.64750

1.04823%

C

0.0 0.1 0.2 0.3 0.4 0.5

0.42575 0.55253 0.61247 0.64797 0.66722 0.67336

0.42932 0.51724 0.57992 0.62180 0.64582 0.65364

0.83852 6.38698 5.31455 4.03877 3.20734 2.92860

3.78579%

Fig. 7. Neutron flux distribution with different n.

average accuracy than scheme B as we can see from Table 2. Besides, the scheme C has the worst solution compare with scheme A and B. Overall, the angular disperses scheme does influence the accuracy, and we prove that the scheme A is the most accurate among the scheme A, B and C after comparing those schemes. More complex disperses schemes will be discussed in further works. Besides n is also an important value which influences the accuracy of results in this problem, thus the neutron flux distribution with different n is shown in Fig. 9 (M = 40), and because of symmetry, only the results from 0 cm to 0.5 cm are put in Table 3. In Table 3, we can see that the average error decrease as the n increases. When n = 5, the errors of the results in 0.25 cm (0.75 cm) are visual, but the result in 0.5 cm has much better accuracy. Besides, in Fig. 9, the results with n greater than 5 fit the curve of exact resolution very well, the errors decrease a little when we strongly increase the value of n from 21 to 101, thus n = 21 is chosen in next problems, in order to save the computational time. Furthermore, during this process, no matter how many nodes are included in the model, only one-order matrixes are used in the problem has at least one vacuum boundary. In order to deeply verify the accuracy and the efficiency of the HBM, as a comparison, some traditional method such as FDM is

applied in this example too. Compare with the exact result, the average errors of the results are 4.13749%, 1.45517%, 0.43344% and 0.20264% under n = 5, n = 11, n = 21 and n = 101 in m = 40 situations, respectively. Overall, the HBM has comparable accuracy with FDM. Example 3.2. In this example, we will solve the neutron problem with two reflective boundaries. Different from the geometric model of example 3.1, there are two parts are involved which are the region I and region II in this example. The thickness of region I and region II are LI = 1 cm and LII = 1 cm respectively, besides, the region I has a uniform fixed neutron source (S = 1 n/cm2 s), and the P constant total cross-section ( t ¼ 1 cm1) is spread all over the region I and II. The boundary of A and B are both reflective boundaries, and the geometric sketch is shown in Fig. 10.

Fig. 8. Neutron flux distribution with scheme A, scheme B, and scheme C.

7

J. Tang et al. / Annals of Nuclear Energy 138 (2020) 107141

Fig. 9. Neutron flux distribution with different n. Table 3 Neutron flux distribution with n = 5, n = 11, n = 21 and n = 101. n

Location (cm)

Exact results

FDM

HBM

Result

Error (%)

Average Error

Result

Error (%)

Average Error

5

0.0 0.25 0.5

0.42575 0.63258 0.67336

0.43457 0.67634 0.69641

2.07164 6.91770 3.42313

4.13749%

0.42733 0.66269 0.67334

0.37111 4.75987 0.00297

1.71132%

11

0.0 0.1 0.2 0.3 0.4 0.5

0.42575 0.55253 0.61247 0.64797 0.66722 0.67336

0.42897 0.5712 0.62147 0.65879 0.67125 0.6791

0.75631 3.37900 1.46946 1.66983 0.60400 0.85244

1.45517%

0.42593 0.55507 0.61610 0.65034 0.66874 0.67461

0.04228 0.45970 0.59268 0.36576 0.22781 0.18564

0.31231%

21

0.0 0.1 0.2 0.3 0.4 0.5

0.42575 0.55253 0.61247 0.64797 0.66722 0.67336

0.42955 0.5515 0.61475 0.651 0.66951 0.67564

0.89254 0.18642 0.37226 0.46761 0.34322 0.33860

0.43344%

0.4257 0.55089 0.61354 0.6488 0.66759 0.67356

0.01174 0.29682 0.17470 0.12809 0.05545 0.02970

0.11609%

101

0.0 0.1 0.2 0.3 0.4 0.5

0.42575 0.55253 0.61247 0.64797 0.66722 0.67336

0.42585 0.55437 0.61472 0.64952 0.66743 0.67485

0.02349 0.33301 0.36736 0.23921 0.03147 0.22128

0.20264%

0.42563 0.54974 0.61277 0.64830 0.66722 0.67321

0.02819 0.50495 0.04898 0.05093 0.00000 0.02228

0.10922%

and the neutron distribution is:

Fig. 10. Geometric sketch for examples 3.2.

The mathematical description of this example 3.2 is as follow:

0 8 P  P 1 > tL > sinh l 2 > txC  > S B P > @1  P  e l A; 0 < x < 2L > > > t tL sinh < l 0 1 /ðx; lÞ ¼ P  P > > tL > sinh l 2 > t ðLxÞ C B  > S P > @1  P  e l A; 2L < x < L > > : t t sinh l L

0 P 1 1 0 8  P > t ðL L Þ sinh > II I R l > 1B C S B > P P  cosh l t xC > AdlA; 0 < x < LI @1  0 @ > > > t sinh < t l LII 0 P  1 /ðxÞ ¼ > P > tL > sinh R1 B > l I > S P > P  cosh l t ðLII  xÞC Adl; LI < x < LII > 0 @ > : t t sinh

ð16Þ

l LII

ð17Þ Different from the examples above, we do not know /ðx; lÞ in any locations, what we can see is that /ð0; l > 0Þ ¼ /ð0; l < 0Þ and /ðL; l > 0Þ ¼ /ðL; l < 0Þ. Thus we need to solve the

8

J. Tang et al. / Annals of Nuclear Energy 138 (2020) 107141

Fig.11. Neutron flux distribution of example 3.2.

Table 4 Neutron flux distribution of example 3.2. Parameters

n = 21, m = 40

Table 5 Order number of the matrices in HBM and FDM of example 3.2.

Location (cm)

Exact results

HBM Result

Error (%)

Average Error

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

0.86124 0.85407 0.83096 0.78608 0.70368 0.50000 0.29632 0.21392 0.16904 0.14593 0.13876

0.86117 0.85399 0.83093 0.78646 0.70480 0.50000 0.29520 0.21354 0.16907 0.14600 0.13877

0.00813 0.00937 0.00361 0.04834 0.15916 0.00000 0.37797 0.17764 0.01775 0.04797 0.00721

0.07792%

/ðBoundary; lÞ first. Flowing the leading of Fig. 4, we can use Eq. (7) to solve this problem, then Eq. (11) is used to solve /ðBoundary; lÞ. Finally, similar to example 3.1, we obtain the neutron flux we want. The results are shown in Fig. 11 and datum are put in Table 4. Fig. 11 and Table 4 show us that the largest neutron flux is 0.86124 where x = 0.0 cm which is the reflective boundary, and the neutron flux decrease as the x increase, finally reduces to the minimum which is 0.13876 in the B surface which is also a reflective boundary. The results from HBM could fit the curve of exact resolution very well, compared with the exact results, the maximum error of HBM is 0.37797% where x = 1.2 cm, and the average error 0.07792%. Those results can prove that the HBM can solve the neutron transfer problems with two reflective boundaries very well, and the accuracy of HBM for those problems is acceptable. In addition, the HBM not only has good accuracy for solving neutron problems, it is also efficient, particularly in tow reflective boundaries problems. The order number of matrices, which are involved in the calculation processes with different numbers of n in FDM and HBM, are put in Table 5.

Number of n, m = 40

Method

Number of the matrices

iterations

Time ratio (HBM/FDM)

5

HBM FDM HBM FDM HBM FDM HBM FDM

2 5 2 11 2 21 2 101

0 7 0 11 0 12 0 12

20.58%

11 21 101

12.68% 5.87% 1.85%

As we can see from Table 5, the order number of the matrices in the computational process in HBM is constant, no matter how many nodes in the numerical model. On the contrary, the order number in FDM is increasing as the node number raises. Besides, there are iterations in FDM, however, the HBM could solve these problems directly, and thus HBM could be more efficient. Overall, the HBM could save much more computational time than FDM, particularly in two reflective boundaries problems. Need to be noted that the iterations and time ratio are variable according to the different coding schemes and convergence accuracy. Besides, only two-order matrixes are used in this problem, no matter how many segments are divided. Example 3.3. After verifying the accuracy of the application of the HBM in the neutron transfer problems with different boundaries. An example with the complex geometric model is simulated in this section. In the actual nuclear power plant, particularly, in platetype nuclear fuel assembly (Kim et al., 2013). The nuclear fuel plates are divided by coolant sections. Therefore, the HBM could solve the neutron transfer problems with multiple neutron source regions. The geometrical model is shown in Fig. 12, there are 3 coolant regions (I region has the constant total cross-section P 1 and without neutron source) and 2 neutron source t ¼ 1 cm regions (II regions have uniform fixed neutron source S = 1 n/cm2 s

9

J. Tang et al. / Annals of Nuclear Energy 138 (2020) 107141 Table 6 Material properties of example 3.4. Material properties Cross-section Fig. 12. Geometric sketch for example 3.3 and example 3.4.

P and the constant total cross-section t ¼ 1 cm1) in this model. In addition, the thicknesses of the I region and the II region are 1 cm and 1 cm, respectively. The results of neutron distributions of example 3.3 with two reflective boundaries and two vacuum boundaries are shown in Fig. 13, respectively. As it is shown in Fig. 13, both the neutron fluxes with two vacuum boundaries and two reflective boundaries fluctuate along the x-direction. The neutron fluxes increase in region II which include a fixed neutron source, and the results decrease in I regions which don’t have neutron source. Besides, compare with the results with two vacuum boundaries, the results with two reflective boundaries are larger, particularly in I regions near the boundary. However, in two II regions and the I region which in the middle, the difference between the results with two vacuum and reflective boundaries are small. Example 3.4. As one of the most important parameters, temperature could influent the material properties, deeply influent the neutron distribution. At the same time, the neutron could affect the heat source as feedback. Therefore, it is more accurate if we couple the temperature and the neutron into the numerical simulation. In previous articles, we have done some researches that HBM is applied to thermal problems. In this example, the temperature and neutron results are coupled to obtain the solution that more approach to reality. The geometrical model is shown in Fig. 12, two vacuum boundaries and two first kind temperature boundaries are applied, besides, the heat sources are only applied in II regions. The

Thermal conductive Thermal generation rate

Equation P 1 t ¼ 1 þ 0:0005T cm K ¼ 0:4 þ 0:01T W=cm  K g ¼ 3:0Neu kW=cm

3

temperatures at A surface and B surface are same, and they are 100 K. Besides, as a fundamental work, we assume that the material properties such as the thermal conductivity, neutron source, and the total macroscopic cross-section are related to the temperature, respectively. In addition, the material properties are shown in Table. 6, where T is the temperature and Neu is the neutron distribution. The flow chart in Fig. 14 illustrates the computational process of example 3.4. The temperature and neutron flux distributions of example 3.4 with and without coupling both the neutron and temperature are shown in Fig. 15. As it shows in Fig. 15 (a), after coupling the thermal problems, the neutron distribution is lower than that without coupling, but it still shows the same trend as that without coupling. However, the difference is obvious, particularly at location x = 1.7 cm and x = 3.3 cm, there are 0.09281 and 0.09287 differences between the neutron flux with and without coupling. Similar to the results of the neutron flux distribution, the temperature distribution has the same status. After coupling the neutron flux distribution, the temperature distribution is reduced, and the biggest difference is 167.50336 K which is about 35.1% difference at x = 3.1 cm. Overall, in example 3.4, the temperature influences the distribution of the neutron flux a lot, simultaneously, the neutron flux distribution also affects the temperature. Example 3.5. To deeply illustrate the application of HBM, it is applied to a realistic nuclear problem (Issa et al., 1986) in the last example. The geometrical model of the realistic example, which has one reflective and one vacuum boundary, is shown in Fig. 16.

Fig. 13. Neutron distribution of example 3.3 with different boundaries.

10

J. Tang et al. / Annals of Nuclear Energy 138 (2020) 107141

Fig. 16. The geometrical model of example 3.5.

Fig. 14. Flow chart of the computational process of example 3.4.

This nuclear fuel model includes 2 cm thickness core and 3 cm thickness shield. In the core fuel section, total cross-section P P 1 1 t ¼ 1:0 cm , scattering cross-section s ¼ 0:5 cm , the average number of neutrons from one fission multiply by fission crossP section is m f ¼ 1:0 cm1 , and the reflective boundary is applied on the left side. While, in the shield region, total cross-section P P 1 1 t ¼ 0:8 cm , scattering cross-section s ¼ 0:4 cm , fission P 1 cross-section is f ¼ 0:0 cm , and the vacuum boundary is applied on the right side. This problem is solved by HBM and MCNP, and two results are shown in Fig. 17. As it is shown in Fig. 17, the value of neutron flux reaches the top at the left side and then decreases along the x-direction. Besides, the results from HBM could fit the MCNP results well.

Fig. 17. The neutron flux of 3.5 by HBM and MCNP.

4. Conclusion In this paper, we first elaborated the derivation process of the HBM for steady neutron transfer problems with constant energy and fixed neutron source. Different from the traditional method such as FEM, no matter how many parts we disperse the geometry model, at most two-order matrixes are involved in calculations of HBM. Since the new method shows great potential in saving computational memory, particularly for neutron problems with huge girds. Then two example with vacuum and reflective boundaries were calculated by HBM, and the accuracy of the HBM was verified by comparing the results with that from the exact solution. Besides, we found that the divided numbers and the schemes in

Fig. 15. The temperature and neutron flux distributions of example 3.4 with and without coupling.

J. Tang et al. / Annals of Nuclear Energy 138 (2020) 107141

angular and space influence the solution accuracy significantly. Different from some traditional method such as FDM, the iterative calculations are avoided during the two reflective boundaries problems in HBM, thus the HBM could save more computational time in these problems. In order to investigate the applicability of HBM, the problem had more complex geometrical modes was simulated. Furthermore, we coupled the temperature into the neutron problem to obtain a more accurate and realistic solution. Finally, HBM solved a realistic problem includes core and shield, and the results were compared with commercial calculational software MCNP, and the accuracy of HBM was verified again. Overall, the HBM has the potential to simulate the neutron transfer problems accurately and efficiently. Acknowledgments The authors would like to appreciate the financial support from the NCEPU ‘‘Double First-Class” Graduate Talent Cultivation Program, the Fundamental Research Funds for the Central Universities (2018ZD10) and the State Key Laboratory of Disaster Prevention and Reduction for Power Grid Transmission and Distribution Equipment (No. B316AF190007) of China. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.anucene.2019.107141. References Aydin, E.D., Katsimichas, S., Oliveira, C.R.E.d., 2005. Time-dependent diffusion and transport calculations using a finite-element-spherical harmonics method. J. Quant. Spectrosc. Radiat. Transfer 95, 349–363. Capilla, M., Talavera, C.F., Ginestar, D., Verdu´, G., 2005. A nodal collocation method for the calculation of the lambda modes of the PL equations. Ann. Nucl. Energy 32, 1825–1853. Capilla, M.T., Talavera, C.F., Ginestar, D., Verdú, G., 2016. Nodal collocation method for the multidimensional PL equations applied to neutron transport source problems. Ann. Nucl. Energy 87, 89–100. Duderstadt, J.J., Hamilton, L.J., 1976. Nuclear reactor analysis. John Wiley & Sons Inc.. Dunn, W.L., Shultis, J.K., 2012. Exploring Monte Carlo methods. Comput. Oper. Res.

11

Ge, J., Zhang, D., Tian, W., Wang, K., Qiu, S., Su, G.H., 2015. Steady and transient solutions of neutronics problems based on finite volume method (FVM) with a CFD code. Prog. Nucl. Energy 85, 366–374. Hébert, A., Sekki, D., Chambon, R., 2016. A user guide for DONJON version4. Institut de Génie Nucléaire, École Polytechnique de Montréal, Montréal. , Technical Report IGE-294. Hu, T., Wu, H., Cao, L., Li, Z., 2017. Finite volume method based neutronics solvers for steady and transient-state analysis of nuclear reactors. Energy Procedia 127, 275–283. Huang, M., Ma, J., Tang, J.N., Yuan, H., 2013. Mixed boundary node method for free vibration analysis of rectangular plates with variable thickness and general boundary conditions. WIT Trans. Model. Simul. 56, 127–137. Huang, M., Tang, J., Zhao, Y., Ouyang, X., 2017. A new procedure for efficient and accurate solving heat conduction problems. Int. J. Heat Mass Transfer 111, 508– 519. Issa, J.G., Riyait, N.S., Goddard, A.J., Stott, G.E., 1986. Multigroup application of the anisotropic FEM code FELTRAN to one, two, three-dimensions and R-Z problems. Prog. Nucl. Energy 18, 251–264. Jamalipour, M., Cammi, A., Ricotti, M.E., 2018. The application of nodal method for dynamic analysis of TRIGA Mark II. Ann. Nucl. Energy 116, 360–375. Kim, H., Yim, J., Lee, B., Oh, J., Tahk, Y., 2013. Drop impact analysis of plate-type fuel assembly in research reactor. Nucl. Eng. Technol. 46, 529–540. Marleau, G., Ebert, A., Roy, R., 2016. A User Guide for DRAGON Version 4. Institut de Génie Nucléaire, École Polytechnique de Montréal, Montréal. Technical report IGE-294. Rhoades, W.A., Childs, R.L., 1991. TORT: a three-dimensional discrete ordinates neutron/photon transport code. Nucl. Sci. Eng. 107 (4), 397–398. Sjoden, G.E., Haghighat, A., 1996. PENTRAN: a three-dimensional scalable transport code with complete phase-space decomposition. Trans. Am. Nucl. Soc. 74. Stacey, W.M., 2001. Nuclear Reactor Physics. Wiley, New York. Tang, J., Huang, M., Zhao, Y., Ouyang, X., Huang, J., 2017. A new procedure for solving steady-state and transient-state nonlinear radial conduction problems of nuclear fuel rods. Ann. Nucl. Energy 110, 492–500. Tang, J., Huang, M., Zhao, Y., Maqsood, S., Ouyang, X., 2018a. Numerical investigations on the melting process of the nuclear fuel rod in RIAs and LOCAs. Int. J. Heat Mass Transfer 124, 990–1002. Tang, J., Huang, M., Yang, M., Zhao, Y., Ouyang, X., 2018b. A procedure for solving transient nonlinear thermal problems of high burnup nuclear fuel rods in a light water reactor. Appl. Therm. Eng. 140, 542–552. Vagheian, M., Vosoughi, N., Gharib, M., 2016. Enhanced finite difference scheme for the neutron diffusion equation using the importance function. Ann. Nucl. Energy 96, 412–421. Yousefi, M., Zolfaghari, A., Minuchehr, A., Abbassi, M.R., 2017. An arbitrary geometry finite element method for the adjoint neutron transport equation. Ann. Nucl. Energy 110, 511–525. Zheng, Z., Wang, M., Li, H., Mei, Q., Deng, L., 2018. Application of a 3D discrete ordinates-Monte Carlo coupling method to deep-penetration shielding calculation. Nucl. Eng. Design 326, 87–96. Zhu, A., Jarrett, M., Xu, Y., Kochunas, B., Larsen, E., Downar, T., 2016. An optimally diffusive Coarse Mesh Finite Difference method to accelerate neutron transport calculations. Ann. Nucl. Energy 96, 116–124.