First Application of Large Reactivity Measurement through Rod Drop Based on Three-Dimensional Space-Time Dynamics
Wencong Wang, Liyuan Huang, Caixue Liu, Han Feng, Jiang Niu, Qidong Dai, Guoen Fu, Linfeng Yang, Mingchang Wu
1 First Application of Large Reactivity Measurement through Rod Drop Based on Three-Dimensional Space-Time Dynamics
Wencong Wang, Liyuan Huang, Caixue Liu, Han Feng, Jiang Niu, Qidong Dai, Linfeng Yang, Mingchang Wu Nuclear Power Institute of China, China
Abstract
Reactivity measurement is an essential part of a zero-power physics test, which is critical to reactor design and development. The rod drop experimental technique is used to measure the control rod worth in a zero-power physics test. The conventional rod drop experimental technique is limited by the spatial effect and the difference between the calculated static reactivity and measured dynamic reactivity; thus, the method must be improved. In this study, a modified rod drop experimental technique that constrains the detector neutron flux shape function based on three-dimensional space-time dynamics to reduce the reactivity perturbation and a new method for calculating the detector neutron flux shape function are proposed. Correction factors were determined using Monte Carlo N-Particle transport code and transient analysis code for a pressurized water reactor at the Ulsan National Institute of Science and Technology and Xi'an Jiaotong University, and a large reactivity of over 2000 pcm was measured using the modified technique. This research evaluated the modified technique accuracy, studied the influence of the correction factors on the modification, and investigated the effect of constraining the shape function on the reactivity perturbation reduction caused by the difference between the calculated neutron flux and true value, using the new method to calculate the shape function of the detector neutron flux and avoiding the neutron detector response function (weighting factor) calculation.
Keywords:
Large reactivity measurement; Rod drop technique; Space-time dynamics; Constrained shape function; Monte Carlo N-Particle
1. Introduction
Reactivity measurement is an essential part of a zero-power physics test, which is vital for reactor design and development. In a zero-power physics test, different core configuration types are investigated, and the reactor power is limited to a low level; thus, the ex-core detector signal level is also low. Therefore, it is critical to determine the reactivity with a low-strength signal, especially when the reactivity is large. The rod swap and boron dilution methods can be used to measure the reactivity in a zero-power physics test. The rod swap method is slow, and it is easily influenced by the control rod shadow effect. The boron dilution method is accurate; however, it takes a long time to dilute the boron to compensate for the reactivity loss. This method is typically applied in commercial nuclear power plants, but it is rarely applied in a zero-power critical reactor because it is difficult to dilute boron in an experimental reactor. In recent decades, dynamic reactivity measurement methods have been widely used in commercial nuclear power plant startups. They separately and independently measure the worth of each control rod bank by inserting each control rod bank into the core at its maximum allowable speed. Dynamic methods, such as the dynamic rod worth measurement(DRWM) method developed by Chao, [1–3] dynamic reactivity measurement of rod worth method developed by Kastanya et al., [4] and the dynamic control rod reactivity measurement method developed by Lee et al., [5,6] have exhibited excellent results for numerous pressurized water reactor (PWR) startups. In these studies, the three-dimensional (3D) space-time kinetics theory is employed to overcome the limitations of a one-point reactor model in dynamic measurement. The reactivities measured by inserting and withdrawing the control rod bank at its maximum allowable speed in a commercial nuclear power plant were all less than 2000 pcm. Recently, Sang Ji Kim et al. [7] presented a preliminary dynamic rod drop simulation study for a transuranic burner core mockup of a sodium-cooled fast reactor, in which the reactivities were approximately 1000 pcm. However, there were several difficulties in applying these methods in a zero-power experimental reactor. (1) The ex-core detector signal level was lower than that of a commercial nuclear power plant due to the small core configurations. Hence, with the control rod inserted at its maximum allowable speed, the detector signal can easily fall below the lower limit of the measurement. (2) In certain cases, the reactivity of a single control rod is much larger than that of a commercial nuclear power plant control rod bank. For such a large reactivity insertion, the neutron flux distributions in the core and ex-core detector change rapidly and substantially, which may increase the difference between the calculated neutron flux and its true value, affecting the accuracy of the correction factors and reactivity results. Due to these difficulties, the DRWM method [1–9] is not always the optimal choice for the zero-power physics test of a zero-power experimental reactor. Thus, in certain cases, the rod drop experimental technique is used to measure the reactivity in a zero-power physics test. [10] The rod drop experimental technique measures the reactivity by dropping the control rod into the core, and the reactivity is analyzed using a one-point reactor model, which assumes that the neutron flux distribution in the reactor core maintains the same shape during the dynamic measurement. However, for large reactivity insertions, the neutron flux distribution changes greatly during the dynamic process, which can significantly affect the detector signal. Thus, the accuracy of a large reactivity measured during this process is unsatisfactory, indicating that the rod drop experimental technique must be modified. In this study, a modified rod drop experimental technique was applied to a large reactivity insertion measurement to reduce the measurement error. We propose a method that constrains the shape function of the detector neutron flux based on 3D space-time dynamics [11] to reduce the reactivity perturbation caused by calculated correction factor errors, which is affected by the difference between the calculated neutron flux and its true value. We also propose a new method to calculate the shape function of the detector neutron flux that avoids the neutron detector response function (weighting factor) calculation. Finally, we evaluate the proposed method accuracy, influence of the correction factors on the modification, and effect of constraining the shape function of the detector neutron flux. Compared to the previous research, our work is innovative in that we are the first to apply the modification of the rod drop experimental technique to a large reactivity measurement of over 2000 pcm in a zero-power experimental reactor, to constrain the shape function to reduce the reactivity perturbation caused by the calculated correction factor errors, and to propose a new method to calculate the detector neutron flux shape function without calculating the neutron detector response function (weighting factor). The difficulties of the conventional rod drop experimental technique and the modified method theory are discussed in Section 2. The proposed method application to large reactivity measurement is introduced in Section 3. The results and are described and discussed in Section 4. The conclusions are summarized in Section 5.
2. Methodology [11]
Thus, in this research, the inverse kinetics based on a one-point reactor model was employed to analyze the data measured by dropping a control rod into a core, which we defined as the conventional rod drop experimental technique to make more easily distinguish the method before and after modification. The measurement procedure is as follows: (1)
Operate the reactor core in the critical state. (2) Drop the control rod into the core while recording the neutron signal with the ex-core neutron detector. (3)
Process the neutron signal using reactivity measurement equipment. (4)
Calculate the reactivity using the experiment data and inverse kinetic equation. The inverse kinetic equation can be described as follows: [12] ( ) 1 ( )( ) ( ) i t ti i ti dN t N e dN t dt N t , (1) where is the reactivity, is the effective fraction of delayed neutron, is the neutron generation time, and is the delayed neutron decay constant. , , and are calculated before the experiment, and they are considered to be known constant parameters in the measurement. The ENDF/B7.0 cross-section libraries were used to calculate and in this research. N ( t ) is the neutron detector signal, which was assumed to be proportional to the amplitude function p ( t ); therefore, the amplitude function p ( t ) was replaced with N ( t ) in Equation (1). There are certain difficulties in using the conventional method, such as the spatial effect of the neutron detector signal and difference between static reactivity and dynamic reactivity. (1) Spatial effect of the neutron detector signal In a one-point reactor model, the time dependence of the flux ( , , ) r E t is assumed to be separable from its space and energy dependences, and the flux can be described as the product of amplitude function p ( t ) and shape function ( , ) r E , assuming that the shape function is unchanged during measurement: ( , , ) ( ) ( , ) r E t p t r E . (2) With this assumption, the shape function ( , ) r E is independent of time, and the neutron flux ( , , ) r E t is proportional to the amplitude function p ( t ) in the measurement; therefore, the detector neutron flux ( , , ) d r E t is proportional to the amplitude function p ( t ). In a real situation, the shape function of the flux is time-dependent, and it changes significantly during the rod drop process. Thus, the detector neutron flux is affected by changes in the neutron flux shape. This violates the assumption that the shape function is unchanged during measurement. For this reason, the detected neutron signal is not proportional to the neutron flux amplitude function, and using it to determine the reactivity leads to measurement inaccuracies. (2) Difference between static reactivity and dynamic reactivity The reactivity determined by the inverse kinetic equation is called the dynamic reactivity. It is definitional different from the reactivity determined by static calculation, which is called the static reactivity. Thus, to obtain superior results, it is important to remedy these difficulties in rod drop reactivity measurement. 2.2 Modified rod drop experimental technique The modified rod drop experimental technique is based on exact point dynamics, which are the 3D reactor dynamics. The modification is focused on two aspects: (1) the difference between the neutron flux amplitude function and neutron detector signal and (2) the difference between the calculated static reactivity and measured dynamic reactivity. 2.2.1 Detector signal correction 2.2.1.1 Detector signal correction factor From the exact point dynamics, the neutron flux in the reactor core ( , , ) r E t can be factorized into a purely time-dependent amplitude function, p ( t ), and a space-, energy-, and time-dependent neutron flux shape function ( , , ) r E t : [11] ( , , ) ( ) ( , , ) r E t p t r E t . (3) Based on the exact point dynamics, the neutron signal ( ) Det
N t is related to the neutron flux in the reactor core ( , , ) r E t : [11] ( ) ( , ) ( , , ) Det V E
N t W r E r E t dEdV , (4) where ( , ) W r E is the weighting factor that denotes the neutron contribution degree at position r in the core region to the neutron detector signal, describing the space and energy dependences of the neutron detector sensitivity. Thus, Equation (5) can be written as follows: ( ) ( , ) ( ) ( , , ) Det V E
N t W r E p t r E t dEdV . (5) The neutron flux amplitude function p ( t ) can be extracted from the neutron signal using the following equation: ( )( ) ( , ) ( , , ) DetV E
N tp t W r E r E t dEdV . (6) Using Equation (6), the neutron signal is converted into the neutron flux amplitude function p ( t ), and the difference between the neutron flux shape and neutron signal can be largely reduced, resulting in a more accurate result. The neutron signal ( ) Det
N t can also be described as follows: ( ) ( ) ( , ) ( ) ( ) ( , )
Det d d
N t t r t t p t r t , (7) where ( ) t is the detector sensitivity. From Equations (6) and (7): ( , ) ( , , ) ( ) ( , ) dV E W r E r E t dEdV t r t . (8) Thus, the denominator on the right side of Equation (6) is considered to be proportional to the shape part of the detector neutron flux ( , ) d r t . To modify the detector signal, one must calculate the denominator on the right side of Equation (6). The neutron detector response function (weighting factor) ( , ) W r E can reduce the accuracy due to the mesh precision and the assumption of no changes under different control rod patterns, and the denominator on the right side of Equation (6) is proportional to the shape part of the detector neutron flux ( , ) d r t . Thus, we use the shape part of the detector neutron flux ( , ) d r t to describe the denominator on the right side of Equation (6) to reduce the inaccuracy induced by the neutron detector response function (weighting factor) ( , ) W r E . The new method for calculating the shape part of the detector neutron flux ( , ) d r t is described in detail in Section 3.2.5. Normalization about the detector signal is performed as follows: det det ( ) ( ) (0) 1 ( ) 1( ) (0) (0) ( ) (0) Det DetDet Det p t N t N tp t p N t C N C . (9) Because the calculated energy spectra of the detector neutron flux stay the same during the rod drop process, ( ) t is assumed to be unchanged in the measurement. det ( , )( , 0) dd r tC r , (10) where C det is the detector signal correction factor, and “0” refers to the critical state at the beginning of the rod drop process. By substituting the normalized amplitude function ( ) p t into Equation (1), instead of the neutron detector signal N ( t ), the reactivity is obtained. 2.2.1.2 Constraining the shape function To modify the detector signal, the shape part of the detector neutron flux ( , ) d r t should be calculated. The purpose of constraining the shape function is to reduce the reactivity perturbation caused by the neutron flux or shape function calculation error. We begin with the definition of ( , , ) r E t , which is the neutron flux shape function. From the 3D neutron diffusion equation: p d r E t F M r E t S r E t S r E tv t , (11) where M is the neutron destruction operator, F p is the prompt neutron production operator, ( , , ) d S r E t is the delayed neutron source, and ( , , )
S r E t is an independent source. By substituting Equation (3) into Equation (11), we obtain the following equation: p d p t r E t p t r E t F M p t r E t S r E t S r E tv . (12) Allowing the shape function to depend on time is a first generalization compared to the use of the time-independent shape in the derivation presented in Equation (2). A second generalization used in the derivation does not remove an approximation, but rather exploits a certain freedom of choice; the neutronics equation is multiplied by a weight function, ( , ) w r E , prior to integration with respect to space and energy. The flux factorization is also introduced into the left side of equation (11). [11] ( ) ( , ) ( , , ) ( , ) ( , , )( )( ) ( )( ) ( ) ( , , ) ( , , ) ( , , ) w wV E V Ew w wp dV E V E V E dp t r E r E t d r E r E tdEdV p t dEdVdt v E dt v EF M p t r E t dEdV S r E t dEdV S r E t dEdV . (13) The second term on the left side of Equation (13) that appears with flux factorization can be eliminated by only using the integral to constrain the time variation of the shape function, thus making the factorization unique: ( , ) ( , , )( ) wV r E r E t dEdV Cv E , (14) where C is an arbitrary constant. Constraining the neutron flux shape using Equation (14) does not introduce an approximation, and the factorization introduces a new degree of freedom. With the constraint defined as Equation (14), Equation (13) performs the same as the one-point model without any assumptions or approximations. Thus, the constraint is vital in determining the neutron flux shape function ( , , ) r E t . The ( , ) w r E choice is also important because it can affect the ( , , ) r E t accuracy and thus affect the detector signal modification. In practice, the shape function ( , , ) r E t cannot be precisely determined (the calculation result is always somewhat different from the true value), and the calculated flux and flux amplitude p ( t ) are therefore only approximate solutions. Both values depend on the weight function. It is thus advantageous to choose a weight function that reduces the error resulting from inaccuracies in the shape function. Because the solution of the point kinetics equation is particularly sensitive to reactivity errors, a weight function should be selected that reduces the effect of shape function inaccuracies on the reactivity. The initial adjoint flux, *0 ( , ) r E , fulfills this objective. [11] This is because, in static perturbation theory, the first order dominates the reactivity perturbation, which is caused by the perturbation in neutron flux, and therefore, the use of *0 ( , ) r E can eliminate the first order in the reactivity perturbation. Thus, the use of *0 ( , ) r E can eliminate most of the reactivity perturbation caused by the difference between the calculated neutron flux and true value, reducing the reactivity perturbation caused by the calculation error of the neutron flux or shape function. It is more precise to use * ( , , ) r E t as the weight function ( , ) w r E . However, ( , ) w r E should remain unchanged in the measurement to maintain the factorization consistency (Equation (3)). Thus, the use of *0 ( , ) r E is the optimal choice because it can eliminate the reactivity perturbation caused by the error in the neutron flux calculation. In summary, to modify the detector signal, it is essential to impose a constraint on the shape function, making the factorization unique with the initial adjoint flux *0 ( , ) r E as the weighting function. *0 00 ( , ) ( , , )( ) V r E r E t dEdV Kv E , (15) where *0 ( , ) r E is the adjoint neutron flux distribution of the critical state at the beginning of the rod drop measurement, v ( E ) is the neutron velocity, and K is an arbitrary constant, which is 1 in this study. 2.2.3 Dynamic reactivity correction The reactivity determined by the inverse kinetic Equation (1) is typically considered the dynamic reactivity, which differs from the static reactivity. To remedy the difference between calculated static reactivity and measured dynamic reactivity, a dynamic reactivity correction is performed as follows: , , st m dyn dyn m C , (16) ,, st cdyn dyn c C , (17) where , dyn m is the measured dynamic reactivity, , st c is the calculated static reactivity, and , dyn c is the calculated dynamic reactivity, which is determined by substituting the calculated neutron flux amplitude function into the inverse kinetic equation. To perform the correction, the difference between the calculated static reactivity and measured dynamic reactivity is determined via theoretical calculation. The static reactivity , st c is determined using a static eigenvalue calculation, and the neutron flux amplitude function is simulated using a transient calculation code. By substituting the simulated neutron flux amplitude function into Equation (1), the calculated dynamic reactivity , dyn c is obtained. Then, the dynamic reactivity correction factor is determined by incorporating the calculated static reactivity and calculated dynamic reactivity into Equation (17). In this research, the static reactivity is calculated using Monte Carlo N-Particle (MCNP) transport code [13] , and the neutron flux amplitude function is simulated transient analysis code for a pressurized water reactor at the Ulsan National Institute of Science and Technology and Xi'an Jiaotong University (TAPUX) [14] .
3. Application B neutron detector AAcoreneutron detector B
Fig. 1. Schematic of experimental core configuration. In this study, the reactivity measurement of a zero-power water-moderated thermal critical reactor was performed. Strong reactivity was locally added into the core by dropping a single control rod cluster into the core from its full out state when the reactor operates in a critical state. Two control rod clusters ( (1) Build 3D MCNP and 3D dynamic calculation models. (2)
Calculate the static physical parameters. (3)
Constrain the detector neutron flux into the shape function, and modify the measured neutron signal. (4)
Calculate the amplitude function p ( t ) of the dynamic rod drop process. (5) Obtain the dynamic reactivity correction factor to determine the measured final reactivity. Fig. 2. Flow diagram of modified rod drop experimental technique. 3.1 Calculation model
Table 1. Results of experimental critical control rod pattern calculations using MCNP
Number Experimental critical control rod pattern
MCNP keff ± error 1 ± ± ± The calculation geometry is based on Fig. 1. Both the MCNP and dynamic calculation models used for modification in this research were verified by experimental critical control rod pattern. The results of the experimental critical control rod pattern calculation by MCNP are presented in Table 1. The ENDF/B-VI cross-section libraries were used in the MCNP code. The detector was located outside the core, and the calculation result of the detector neutron flux energy spectra remained the same before and after the rod drop. Therefore, the neutron detector efficiency was assumed to be constant in the rod drop process. Because the shape function at detector position ( , ) d r t appears a as ratio in Equation (10), the neutron detector in the model was simplified as a cylinder with its shell, and each detector was located in a tube. The detector was located approximately 10 cm above the bottom of the active core in the axial direction. The sensitive height and radius of the detector were approximately 36 cm and 2.5 cm, respectively. In this study, all of the static physical parameters were determined using the MCNP code, such as the detector neutron flux ( , ) d r t , neutron velocity v , neutron flux ( , , ) r E t , and adjoint neutron flux distribution in the core *0 ( , ) r E . The dynamic physical parameter, which is the amplitude function p ( t ) in the rod drop process, was calculated using TAPUX. 3.2 Calculation of static physical parameters 3.2.1 Detector neutron flux The detector neutron flux ( , ) d r t was calculated using an MCNP neutron flux tally card based on the previously mentioned MCNP calculation model. The geometry split technique was used to improve the calculation efficiency. The calculated detector neutron flux ( , ) d r t was constrained into a shape function at detector position ( , ) d r t , which is elaborated on in Section 3.2.5. 3.2.2 Neutron flux distribution The neutron flux distribution calculation in the core ( , , ) r E t was directly performed using the neutron flux tally card. The volume tally mesh was in the XY plane assembly-wise and 10 cm in the Z-direction. 3.2.3 Neutron velocity The calculation of neutron velocity v in the core was realized using the tally energy card and MCNP neutron flux tally card. The volume tally mesh was the same as that of the neutron flux distribution calculation. Using the calculated ( , , ) r E t of each energy bin in each volume mesh and the relationship between energy and velocity, the neutron velocity v was calculated as follows: , , ,1 1 ,1 1
2( )
V E V E n n i up i down i jj i n n i jj i
E EE t , (18) E tv t m . (19) The average neutron energy in core E was determined by the upper and lower energy boundaries and the corresponding neutron flux of each energy bin in the volume tally mesh elements. n E is the number of energy bins in a volume tally mesh element, and it was the same for all of the n V volume tally mesh elements. 3.2.4 Adjoint neutron flux To constrain the detector neutron flux ( , ) d r t into a shape function at the detector position ( , ) d r t , the adjoint neutron flux distribution in the core *0 ( , ) r E was calculated. The proportionality of the iterated fission probability ( I FP ) to the adjoint flux was demonstrated [15] , and I FP was therefore calculated instead of the adjoint neutron flux *0 ( , ) r E , which was difficult to obtain using the MCNP code. I FP was determined using the formulas below [15,16] : (1) (2) ( 1) ( )( ) 0 (0) (1) ( 2) ( 1) ( ) FP s s s sI s s s s , (20) ( ) (1) (2) ( 1) ( )0 ( ) FP I k k k k . (21) As the generation number, , increases, distributions of the neutron flux ( ) and the fission neutron emission (source) produced by the progenies converge to those of the fundamental mode. ( ) i s is the number of fission neutrons of the i-th generation; hence, the corresponding ratio ( ) ( 1) / i i s s can be written as ( ) i k , which can be estimated by the eigenvalue calculation in MCNP. [15] is the initial point source of I FP , and the I FP of ( , , ) r E can be estimated by calculating the initial source points at . The I FP was calculated using Equation (21). In the I FP calculation, ( ) i k was determined using k-collision, which was printed in the MCNP output file. The I FP of every tally mesh element was determined using a corresponding eigenvalue calculation. In every eigenvalue calculation, sufficient neutrons were sampled by locating hundreds of homogenous source points in the corresponding mesh grid. 3.2.5 Shape function at detector position ( , ) d r t The new method for calculating the shape part of the detector neutron flux ( , ) d r t is described in detail in this section. ( , ) d r t was determined using the detector neutron flux ( , ) d r t and constraint factor of the shape function C ( t ). For the static physical parameters calculated using the methods described in Sections 3.2.2–3.2.4, the neutron flux in the core ( , , ) r E t was constrained into shape function in the core ( , , ) r E t as follows: ( , , )( , , ) ( ) r E tr E t C t , (22) where C ( t ) is the constraint factor of the shape function: *00 ( , ) ( , , )( ) ( ) V r E r E tC t dEdVv E . (23) The shape function satisfies the constraint condition Equation (15) when K =1. The detector neutron flux ( , ) d r t was constrained into shape function at detector position ( , ) d r t as follows. From Equation (8): ( , ) ( , ) ( , , ) / ( ) d V E r t W r E r E t dEdV t . (24) By substituting Equation (22) into Equation (24): ( , ) ( , , ) / ( )( , ) ( ) V Ed
W r E r E t dEdV tr t C t . (25) Thus, ( , ) d r t can be described as ( , )( , ) ( ) dd r tr t C t . (26) *00 ( , )( , ) ( , ) ( , , )( ) dd V r tr t r E r E t dEdVv E . (27) As displayed in Equation (27), the shape function at detector position ( , ) d r t can be determined using the detector neutron flux ( , ) d r t and constraint factor of the shape function C ( t ). Hence, the weighting factor ( , ) W r E calculation is avoided. This method appears to be superior to that calculating the neutron detector response function (weighting factor) ( , )
W r E [17,18] because ( , )
W r E can reduce the accuracy due to the mesh precision and assumption that it does not change under different control rod patterns. 3.3 Dynamic physical parameter calculations The amplitude function p ( t ) was determined using TAPUX and the dynamic model mentioned above. The TAPUX code is a recently developed 3D two-group light water reactor core analysis code by the Nuclear Engineering Computational Physics laboratory at Xi’an Jiaotong University. It can be used by utilities to perform transient analysis in neutronics. The code adopts the non-linear coarse-mesh finite difference method based on the nodal methodology in steady-state and transient core calculations. The frequency transform method was applied based on the θ method in time discretization. Thus, the time step can be expanded to enhance the efficiency. The cross-section of TAPUX was derived from the Bamboo lattice code [19–21] , which was based on the ENDF/B7.0 library and 69-group lattice calculation. A two-group homogenized cross-section was generated for each assembly and made into a look-up table for the transient calculation in TAPUX considering the fuel temperature feedback, coolant density, and control rod movements. For the calculation, the control rod was dropped into the core for free fall, which was the same as that in a real experiment situation. Thus, the amplitude function p ( t ) in the rod drop process was determined. 3.4 Correction factor calculations Fig. 3. Detector signal correction factor. The shape function at detector position ( , ) d r t was determined using the static physical parameters, which were calculated using the MCNP code. The detector signal correction factors in multiple rod positions were calculated using Equation (10), and the relationship between the correction factor and rod position was obtained by high-order polynomial fitting. The detector signal correction factor curve is illustrated in Fig. 3. Then, the raw signals N Det ( t ) were normalized and modified using Equation (9). We can see from the figure that the spatial effects of the Table 2. Dynamic reactivity correction factors Measured rod cdyn , (pcm) cst , (pcm) C dyn By substituting the amplitude function p ( t ) calculated with TAPUX into Equation (1), the calculated dynamic reactivity , dyn c was obtained. With the static k eff calculated by MCNP, the static reactivity , st c was determined. In this research, we focused on the integral rod worth; hence, the corresponding dynamic reactivity correction factor was calculated using Equation (17), and the results are listed in Table 2. The dynamic reactivity correction factors are close to 1, which means that the dynamic effect is small during the process. 3.5 Uncertainty Analysis In references [22] and [23], the uncertainty was estimated using the standard deviation, and a first-order Taylor formula was used to linearize and calculate an approximation of the uncertainty. Based on the analytical method of uncertainty in references [22] and [23], the uncertainty of the detector signal correction factor C det and the dynamic reactivity correction factor C dyn were determined. Then, the uncertainties were propagated to calculate the final reactivity. 3.5.1 Detector signal correction factor C det uncertainty According to the uncertainty definition in reference [23], the relative uncertainty was evaluated as the relative standard deviation. For the MCNP output, the relative error was evaluated as the relative standard deviation. Thus, the relative uncertainties of the detector neutron flux ( , ) d r t and neutron flux distribution ( , , ) r E t were evaluated using the relative error in the MCNP output file. The relative uncertainties of the physical parameters were propagated to the detector signal correction factor C det . Using the analytical method of the standard deviation [22,23] , the relative uncertainty of the adjoint neutron flux distribution *0 ( , ) r E was determined.
2( )* ( ') ' ,0, ,* ( ') ( )10, , , ( )( ( )) ( )( ) imm FP m iim FP m m u ku E u IE I k . (28) For each criticality calculation, we assumed that the relative uncertainties of ( ) i k from different cycles were the same because the calculation condition was unchanged even though the source points changed across cycles. Here, ( ) ( ) ( ) / i i u k k was assumed to be the same in one critical calculation. Then, we obtained the relative uncertainty of I FP , [16] ( )( ') ,, ( ') ( ), , ( )( ) ' imFP m iFP m m u ku II k . (29) As the neutrons propagated, the fission neutron emission (source) produced by the progenies converged to the fundamental mode. The k calculated in these fundamental mode source condition was considered to be under the same simulation conditions, and it was used to determine the relative uncertainty of k . [16] Therefore, we skipped 50 cycles and used the k of the last 50 cycles to calculate the relative uncertainty of *0 ( , ) r E . 3.5.2 Uncertainty of dynamic reactivity correction factor C dyn According to the uncertainty of k eff , which is illustrated by the MCNP results, the uncertainty of the dynamic reactivity correction factor C dyn was determined. Using Equation (17), the relative uncertainty of C dyn can be described as follows: ( ) ( )( ) dyn effrel dyn dyn eff u C u ku C C k . (30) Here, the uncertainty of TAPUX was ignored because it is a deterministic code. 3.5.3 Uncertainty of final reactivity modification Based on the above calculations, the relative uncertainty of the detector signal correction factor C det was obtained, and the relative uncertainty curve of the detector signal correction factor was obtained with the relative uncertainties at different rod positions by fitting. Then, the relative uncertainty of the detector signal correction factor C det and the dynamic reactivity correction factor C dyn were propagated to the reactivity. Using Equations (1), (9), and (16), the relative uncertainty caused by the modification in the reactivity was determined, as illustrated in Table 3. Fig. 4. Normalized detector signals before and after modification. According to the measurement procedure described in Section 2.1, the detector signal was obtained through experimentation. Because the control rod position changes rapidly during the rod drop [24–26] , the sample frequency was set to 100 Hz in the measurement. Neutron detector A measured the neutron signal of control rod cluster N (0), and the data were transformed into the relationship between the signal and relative control rod position using the relationship between time and the relative control rod position during the free fall rod drop process. The normalized detector signal N ( t )/ N (0) is displayed in Fig. 4 (raw signal), where the relative control rod position ‘1’ means that the rod is at the top of the core, and the relative control rod position ‘0’ means that the rod is at the bottom of the core. Using the signal modification procedure in Section 3, the detector signal correction factor was determined, and the modified detector signal was obtained. The normalized modified detector signal N c ( t )/ N c (0) is displayed in Fig. 4 (modified signal). As illustrated in this figure, the detector signal falls rapidly during the rod drop process, which is less than 1 s, and the signal decreases by approximately one order of magnitude for the reactivity was obtained. The integral rod worth calculated with both raw and modified data is listed in Table 3. 4.1.2 Results based on dynamic reactivity modification The integral rod worth both with and without signal modification were modified with dynamic reactivity correction factor C dyn using Equation (16) are also displayed in Table 3. 4.1.3 Final results With both detector signal and dynamic reactivity modifications, the final results were obtained, as exhibited in Table 3. The measurement result is compared with the MCNP result, and the relative difference with the MCNP result is listed in the table. Table 3. Comparison of modified integral rod worth Results Measured rod ± Uncertainty (pcm) −4874.9 ± −2664.2 ± Conventional rod drop experimental technique (with raw data) Reactivity (pcm) −5604.4 −2803.3 Relative difference 14.96% 5.22% Modified results 1 (with signal modification only) Reactivity ± Uncertainty (pcm) −4288.1 ± −2733.1 ± Relative difference −12.04% 2.59% Modified results 2 (with dynamic reactivity modification only) Reactivity ± Uncertainty (pcm) −5793.2 ± −2655.6 ± Relative difference 18.84% −0.32% Final reactivity (with signal and dynamic reactivity modifications) Reactivity ± Uncertainty (pcm) −4432.5 ± −2589.1 ± Relative difference −9.07% −2.82%
As the results demonstrate, the discrepancy between the conventional rod drop experimental technique and MCNP is large because of the spatial and dynamic effects in the conventional measurement, while most of the modified results are improved. As we discussed in Section 2.2.1.2, the purpose of constraining the shape function is to reduce the reactivity perturbation caused by the calculation error of the neutron flux or shape function. The raw signal was modified using the detector neutron flux ( , ) d r t , which was calculated by MCNP without constraining the shape function, and the rest of the modifications are the same as those previously mentioned. Then, the final reactivities both with and without constraining the shape function were compared. The investigation was based on Table 4. Effect of constraining shape function in ± Uncertainty (pcm) −4874.9 ± −5447.8 ± Conventional rod drop experimental technique Reactivity (pcm) −5604.4 −5869.5 Relative difference 14.96% 7.74% Final reactivity (without constraining the shape function) Reactivity (pcm) ± Uncertainty (pcm) −4373.4 ± −5231.0 ± Relative difference −10.29% −3.98% Final reactivity Reactivity (pcm) −4432.5 −5301.8 (with constraining the shape function) ± Uncertainty (pcm) ± ± Relative difference −9.07% −2.68%
The results of the clear improvement compared to the results of single modification (modified results 1 and 2). The final reactivity of the In this study, a modified rod drop experimental technique was proposed, and it was applied to large reactivity measurement. In the modification, static physical parameters were calculated using the MCNP code, and dynamic physical parameters were calculated using a transient code. The reactivity of the rod drop process, in which large reactivity (approximately −5500 pcm) is locally inserted, was measured using the modified rod drop experimental technique. The primary conclusions can be drawn from the results as follows: (1)
When the large reactivity is locally inserted in the rod drop, the conventional rod drop experimental technique accuracy is limited by the spatial effect and the difference between the static reactivity and dynamic reactivity. The modified rod drop experimental technique can reduce the spatial and dynamic effects in the measurement, and it is more accurate and valid. (2)
The detector signal modification is important and essential, and it dominates the modification of large reactivity measurement. Modifications based on MCNP exhibit satisfactory results. (3)
The dynamic reactivity modification is necessary for large reactivity measurement. (4)
Constraining the shape function can reduce the reactivity perturbation caused by the difference between the calculated neutron flux and its true value, and the results suggest that the modification can improve the results.
Acknowledgments The authors would like to thank Guoen Fu, Dong Wei, Yongmu Yang, Shibiao Pan, Cuiyun Peng, Yonglin Zhang, Yan Guo, Jie He, Hang Zhou, Aining Deng, and Yifan Zhang of the Reactor Engineering Research Sub-institute of the Nuclear Power Institute of China for their help in the experiment related to this work.
References Y. A. Chao, D. M. Chapman, L.R. Grobmyer, et al., Dynamic rod worth measurement (DRWM) next generation rod worth measurement method. In: Proc. Topl. Mtg. Advances in Nuclear Fuel Management II, Myrtle Beach, South Carolina, USA, March 23–26, 14–53(1997), ANS. Y.A. Chao, D.M. Chapman, D. J. Hill, et al., Dynamic rod worth measurement. Nucl. Tech. 132, 403(2000). DOI: 10.13182/NT00-A3153 Kastanya, D.F., Ariani, I., Turinsky, P.J., Verification of dynamic rod worth measurement calculation methodology. In: Proc. Topl. Mtg. Advances in Nuclear Fuel Management II, Myrtle Beach, South Carolina, USA, 23–26 March, 14–65(1997), ANS. http://refhub.elsevier.com/S0029-5493(15)00456-2/sbref0035 E. K. Lee, H. C. Shin, S. M. Bae, et al, Application of the dynamic control rod reactivity measurement method to Korea Standard Nuclear Power Plants. In: Proceedings of the PHYSOR 2004, Chicago, IL, April 25–29(2004), ANS. E.K. Lee, H.C. Shin, S.M. Bae, New dynamic method to measure rod worths in zero power physics test at PWR startup. Annals of Nuclear Energy, 32, 1457-1475(2005). doi:10.1016/j.anucene.2005.04.008 S.J. Kim, Ha, PNV, Lee, MJ, et al, Dynamic rod worth simulation study for a sodium-cooled TRU burner. Nuclear Engineering and Design, 295, p192-203(2015). DOI: 10.1016/j.nucengdes.2015.10.007 Hong, S.G., Song, J.S., A preliminary simulation study of dynamic rod worth for the SMART (System-integrated Modular Advanced ReacTor) reactor. Ann. Nucl. Energy, 60, 350-356(2013). DOI: 10.1016/j.anucene.2013.05.014 A. A. Bobrov, G. V. Lebedev, Yu. A. Nechaev, Measurements of Control Rod Worth by Modified Inverse Kinetic Method. Physics of Atomic Nuclei, 74(14), 1917–1920(2011). DOI: 10.1134/S106377881114002X Krell, J, Mangel, B, Rebohm, H, Analysis of Rod-Drop Measurements at VVER-440 Reactor. Nuclear Engineering and Design, 159, 265-271(1995). DOI:10.1016/0029-5493(95)01090-5 Karl O. Ott, Robert J. Neuhold. Introductory Nuclear Reactor Dynamics. (American Nuclear Society, Illinois, 1985), pp 55-58. Shi, YQ.
Neutronics Experimental Technology in Nuclear Reactor . (Atomic Energy Press, Beijing, 2011). X-5 Monte Carlo Team, MCNP-A General Monte Carlo N-Particle Transport Code, Version 5, LA-UR-03-1987, Los Alamos National Laboratory (LANL), 2003. Zheng, YQ, Lee, DK, A nodal code development for PWR transient analysis based on non-linear iteration method. High Power Laser and Particle Beams, 29, 036001(2017). DOI:10.11884/HPLPB201729.160297 Y.Nauchi, T.Kameyama, Development of Calculation Technique for Iterated Fission Probability and Reactor Kinetic Parameters Using Continuous-Energy Monte Carlo Method. Journal of Nuclear Science and Technology, 47(11), 977-990(2010). DOI:10.3327/jnst.47.977 Wang, WC, Liu, CX, Huang, LY, The first application of modified neutron source multiplication method in subcriticality monitoring based on Monte Carlo, Nuclear Engineering and Technology, 52 , 477-484(2020). DOI:10.1016/j.net.2019.08.014 Farkas,G., Lipka,J., Hascik,J., et al., Computation of ex-core detector weighting functions for VVER-440 using MCNP5. Nucl. Eng. Des. 261, 226–231(2013). DOI:10.1016/j.nucengdes.2013.01.026 Crump, M.W., Lee, J.C., Calculation of spatial weighting functions for ex-core detectors. Nucl. Technol. 41, 87–96(1978). DOI:10.13182/NT78-A32135 Li Y , Zhang B , He Q , et al, Development and verification of PWR-core fuel management calculation code system NECP-Bamboo: Part I Bamboo-Lattice. Nuclear Engineering and Design, 335, 432-440(2018). DOI:10.1016/j.nucengdes.2018.05.030 Yang W , Wu H , Li Y , et al, Development and verification of PWR-core fuel management calculation code system NECP-Bamboo: Part II Bamboo-Core. Nuclear Engineering and Design, 337, 279-290(2018). DOI:10.1016/j.nucengdes.2018.07.017 Li Y , He T , Liang B , et al, Development and verification of PWR-core nuclear design code system NECP-Bamboo: Part III: Bamboo-Transient. Nuclear Engineering and Design, 359, 110462(2020). DOI:10.1016/j.nucengdes.2019.110462 ISO/IEC GUIDE 98-3:2008(Uncertainty of measurement - Part 3: Guide to the expression of uncertainty in measurements (GUM)). http://refhub.elsevier.com/S1738-5733(19)30063-4/sref8 ISO/IEC GUIDE 98-3:2008(Uncertainty of measurement - Part 3: Guide to the expression of uncertainty in measurements (GUM)) Supplement 1: Propagation of distributions using a Monte Carlo method. http://refhub.elsevier.com/S1738-5733(19)30063-4/sref9 Kai WANG, Xiaowei JIAO, Qun YANG, Yanhua WANG Chaoqun WU, Zhaozhong HE. The effect of scram rod drop time on the consequences of molten salt reactor reactivity insertion transient[J]. Nuclear Techniques, 2020, 43(9): 90606-090606. DOI: 10.11889/j.0253-3219.2020.hjs.43.090606 Xungang CHEN, Jintao MO, Ying LUO, Dapeng YAN, Yueyin SHEN, Haoxuan NIU. Simulation study on the characteristics of driving line rod drop based on passive acceleration[J]. Nuclear Techniques, 2020, 43(7): 70604-070604. DOI: 10.11889/j.0253-3219.2020.hjs.43.07060426