- Email: [email protected]
Determination of the effective delayed neutron fraction for training reactor VR-1
Progress in Nuclear Energy 52 (2010) 735e742
Contents lists available at ScienceDirect
Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene
Determination of the effective delayed neutron fraction for training reactor VR-1 b a, Ján Rataj c, Svetozár Michálek a, *, Stanislav Stevo , Gabriel Farkas a, Ján Has cík a, Vladimír Slugen c Antonín Kolros Slovak University of Technology, Faculty of Electrical Engineering and Information Technology, Department of Nuclear Physics and Technology, Ilkovicova 3, 81219 Bratislava, Slovakia b Slovak University of Technology, Faculty of Electrical Engineering and Information Technology, Institute of Control and Industrial Informatics, Ilkovicova 3, 81219 Bratislava, Slovakia c Czech Technical University in Prague, Faculty of Nuclear Sciences and Physical Engineering, Department of Nuclear Reactors, V Holesovickách 2, 18000, Prague 8, Czech Republic a
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 January 2010 Received in revised form 6 May 2010 Accepted 14 May 2010
This paper describes a methodology of effective delayed neutron fraction (beff) determination for the VR1 training reactor located at the Czech Technical University (CTU) in Prague. The methodology contains both experimental and computational elements. The experimental determination is based on a measurement of critical, zero-power reactor response to a periodical, low-worth reactivity insertion induced by fast periodical movement of a cadmium absorber in a special pneumatic system. From the power oscillations measurement, the reactor transfer function is determined using a genetic algorithm in the numerical program MATLAB. A final value of delayed neutron fraction is then determined directly from the amplitude of the transfer function. The measurements are verified by calculations using the stochastic MCNP5 code. The results of this calculation and measurement differ only by approximately 1.76% and indicate that the currently used value of beff on VR-1 has been underestimated by more than 8%. Ó 2010 Elsevier Ltd. All rights reserved.
Keywords: Effective delayed neutron fraction (beff) Reactivity Transfer function Genetic algorithm Monte Carlo method
1. Introduction Despite the fact that delayed neutrons constitute only a very small fraction of the total number of neutrons generated from fission, they play a dominant role in fission chain reaction control. If only the prompt neutrons existed, reactor operation would become impossible, due to rapid reactor power changes. Therefore, the exact determination of the value of the effective delayed neutron fraction (beff) is one of the main requirements in reactor physics. The value of beff for the VR-1 training reactor has never been determined experimentally. beff has only been estimated empirically considering the physical properties of the core and the used jka et al., 2003). Since 2005, a new type fuel, as beff ¼ 0.00714 (Mate of fuel with lower enrichment of 19.7% is being used in the VR-1 reactor. This has influenced the currently used beff value and has served as a basis for research into a methodology for its experimental determination. This process is based on the theory of the linearized reactor of zero-power. It begins with development of a methodology of beff measurement, continues with performance of
* Corresponding author. Tel.: þ421 908 500 601; fax: þ421 2 6542 7207. E-mail address: [email protected] (S. Michálek). 0149-1970/$ e see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.pnucene.2010.05.002
several reactor measurement series, and ends with verification of the achieved results by calculations, performed by stochastic MCNP5 code. The method of our measurement of beff on the VR-1 training reactor is the in-pile kinetic technique, which uses the reactor itself as a measurement system. This method is based on periodical, lowworth reactivity insertion to a critical zero-power reactor by periodical movement of control rod and following power oscillations measurement. From this measurement, the reactor transfer function is determined by fitting the reactor response function, using Genetic Algorithms (GA) in the numerical program MATLAB. A final value of effective delayed neutron fraction is then determined directly from the amplitude of the reactor transfer function. A similar approach can be found in contribution of Yedvab et al. (2006) and another method of transfer function-based beff determination in Khamis et al. (2004). The verification, by stochastic Monte Carlo method calculations, was performed using the prompt method published by Bretscher (1997). MCNP5 code was selected for this purpose. Other Monte Carlo approaches of beff calculation can be found in contributions from Spriggs et al. (2001), Nauchi and Kameyama (2005), Van der Mark and Klein Meulekamp (2006), Nagaya and Nakajima (2008), Chiba (2009), and Irwanto et al. (2009).
736
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742
Fig. 1. Reactor transfer function amplitude response for
235
U and
2. Theoretical analysis The transfer of linearized model of zero-power reactor can be expressed as follows (Hermanský, 1987):
" #1 m ~ X bi 1 DPðsÞ 1 L þ ; ¼ G0 ðsÞ ¼ rðsÞ P0 ~ s s þ li
238
U systems. Parameter [p e prompt neutron lifetime (Keepin, 1965).
small reactivity insertions, does not behave as a harmonic oscillator, due to the increase of reactor power mean value. Substituting s ¼ ju in Eq. (1) a frequency response is obtained as:
" (1)
G0 ðjuÞ ¼ juL þ
i¼1
where G0(s) is linearized zero-power reactor transfer function, L e prompt neutron generation time, bi, li e delayed neutron fraction, decay constant of ith-group of delayed neutron precursors, s e parameter of Laplace-transform. The graphical interpretation of Eq. (1) is expressed as:
r/½G0 /DP=P:
(2)
A basic characteristic of transfer G0 is its dependence on the initial power level P0. If the amplitude of the harmonic reactivity changes is constant, the amplitude of the power deviation is proportional to the initial (steady state) power level P0. Hence, the gain is based on reactor power. The zero-power reactor, even at
m X jubi u j þ li i¼1
#1 (3)
Fig. 1 shows the amplitude frequency response of linearized transfer G0 for delayed neutron parameters form 235U and 238U fission. It successfully corresponds with the measured values. Therefore usage of a transfer in form of G0 is widely applied, e.g. in reactor control systems design (Fig. 2). This logarithmic frequency response of transfer G0(s) has three areas: 1. Slow reactivity changes
limin
s
0 G0 ðsÞ ¼
1 1
sbl
0 jG0 ðjuÞj ¼
Fig. 2. Measurement (in critical state) of neutron source (S) contribution to the increase of the VR-1 reactor power level.
1
ubl1
(4)
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742 Table 1 Results of the non-dimensional reactivity measurement (parameters: N0 ¼ 5E5 s1; S ¼ 101.4 s1; [ ¼ 0.0942 s*; beff ¼ 0.00714**).
h [mm]
Nh [s1] keffh
r1 []
400 525 650
9564 12 388 16 755
0.000228 0.000201 0.031933 0.028151
0.99900 0.99923 0.99943
r2 []
r1 [$]
2. Intermediate reactivity changes
s^L
s^L
m X bi 1 1 0 G0 ðsÞ ¼ 0 jG0 ðjuÞj ¼ u l [p s[ s þ p i i¼1
(6)
r2 [$]
* Value for 235U. ** Value for approx. determination of r [$] as r½$ ¼ r½ =beff .
limax
3. Fast reactivity changes
limax
Absorber Reactor Effective Reactivity position power multiple factor
737
m X bi 1 1 0 G0 ðsÞ ¼ 0 jG0 ðjuÞj ¼ l b b s þ i i¼1
(5)
3. Measurements on VR-1 reactor The development of the methodology of beff measurement on VR-1 training reactor was done over a period, a year and a half and consisted of 6 measurement series. Results of the last measurement series are presented hereinafter. As previously mentioned, the trapezoidal reactivity perturbations were caused by periodical movement of a cadmium strip embedded in a polyethylene capsule oscillating in a special pneumatic system. Quantification of the reactivity was done based on additional
Fig. 3. Measurement in subcritical state performed on VR-1 reactor.
Fig. 4. Reactor response to the periodical changes of reactivity e v ¼ 999 mm/s, oscillation period 2.1667 s.
738
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742
Fig. 5. Approximation of the VR-1 reactor power response by ARMAX e velocity 999 mm/s, oscillation period 2.1944 s, TF order ¼ 4, LCM.
measurement. For evaluation of the measurement it was required to set non-dimensionally (r[]). Therefore, a measurement based on the theory of source multiplication was performed. According to this theory, the reactivity can be determined from following relation:
r¼
keff2 keff1 1 1 zDkeff ¼ S[ N N keff1 N1 N2
(7)
This determination has disadvantage in the necessity of knowledge of the total average neutron lifetime [. To determine the
Fig. 6. Approximation of the VR-1 reactor power response by GA e velocity 900 mm/s, NG ¼ 6000, oscillation period 2.1711 s, TF order ¼ 3, LCM.
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742 Table 2 Results of GA fitting for the case of LCM310 detection system, NG ¼ 6000. v [mm/s]
T [s]
u [rad/s]
jGðsÞj
jGðsÞj1
100 200 300 400 500 600 700 800 900 999
6.3925 3.8267 3.0378 2.6589 2.3978 2.3533 2.2033 2.1911 2.1711 2.1667
0.982900 1.641933 2.068334 2.363077 2.620396 2.669947 2.851716 2.867594 2.894010 2.899887
140.6340* 132.4967* 130.6938 128.0086 127.0678 126.2476 125.0249 124.6390 127.1129 127.2720
0.007111 0.007547 0.007651 0.007812 0.007870 0.007921 0.007998 0.008023 0.007867 0.007857
* Calculated using NG ¼ 12 000 to reach the optimal fitness.
influence of the selection of its value on the result, a sensitivity analysis was conducted in numerical program Bokin20001. In this program, it is possible to simulate the reactor response to various periodical- or step-changes of inserted reactivity. These changes can be of different shapes, e.g. sinusoidal, trapezoidal, saw-shaped, or various step changes. This program requires entering a few parameters such as: delayed neutron group parameters, initial conditions, shape of the inserted reactivity, source position, etc., to simulate the real reactor behavior during reactivity transient. The aim of the program is to determine the solution of one-point kinetic equations for six groups of delayed neutrons, external source count and reactivity variable in time. Simulations of the reactor response to sinusoidal reactivity insertion were performed. The worth of the reactivity was set nondimensionally proceeding from the saturated levels of counts (power levels) and keff parameters according to Table 1. Its value was changed according to the total average neutron lifetime [ as a parameter in the interval of <0.09 O 0.1> s sampled by 0.001 s. The standard MATLAB identificators ARX, ARMAX and PEM were chosen for the reactor response function approximation. These were compared with Genetic Algorithm (GA) with input parameters: NG ¼ 150 and TF order ¼ 2 providing almost identical results. As a conclusion to this analysis it can be stated that the value of parameter [ in the investigated range influences the value of beff only slightly. Relative uncertainty, according to the reference value set as a required program input beff ¼ 0.006885, lies in the interval <1.05 O 1.12> %. According to this conclusion, in the measurement of the non-dimensional value of reactivity on VR-1 reactor, the value of total delay time [d y[ ¼ 0:0942 s for 235U is used for beff experimental determination. Results of the reactivity measurement are summarized in Table 1. Results of the non-dimensional reactivity measurement showed that the positive reactivity insertion is slightly less than negative (w0.38 ¢) and that the condition of symmetry was not totally fulfilled. As shown later on, the consequences were not as significant as expected. The reactor response was measured by Low Current Monitor (LCM310), connected to one of the four fission chambers RJ1300, providing the Operational Power Measurement (Fig. 3). Ten oscillation cycles of 10 periods were performed increasing the velocity of absorber movement in the pneumatic system by 100 mm/s, from 100 to 999 mm/s. The measured response for the case of the highest velocity 999 mm/s is shown in Fig. 4. Since the pneumatic system did not provide an independent, real-time measurement of the absorber position, its transition from one
1 Program developed at Dept. of Nuclear Reactors of CTU in Prague by J. Krepel, MSc.
739
position to another had to be reconstructed from the measured response. Following the measurement of the reactor response, the transfer function was determined. With the use of standard MATLAB identificators, the quality of the reactor response approximation is inadequate, which consequently influences the precision of the final result of beff. ARX, ARMAX and PEM identificators were compared in terms of fitting the measured reactor response and following transfer function determination. Approximation of the data measured by LCM and their approximation by the ARMAX identificator for the oscillation cycle with a velocity of 999 mm/s (T ¼ 2.1944 s) is depicted in Fig. 5. The resulting transfer function obtained was:
GðsÞ ¼
s3
1318s2 þ 24750s þ 10800 þ 10:85s2 þ 201:8s þ 25:37
(8)
with the inverse amplitude of 0.007870. The quality of the fit, not the result was an important indicator of credibility of the result. In the case of the remaining identificators, the quality of the fit was far lower than the ARMAX one; therefore other fitting technique had to be found to obtain more precise approximation and more credible result. GA was selected for this purpose. GA fitting for the case of the 900 mm/s fast reactivity oscillations is depicted in Fig. 6. Fig. 6 illustrates that GA is a more suitable approximation technique for final determination of beff. Results of the GA approximation of the measured reactor response considering 6000 generations are summarized in Table 2. The fitness function, representing the absolute difference between the real and by GAapproximated response is depicted in Fig. 7. Constant value of the fitness function during the defined number of generations (NG) in this case represents reaching the global minimum of the function. From Fig. 7 it is obvious that the global minimum has been reached after less than 1000 generations and it remained constant until reaching the ending condition of 6000 generations, which serves as the indicator of reaching the optimal (best possible) quality of the fit (Table 3). Figs. 8 and 9 depict the part of VR-1 transfer function amplitude response close to the boundary period of 1.5 s (u ¼ 4.18 rad/s). Further decrease of the period (increase in angular frequency) should represent approaching the area of intermediate reactivity changes where the amplitude response remains constant until it reaches the other period boundary of approximately 0.06 s (u ¼ 104.72 rad/s). Fig. 9 depicts the same dependence in the case of beff and shows its convergence to the real value. The true value of beff for VR-1 reactor has not yet been discussed. The exact value has not been found by these aforementioned measurements. Further research has to be done in this field. The pneumatic system must be made capable of performing the required harmonic oscillations with a period from the area of the intermediate reactivity changes, i.e. in the interval <0.06 O 1.57> s (Fig. 1). Then, it can be reliably stated that the value of beff for the training reactor VR-1 is higher than the currently used value of 0.00714 and it asymptotically reaches the value, which lies in the interval of <0.0078 O 0.0081>. 4. Verification by Monte Carlo calculation The calculation was performed using MCNP5 code. The model of VR-1 reactor, containing both reactor tanks with all inner reactor components, excluding those summarized below and including the heavy-concrete biological shielding together with surrounding air. Fig. 10 depicts the VR-1 actual reactor core composition and Figs. 11 and 12 the MCNP plots of the reactor.
740
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742
Fig. 7. Fitness function for GA fitting e velocity 900 mm/s, LCM.
In the model, several simplifications were considered: Table 3 Influence of XS library selection in MCNP criticality calculations for beff determination. XS library JENDL-3.3 kp kt
beff
ENDF/B-VI
ENDF/B-VII
0.99244 0.00005 0.99792 0.00005 0.99264 0.00005 1.00016 0.00005 1.00636 0.00005 1.00036 0.00005 0.007719 0.000100 0.008387 0.000100 0.007717 0.000100
a) Pneumatic system (“Hopik”) is not included e model of fuel assembly is used instead, b) Fission chambers RJ1300 (M1-4) for power measurement during operation are not considered in the model, c) SNM detectors are not considered in the model. The calculation was done using Bretscher’s prompt method (Bretscher, 1997). This method consists of two steps. In the first
Fig. 8. Part of the VR-1 training reactor transfer function amplitude response measured by LCM310 detection system.
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742
741
Fig. 9. Dependence of beff on the oscillation period measured by LCM310 detection system. (beff values depicted here were calculated as ßeff ¼ jGðjuÞj1 , although the oscillation periods do not come under the area of intermediate changes.)
step, KCODE calculation of the effective multiplication factor for prompt neutrons (kp) was determined as kp ¼ 0.99244 0.00005. In the second, the value of total effective multiplication factor (kt) was calculated as: kt ¼ 1.00016 0.00005. The final effective delayed neutron fraction for VR-1 reactor was determined from Eq. (9):
beff y1
kp y0:007719 0:000100 kt
(9)
For this calculation, Japanese JENDL-3.3 XS library was chosen as a reference. The selection of cross-section library influences the results of the MCNP criticality calculations. To find the most suitable XS library for the calculation, the influence of XS data has been investigated and is described hereinafter. Also, calculations with the U.S. ENDF/B-VI and B-VII were made. The difference between
beff values based on the use of ENDF/B.VI or JENDL-3.3 libraries, respectively, calculated as 7.3%:
sbeff ¼
jbeff ENDF beff JENDL j
beff JENDL
$100:
(10)
This analysis confirms the results of the benchmark presented by Dos Santos et al. (2006), according to which the ENDF/B-VI library overpredicts beff by more than 7% in beff, beff/L and L calculations in comparison with JENDL-3.3. According to Mosteller (2008), “ENDF/B-VII.0 produces significantly better overall results for the cases in the MCNP criticality validation suite than do ENDF/B-VI, JENDL-3.3, or JEFF-3.1. However, further improvements still are needed in some areas.”, it has been concluded, that the level of ENDF/B-VII accuracy exceeded the JENDL-3.3 library. In our calculations for the case of
Fig. 10. VR-1 reactor core composition C1.
742
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742
MCNP calculated beff ¼ 0.007719 0.000100 overvalues the currently used value by approximately 8.1%. In comparison with the MCNP value, the measurement provides a value overestimated by 1.76%. The 3% difference would be reached by beff ¼ 0.007957. The true value of beff must be determined by further measurements with oscillation periods of approximately 1 s. The performed measurements presented in this work point to the fact that there is a very high probability, that the true value of beff lies in the interval of <0.0078 O 0.0081>. Together with computational verification, they indicate that the currently used value has been underestimated and should be changed. Since the currently used value is based on a higher level of conservatism; this does not require immediate action in terms of reactor operation and safety. Any action should be postponed until measurements with lower reactivity oscillation periods are performed. More accurate results would be achieved by substituting the pneumatically realized reactivity insertions for a linear stepper motor. This would allow harmonic sinusoidal oscillations instead of trapezoidal with inconstant time delay of the absorber in the boundary positions. This was the limitation of the current pneumatic system-driven reactivity oscillations.
Fig. 11. MCNP5 code axial plot of whole reactor.
Acknowledgements The authors would like to thank the VR-1 reactor operational staff from Department of Nuclear Reactors of CTU in Prague for their help and valuable advice during the experiments.
References
Fig. 12. MCNP5 code radial plot of whole reactor.
VR-1 reactor, the difference between the final beff values based on JENDL-3.3 and ENDF/B-VII libraries, is only w0.03%, which leads to the conclusion, that both XS libraries can be used for these calculations, providing almost the same results. 5. Conclusion In this paper, the methodology of the effective delayed neutron fraction (beff) determination in reactor VR-1 has been presented. This methodology serves as the first experimental determination of this important reactor-safety parameter. The contemporary estimated value of beff ¼ 0.00714 is used in reactor operation. If the value for the fastest oscillation cycle performed, beff ¼ 0.007857, is considered, the condition of Eq. (11) published by Dos Santos et al. (2006), describing the difference between calculated (C) and experimental (E) value, is reliably fulfilled.
jC Ej $100 < 3% E
(11)
For the calculation, stochastic Monte Carlo method, realized by means of the MCNP5 code was applied for its verification. The
Bretscher, M.M., 1997. Evaluation of reactor kinetic parameters without the need for perturbation codes. In: Proc. Int. Meeting on Reduced Enrichment for Research and Test Reactors, Wyoming, USA, 5e10 October. Chiba, G., 2009. Calculation of effective delayed neutron fraction using a modified K-ratio method. Journal of Nuclear Science and Technology 46 (5), 399. Dos Santos, A., Diniz, R., Fanaro, L., Jerez, R., Silva, G., Yamaguchi, M., 2006. A proposal of a benchmark for beff, beff/ L, and L of thermal reactors fueled with slightly enriched uranium. Annals of Nuclear Energy 33 (9), 848e855. CSR Hermanský, B., 1987. Dynamics nuclear reactors. MS (in Czech). Irwanto, D., Obara, T., Chiba, G., Nagaya, Y., 2009. Study on calculation methods for effective delayed neutron fraction. In: Proc. 2009 Annual Meeting of the Atomic Energy Society of Japan, Tokyo, Japan, 23e25 March. Keepin, G.R., 1965. Physics of Nuclear Kinetics. Addison Wesley, Reading, Massachusetts, USA. Khamis, I., Hainoun, A., Suleiman, W., 2004. Measurement of the Syrian MNSR delayed neutron fraction and neutron generation time by noise analysis. Annals of Nuclear Energy 31 (3), 331e341. jka, K., Sklenka, L., Kropík, M., Kolros, A., 2003. Eugene Wiegner Training Mate Course at VR-1 Reactor (May 2003). Textbook for practical Course on Reactor Physics in Framework of European Nuclear Engineering Network, DNR FNSPE CTU Prague. p. 11. Mosteller, R.D., 2008. ENDF/B-VII.0, ENDF/B-VI, JEFF-3.1, and JENDL-3.3 results for the MCNP criticality validation suite and other criticality benchmarks. In: Proc. Int. Conf. Reactor Physics, PHYSOR2008, Interlaken, Switzerland, 14e19 September, LA-UR-07-6284. Nagaya, Y., Nakajima, K., 2008. Approximate estimation of effective delayed neutron fraction with correlated sampling method. In: Proc. Int. Conf. Reactor Physics, PHYSOR2008, Interlaken, Switzerland, 14e19 September. Nauchi, Y., Kameyama, T., 2005. Proposal of direct calculation of kinetic parameters beff and L based on continuous energy Monte Carlo. Journal of Nuclear Science and Technology 42, 503. Spriggs, G.D., Busch, R.D., Campbell, J.M., 2001. Calculation of the delayed neutron effectiveness factor using rations of k-eigenvalues. Annals of Nuclear Energy 28, 477. Van der Mark, S.C., Klein Meulekamp, R., 2006. Calculation of the effective delayed neutron fraction by Monte Carlo. Nuclear Science and Engineering 152, 142. Yedvab, Y., Grober, A., Ettedgui, H., Harari, R., Shvetz, I., Caspi, E., 2006. Determination of delayed neutrons source in the frequency domain based on in-pile oscillation measurements. In: Proceedings of PHYSOR-2006, Canadian Nuclear Society, Vancouver, Canada, 10e14 September.
Contents lists available at ScienceDirect
Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene
Determination of the effective delayed neutron fraction for training reactor VR-1 b a, Ján Rataj c, Svetozár Michálek a, *, Stanislav Stevo , Gabriel Farkas a, Ján Has cík a, Vladimír Slugen c Antonín Kolros Slovak University of Technology, Faculty of Electrical Engineering and Information Technology, Department of Nuclear Physics and Technology, Ilkovicova 3, 81219 Bratislava, Slovakia b Slovak University of Technology, Faculty of Electrical Engineering and Information Technology, Institute of Control and Industrial Informatics, Ilkovicova 3, 81219 Bratislava, Slovakia c Czech Technical University in Prague, Faculty of Nuclear Sciences and Physical Engineering, Department of Nuclear Reactors, V Holesovickách 2, 18000, Prague 8, Czech Republic a
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 January 2010 Received in revised form 6 May 2010 Accepted 14 May 2010
This paper describes a methodology of effective delayed neutron fraction (beff) determination for the VR1 training reactor located at the Czech Technical University (CTU) in Prague. The methodology contains both experimental and computational elements. The experimental determination is based on a measurement of critical, zero-power reactor response to a periodical, low-worth reactivity insertion induced by fast periodical movement of a cadmium absorber in a special pneumatic system. From the power oscillations measurement, the reactor transfer function is determined using a genetic algorithm in the numerical program MATLAB. A final value of delayed neutron fraction is then determined directly from the amplitude of the transfer function. The measurements are verified by calculations using the stochastic MCNP5 code. The results of this calculation and measurement differ only by approximately 1.76% and indicate that the currently used value of beff on VR-1 has been underestimated by more than 8%. Ó 2010 Elsevier Ltd. All rights reserved.
Keywords: Effective delayed neutron fraction (beff) Reactivity Transfer function Genetic algorithm Monte Carlo method
1. Introduction Despite the fact that delayed neutrons constitute only a very small fraction of the total number of neutrons generated from fission, they play a dominant role in fission chain reaction control. If only the prompt neutrons existed, reactor operation would become impossible, due to rapid reactor power changes. Therefore, the exact determination of the value of the effective delayed neutron fraction (beff) is one of the main requirements in reactor physics. The value of beff for the VR-1 training reactor has never been determined experimentally. beff has only been estimated empirically considering the physical properties of the core and the used jka et al., 2003). Since 2005, a new type fuel, as beff ¼ 0.00714 (Mate of fuel with lower enrichment of 19.7% is being used in the VR-1 reactor. This has influenced the currently used beff value and has served as a basis for research into a methodology for its experimental determination. This process is based on the theory of the linearized reactor of zero-power. It begins with development of a methodology of beff measurement, continues with performance of
* Corresponding author. Tel.: þ421 908 500 601; fax: þ421 2 6542 7207. E-mail address: [email protected] (S. Michálek). 0149-1970/$ e see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.pnucene.2010.05.002
several reactor measurement series, and ends with verification of the achieved results by calculations, performed by stochastic MCNP5 code. The method of our measurement of beff on the VR-1 training reactor is the in-pile kinetic technique, which uses the reactor itself as a measurement system. This method is based on periodical, lowworth reactivity insertion to a critical zero-power reactor by periodical movement of control rod and following power oscillations measurement. From this measurement, the reactor transfer function is determined by fitting the reactor response function, using Genetic Algorithms (GA) in the numerical program MATLAB. A final value of effective delayed neutron fraction is then determined directly from the amplitude of the reactor transfer function. A similar approach can be found in contribution of Yedvab et al. (2006) and another method of transfer function-based beff determination in Khamis et al. (2004). The verification, by stochastic Monte Carlo method calculations, was performed using the prompt method published by Bretscher (1997). MCNP5 code was selected for this purpose. Other Monte Carlo approaches of beff calculation can be found in contributions from Spriggs et al. (2001), Nauchi and Kameyama (2005), Van der Mark and Klein Meulekamp (2006), Nagaya and Nakajima (2008), Chiba (2009), and Irwanto et al. (2009).
736
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742
Fig. 1. Reactor transfer function amplitude response for
235
U and
2. Theoretical analysis The transfer of linearized model of zero-power reactor can be expressed as follows (Hermanský, 1987):
" #1 m ~ X bi 1 DPðsÞ 1 L þ ; ¼ G0 ðsÞ ¼ rðsÞ P0 ~ s s þ li
238
U systems. Parameter [p e prompt neutron lifetime (Keepin, 1965).
small reactivity insertions, does not behave as a harmonic oscillator, due to the increase of reactor power mean value. Substituting s ¼ ju in Eq. (1) a frequency response is obtained as:
" (1)
G0 ðjuÞ ¼ juL þ
i¼1
where G0(s) is linearized zero-power reactor transfer function, L e prompt neutron generation time, bi, li e delayed neutron fraction, decay constant of ith-group of delayed neutron precursors, s e parameter of Laplace-transform. The graphical interpretation of Eq. (1) is expressed as:
r/½G0 /DP=P:
(2)
A basic characteristic of transfer G0 is its dependence on the initial power level P0. If the amplitude of the harmonic reactivity changes is constant, the amplitude of the power deviation is proportional to the initial (steady state) power level P0. Hence, the gain is based on reactor power. The zero-power reactor, even at
m X jubi u j þ li i¼1
#1 (3)
Fig. 1 shows the amplitude frequency response of linearized transfer G0 for delayed neutron parameters form 235U and 238U fission. It successfully corresponds with the measured values. Therefore usage of a transfer in form of G0 is widely applied, e.g. in reactor control systems design (Fig. 2). This logarithmic frequency response of transfer G0(s) has three areas: 1. Slow reactivity changes
limin
s
0 G0 ðsÞ ¼
1 1
sbl
0 jG0 ðjuÞj ¼
Fig. 2. Measurement (in critical state) of neutron source (S) contribution to the increase of the VR-1 reactor power level.
1
ubl1
(4)
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742 Table 1 Results of the non-dimensional reactivity measurement (parameters: N0 ¼ 5E5 s1; S ¼ 101.4 s1; [ ¼ 0.0942 s*; beff ¼ 0.00714**).
h [mm]
Nh [s1] keffh
r1 []
400 525 650
9564 12 388 16 755
0.000228 0.000201 0.031933 0.028151
0.99900 0.99923 0.99943
r2 []
r1 [$]
2. Intermediate reactivity changes
s^L
s^L
m X bi 1 1 0 G0 ðsÞ ¼ 0 jG0 ðjuÞj ¼ u l [p s[ s þ p i i¼1
(6)
r2 [$]
* Value for 235U. ** Value for approx. determination of r [$] as r½$ ¼ r½ =beff .
limax
3. Fast reactivity changes
limax
Absorber Reactor Effective Reactivity position power multiple factor
737
m X bi 1 1 0 G0 ðsÞ ¼ 0 jG0 ðjuÞj ¼ l b b s þ i i¼1
(5)
3. Measurements on VR-1 reactor The development of the methodology of beff measurement on VR-1 training reactor was done over a period, a year and a half and consisted of 6 measurement series. Results of the last measurement series are presented hereinafter. As previously mentioned, the trapezoidal reactivity perturbations were caused by periodical movement of a cadmium strip embedded in a polyethylene capsule oscillating in a special pneumatic system. Quantification of the reactivity was done based on additional
Fig. 3. Measurement in subcritical state performed on VR-1 reactor.
Fig. 4. Reactor response to the periodical changes of reactivity e v ¼ 999 mm/s, oscillation period 2.1667 s.
738
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742
Fig. 5. Approximation of the VR-1 reactor power response by ARMAX e velocity 999 mm/s, oscillation period 2.1944 s, TF order ¼ 4, LCM.
measurement. For evaluation of the measurement it was required to set non-dimensionally (r[]). Therefore, a measurement based on the theory of source multiplication was performed. According to this theory, the reactivity can be determined from following relation:
r¼
keff2 keff1 1 1 zDkeff ¼ S[ N N keff1 N1 N2
(7)
This determination has disadvantage in the necessity of knowledge of the total average neutron lifetime [. To determine the
Fig. 6. Approximation of the VR-1 reactor power response by GA e velocity 900 mm/s, NG ¼ 6000, oscillation period 2.1711 s, TF order ¼ 3, LCM.
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742 Table 2 Results of GA fitting for the case of LCM310 detection system, NG ¼ 6000. v [mm/s]
T [s]
u [rad/s]
jGðsÞj
jGðsÞj1
100 200 300 400 500 600 700 800 900 999
6.3925 3.8267 3.0378 2.6589 2.3978 2.3533 2.2033 2.1911 2.1711 2.1667
0.982900 1.641933 2.068334 2.363077 2.620396 2.669947 2.851716 2.867594 2.894010 2.899887
140.6340* 132.4967* 130.6938 128.0086 127.0678 126.2476 125.0249 124.6390 127.1129 127.2720
0.007111 0.007547 0.007651 0.007812 0.007870 0.007921 0.007998 0.008023 0.007867 0.007857
* Calculated using NG ¼ 12 000 to reach the optimal fitness.
influence of the selection of its value on the result, a sensitivity analysis was conducted in numerical program Bokin20001. In this program, it is possible to simulate the reactor response to various periodical- or step-changes of inserted reactivity. These changes can be of different shapes, e.g. sinusoidal, trapezoidal, saw-shaped, or various step changes. This program requires entering a few parameters such as: delayed neutron group parameters, initial conditions, shape of the inserted reactivity, source position, etc., to simulate the real reactor behavior during reactivity transient. The aim of the program is to determine the solution of one-point kinetic equations for six groups of delayed neutrons, external source count and reactivity variable in time. Simulations of the reactor response to sinusoidal reactivity insertion were performed. The worth of the reactivity was set nondimensionally proceeding from the saturated levels of counts (power levels) and keff parameters according to Table 1. Its value was changed according to the total average neutron lifetime [ as a parameter in the interval of <0.09 O 0.1> s sampled by 0.001 s. The standard MATLAB identificators ARX, ARMAX and PEM were chosen for the reactor response function approximation. These were compared with Genetic Algorithm (GA) with input parameters: NG ¼ 150 and TF order ¼ 2 providing almost identical results. As a conclusion to this analysis it can be stated that the value of parameter [ in the investigated range influences the value of beff only slightly. Relative uncertainty, according to the reference value set as a required program input beff ¼ 0.006885, lies in the interval <1.05 O 1.12> %. According to this conclusion, in the measurement of the non-dimensional value of reactivity on VR-1 reactor, the value of total delay time [d y[ ¼ 0:0942 s for 235U is used for beff experimental determination. Results of the reactivity measurement are summarized in Table 1. Results of the non-dimensional reactivity measurement showed that the positive reactivity insertion is slightly less than negative (w0.38 ¢) and that the condition of symmetry was not totally fulfilled. As shown later on, the consequences were not as significant as expected. The reactor response was measured by Low Current Monitor (LCM310), connected to one of the four fission chambers RJ1300, providing the Operational Power Measurement (Fig. 3). Ten oscillation cycles of 10 periods were performed increasing the velocity of absorber movement in the pneumatic system by 100 mm/s, from 100 to 999 mm/s. The measured response for the case of the highest velocity 999 mm/s is shown in Fig. 4. Since the pneumatic system did not provide an independent, real-time measurement of the absorber position, its transition from one
1 Program developed at Dept. of Nuclear Reactors of CTU in Prague by J. Krepel, MSc.
739
position to another had to be reconstructed from the measured response. Following the measurement of the reactor response, the transfer function was determined. With the use of standard MATLAB identificators, the quality of the reactor response approximation is inadequate, which consequently influences the precision of the final result of beff. ARX, ARMAX and PEM identificators were compared in terms of fitting the measured reactor response and following transfer function determination. Approximation of the data measured by LCM and their approximation by the ARMAX identificator for the oscillation cycle with a velocity of 999 mm/s (T ¼ 2.1944 s) is depicted in Fig. 5. The resulting transfer function obtained was:
GðsÞ ¼
s3
1318s2 þ 24750s þ 10800 þ 10:85s2 þ 201:8s þ 25:37
(8)
with the inverse amplitude of 0.007870. The quality of the fit, not the result was an important indicator of credibility of the result. In the case of the remaining identificators, the quality of the fit was far lower than the ARMAX one; therefore other fitting technique had to be found to obtain more precise approximation and more credible result. GA was selected for this purpose. GA fitting for the case of the 900 mm/s fast reactivity oscillations is depicted in Fig. 6. Fig. 6 illustrates that GA is a more suitable approximation technique for final determination of beff. Results of the GA approximation of the measured reactor response considering 6000 generations are summarized in Table 2. The fitness function, representing the absolute difference between the real and by GAapproximated response is depicted in Fig. 7. Constant value of the fitness function during the defined number of generations (NG) in this case represents reaching the global minimum of the function. From Fig. 7 it is obvious that the global minimum has been reached after less than 1000 generations and it remained constant until reaching the ending condition of 6000 generations, which serves as the indicator of reaching the optimal (best possible) quality of the fit (Table 3). Figs. 8 and 9 depict the part of VR-1 transfer function amplitude response close to the boundary period of 1.5 s (u ¼ 4.18 rad/s). Further decrease of the period (increase in angular frequency) should represent approaching the area of intermediate reactivity changes where the amplitude response remains constant until it reaches the other period boundary of approximately 0.06 s (u ¼ 104.72 rad/s). Fig. 9 depicts the same dependence in the case of beff and shows its convergence to the real value. The true value of beff for VR-1 reactor has not yet been discussed. The exact value has not been found by these aforementioned measurements. Further research has to be done in this field. The pneumatic system must be made capable of performing the required harmonic oscillations with a period from the area of the intermediate reactivity changes, i.e. in the interval <0.06 O 1.57> s (Fig. 1). Then, it can be reliably stated that the value of beff for the training reactor VR-1 is higher than the currently used value of 0.00714 and it asymptotically reaches the value, which lies in the interval of <0.0078 O 0.0081>. 4. Verification by Monte Carlo calculation The calculation was performed using MCNP5 code. The model of VR-1 reactor, containing both reactor tanks with all inner reactor components, excluding those summarized below and including the heavy-concrete biological shielding together with surrounding air. Fig. 10 depicts the VR-1 actual reactor core composition and Figs. 11 and 12 the MCNP plots of the reactor.
740
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742
Fig. 7. Fitness function for GA fitting e velocity 900 mm/s, LCM.
In the model, several simplifications were considered: Table 3 Influence of XS library selection in MCNP criticality calculations for beff determination. XS library JENDL-3.3 kp kt
beff
ENDF/B-VI
ENDF/B-VII
0.99244 0.00005 0.99792 0.00005 0.99264 0.00005 1.00016 0.00005 1.00636 0.00005 1.00036 0.00005 0.007719 0.000100 0.008387 0.000100 0.007717 0.000100
a) Pneumatic system (“Hopik”) is not included e model of fuel assembly is used instead, b) Fission chambers RJ1300 (M1-4) for power measurement during operation are not considered in the model, c) SNM detectors are not considered in the model. The calculation was done using Bretscher’s prompt method (Bretscher, 1997). This method consists of two steps. In the first
Fig. 8. Part of the VR-1 training reactor transfer function amplitude response measured by LCM310 detection system.
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742
741
Fig. 9. Dependence of beff on the oscillation period measured by LCM310 detection system. (beff values depicted here were calculated as ßeff ¼ jGðjuÞj1 , although the oscillation periods do not come under the area of intermediate changes.)
step, KCODE calculation of the effective multiplication factor for prompt neutrons (kp) was determined as kp ¼ 0.99244 0.00005. In the second, the value of total effective multiplication factor (kt) was calculated as: kt ¼ 1.00016 0.00005. The final effective delayed neutron fraction for VR-1 reactor was determined from Eq. (9):
beff y1
kp y0:007719 0:000100 kt
(9)
For this calculation, Japanese JENDL-3.3 XS library was chosen as a reference. The selection of cross-section library influences the results of the MCNP criticality calculations. To find the most suitable XS library for the calculation, the influence of XS data has been investigated and is described hereinafter. Also, calculations with the U.S. ENDF/B-VI and B-VII were made. The difference between
beff values based on the use of ENDF/B.VI or JENDL-3.3 libraries, respectively, calculated as 7.3%:
sbeff ¼
jbeff ENDF beff JENDL j
beff JENDL
$100:
(10)
This analysis confirms the results of the benchmark presented by Dos Santos et al. (2006), according to which the ENDF/B-VI library overpredicts beff by more than 7% in beff, beff/L and L calculations in comparison with JENDL-3.3. According to Mosteller (2008), “ENDF/B-VII.0 produces significantly better overall results for the cases in the MCNP criticality validation suite than do ENDF/B-VI, JENDL-3.3, or JEFF-3.1. However, further improvements still are needed in some areas.”, it has been concluded, that the level of ENDF/B-VII accuracy exceeded the JENDL-3.3 library. In our calculations for the case of
Fig. 10. VR-1 reactor core composition C1.
742
S. Michálek et al. / Progress in Nuclear Energy 52 (2010) 735e742
MCNP calculated beff ¼ 0.007719 0.000100 overvalues the currently used value by approximately 8.1%. In comparison with the MCNP value, the measurement provides a value overestimated by 1.76%. The 3% difference would be reached by beff ¼ 0.007957. The true value of beff must be determined by further measurements with oscillation periods of approximately 1 s. The performed measurements presented in this work point to the fact that there is a very high probability, that the true value of beff lies in the interval of <0.0078 O 0.0081>. Together with computational verification, they indicate that the currently used value has been underestimated and should be changed. Since the currently used value is based on a higher level of conservatism; this does not require immediate action in terms of reactor operation and safety. Any action should be postponed until measurements with lower reactivity oscillation periods are performed. More accurate results would be achieved by substituting the pneumatically realized reactivity insertions for a linear stepper motor. This would allow harmonic sinusoidal oscillations instead of trapezoidal with inconstant time delay of the absorber in the boundary positions. This was the limitation of the current pneumatic system-driven reactivity oscillations.
Fig. 11. MCNP5 code axial plot of whole reactor.
Acknowledgements The authors would like to thank the VR-1 reactor operational staff from Department of Nuclear Reactors of CTU in Prague for their help and valuable advice during the experiments.
References
Fig. 12. MCNP5 code radial plot of whole reactor.
VR-1 reactor, the difference between the final beff values based on JENDL-3.3 and ENDF/B-VII libraries, is only w0.03%, which leads to the conclusion, that both XS libraries can be used for these calculations, providing almost the same results. 5. Conclusion In this paper, the methodology of the effective delayed neutron fraction (beff) determination in reactor VR-1 has been presented. This methodology serves as the first experimental determination of this important reactor-safety parameter. The contemporary estimated value of beff ¼ 0.00714 is used in reactor operation. If the value for the fastest oscillation cycle performed, beff ¼ 0.007857, is considered, the condition of Eq. (11) published by Dos Santos et al. (2006), describing the difference between calculated (C) and experimental (E) value, is reliably fulfilled.
jC Ej $100 < 3% E
(11)
For the calculation, stochastic Monte Carlo method, realized by means of the MCNP5 code was applied for its verification. The
Bretscher, M.M., 1997. Evaluation of reactor kinetic parameters without the need for perturbation codes. In: Proc. Int. Meeting on Reduced Enrichment for Research and Test Reactors, Wyoming, USA, 5e10 October. Chiba, G., 2009. Calculation of effective delayed neutron fraction using a modified K-ratio method. Journal of Nuclear Science and Technology 46 (5), 399. Dos Santos, A., Diniz, R., Fanaro, L., Jerez, R., Silva, G., Yamaguchi, M., 2006. A proposal of a benchmark for beff, beff/ L, and L of thermal reactors fueled with slightly enriched uranium. Annals of Nuclear Energy 33 (9), 848e855. CSR Hermanský, B., 1987. Dynamics nuclear reactors. MS (in Czech). Irwanto, D., Obara, T., Chiba, G., Nagaya, Y., 2009. Study on calculation methods for effective delayed neutron fraction. In: Proc. 2009 Annual Meeting of the Atomic Energy Society of Japan, Tokyo, Japan, 23e25 March. Keepin, G.R., 1965. Physics of Nuclear Kinetics. Addison Wesley, Reading, Massachusetts, USA. Khamis, I., Hainoun, A., Suleiman, W., 2004. Measurement of the Syrian MNSR delayed neutron fraction and neutron generation time by noise analysis. Annals of Nuclear Energy 31 (3), 331e341. jka, K., Sklenka, L., Kropík, M., Kolros, A., 2003. Eugene Wiegner Training Mate Course at VR-1 Reactor (May 2003). Textbook for practical Course on Reactor Physics in Framework of European Nuclear Engineering Network, DNR FNSPE CTU Prague. p. 11. Mosteller, R.D., 2008. ENDF/B-VII.0, ENDF/B-VI, JEFF-3.1, and JENDL-3.3 results for the MCNP criticality validation suite and other criticality benchmarks. In: Proc. Int. Conf. Reactor Physics, PHYSOR2008, Interlaken, Switzerland, 14e19 September, LA-UR-07-6284. Nagaya, Y., Nakajima, K., 2008. Approximate estimation of effective delayed neutron fraction with correlated sampling method. In: Proc. Int. Conf. Reactor Physics, PHYSOR2008, Interlaken, Switzerland, 14e19 September. Nauchi, Y., Kameyama, T., 2005. Proposal of direct calculation of kinetic parameters beff and L based on continuous energy Monte Carlo. Journal of Nuclear Science and Technology 42, 503. Spriggs, G.D., Busch, R.D., Campbell, J.M., 2001. Calculation of the delayed neutron effectiveness factor using rations of k-eigenvalues. Annals of Nuclear Energy 28, 477. Van der Mark, S.C., Klein Meulekamp, R., 2006. Calculation of the effective delayed neutron fraction by Monte Carlo. Nuclear Science and Engineering 152, 142. Yedvab, Y., Grober, A., Ettedgui, H., Harari, R., Shvetz, I., Caspi, E., 2006. Determination of delayed neutrons source in the frequency domain based on in-pile oscillation measurements. In: Proceedings of PHYSOR-2006, Canadian Nuclear Society, Vancouver, Canada, 10e14 September.