Development of a Molten Salt Reactor specific depletion code MODEC

Annals of Nuclear Energy 124 (2019) 88–97

Contents lists available at ScienceDirect

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

Development of a Molten Salt Reactor specific depletion code MODEC Shaopeng Xia a,b,c, Jingen Chen a,b,c, Wei Guo a,b,c, Deyang Cui a,b, Jianlong Han a,b,c, Jianhui Wu a,b,⇑, Xiangzhou Cai a,b,c,* a

Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China Innovative Academies in TMSR Energy System, Chinese Academy of Sciences, Shanghai 201800, China c University of Chinese Academy of Sciences, Beijing 100049, China b

a r t i c l e

i n f o

Article history: Received 24 April 2018 Received in revised form 15 August 2018 Accepted 23 September 2018

Keywords: Molten Salt Reactor Burnup calculation Nonhomogeneous term TTA CRAM Error analysis of ORIGEN-S

a b s t r a c t Molten Salt Reactor (MSR) has the characteristics of on-line reprocessing and continuously refueling, which bring significant differences from traditional reactors in burnup calculations. To handle its specific burnup characteristics, a Molten Salt Reactor specific depletion code - MODEC has been newly developed. MODEC implements two depletion algorithms to solve the basic burnup equations: the transmutation trajectory analysis (TTA) and the Chebyshev rational approximation method (CRAM). To simulate the on-line reprocessing, the fictive decay constant method is applied. And three different methods are implemented in MODEC to solve the nonhomogeneous burnup equations in the continuously refueling problems. Moreover, it can trace in real time the evolution of the in-stockpile nuclides which is extracted by on-line reprocessing. MODEC has three calculating modes for decay, constant flux and constant power calculations. By comparing with ORIGEN-S, the validity of performance of MODEC in conventional burnup calculations, burnup calculations with on-line reprocessing and burnup calculations with continuously refueling is proved. And a comparison of the three methods of solving nonhomogeneous burnup equations is presented and discussed. Additionally, a detailed analysis of error sources in ORIGEN-S is applied and an unpublished error source is found. Finally, a specific Monte Carlo burnup procedure for actual MSR burnup calculations is developed by coupling KENO-VI with MODEC, and the burnup benchmark of Molten Salt Fast Reactor (MSFR) is calculated to validate the specific Monte Carlo burnup procedure. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction Molten Salt Reactor (MSR) is one of the six candidates of Gen-IV advanced nuclear power system (U.S. DOE, 2002). The fuel in the MSR is dissolved in molten salt carrier which acts as coolant as well. The molten salt mixture circulates through the primary loop continuously, producing heat in the core and removing the heat by the external heat exchangers. Because of the unique fluid fuel form, the composition of fuel salt can be readily altered by on-line fuel salt reprocessing and continuously refueling during operation. The detailed reprocessing diagram referred to MSBR and MSFR (Nuttin et al., 2005; Mathieu et al., 2009) is shown in Fig. 1. In this reprocessing scheme, the gaseous and noble metal fission products are removed by helium bubbling, and the soluble fission products are removed by on-line reductive extraction. Besides, for a thorium-based MSR, Pa extracted from the reactor core by ⇑ Corresponding authors at: Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China. E-mail addresses: [email protected] (J. Wu), [email protected] (X. Cai). https://doi.org/10.1016/j.anucene.2018.09.032 0306-4549/Ó 2018 Elsevier Ltd. All rights reserved.

on-line reprocessing is stored in the stockpile for some time to let 233Pa decay into 233U, and the decay product 233U is then reinjected into the core. Additionally, the fresh fuel should be fed continuously into the core to maintain critical operation. The on-line reprocessing and refueling system can bring good neutron economy, which makes it possible for Th/233U breeding in thermal MSR. However, it causes significant differences in burnup calculation from the traditional solid fueled reactor. The biggest difference in mathematics is that the characteristic of continuously refueling would bring a nonhomogeneous term into the conventional burnup equation. The depletion codes which can only solve the conventional burnup equations, such as WIMS (Lindley et al., 2017), DRAGON (Marleau et al., 2011) for simplified nuclides system, and CINDER (Wilson et al., 2008), DEPTH (She et al., 2013) for complicated nuclides system, cannot be applied directly to the MSR. ORIGEN2 (Croff, 1983) and ORIGEN-S (Gauld, 2011), which can solve the nonhomogeneous burnup equations, are suitable for MSR depletion calculation, while the depletion algorithm implemented in the two codes of truncation Taylor series expansion method with the separate treatment of short-lived

89

S. Xia et al. / Annals of Nuclear Energy 124 (2019) 88–97

(1) For the first node:

a1 ¼ 1 b1;1 ¼ 0 k~j ¼ k1 c1;1;0 ¼ n1 ð0Þ

ð3Þ

(2) For all j where ~ kj – kiþ1 :

biþ1;j ¼ bi;j kiþ1;j;k

ciþ1;j;biþ1;j ¼

kiþ1 ci;j;bi;j

kiþ1  kj kiþ1 ci;j;k  ðk þ 1Þciþ1;j;kþ1 ¼ kiþ1  cj



0 6 k 6 biþ1;j



ð4Þ

(3) If 9^j 2 f1; 2; . . . ; ai g; ~ k^j ¼ kiþ1 ,then:

aiþ1 ¼ ai

biþ1;^j ¼ bi;^j þ 1   ciþ1;^j;k ¼ 0 6 k 6 biþ1;^j k aiþ1 X ciþ1;^j;0 ¼ niþ1 ð0Þ  ciþ1;j;0 kiþ1;i ci;^j;k1

Fig. 1. Reprocessing diagram of MSR.

nuclides causes a loss of accuracy (Songqian et al., 2016). SERPENT2 (Aufiero et al., 2013) has been recently extended and employed to investigate burnup evolution of the Molten Salt Fast Reactor (MSFR) (Merle-Lucotte et al., 2011). However, the extended SERPENT-2 cannot directly solve the nonhomogeneous burnup equations but convert the nonhomogeneous term to the homogeneous term based on the conservation of the total molar fraction of heavy metals in the fuel salt. In order to address the unique burnup features of MSRs and meanwhile, maintain a high precision in depletion calculations, a specific depletion code MODEC for MSRs is developed. In MODEC, several advanced and high-precision algorithms for burnup equations and nonhomogeneous systems are employed, which are provided in detail in Sections 2 and 3. Then, its reliability is verified in Section 4, and finally, in Section 5, conclusions are given.

ð5Þ

j¼1;j–^j

Otherwise 8^j 2 f1; 2; . . . ; ai g; ~ k^j – kiþ1 , then:

aiþ1 ¼ ai þ 1 biþ1;aiþ1 ¼ 0 ~ka ¼ kiþ1 iþ1

ciþ1;aiþ1 ;0 ¼ niþ1 ð0Þ 

ai X

ciþ1;j;0

ð6Þ

j¼1;j–^j

To avoid the unnecessary long chains, the trajectory passage (Cetnar, 2006) is introduced to cut off the unimportant chains:

Pn ¼

niþ1 ðtÞ
ð7Þ

The second method implemented in MODEC is Chebyshev rational approximation method (CRAM) (Pusa and Leppänen, 2010), which is one of the matrix exponential methods. The burnup equations described by Eq. (1) can be written in a matrix form:

2. Mathematical model and method 2.1. Algorithms of conventional burnup equation The conventional burnup equations consist of the following set of first-order liner differential governing equations:

d~ n ¼ A~ n dt

8 X dni ðt Þ > kj;i nj  ki ni ; < dt ¼

where ~ n is the vector of nuclide concentrations, and A is the transition matrix containing the rate coefficients of nuclide transmutation by decay and/or neutron reaction. The solution of the matrixform burnup equations has an exponential form:

j

> : k ¼ f rtot / þ c kdecay ; j;i j!i j j!i j

decay ki ¼ rtot : i / þ ki

ð1Þ

where ni is the concentration of nuclide i; f j!i is the probability of a neutron reaction from nuclide j into nuclide i; rtot and rtot are the j i microscopic one-group total cross sections of nuclide j and nuclide i, respectively; / is the average neutron flux; cj!i is the branching ratio for the decay of nuclide j into nuclide i; cdecay is the decay coni stant of nuclide i. To solve the burnup system described by Eq. (1), two basic methods are implemented in MODEC. The first is transmutation trajectory analysis (TTA). The basic idea of TTA is to decompose the complex burnup web into a set of linear nuclide-chains, and then solve them one by one. A recursive form of the TTA (Huang et al., 2016) is implemented in MODEC. The solution of the ith node takes the following form:

ni ðt Þ ¼

bi;j ai X X

ci;j;k tk e~kj t

ð2Þ

j¼1 k¼0

And the recursion formulas of the associated variables are given as follows:

ð8Þ

~ nðtÞ ¼ eAt~ n0

ð9Þ

where eAt is the matrix exponential, and ~ n0 is the vector of initial nuclide concentrations. The CRAM method approximates the exponential by a rational function for the interval ð1; 0. And the rational function is then expressed in a pole-residue form:

e z  r k ð zÞ ¼

k=2 k X X P k ð zÞ ai ai ¼ a0 þ 2Re ¼ a0 þ Q k ð zÞ z  hi z  hi i¼1 i¼1

ð10Þ

where P k and Q k are polynomials of order k; ai and hi are the residues and poles. Since all eigenvalues of the matrix At are negative (Pusa and Leppänen, 2010), the matrix exponential eAt can be computed by CRAM. Thus, Eq. (9) can be written in CRAM as:

~ n0 nðt Þ ¼ eAt~ n0  ^r k;k ðAt Þ~ ¼ a0~ n0 þ 2Re

k=2 X ðhi I þ AtÞ1 ai~ n0 i¼1

!

ð11Þ

90

S. Xia et al. / Annals of Nuclear Energy 124 (2019) 88–97

where I is the identity matrix. The precision of the solution can be improved by increasing the order of approximation. In the MODEC code, a 16-order CRAM, which shows high enough accuracy in burnup system (Isotalo and Aarnio, 2011), is employed. 2.2. Treatment of on-line reprocessing

kci ¼

e

ekp t ¼ 1  e 1 kp ¼  lnð1  eÞ t

where e is extraction efficiency of nuclide i; T r is time required for reprocessing the entire fuel salt in the MSR system. For the TTA method, the effective decay constants ki in Eqs. (2)– (6) should be modified by adding the fictive decay constant: e i

c km i ¼ ki þ ki

ð13Þ

Regarding to the CRAM method, the matrix form of burnup equation can be modified as:

  d~ n  ¼ A  diag kci n dt

ð14Þ

  where diag kci is the diagonal matrix composed of fictive decay constants. 2.3. Treatment of continuously refueling The continuously refueling introduces a nonhomogeneous term into the burnup equation, as indicated in the following equation:

dni ðtÞ X ¼ kj;i nj  ki ni þ f i dt j

P

ð15Þ

ð16Þ



t ð1  eÞlnð1  eÞ

2.3.1. Treatment based on the TTA method For the TTA method, the treatment of the nonhomogeneous term is to represent the feed rate f i by the decay of a pseudo nuclide (Cetnar, 2006; Li et al., 2017):

ð17Þ

where kp;i is the transfer coefficient from the pseudo nuclide to the nuclide i; and np is the nuclide concentration of the pseudo nuclide. Assuming there are k nuclides being fed into the burnup system, kp;i can be expressed as:

ð22Þ

Utilizing Eqs. 17,18 and 21,22, the nonhomogeneous burnup equation in the form of Eq. (15) can be converted into the homogeneous burnup equation in the form of Eq. (1). So the TTA method described in Section 2.1 can be directly used without any modifications. 2.3.2. Treatment based on the CRAM method The matrix form of the nonhomogeneous burnup equations described by Eq. (16) cannot be directly computed by CRAM. In order to solve the nonhomogeneous burnup equations, two different methods are implemented in MODEC. The first method is based on the numerical integral. The solution of Eq. (16) has been proved to be the sum of the general solution of the related homogeneous equations and the particular integral (Ball and Adams, 1967), which can be expressed as:

~ n ðt Þ ¼ eAt~ n0 þ eAt

Z

t

f eAs ds~

ð23Þ

The analytical solution of the integral of matrix exponential in Eq. (23) cannot be given since A is a singular matrix. Thus, the numerical integral method must be used to compute the integral. The integral calculated by the composite trapezoidal quadrature formula is given by Li et al. (2017). The MODEC code implements Gauss-Legendre quadrature to calculate the integral for higher algebraic accuracy:

Z

where ~ f is the vector of feed rate.

f i ¼ kp;i np

np ð0Þ ¼ 

kf k

0

where f i is the feed rate of nuclide i. The matrix form of the nonhomogeneous burnup equations is described as:

d~ n ¼ A~ n þ~ f dt

ð21Þ

According to Eqs. (17)–(21), the initial concentration of the pseudo nuclide is calculated as:

ð12Þ

Tr

ð20Þ

and therefore, the decay constant of the pseudo nuclide is

A fictive decay constant method (Nuttin et al., 2005) is introduced for simulating the on-line reprocessing. The fictive decay constant is described as: e i

The item ekp t must be approach to 1 to ensure the feed rate almost unchanged. Introducing a positive constant e which is very close to 0, the item ekp t can be expressed as:

0

t

eAs ds ¼

n h i tX t xi  eA2ð1þxi Þ 2 i¼1

ð24Þ

where xi ranging from 1 to 1 is the quadrature point for n-order Gauss-Legendre quadrature; xi is the corresponding weight. Thus, Eq. (23) turns into:

nðt Þ ¼ eAt~ n0 þ

n   t X t  xi  eA2ð1xi Þ  ~f 2 i¼1

ð25Þ

where Pf i f represents the branching ratio of nuclide i; kp is the

Since ð1  xi ) is greater than 1, all the eigenvalues of matrix A 2t ð1  xi Þ are still strictly negative. Thus, the series of discrete matrix exponentials at all quadrature points can be accurately computed by CRAM. Another method is to convert the original nonhomogeneous burnup system to a homogeneous one by using the augment matrix (Isotalo and Wieselquist, 2015):

decay constant of the pseudo nuclide. Because the pseudo nuclide is the first-node nuclide without any precursors, the concentration of the pseudo nuclide can be easily obtained according to Eqs. (2) and (3), which is described as:

a1;1 6 . 6 .. e ~ A~ n þ~ f ¼ ^ A n~ ¼ 6 6 4 an;1

f kp;i ¼ P i

kf k

kp

ð18Þ

k k

np ðt Þ ¼ np ð0Þekp t

ð19Þ

2

0

   an;n

32 3 f1 n1 6 . 7 .. 7 6 . 7 . 7 76 . 7 76 7 5 f n 4 nn 5

0

0

   a1;n .. .. . . 0

d0

ð26Þ

S. Xia et al. / Annals of Nuclear Energy 124 (2019) 88–97

with the initial condition

~ ne0 ¼ ½ n1 ð0Þ n2 ð0Þ    nn ð0Þ 1 

T

where d0 is a pseudo nuclide for constructing the augment matrix. Thus, Eq. (16) can be turned into homogeneous:

e ~ e~ ~ n~ðtÞ ¼ A n~ ) ~ n~ðt Þ ¼ e At~ n0

ð27Þ

which can also be directly computed by CRAM. 2.4. Real-time tracing of nuclides evolution in the stockpile The nuclides extracted from the reactor core by on-line reprocessing are firstly stored in the stockpile for some time, and then a portion of the nuclides, such as 233U from decay of extracted 233 Pa, can be reinjected into the core, and the others, mostly fission products, will be further stored for the final geological disposal. Because of the continuously on-line reprocessing, tracing the decay evolution of the stockpiled nuclides must be in real time during the burnup calculation of the core. The evolution equations of the extracted nuclides can be described as:

X dni cj!i kj nout  ki nout þ kci nin ¼ j i i dt j out

ð28Þ

where nout is the concentration of extracted nuclide i; nin i i is the concentration of in-core nuclide i. It should be noticed that the last term of Eq. (28) represents the source from the in-core extraction. Thus, it is necessary to couple the burnup equations of in-core nuclides and the evolution equations of extracted nuclides. And the coupled equations can be written in matrix form:

"

~ nin d out ~ n dt

#

" ¼

# " #   ~ 0 A  diag kci nin 0    ~ D diag kci nout 0 )

ð29Þ

d~ n ¼ A~ n dt

where 0 is zero matrix; D is transition matrix containing only decay coefficients of the extracted nuclides. Eq. (29) has the same form as Eq. (14) and can be directly computed by CRAM method, and the block matrix method can be easily employed to improve computation efficiency. 3. MODEC features 3.1. Data library The data library adopted by the MODEC code is identical to that of ORIGEN-S developed by ORNL, which contains 1693 nonrepeated nuclides, 11 decay modes, 23 neutron reactions and 30 nuclides assigned mean fission yields. The MODEC code can convert the ORIGEN-S data to its own or directly read ORIGEN-S data library made by COUPLE (Gauld and Wiarda, 2011) which is a nuclear decay and cross-section data processing code used for creating ORIGEN-S libraries. Additionally, the MODEC code can also read or convert the data library of ORIGEN-2. It should be noticed that there are only 8 nuclides assigned fission yields in ORIGEN-2, thus the MODEC code adopts the nearest neighbor approximation method employed by ORIGEN-2 itself to treat the fissionable nuclides without explicit fission yields. Besides, recoverable energy is another important data which affects the accuracy of power-to-flux conversion. The MODEC code adopts the same recoverable energy data as ORIGEN-S. There are 24 fission nuclides and 32 important neutron-absorbing nuclides

91

which own specific values of the recoverable energy for fission and radiative capture based on ENDF/B data, while an average recoverable energy of 5 MeV per radiative neutron capture reaction and an average recoverable energy of 200 MeV per fission are applied for all other nuclides. 3.2. Functions and implementation MODEC is written in the C++ programming language and its CRAM solver is based on a homemade sparse matrix library, where the LU factorization and Gauss elimination are implemented to compute the matrix inverse. The MODEC code has three calculating modes including the decay mode, the constant flux mode and the constant power mode. MODEC can also calculate instantaneous physical quantities related to burnup calculations, such as radioactivity, decay heat and radiotoxicity. The same multiplication factor approximation method as ORIGEN-S is implemented in MODEC, which is expressed as (Lovecky` et al., 2014):

P P i Ni j mj ri;j k¼ P P N i i j ri;j

ð30Þ

where mj is neutron multiplication of the nuclear reaction j (0 for absorption reactions, x for ðn; xnÞ reactions, m for fission). Besides, MODEC uses Extensible Markup Language (XML) (Bray et al., 1997) for user input and output files. The XML-format file is easy to read and set by users. Meanwhile, the XML format makes it normalized to exchange data between MODEC and the Monte Carlo transport code for actual burnup calculations. Additionally, an alternative version of the TTA solver utilizing the MPFR multiple precision library (Fousse et al., 2007) is implemented in MODEC. This version of the TTA solver can be applied to obtain exact solutions without numerical precision loss under the preset cutoff, which would be used to analyze the errors of other depletion algorithms. 4. Validation and discussion In this section, three cases will be firstly implemented by comparing with ORIGEN-S to validate the correctness of MODEC in general burnup calculation, burnup calculation with on-line reprocessing and continuously refueling, respectively. All three cases are tested in the constant flux mode. Then, a detailed analysis of the error sources in ORIGEN-S will be shown in Section 4.4. Finally, a special Monte Carlo burnup procedure is established and the performance of the actual burnup calculation for MSR is illustrated. It should be noted that unless explicitly stated, the ORIGEN-S applied in this paper is from the SCALE6.1 code package (Bowman, 2011). 4.1. Case 1: depletion of Thorium fuel salt The first case was depletion calculation of Thorium fuel salt. The compositions of the selected Thorium fuel salt are similar to that of Molten Salt Fast Reactor (MSFR): 77.5% mol LiF – 22.5% mol ThF4. It is irradiated by 12 months through 10 steps with constant flux 1:0  1014 cm2  s1 . To avoid the precision errors from data converting, the MODEC code reads the data library generated by COUPLE code directly. A cutoff density of 1:0  1020 atom=ðbarn  cm Þ is applied to the solution. Table 1 gives the concentrations of several main nuclides at the end of the irradiation. The differences in the table are between ORIGEN-S and MODEC-CRAM. It is shown that the relative errors for both fission products and actinides are generally below 1.0%, which indicates the good agreement between ORIGEN-S and

92

S. Xia et al. / Annals of Nuclear Energy 124 (2019) 88–97

Table 1 Calculation results of case 1. Nuclide Concentrations (mol) Nuclide

79

Se Sr 91 Y 95 Zr 95 Nb 103 Ru 103m Rh 131 I 135 Xe 137 Cs 140 La 144 Ce 143 Pr 147 Pm 151 Sm 156 Eu 232 Th 233 Th 233 Pa 232 U 233 U 234 U 235 U 236 U 237 Np 238 Pu 239 Pu 241 Am 91

MODEC CRAM

TTA(1020 )

2.9794E05 4.4464E06 5.0116E04 5.1810E04 2.3170E04 8.5918E05 8.4297E08 4.6758E05 1.5971E06 1.4094E03 1.7549E05 7.9747E04 1.3037E04 3.0517E04 2.6358E05 6.4236E07 2.2278E+01 1.3409E05 2.3172E02 3.4531E05 1.7352E01 3.8715E03 1.1111E04 1.2914E06 6.7020E09 1.3445E10 3.3281E12 3.3754E17

2.9794E05 4.4464E06 5.0116E04 5.1810E04 2.3170E04 8.5918E05 8.4297E08 4.6758E05 1.5971E06 1.4094E03 1.7549E05 7.9747E04 1.3037E04 3.0517E04 2.6358E05 6.4236E07 2.2278E+01 1.3409E05 2.3172E02 3.4531E05 1.7352E01 3.8715E03 1.1111E04 1.2914E06 6.7020E09 1.3445E10 3.3280E12 3.3165E17

ORIGEN-S

Diff.

2.9802E05 4.4458E06 5.0510E04 5.1815E04 2.3175E04 8.5928E05 8.4401E08 4.6726E05 1.6000E06 1.4097E03 1.7551E05 7.9762E04 1.3037E04 3.0526E04 2.6363E05 6.4249E07 2.2278E+01 1.3409E05 2.3165E02 3.4536E05 1.7354E01 3.8767E03 1.1134E04 1.2951E06 6.7260E09 1.3503E10 3.3442E12 3.3980E17

0.025% 0.013% 0.786% 0.011% 0.022% 0.012% 0.122% 0.069% 0.183% 0.022% 0.011% 0.018% 0.004% 0.028% 0.020% 0.020% 0.000% 0.000% 0.028% 0.015% 0.010% 0.135% 0.212% 0.286% 0.358% 0.427% 0.485% 0.671%

MODEC. Additionally, it is seen from the table that the results from the two methods implemented in MODEC are almost the same, which further validates the correctness of MODEC by selfverification. 4.2. Case 2: depletion of Thorium fuel salt with on-line reprocessing An extended case 2 is designed to validate the method of treating on-line reprocessing: during the period of depletion of Thorium fuel salt in case 1, 233Pa from irradiated 232Th is on-line extracted with a fictive decay constant of 1:15741  106 s1 and then stored until it decays to 233U. The time evolutions of 233Pa nuclide density and its daughter nuclide 233U are shown in Fig. 1 and Fig. 2(b), respectively. The results of the CRAM method and of the TTA method are completely the same. Therefore, only one set of the results from MODEC is shown in each figure. It can be seen that the MODEC code gives the consistent results with OIRGEN-S. Additionally, because ORIGEN-S is not able to trace the evolution of nuclides extracted by on-line reprocessing in real time, an indirect way is taken to validate the results of the extracted nuclides computed by the MODEC code. In this case, the nuclide densities of extracted 233Pa and its daughter nuclide 233U in the stockpile are calculated. The governing equation of evolution of 233 Pa in the stockpile can be expressed as: out

dnp3 c in ¼ kp3 nout p3 þ kp3 np3 dt

ð31Þ

where the last term of Eq. (31) is the external source of the stockpiled 233Pa. Thus, the nuclide density of 233Pa at time t1 can be obtained analytically by Eq. (31), and the expression of the result is shown as:

Fig. 2. The time evolution of (a)

233

Pa density and (b)

233

U density.

93

S. Xia et al. / Annals of Nuclear Energy 124 (2019) 88–97

out kp3 t1 nout þ kcp3 ekp3 t1 p3 ðt 1 Þ ¼ np3 ðt 0 Þ  e

Z

Besides, the nuclide density of obtained:

Z c nout u3 ðt 1 Þ ¼ kp3

t1

t0

t1

kp3 s nin ds p3 ðsÞ  e

t0

ð32Þ

233

U at time t 1 can be also

Table 3 Calculation results of case 3. Nuclide

MODEC (mol)

ORIGEN-S (mol)

Diff.

79

1.8022E06 1.0090E06 5.1027E05 6.4678E10 2.1237E15 1.6876E06 9.4756E05 4.5843E07 6.2577E05 1.7874E05 2.7848E06 7.1646E08 2.2304E+01 1.3425E05 4.7724E03 1.7216E06 3.8721E02 9.0398E04 2.7586E05 3.4045E07 1.8749E09 3.9590E11 1.0207E12 1.1582E17

1.8022E06 1.0088E06 5.1451E05 6.4685E10 2.1239E15 1.6927E06 9.4769E05 4.5870E07 6.2582E05 1.7875E05 2.7850E06 7.1654E08 2.2304E+01 1.3425E05 4.7696E03 1.7279E06 3.8724E02 9.0782E04 2.7734E05 3.4255E07 1.8878E09 3.9884E11 1.0287E12 1.1689E17

0.000% 0.017% 0.831% 0.011% 0.009% 0.306% 0.013% 0.059% 0.008% 0.005% 0.006% 0.011% 0.000% 0.000% 0.058% 0.361% 0.009% 0.424% 0.534% 0.619% 0.688% 0.743% 0.787% 0.917%

Se Sr 91 Y 95 Zr 99 Tc 133 I 137 Cs 141 La 144 Ce 147 Pm 151 Sm 156 Eu 232 Th 233 Th 233 Pa 232 U 233 U 234 U 235 U 236 U 237 Np 238 Pu 239 Pu 241 Am 91

out nin p3 ðsÞds  np3 ðt 1 Þ

ð33Þ

It should be emphasized that Eq. (33) is held based on the assumption that 233U is a stable nuclide, and the error of the assumption is negligible because the evolution time is much shorter compared to the half-life of 233U in this case. The integrals in Eq. (32) and in Eq. (33) are calculated by numerical quadrature. The time point of each sub step is selected as quadrature point, thus nin p3 ðsÞ can be easily obtained. Three sets of results based on trapezoidal quadrature with 10, 100, 1000 quadrature points are calculated, and the last set is chosen as the reference. As can be seen in Table 2, the direct results calculated by the MODEC code agree well with the reference, which demonstrates the validity of the MODEC code in tracing the evolution of nuclides extracted by on-line reprocessing.

4.3. Case 3: depletion of Thorium fuel salt with continuously refueling Test case 3 is a further extension from case 2: the Thorium fuel salt with the composition mentioned in case 1 is continuously fed

4.4. Error analysis of ORIGEN-S

with a constant of 1:0  109 mol  s1 feed rate during the irradiated period, which means the feed rate of 7Li, 19F and 232Th are

From the above three cases, the correctness of MODEC has been validated by comparing with ORIGEN-S. Some studies (Songqian et al., 2016; Isotalo and Aarnio, 2011) have pointed out that ORIGEN-S have relatively low accuracy and are sensitive to the time interval in burnup calculations. According to these studies, the cause of the relatively low accuracy in ORIGEN-S is simply attributed to the treatment of short-lived nuclides and is not fully discussed yet. Therefore, a further study of error sources in ORIGEN-S has been implemented and will be discussed in detail in this section. Fig. 4, where each point represents a nuclide, shows the relative errors of ORIGEN-S compared to the MPFR-version TTA solver in MODEC following a sub-step of 0.1 year and 0.01 year in the case

2:897  1010 mol  s1 ; 6:262  1010 mol  s1 and 8:411  1011 mol  s1 , respectively. Similar to case 1, the results of the concentrations of several main nuclides are presented in Table 3. The nuclide densities calculated by MODEC are based on augment matrix method of CRAM. Still, it can be seen that MODEC and ORIGEN-S agree well in both fission products and actinides. Besides, the three methods of treating constant feed rate mentioned in Section 2.3 are compared. Fig. 3 shows the relative errors between these methods. The reference results are computed by the augment matrix method of CRAM. For the 30-order GaussLegendre quadrature method, the relative errors of most nuclides are below 109 , and the maximum relative error is about 7

5  10 . For the pseudo nuclide method of TTA, the relative errors increase with the reduction of concentrations of nuclides, which is due to the numerical errors and the chain cutoff (Isotalo and Aarnio, 2011). However, for most nuclides, the pseudo nuclide method is precise enough. Table 4 gives the total computing time in this case of the three methods executed on the same personal computer. Augmented matrix method of CRAM is highly efficient. While the numerical integral method of CRAM is much slower because a series of matrix exponential at quadrature points must be computed by CRAM. Thanks to the efficient CRAM solver, the efficiencies of the both methods are acceptable. However, the pseudo nuclide method of TTA has the worst computing efficiency, which is mainly because searching and calculating the transmutation chains one by one will be very time consuming (Isotalo and Aarnio, 2011).

Table 2 Comparisons of nuclide densities of

233 233

Pa U

233

Pa and

233

1. Cutoff 1030 and 256-bit precision are used in the TTA solver, which ensures at least nine correct digits for nuclides with final concentration fraction of 1020 . The horizontal axis of the chart is the ratio of the nuclide effective half-life to the burnup sub-step time interval, which will be called ”Ratio” for short thereafter. And for more clear explanation about the impact for the nuclides with different effective half-lives, three zones, which is shown in Fig. 4, are divided according to the nuclide effective half-lives.

4.4.1. Error distribution of nuclides in ORIGEN-S In ORIGEN-S, Ratio is the criterion for considering whether the nuclides are short-lived or not. The nuclides with Ratio less than 0.1 are regarded as short-lived, and the other nuclides are regarded as long-lived. So, the line Ratio ¼ 0:1 in Fig. 4 is the dividing line between the long-lived nuclides and the short-lived nuclides.

U in the stockpile at 1 year.

References

MODEC

100 sub steps

10 sub step

1.8569E02 1.5257E01

0.00073% 0.00020%

0.073% 0.051%

7.16% 4.30%

94

S. Xia et al. / Annals of Nuclear Energy 124 (2019) 88–97

Fig. 3. Relative errors of nuclide densities between augmented matrix method of CRAM and the other two methods. (a) 30-order Gauss-Legendre quadrature method of CRAM. (b) pseudo nuclide method of TTA (1.0E-20).

Table 4 Computing time of the three methods.

Computing time (s)

Method I⁄

Method II⁄⁄

Method III⁄⁄⁄

0.336

3.658

 5  103

⁄

Method I is the augmented matrix method of CRAM. Method II is the 30-order Gauss quadrature of CRAM. ⁄⁄⁄ Method III is the pseudo nuclide method of TTA.

Fig. 4. Relative errors of nuclide densities between ORIGEN-S and MODEC in (a) 0.1-year case and (b) 0.01-year case, where ‘Ratio’ is the ratio of the nuclide effective half-life to the burnup sub-step time interval.

⁄⁄

It can be seen from Fig. 4 that the distribution regularities of the errors are significantly different between two sides of the line Ratio ¼ 0:1, which is because that ORIGEN-S implements two different depletion methods for the long-lived nuclides and for the short-lived nuclides, respectively. In ORIGEN-S, the power series approximation of the matrix exponential is applied to the longlived nuclides. While for the short-lived nuclides, the generalized solution form of Bateman equations with secular equilibrium assumptions is utilized. And secular equilibrium assumptions would become more precise as the effective half-lives of the nuclides decrease. So, the relative errors of the short-lived nuclides are negatively correlated to the effective half-lives, which is shown in Fig. 4. Additionally, when the time interval of one sub-step is altered, the categories of short-lived or long-lived of some nuclides would also be changed in ORIGEN-S because of the invariable criterion. For example, the nuclides in zone II of Fig. 4 are treated as the

short-lived nuclides in the 0.1-year case and then converted to the long-lived nuclides in the 0.01-year case. For the long-lived nuclides, the relative errors are distributed almost stochastically, which is shown in zone III, for instance. Besides, the mean and the deviation of relative errors decrease with the reduction of the time interval. It is because that the new transition matrix elements, which are determined for the new transitions after removing the short-lived nuclides, asymptotically approach the correct value with time as successive time intervals are executed. Therefore, the results in the 0.01-year case are more precise because of the more time intervals. Moreover, for the short-lived nuclides, the relative errors are positively correlated with the effective half-lives, which is indicated in zone I, for instance. It is because that secular equilibrium assumptions, which are the main error sources for the short-lived nuclides, would become more real and therefore more accurate as the effective half-lives of the nuclides decrease. Meanwhile, the long-lived nuclides would propagate errors to the short-lived nuclides through secular equilibrium assumptions. So, when the time interval is reduced, the relative errors of some short-lived

95

S. Xia et al. / Annals of Nuclear Energy 124 (2019) 88–97

nuclides would decrease due to the more accurate results of their long-lived precursors. And, it can be seen that the relative errors of most short-lived nuclides are non time-sensitive, which is mainly because that the relative errors of their long-lived precursors, for example 232 Th in these cases, are non time-sensitive. 4.4.2. Incomplete nuclides classification in ORIGEN-S It is also noticed that some nuclides with over 10% and even 100% relative errors are negligibly affected by the reduction of sub-step time interval. These nuclides are listed in Table 5. As will be explained later, it is mostly due to the incomplete nuclides classification in ORIGEN-S. In ORIGEN-S, the number of nuclides included in the libraries is expanded to 910 activation nuclides, 176 actinides, and 1151 fission products. The activation nuclides construct a burnup space, and the actinides together with the fission products construct the other burnup space. These two burnup spaces are calculated separately and independently. Considering the total number of nonrepeated nuclides in ORIGEN-S is about 1700, hundreds of nuclides must be repeatedly classified and included in the both burnup spaces. However, the precursors and the daughters of the repeatedly-classified nuclides may be included in only one burnup space, which would cause the incomplete burnup chains for these nuclides and therefore bring computing errors in ORIGEN-S. For instance, Fig. 5 shows the main transmutation paths for 133 Ba and 133mBa. In ORIGEN-S, 133Ba is considered as both the activation nuclide and the fission product, and 133mBa, which can decay to 133Ba, is only classified as the activation nuclides. Their common parent nuclides 132Ba and 134Ba, which are considered as both the activation nuclide and the fission product, however are produced only as the fission products from 232Th in the both subcases. Since 133mBa is only the activation nuclide, the paths from the fission-product 132Ba and 134Ba to 133mBa are ignored by ORIGEN-S and so the concentration of the 133mBa is always the incorrect 0. Furthermore, the contribution from 133mBa to 133Ba is therefore 0 as well, which causes the incorrect lower value of 133 Ba compared to the reference concentration. To validate the above analysis, the three reaction channels of Baðn; cÞ133m Ba, 132mBaðn; cÞ133 Ba and 134Baðn; 2nÞ133m Ba are artificially closed in the MODEC code and the recalculated results of 133 Ba and 133mBa are shown in Table 6. It can be seen that the recalculated results are very close to the results from ORIGEN-S, which indirectly validates that the incomplete burnup chain is the main error source of 133Ba and 133mBa in ORIGEN-S. Other nuclides listed in Table 5 have the same error reasons as 133 Ba and 133mBa. And some of them have relatively large nuclide 132

densities (> 1:015 ), which cannot be easily ignored.

Fig. 5. The main production and decay paths for 133Ba and 133mBa in the burnup system, where ‘‘Burnup Space I” contains all the actinides and fission products and ‘‘Burnup Space II” contains all the activation nuclides.

Table 6 Recalculated results of

133

Ba and

133m

Ba.

Nuclides

Recalculated results from MODEC

0.01-year-case results from ORIGEN-S

133

2.3408E13 0.0(4.0706E22)

2.3429E13 0.0

Ba Ba

133m

oped by coupling KENO-VI with MODEC. The calculation procedure starts from the preparation of the core geometry and of the isotopic composition of the fuel salt. And in each time step, the feed rate of fuel salt is searched based on secant method to ensure the criticality of the MSR. If the keff does not meet the requirement, the feed rate is updated, and the time step is repeated. Once the keff is within the pre-set margins, the calculation results of the time step are saved, and the next time step is started. The main procedure described above is shown in Fig. 6. In this subsection, Molten Salt Fast Reactor (MSFR) is chosen as the actual burnup calculation test case (Brovchenko et al., 2013). The time evolutions up to equilibrium of main actinides are shown in Fig. 7. In the figure, the results of reference were from LPSC/ IN2P3/CNRS (Heuer et al., 2014) relying on coupling MCNP (Brown et al., 2003) with a home-made materials evolution code REM (Nuttin et al., 2005). As can be seen in the figure, the results simulated by the specific procedure are in good agreement with that of the reference both for 233U-started and TRU-started MSFR cases, which validates the correctness of the specific procedure for actual burnup calculation of MSR.

4.5. Burnup calculation of a reference MSR concept

5. Conclusion

To validate the MODEC code performance in actual burnup calculation of MSR, a specific Monte Carlo burnup procedure is devel-

In this work, a Molten Salt Reactor specific depletion code – MODEC is developed. The two depletion algorithms of TTA and

Table 5 List of the nuclides with large and not-time-sensitive relative errors in ORIGEN-S. Nuclide

Reference

133

Ba Ba 205 Tl 131

133m

Ba,150Eu,

150

62

63

63

64

Gd,173Yb

Ni, Ni, Ni, Cu, 64Zn,65Zn,150mEu, 173 Lu,175Lu,204Hg

173

Tm,174Yb,175Yb,

0.01-year-case results from ORIGEN-S

0.1-year-case results from ORIGEN-S

Concentrations

Relative errors

Concentrations

Relative errors

3.4625E13 3.8162E19 1.8851E18

2.3429E13 1.7233E19 5.3430E22

32.34% 54.84% 99.97%

2.3716E13 1.7294E19 5.5613E22

31.51% 54.33% 99.97%

>1.0E15

0.0

100%

0.0

100%

<1.0E15

0.0

100%

0.0

100%

96

S. Xia et al. / Annals of Nuclear Energy 124 (2019) 88–97

Fig. 6. Flow chart of the specific Monte Carlo burnup procedure for MSRs.

CRAM are implemented in the MODEC code. Then the fictive decay constant method is introduced in the MODEC code to simulate online reprocessing. And the continuously refueling characteristic introduces the nonhomogeneous term into the conventional burnup equations. To solve the nonhomogeneous burnup equations, the MODEC code implemented three methods: the pseudo nuclide method of TTA, the numerical integral method of CRAM and the augmented matrix method of CRAM. Besides, the evolution of nuclides extracted by on-line reprocessing and stored in the stockpile can be traced in real time during the burnup calculation of the core. Then, the correctness of the MODEC code is validated through the three cases by comparing with ORIGEN-S. Moreover, the evolution of the on-line reprocessed nuclides is validated by an indirect approach, and the three methods of solving nonhomogeneous burnup equations are validated and compared. In addition, a detailed error analysis of ORIGEN-S is implemented. The separate treatment of short-lived nuclides and the incomplete nuclides classification are found out to be the two main error sources of ORIGEN-S. And the results calculated from ORIGEN-S are sensitive to the burnup time interval, especially for the long-lived nuclides and a part of the short-lived nuclides that would be converted to long-lived nuclides with the decrease of the burnup time interval. In SCALE6.2, the CRAM method has already been applied in the ORIGEN-S code (Isotalo and Wieselquist, 2015), therefore the short-lived nuclides no longer need to be treated separately. However, it seems that the incomplete nuclides classification in the latest ORIGEN-S have not been amended yet. Moreover, a specific Monte Carlo burnup procedure for the Molten Salt Reactor based on coupling MODEC with KENO-VI is developed, and the MSFR test case is calculated to validate MODEC code in actual burnup calculation of MSR.

Fig. 7. Comparison of time evolution of the actinides elements between the reference and KENOVI-MODEC for (a) 233U started and (b) TRU started.

Acknowledgements This work is supported by the Chinese TMSR Strategic Pioneer Science and Technology Project under Grant No. XDA02010000, the National Natural Science Foundation of China under Grant No. 91326201 and the Frontier Science Key Program of the Chinese Academy of Sciences under Grant No. QYZDY-SSW-JSC016. References Aufiero, M., Cammi, A., Fiorina, C., Leppänen, J., Luzzi, L., Ricotti, M., 2013. An extended version of the serpent-2 code to investigate fuel burn-up and core material evolution of the molten salt fast reactor. J. Nucl. Mater. 441 (1–3), 473– 486. Ball, S., Adams, R., 1967. matexp, a general purpose digital computer program for solving ordinary differential equations by the matrix exponential method (Tech. rep.). Oak Ridge National Lab., Tenn. Bowman, S.M., 2011. Scale 6: comprehensive nuclear safety analysis code system. Nucl. Technol. 174 (2), 126–148. Bray, T., Paoli, J., Sperberg-McQueen, C.M., Maler, E., Yergeau, F., 1997. Extensible markup language (xml). World Wide Web J. 2 (4), 27–66. Brovchenko, M., Merle-Lucotte, E., Rouch, H., Alcaro, F., Allibert, M., Aufiero, M., Cammi, A., Dulla, S., Feynberg, O., Frima, L., et al., 2013. Optimization of the preconceptual design of the msfr. Brown, F.B. et al., 2003. Mcnp – A General Monte Carlo N-particle Transport Code, Version 5. Los Alamos National Laboratory, Oak Ridge, TN. Cetnar, J., 2006. General solution of bateman equations for nuclear transmutations. Ann. Nucl. Energy 33 (7), 640–645. Croff, A.G., 1983. Origen2: a versatile computer code for calculating the nuclide compositions and characteristics of nuclear materials. Nucl. Technol. 62 (3), 335–352.

S. Xia et al. / Annals of Nuclear Energy 124 (2019) 88–97 Fousse, L., Hanrot, G., Lefèvre, V., Pélissier, P., Zimmermann, P., 2007. Mpfr: a multiple-precision binary floating-point library with correct rounding. ACM Trans. Math. Software 33 (2), 13. Gauld, I., 2011. Origen-s: Depletion Module to Calculate Neutron Activation, Actinide Transmutation, Fission Product Generation, and Radiation Source Terms. Oak Ridge National Laboratory, Oak Ridge, Tennessee. 37831-36170. Gauld, I., Wiarda, D., 2011. Couple: A Nuclear Decay and Cross Section Data Processing Code for Creating Origen-s Libraries. Oak Ridge National Laboratory, Oak Ridge. Heuer, D., Merle-Lucotte, E., Allibert, M., Brovchenko, M., Ghetta, V., Rubiolo, P., 2014. Towards the thorium fuel cycle with molten salt fast reactors. Ann. Nucl. Energy 64, 421–429. Huang, K., Wu, H., Cao, L., Li, Y., Shen, W., 2016. Improvements to the transmutation trajectory analysis of depletion evaluation. Ann. Nucl. Energy 87, 637–647. Isotalo, A., Aarnio, P., 2011. Comparison of depletion algorithms for large systems of nuclides. Ann. Nucl. Energy 38 (2–3), 261–268. Isotalo, A.E., Wieselquist, W.A., 2015. A method for including external feed in depletion calculations with cram and implementation into origen. Ann. Nucl. Energy 85, 68–77. Li, J., She, D., Shi, L., 2017. Algorithms of solving the burnup equation with external feed. Ann. Nucl. Energy 110, 412–417. Lindley, B., Hosking, J., Smith, P., Powney, D., Tollit, B., Newton, T., Perry, R., Ware, T., Smith, P., 2017. Current status of the reactor physics code wims and recent developments. Ann. Nucl. Energy 102, 148–157.

97

Lovecky`, M., Piterka, L., Prehradny`, J., Škoda, R., 2014. Uwb1-fast nuclear fuel depletion code. Ann. Nucl. Energy 71, 333–339. Marleau, G., Hébert, A., Roy, R., 2011. A user guide for dragon version 4. IGE-294, Ecole Polytechnique de Montréal, Institut de génie nucléaire Département de génie mécanique. Mathieu, L., Heuer, D., Merle-Lucotte, E., Brissot, R., Le Brun, C., Liatard, E., Loiseaux, J.-M., Méplan, O., Nuttin, A., Lecarpentier, D., 2009. Possible configurations for the thorium molten salt reactor and advantages of the fast nonmoderated version. Nucl. Sci. Eng. 161 (1), 78–89. Merle-Lucotte, E., Heuer, D., Allibert, M., et al., 2011. The thorium molten salt reactor: launching the thorium fuel cycle with the molten salt fast reactor. Proc. ICAPP, 842. Nuttin, A., Heuer, D., Billebaud, A., Brissot, R., Le Brun, C., Liatard, E., Loiseaux, J.-M., Mathieu, L., Meplan, O., Merle-Lucotte, E., et al., 2005. Potential of thorium molten salt reactorsdetailed calculations and concept evolution with a view to large scale energy production. Prog. Nucl. Energy 46 (1), 77–99. Pusa, M., Leppänen, J., 2010. Computing the matrix exponential in burnup calculations. Nucl. Sci. Eng. 164 (2), 140–150. She, D., Wang, K., Yu, G., 2013. Development of the point-depletion code depth. Nucl. Eng. Des. 258, 235–240. Songqian, T., Feng, X., Yi, T., Shuping, W., 2016. Analysis of complex nuclide chain equation calculation based on different method. Nucl. Power Eng. 37 (1), 8–12. U.S. DOE, 2002. A technology roadmap for generation IV nuclear energy systems. Wilson, W., Cowell, S., England, T., Hayes, A., Moller, P., 2008. A Manual for Cinder 90 Version 07.4 Codes and Data. Los Alamos National Laboratory.