One million neutrons were simulated, with a packing factor of 0.093, to compare the linear and segmented scripts with the projectile collision method. The results are listed in table 6.5. Additionally, the non-leaked free path length distribution can be seen from figure 6.3.
Table 6.5. Projectile collision results, using afp of 0.093, compared to the reference model.
Description Symbol Unit Projectile collision Reference model
Neutrons nn [−] 1.0·106 1.0·106
Leaked neutrons nl [−] 39 −
Mean free path
¯l [cm] 0.587
0.595
Leakless mean free path 0.587
Error a [cm] ±0.003 ±0.003
r [%] ±0.470 ±0.475
Leakless error a [cm] ±0.003 ±0.003
r [%] ±0.468 ±0.475
Significance α [−] 4.750σ 4.750σ
Median
˜l [cm] 0.409
0.413
Leakless median [cm] 0.408
Min lmin [cm] 1.151·10−4
2.641·10−7
Leakless min 1.151·10−4
Max lmax [cm] 10.470
8.163
Leakless max 7.491
Figure 6.3. Projectile collision method free path length distribution with fp of 0.093.
The projectile collision detection method was also tested with a packing factor of 17 % using 1,000,000 neutrons. These results were then compared to an equal amount of ran-domly generated mean free paths using equation (4.3). The results are shown in table 6.6. The comparison plot is shown in figure 6.4. The packing factor experienced by the neutron was sampled 400,000 times and these results are included in table 6.6.
Table 6.6. Projectile collision results, using afp of 0.17, compared to the reference model.
Description Symbol Unit Projectile collision Reference model
Neutrons nn [−] 1.4·106 1.4·106
Mean free path ¯l [cm] 0.292 0.299
Error a [cm] ±0.001 ±0.001
r [%] ±0.205 ±0.210
Significance α [−] 4.750σ 4.750σ
Observed packing factor fp,o [−] 0.1729 0.1700
Median ˜l [cm] 0.204 0.207
Min lmin [cm] 3.509·10−7 2.492·10−8
Max lmax [cm] 4.279 4.407
Figure 6.4. Projectile collision method and reference model results with a packing factor of 0.17.
In order to determine possible differences between the distributions of neutrons leaving from the surfaces of fuel particles and those generated in the graphite matrix, both cases were simulated 400,000 times. Based on these results, the combined mean free path was calculated. The results are listed in table 6.7 and they are plotted in figure 6.5.
Table 6.7. A comparison of results between generating the neutrons on the surface and off the surface of fuel particles, with a fp of 0.017, using
the projectile collision method.
Description Symbol Unit On surface In matrix Combined
Neutrons nn [−] 4·105 4·105 8·105
Mean free path ¯l [cm] 0.293 0.284 0.288
Error a [cm] ±0.002 ±0.002 ±0.002
r [%] ±0.733 ±0.755 ±0.526
Significance α [−] 4.750σ 4.750σ 4.750σ
Median ˜l [cm] 0.205 0.196 0.201
Min lmin [cm] 1.753·10−4 5.426·10−8 5.426·10−8
Max lmax [cm] 3.883 3.679 3.883
Figure6.5.Projectilecollisionmethodfreepathlengthdistributioncomparisonbetweenneutronsgeneratedoff(top)oron(middle) thesurfaceofaparticle,andthecombineddistribution(bottom),withfpof0.17.Therelativeproportionsoftheoffandon distributionsisnotknownforarealcaseandtheycanthusbedifferent.
7 EVALUATION
The reference results received with Serpent, seen from the table 6.1, using the explicit model show an initial increase in the multiplication factor as a function of the packing factor. It peaks at a packing factor of 0.04, after which starts to constantly decrease. This peak is caused by the presence of an ideal mixture of fuel and moderator. Below 0.04 the amount of moderator present is larger than necessary, reducing the amount of inter-actions between neutrons and fuel. Above it the increased presence of U-238, combined with the reduction of the amount of moderator, causes more neutrons to be captured and less thermal neutrons to be produced and viable for fission.
A similar trend can be seen with the use of the implicit model. The initial increase of the neutron multiplication factor peaks at 4 %, after which it declines. The results differ noticeably from the explicit models from a packing factor of 7 % upwards. The error relative to the explicit results increase on average by 0.06 % for every 1 % of pack-ing factor, rispack-ing from 0.11 % at a 7 % packpack-ing factor to 0.87 % at 20 %. These errors correspond to reactivity differences of 97 pcm and 977 pcm, respectively. The decreasing trend in the multiplication factors after the peak suggests that an increase in the packing factor causes a decrease in the neutron multiplication factor and, conversely, a decrease in the packing factor leads to an increase in the neutron multiplication factor. This implies that the implicit model observes a lower packing factor than is actually used. This would lead to a conclusion that the current model’s interaction mean free path, or the mean particle distance, given by the equation (4.2), is too large for packing factors above 7 %.
The results from the linear and segmented programs, seen from table 6.3 are greater than the reference model’s interaction mean free path but lie within the confidence inter-val. The neutrons were generated outside of the fuel particles and were found to follow an exponential distribution similar to the reference model, as can be seen from figures 6.1 and 6.2. This would suggest that the current model is still accurate enough around a packing factor of 0.093. The comparison between the reference results and the linear program, seen from table 6.4, suggests that the model produces too short interaction mean free paths up to a packing factor of 0.12, after which they are too large. This would, given the previously suggested explanation, mean that the neutron multiplication factors below a packing factor of 0.04 should be lower for the implicit method, due to the apparent lower packing density. Above 0.04 the multiplication factor would thus also be lower up to 0.12 due to an apparent increase in packing factor. Past 0.12, the neutron multiplication factor would be greater for the implicit model due to the lessened effect of the increase of
packing factor on the neutron multiplication factor. The latter is found to be true in the aforementioned results from Serpent.
In the case of the projectile collision detection method, the interaction mean free path received with a packing factor of 0.093, seen from table 6.5, was approximately 1.34 % lower than that of the reference model and thus the reference value was outside of the confidence interval of the simulated value. 39 neutrons leaked from the simulated space, thus causing an increased bias towards the higher values of the mean free paths. This did not, however have a noticeable effect on the mean. This result suggests that the neutron multiplication factor at 0.093, given by the reference model, is higher than the real value.
The free path length distribution does not follow an exponential distribution, however, as can be seen from figure 6.3, due to the sudden drop in the probability near the origin.
Given the Serpent reference results, it can be deduced that the reference model’s interac-tion mean free path at a packing factor of 0.093 is too high. The simplicity of this method implies that the results received from the linear and segmented programs are not accu-rate. This inaccuracy can be due to mistakes in programming or intrinsic flaws in its logic.
The projectile collision method was also used with a higher packing factor of 0.17, seen from table 6.6, in order to further determine the validity of the suggested explanation.
In the first simulation one million neutrons were simulated to leave from the surfaces of fuel particles and 0.4 million were generated on the boundaries to determine the observed packing factors in addition to the free path length. The resulting interaction mean free path was approximately 2.46 % lower than the reference result. This result follows the explanation as a lower mean free path means that the packing factor observed by the neutrons is higher than what is used in the model. This is confirmed by the simulation’s resulting observed packing factor, which is found to be 1.71 % higher than the actual packing factor of the fuel particles. The free path length distribution, as seen from figure 6.4, is again similar to the exponential distribution with the exception of the free path lengths close to 0, the probabilities of which increase steeply to the peak probability, after which it follows an exponential curve.
The final simulation was performed to determine the possible differences between the interaction free path length distributions between neutrons that are generated on the sur-faces of the fuel particles and those generated in the graphite matrix. It can be seen from tables 6.7 and 6.6 that the neutrons generated on the surface appear to have a mean free path between 0.292 and 0.293 cm. The neutrons generated inside the matrix, however, have a shorter mean free path of only 0.284 cm. The average mean free path, when the results are combined, is approximately 0.288 cm, that is 0.011 cm or 3.68 % shorter than
that of the reference model. The short mean free path of the neutrons inside the graphite matrix can be explained by a greater probability of the neutron being located near the surface of a fuel particle. The difference can also be seen in the free path length distribu-tions shown in figure 6.5. The neutrons generated on the surfaces still follow the previous trend of a decreasing probability near the source. The path lengths of the neutrons inside the matrix, however, appear to follow the exponential distribution similarly to the linear and segmented models’ results. This disparity between the distributions would suggest that a single exponential distribution is not enough in the determination of the interaction mean free paths of neutrons.
Following the conclusion that the exponential distribution is not fit for the data, an-other type of a distribution would be necessary. Multiple alternative distribution models, including Erlang and Gamma distributions, were tested with a nonlinear fitting for their parameters. The Gamma distribution was the closest fit but was not able to follow the distribution of the data in the vicinity of its peak.
8 CONCLUSIONS
The error apparent in the results of pebble bed reactor simulations with a packing factor higher than 7 % in Serpent was expected to stem from the current model used to estimate the mean distances between particles. This model gives, as the interaction mean free path with a packing factor of about 9.34 %, a length of 0.595 cm. This was tested with three different programs with results ranging from 0.587 to 0.598. The increase of the neutron multiplication factor as a function of decreasing packing factor lead to the deduction that the interaction mean free path given by the model after 7 % is too large.
It was found that the free path lengths of neutrons generated inside the graphite mod-erator matrix, between the fuel particles, follow an exponential distribution similar to the current model. The neutrons generated on the outer surfaces of the fuel particles, however, follow a distribution nearly identical to the exponential distribution, with the exception of the probability of the free path length dropping close to the origin. That is, the probability of a given free path length for a neutron leaving the surface of a fuel particle is initially zero and increases steeply to maximum shortly after, after which it follows an exponential curve. These results are found to be similar to those received by Torquato et al. [14]
The implicit model was found to produce a neutron multiplication factor that increases with the packing factor, that is the opposite of the referential explicit model’s results, where the packing factor’s increase decreases the multiplication factor. As an explanation for this it was suggested that the packing factor observed by neutrons differs from the actual packing factor of the fuel-filled space. The results from the simulations agreed with this suggestion. The simulated mean free paths were lower than that given by the model, suggesting a higher observed packing factor. This was confirmed with a ray tracing sim-ulation performed inside the same fuel space. Thus the data strongly suggests that the exponential distribution used to model the interaction mean free path of neutrons is not valid.
REFERENCES
[1] J.J. Duderstadt and L.J. Hamilton. 1976. Nuclear reactor analysis. John Wiley &
Sons.
[2] IAEA. 1997. Fuel performance and fission product behaviour in gas cooled reac-tors. [E-document]. From: http://www-pub.iaea.org/books/iaeabooks/5633/Fuel-Performance-and-Fission-Product-Behaviour-in-Gas-Cooled- Reactors. [Retrieved 5th June, 2012].
[3] IAEA. 2001. Current status and future development of modular high tem-perature gas cooled reactor technology. [E-document]. From: http://www- pub.iaea.org/books/IAEABooks/6124/Current-Status-and-Future-Development-of-Modular-High- Temperature-Gas-Cooled-Reactor-Technology. [Retrieved 4th June, 2012].
[4] A.C. Kadak. 2005. A future for nuclear energy: pebble bed reactors. Int. J. Critical Infrastructures, Vol. 1, No. 4, pp. 330–345. [E-document]. From: http://web.mit.edu/pebble-bed/papers1 files/Future%20for%20Nuclear%20Energy.pdf. [Retrieved 4th June, 2012].
[5] B. Lapeyre, ´E. Pardoux, R. Sentis. 2003. Introduction to Monte-Carlo Methods for Transport and Diffusion Equations. Oxford University Press: New York.
[6] J. Lepp¨anen. 2007. Development of a New Monte Carlo Reactor Physics Code. [E-document]. From: http://montecarlo.vtt.fi/download/P640.pdf. [Retrieved 4th June, 2012].
[7] J. Lepp¨anen. 2007.Randomly Dispersed Particle Fuel Model in the PSG Monte Carlo Neutron Transport Code. Joint International Topical Meeting on Mathematics &
Computation and Supercomputing in Nuclear Applications, Monterey, California, 15th–19th April, 2007.
[8] C. Liang, W. Ji and F.B. Brown. 2013. Chord length sampling method for analyzing stochastic distribution of fuel particles in continuous energy simulations. Annals of Nuclear Energy, Vol. 53, pp. 140–146.
[9] B. Lu and S. Torquato. 1993. Chord-length and free-path distribution functions for many-body systems. Journal of Chemical Physics, Vol. 98, No. 8, pp. 6472–6482.
[10] I. Murata, T. Mori and M. Nakagawa. 1995. Continuous Energy Monte Carlo Cal-culations of Randomly Distributed Spherical Fuels in High-Temperature Gas- Cooled Reactors Based on a Statistical Geometry Model. Nuclear science and engineering, Vol. 123, No. 1, pp. 96–109.
[11] C.H. Rycroft, T. Lind, S. G¨untay and A. Dehbi. 2012 Granular flow in pebble bed reactors: dust generation and scaling. 2012 International Congress on Advances in National Power Plants, Chicago, Illinois, 24th–28th June, 2012.
[12] Serpent - A Monte Carlo Reactor Physics Burnup Calculation Code. [VTT web page].
From: http://montecarlo.vtt.fi/. [Retrieved 25th January, 2013].
[13] N.G. Sj¨ostrand. 2002. What is the average chord length?. Annals of Nuclear Energy, Vol. 29, No. 13, pp. 1607–1608.
[14] S. Torquato, B. Lu and J. Rubinstein. 1990.Nearest-neighbor distribution functions in many-body systems. Physical Review A, Vol. 41, No. 4, pp. 2059–2075.
[15] World Nuclear Association. 2012.Nuclear Power in the World Today. [E-document].
From: http://www.world-nuclear.org/info/inf01.html. [Retrieved 7th November, 2012].
A LINEAR PROGRAM
43 w h i l e (any(ID( : , 1 ) ) | | any(ID( : , 2 ) ) ) && ITER <= 1e2
87 %−−− G e n e r a t e c o o r d i n a t e s i t e r a t i v e l y
133 OUT = nn−sum(id) ;
179 B( : ,i) = b;
224 case {’ f ’,’ b ’}
270
271 % Update t h e n e u t r o n c o o r d i n a t e s .
272 i f strcmp(TRIG, ’ f ’) | | strcmp(TRIG, ’ b ’)
273 NC(r, : ) = NC(r, : ) +t(r)∗ones( 1 , 3 ) .∗ND(r, : ) ; 274 e l s e i f strcmp(TRIG, ’ b f ’)
275 NC(r, : ) = NC(r, : )−t(r)∗ones( 1 , 3 ) .∗ND(r, : ) ; 276 end
277
278 % Move t h e n e u t r o n s t o a random l o c a t i o n on t h e e d g e o f t h e f u e l p e b b l e . 279 f o r i = 1 :l e n g t h(r)
280 [x,y,z] = s p h 2 c a r t( 2∗p i∗rand( 1 ) , ( 2∗rand( 1 )−1)∗p i,BR) ; 281 NC(r(i) , : ) = [x y z] ;
282
283 [ox,oy,oz] = s p h 2 c a r t( 2∗p i∗rand( 1 ) , ( 2∗rand( 1 )−1)∗p i, 1 ) ; 284 ND(r(i) , : ) = [ox oy oz] ;
285 end
B CELL TRACKING PROGRAM
41 %−−− I n d e x i n g o f t h e p a r t i c l e s and n e u t r o n s . Space d i v i d e d i n t o
ismember(iPxy,ind, ’ rows ’) | ismember(iPxz,ind, ’ rows ’) |
ismember(iPxym,indexB, ’ rows ’) | ismember(iPxzm,
142 i f id( 1 ) , [tfp,ncF,ndF] = sub_trs(’ f ’,tfp,id,BR,ncF
176 i f norm(PC( 1 , : ) ) <= BR−PR, D = 0 ; end
222 id(i) = true;
268 switch TRIG
314 cB = sum(NCB. ˆ 2 , 2 )−BRˆ 2 ;
B.2 Cell division
44 % Track t h e n e u t r o n .
86 xyzB = f l i p u d(xyzB) ;
C PROJECTILE COLLISION
41 % PARTICLES IN WAY
22 0 1 0
68 count = 0 ;
111 end
112 end
113
114 % Add t h e p a r t i c l e c o o r d i n a t e s t o t h e combined c o o r d i n a t e 115 % m a t r i x o f t h e s p a c e .
116 c_spa( (bid−1)∗n_par+1:bid∗n_par, : ) = c_box;
117 end
118 end
119 end
D EXPERIENCED PACKING FACTOR
43
89
134 i f s i g n(tb) == −1
135 t = s o r t(−t(t < −tol) ) ;
136 tb =−tb;
137 e l s e i f s i g n(tb) == 1 138 t = s o r t(t(tol < t) ) ;
139 end
140
141 % Form t h e v e c t o r f o r t h e t i m e s o f c o l l i s i o n , i n c l u d i n g t h e i n i t i a l and 142 % f i n a l b o u n d a r i e s h i t and d e t e r m i n e t h e t i m e d i f f e r e n c e s .
143 T = [ 0 ;t;tb] ; 144 dT = d i f f(T) ; 145
146 % Determine t h e d i s t a n c e s t h r o u g h p a r t i c l e s ( dp ) and v o i d ( dv ) 147 % t r a v e l l e d .
148 i f piny
149 dp = dT( 1 : 2 :end) ; 150 dv = dT( 2 : 2 :end) ; 151 e l s e
152 dp = dT( 2 : 2 :end) ; 153 dv = dT( 1 : 2 :end) ;
154 end
155
156 % C a l c u l a t e t h e p a c k i n g f a c t o r e x p e r i e n c e d ( d f ) and t h e f i r s t f r e e path 157 % between p a r t i c l e s ( f p ) .
158 df(i) = sum(dp) / (sum(dp)+sum(dv) ) ; 159 fp(i) = dv( 1 ) ;
160 f p r i n t f(’ %1.0 f . mean p f | f p : %1.5 f | %1.5 f cm\n ’,i,mean(df( 1 :i) ) ,mean( fp( 1 :i) ) ) ;
161 end 162
163 % C a l c u l a t e t h e means . 164 DF = mean(df) ;
165 BFP = mean(fp) ; 166
167 display(’ ’) ;
168 f p r i n t f(’ Mean v i r t u a l p f : %1.7 f\n ’,DF) ; 169 f p r i n t f(’ Mean f r e e path : %1.7 f cm\n ’,BFP) ;
E SHARED
E.1 A faster dsearchn
1 f u n c t i o n d = dsrchn(x,xi) 2
3 % A s i m p l i f i e d , f a s t e r v e r s i o n o f d s e a r c h n f o r t h e p u r p o s e s needed . 4
5 % I n i t i a l i z e t h e d i s t a n c e v e c t o r . 6 d = z e r o s(s i z e(xi, 1 ) , 1 ) ;
7
8 % C a l c u l a t e t h e s q u a r e s o f minimum d i s t a n c e s t o x f o r e a c h s e t o f 9 % c o o r d i n a t e s o f x i .
10 f o r i = 1 :s i z e(xi, 1 )
11 yi = ones(s i z e(x, 1 ) , 1 )∗xi(i, : ) ; 12 d(i) = min(sum( (x−yi) . ˆ 2 , 2 ) ) ; 13 end
14
15 % Take t h e s q u a r e r o o t s o f t h e s q u a r e s . 16 d = s q r t(d) ;