Subcriticality measurement using time-domain decomposition-based integral method for simultaneous reactivity and source changes

ABSTRACT

To estimate the subcriticality in dollar units for an arbitrary state-change, the time-domain decomposition-based integral method (TDDI) is proposed using the point kinetics theory based on the fundamental mode approximation. In a general transient subcritical system, reactivity, neutron source intensity, and point kinetics parameters can vary simultaneously. Furthermore, the state-change may not necessarily be a stepwise change. For such a transient, the TDDI method can estimate the subcriticality after the transient using only the time variation of the neutron count rate. Therefore, the proposed method is useful to approximately estimate the subcriticality in a system where a detailed core configuration is unknown. To investigate the applicability of the TDDI method, transient experiments with simultaneous reactivity and source changes or to two successive safety rods dropping were performed at the Kindai University Training and Research Reactor (UTR-KINKI). By comparing with reference values using excess reactivity and control rod worth, it was validated that the subcriticality values obtained by the TDDI method better agree with the reference values than the previous integral method.

KEYWORDS:

• UTR-KINKI
• subcriticality
• transient measurement
• control rod
• neutron source
• point kinetics theory
• integral method
• random sampling method

1. Introduction

The subcriticality measurement technique is important to assure the safety margin from the critical state. Although various measurement techniques have been proposed to measure the subcriticality (or negative reactivity $-\rho \equiv \left(1-k_{\mathrm{e}\mathrm{f}\mathrm{f}}\right)/k_{\mathrm{e}\mathrm{f}\mathrm{f}}$), there are still problems to be addressed in applying the technique to a general transient situation, where not only $-\rho$ but also the neutron source intensity $S$ and point kinetics parameters (neutron generation time $\mathrm{\Lambda }$ and the effective delayed neutron fraction ${\beta }_{\mathrm{e}\mathrm{f}\mathrm{f}}$) can be simultaneously varied.

For example, consider the application of the subcriticality monitoring to the refuelling process at a boiling water reactor (BWR) in shutdown state [1,2]. During refuelling or handling of a control rod in BWR, operators confirm the magnitude of neutron flux in the low-power (startup) condition by the startup range neutron monitor (SRNM) [3]. Because the SRNM requires minimum neutron count rate to assure its intactness, burnt fuel assemblies containing inherent neutron sources (e.g. spontaneous fissions of ²⁴²Cm and ²⁴⁴Cm) or an external neutron source is loaded at first around a detector during refuelling. Thereby, $S$ and $-\rho$ can change simultaneously during the refuelling. In addition, the point kinetics parameters can also vary, because the change of fuel loading pattern yields the spatial variations in fission rate and importance function, and thereby affects the point kinetics parameters.

Another example is the subcriticality measurement in the removal of fuel debris at the Fukushima Daiichi nuclear power plants. The fuel debris includes the spontaneous fission nuclides such as ²⁴²Cm and ²⁴⁴Cm, which are also contributors for ¹³⁵Xe and ⁸⁸Kr detected by the primary containment vessel (PCV) gas monitoring system [4]. The removal process has an impact on both $\rho$ and $S$, i.e. fissile and neutron-source nuclides are simultaneously removed. If the moderation ratio (or the H/U ratio) varies owing to the submerging and the drilling operation, this variation also results in the simultaneous changes of $\rho$ and point kinetics parameters. In addition, because the detailed condition is not well known, there are large uncertainties in the preliminary analysis of $S$ and point kinetics parameters. Although the detailed material compositions and spatial distribution are unknown, the major fissile nuclide in fuel debris seems to be ²³⁵U. The relative delayed neutron yields $a_i={\beta }_i/{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f}}$ of ²³⁵U [5] is reasonable as the representative value used in the subcritical measurement.

In the above-mentioned simultaneous transient, time variation in neutron count rate cannot be distinguished into each of the components due to $-\rho$, $S$, and point kinetics parameters, meaning that the subcriticality measurement by the conventional neutron source multiplication (NSM) method [6–8] or inverse kinetics method [9–11] is difficult because $S$ and point kinetics parameters are generally assumed to be constant. Although the continuous subcriticality monitoring during the transient is a challenging problem, if steady-state values of neutron count rates can be measured under both steady-states before and after the transient, other measurement techniques have applicability to evaluate the subcriticality after the transient. For this investigation, we refocused on the integral method, which can estimate the subcriticality in dollar units after a stepwise transient, using only the time variation in the neutron count rate [12]. Note that the previous integral method had an issue with its applicability to a more general event where a stepwise state-change does not necessarily occur. Thus, in the previous study [13,14], ‘time-domain decomposition-based integral (TDDI) method’ was proposed using the point kinetics theory based on the fundamental mode approximation. Similar to the integral method, the TDDI method requires only the time variation in the neutron count rate to estimate the subcriticality. Furthermore, when an arbitrary state-change (e.g. step- or ramp-wise variations in $-\rho$, $S$, and kinetics parameters) occurs, the TDDI method can approximately estimate the subcriticality without the following information: (1) the reference subcriticality $-{\rho }_0$ before the transient, (2) the absolute values of neutron source intensity and the point kinetics parameters before and after the transient, and (3) the time variation in them. The above-mentioned key points are important advantages in the use of the TDDI method, especially for an unknown system. Note that, since the TDDI method is based on the point kinetics theory, there is a measurement error due to this approximation. To reduce the measurement error due to the higher-order modes of neutron flux, neutron detector position should be carefully chosen.

The goal of our research is to develop the subcriticality measurement technique after the transient, when the neutron count during a certain counting gate width (typical width between 0.1 and 1 s) is successively measured in the subcritical system. In the previous study, the verification of the integral and TDDI method was performed via virtual numerical experiments [12,13]. Subsequently, the applicability of the TDDI method to an actual experimental facility was investigated through transient experiments in a source-driven subcritical system conducted at the Kyoto University Critical Assembly (KUCA) [14]. However, a transient experiment with a simultaneous change of $-\rho$ and $S$ was not sufficiently investigated because a startup source of Am-Be was utilized as the steady external neutron source. This study aims to investigate the further applicability of the TDDI method to a more complicated transient caused by the simultaneous change. For this purpose, transient experiments with the operation of control rods and an external neutron source were conducted at the Kindai (formerly named Kinki) University Training and Research Reactor (UTR-KINKI).

Subsequent sections of the paper are structured as follows. In Section 2, the theory of the TDDI method is explained. The experimental conditions and results are described in Section 3. Finally, conclusions are presented in Section 4.

2. Time-domain decomposition-based integral method

2.1. Theory

In this section, the theory of the measurement method for the subcriticality in dollar units using the time-domain decomposition-based integral method (TDDI) is explained. The proposed methodology is based on the ‘point kinetics theory,’ which can be derived by the fundamental mode approximation for space and time-dependent neutron transport equation. By assuming time-dependent neutron flux $\varphi \left(\vec{r},E,t\right)$ can be expanded by forward $k_{\mathrm{e}\mathrm{f}\mathrm{f}}$-eigenfunctions and the fundamental mode component is dominant, $\varphi \left(\vec{r},E,t\right)$ is approximated by a product of amplitude function $f\left(t\right)$ and the fundamental eigenfunction ${\psi }_0\left(\vec{r},E\right)$:

(1) $\varphi \left(\overrightarrow{r},E,t\right)\approx f\left(t\right){\psi }_0\left(\overrightarrow{r},E\right),$(1)

where ${\psi }_0\left(\vec{r},E\right)$ satisfies the following $k_{\mathrm{e}\mathrm{f}\mathrm{f}}$-eigenvalue equation:

(2) $\mathbf{A}{\psi }_0=\frac{1}{k_{\mathrm{e}\mathrm{f}\mathrm{f}}}\mathbf{F}{\psi }_0,$(2)

where $\mathbf{A}$ and $\mathbf{F}$ represent net neutron loss and total neutron production operators, respectively. Based on this assumption, time-dependent neutron count rate $n\left(t\right)$ is approximately proportional to the amplitude function $f\left(t\right)$ related to the fundamental mode:

(3) $n\left(t\right)={⟨\mathrm{\Sigma }}_{\mathrm{d}}\varphi ⟩\approx \epsilon f\left(t\right),$(3)

(4) $\epsilon \equiv {⟨\mathrm{\Sigma }}_{\mathrm{d}}{\psi }_0⟩,$(4)

where ${\mathrm{\Sigma }}_{\mathrm{d}}\left(\vec{r},E\right)$ is macroscopic neutron detection cross-section; the brackets $⟨⟩$ indicate the integral over all phase space; and the factor $\epsilon$ corresponds to the detection efficiency, which is assumed to be constant in this study. The approximation of Equation (3)(3) $n\left(t\right)={⟨\mathrm{\Sigma }}_{\mathrm{d}}\varphi ⟩\approx \epsilon f\left(t\right),$(3) is appropriate for a shallow subcritical system. In the case of a deeper subcritical system, the neutron detector should be carefully located to reduce the higher-order mode components because the higher mode effect increases as the subcriticality deepens.

By applying the approximation of Equation (3)(3) $n\left(t\right)={⟨\mathrm{\Sigma }}_{\mathrm{d}}\varphi ⟩\approx \epsilon f\left(t\right),$(3) to the time-dependent neutron transport equation in a subcritical system with an external neutron source, the point kinetics equation with six delayed neutron precursor groups can be obtained as the following simultaneous differential equations [15]:

(5) $\frac{dn\left(t\right)}{dt}=\frac{\rho \left(t\right)-{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f}}\left(t\right)}{\mathrm{\Lambda }\left(t\right)}n\left(t\right)+\sum _{i=1}^6{\lambda }_i{C}_i\left(t\right)+S\left(t\right),$(5)

(6) $\frac{d{C}_i\left(t\right)}{dt}=\frac{{\beta }_i\left(t\right)}{\mathrm{\Lambda }\left(t\right)}n\left(t\right)-{\lambda }_i{C}_i\left(t\right),$(6)

(7) ${\beta }_{\mathrm{e}\mathrm{f}\mathrm{f}}\left(t\right)\equiv \sum _{i=1}^6{\beta }_i\left(t\right),$(7)

where the notations are conventional meanings in the reactor physics. Since $n\left(t\right)$ corresponds to the product $\epsilon$ and $f\left(t\right)$, the notations $C_i\left(t\right)$ and $S\left(t\right)$ are also defined by the products $\epsilon$ and these time-dependent fundamental modes of amplitude functions for the delayed neutron precursors and the external neutron source, respectively. It is noted that the point kinetics parameters (effective delayed neutron fraction ${\beta }_{\mathrm{e}\mathrm{f}\mathrm{f}}\left(t\right)$ and neutron generation time $\mathrm{\Lambda }\left(t\right)$) in Equations (5)(5) $\frac{dn\left(t\right)}{dt}=\frac{\rho \left(t\right)-{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f}}\left(t\right)}{\mathrm{\Lambda }\left(t\right)}n\left(t\right)+\sum _{i=1}^6{\lambda }_i{C}_i\left(t\right)+S\left(t\right),$(5) –(7(7) ${\beta }_{\mathrm{e}\mathrm{f}\mathrm{f}}\left(t\right)\equiv \sum _{i=1}^6{\beta }_i\left(t\right),$(7) ) are defined by the forward and adjoint fundamental $k_{\mathrm{e}\mathrm{f}\mathrm{f}}$-eigenfunctions, which depend on a state of core configuration at time $t$. The neutron count rate $n\left(t\right)$ and the amplitudes of delayed neutron precursors $C_i\left(t\right)$ vary over time, according to the temporal changes of not only the reactivity $\rho \left(t\right)$ but also the amplitude of external neutron source $S\left(t\right)$ and the point kinetics parameters.

Considering the steady-state before a transient (e.g. the ramp-wise state-change), the initial values of $n_0$ and $C_{i,0}$ before the transient can be obtained as follows:

(8) $n_0=\frac{S_0{\mathrm{\Lambda }}_0}{-{\rho }_0},$(8)

(9) $C_{i,0}=\frac{{\beta }_{i,0}}{{\lambda }_i{\mathrm{\Lambda }}_0}{n}_0,$(9)

where the subscript ‘0ʹ indicates the variable at steady-state before the transient. It is assumed that the transient of state-change occurs over a certain time interval from $t_0$ to $t_1$. An example of a ramp-wise transient of state-change is shown in Figure 1. It is noted that the transient of state-change is not necessarily linear ramp-wise as shown in Figure 1, although the variation in $n\left(t\right)$ is assumed to be monotonically increasing or decreasing.

After the transient of state-change terminates at $t=t_1$, $n\left(t\right)$ and $C_i\left(t\right)$ converge to the saturation values $n_{\mathrm{\infty }}$ and $C_{i,\mathrm{\infty }}$, respectively, in the limit of $t\to \mathrm{\infty }$. After the ending time $t_1$, the point kinetics equations are expressed as:

(10) $\frac{dn\left(t\right)}{dt}=\frac{{\rho }_1-{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}{{\mathrm{\Lambda }}_1}n\left(t\right)+\sum _{i=1}^6{\lambda }_i{C}_i\left(t\right)+S_1,$(10)

(11) $\frac{d{C}_i\left(t\right)}{dt}=\frac{{\beta }_{i,1}}{{\mathrm{\Lambda }}_1}n\left(t\right)-{\lambda }_i{C}_i\left(t\right),$(11)

where the subscript ‘1ʹ indicates variable after the transient ($t\ge t_1$). Based on the fundamental mode approximation, ${\mathrm{\Lambda }}_1$ and ${\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}$ are assumed to be constant because the forward and adjoint fundamental $k_{\mathrm{e}\mathrm{f}\mathrm{f}}$-eigenfunctions are unchanged after the transient. In Equations (10)(10) $\frac{dn\left(t\right)}{dt}=\frac{{\rho }_1-{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}{{\mathrm{\Lambda }}_1}n\left(t\right)+\sum _{i=1}^6{\lambda }_i{C}_i\left(t\right)+S_1,$(10) and (11(11) $\frac{d{C}_i\left(t\right)}{dt}=\frac{{\beta }_{i,1}}{{\mathrm{\Lambda }}_1}n\left(t\right)-{\lambda }_i{C}_i\left(t\right),$(11) ), the time derivatives become zero at steady-state after a sufficiently long time from $t_1$, if the system is still in the subcritical state after the transient. From Equations (10)(10) $\frac{dn\left(t\right)}{dt}=\frac{{\rho }_1-{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}{{\mathrm{\Lambda }}_1}n\left(t\right)+\sum _{i=1}^6{\lambda }_i{C}_i\left(t\right)+S_1,$(10) and (11(11) $\frac{d{C}_i\left(t\right)}{dt}=\frac{{\beta }_{i,1}}{{\mathrm{\Lambda }}_1}n\left(t\right)-{\lambda }_i{C}_i\left(t\right),$(11) ), the saturation values $n_{\mathrm{\infty }}$ and $C_{i,\mathrm{\infty }}$ satisfy the following equations:

(12) $0=\frac{{\rho }_1-{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}{{\mathrm{\Lambda }}_1}{n}_{\mathrm{\infty }}+\sum _{i=1}^6{\lambda }_i{C}_{i,\mathrm{\infty }}+S_1,$(12)

(13) $0=\frac{{\beta }_{i,1}}{{\mathrm{\Lambda }}_1}{n}_{\mathrm{\infty }}-{\lambda }_i{C}_{i,\mathrm{\infty }}.$(13)

By subtracting Equations (12)(12) $0=\frac{{\rho }_1-{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}{{\mathrm{\Lambda }}_1}{n}_{\mathrm{\infty }}+\sum _{i=1}^6{\lambda }_i{C}_{i,\mathrm{\infty }}+S_1,$(12) and (13(13) $0=\frac{{\beta }_{i,1}}{{\mathrm{\Lambda }}_1}{n}_{\mathrm{\infty }}-{\lambda }_i{C}_{i,\mathrm{\infty }}.$(13) ) from Equations (10)(10) $\frac{dn\left(t\right)}{dt}=\frac{{\rho }_1-{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}{{\mathrm{\Lambda }}_1}n\left(t\right)+\sum _{i=1}^6{\lambda }_i{C}_i\left(t\right)+S_1,$(10) and (11(11) $\frac{d{C}_i\left(t\right)}{dt}=\frac{{\beta }_{i,1}}{{\mathrm{\Lambda }}_1}n\left(t\right)-{\lambda }_i{C}_i\left(t\right),$(11) ) respectively, the following two equations can be obtained to eliminate the term of $S_1$:

(14) $\frac{dn\left(t\right)}{dt}=\frac{{\rho }_1-{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}{{\mathrm{\Lambda }}_1}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)+\sum _{i=1}^6{\lambda }_i\left(C_i\left(t\right)-C_{i,\mathrm{\infty }}\right),$(14)

(15) $\frac{d{C}_i\left(t\right)}{dt}=\frac{{\beta }_{i,1}}{{\mathrm{\Lambda }}_1}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)-{\lambda }_i\left(C_i\left(t\right)-C_{i,\mathrm{\infty }}\right).$(15)

By integrating Equations (14)(14) $\frac{dn\left(t\right)}{dt}=\frac{{\rho }_1-{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}{{\mathrm{\Lambda }}_1}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)+\sum _{i=1}^6{\lambda }_i\left(C_i\left(t\right)-C_{i,\mathrm{\infty }}\right),$(14) and (15(15) $\frac{d{C}_i\left(t\right)}{dt}=\frac{{\beta }_{i,1}}{{\mathrm{\Lambda }}_1}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)-{\lambda }_i\left(C_i\left(t\right)-C_{i,\mathrm{\infty }}\right).$(15) ) within the range of $t_1\le t<\mathrm{\infty }$, these equations are transformed into:

(16) $n_{\mathrm{\infty }}-n_1=\frac{{\rho }_1-{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}{{\mathrm{\Lambda }}_1}{\int }_{t_1}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt+\sum _{i=1}^6{\lambda }_i{\int }_{t_1}^{\mathrm{\infty }}\left(C_i\left(t\right)-C_{i,\mathrm{\infty }}\right)dt,$(16)

(17) $C_{i,\mathrm{\infty }}-C_{i,1}=\frac{{\beta }_{i,1}}{{\mathrm{\Lambda }}_1}{\int }_{t_1}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt-{\lambda }_i{\int }_{t_1}^{\mathrm{\infty }}\left(C_i\left(t\right)-C_{i,\mathrm{\infty }}\right)dt,$(17)

where $n_1$ and $C_{i,1}$ represent the variables just at $t=t_1$. By substituting Equation (17)(17) $C_{i,\mathrm{\infty }}-C_{i,1}=\frac{{\beta }_{i,1}}{{\mathrm{\Lambda }}_1}{\int }_{t_1}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt-{\lambda }_i{\int }_{t_1}^{\mathrm{\infty }}\left(C_i\left(t\right)-C_{i,\mathrm{\infty }}\right)dt,$(17) into Equation (16)(16) $n_{\mathrm{\infty }}-n_1=\frac{{\rho }_1-{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}{{\mathrm{\Lambda }}_1}{\int }_{t_1}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt+\sum _{i=1}^6{\lambda }_i{\int }_{t_1}^{\mathrm{\infty }}\left(C_i\left(t\right)-C_{i,\mathrm{\infty }}\right)dt,$(16) and utilizing Equation (7)(7) ${\beta }_{\mathrm{e}\mathrm{f}\mathrm{f}}\left(t\right)\equiv \sum _{i=1}^6{\beta }_i\left(t\right),$(7) , the following equation can be obtained:

(18) $n_{\mathrm{\infty }}-n_1=\frac{{\rho }_1}{{\mathrm{\Lambda }}_1}{\int }_{t_1}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt+\sum _{i=1}^6\left(C_{i,1}-C_{i,\mathrm{\infty }}\right).$(18)

The purpose of this study is to approximately estimate the subcriticality $-{\rho }_1/{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}$ using only the time variation in $n\left(t\right)$, without any information on the parameters after the state-change, i.e. it is required to eliminate $C_{i,1}$, $C_{i,\mathrm{\infty }}$ and ${\mathrm{\Lambda }}_1$ in Equation (18)(18) $n_{\mathrm{\infty }}-n_1=\frac{{\rho }_1}{{\mathrm{\Lambda }}_1}{\int }_{t_1}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt+\sum _{i=1}^6\left(C_{i,1}-C_{i,\mathrm{\infty }}\right).$(18) . For this purpose, firstly, the saturation value of $C_{i,\mathrm{\infty }}$ can be derived from Equation (13)(13) $0=\frac{{\beta }_{i,1}}{{\mathrm{\Lambda }}_1}{n}_{\mathrm{\infty }}-{\lambda }_i{C}_{i,\mathrm{\infty }}.$(13) :

(19) $C_{i,\mathrm{\infty }}=\frac{{\beta }_{i,1}}{{\lambda }_i{\mathrm{\Lambda }}_1}{n}_{\mathrm{\infty }}.$(19)

Secondly, to evaluate $C_{i,1}$ using only the time variation in $n\left(t\right)$, both sides of Equation (6)(6) $\frac{d{C}_i\left(t\right)}{dt}=\frac{{\beta }_i\left(t\right)}{\mathrm{\Lambda }\left(t\right)}n\left(t\right)-{\lambda }_i{C}_i\left(t\right),$(6) are multiplied by $\mathrm{\Lambda }\left(t\right)/{\beta }_i\left(t\right)$ and transformed into the following time differential equation:

(20) $\frac{d}{dt}\left(\frac{\mathrm{\Lambda }(t)C_i(t)}{{\beta }_i(t)}\right)=\left(\frac{d}{dt}\left(log\frac{\mathrm{\Lambda }(t)}{{\beta }_i(t)}\right)-{\lambda }_{\mathrm{i}}\right)\left(\frac{\mathrm{\Lambda }(t)C_i(t)}{{\beta }_i(t)}\right)+n(t).$(20)

For simplicity, the proposed method assumes that the decay of delayed neutron precursor is dominant in the time variation in $\mathrm{\Lambda }\left(t\right)C_i\left(t\right)/{\beta }_i\left(t\right)$, i.e., $\left|\frac{d}{dt}\left(log\frac{\mathrm{\Lambda }\left(t\right)}{{\beta }_i\left(t\right)}\right)\right|$ during the transient is smaller than ${\lambda }_i$. Although the assumption of $\left|\frac{d}{dt}\left(log\frac{\mathrm{\Lambda }\left(t\right)}{{\beta }_i\left(t\right)}\right)\right|\ll {\lambda }_i$ yields the approximation error, this assumption is reasonable if the transient interval ($t_1-t_0$) is relatively long (e.g., several tens of seconds) for the change ratio of $\left({\mathrm{\Lambda }}_1{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},0}\right)/\left({\mathrm{\Lambda }}_0{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}\right)$. Based on this approximation, Equation (20)(20) $\frac{d}{dt}\left(\frac{\mathrm{\Lambda }(t)C_i(t)}{{\beta }_i(t)}\right)=\left(\frac{d}{dt}\left(log\frac{\mathrm{\Lambda }(t)}{{\beta }_i(t)}\right)-{\lambda }_{\mathrm{i}}\right)\left(\frac{\mathrm{\Lambda }(t)C_i(t)}{{\beta }_i(t)}\right)+n(t).$(20) can be approximately solved to obtain the solution of $\mathrm{\Lambda }\left(t\right)C_i\left(t\right)/{\beta }_i\left(t\right)$ at $t=t_1$:

(21) $\frac{{\mathrm{\Lambda }}_1{C}_{i,1}}{{\beta }_{i,1}}\approx \frac{{\mathrm{\Lambda }}_0{C}_{i,0}}{{\beta }_{i,0}}{e}^{-{\lambda }_i\left(t_1-t_0\right)}+{\int }_{t_0}^{t_1}n\left(t\right)e^{-{\lambda }_i\left(t_1-t\right)}dt=\frac{n_0}{{\lambda }_i}{e}^{-{\lambda }_i\left(t_1-t_0\right)}+{\int }_{t_0}^{t_1}n\left(t\right)e^{-{\lambda }_i\left(t_1-t\right)}dt,$(21)

where the relationship in Equation (9)(9) $C_{i,0}=\frac{{\beta }_{i,0}}{{\lambda }_i{\mathrm{\Lambda }}_0}{n}_0,$(9) is utilized to express ${\mathrm{\Lambda }}_1{C}_{i,1}/{\beta }_{i,1}$ using $n_0$ and $n\left(t\right)$. Thirdly, by substituting Equations (19)(19) $C_{i,\mathrm{\infty }}=\frac{{\beta }_{i,1}}{{\lambda }_i{\mathrm{\Lambda }}_1}{n}_{\mathrm{\infty }}.$(19) and (21(21) $\frac{{\mathrm{\Lambda }}_1{C}_{i,1}}{{\beta }_{i,1}}\approx \frac{{\mathrm{\Lambda }}_0{C}_{i,0}}{{\beta }_{i,0}}{e}^{-{\lambda }_i\left(t_1-t_0\right)}+{\int }_{t_0}^{t_1}n\left(t\right)e^{-{\lambda }_i\left(t_1-t\right)}dt=\frac{n_0}{{\lambda }_i}{e}^{-{\lambda }_i\left(t_1-t_0\right)}+{\int }_{t_0}^{t_1}n\left(t\right)e^{-{\lambda }_i\left(t_1-t\right)}dt,$(21) ) into Equation (18)(18) $n_{\mathrm{\infty }}-n_1=\frac{{\rho }_1}{{\mathrm{\Lambda }}_1}{\int }_{t_1}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt+\sum _{i=1}^6\left(C_{i,1}-C_{i,\mathrm{\infty }}\right).$(18) , the subcriticality in dollar units $-{\rho }_1/{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}$ after the transient of state-change can be expressed as follows:

(22) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}=\frac{\frac{{\mathrm{\Lambda }}_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\left(n_1-n_{\mathrm{\infty }}\right)+{\sum }_{i=1}^6\frac{1}{{\lambda }_i}\frac{{\beta }_{i,1}}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\left(n_0{e}^{-{\lambda }_i\left(t_1-t_0\right)}-n_{\mathrm{\infty }}+{\lambda }_i{\int }_{t_0}^{t_1}n\left(t\right)e^{-{\lambda }_i\left(t_1-t\right)}dt\right)}{{\int }_{t_1}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt}.$(22)

In Equation (22)(22) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}=\frac{\frac{{\mathrm{\Lambda }}_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\left(n_1-n_{\mathrm{\infty }}\right)+{\sum }_{i=1}^6\frac{1}{{\lambda }_i}\frac{{\beta }_{i,1}}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\left(n_0{e}^{-{\lambda }_i\left(t_1-t_0\right)}-n_{\mathrm{\infty }}+{\lambda }_i{\int }_{t_0}^{t_1}n\left(t\right)e^{-{\lambda }_i\left(t_1-t\right)}dt\right)}{{\int }_{t_1}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt}.$(22) , the following approximation is assumed to be applicable. First, relative delayed neutron yields are known as constant, i.e. ${\beta }_{i,1}/{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}\approx a_i$, where $a_i$ is the $i$th group relative delayed neutron yield [5]. Second, $\left|\frac{{\mathrm{\Lambda }}_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\left(n_1-n_{\mathrm{\infty }}\right)\right|$ is negligibly small compared to the second term of the numerator, because the magnitude of ${\mathrm{\Lambda }}_1/{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}\approx 0.001$ s is much smaller than $\sum _{i=1}^6\text{break}{\beta }_{i,1}/\left({\lambda }_i{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}\right)\approx \sum _{i=1}^6{a}_i/{\lambda }_i\approx 10$ s and the following magnitude correlation holds true if the time variation in $n\left(t\right)$ is monotonically increasing ($n_0<n_1<n_{\mathrm{\infty }}$) or decreasing ($n_0>n_1>n_{\mathrm{\infty }}$):

(23) $\begin{array}{rl}& \left|n_{\mathrm{\infty }}-n_0{e}^{-{\lambda }_i\left(t_1-t_0\right)}-{\lambda }_i{\int }_{t_0}^{t_1}n\left(t\right)e^{-{\lambda }_i\left(t_1-t\right)}dt\right|\\ & >\left|n_{\mathrm{\infty }}-n_0{e}^{-{\lambda }_i\left(t_1-t_0\right)}-{\lambda }_i{\int }_{t_0}^{t_1}{n}_1{e}^{-{\lambda }_i\left(t_1-t\right)}dt\right|\\ & =\left|n_{\mathrm{\infty }}-n_1+(n_1-n_0)e^{-{\lambda }_i(t_1-t_0)}\right|>|n_{\mathrm{\infty }}-n_1|.\end{array}$(23)

Finally, the TDDI method can be obtained to estimate the subcriticality in dollar units for the transient as shown in Equation (24)(24) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\approx \frac{{\sum }_{i=1}^6\frac{a_i}{{\lambda }_i}\left(n_0{e}^{-{\lambda }_i(t_1-t_0)}-n_{\mathrm{\infty }}+{\lambda }_i{\int }_{t_0}^{t_1}n(t)e^{-{\lambda }_i(t_1-t_0)}dt\right)}{{\int }_{t_1}^{\mathrm{\infty }}\left(n(t)-n_{\mathrm{\infty }}\right)dt}.$(24) :

(24) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\approx \frac{{\sum }_{i=1}^6\frac{a_i}{{\lambda }_i}\left(n_0{e}^{-{\lambda }_i(t_1-t_0)}-n_{\mathrm{\infty }}+{\lambda }_i{\int }_{t_0}^{t_1}n(t)e^{-{\lambda }_i(t_1-t_0)}dt\right)}{{\int }_{t_1}^{\mathrm{\infty }}\left(n(t)-n_{\mathrm{\infty }}\right)dt}.$(24)

Equation (24)(24) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\approx \frac{{\sum }_{i=1}^6\frac{a_i}{{\lambda }_i}\left(n_0{e}^{-{\lambda }_i(t_1-t_0)}-n_{\mathrm{\infty }}+{\lambda }_i{\int }_{t_0}^{t_1}n(t)e^{-{\lambda }_i(t_1-t_0)}dt\right)}{{\int }_{t_1}^{\mathrm{\infty }}\left(n(t)-n_{\mathrm{\infty }}\right)dt}.$(24) is the measurement principle of the TDDI method that can estimate the subcriticality in dollar units using only the time variation in $n\left(t\right)$ for an arbitrary transient of state-change in a subcritical system. If a transient of state-change is as quick as a stepwise change, i.e. $t_0\approx t_1$, the TDDI method of Equation (24)(24) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\approx \frac{{\sum }_{i=1}^6\frac{a_i}{{\lambda }_i}\left(n_0{e}^{-{\lambda }_i(t_1-t_0)}-n_{\mathrm{\infty }}+{\lambda }_i{\int }_{t_0}^{t_1}n(t)e^{-{\lambda }_i(t_1-t_0)}dt\right)}{{\int }_{t_1}^{\mathrm{\infty }}\left(n(t)-n_{\mathrm{\infty }}\right)dt}.$(24) can be further simplified into the integral method [12], which is essentially the same as the extended integral count technique [16]:

(25) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\approx \frac{{\sum }_{i=1}^6\frac{a_i}{{\lambda }_i}\left(n_0-n_{\mathrm{\infty }}\right)}{{\int }_{t_0}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt}.$(25)

Compared with Equation (25)(25) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\approx \frac{{\sum }_{i=1}^6\frac{a_i}{{\lambda }_i}\left(n_0-n_{\mathrm{\infty }}\right)}{{\int }_{t_0}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt}.$(25) , the subcriticality estimation using the TDDI method requires an additional numerical calculation for a convolution of ${\int }_{t_0}^{t_1}n\left(t\right)e^{\text{break}-{\lambda }_i\left(t_1-t\right)}dt$.

2.2. Actual experimental procedures

To measure the subcriticality using the TDDI method based on Equation (24)(24) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\approx \frac{{\sum }_{i=1}^6\frac{a_i}{{\lambda }_i}\left(n_0{e}^{-{\lambda }_i(t_1-t_0)}-n_{\mathrm{\infty }}+{\lambda }_i{\int }_{t_0}^{t_1}n(t)e^{-{\lambda }_i(t_1-t_0)}dt\right)}{{\int }_{t_1}^{\mathrm{\infty }}\left(n(t)-n_{\mathrm{\infty }}\right)dt}.$(24) , the initial value $n_0$, saturation value $n_{\mathrm{\infty }}$, convolution ${\int }_{t_0}^{t_1}n\left(t\right)e^{-{\lambda }_i\left(t_1-t\right)}dt$, and integral ${\int }_{t_1}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt$ are necessary. In the case of an actual experiment, a neutron count rate at $t=t_j$ is obtained from a neutron count $N_j$ during a certain counting gate width $\mathrm{\Delta }t$:

(26) $n\left(t_j\right)\approx \frac{{\int }_{t_j-\mathrm{\Delta }t/2}^{t_j+\mathrm{\Delta }t/2}n\left(t\right)dt}{{\int }_{t_j-\mathrm{\Delta }t/2}^{t_j+\mathrm{\Delta }t/2}dt}=\frac{N_j}{\mathrm{\Delta }t}.$(26)

Since each of neutron count $N_j$ statistically fluctuates, the time variation in the neutron count rate has a statistical error. If instantaneous values of neutron count rates at certain times are simply chosen as $n_0$ and $n_{\mathrm{\infty }}$, the subcriticality obtained by the TDDI method will have a large statistical error owing to the statistical fluctuations of $n_0$ and $n_{\mathrm{\infty }}$.

Figure 1. Example of a ramp-wise transient.

Figure 2. Time intervals for initial value, convolution, integral, and saturation value of neutron counts.

Figure 3. Top view of experimental core.

To address this problem, the mean values of the neutron count rates within two intervals are used as $n_0$ and $n_{\mathrm{\infty }}$ as shown in Figure 2. Namely, the moving average within the range of $t_{0,a}\le t<t_0$ is used as the initial value $n\left(t_0\right)=n_0$ in the steady-state. Similarly, the moving average within the range of $t_{\mathrm{\infty },a}\le t<t_{\mathrm{\infty },b}$ is regarded as the saturation value $n\left(t_{\mathrm{\infty }}\right)=n_{\mathrm{\infty }}$. Note that the neutron count rates within the range of $t_{\mathrm{\infty },a}$ and $t_{\mathrm{\infty },b}$ gradually approaches the ideal saturation value. Therefore, the average time $\left(t_{\mathrm{\infty },a}+t_{\mathrm{\infty },b}\right)/2$ is set as $t_{\mathrm{\infty }}$ in this experimental analysis. The convolution and integral in Equation (24)(24) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\approx \frac{{\sum }_{i=1}^6\frac{a_i}{{\lambda }_i}\left(n_0{e}^{-{\lambda }_i(t_1-t_0)}-n_{\mathrm{\infty }}+{\lambda }_i{\int }_{t_0}^{t_1}n(t)e^{-{\lambda }_i(t_1-t_0)}dt\right)}{{\int }_{t_1}^{\mathrm{\infty }}\left(n(t)-n_{\mathrm{\infty }}\right)dt}.$(24) are numerically calculated using the trapezoidal rule.

2.3. Statistical error estimation using random sampling method

In an actual experiment, there are statistical fluctuations in measured time-series data of neutron counts $N_j$, which result in the statistical error of estimated subcriticality using the TDDI method. If the probability distribution of $N_j$ follows the Poisson distribution, the statistical error ${\sigma }_{N_j}$ can be estimated by the standard deviation of the Poisson distribution, because the variance is equal to the mean:

(27) ${\sigma }_{N_j}\approx \sqrt{N_j}.$(27)

Note that, in a neutron multiplication system, the variance of neutron count is generally larger than the mean, i.e. the second-order neutron correlation factor $Y\equiv \left(\mathrm{v}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e}/\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}-1\right)$ is larger than zero owing to the fission chain reactions [17]. The $Y$ value represents the relative difference between the mean and the variance from the Poisson distribution, thus the statistical error of Equation (27)(27) ${\sigma }_{N_j}\approx \sqrt{N_j}.$(27) can be corrected using the $Y$ value:

(28) ${\sigma }_{N_j}\approx \sqrt{\left(1+Y\right)N_j}.$(28)

Here, the magnitude of $Y$ is proportional to neutron detection efficiency. If the neutron detection efficiency is relatively low or $Y<1$, the variance of neutron count is nearly equal to the mean because the Poisson process of neutron-detection is dominant. Then, the statistical error of neutron count for each counting gate can be roughly estimated by Equation (27)(27) ${\sigma }_{N_j}\approx \sqrt{N_j}.$(27) based on approximation of the Poisson distribution.

In addition, the statistical error of the subcriticality can be approximately estimated using the propagation of uncertainty:

(29) ${\sigma }_f\approx \sqrt{\sum _{j=1}^J{\left(\frac{\mathrm{\partial }f}{\mathrm{\partial }{N}_j}{\sigma }_{N_j}\right)}^2},$(29)

where $J$ is the total number of neutron-count data and $f\left(N_1,N_2\cdots ,N_J\right)$ is a function to evaluate the subcriticality from input variables (i.e., time-series data of neutron counts $N_j$), which corresponds to Equation (24)(24) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\approx \frac{{\sum }_{i=1}^6\frac{a_i}{{\lambda }_i}\left(n_0{e}^{-{\lambda }_i(t_1-t_0)}-n_{\mathrm{\infty }}+{\lambda }_i{\int }_{t_0}^{t_1}n(t)e^{-{\lambda }_i(t_1-t_0)}dt\right)}{{\int }_{t_1}^{\mathrm{\infty }}\left(n(t)-n_{\mathrm{\infty }}\right)dt}.$(24) . However, even if the statistical error of neutron count ${\sigma }_{N_j}$ is given by Equation (27)(27) ${\sigma }_{N_j}\approx \sqrt{N_j}.$(27) , it is cumbersome to estimate the statistical error of the subcriticality using the propagation of uncertainty. This is because the estimation formula using Equation (29)(29) ${\sigma }_f\approx \sqrt{\sum _{j=1}^J{\left(\frac{\mathrm{\partial }f}{\mathrm{\partial }{N}_j}{\sigma }_{N_j}\right)}^2},$(29) with Equation (24)(24) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\approx \frac{{\sum }_{i=1}^6\frac{a_i}{{\lambda }_i}\left(n_0{e}^{-{\lambda }_i(t_1-t_0)}-n_{\mathrm{\infty }}+{\lambda }_i{\int }_{t_0}^{t_1}n(t)e^{-{\lambda }_i(t_1-t_0)}dt\right)}{{\int }_{t_1}^{\mathrm{\infty }}\left(n(t)-n_{\mathrm{\infty }}\right)dt}.$(24) is complicated due to the numerical convolution and integral.

To address the cumbersome procedure of error estimation, the uncertainty quantification based on the random sampling method is applied in this study [18–20]. In other words, the statistical error of the TDDI method is evaluated using the random sampling method, instead of Equation (29)(29) ${\sigma }_f\approx \sqrt{\sum _{j=1}^J{\left(\frac{\mathrm{\partial }f}{\mathrm{\partial }{N}_j}{\sigma }_{N_j}\right)}^2},$(29) based on the propagation of uncertainty. The detailed procedures are presented below:

1. Original time-series data of neutron counts $\vec{N}=\left\{N_1,N_2,\cdots ,N_J\right\}$ are provided by arranging $N_j$ in a time-series order.

2. A vector of the virtual time-series data ${\vec{N}}^{\ast b}=\left\{N_1^{\ast b},N_2^{\ast b},\cdots ,N_J^{\ast b}\right\}$ is stochastically generated using random numbers from the Poisson distribution, to add the statistical error to the original time-series data $\vec{N}$. The superscript ‘$\ast b$’ indicates the $b$th virtual sample. Here, each count $N_j^{\ast b}$ is independently sampled based on the Poisson distribution where the parameter of mean is $\lambda =N_j$. Note that the correlation of statistical error between different counting gates is neglected in this study. For each sample $N_j^{\ast b}$, expected values of mean and variance are equal to the original $N_j$.

3. Using equation (24)(24) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\approx \frac{{\sum }_{i=1}^6\frac{a_i}{{\lambda }_i}\left(n_0{e}^{-{\lambda }_i(t_1-t_0)}-n_{\mathrm{\infty }}+{\lambda }_i{\int }_{t_0}^{t_1}n(t)e^{-{\lambda }_i(t_1-t_0)}dt\right)}{{\int }_{t_1}^{\mathrm{\infty }}\left(n(t)-n_{\mathrm{\infty }}\right)dt}.$(24) with ${\vec{N}}^{\ast b}$, the corresponding output parameter $(-{\rho }_1/{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}{)}^{\ast b}$ is evaluated

4. To estimate the statistical error of subcriticality, steps 2–3 are repeated $B$ times. Consequently, a set of output parameters $(-{\rho }_1/{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}{)}^{\ast b}$ are obtained for $b=1,2,\cdots ,B$, where $B$ is the sample size of output parameters and $B=1000$ in this study.

5. Finally, the statistical error of the subcriticality by the TDDI method can be simply estimated by the 95% confidence interval (CI) of $(-{\rho }_1/{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}{)}^{\ast b}$. In particular, $B$ output parameters $(-{\rho }_1/{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}{)}^{\ast b}$ are sorted in ascending order. The lower and upper limits of the 95% CI are simply estimated from the $\left(0.025\times B\right)$th and $\left(0.975\times B\right)$th smallest values of sorted $(-{\rho }_1/{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}{)}^{\ast b}$, respectively.

3. Experimental conditions of subcritical transient experiments

To confirm the applicability of the TDDI method to an actual subcritical system with simultaneous reactivity and source changes, subcritical transient experiments were conducted in the UTR-KINKI [21].

The UTR-KINKI reactor is a light-water-moderated and graphite-reflected two-core reactor. Figure 3 illustrates the top view of the experimental core with positions of a neutron detector and an external source. The reactor has four control rods in symmetrical positions. Two of them are assigned as safety rods (#1 and #2) for reactor scram and manual shutdown. Others are a ‘shim rod’ for scram and coarse adjustment of reactivity, and a ‘regulating rod’ for auto operation and fine adjustment of reactivity. Positive or negative reactivity insertion was given by withdrawing or inserting the shim and regulating rods. To change the effective neutron source strength, an external neutron source of Pu-Be was jerked or dropped through the source hole in the experimental system. Although the detection efficiency of neutron count decreases, a BF₃ neutron detector was located as far from the neutron source as possible to reduce the detection of direct neutrons from the Pu-Be source. Thereby, we tried reducing to the higher mode effect due to the external neutron source. Using the detector with a list-mode data acquisition system (ANSeeN, HSDMCA), time-series data of neutron counts were successively measured during the counting gate width $\mathrm{\Delta }t=0.1$ s, of which the discretization error seemed to be smaller than a typical statistical error of the neutron count during $\mathrm{\Delta }t$ in the present experiments (e.g., the relative statistical error of the neutron count is larger than 3%, as presented later in Table 3).

Table 1. Time configurations for TDDI method.

Case	$t_{0,a}$ (s)	$t_0$ (s)	$t_1$ (s)	$t_{\mathrm{\infty },a}$ (s)	$t_{\mathrm{\infty }}$ (s)	$t_{\mathrm{\infty },b}$ (s)
1	0	300	332	532	632	732
2	0	300	332	532	632	732
3	0	300	332	532	632	732
4	0	300	332	532	632	732
5	0	300	305	505	605	705

Table 2. Measurement results of subcriticality.

Case	reference ($)	NSM method ($)	integral method ($)	TDDI method ($)
1	0.463	0.0095 [0.0094,0.0096]^†	0.395 [0.391,0.399]	0.447 [0.439,0.455]
2	0.638	83.9 [81.5,86.7]	0.566 [0.562,0.570]	0.660 [0.653,0.668]
3	0.463	37.3 [36.6,38.1]	0.512 [0.507,0.517]	0.494 [0.487,0.500]
4	0.638	0.0443 [0.0440,0.0447]	0.831 [0.789,0.879]	0.738 [0.676,0.813]
5	2.093	2.194 [2.186,2.202]	1.744 [1.651,1.852]	2.245 [2.072,2.450]

^†95% CI.

Table 3. Mean, variance, and $Y$ value for initial and saturation values of neutron counts.

	$t_{0,a}\le t<t_0$			$t_{\mathrm{\infty },a}\le t<t_{\mathrm{\infty },b}$
Case	mean	variance	$Y$	mean	variance	$Y$
1	22	36	0.59	1506	2594	0.72
2	369	672	0.82	2.0	4.0	0.95
3	250	403	0.62	4.3	8.2	0.92
4	28	54	0.89	296	468	0.58
5	602	941	0.56	175	212	0.21

Figure 4. Time variations in neutron counts after transients.

When time-series data of neutron counts were measured from $t=0$ s, the following various transient experiments were carried out by operations of the shim and regulating rods, and the external neutron source in the subcritical core:

• Case 1: Positive reactivity insertion and source increase. Initially, two safety rods were fully withdrawn, the shim and regulating rods were fully inserted (initial subcriticality 0.638$), and the neutron source was partially inserted. At $t=300$ s, the partially inserted neutron source was dropped (or fully inserted in a moment) and withdrawing of regulating rod was started. The regulating rod was fully withdrawn at $t=332$ s.

• Case 2: Negative reactivity insertion and source reduction. The initial condition of the safety, shim, and regulating rods were the same as the last condition of Case 1, i.e. two safety rods and regulating rod were fully withdrawn, the shim rod was fully inserted (initial subcriticality 0.463$), while the neutron source was fully inserted. At $t=300$ s, the fully inserted neutron source was quickly jerked for approximately 0.5 s or less and inserting of regulating rod was started. The regulating rod was fully inserted at $t=332$ s.

• Case 3: Positive reactivity insertion and source reduction. The initial condition of the safety, shim, and regulating rods were the same as the last condition of Case 2 (initial subcriticality 0.638$), while the neutron source was fully inserted. At $t=300$ s, the fully inserted neutron source was quickly jerked and withdrawing of regulating rod was started. The regulating rod was fully withdrawn at $t=332$ s.

• Case 4: Negative reactivity insertion and source increase. The initial condition of the safety, shim, and regulating rods were the same as the last condition of Case 3 (initial subcriticality 0.463$), while the neutron source was partially inserted. At $t=300$ s, the partially inserted neutron source was dropped and inserting of the regulating rod was started. The regulating rod was fully inserted at $t=332$ s.

• Case 5: Two successive events of negative reactivity insertion. The initial condition of the safety, shim, and regulating rods were the same as the last condition of Case 4 (initial subcriticality 0.638$), and the neutron source was completely inserted. Initially, the safety rod #1 was dropped at $t=300$ s. Subsequently, the safety rod #2 was dropped at $t=305$ s.

Figure 4 shows the time variation in neutron counts for Cases 1–5. For Case 5, two figures with different time ranges are presented to show the overall trend and the short-term variation after two successive rod drops. As can be seen from Figure 4, the neutron counts after the transient in Case 3 decrease in spite of the positive reactivity insertion by withdrawing the regulating rod, because the external source was jerked. Contrarily, the neutron counts after the transient in Case 4 gradually saturates owing to the external source drop, although the negative reactivity insertion was given by the regulating rod. In Case 5, the transient behaviour was more complicated due to the two successive safety rod drops.

Table 1 lists the time configurations of $t_{0,a}$, $t_0$, $t_1$, $t_{\mathrm{\infty },a}$, $t_{\mathrm{\infty }}$, and $t_{\mathrm{\infty },b}$ for the TDDI method. To confirm the improvement of the measurement accuracy by the TDDI method, the estimation using the integral method approximated by Equation (25)(25) $\frac{-{\rho }_1}{{\beta }_{\mathrm{e}\mathrm{f}\mathrm{f},1}}\approx \frac{{\sum }_{i=1}^6\frac{a_i}{{\lambda }_i}\left(n_0-n_{\mathrm{\infty }}\right)}{{\int }_{t_0}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt}.$(25) was also performed in almost the same conditions except for $t_1=t_0$. In the integral method, the decomposition of the time-domain for the convolution and integral is not necessary, thus only one integral-interval was set from $t_0$ to $t_{\mathrm{\infty }}$. Using these conditions, evaluation of the subcriticality in dollar units by both the TDDI and integral method was performed.

Before the above-mentioned transient experiments, the excess reactivity and control rod worth were preliminarily measured by the positive period method and the rod drop method, respectively. These measurement values were utilized to evaluate reference values of the subcriticality in dollar units for each of Cases 1–5. For supplemental information, the subcriticality values after the transients were estimated by the conventional NSM method [6–8] using the initial subcriticality and the count ratio of $n_0$ to $n_{\mathrm{\infty }}$.

4. Experimental results

Table 2 lists the subcritical measurement results using the NSM, integral, and TDDI methods with these statistical errors (95% CIs). As can be seen from Table 2, the conventional NSM method for Cases 1–4 had larger discrepancies from the reference values. Especially in Cases 3 and 4, the estimated subcriticality is qualitatively inappropriate. This is because the conventional NSM method neglects the neutron-source change. Since the neutron source was unchanged in Case 5, the estimated subcriticality using the NSM method well agree with the reference value after the transient.

Interestingly, note that the integral and TDDI methods enabled estimation of subcriticality even when the increasing or decreasing trend of neutron count rate did not correspond to that of reactivity. In addition, it was confirmed that the measurement results of the TDDI method were closer to reference values than those of the integral method, because the integral method assumes a transient as a stepwise change. Although the statistical errors of estimated subcriticality using the integral method were smaller than those of the TDDI method, the difference between the estimated subcriticality using the integral method and reference values were quite larger than these statistical errors of 95% CIs.

The statistical error of the TDDI method seemed to increase as the magnitude of the integral ${\int }_{t_1}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt$ decreases according to the subcriticality and source strength after the transient, because $1/\sqrt{{\int }_{t_1}^{\mathrm{\infty }}\left(n\left(t\right)-n_{\mathrm{\infty }}\right)dt}$ could be a major contribution to the relative statistical error based on the propagation of uncertainty. Compared with the difference between the measurement results by the TDDI method and reference values, the estimated statistical errors by the 95% CIs based on the Poisson distribution seemed to be slightly underestimated. To discuss this reason, Table 3 summarizes the mean, variance, and second-order neutron correlation factor $Y$ for initial and saturation values of neutron counts. Here, the values of the mean, variance, and $Y$ were estimated from the time-series data of neutron counts within the range of steady-states ($t_{0,a}\le t<t_0$ and $t_{\mathrm{\infty },a}\le t<t_{\mathrm{\infty },b}$), therefore the estimated $Y$ value could be regarded as the correction factor for the statistical error of neutron counts in Equation (28)(28) ${\sigma }_{N_j}\approx \sqrt{\left(1+Y\right)N_j}.$(28) . Although the $Y$ values were roughly less than unity, the $Y$ value tended to be larger than zero as the subcriticality became shallow. Namely, since the neutron-correlation due to the fission chain reactions in the subcritical systems was detectable, the probability distributions of neutron counts slightly differed from the Poisson distribution. If the neutron-correlation effect is considered, the statistical errors by the random sampling method seem to be underestimated by a factor of $1/\sqrt{1+Y}$ based on Equation (28)(28) ${\sigma }_{N_j}\approx \sqrt{\left(1+Y\right)N_j}.$(28) . This correction method for the statistical error by the random sampling method is future research because the issue is more difficult for time-series data of neutron counts under a transient from $t_0$ to $t_1$.

In the present transient experiments, the starting time $t_0$ and ending time $t_1$ of transient were known beforehand by the operation records. The saturation time $t_{\mathrm{\infty }}$ was determined after a sufficiently long time had passed from $t_1$. If the operation information is unknown, the automatic determination of $t_0$ and $t_1$ from the measured time-series data of neutron counts will be an additional issue for the present TDDI method. Furthermore, the subcriticality in the present experiments was relatively shallow; therefore, the point kinetics parameters were nearly equal to those at critical state and time variations in the point kinetics parameters seemed to be small. To investigate the applicability of the TDDI method to a situation with time variation of point kinetics parameters, for example, the actual measurement for a deeper subcritical system is important research. In the case of a deeper subcritical system, detected neutron counts tend to be generally small, which results in a large statistical error of the TDDI method. Furthermore, the spatial higher mode effect becomes larger as the subcriticality deepens. Thus, the reduction techniques of the statistical error and the higher mode effect are necessary to be developed.

5. Conclusion

For further investigation about the applicability of the TDDI method to actual subcritical measurement, the transient experiments with simultaneous changes of subcriticality and external neutron source or to two successive safety rods dropping were conducted at UTR-KINKI. To evaluate the statistical error of subcriticality using the TDDI method, the uncertainty quantification using the random sampling method was utilized based on the approximation of the Poisson distribution of neutron counts. Consequently, it was validated that the subcriticality in dollar units by the TDDI method better agreed with the reference values than the previous integral method. Thus, the superiority of the TDDI method was demonstrated for the non-stepwise transients in the present experiments.

For further improvements of the TDDI method, the following future researches are necessary: (1) the automatic determination for time configurations of $t_0$, $t_1$, and $t_{\mathrm{\infty }}$; (2) the reduction techniques of the statistical error and the higher mode effect for a deeper subcritical system; and (3) the correction method for underestimation of the statistical error by the random sampling method based on the Poisson distribution.

Acknowledgments

This work was carried out as a joint research program of the UTR-KINKI under the support of the Graduate School of Engineering, Osaka University. The authors are grateful to all the staff at UTR-KINKI for their assistance during the experiments.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for Scientific Research (C) [Grant Number 19K05328].