- Email: [email protected]
An alternative stochastic formulation for the point reactor
Annals of Nuclear Energy 63 (2014) 691–695
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
An alternative stochastic formulation for the point reactor Seyed Mohsen Ayyoubzadeh ⇑, Naser Vosoughi Department of Energy Engineering, Sharif University of Technology, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 23 August 2013 Received in revised form 3 September 2013 Accepted 4 September 2013
Keywords: Stochastic differential equation Point kinetics Neutron density variance Precursor concentration variance
a b s t r a c t The stochastic behavior of a point reactor is modeled with a system of Ito stochastic differential equations. This new approach does not require computing the square root of a matrix which is a great computational advantage. Moreover, the derivation procedure clearly demonstrates the mathematical approximations involved in the final formulation. Three numerical benchmarks show the accuracy of this model in predicting the mean and variance of the neutron and precursor population in a point reactor. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction It is well known that the reactions in a nuclear system are not fully describable by deterministic laws. This fact, at the most fundamental level, is due to the laws of quantum mechanics, which only give probabilities of various interactions for a neutron, which are manifest in the interaction cross sections of atoms with neutrons. There are various situations which this probabilistic behavior could be readily observed for a nuclear system, e.g. in the startup of the reactors, in zero power reactors, in most laboratory source – detector configurations, etc. There has been an extensive research effort to model this stochastic behavior. Some researchers have attempted to give a description based on the quantum mechanical laws of evolution in such a system (Osborn, 1969; Osborn and Yip, 1966). Others have tried to obtain suitable master equations for the system under study, which embed the stochastic information of the system (Degweker, 1994; Lewins, 1978; Pázsit and Pál, 2008). Some have tried to attribute this stochastic behavior to an (imaginary) fluctuation source, and have developed the so-called Langevinian description (Akcasu and Stolle, 1989; Gençay and Akcasu, 1981; Williams, 1974). These models generally trade simplicity for accuracy. However, the most accurate formulation of the stochastic processes involved may not always be useful. Measuring higher order moments requires more data from the system for a given accuracy. Actually, most of the times in practice, one only measures the first and second order moments in a system, i.e. the mean and variance.
⇑ Corresponding author. Tel.: +98 21 66166130; fax: +98 21 66081723. E-mail address: [email protected] (S.M. Ayyoubzadeh). 0306-4549/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2013.09.013
The Central Limit Theorem (CLT) assures the accuracy of such an examination from a mathematical perspective (Eagleson, 1975). This fact has manifest itself in recent stochastic research trends, namely the Stochastic Differential Equation (SDE) formulation (Allen, 2007; Ikeda and Watanabe, 1989; Øksendal, 2003). Formulations based on this method are capable of describing stochastic processes up to the second moment, which is enough in most applications. Most of these formulations have two main disadvantages. First, they may not be quite satisfactory from a physicist’s perspective, in the rather sudden appearance of Weiner processes. Second, the need to compute the square root of the covariance matrix in these formulations makes these undesirable from a computational point of view. This research describes an alternative SDE developed for a point reactor system, namely the Simplified Stochastic Point Kinetics equation (SSPK), based on (Gillespie, 2000) methodology. The current research is organized in six sections. In the next section, the stochastic processes involved in a point reactor are described. In the third section, the methodology for obtaining an SDE, assuming that the reactions of the system are known, is described. The fourth section deals with applying this methodology to a point reactor system with multiple precursor groups. In the fifth section, numerical benchmarks are given which show the accuracy of this method. Finally, the results are discussed and further improvements are suggested.
2. The stochastic point reactor A point reactor is a reactor in which the spatial effects have been eliminated. This is obviously possible if the reactors length is infinite in all spatial dimensions. However definitions exist
692
S.M. Ayyoubzadeh, N. Vosoughi / Annals of Nuclear Energy 63 (2014) 691–695
which eliminate the need for such an (unphysical) requirement (Henry, 1975). Fortunately, the final form in all of these formulations is alike. Namely, for a reactor with J precursor groups, one may write
2
J X
2
3 6 q bj n 6 j¼1 6 c 7 6 6 K 17 d6 6 . 7¼6 6 b1 6 7 dt 4 .. 5 6 K 6 . 6 . cJ 4 . bJ
K
k1 k1
3 3 2 3 72 n q 7 7 6 . . . kJ 7 c 7 6 0 7 17 6 7 76 6 7 76 .. 7 7 þ 6 .. 7; 76 4 5 4.5 7 . 7 5 cJ 0
ð1Þ
0
where in Eq. (1), n and cj are the neutron and jth precursor group populations, q is the strength of an external neutron source, q is the reactors reactivity, bj and kj are the delayed neutron fission fraction and decay constant for the jth precursor group, K is the prompt neutron lifetime (Hetrick, 1993), and the bar indicates the ensemble averaging operator. Study of a point reactor, i.e. studying the properties of Eq. (1), is desirable in the sense that it captures some of the most essential features of the reactor dynamics without involving into the complexities of integro-differential equations, i.e. the transport equation, or partial differential equations, i.e. the diffusion equation. However, as it is evident from Eq. (1), this formulation does not contain any stochastic terms, hence it could not describe the fluctuations in neutron and precursor populations. As a matter of fact, this equation should be interpreted as the ensemble averaged quantities in infinite alike reactors (Bell, 1969). The starting point in accurately modeling the stochastic processes in a point reactor is the consideration of basic processes involved in a small time step dt. These processes are listed in Table 1. In this table, v is the mean velocity of neutrons, Rc is the probability per unit length of neutron capture and Rf,i1, i2, . . . , iJ+1 is the probability per unit length of a neutron to induce a fission which yields i1 prompt neutrons and ij+1 precursors for the jth precursor group. Using the data embedded in Table 1, the stochastic point reactor could be readily simulated by the Monte Carlo (MC) method, which the results may then be used to extract the stochastic properties of such a reactor. Unfortunately in this method, one relies heavily on simulations and could gain little information a priori on the stochastic properties. In other words, to obtain the moments evolution, even the mean behavior, needs a full MC simulation. A circumvent to these drawbacks is using an SDE model. Usually, in this method, one uses the model (Allen, 2007) 1
dx ¼ mdt þ C 2 dw;
ð2Þ
where in Eq. (2), w is a J + 1 dimensional Wiener process, m is the mean vector, and C is the covariance matrix, for reasons which are seen in the following. From the properties of Wiener processes,z namely
EðdxÞ ¼ mdt;
ð4Þ
EðdxdxT Þ ¼ Cdt:
Eq. (4) readily justifies the terms chosen for m and C. Using the data in Table 1, one could use Eq. (4) to readily compute m and C, as in (Allen, 2007; Hayes and Allen, 2005; Saha Ray, 2012). Such a computation leads to (Hayes and Allen, 2005)
2
3
J X
3
2
6 q bj n 6 j¼1 6 7 6 K d 6 c1 7 6 b1 6 .. 7 ¼ 6 dt 4 . 5 6 6 K 6 .. 4 . cJ bJ
0
K
2
f 6 a1 6 þ 6 .. 4 .
a1 r1
aJ
bJ;2
f ¼ cn þ
and by noting the self-adjointness of the covariance matrix, it could be seen that
7 7 7 dw; 5
ð5Þ
J X kj cj þ q; j¼1
1 q þ 2
c¼
J P
bj þ 1
j¼1
J P
!2 bj
m
j¼1
;
K ! ! J P bi 1 þ 1 bj m j¼1
ai ¼
ð6Þ n ki ci ;
K
bi1 bj1 m n; K b2 m r i ¼ i n þ ki c i : K bi;j ¼
While this model takes account of the stochastic process involved, using a Wiener process term, routine application of this model is difficult due to the need to compute the square root of the covariance matrix in each step. Note that this process is computationally expensive and may result in non-real values for the physical quantities due to numerical round offs. Also, the reason for introducing the Wiener processes in this model may not be readily clear and be a source for confusion. In the next section we aim at eliminating these drawbacks. 3. An alternative SDE modeling Assume that the stochastic quantities in a system are shown by the vector x. Assuming the system to follow a Markov model, for the ith component of such a system, one may write
X
Rr ðxðt0 Þ; sÞmr;i ;
ð7Þ
r2All reactions
ð3Þ
EðdwdwT Þ ¼ Idt;
312
. . . aJ
where in Eq. (5)
xi ðt0 þ sÞ ¼ xi ðt 0 Þ þ
EðdwÞ ¼ 0;
k1 k1
72 n 3 2 3 q 7 607 . . . kJ 7 c1 7 76 7 6 6 76 . 7 þ 6 . 7 74 . 5 4 .. 7 5 7 . 7 0 5 cJ
where in Eq. (7), we have represented the total number of r-type reactions starting from the state x(t0) after a time s by Rr(x(t0), s). Also, the number of i-type particles generated from a r-type
Table 1 Basic processes in a point reactor in an infinitesimal time interval. Probability rate of occurrence
Change in neutron population
Change in jth precursor group population
Radiative capture Fission that results in i1 neutrons and i2, . . . , iJ+1 precursors
nvRc nv Rf;i1 ;i2 ;...;iJþ1
1 i1 1
0 iJ+1
Decay for a precursor in group j External source emission
c j kj q
+1 +1
1 0
693
S.M. Ayyoubzadeh, N. Vosoughi / Annals of Nuclear Energy 63 (2014) 691–695
reaction is denoted by mr,i. Note that obtaining the evolution kernel, Rr(x(t0), s), for an arbitrary time interval requires detailed knowledge of the values of x in that interval. However, if s is chosen to be sufficiently small, such that the population change in this time interval is negligible, the probability rate of an r-type reaction could be assumed to be constant which we may represent by kr ðxðt0 ÞÞ. Hence Rr(x(t0), s) would follow a Poisson distribution, and one may write
Note that in Eq. (16), there is no need for operation on matrices. Also, the linear increment of the deterministic part and the superlinear increment, namely square root, of the stochastic part with time have appeared naturally as a consequence of the addition theorem for stable distributions.
xi ðt 0 þ sÞ ¼ xi ðt0 Þ þ
Using the foundation laid in the previous section, we shall derive the stochastic equations of evolution in a point reactor in this section. For obtaining the stochastic dynamics equations in a point reactor, one may replace the parameters in Eq. (16) with the data in Table 1.0Such replacement yields 1
X
Pr ðkr ðxðt 0 ÞÞ; sÞmr;i ;
ð8Þ
r2All reactions
where we have shown the Poisson type random variable for the rth reaction in Eq. (8) by Pr . Note that the probability for a change of m particles according to this distribution is equal to
PðPr ðkr ðxðt 0 ÞÞ; sÞ ¼ mÞ ¼
ðkr ðxðt0 ÞÞsÞm ekr ðxðt0 ÞÞs : m!
ð9Þ
4. The alternative SDE for a point reactor
Eq. (10) shows that the probability of a change of m particles in the population is zero for all values except those sufficiently close to the mean, namely those which may be expressed as (Evans, 2002)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m ¼ kr ðxðt 0 ÞÞs þ e kr ðxðt 0 ÞÞs;
j
ð17Þ and
0
dcj ¼ @nv 0
e
ð12Þ
Note that on may further use Eq. (11) into Eq. (12) to eliminate e. With such a replacement one may arrive at
ðm kr ðxðt 0 ÞÞsÞ2 2kr ðxðt 0 ÞÞ e PðBr ðkr ðxðt 0 ÞÞ; sÞ ¼ mÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2pkr ðxðt 0 ÞÞs
xi ðt 0 þ sÞ ¼ xi ðt0 Þ þ
Nr ðkr ðxðt 0 ÞÞs; kr ðxðt0 ÞÞsÞmr;i :
ð18Þ Note that the complete stochastic data for Rf,i1, i2, . . . , iJ+1 is usually not present. In this case, one assumes that a fission always results in ! X bj m neutrons and bjm precursors jth precursor group. 1 j
Hence, Eqs. (17) and (18) could be rewritten as
ð14Þ
ð15Þ
where nr is a standard Normal variable. Replacing Eq. (15) into Eq. (14) and replacing s by dt, yields
dxi ¼
X
! kr ðxÞmr;i dt þ
r2All reactions
X
! pffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffi kr ðxÞmr;i nr dt : ð16Þ
r2All reactions
X 1 bj
!
!
! X m 1 þ cj kj þ q dt
j
ð13Þ
One may note that in this stage, x has lost the property of being an integer and has transformed into a real-valued variable. This is not a source of concern if the population is large enough, as assumed formerly. To proceed, properties of Normal distributions should be exploited. One core property of such distributions is that they belong to the class of stable distributions, i.e. distributions in which linear combinations belong to the same distribution. Hence, one may write (Zolotarev, 1986)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kr ðxðt 0 ÞÞsnr ;
ijþ1 Rf;i1 ;i2 ;...;iJþ1 cj kj Adt
i1 ;i2 ;...;iJþ1
dn ¼ nv Rc þ nv Rf
r2All reactions
Nr ðkr ðxðt 0 ÞÞs; kr ðxðt0 ÞÞsÞ ¼ kr ðxðt 0 ÞÞs þ
1
i1 ;i2 ;...;iJþ1
X bj 1
pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffi þ nv Rc nc þ nv Rf
Note that the R.H.S. of Eq. (13) represents the probability distribution function for a Normally distributed random variable whose mean and variance are equal to kr ðxðt 0 ÞÞs. Hence Eq. (8) could be re-written as
X
X
1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffi pffiffiffiffiffiffiffiffi pffiffiffiffiffiffi X ijþ1 Rf;i1 ;i2 ;...;iJþ1 nf;i1 ;i2 ;...;iJþ1 cj kj nd;j A dt : þ @ nv
where in Eq. (11), e is a sufficiently small constant. Replacing Eq. (11) into (10) and using the smallness assumption enforced on e yields
e 2 PðBr ðkr ðxðt 0 ÞÞ; sÞ ¼ mÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2pkr ðxðt 0 ÞÞs
i1 ;i2 ;...;iJþ1
Xpffiffiffiffiffiffiffiffi pffiffiffiffiffi pffiffiffi cj kj nd;j þ qnq þ dt;
ð11Þ
2
j
!
:
ð10Þ
X cj kj þ qAdt
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi X þ @ nv Rc nc þ nv ði1 1Þ Rf;i1 ;i2 ;...;iJþ1 nf;i1 ;i2 ;...;iJþ1
; small m ; large m
ði1 1ÞRf;i1 ;i2 ;...;iJþ1 þ
i1 ;i2 ;...;iJþ1
0
Now, if one assumes kr ðxðt0 ÞÞs to be sufficiently large, the term ekr ðxðt0 ÞÞs makes the probability of small changes negligibly small, hence one may consider Eq. (9) only for large m. Using Stirling’s approximation one may write
8 <0 m kr ðxðt ÞÞs 0 PðPr ðkr ðxðt 0 ÞÞ; sÞ ¼ mÞ k ðxðt 0 ÞÞse e : r pffiffiffiffiffiffiffiffiffiffiffi m 2p m
X
dn ¼@nv Rc þ nv
!
j
!
m 1 nf
ð19Þ
j
! Xpffiffiffiffiffiffiffiffi pffiffiffiffiffi pffiffiffi cj kj nd;j þ qnq þ dt; j
and
pffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi pffiffiffiffiffi dcj ¼ nv Rf bj m cj kj dt þ nv Rf bj mnf cj kj nd;j dt :
ð20Þ
The obtained model shown by Eqs. (19) and (20) could be readily simulated. The only disadvantage of this model is that it requires J + 3 Wiener processes in comparison with the J + 1 processes needed in Eq. (2). However, further simplification is possible by noting that seldom the covariance is needed in practical applications, so that one may use the addition theorem for independent Normally distributed random variables, namely qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi
aN1 ð0; 1Þ þ bN2 ð0; 1Þ ¼
2
a2 þ b N3 ð0; 1Þ:
Using Eq. (21) into Eq. (19) yields ! dn ¼
nv Rc þ nv Rf
1
X bj j
!
m1 þ
ð21Þ
! X cj kj þ q dt j
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! !2 u u X X pffiffiffiffiffi þ tnv Rc þ nv Rf 1 bj m 1 þ cj kj þ qnn dt : j
j
ð22Þ
694
S.M. Ayyoubzadeh, N. Vosoughi / Annals of Nuclear Energy 63 (2014) 691–695
Also, use of Eq. (21) into Eq. (20) yields
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffi dcj ¼ ðnv Rf bj m cj kj Þdt þ nv Rf ðbj mÞ2 þ cj kj ncj dt :
ð23Þ
Note that use of Eq. (21) has coupled the random variables nn and ncj in a way that they are no longer independent. Hence, Eqs. (22) and (23) should only be used when correlations between n and cj are not of concern. Note that
k1 ¼ mRRaf ;
q ¼ 1 k11 ; K¼
1 v Ra k 1
ð24Þ
:
Using Eq. (24)into Eqs. (22) and (23), and casting the resulting equations into a matrix form yields
2 6 6 6 6c 7 6 6 17 6 7 6 d6 6 .. 7 ¼ 6 4 . 5 6 6 6 cJ 4 2
n
3
q
3 2 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X uX u2 bj þð1 bj Þ2 mð1þqÞ 7 6u bj X 7 6u j X j 7 7 6t 0 ... 0 nþ kj c j þ q nþ kj cj þ q 7 K 7 6 7 7 6 j 7 j 7 6 7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 7 7dt þ 6 b1 2 7dw 6 b1 m n k c 7 1 1 7 6 K n þ k 0 1 c1 7 K 7 6 7 .. 7 6 7 . 7 6 . .. 5 7 6 bJ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 5 4 n kJ cJ b2J m K n þ kJ cJ 0 0 K 3
X j K
As the first benchmark, we shall consider the one group delayed neutron point kinetics model of a TRIGA reactor. The parameters of this reactor in clean startup conditions, i.e. n(0) = c(0) = 0, are k ¼ 0:077 ðs1 Þ, b = 0.0079, m = 2.432, K = 0.0001(s), q = 0.003 (Hetrick, 1993) and an external neutron source of q = 10,000 (s1) has been assumed (Allen, 2007). This reactor has been simulated for 0.1 (s) with 106 computational time steps for 10,000 realizations. The mean and standard deviation of the neutron and precursor populations using SSPK have been compared with the exact MC simulation and Stochastic PCA (Allen, 2007) in Table 2. Also the time behavior of the mean neutron population is compared with a sample in Fig. 1. Note that while the population of neutrons and precursors are zero at startup, the final results are close to the MC simulations. This indicates that while the requirements in Section 3 are sufficient for a rigorous mathematical derivation
ð25Þ
we shall refer to Eq. (25) by SSPK equations. Table 3 Computational times for the first benchmark problem.
5. Numerical results Here, we shall measure the accuracy of SSPK by comparing the results of Eq. (25) with published benchmarks for the mean and standard deviations in the point kinetics equation.
Time (s)
MC
PCA
SSPK
Deterministic
6.71
414.34
5.15
4.23
Table 2 Mean and standard deviation of neutron and precursor population for the first benchmark problem.
E[n(0.1)] r[n(0.1)] E[c(0.1)] r[c(0.1)]
Monte Carlo
Stochastic PCA
SSPK
199.15 152.63 1254.5 613.94
204.52 174.03 1294.0 620.68
208.14 174.30 1293.22 622.12
Fig. 2. Time distribution to reach a neutron level of 4000.
Table 4 The mean and standard deviation of time needed to reach 4000 neutrons.
Fig. 1. Sample and average neutron population for first benchmark problem.
E[t](s) r[t](s)
Monte Carlo
Stochastic PCA
SSPK
33.136 2.0886
33.157 2.5772
33.125 2.6731
695
S.M. Ayyoubzadeh, N. Vosoughi / Annals of Nuclear Energy 63 (2014) 691–695 Table 5 The mean and standard deviation for the third benchmark problem.
E[n(0.1)]
r[n(0.1)]
P E½ 6j¼1 cj ð0:1Þ P r½ 6j¼1 cj ð0:1Þ
Monte Carlo
Stochastic PCA
Euler–Maruyama approximation
SSPK
183.04 168.79 4.478 105
186.31 164.16 4.491 105
208.6 255.95 4.498 105
184.8 186.96 4.489 105
1495.7
1917.2
1233.38
982.64
Table 6 Computational times for the third benchmark problem.
Time (s)
MC
PCA
SSPK
Deterministic
7.34
4610.4
11.50
4.67
of a system of SDEs, they may not always be necessary. Obtaining the full range of validity for SSPK is beyond the scope of this paper. In Table 3, we have compared the computational time needed to obtain a realization for this benchmark problem using the stochastic equations with the time needed to solve the forward finite differenced, deterministic PK equations. The implementation in (Hayes, 2005) has been used for the PCA method and the computations have been performed on a 2.2 GHz Core i7 CPU. In the second example, we examine the accuracy of the model by calculating the distribution of the time necessary for a reactor to reach a predefined neutron population. We shall use the (non-physical) parameters used by (Hayes and Allen, 2005) as b = 0.05, k ¼ 0:1 ðs1 Þ, m = 2.5, q = 200 (s1), K = 0.499002 (s), and q = 0.001996 for a clean reactor. The distribution of the time to reach 4000 neutrons is shown in Fig. 2. The mean and standard deviation of this distribution for 5000 realizations is obtained using SSPK and compared with the values obtained from MC and Stochastic PCA (Hayes and Allen, 2005) methods in Table 4. In the last benchmark problem, we consider the 6-group delayed neutron model of a pool-type light-water research reactor (Chao and Attard, 1985; Kinard and Allen, 2004; Saha Ray, 2012; Yang and Jevremovic´, 2009). The parameters of this reactor are that of the VR-1 research reactor and equal to b = [266, 1491, 1316, 2849, 896, 182] 106, k ¼ ½127; 317; 1150; 3110; 14000; 38700 104 ðs1 Þ, m = 2.5, q = 0, K = 2 105 (s), and q = 0.003 (Matejka et al., 2003). The equilibrium conditions are considered at the time origin, namely ½ nð0Þ c1 ð0Þ . . . c6 ð0Þ ¼ 100 ½ 1 kb11K . . . kb66K . The result has been compared at t = 0.1(s) for 5000 realizations with (Saha Ray, 2012) in Table 5. The computational times required for obtaining a realization using the MC, PCA, SSPK and deterministic methods for this reactor have been summarized in Table 6. As it is seen, eliminating the computation of the square root of a matrix has reduced the computational time required for SSPK to solve this problem nearly 400 times in comparison with PCA (Hayes and Allen, 2005). 6. Conclusion In this paper, an alternative derivation of the stochastic differential equations that describe a point reactor has been presented. A few approximations have been used to simplify the final form
of the equations. The validity of the resultant equations, namely the Simplified Stochastic Point Kinetics equations (SSPK), has been confirmed by numerical benchmarks which show good agreement with published results. Possible future works include improving the current work by advising advanced numerical SDE algorithms and generalization of these equations to include spatial and multigroup neutron effects. References Akcasu, A.Z., Stolle, A., 1989. Comments on the noise equivalent source in the Langevin technique. Annals of Nuclear Energy 16, 493–508. Allen, E., 2007. Modeling with Itô Stochastic Differential Equations. Springer. Bell, G.I., 1969. Stochastic formulations of neutron transport. SIAM-AMS Proceedings. American Mathematical Society, 181–197. Chao, Y.-A., Attard, A., 1985. A resolution of the stiffness problem of reactor kinetics. Nuclear Science and Engineering; (United States) 90. Degweker, S., 1994. A forward equation for stochastic neutron transport. Annals of Nuclear Energy 21, 531–539. Eagleson, G., 1975. On Gordin’s central limit theorem for stationary processes. Journal of Applied Probability, 176–179. Evans, L.C., 2002. An Introduction to Stochastic Differential Equations. UC Berkley. Gençay, Sß., Akcasu, Z., 1981. Noise analysis of a zero power critical line reactor using Langevin’s technique. Annals of Nuclear Energy 8, 65–77. Gillespie, D.T., 2000. The chemical Langevin equation. The Journal of Chemical Physics 113, 297. Hayes, J.G., 2005. Stochastic Point Kinetics Equations in Nuclear Reactor Dynamics, Department of Mathematics. Texas Tech University. Hayes, J.G., Allen, E., 2005. Stochastic point-kinetics equations in nuclear reactor dynamics. Annals of Nuclear Energy 32, 572–587. Henry, A.F., 1975. Nuclear-Reactor Analysis. MIT press Cambridge. Hetrick, D.L., 1993. Dynamics of Nuclear Reactors. American Nuclear Society (La Grange Park, Ill., USA). Ikeda, N., Watanabe, S., 1989. Stochastic Differential Equations and Diffusion Processes. Kinard, M., Allen, E., 2004. Efficient numerical solution of the point kinetics equations in nuclear reactor dynamics. Annals of Nuclear Energy 31, 1039– 1051. Lewins, J., 1978. Linear stochastic neutron transport theory. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 362, 537–558. Matejka, K., Kolros, A., Kropik, M., Sklenka, L., 2003. Eugene Wigner Training Course at VR-1 Reactor. Department of Nuclear Reactors, CTU in Prague. Øksendal, B., 2003. Stochastic Differential Equations. Springer. Osborn, R., 1969. Kinetic theory of transport and fluctuation phenomena. SIAM-AMS Proceedings. American Mathematical Society, p. 213. Osborn, R.K., Yip, S., 1966. The foundations of neutron transport theory. In: Osborn, Richard K., (Ed.), The foundations of neutron transport theory. Yip, Sidney. New York, Gordon and Breach [1966]. Monograph Series on Nuclear, Science and Technology 1. Pázsit, I., Pál, L., 2008. Neutron Fluctuations: A Treatise on the Physics on Branching Processes. Elsevier Science Limited. Saha Ray, S., 2012. Numerical simulation of stochastic point kinetic equation in the dynamical system of nuclear reactor. Annals of Nuclear Energy 49, 154–159. Williams, M.M.R., 1974. Random Processes in Nuclear Reactors. Yang, X., Jevremovic´, T., 2009. Revisiting the Rosenbrock numerical solutions of the reactor point kinetics equation with numerous examples. Nuclear Technology and Radiation Protection 24, 3–12. Zolotarev, V.M., 1986. One-dimensional stable distributions. American Mathematical Soc.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
An alternative stochastic formulation for the point reactor Seyed Mohsen Ayyoubzadeh ⇑, Naser Vosoughi Department of Energy Engineering, Sharif University of Technology, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 23 August 2013 Received in revised form 3 September 2013 Accepted 4 September 2013
Keywords: Stochastic differential equation Point kinetics Neutron density variance Precursor concentration variance
a b s t r a c t The stochastic behavior of a point reactor is modeled with a system of Ito stochastic differential equations. This new approach does not require computing the square root of a matrix which is a great computational advantage. Moreover, the derivation procedure clearly demonstrates the mathematical approximations involved in the final formulation. Three numerical benchmarks show the accuracy of this model in predicting the mean and variance of the neutron and precursor population in a point reactor. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction It is well known that the reactions in a nuclear system are not fully describable by deterministic laws. This fact, at the most fundamental level, is due to the laws of quantum mechanics, which only give probabilities of various interactions for a neutron, which are manifest in the interaction cross sections of atoms with neutrons. There are various situations which this probabilistic behavior could be readily observed for a nuclear system, e.g. in the startup of the reactors, in zero power reactors, in most laboratory source – detector configurations, etc. There has been an extensive research effort to model this stochastic behavior. Some researchers have attempted to give a description based on the quantum mechanical laws of evolution in such a system (Osborn, 1969; Osborn and Yip, 1966). Others have tried to obtain suitable master equations for the system under study, which embed the stochastic information of the system (Degweker, 1994; Lewins, 1978; Pázsit and Pál, 2008). Some have tried to attribute this stochastic behavior to an (imaginary) fluctuation source, and have developed the so-called Langevinian description (Akcasu and Stolle, 1989; Gençay and Akcasu, 1981; Williams, 1974). These models generally trade simplicity for accuracy. However, the most accurate formulation of the stochastic processes involved may not always be useful. Measuring higher order moments requires more data from the system for a given accuracy. Actually, most of the times in practice, one only measures the first and second order moments in a system, i.e. the mean and variance.
⇑ Corresponding author. Tel.: +98 21 66166130; fax: +98 21 66081723. E-mail address: [email protected] (S.M. Ayyoubzadeh). 0306-4549/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2013.09.013
The Central Limit Theorem (CLT) assures the accuracy of such an examination from a mathematical perspective (Eagleson, 1975). This fact has manifest itself in recent stochastic research trends, namely the Stochastic Differential Equation (SDE) formulation (Allen, 2007; Ikeda and Watanabe, 1989; Øksendal, 2003). Formulations based on this method are capable of describing stochastic processes up to the second moment, which is enough in most applications. Most of these formulations have two main disadvantages. First, they may not be quite satisfactory from a physicist’s perspective, in the rather sudden appearance of Weiner processes. Second, the need to compute the square root of the covariance matrix in these formulations makes these undesirable from a computational point of view. This research describes an alternative SDE developed for a point reactor system, namely the Simplified Stochastic Point Kinetics equation (SSPK), based on (Gillespie, 2000) methodology. The current research is organized in six sections. In the next section, the stochastic processes involved in a point reactor are described. In the third section, the methodology for obtaining an SDE, assuming that the reactions of the system are known, is described. The fourth section deals with applying this methodology to a point reactor system with multiple precursor groups. In the fifth section, numerical benchmarks are given which show the accuracy of this method. Finally, the results are discussed and further improvements are suggested.
2. The stochastic point reactor A point reactor is a reactor in which the spatial effects have been eliminated. This is obviously possible if the reactors length is infinite in all spatial dimensions. However definitions exist
692
S.M. Ayyoubzadeh, N. Vosoughi / Annals of Nuclear Energy 63 (2014) 691–695
which eliminate the need for such an (unphysical) requirement (Henry, 1975). Fortunately, the final form in all of these formulations is alike. Namely, for a reactor with J precursor groups, one may write
2
J X
2
3 6 q bj n 6 j¼1 6 c 7 6 6 K 17 d6 6 . 7¼6 6 b1 6 7 dt 4 .. 5 6 K 6 . 6 . cJ 4 . bJ
K
k1 k1
3 3 2 3 72 n q 7 7 6 . . . kJ 7 c 7 6 0 7 17 6 7 76 6 7 76 .. 7 7 þ 6 .. 7; 76 4 5 4.5 7 . 7 5 cJ 0
ð1Þ
0
where in Eq. (1), n and cj are the neutron and jth precursor group populations, q is the strength of an external neutron source, q is the reactors reactivity, bj and kj are the delayed neutron fission fraction and decay constant for the jth precursor group, K is the prompt neutron lifetime (Hetrick, 1993), and the bar indicates the ensemble averaging operator. Study of a point reactor, i.e. studying the properties of Eq. (1), is desirable in the sense that it captures some of the most essential features of the reactor dynamics without involving into the complexities of integro-differential equations, i.e. the transport equation, or partial differential equations, i.e. the diffusion equation. However, as it is evident from Eq. (1), this formulation does not contain any stochastic terms, hence it could not describe the fluctuations in neutron and precursor populations. As a matter of fact, this equation should be interpreted as the ensemble averaged quantities in infinite alike reactors (Bell, 1969). The starting point in accurately modeling the stochastic processes in a point reactor is the consideration of basic processes involved in a small time step dt. These processes are listed in Table 1. In this table, v is the mean velocity of neutrons, Rc is the probability per unit length of neutron capture and Rf,i1, i2, . . . , iJ+1 is the probability per unit length of a neutron to induce a fission which yields i1 prompt neutrons and ij+1 precursors for the jth precursor group. Using the data embedded in Table 1, the stochastic point reactor could be readily simulated by the Monte Carlo (MC) method, which the results may then be used to extract the stochastic properties of such a reactor. Unfortunately in this method, one relies heavily on simulations and could gain little information a priori on the stochastic properties. In other words, to obtain the moments evolution, even the mean behavior, needs a full MC simulation. A circumvent to these drawbacks is using an SDE model. Usually, in this method, one uses the model (Allen, 2007) 1
dx ¼ mdt þ C 2 dw;
ð2Þ
where in Eq. (2), w is a J + 1 dimensional Wiener process, m is the mean vector, and C is the covariance matrix, for reasons which are seen in the following. From the properties of Wiener processes,z namely
EðdxÞ ¼ mdt;
ð4Þ
EðdxdxT Þ ¼ Cdt:
Eq. (4) readily justifies the terms chosen for m and C. Using the data in Table 1, one could use Eq. (4) to readily compute m and C, as in (Allen, 2007; Hayes and Allen, 2005; Saha Ray, 2012). Such a computation leads to (Hayes and Allen, 2005)
2
3
J X
3
2
6 q bj n 6 j¼1 6 7 6 K d 6 c1 7 6 b1 6 .. 7 ¼ 6 dt 4 . 5 6 6 K 6 .. 4 . cJ bJ
0
K
2
f 6 a1 6 þ 6 .. 4 .
a1 r1
aJ
bJ;2
f ¼ cn þ
and by noting the self-adjointness of the covariance matrix, it could be seen that
7 7 7 dw; 5
ð5Þ
J X kj cj þ q; j¼1
1 q þ 2
c¼
J P
bj þ 1
j¼1
J P
!2 bj
m
j¼1
;
K ! ! J P bi 1 þ 1 bj m j¼1
ai ¼
ð6Þ n ki ci ;
K
bi1 bj1 m n; K b2 m r i ¼ i n þ ki c i : K bi;j ¼
While this model takes account of the stochastic process involved, using a Wiener process term, routine application of this model is difficult due to the need to compute the square root of the covariance matrix in each step. Note that this process is computationally expensive and may result in non-real values for the physical quantities due to numerical round offs. Also, the reason for introducing the Wiener processes in this model may not be readily clear and be a source for confusion. In the next section we aim at eliminating these drawbacks. 3. An alternative SDE modeling Assume that the stochastic quantities in a system are shown by the vector x. Assuming the system to follow a Markov model, for the ith component of such a system, one may write
X
Rr ðxðt0 Þ; sÞmr;i ;
ð7Þ
r2All reactions
ð3Þ
EðdwdwT Þ ¼ Idt;
312
. . . aJ
where in Eq. (5)
xi ðt0 þ sÞ ¼ xi ðt 0 Þ þ
EðdwÞ ¼ 0;
k1 k1
72 n 3 2 3 q 7 607 . . . kJ 7 c1 7 76 7 6 6 76 . 7 þ 6 . 7 74 . 5 4 .. 7 5 7 . 7 0 5 cJ
where in Eq. (7), we have represented the total number of r-type reactions starting from the state x(t0) after a time s by Rr(x(t0), s). Also, the number of i-type particles generated from a r-type
Table 1 Basic processes in a point reactor in an infinitesimal time interval. Probability rate of occurrence
Change in neutron population
Change in jth precursor group population
Radiative capture Fission that results in i1 neutrons and i2, . . . , iJ+1 precursors
nvRc nv Rf;i1 ;i2 ;...;iJþ1
1 i1 1
0 iJ+1
Decay for a precursor in group j External source emission
c j kj q
+1 +1
1 0
693
S.M. Ayyoubzadeh, N. Vosoughi / Annals of Nuclear Energy 63 (2014) 691–695
reaction is denoted by mr,i. Note that obtaining the evolution kernel, Rr(x(t0), s), for an arbitrary time interval requires detailed knowledge of the values of x in that interval. However, if s is chosen to be sufficiently small, such that the population change in this time interval is negligible, the probability rate of an r-type reaction could be assumed to be constant which we may represent by kr ðxðt0 ÞÞ. Hence Rr(x(t0), s) would follow a Poisson distribution, and one may write
Note that in Eq. (16), there is no need for operation on matrices. Also, the linear increment of the deterministic part and the superlinear increment, namely square root, of the stochastic part with time have appeared naturally as a consequence of the addition theorem for stable distributions.
xi ðt 0 þ sÞ ¼ xi ðt0 Þ þ
Using the foundation laid in the previous section, we shall derive the stochastic equations of evolution in a point reactor in this section. For obtaining the stochastic dynamics equations in a point reactor, one may replace the parameters in Eq. (16) with the data in Table 1.0Such replacement yields 1
X
Pr ðkr ðxðt 0 ÞÞ; sÞmr;i ;
ð8Þ
r2All reactions
where we have shown the Poisson type random variable for the rth reaction in Eq. (8) by Pr . Note that the probability for a change of m particles according to this distribution is equal to
PðPr ðkr ðxðt 0 ÞÞ; sÞ ¼ mÞ ¼
ðkr ðxðt0 ÞÞsÞm ekr ðxðt0 ÞÞs : m!
ð9Þ
4. The alternative SDE for a point reactor
Eq. (10) shows that the probability of a change of m particles in the population is zero for all values except those sufficiently close to the mean, namely those which may be expressed as (Evans, 2002)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m ¼ kr ðxðt 0 ÞÞs þ e kr ðxðt 0 ÞÞs;
j
ð17Þ and
0
dcj ¼ @nv 0
e
ð12Þ
Note that on may further use Eq. (11) into Eq. (12) to eliminate e. With such a replacement one may arrive at
ðm kr ðxðt 0 ÞÞsÞ2 2kr ðxðt 0 ÞÞ e PðBr ðkr ðxðt 0 ÞÞ; sÞ ¼ mÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2pkr ðxðt 0 ÞÞs
xi ðt 0 þ sÞ ¼ xi ðt0 Þ þ
Nr ðkr ðxðt 0 ÞÞs; kr ðxðt0 ÞÞsÞmr;i :
ð18Þ Note that the complete stochastic data for Rf,i1, i2, . . . , iJ+1 is usually not present. In this case, one assumes that a fission always results in ! X bj m neutrons and bjm precursors jth precursor group. 1 j
Hence, Eqs. (17) and (18) could be rewritten as
ð14Þ
ð15Þ
where nr is a standard Normal variable. Replacing Eq. (15) into Eq. (14) and replacing s by dt, yields
dxi ¼
X
! kr ðxÞmr;i dt þ
r2All reactions
X
! pffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffi kr ðxÞmr;i nr dt : ð16Þ
r2All reactions
X 1 bj
!
!
! X m 1 þ cj kj þ q dt
j
ð13Þ
One may note that in this stage, x has lost the property of being an integer and has transformed into a real-valued variable. This is not a source of concern if the population is large enough, as assumed formerly. To proceed, properties of Normal distributions should be exploited. One core property of such distributions is that they belong to the class of stable distributions, i.e. distributions in which linear combinations belong to the same distribution. Hence, one may write (Zolotarev, 1986)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kr ðxðt 0 ÞÞsnr ;
ijþ1 Rf;i1 ;i2 ;...;iJþ1 cj kj Adt
i1 ;i2 ;...;iJþ1
dn ¼ nv Rc þ nv Rf
r2All reactions
Nr ðkr ðxðt 0 ÞÞs; kr ðxðt0 ÞÞsÞ ¼ kr ðxðt 0 ÞÞs þ
1
i1 ;i2 ;...;iJþ1
X bj 1
pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffi þ nv Rc nc þ nv Rf
Note that the R.H.S. of Eq. (13) represents the probability distribution function for a Normally distributed random variable whose mean and variance are equal to kr ðxðt 0 ÞÞs. Hence Eq. (8) could be re-written as
X
X
1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffi pffiffiffiffiffiffiffiffi pffiffiffiffiffiffi X ijþ1 Rf;i1 ;i2 ;...;iJþ1 nf;i1 ;i2 ;...;iJþ1 cj kj nd;j A dt : þ @ nv
where in Eq. (11), e is a sufficiently small constant. Replacing Eq. (11) into (10) and using the smallness assumption enforced on e yields
e 2 PðBr ðkr ðxðt 0 ÞÞ; sÞ ¼ mÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2pkr ðxðt 0 ÞÞs
i1 ;i2 ;...;iJþ1
Xpffiffiffiffiffiffiffiffi pffiffiffiffiffi pffiffiffi cj kj nd;j þ qnq þ dt;
ð11Þ
2
j
!
:
ð10Þ
X cj kj þ qAdt
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi X þ @ nv Rc nc þ nv ði1 1Þ Rf;i1 ;i2 ;...;iJþ1 nf;i1 ;i2 ;...;iJþ1
; small m ; large m
ði1 1ÞRf;i1 ;i2 ;...;iJþ1 þ
i1 ;i2 ;...;iJþ1
0
Now, if one assumes kr ðxðt0 ÞÞs to be sufficiently large, the term ekr ðxðt0 ÞÞs makes the probability of small changes negligibly small, hence one may consider Eq. (9) only for large m. Using Stirling’s approximation one may write
8 <0 m kr ðxðt ÞÞs 0 PðPr ðkr ðxðt 0 ÞÞ; sÞ ¼ mÞ k ðxðt 0 ÞÞse e : r pffiffiffiffiffiffiffiffiffiffiffi m 2p m
X
dn ¼@nv Rc þ nv
!
j
!
m 1 nf
ð19Þ
j
! Xpffiffiffiffiffiffiffiffi pffiffiffiffiffi pffiffiffi cj kj nd;j þ qnq þ dt; j
and
pffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi pffiffiffiffiffi dcj ¼ nv Rf bj m cj kj dt þ nv Rf bj mnf cj kj nd;j dt :
ð20Þ
The obtained model shown by Eqs. (19) and (20) could be readily simulated. The only disadvantage of this model is that it requires J + 3 Wiener processes in comparison with the J + 1 processes needed in Eq. (2). However, further simplification is possible by noting that seldom the covariance is needed in practical applications, so that one may use the addition theorem for independent Normally distributed random variables, namely qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi
aN1 ð0; 1Þ þ bN2 ð0; 1Þ ¼
2
a2 þ b N3 ð0; 1Þ:
Using Eq. (21) into Eq. (19) yields ! dn ¼
nv Rc þ nv Rf
1
X bj j
!
m1 þ
ð21Þ
! X cj kj þ q dt j
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! !2 u u X X pffiffiffiffiffi þ tnv Rc þ nv Rf 1 bj m 1 þ cj kj þ qnn dt : j
j
ð22Þ
694
S.M. Ayyoubzadeh, N. Vosoughi / Annals of Nuclear Energy 63 (2014) 691–695
Also, use of Eq. (21) into Eq. (20) yields
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffi dcj ¼ ðnv Rf bj m cj kj Þdt þ nv Rf ðbj mÞ2 þ cj kj ncj dt :
ð23Þ
Note that use of Eq. (21) has coupled the random variables nn and ncj in a way that they are no longer independent. Hence, Eqs. (22) and (23) should only be used when correlations between n and cj are not of concern. Note that
k1 ¼ mRRaf ;
q ¼ 1 k11 ; K¼
1 v Ra k 1
ð24Þ
:
Using Eq. (24)into Eqs. (22) and (23), and casting the resulting equations into a matrix form yields
2 6 6 6 6c 7 6 6 17 6 7 6 d6 6 .. 7 ¼ 6 4 . 5 6 6 6 cJ 4 2
n
3
q
3 2 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X uX u2 bj þð1 bj Þ2 mð1þqÞ 7 6u bj X 7 6u j X j 7 7 6t 0 ... 0 nþ kj c j þ q nþ kj cj þ q 7 K 7 6 7 7 6 j 7 j 7 6 7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 7 7dt þ 6 b1 2 7dw 6 b1 m n k c 7 1 1 7 6 K n þ k 0 1 c1 7 K 7 6 7 .. 7 6 7 . 7 6 . .. 5 7 6 bJ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 5 4 n kJ cJ b2J m K n þ kJ cJ 0 0 K 3
X j K
As the first benchmark, we shall consider the one group delayed neutron point kinetics model of a TRIGA reactor. The parameters of this reactor in clean startup conditions, i.e. n(0) = c(0) = 0, are k ¼ 0:077 ðs1 Þ, b = 0.0079, m = 2.432, K = 0.0001(s), q = 0.003 (Hetrick, 1993) and an external neutron source of q = 10,000 (s1) has been assumed (Allen, 2007). This reactor has been simulated for 0.1 (s) with 106 computational time steps for 10,000 realizations. The mean and standard deviation of the neutron and precursor populations using SSPK have been compared with the exact MC simulation and Stochastic PCA (Allen, 2007) in Table 2. Also the time behavior of the mean neutron population is compared with a sample in Fig. 1. Note that while the population of neutrons and precursors are zero at startup, the final results are close to the MC simulations. This indicates that while the requirements in Section 3 are sufficient for a rigorous mathematical derivation
ð25Þ
we shall refer to Eq. (25) by SSPK equations. Table 3 Computational times for the first benchmark problem.
5. Numerical results Here, we shall measure the accuracy of SSPK by comparing the results of Eq. (25) with published benchmarks for the mean and standard deviations in the point kinetics equation.
Time (s)
MC
PCA
SSPK
Deterministic
6.71
414.34
5.15
4.23
Table 2 Mean and standard deviation of neutron and precursor population for the first benchmark problem.
E[n(0.1)] r[n(0.1)] E[c(0.1)] r[c(0.1)]
Monte Carlo
Stochastic PCA
SSPK
199.15 152.63 1254.5 613.94
204.52 174.03 1294.0 620.68
208.14 174.30 1293.22 622.12
Fig. 2. Time distribution to reach a neutron level of 4000.
Table 4 The mean and standard deviation of time needed to reach 4000 neutrons.
Fig. 1. Sample and average neutron population for first benchmark problem.
E[t](s) r[t](s)
Monte Carlo
Stochastic PCA
SSPK
33.136 2.0886
33.157 2.5772
33.125 2.6731
695
S.M. Ayyoubzadeh, N. Vosoughi / Annals of Nuclear Energy 63 (2014) 691–695 Table 5 The mean and standard deviation for the third benchmark problem.
E[n(0.1)]
r[n(0.1)]
P E½ 6j¼1 cj ð0:1Þ P r½ 6j¼1 cj ð0:1Þ
Monte Carlo
Stochastic PCA
Euler–Maruyama approximation
SSPK
183.04 168.79 4.478 105
186.31 164.16 4.491 105
208.6 255.95 4.498 105
184.8 186.96 4.489 105
1495.7
1917.2
1233.38
982.64
Table 6 Computational times for the third benchmark problem.
Time (s)
MC
PCA
SSPK
Deterministic
7.34
4610.4
11.50
4.67
of a system of SDEs, they may not always be necessary. Obtaining the full range of validity for SSPK is beyond the scope of this paper. In Table 3, we have compared the computational time needed to obtain a realization for this benchmark problem using the stochastic equations with the time needed to solve the forward finite differenced, deterministic PK equations. The implementation in (Hayes, 2005) has been used for the PCA method and the computations have been performed on a 2.2 GHz Core i7 CPU. In the second example, we examine the accuracy of the model by calculating the distribution of the time necessary for a reactor to reach a predefined neutron population. We shall use the (non-physical) parameters used by (Hayes and Allen, 2005) as b = 0.05, k ¼ 0:1 ðs1 Þ, m = 2.5, q = 200 (s1), K = 0.499002 (s), and q = 0.001996 for a clean reactor. The distribution of the time to reach 4000 neutrons is shown in Fig. 2. The mean and standard deviation of this distribution for 5000 realizations is obtained using SSPK and compared with the values obtained from MC and Stochastic PCA (Hayes and Allen, 2005) methods in Table 4. In the last benchmark problem, we consider the 6-group delayed neutron model of a pool-type light-water research reactor (Chao and Attard, 1985; Kinard and Allen, 2004; Saha Ray, 2012; Yang and Jevremovic´, 2009). The parameters of this reactor are that of the VR-1 research reactor and equal to b = [266, 1491, 1316, 2849, 896, 182] 106, k ¼ ½127; 317; 1150; 3110; 14000; 38700 104 ðs1 Þ, m = 2.5, q = 0, K = 2 105 (s), and q = 0.003 (Matejka et al., 2003). The equilibrium conditions are considered at the time origin, namely ½ nð0Þ c1 ð0Þ . . . c6 ð0Þ ¼ 100 ½ 1 kb11K . . . kb66K . The result has been compared at t = 0.1(s) for 5000 realizations with (Saha Ray, 2012) in Table 5. The computational times required for obtaining a realization using the MC, PCA, SSPK and deterministic methods for this reactor have been summarized in Table 6. As it is seen, eliminating the computation of the square root of a matrix has reduced the computational time required for SSPK to solve this problem nearly 400 times in comparison with PCA (Hayes and Allen, 2005). 6. Conclusion In this paper, an alternative derivation of the stochastic differential equations that describe a point reactor has been presented. A few approximations have been used to simplify the final form
of the equations. The validity of the resultant equations, namely the Simplified Stochastic Point Kinetics equations (SSPK), has been confirmed by numerical benchmarks which show good agreement with published results. Possible future works include improving the current work by advising advanced numerical SDE algorithms and generalization of these equations to include spatial and multigroup neutron effects. References Akcasu, A.Z., Stolle, A., 1989. Comments on the noise equivalent source in the Langevin technique. Annals of Nuclear Energy 16, 493–508. Allen, E., 2007. Modeling with Itô Stochastic Differential Equations. Springer. Bell, G.I., 1969. Stochastic formulations of neutron transport. SIAM-AMS Proceedings. American Mathematical Society, 181–197. Chao, Y.-A., Attard, A., 1985. A resolution of the stiffness problem of reactor kinetics. Nuclear Science and Engineering; (United States) 90. Degweker, S., 1994. A forward equation for stochastic neutron transport. Annals of Nuclear Energy 21, 531–539. Eagleson, G., 1975. On Gordin’s central limit theorem for stationary processes. Journal of Applied Probability, 176–179. Evans, L.C., 2002. An Introduction to Stochastic Differential Equations. UC Berkley. Gençay, Sß., Akcasu, Z., 1981. Noise analysis of a zero power critical line reactor using Langevin’s technique. Annals of Nuclear Energy 8, 65–77. Gillespie, D.T., 2000. The chemical Langevin equation. The Journal of Chemical Physics 113, 297. Hayes, J.G., 2005. Stochastic Point Kinetics Equations in Nuclear Reactor Dynamics, Department of Mathematics. Texas Tech University. Hayes, J.G., Allen, E., 2005. Stochastic point-kinetics equations in nuclear reactor dynamics. Annals of Nuclear Energy 32, 572–587. Henry, A.F., 1975. Nuclear-Reactor Analysis. MIT press Cambridge. Hetrick, D.L., 1993. Dynamics of Nuclear Reactors. American Nuclear Society (La Grange Park, Ill., USA). Ikeda, N., Watanabe, S., 1989. Stochastic Differential Equations and Diffusion Processes. Kinard, M., Allen, E., 2004. Efficient numerical solution of the point kinetics equations in nuclear reactor dynamics. Annals of Nuclear Energy 31, 1039– 1051. Lewins, J., 1978. Linear stochastic neutron transport theory. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 362, 537–558. Matejka, K., Kolros, A., Kropik, M., Sklenka, L., 2003. Eugene Wigner Training Course at VR-1 Reactor. Department of Nuclear Reactors, CTU in Prague. Øksendal, B., 2003. Stochastic Differential Equations. Springer. Osborn, R., 1969. Kinetic theory of transport and fluctuation phenomena. SIAM-AMS Proceedings. American Mathematical Society, p. 213. Osborn, R.K., Yip, S., 1966. The foundations of neutron transport theory. In: Osborn, Richard K., (Ed.), The foundations of neutron transport theory. Yip, Sidney. New York, Gordon and Breach [1966]. Monograph Series on Nuclear, Science and Technology 1. Pázsit, I., Pál, L., 2008. Neutron Fluctuations: A Treatise on the Physics on Branching Processes. Elsevier Science Limited. Saha Ray, S., 2012. Numerical simulation of stochastic point kinetic equation in the dynamical system of nuclear reactor. Annals of Nuclear Energy 49, 154–159. Williams, M.M.R., 1974. Random Processes in Nuclear Reactors. Yang, X., Jevremovic´, T., 2009. Revisiting the Rosenbrock numerical solutions of the reactor point kinetics equation with numerous examples. Nuclear Technology and Radiation Protection 24, 3–12. Zolotarev, V.M., 1986. One-dimensional stable distributions. American Mathematical Soc.