The time-dependent neutron transport equation and the precursor concentration balance equation need to be simultaneously solved for transient reactor analysis. 1 v E ∂ ψ r → , E , Ω → , t ∂ t = − L ψ r → , E , Ω → , t − T ψ r → , E , Ω → , t + S ψ r → , E , Ω → , t + ∑ d = 1 G d χ d E 4 π λ d C d r → , t + 1 k 0 χ p E 4 π 1 − β F ψ r → , E , Ω → , t , (1) ∂ C d r → , t ∂ t = 1 k 0 β d F ψ r → , E , Ω → , t − λ d C d r → , t , (2) where k 0 is the multiplication factor at the initial steady state, ψ r → , E , Ω → , t is the angular flux and notations L, T, S, and F represent the standard leakage, transport, scattering, and fission operators respectively. The number of precursor groups is represented as G _d and all the other notations are those of the convention (Jo et al., 2016). L ψ r → , E , Ω → , t = Ω → ⋅ ∇ ψ r → , E , Ω → , t , (3) T ψ r → , E , Ω → , t = σ t r → , E , t ψ r → , E , Ω → , t , (4) S ψ r → , E , Ω → , t = ∫ d Ω → ′ ∫ d E ′ σ s r → , Ω → ′ ⋅ Ω → , E ′ → E , t ψ r → , E ′ , Ω → ′ , t , (5) F ψ r → , E , Ω → , t = ∫ d Ω → ′ ∫ d E ′ ν σ f r → , E ′ , t ψ r → , E ′ , Ω → ′ , t . (6)

Implementation of implicit Euler method with a time-step of Δ t S , the aforementioned set of equations reduces into the following form: L ψ r → , E , Ω → , t s + T P C Q S ψ r → , E , Ω → , t s = S ψ r → , E , Ω → , t s + 1 k 0 χ p E 4 π 1 − β F ψ r → , E , Ω → , t s + ∑ d = 1 G d χ d E 4 π λ d C d r → , t s , (7) where T _PCQS is defined as T P C Q S ≡ σ t r → , E , t s + 1 v E Δ t s ψ r → , E , Ω → , t s . (8)

The following linear variation assumption for fission source and exponential transformation are then applied.

1) Linear Variation of Fission Source

Substituting Eqs. 9,10 into Eq. 7, one yields the following expression for the time-dependent transport equation. L + T ∼ P C Q S − S ψ r → , E , Ω → , t s = ψ r → , E , Ω → , t s − 1 e γ s Δ t v E Δ t s + ∑ d = 1 G d χ d E 4 π λ d C d r → , t s − 1 f 1 , d + ∑ d = 1 G d χ d E 4 π β d F ψ r → , E , Ω → , t s − 1 k 0 f 2 , d + ∑ d = 1 G d χ d E 4 π β d F ψ r → , E , Ω → , t s k 0 f 3 , d + 1 k 0 χ p E 4 π 1 − β F ψ r → , E , Ω → , t s , (11) T ∼ P C Q S ψ r → , E , Ω → , t s = σ t r → , E , t s + 1 v E Δ t s + γ s v E ψ r → , E , Ω → , t s , (12) where f 1 , k = e − λ k Δ t s , f 2 , k = 1 − e − λ k Δ t s − λ k Δ t s e − λ k Δ t s λ k Δ t s , f 3 , k = e − λ k Δ t 1 − e λ k Δ t s + λ k Δ t s e λ k Δ t s λ k Δ t s . (13)

Similarly, the balance equation for precursor concertation can be expressed as C d r → , t s = C d r → , t s − 1 f 1 , d + β d F ψ r → , E , Ω → , t s − 1 λ d k 0 f 2 , d + β d F ψ r → , E , Ω → , t s λ d k 0 f 3 , d . (14)

One could recognize that Eq. 11 is in the form of a transient fixed-source problem (TFSP) and its calculation result directly determines the precursor concentration according to Eq. 14. The MC method applied to Eq. 11 samples the fission source for the updated time-step ( t S ) based on the banked fission source and precursor concentration from the previous time-step ( t S − 1 ), where the increment of time ( Δ t S ) for solving TFSP will be referred to as a macro time-step throughout the manuscript. Note that contribution from t S dependent terms in Eq. 11 should be determined for every cycle whilst updating the fission source.

As aforementioned, it is the correction based on the point-kinetics (PK) equation that determines the final solution for the PCQS-MC method. For derivation of PK equation, the angular flux is factorized into the amplitude function n(t) and the shape function as: ψ r → , E , Ω → , t = n t φ r → , E , Ω → , t . (15)

Note that such a factorization is not an assumption, but rather demands an additional equation to render such a factorization to be unique. Mathematically, the following condition is imposed, ∫ d V ∫ d Ω → ∫ d E W r → , E , Ω → φ r → , E , Ω → , t v E = ∫ d V ∫ d Ω → ∫ d E W r → , E , Ω → ψ r → , E , Ω → , t 0 v E , (16) where W r → , E , Ω → denotes a weighting function. It is the consistent usage of weighting function, not its value, that guarantees reliable factorization, although exploitation of adjoint flux as a weighting function is recommended for stifling 1st order errors (Ott and Neuhold, 1985).

By multiplying the weighting function to the original time-dependent balance equations and integrating over the space, angle, and energy, the following point-kinetics (PK) time-dependent differential equations can be obtained: d n t d t = ρ t − β t Λ t n t + ∑ d = 1 G d λ d c d t d c d t d t = − λ d c d t + β d t Λ t n t (17) where n(t) is the amplitude function and c _d (t) is the weighted integrals of precursor concentration. Note that Equation 17 will be referred to as the point-kinetics equations (PKE) throughout this manuscript. All the other notations are defined as below: ρ t = 1 − k 0 k t , w i t h k t = W r → , E , Ω → , χ E 4 π F φ r → , E , Ω → , t W r → , E , Ω → , L + T − S φ r → , E , Ω → , t , (18) β t = ∑ d = 1 G d β d t , w i t h β d t = W r → , E , Ω → , χ d E 4 π β d F φ r → , E , Ω → , t W r → , E , Ω → , χ E 4 π F φ r → , E , Ω → , t , (19) Λ t = W r → , E , Ω → , 1 v E φ r → , E , Ω → , t W r → , E , Ω → , 1 k 0 χ E 4 π F φ r → , E , Ω → , t (20) c t = W r → , E , Ω → , χ d E 4 π C d r → , t W r → , E , Ω → , 1 v E φ r → , E , Ω → , t (21) α p t = ρ t − β t Λ t (22)

The PK parameters are evaluated (tallied) based on Eqs 18–21 for every cycle, thus insinuating stochasticity in nature. From the cycle averaged values; Eq. 17 is solved using the implicit Euler method with a finer time-step difference of ∆ t K , where a linear interpolation of α p t is postulated between the consecutive macro time-steps ( ∆ t S ) of interest.

The total spectrum χ E that appears in the above equations is defined to satisfy the following condition: χ E F φ r → , E , Ω → , t = χ p E 1 − β F φ r → , E , Ω → , t + ∑ d χ d E β d F φ r → , E , Ω → , t , (23)

which is necessary for the preservation of Eq. 1 (Ott and Neuhold, 1985). For the steady-state MC calculation, by sampling the quasi-static precursor concentration, the expression given in Eq. 23 can be met, which allows the evaluation of effective (adjoint-weighted) kinetic parameters according to Eqs 19, 20 (Leppänen et al., 2014). Considering the TFSP for PCQS-MC calculation, Eq. 23 can be met by assessing the PK parameters based on particles originating from t _s related source terms in Eq. 11. As aforementioned, the t _s-1 related source terms are pre-evaluated from the previous time-step whereas the t _s related source terms are re-evaluated for every cycle.

It should be noted that variation exists for the definitions of PK parameters concerning the usage of proper spectrum, which also have been successfully applied for solving the time-dependent transport problem under the PCQS framework (Shang et al., 2018). Nevertheless, it must be articulated that the definitions included in this manuscript is also consistently formulated, where its applicability for PCQS-MC calculation have been demonstrated several times (Jo et al., 2016; Kim and Kim, 2021; Kim and Kim, 2022b; Kim, 2022).

The updated amplitude function n t s is then multiplied to the original solution ψ p r e d i c t e d r → , E , Ω → , t s of TFSP for correction as below: ψ c o r r e c t e d r → , E , Ω → , t S = ψ p r e d i c t e d r → , E , Ω → , t S n t S Z t S , (24) where Z t s is the normalization factor introduced for the consistent factorization according to Eq. 16. Note that the corrected angular flux ψ c o r r e c t e d r → , E , Ω → , t s is employed whilst solving Eq. 11 for the next time-step. Figure 1 summarizes the overall flowchart of PCQS-MC calculation, and Figure 2 depicts the comparison between macro and micro time-steps. Z t s ≡ W r → , E , Ω → ,   ψ p r e d i c t o r r → , E , Ω → , t s v E W r → , E , Ω → ,   ψ r → , E , Ω → , t 0 v E . (25)

Unlike the conventional steady-state MC calculation, not only the cycle-wise source correlation but also the innate uncertainty of tallied PK parameters contributes to the stochastic uncertainty of the PCQS-MC calculation. For a quantitative assessment of each uncertainty source, a null-transient simulation using the critical GODIVA benchmark has been considered. Note that ‘null-transient’ refers to a transient calculation without invoking any perturbation, hence any deviation from the original steady-state calculation result can be perceived as an extent of uncertainty.

The cycle-wise correlation-induced uncertainty can be easily captured by employing the conventional MC uncertainty analysis scheme. As shown in Figure 1, during the active cycle, a predictor transport is performed, where the power level of the predicted solution can be scored for each cycle. From the deviation of cycle-wise predicted power, its uncertainty can be obtained and can be scaled accordingly to Eq. 24 using the correction factor. However, such an approach cannot accommodate the uncertainty stemming from deducing PK parameters.

For the null-transient test using the critical GODIVA benchmark, three different cases have been envisioned. Case 1 is a non-controlled case which represents a normal PCQS-MC calculation condition. For case 2, the dynamic reactivity (ρ) used in the PK equation calculation is set to zero, which signifies the exclusion of PK parameter-associated uncertainty. In the last case, variance due to cycle-wise source correlation is suppressed through a large number of histories per cycle. Table 1 enumerates a detailed calculation condition for each case with a unit vector being universally employed as a weighting function in this work, which alleviates any concerns related to the proper usage of spectrum whilst tallying the PK parameters, i.e., energy dependency of χ E F φ r → , E , Ω → , t is solely determined by the spectrum.

Figure 3 depicts the calculated time-wise power profiles for each case from the null-transient simulation, where 500 micro time-steps were considered for each macro time-step interval. The uncertainty (error bar) corresponds to the 1-sigma range, and has been obtained from the cycle-wise tallied predicted power scaled with a correction factor. One can notice that the estimated PCQS-MC error bar cannot cover the extent of stochastic power fluctuation in the time domain for case 1. In addition, the effect of source correlation for the magnitude of fluctuation was marginal which can be seen from case 3.

However, the fluctuation of power in time was noticeably suppressed, being comparable to the range of the error bar, for case 2. Note that dynamic reactivity was intentionally set to be zero for such a case, i.e., zero uncertainty, implying that PK parameter-related uncertainty is the major source for fluctuation of the PCQS-MC calculation, which must be taken into account for the proper estimation of uncertainty.

From the null-transient study, it is plain that the uncertainty of tallied PK parameters, especially the dynamic reactivity, has the most noticeable effect on the fluctuation of PCQS-MC power estimation. As discussed, the aforementioned conventional cycle-wise uncertainty assessment severely underestimates such an inherent uncertainty of the PCQS-MC framework. To circumvent such a limitation of the conventional approach, a modified simulation scheme referred to as the point-kinetics (PK) sampling method is proposed in this work.

In comparison with the conventional PCQS-MC framework, which solves the PK equation based on the cycle-averaged tallied information, i.e., PK parameters, the PK sampling approach additionally solves the PK equation based on each cycle-wise tallied information. Note that the extra computing burden introduced by solving cycle-wise PK equation is insignificant compared to the MC transport calculation. The cycle-wise PK equation solution is then applied for each cycle-wise predicted power resulting in an accumulation of cycle-wise corrected power, in which the PK parameter-driven stochasticity is embedded. The overall procedure of the PK sampling scheme is shown in Figure 4.

Nevertheless, caution is needed whilst deducing the uncertainty from the heap of cycle-wise predicted power since solving the PK equation introduces non-linearity. It is worthwhile to explicitly mention that we are interested in the uncertainty of cycle-averaged PK parameter-based corrected power, not the average of cycle-wise corrected power. In fact, due to the non-linearity of the PK equation, one cannot guarantee equivalence between the cycle-averaged PK information-based corrected power and the average of cycle-wise PK information-based corrected power. Mathematically, one can express such an issue as below: f P K E ρ c y c ≈ E f P K ρ c y c , (26) where f _PK denotes the operation of solving PK equation, E represents an average, and ρ c y c is the cycle-wise PK information. In addition, it is expected that deviation between f P K E ρ c y c and E f P K ρ c y c will be strengthened with a smaller number of histories per cycle being sampled.

To validate the above speculation, a numerical test which mimics the null-transient PCQS-MC calculation has been envisioned. The PK calculation based on typical PWR kinetic parameters listed in Table 2 has been performed whilst perturbing the dynamic reactivity based on Eq. 18, where the multiplication factor is sampled from a normal distribution having a mean value of unity with a prescribed uncertainty ( σ k ) of 50 pcm, 100 pcm, and 150 pcm. By enlarging the uncertainty, one can mimic the case of having a reduced number of histories per cycle. The macro and micro time-steps for PK calculation were set to be 0.1 [s] and 0.2 [ms] respectively, and a set of 100 independent samples (cycles) were repeatedly calculated 25 times, i.e., a total of 25 batches, to acquire the uncertainty of cycle-averaged PK parameter-based corrected power.

Figure 5 illustrates the scatter plots of cycle-wise PKE calculated null-transient powers where the presence of seemingly outliers prevails with an increase of uncertainty. For comparison, the sample mean and 3-sigma range based on sample variance are marked for each case, where an increase in the multiplication uncertainty renders the sampled null-transient power to have a skewed distribution. Such an observation highlights the effect of non-linearity stemming from solving the PKE. As expected, the deviation between the average value of the cycle-wise PK calculation and PK calculation based on the cycle-averaged multiplication factor escalates with an increase in the uncertainty ( σ k ) as shown in Table 3, where the former tends to overestimate the reference null-transient calculation result which is unity by definition. Note that uncertainty for cycle-wise PK calculation is obtained from 100 samples, whereas the uncertainty for cycle-averaged PK calculation is estimated based on 25 independent batch calculations.

The overestimation of the average value of the cycle-wise PK calculation, i.e., the right-skewed distribution of cycle-wise PK calculation result as shown in Figure 5, can be explained from Eq. 17. The point-kinetics equation is in the form of 1^st order differential equation, implying exponential growth term in its solution. Even with the same magnitude of uncertainty for k(t) in Eq. 18, the extent of resulting deviation in the calculated corrected power will vary according to the sign of ρ t . To visualize such disparity, the PKE calculated result with an increase in the magnitude of uncertainty ( σ k ) has been performed as shown in Figure 6, where positive means + σ k has been applied during calculation whereas negative can be understood as vice versa. Note that exactly the prescribed value of uncertainty was imposed rather than sampling from a normal distribution, i.e., k = 1.0 ± σ k .

Figure 6 exhibits that the extent of change in the PKE calculated power becomes relatively larger for having a positive perturbation (+ σ k ) as the magnitude of the imposed uncertainty increases. Hence, the average value of the cycle-wise PK calculation will tend to overestimate the reference value as discussed, i.e., the non-linearity issue. From such an understanding, it can be argued that the applicability of Eq. 26 can be only met by having enough number of histories being sampled, i.e., reduced uncertainty in the tallied PK parameter. For such a case, the uncertainty of the PCQS-MC power can be directly calculated from the PK sampling scheme. However, uncertainty based on sample variance of cycle-wise corrected power will be inappropriate to be applied as less number of histories are involved in the calculation because the set of sampled cycle-wise corrected power cannot be interpreted as a normal distribution. For the proposed PK sampling scheme to be generally reliable, an additional systematic measure that overcomes non-linearity issues is required that enables proper estimation of uncertainty.

When a noticeable non-linearity effect exists, the sampled cycle-wise corrected power deviates from a normal distribution as discussed in Section 3.2. Especially, the presence of outlier-like PK corrected power manifests as an overestimation in the sample variance for null-transient calculation. From such observations, the authors conjectured that through a proper screening process which excludes outlier-like sample data, the remaining heap of cycle-wise sample data could be perceived as following a normal distribution and a reasonable uncertainty estimation can be made. Obviously, a mathematical standard must be employed for determining normality, i.e., the degree of following a normal distribution, and screening of the sample data.

To verify the normality of a certain sample, the Shapiro-Wilk (SW) normality test (Shapiro and Wilk, 1965), which is a null hypothesis testing, has been exploited. The p-value calculated from the SW normality test is compared with a prescribed significance level α, and the sample is regarded as normal when the p-value is bigger than alpha. When the given sample does not suffice such a condition, the data residing out of μ s ±   z α / 2 σ s range are excluded, where μ s and σ s are the sample mean and standard deviation of the sample, and ± z α / 2 are the critical values corresponding to the confidence interval of (1-α). The heap of data after the screening process undergoes the SW normality test again, and the overall procedure is repeated until the Shapiro-Wilk normality test is satisfied or all the elements reside within the range of screening, i.e., x i ∈ μ s − z α / 2 σ s , μ s + z α / 2 σ s . A more detailed description of the screening process is enumerated below:

Step1) Check whether the given set of sample data satisfies the Shapiro-Wilk normality test

The screening process has been applied to a null-transient simulation result with an uncertainty of 150 pcm. As shown in Table 3, the direct calculation of sample standard deviation from the heap of data overestimates the real uncertainty evaluated from batch-wise calculation. However, the uncertainty after the screening process with a significance level of 0.05 and z α / 2 of 1.96 becomes about 0.01810, which accurately resembles the real uncertainty. Note that a significance level of 0.05 is the standard value often imposed for a null hypothesis and z 0.05 / 2 equals 1.96. Figure 7 illustrates the Q-Q (Quantile-to-Quantile) plots for sample data before and after the screening process alongside a scatter plot denoting screened-out data with a blue colour.

As aforementioned, the fluctuation of power in time for a null-transient simulation can be directly regarded as an extent of PCQS-MC uncertainty. For the verification of the proposed PK sampling scheme, the uncertainty based on PK sampling including the screening process is compared with the real variance acquired from batch-wise calculation. The time-step was set to be 0.1 s with 150,000 histories per cycle, having 100 inactive and 100 active cycles. Each consecutive macro time interval consists of 500 micro time-steps for PK calculation, and a total number of 40 independent batch runs were considered for the assessment of real uncertainty. Regarding the screening process, a significance level of 0.05 with a z α / 2 of 1.96 has been imposed. Figure 9 depicts both the real variance and PK sampling-based uncertainty for a time duration of 1.0 s where each error bar corresponds to the 1-sigma range.

One could observe that extent of real and PK sampling-based uncertainties are consistent throughout the simulation, where the latter always underestimates the real variance due to the presence of cycle-wise correlation. Such an attribute has been reported before during the preliminary assessment of the PK sampling strategy (Kim and Kim, 2022b; Kim, 2022; Oh et al., 2022). Nevertheless, the presented null-transient result plainly attests to the effectiveness of the proposed scheme for deducing proper uncertainty for PCQS-MC transient calculation.

It is the screening process that guarantees the set of samples, i.e., cycle-wise PK corrected power, to follow a normal distribution. Figure 10 illustrates the Q-Q plot and the raw data that have been obtained at a time-step of 0.1 s, where the screened-out data are marked with a blue colour. Identical magnitudes of uncertainties were obtained with different values of significance levels: α = 0.046 ( z α / 2 = 2.00) and α = 0.055 ( z α / 2 = 1.920), which are conventional values often applied for null-hypothesis testing.

As discussed, the proposed PK sampling scheme postulates near equivalence concerning Eq. 26 for deducing proper uncertainty estimation from a heap of PK-corrected power. Although inherent non-linearity forever exists, it is expected that sampled cycle-wise data will follow a normal distribution when enough number of particles are considered. Such speculation has been partially tested and confirmed using a numerical test that mimics a null-transient, but it is essential for it to be appraised based on a real PCQS-MC calculation. Hence, the same null-transient test for the C5G7-TD benchmark has been performed with an enlarged number of particles, e.g., one-million histories per cycle with 150 active and inactive cycles. It must be articulated that the number of cycles is irrelevant to the normality issue, which only determines the number of sampled cycle-wise corrected power to undergo the screening process. It is the number of histories per cycle that dictates the extent of attaining normality. Figure 11 depicts the Q-Q plot for the raw data, where the cycle-wise corrected power satisfies the Shapiro-Wilk test without any screening process, i.e., a p-value of 0.2132 was obtained. Note that it is mathematically valid to regard the presented data to be sampled from a normal distribution considering such a large p-value.

The C5G7-TD0 benchmark transient calculation is induced with a step reactivity insertion, where the control rods abruptly plunge into 10% of the core and reside for one second. Then, the rods are lifted to half of the originally inserted depth, i.e., 5% of the core, and remain for another second. Finally, the control rods are fully withdrawn out of the core. With respect to the number of banks involved, several test cases are considered in the benchmark, where the presented manuscript postulates all the banks to be subjected to movement. The analytic representation of the guide thimble cross-section in time, where the control rods are inserted, is written as Σ x t = Σ x G T ,   t < 0 and t ≥ 2 s Σ x t = Σ x G T + 0.1 Σ x C R − Σ x G T ,   0 < t ≤ 1 s Σ x t = Σ x G T + 0.05 Σ x C R − Σ x G T ,   1 s < t ≤ 2 s , (27) where Σ x G T and Σ x C R are macroscopic cross-sections for the empty guide thimble and control rod loaded guide thimble respectively. For PCQS-MC calculation using iMC, 150,000 histories per cycle with 150 inactive and 150 active cycles were considered. To exhibit the inherent uncertainty of PCQS-MC calculation, a null-transient for 0.5 s has been performed before the initiation of transient, i.e., control rod movements. The macro time-step was set to be 0.1 s and was divided into 500 equally spaced micro time-steps for PK calculation. Figure 12 shows the evolution of normalized power from the PCQS-MC calculation along with its 2-sigma uncertainty evaluated using the PK sampling scheme, where a comparison is made with the deterministic transport code PANDAS-MOC result (Tao and Xu, 2022). Note that the red coloured error bars correspond to the case for including the screening process whereas the blue coloured ones do not.

Although seemingly no significant difference exists for the two evaluated uncertainties for Figure 12, it is philosophically imperative to include the screening process. As aforementioned, the sampled cycle-wise PK corrected power is subjected to a non-linearity issue, where one cannot guarantee the normality of the sampled data. Whereupon, the calculated sample variance without the screening process not only tends to overestimate the true uncertainty but its appropriateness for representing the uncertainty of the heap of cycle-wise corrected power also becomes vague. Figure 13 depicts the Q-Q plot and the raw data obtained at a time step of 4.0 s, where the screened-out data are marked with a blue colour. Mathematically, it is the data after the screening process that can be regarded as normal, which renders the calculated sample variance to be applicable for representing the uncertainty.

As discussed in Section 2.2, only the neutrons born from t _s related source terms should be included whilst tallying the PK parameters for retaining consistency. However, during the actual PCQS-MC calculation implemented in the iMC, all the neutrons are tracked for tallying the PK parameters regardless of its origin. Nevertheless, such difference does not impair the realization of Eq. 23 since predominant contribution originates from the t _s related source terms as shown in Figure 14.

Unlike the t _s-1 related source terms, which are predetermined from the previous time step, the t _s related source terms are updated for every active cycle. Hence, it can be envisaged that the uncertainty stemming from the cycle wisely evaluated t _s related source terms determines the uncertainty of tallied PK parameters, which manifests into the overall uncertainty of PCQS-MC calculation. Such an interpretation could explain the noticeable decrease in the PCQS-MC uncertainty exhibited in Figure 12, where the time intervals for reduced uncertainty and the summation of weights concur with each other as shown in Figure 14.

Benchmarks TD1 and TD2 postulate continuous (constant) movement of control rods which results in a ramp reactivity insertion. The duration of insertion is set to be one second, and the fractional depth of the rods after 1 s of insertion are 1% and 10% for TD1 and TD2 respectively. The rods are then withdrawn with the same speed, being fully withdrawn after 2 s after the initiation of a transient. The analytic representation analogous to Section 4.2 is given below where w denotes the fractional depth of insertion. Σ x t = Σ x G T ,   t < 0 and t ≥ 2 s Σ x t = Σ x G T + w Σ x C R − Σ x G T t ,   0 < t ≤ 1 s Σ x t = Σ x G T + w Σ x C R − Σ x G T 2 − t ,   1 s < t ≤ 2 s . (28)

The same calculation conditions for C5G7-TD0 were imposed for TD1 and TD2 simulation and the iMC PCQS-MC result are depicted in Figs. 15 and 16, where only the uncertainty after the screening process is included. Note that PCQS-MC estimated power and MOC-based power reside within the range of uncertainty, i.e., error bar, whilst accommodating the inherent fluctuation of the PCQS-MC method.