On the transient models of the VITAS code: applications to the C5G7-TD pin-resolved benchmark problem

On the transient models of the VITAS code: applications to the C5G7-TD pin-resolved benchmark problem 1 first-author https://orcid.org/0000-0002-6623-5174 1 1 1 1 corresp † zhangtengfei@sjtu.edu.cn † zhangtengfei@sjtu.edu.cn zhangtengfei@sjtu.edu.cn https://orcid.org/0000-0002-6763-2534 School of Nuclear Science and Engineering, Shanghai Jiao Tong University, Shanghai 2002040, China School of Nuclear Science and Engineering, Shanghai Jiao Tong University, Shanghai 2002040, China This article describes the transient models of the neutronics code VITAS that are used for solving time-dependent, pin-resolved neutron transport equations. VITAS uses the stiffness confinement method (SCM) for temporal discretization to transform the transient equation into the corresponding transient eigenvalue problem (TEVP). To solve the pin-resolved TEVP, VITAS uses a heterogeneous variational nodal method (VNM). The spatial flux is approximated at each Cartesian node using finite elements in the x-y plane and orthogonal polynomials along the z -axis. Angular discretization utilizes the even-parity integral approach at the nodes and spherical harmonic expansions at the interfaces. To further lower the computational cost, a predictor-corrector quasi-static SCM (PCQ-SCM) was developed. Within the VNM framework, computational models for the adjoint neutron flux and kinetic parameters are presented. The direct-SCM and PCQ-SCM were implemented in VITAS and verified using the two-dimensional (2D) and three-dimensional (3D) exercises on the OECD/NEA C5G7-TD benchmark. In the 2D and 3D problems, the discrepancy between the direct-SCM solver’s results and those reported by MPACT and PANDAS-MOC was under 0.97% and 1.57%, respectively. In addition, numerical studies comparing the PCQ-SCM solver to the direct-SCM solver demonstrated that the PCQ-SCM enabled substantially larger time steps, thereby reducing the computational cost 100-fold, without compromising numerical accuracy. Stiffness confinement method Quasi-static method C5G7-TD benchmark Pin-resolved transient analysis Introduction 1 Neutronic calculations are critical for analyzing different types of nuclear systems[ 1 1 - 3 3 1-3 ]. Without spatial homogenization and lower-order approximations[ 4 4 ], high-fidelity neutron transport methods can accurately predict the steady-state and transient behaviors of neutron flux and thermal power during operation or severe accidents in nuclear reactors. Solving the transient pin-resolved equation is particularly challenging because the calculation is computationally expensive and it is difficult to achieve a balance between acceptable numerical accuracy and computational cost. The method of characteristics (MOC) and the finite element method (FEM) have shown promise in this area for high-fidelity calculations. In recent years, several high-fidelity neutronics codes have been developed. The MOC-based codes MPACT[ 4 4 ], PROTEUS-MOC[ 5 5 ], PANDAS-MOC[ 6 6 ], NECP-X[ 7 7 , 8 8 ], and the FEM-based code Rattlesnake[ 9 9 ] constitute some typical examples. In addition, the variational nodal method (VNM) has shown tremendous promise for solving high-fidelity problems[ 10 10 ]. The VNM uses the functional form of the second-order even-parity transport equation, with odd-parity Lagrange multipliers for enforcing nodal balance[ 11 11 , 12 12 ]. It has been demonstrated that the VNM is compatible with various discretization techniques and node geometries. Owing to its adaptability and accuracy, the VNM has been successfully implemented in neutronics codes, including VARIANT[ 13 13 ], VIOLET[ 14 14 ], and VITAS[ 15 15 ]. In particular, VITAS is a multi-purpose neutron transport solver with high-fidelity modeling capability (i.e., a heterogeneous VNM). Therefore, the spatial flux distribution is approximated using quadratic finite elements (FEs) and orthogonal polynomials. Within the node, the integral approach for angular discretization is implemented, whereas spherical harmonics functions are used at the interfaces. These discretization approaches enable VITAS to accomplish refinements in both space and angle. Efficient and accurate temporal discretization schemes are essential[ 16 16 ]. The quasi-static approach is the most frequently adopted temporal discretization scheme among high-fidelity methods. It consists of an improved quasi-static method and a predictor-corrector quasi-static (PCQ) approach[ 17 17 , 18 18 ]. In contrast to direct methods, the computational efficiency of the PCQ method is accomplished by coupling the transport equation and the exact point kinetics equation (EPKE) based on the observation that the flux shape varies more slowly with time than the flux amplitude. The shape function is computed by solving the transport equation with coarse time steps, whereas the amplitude function is determined by solving the EPKE equation with fine time steps. Consequently, high-fidelity transport codes that utilize the PCQ approach can provide time-dependent solutions at an acceptable computational cost. In the PCQ method, the transient problem is typically transformed into a transient fixed-source problem (TFSP) based on the backward difference method (BDM), which is then solved using TFSP solvers. However, in the absence of an efficient acceleration method, TFSP solvers converge slowly[ 5 5 ]. The stiffness confinement technique (SCM) is another method that can be used to address temporal dependence. In contrast to the BDM, the SCM transforms the transient problem into the corresponding transient eigenvalue problem (TEVP). The method was originally proposed to alleviate the stiffness in the point kinetics equation by incorporating dynamic frequencies[ 19 19 ], and subsequently extended for use in time-dependent diffusion and transport problems[ 20 20 - 22 22 20-22 ]. Recent studies on the use of the VNM in time-dependent problems have confirmed its promise. For instance, the VNM with the SCM has been examined in whole-core transport calculations with spatial homogenization[ 23 23 ]. It was shown that the VNM with the SCM leveraged coarse time steps, thus drastically minimizing the computation cost. To obtain an accurate flux amplitude, however, the traditional SCM (denoted as the direct-SCM) still requires finer time steps than the PCQ method. High-fidelity transient problems may impose large computing overheads. To extend the accurate and efficient application of the VNM and SCM to high-fidelity-transient situations, it is necessary to develop an improved VNM-SCM framework that incorporates the advantages of other temporal discretization techniques, such as the PCQ method. In this study, direct-SCM is used in conjunction with the VITAS code to solve time-dependent high-fidelity-transport problems. The PCQ-SCM is proposed and implemented in VITAS to reduce the computational cost, which incorporates the fundamental idea of the PCQ but obtains the flux shape by solving the TEVP, which is transformed by the SCM. The objective of this study was to demonstrate the capability of a heterogeneous VNM for solving transient pin-resolved problems. The remainder of the article is structured as follows. Methods for the direct-SCM, heterogeneous VNM, and PCQ-SCM are introduced and derived in Sect. 2. Section 3 describes the OECD/NEA C5G7-TD benchmark and the associated two-dimensional (2D) and three-dimensional (3D) exercises. The results obtained using the direct-SCM and PCQ-SCM approaches are compared in detail in Sect. 4. Moreover, comparisons with other codes are provided. Section 5 presents the summary and conclusions of the study. Theory 1 In the SCM, the time-dependent transport equation is first transformed into the corresponding TEVP using the frequency transformation. Then, by solving the TEVP using the steady solver iteratively, the necessary quantities for updating the frequencies are obtained. In VITAS, a steady solver based on a heterogeneous solver was used for solving the TEVP. In the following sections, we derive the frequency-transformed equation, which yields the form of the TEVP. Subsequently, based on the heterogeneous VNM, the nodal functional of the TEVP is presented and discretized to obtain the response matrix equations. Next, a flowchart of the essential iterations for updating frequencies in the direct-SCM approach is presented. Finally, the PCQ-SCM is formulated based on the direct-SCM and EPKE. Frequency-transformed equation 2 The dynamic frequency of the scalar flux is defined as <math display="block" id="M1"><mrow><msub><mi>ω</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>≡</mo><mfrac><mn>1</mn><mrow><msub><mi>ϕ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mfrac><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><msub><mi>ϕ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where the scalar flux <math id="M2"><mrow><msub><mi>ϕ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></math> of group g is given by <math display="block" id="M3"><mrow><msub><mi>ϕ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M4"><mrow><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mtext>r</mtext><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></math> is the angular flux at position 𝒓 in direction <math id="M5"><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle></math> . Note that to remove factors of π from these equations, dΩ is normalized such that <math id="M6"><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mtext>Ω</mtext><mo>=</mo><mn>1</mn></mrow></math> . Generally, it is assumed that the dynamic frequency of the angular flux is isotropic, such that <math display="block" id="M7"><mrow><msub><mi>ω</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mfrac><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>.</mo></mrow></math> Then, the time-variant angular flux can be expressed as <math display="block" id="M8"><mrow><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><msup><mi>e</mi><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mi>t</mi></msubsup></mrow></mstyle><mi>d</mi><msup><mi>t</mi><mo>'</mo></msup><msub><mi>ω</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mtext>r</mtext></mstyle><mn>,</mn><msup><mi>t</mi><mo>'</mo></msup></mrow><mo>)</mo></mrow></mrow></msup><mo>.</mo></mrow></math> The dynamic frequency can be further decomposed into <math id="M9"><mrow><msub><mi>ω</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><msub><mi>ω</mi><mrow><mi>S</mi><mn>,</mn><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><msub><mi>ω</mi><mtext>T</mtext></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> [ 21 21 ]. The flux amplitude frequency ω T ( t ) represents a global quantity and depends only on time, and the flux shape frequency <math id="M10"><mrow><msub><mi>ω</mi><mrow><mtext>S</mtext><mn>,</mn><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></math> depends on space, time, and energy. With this approximation, the angular flux can be decomposed using the amplitude and shape terms as <math display="block" id="M11"><mrow><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><mi>P</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mover accent="true"><mi>ψ</mi><mo>^</mo></mover><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where P ( t ) is the flux amplitude and <math id="M12"><mrow><msub><mover accent="true"><mi>ψ</mi><mo>^</mo></mover><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></math> is the flux shape. The time-variant flux amplitude and flux shape are evaluated as <math display="block" id="M13"><mrow><mi>P</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mi>P</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo><msup><mi>e</mi><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msup><mi>t</mi><mn>0</mn></msup></mrow><mi>t</mi></msubsup></mrow></mstyle><mtext>d</mtext><msup><mi>t</mi><mo>'</mo></msup><msub><mi>ω</mi><mtext>T</mtext></msub><mo stretchy="false">(</mo><msup><mi>t</mi><mo>'</mo></msup><mo stretchy="false">)</mo></mrow></msup><mo>,</mo></mrow></math> <math display="block" id="M14"><mrow><msub><mover accent="true"><mi>ψ</mi><mo>^</mo></mover><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><msub><mover accent="true"><mi>ψ</mi><mo>^</mo></mover><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><msup><mi>e</mi><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mi>t</mi></msubsup></mrow></mstyle><mi>d</mi><msup><mi>t</mi><mo>'</mo></msup><msub><mi>ω</mi><mrow><mi>S</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mtext>r</mtext></mstyle><mn>,</mn><msup><mi>t</mi><mo>'</mo></msup></mrow><mo>)</mo></mrow></mrow></msup><mo>.</mo></mrow></math> To make the decomposition unique[ 21 21 ], a physics-based constraint on the shape frequency is introduced, as follows: <math display="block" id="M15"><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub></mrow></mstyle><mtext>d</mtext><mi>r</mi></mrow></mstyle><mtext> </mtext><mi>ν</mi><msub><mi>Σ</mi><mrow><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mn>,</mn><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mi>d</mi><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><msub><mover accent="true"><mi>ψ</mi><mo>^</mo></mover><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><mi>P</mi><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where ν is the number of released neutrons per fission reaction and ∑ f g ( r ,t ) is the macroscopic fission cross-section. This constraint guarantees that the shape frequency affects only the flux shape and not the amplitude. Similarly, the precursor concentration frequency of delayed group i is defined as <math display="block" id="M16"><mrow><msub><mi>μ</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>≡</mo><mfrac><mn>1</mn><mrow><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mfrac><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where C i ( r ,t ) is the precursor concentration at position 𝒓 , which is commonly assumed to be isotropic. Introducing Eqs. 3 and 9 into the time-dependent transport equation [ 23 23 ] with isotropic scattering, the following frequency-transformed neutron transport equations are obtained: <math display="block" id="M17"><mtable columnalign="left"><mtr><mtd><mfrac><mrow><msub><mi>ω</mi><mtext>T</mtext></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>+</mo><msub><mi>ω</mi><mrow><mtext>S</mtext><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>v</mi><mi>g</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></mfrac><msub><mi>ψ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>⋅</mo><mo>∇</mo><msub><mi>ψ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mo>+</mo><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mtext>t</mtext><mi>g</mi></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle><msub><mi>ψ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mo>=</mo><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mtext>s</mtext><mi>g</mi></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle><msub><mi>ϕ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><msub><mi>q</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mtext>   </mtext><mi>g</mi><mo>=</mo><mn>1,</mn><mo>⋯</mo><mn>,</mn><mi>G</mi></mtd></mtr></mtable></math> where <math id="M18"><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mtext>t</mtext><mi>g</mi></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></math> and <math id="M19"><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mtext>s</mtext><mi>g</mi></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></math> are the macroscopic total and scattering cross sections, respectively. The group source q g ( r ,t ) is the contribution from scattering, prompt, and delayed neutrons: <math display="block" id="M20"><mtable><mtr><mtd><msub><mi>q</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><msup><mi>g</mi><mo>'</mo></msup><mo>≠</mo><mi>g</mi></mrow></munder><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mtext>s</mtext><mi>g</mi><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></mstyle></mtd></mtr><mtr><mtd><mo>+</mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>χ</mi><mi>g</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow><mrow><mo>[</mo><mrow><mn>1</mn><mo>−</mo><mi>β</mi><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow><mo>]</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>+</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></msubsup><mrow><mfrac><mrow><msub><mi>χ</mi><mrow><mi>i</mi><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>β</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>λ</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow><mrow><msub><mi>μ</mi><mi>i</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>λ</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></mfrac></mrow></mstyle></mrow></mtd></mtr></mtable></mrow><mo>}</mo></mrow></mtd></mtr><mtr><mtd><mo>×</mo><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><mi>ν</mi><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></mstyle></mtd></mtr></mtable></math> where <math id="M21"><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mtext>s</mtext><mi>g</mi><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mrow><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mstyle></mrow></math> is the macroscopic scattering cross-section, <math id="M22"><mrow><msub><mi>χ</mi><mi>g</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> and <math id="M23"><mrow><msub><mi>χ</mi><mrow><mi>i</mi><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> are the fission spectra of prompt and delayed neutrons, <math id="M24"><mrow><msub><mi>β</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> and <math id="M25"><mrow><mi>β</mi><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> are the fractions of delayed neutrons, and <math id="M26"><mrow><msub><mi>λ</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></math> is the delayed constant of the precursor. The frequency-transformed precursor equations are given as <math display="block" id="M27"><mtable columnalign="left"><mtr><mtd><msub><mi>μ</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><msub><mi>β</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mi>ν</mi></mstyle><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mtext>                         </mtext><mo>−</mo><msub><mi>λ</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mtext>   </mtext><mi>i</mi><mo>=</mo><mn>1,...,</mn><mi>I</mi></mtd></mtr></mtable></math> By introducing the dynamic eigenvalue k D into Eq. 10 and moving the frequency term to the right-hand side, the equation can be transformed into an eigenvalue problem (EVP) as <math display="block" id="M28"><mtable columnalign="left"><mtr><mtd><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>⋅</mo><mo>∇</mo><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>+</mo><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mtext>t</mtext><mi>g</mi></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle><msub><mi>ψ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mtext>    </mtext><mo>=</mo><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mtext>s</mtext><mi>g</mi></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle><msub><mi>ϕ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><msubsup><mi>q</mi><mi>g</mi><mo>'</mo></msubsup><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>,</mo></mtd></mtr></mtable></math> <math display="block" id="M29"><mtable columnalign="left"><mtr><mtd><msub><msup><mi>q</mi><mo>′</mo></msup><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><msup><mi>g</mi><mo>'</mo></msup><mo>≠</mo><mi>g</mi></mrow></munder><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mtext>s</mtext><mi>g</mi><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></mstyle><msub><mi>ϕ</mi><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mo>−</mo><mfrac><mrow><msub><mi>ω</mi><mtext>T</mtext></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>+</mo><msub><mi>ω</mi><mrow><mtext>S</mtext><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>v</mi><mi>g</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></mfrac><msub><mi>ϕ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mo>+</mo><msubsup><mi>k</mi><mtext>D</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msubsup><msub><msup><mi>χ</mi><mo>′</mo></msup><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mi>ν</mi></mstyle><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>,</mo></mtd></mtr></mtable></math> where <math id="M30"><mrow><msubsup><mi>χ</mi><mi>g</mi><mo>'</mo></msubsup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></math> is the dynamic fission spectrum, defined as <math display="block" id="M31"><mrow><msubsup><mi>χ</mi><mi>g</mi><mo>'</mo></msubsup><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>≡</mo><msub><mi>χ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo stretchy="false">)</mo><mrow><mo>[</mo><mrow><mn>1</mn><mo>−</mo><mi>β</mi><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover><mrow><mfrac><mrow><msub><mi>χ</mi><mrow><mi>i</mi><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>β</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>λ</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow><mrow><msub><mi>μ</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><msub><mi>λ</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></mfrac></mrow></mstyle><mo>.</mo></mrow></math> After the transformation into an EVP problem, Eq. 13 can be solved using existing neutron transport solvers. The dynamic frequencies corresponding to the dominant eigenvalue of 1 are the solutions to Eq. 13. Thus, the dynamic frequencies are solved iteratively using the power iteration and secant method for nonlinear equations. This nonlinear iteration algorithm for frequency-transformed equations is called the k-ω iteration [ 23 23 ], and is illustrated in Fig. 1 Fig. 1 . Variational nodal formulation 2 In this section, the discretized functional of the EVP is derived for the VNM. Finally, the corresponding response matrix equations are presented in a multi-group framework. The discretized functional 3 In the VNM, the even- and odd-parity fluxes ψ + and ψ - are defined, and the corresponding functional for node ν can be written with explicit separation of the radial and axial interfaces as [ 11 11 ] <math display="block" id="M32"><mtable><mtr><mtd><msub><mi>F</mi><mi>ν</mi></msub><mrow><mo>[</mo><mrow><msup><mi>ψ</mi><mo>+</mo></msup><mn>,</mn><msup><mi>ψ</mi><mo>−</mo></msup></mrow><mo>]</mo></mrow><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>ν</mi></msub></mrow></mstyle><mtext>d</mtext><mi>r</mi><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mtext>Ω</mtext><mrow><mo>[</mo><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mtext>t</mtext><mrow><mo>'</mo><mo>−</mo><mn>1</mn></mrow></msubsup><mrow><msup><mrow><mrow><mo>(</mo><mrow><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>⋅</mo><mo>∇</mo><msup><mi>ψ</mi><mo>+</mo></msup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></mstyle><mo>+</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><mtext>t</mtext><mo>'</mo></msubsup><mrow><msup><mi>ψ</mi><mrow><mo>+</mo><mn>2</mn></mrow></msup></mrow></mstyle></mrow><mo>]</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>−</mo><msub><mi>Σ</mi><mtext>s</mtext></msub><msup><mi>ϕ</mi><mn>2</mn></msup><mo>−</mo><mn>2</mn><mi>ϕ</mi><msup><mi>q</mi><mo>′</mo></msup></mrow></mtd></mtr></mtable></mrow><mo>}</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><mn>2</mn><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mi>z</mi><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>Γ</mi></msub></mrow></mstyle><mtext>d</mtext><mi>Γ</mi><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mtext>Ω</mtext><msub><mi>n</mi><mi>p</mi></msub><mo>⋅</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><msup><mi>ψ</mi><mo>+</mo></msup><msup><mi>ψ</mi><mo>−</mo></msup></mtd></mtr><mtr><mtd><mo>+</mo><mn>2</mn><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>A</mi></msub></mrow></mstyle><mtext>d</mtext><mi>A</mi><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mtext>Ω</mtext><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><msub><mstyle mathvariant="bold-italic"><mi>n</mi></mstyle><mrow><msup><mi>z</mi><mo>+</mo></msup></mrow></msub><mo>⋅</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><msup><mi>ψ</mi><mo>+</mo></msup><msup><mi>ψ</mi><mo>−</mo></msup><msub><mo>|</mo><mrow><msup><mi>z</mi><mo>+</mo></msup></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><mo>+</mo><msub><mstyle mathvariant="bold-italic"><mi>n</mi></mstyle><mrow><msup><mi>z</mi><mo>−</mo></msup></mrow></msub><mo>⋅</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><msup><mi>ψ</mi><mo>+</mo></msup><msup><mi>ψ</mi><mo>−</mo></msup><msub><mo>|</mo><mrow><msup><mi>z</mi><mo>−</mo></msup></mrow></msub></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mo>.</mo></mtd></mtr></mtable></math> For brevity, 𝒓 , Ω and t are suppressed in the unknowns. In the local system of coordinates, the node is defined as -Δ z /2≤ z ≤Δ z /2, the planar area is A =Δ x Δ y , 𝒏 p is the outward normal to the lateral interfaces extending over the periphery Γ, while <math id="M33"><mrow><msub><mstyle mathvariant="bold-italic"><mi>n</mi></mstyle><mrow><msup><mi>z</mi><mo>+</mo></msup></mrow></msub></mrow></math> and <math id="M34"><mrow><msub><mstyle mathvariant="bold-italic"><mi>n</mi></mstyle><mrow><msup><mi>z</mi><mo>−</mo></msup></mrow></msub></mrow></math> are the outward normals to the top and bottom axial interfaces, respectively. The functional in the entire problem domain is the superposition of the functionals in each nodal volume, given as <math id="M35"><mrow><mi>F</mi><mrow><mo>[</mo><mrow><msup><mi>ψ</mi><mo>+</mo></msup><mn>,</mn><msup><mi>ψ</mi><mo>−</mo></msup></mrow><mo>]</mo></mrow><mo>=</mo><mstyle displaystyle="true"><msub><mo>∑</mo><mi>ν</mi></msub><mrow><msub><mi>F</mi><mi>ν</mi></msub></mrow></mstyle><mrow><mo>[</mo><mrow><msup><mi>ψ</mi><mo>+</mo></msup><mn>,</mn><msup><mi>ψ</mi><mo>−</mo></msup></mrow><mo>]</mo></mrow></mrow></math> . Within the node, the spatial distribution of the even-parity flux is approximated by <math display="block" id="M36"><mrow><msup><mi>ψ</mi><mo>+</mo></msup><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>≈</mo><msubsup><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>z</mi><mtext>T</mtext></msubsup><mo stretchy="false">(</mo><mi>z</mi><mo stretchy="false">)</mo><mo>⊗</mo><msup><mstyle mathvariant="bold-italic"><mtext>g</mtext></mstyle><mtext>T</mtext></msup><mrow><mo>(</mo><mrow><mi>x</mi><mn>,</mn><mi>y</mi></mrow><mo>)</mo></mrow><mi>ψ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mn>.</mn></mrow></math> The scalar flux can be expanded as <math display="block" id="M37"><mrow><mi>ϕ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><msubsup><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>z</mi><mi>T</mi></msubsup><mo stretchy="false">(</mo><mi>z</mi><mo stretchy="false">)</mo><mo>⊗</mo><msup><mstyle mathvariant="bold-italic"><mtext>g</mtext></mstyle><mi>T</mi></msup><mrow><mo>(</mo><mrow><mi>x</mi><mn>,</mn><mi>y</mi></mrow><mo>)</mo></mrow><mi>ϕ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mn>.</mn></mrow></math> where the time-dependent scalar flux moments are defined as <math id="M38"><mrow><mi>ϕ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mi>ψ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></math> . <math id="M39"><mrow><msub><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>z</mi></msub><mrow><mo>(</mo><mi>z</mi><mo>)</mo></mrow></mrow></math> is a vector of orthogonal polynomials defined at node ν . The polynomials are governed by <math display="block" id="M40"><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mi>z</mi><msubsup><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>z</mi><mtext>T</mtext></msubsup><mrow><mo>(</mo><mi>z</mi><mo>)</mo></mrow><msub><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>z</mi></msub><mrow><mo>(</mo><mi>z</mi><mo>)</mo></mrow><mo>=</mo><mi>Δ</mi><mi>z</mi><mstyle mathvariant="bold-italic"><mi>I</mi></mstyle><mo>,</mo></mrow></math> where 𝑰 denotes the identity matrix. <math id="M41"><mrow><mi>g</mi><mrow><mo>(</mo><mrow><mi>x</mi><mn>,</mn><mi>y</mi></mrow><mo>)</mo></mrow></mrow></math> is a vector of continuous FE basis functions. Triangular and quadrilateral iso-parametric quadratic FEs [ 11 11 ] were employed to map the curved interfaces between materials. Figure 2 Figure 2 shows the FE mesh with 32 elements, used for describing a fuel pin cell in this paper. The odd-parity flux at the axial interface is approximated as <math display="block" id="M42"><mtable><mtr><mtd><msup><mi>ψ</mi><mo>−</mo></msup><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>≈</mo><msubsup><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>z</mi><mtext>T</mtext></msubsup><mrow><mo>(</mo><mrow><mo>±</mo><mi>Δ</mi><mi>z</mi><mo>/</mo><mn>2</mn></mrow><mo>)</mo></mrow><msubsup><mstyle mathvariant="bold-italic"><mtext>y</mtext></mstyle><mi>z</mi><mtext>T</mtext></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>)</mo></mrow><mo>⊗</mo><msup><mstyle mathvariant="bold-italic"><mi>h</mi></mstyle><mtext>T</mtext></msup><mrow><mo>(</mo><mrow><mi>x</mi><mn>,</mn><mstyle mathvariant="bold-italic"><mi>y</mi></mstyle></mrow><mo>)</mo></mrow><msub><mi>χ</mi><mi>z</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><msubsup><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>z</mi><mtext>T</mtext></msubsup><mrow><mo>(</mo><mrow><mo>±</mo><mi>Δ</mi><mi>z</mi><mo>/</mo><mn>2</mn></mrow><mo>)</mo></mrow><msubsup><mstyle mathvariant="bold-italic"><mtext>y</mtext></mstyle><mrow><mi>z</mi><mi>o</mi></mrow><mtext>T</mtext></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>)</mo></mrow><mo>⊗</mo><msup><mstyle mathvariant="bold-italic"><mi>h</mi></mstyle><mtext>T</mtext></msup><mrow><mo>(</mo><mrow><mi>x</mi><mn>,</mn><mstyle mathvariant="bold-italic"><mi>y</mi></mstyle></mrow><mo>)</mo></mrow><msub><mi>χ</mi><mrow><mi>z</mi><mi>o</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>,</mo></mtd></mtr></mtable></math> where <math id="M43"><mrow><mstyle mathvariant="bold-italic"><mi>h</mi></mstyle><mrow><mo>(</mo><mrow><mi>x</mi><mn>,</mn><mi>y</mi></mrow><mo>)</mo></mrow></mrow></math> denotes a vector of piecewise constants. <math id="M44"><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>y</mi></mstyle><mi>z</mi><mtext>T</mtext></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>)</mo></mrow></mrow></math> and <math id="M45"><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>y</mi></mstyle><mrow><mi>z</mi><mi>o</mi></mrow><mtext>T</mtext></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>)</mo></mrow></mrow></math> are odd-parity spherical harmonic vectors comprising both sine and cosine functions[ 24 24 ]. The vectors <math id="M46"><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>y</mi></mstyle><mi>z</mi><mtext>T</mtext></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>)</mo></mrow></mrow></math> are low-order P n approximations, whereas <math id="M47"><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>y</mi></mstyle><mrow><mi>z</mi><mi>o</mi></mrow><mtext>T</mtext></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>)</mo></mrow></mrow></math> and <math id="M48"><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>y</mi></mstyle><mrow><mi>γ</mi><mi>o</mi></mrow><mtext>T</mtext></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>)</mo></mrow></mrow></math> contain higher-order terms ranging from n +2 to a larger number N . In this study, we refer to these as P N P n approximations. <math id="M49"><mrow><msub><mi>χ</mi><mi>z</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> and <math id="M50"><mrow><msub><mi>χ</mi><mrow><mi>z</mi><mi>o</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> are the corresponding low- and high-order moments of the odd-parity flux on the axial interfaces, respectively. The odd-parity flux at the lateral interfaces is approximated as <math display="block" id="M51"><mtable><mtr><mtd><msup><mi>ψ</mi><mo>−</mo></msup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><msubsup><mstyle mathvariant="bold-italic"><mtext>y</mtext></mstyle><mi>γ</mi><mtext>T</mtext></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>)</mo></mrow><mo>⊗</mo><msubsup><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>γ</mi><mtext>T</mtext></msubsup><mrow><mo>(</mo><mi>ζ</mi><mo>)</mo></mrow><msub><mi>χ</mi><mi>γ</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><msubsup><mi>y</mi><mrow><mi>γ</mi><mi>o</mi></mrow><mtext>T</mtext></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>)</mo></mrow><mo>⊗</mo><msubsup><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mrow><mi>γ</mi><mi>o</mi></mrow><mtext>T</mtext></msubsup><mrow><mo>(</mo><mi>ζ</mi><mo>)</mo></mrow><msub><mi>χ</mi><mrow><mi>γ</mi><mi>o</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mi>γ</mi><mo>=</mo><mi>x</mi><mn>,</mn><mi>y</mi><mtext> </mtext><mi>ζ</mi><mo>=</mo><mi>y</mi><mn>,</mn><mi>x</mi></mtd></mtr></mtable></math> where γ =1,2,3,4 corresponds to the γ = x - , x + , y - , y + lateral interface. In VITAS, <math id="M52"><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>γ</mi><mtext>T</mtext></msubsup><mrow><mo>(</mo><mi>ζ</mi><mo>)</mo></mrow></mrow></math> and <math id="M53"><mrow><msubsup><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mrow><mi>γ</mi><mi>o</mi></mrow><mtext>T</mtext></msubsup><mrow><mo>(</mo><mi>ζ</mi><mo>)</mo></mrow></mrow></math> are chosen as fourth-order polynomials. Similarly, <math id="M54"><mrow><msubsup><mi>y</mi><mi>γ</mi><mtext>T</mtext></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>)</mo></mrow></mrow></math> and <math id="M55"><mrow><msubsup><mi>y</mi><mrow><mi>γ</mi><mi>o</mi></mrow><mtext>T</mtext></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>)</mo></mrow></mrow></math> are low- and high-order P n approximations, while <math id="M56"><mrow><msub><mi>χ</mi><mi>γ</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> and <math id="M57"><mrow><msub><mi>χ</mi><mrow><mi>γ</mi><mi>o</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> are the corresponding low- and high-order moments of the odd-parity flux on the lateral interfaces. Inserting the trial functions into Eqs. 17, 18, 20 and 21 yields the following discretized functional: <math display="block" id="M58"><mtable><mtr><mtd><mi>F</mi><mrow><mo>[</mo><mrow><mi>ψ</mi><mn>,</mn><msub><mi>χ</mi><mi>γ</mi></msub><mn>,</mn><msub><mi>χ</mi><mi>z</mi></msub></mrow><mo>]</mo></mrow><mo>=</mo><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><msup><mi>ψ</mi><mtext>T</mtext></msup><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mtext>A</mtext><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mi>ψ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mo>−</mo><msup><mi>ϕ</mi><mtext>T</mtext></msup><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mtext>I</mtext><mi>z</mi></msub><mo>⊗</mo><msub><mtext>F</mtext><mi>s</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mi>ϕ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><mn>2</mn><mi>ϕ</mi><msup><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mtext>T</mtext></msup><msup><mi>q</mi><mo>′</mo></msup><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><mn>2</mn><mstyle displaystyle="true"><munder><mo>∑</mo><mi>γ</mi></munder><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><msup><mi>ψ</mi><mtext>T</mtext></msup></mrow></mstyle><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><msub><mi>E</mi><mi>γ</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo stretchy="false">)</mo><msub><mi>χ</mi><mi>γ</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><mn>2</mn><mstyle displaystyle="true"><munder><mo>∑</mo><mi>γ</mi></munder><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><msup><mi>ψ</mi><mtext>T</mtext></msup></mrow></mstyle><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><msub><mi>E</mi><mi>z</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo stretchy="false">)</mo><msub><mi>χ</mi><mi>z</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><mn>2</mn><mstyle displaystyle="true"><munder><mo>∑</mo><mi>γ</mi></munder><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><msup><mi>ψ</mi><mtext>T</mtext></msup></mrow></mstyle><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><msub><mi>E</mi><mrow><mi>γ</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo stretchy="false">)</mo><msub><mi>χ</mi><mrow><mi>γ</mi><mi>o</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><mn>2</mn><mstyle displaystyle="true"><munder><mo>∑</mo><mi>γ</mi></munder><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><msup><mi>ψ</mi><mtext>T</mtext></msup></mrow></mstyle><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><msub><mi>E</mi><mrow><mi>z</mi><mi>o</mi></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo stretchy="false">)</mo><msub><mi>χ</mi><mrow><mi>z</mi><mi>o</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>,</mo></mtd></mtr></mtable></math> <math display="block" id="M59"><mrow><msup><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle><mo>′</mo></msup><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mi>r</mi><msub><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>z</mi></msub><mrow><mo>(</mo><mi>z</mi><mo>)</mo></mrow><mo>⊗</mo><mi>g</mi><mrow><mo>(</mo><mrow><mi>x</mi><mn>,</mn><mi>y</mi></mrow><mo>)</mo></mrow><msup><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle><mo>'</mo></msup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where 𝑨 and 𝑬 are the coefficient matrices. Matrices containing integrals over the spatial trial functions were numerically evaluated using standard FE techniques. The definitions of these matrices were detailed in a previous study[ 11 11 , 12 12 ]. Requiring the discretized functional in Eq. 22 to be stationary with respect to variations in <math id="M60"><mrow><mi>ψ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></math> yields: <math display="block" id="M61"><mtable columnalign="left"><mtr><mtd><mstyle mathvariant="bold-italic"><mi>A</mi></mstyle><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mi>ψ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mstyle mathvariant="bold-italic"><mi>I</mi></mstyle><mi>z</mi></msub><mo>⊗</mo><msub><mstyle mathvariant="bold-italic"><mi>F</mi></mstyle><mi>s</mi></msub><mi>ϕ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mtext>   </mtext><mo>=</mo><msup><mstyle mathvariant="bold-italic"><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle></mstyle><mo>'</mo></msup><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><msub><mstyle mathvariant="bold-italic"><mi>E</mi></mstyle><mi>l</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo stretchy="false">)</mo><msub><mi>χ</mi><mi>l</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><msub><mstyle mathvariant="bold-italic"><mi>E</mi></mstyle><mi>o</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo stretchy="false">)</mo><msub><mi>χ</mi><mi>o</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>.</mo></mtd></mtr></mtable></math> Note that the terms in Eq. 24 are re-grouped by low-order (with subscript l ) and high-order (with subscript o ) terms. Nodal response matrix equations 3 The odd-parity moments <math id="M62"><mrow><msub><mi>χ</mi><mi>l</mi></msub></mrow></math> and <math id="M63"><mrow><msub><mi>χ</mi><mi>o</mi></msub></mrow></math> are defined as continuous across interfaces. Also, requiring Eq. 22 to be stationary with respect to variations in <math id="M64"><mrow><msub><mi>χ</mi><mi>l</mi></msub></mrow></math> and <math id="M65"><mrow><msub><mi>χ</mi><mi>o</mi></msub></mrow></math> implies that the even-parity flux is continuous across interfaces: <math display="block" id="M66"><mrow><msub><mstyle mathvariant="bold-italic"><mi>φ</mi></mstyle><mi>l</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><msup><mstyle mathvariant="bold"><mtext>Λ</mtext></mstyle><mrow><mo>−</mo><mn>1</mn></mrow></msup><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><msubsup><mstyle mathvariant="bold-italic"><mtext>E</mtext></mstyle><mi>l</mi><mi>T</mi></msubsup><mrow><mo>(</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>)</mo></mrow><mi>ψ</mi><mrow><mo>(</mo><mrow><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>.</mo></mrow></math> By integrating the angle and matrix operations[ 11 11 ] into Eq. 24, the scalar flux can be determined in terms of the response matrix equations, given as <math display="block" id="M67"><mrow><mstyle mathvariant="bold-italic"><mi>ϕ</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mi>V</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msup><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle><mo>'</mo></msup><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><mstyle mathvariant="bold-italic"><mi>C</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mrow><mo>[</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>j</mi></mstyle><mo>+</mo></msup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>−</mo><msup><mstyle mathvariant="bold-italic"><mi>j</mi></mstyle><mo>−</mo></msup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow><mo>.</mo></mrow></math> Coefficient matrices can also be found in [ 11 11 , 12 12 ]. <math id="M68"><mrow><msup><mstyle mathvariant="bold-italic"><mi>j</mi></mstyle><mo>+</mo></msup></mrow></math> and <math id="M69"><mrow><msup><mstyle mathvariant="bold-italic"><mi>j</mi></mstyle><mo>−</mo></msup></mrow></math> correspond to the interfacial expansion moments of the outgoing and incoming neutron currents, respectively. They are defined as <math display="block" id="M70"><mrow><msup><mstyle mathvariant="bold-italic"><mi>j</mi></mstyle><mo>±</mo></msup><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mn>4</mn></mfrac><msub><mi>φ</mi><mi>l</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>±</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><msub><mi>χ</mi><mi>l</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mn>.</mn></mrow></math> The outgoing current moments <math id="M71"><mrow><msup><mstyle mathvariant="bold-italic"><mi>j</mi></mstyle><mo>+</mo></msup></mrow></math> satisfy <math display="block" id="M72"><mrow><msup><mstyle mathvariant="bold-italic"><mi>j</mi></mstyle><mo>+</mo></msup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle mathvariant="bold-italic"><mi>B</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msup><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle><mo>'</mo></msup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><mstyle mathvariant="bold-italic"><mi>R</mi></mstyle><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msup><mstyle mathvariant="bold-italic"><mi>j</mi></mstyle><mo>−</mo></msup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mn>.</mn></mrow></math> In summary, the flowchart of the VNM method is shown in Fig. 3 Fig. 3 . Direct-SCM in the VNM 2 The application of high-order iso-parametric FEs incurs significant computational cost and memory usage. To alleviate this burden, the element-wise flux values are stored in the code rather than the vertex-wise values generally adopted in the FEM framework. The element-wise average scalar flux is defined as <math display="block" id="M73"><mrow><mover accent="true"><mstyle mathvariant="bold-italic"><mi>ϕ</mi></mstyle><mo stretchy="true">¯</mo></mover><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><msup><mstyle mathvariant="bold"><mtext>Ξ</mtext></mstyle><mrow><mo>−</mo><mn>1</mn></mrow></msup><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mi>x</mi><mtext>d</mtext><mi>y</mi><mstyle mathvariant="bold-italic"><mi>h</mi></mstyle><mrow><mo>(</mo><mrow><mi>x</mi><mn>,</mn><mi>y</mi></mrow><mo>)</mo></mrow><msup><mstyle mathvariant="bold-italic"><mi>g</mi></mstyle><mi>T</mi></msup><mrow><mo>(</mo><mrow><mi>x</mi><mn>,</mn><mi>y</mi></mrow><mo>)</mo></mrow><msub><mstyle mathvariant="bold-italic"><mi>ϕ</mi></mstyle><mn>0</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where <math id="M74"><mrow><msub><mstyle mathvariant="bold-italic"><mi>ϕ</mi></mstyle><mn>0</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> is the segment of <math id="M75"><mrow><mstyle mathvariant="bold-italic"><mi>ϕ</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> corresponding to the components of the vector of orthogonal polynomials <math id="M76"><mrow><msub><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>z</mi></msub><mrow><mo>(</mo><mi>z</mi><mo>)</mo></mrow></mrow></math> of order 0, and <math id="M77"><mstyle mathvariant="bold"><mi>Ξ</mi></mstyle></math> is a diagonal matrix composed of the FE areas. <math id="M78"><mrow><mstyle mathvariant="bold-italic"><mi>h</mi></mstyle><mrow><mo>(</mo><mrow><mi>x</mi><mn>,</mn><mi>y</mi></mrow><mo>)</mo></mrow></mrow></math> is a vector of piecewise functions, which are equal to one in the domain of the corresponding FE and zero elsewhere. By averaging over each FE, the number of stored values is significantly reduced. Correspondingly, the flux shape frequency and precursor concentration are also defined as element-wise quantities. For brevity, the subscripts ν and k denote the average quantity of element k at node ν in Sect. 2.3. For example, ϕ ν k,g ( t ) is the element-wise average scalar flux for group g of element k at node ν . The calculation procedure for the direct-SCM within the framework of the VNM is shown in Fig. 4 Fig. 4 . In the k-ω iteration, the secant method [ 21 21 ] can be applied to update the amplitude frequency. The updated formula of the amplitude frequency based on the secant method is <math display="block" id="M79"><mtable><mtr><mtd><msubsup><mi>ω</mi><mtext>T</mtext><mrow><mo stretchy="false">(</mo><mi>m</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></msubsup><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><msubsup><mi>ω</mi><mtext>T</mtext><mrow><mo stretchy="false">(</mo><mi>m</mi><mo stretchy="false">)</mo></mrow></msubsup><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><mrow><mo>[</mo><mrow><msubsup><mi>ω</mi><mtext>T</mtext><mrow><mo stretchy="false">(</mo><mi>m</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></msubsup><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>−</mo><msubsup><mi>ω</mi><mtext>T</mtext><mrow><mo stretchy="false">(</mo><mi>m</mi><mo stretchy="false">)</mo></mrow></msubsup><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mfrac><mrow><mn>1</mn><mo>−</mo><msubsup><mi>k</mi><mtext>D</mtext><mrow><mo stretchy="false">(</mo><mi>m</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow><mrow><msubsup><mi>k</mi><mtext>D</mtext><mrow><mo stretchy="false">(</mo><mi>m</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></msubsup><mo>−</mo><msubsup><mi>k</mi><mtext>D</mtext><mrow><mo stretchy="false">(</mo><mi>m</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow></mfrac><mo>,</mo></mtd></mtr></mtable></math> where m denotes the iteration step of the k-ω iteration. The normalized element-wise average scalar flux <math id="M80"><mrow><msub><mover accent="true"><mi>ϕ</mi><mo>^</mo></mover><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow></math> is <math display="block" id="M81"><mtable><mtr><mtd><msub><mover accent="true"><mi>ϕ</mi><mo>^</mo></mover><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>×</mo><mfrac><mrow><mi>P</mi><mrow><mo>(</mo><mrow><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow></mrow><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mi>ν</mi></msub><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mi>k</mi></msub><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></msubsup><mi>ν</mi></mstyle></mrow></mstyle></mrow></mstyle><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>g</mi></mrow></msub><mrow><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow></mstyle><msub><mi>A</mi><mrow><mi>ν</mi><mi>k</mi></mrow></msub><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow></mfrac><mo>,</mo></mtd></mtr></mtable></math> in which A νk is the element area and <math id="M82"><mrow><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow></math> is the element-wise average scalar flux obtained from the EVP solver directly before the normalization. This normalization enforces the total power contributed by the normalized flux to equal the power of the previous time step. Therefore, the shape frequency obtained from the normalized flux is independent of the flux amplitude. In addition, it is assumed that the shape frequency is homogeneous within the element. Thus, based on the isotropic and homogeneous approximation, the element-wise shape frequency is <math display="block" id="M83"><mrow><msub><mi>ω</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>S</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><mi>Δ</mi><msub><mi>t</mi><mi>n</mi></msub></mrow></mfrac><mi>ln</mi><mrow><mo>[</mo><mrow><mfrac><mrow><msub><mover accent="true"><mi>ϕ</mi><mo>^</mo></mover><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>ϕ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow></mrow></mfrac></mrow><mo>]</mo></mrow><mo>,</mo></mrow></math> where Δ t n = t n - t n -1 , and ω ν k,S,g ( t n ) is the element-wise shape frequency within [ t , n -1 , t n ]. It is assumed that the amplitude frequency varies linearly within [ t n -1 , t n ]. <math display="block" id="M84"><mrow><msub><mi>ω</mi><mtext>T</mtext></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><msub><mi>ω</mi><mtext>T</mtext></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><mfrac><mrow><msub><mi>ω</mi><mtext>T</mtext></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo><mo>−</mo><msub><mi>ω</mi><mtext>T</mtext></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mi>t</mi><mo>−</mo><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow></mrow></math> Then, according to Eq. 5, the actual element-wise average scalar flux is <math display="block" id="M85"><mtable><mtr><mtd><msub><mi>ϕ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><msub><mi>ϕ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow><msup><mi>e</mi><mrow><mfrac><mrow><msub><mi>ω</mi><mtext>T</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>ω</mi><mtext>T</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow></mrow><mn>2</mn></mfrac><mi>Δ</mi><msub><mi>t</mi><mi>n</mi></msub><mo>+</mo><msub><mi>ω</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mtext>S</mtext><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mi>Δ</mi><msub><mi>t</mi><mi>n</mi></msub></mrow></msup></mtd></mtr><mtr><mtd><mo>=</mo><msub><mover accent="true"><mi>ϕ</mi><mo>^</mo></mover><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><msup><mi>e</mi><mrow><mfrac><mrow><msub><mi>ω</mi><mtext>T</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>ω</mi><mi>T</mi></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow></mrow><mn>2</mn></mfrac><mi>Δ</mi><msub><mi>t</mi><mi>n</mi></msub></mrow></msup><mo>.</mo></mtd></mtr></mtable></math> The first-order analytical precursor integration method is implemented to calculate the precursor concentration, where the fission source is assumed to vary linearly with time [ t n -1 , t n ]. The analytical solution for the element-wise precursor concentration is <math display="block" id="M86"><mtable><mtr><mtd><msub><mi>C</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><msub><mi>C</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow><msup><mi>e</mi><mrow><mo>−</mo><msub><mi>λ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><mi>Δ</mi><msub><mi>t</mi><mi>n</mi></msub></mrow></msup></mtd></mtr><mtr><mtd><mo>+</mo><msub><mi>β</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><msup><mi>e</mi><mrow><mo>−</mo><msub><mi>λ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><mi>Δ</mi><msub><mi>t</mi><mi>n</mi></msub></mrow></msup><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow></msubsup></mrow></mstyle><msub><mi>Q</mi><mrow><mi>ν</mi><mi>k</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msup><mi>e</mi><mrow><msub><mi>λ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><mi>t</mi></mrow></msup><mtext>d</mtext><mi>t</mi><mo>,</mo></mtd></mtr></mtable></math> in which the element-wise fission source <math id="M87"><mrow><msub><mi>Q</mi><mrow><mi>ν</mi><mi>k</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> is expressed as <math display="block" id="M88"><mrow><msub><mi>Q</mi><mrow><mi>ν</mi><mi>k</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><msub><mi>Q</mi><mrow><mi>ν</mi><mi>k</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><mfrac><mrow><msub><mi>Q</mi><mrow><mi>ν</mi><mi>k</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo><mo>−</mo><msub><mi>Q</mi><mrow><mi>ν</mi><mi>k</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow><mrow><mi>Δ</mi><msub><mi>t</mi><mi>n</mi></msub></mrow></mfrac><mrow><mo>(</mo><mrow><mi>t</mi><mo>−</mo><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M89"><mrow><msub><mi>Q</mi><mrow><mi>ν</mi><mi>k</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mi>g</mi></munder><mi>κ</mi></mstyle><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>g</mi></mrow></msub><mrow><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow></mstyle><msub><mi>ϕ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mn>,</mn></mrow></math> The updated formula for the element-wise precursor frequency is <math display="block" id="M90"><mrow><msub><mi>μ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>β</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><mfrac><mrow><msub><mi>Q</mi><mi>ν</mi></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo></mrow><mrow><msub><mi>C</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo></mrow></mfrac><mo>−</mo><msub><mi>λ</mi><mrow><msub><mrow/><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>C</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo><mo>≠</mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mrow><msub><mi>C</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo><mo>=</mo><mn>0</mn></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> PCQ-SCM in the VNM 2 To reduce the computational cost without a significant loss of accuracy, the PCQ-SCM is proposed. The main idea of the PCQ-SCM is to first solve the TEVP with coarse time steps, to obtain the predicted neutron flux. Then, by determining the kinetic parameters (KPs) with the predicted flux, the EPKE is solved with fine time steps, to obtain the amplitude function. The time step specifications are shown in Fig. 5 Fig. 5 . The error in the neutron flux resulting from the coarse time step is then corrected using the amplitude function. In the following sections, the PCQ-SCM is derived from the solution to the adjoint equation using the VNM. Then, the approach for evaluating the reactivity and kinetic parameters (KPs) with the VNM-based adjoint flux is discussed. Adjoint solution using the VNM 3 The solution to the adjoint steady transport equation, that is, the adjoint neutron flux, is adopted as a weighting function in perturbation problems. For example, the EPKE is derived from the first-order perturbation theory with the adjoint flux as the weighting function[ 25 25 , 26 26 ]. To implement the PCQ method in the transient calculation, the VNM-based adjoint flux should be used as the weighting function. The VNM transforms the original transport equation into a second-order even-parity equation. As the leakage operator in the second-order even-parity equation is self-adjoint[ 27 27 ], the adjoint calculation is significantly simplified. The forward and adjoint response matrix equations differ only in their source terms. Consequently, forward and adjoint calculations can be performed using the same within-group algorithms and response matrices shown in Fig. 3 Fig. 3 . The adjoint source term can be constructed by switching the group number of the scattering and fission cross-sections. In the following derivations, an asterisk denotes the solution of the adjoint equation. The adjoint response matrix equations for the initial condition are <math display="block" id="M91"><mrow><msup><mstyle mathvariant="bold-italic"><mi>ϕ</mi></mstyle><mo>*</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>=</mo><mstyle mathvariant="bold-italic"><mi>V</mi></mstyle><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><msup><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle><mo>*</mo></msup><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mo>−</mo><mstyle mathvariant="bold-italic"><mi>C</mi></mstyle><mo stretchy="false">(</mo><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mrow><mo>[</mo><mrow><msup><mstyle mathvariant="bold-italic"><mi>j</mi></mstyle><mrow><mo>*</mo><mo>+</mo></mrow></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>−</mo><msup><mstyle mathvariant="bold-italic"><mi>j</mi></mstyle><mrow><mo>*</mo><mo>−</mo></mrow></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M92"><mrow><msup><mstyle mathvariant="bold-italic"><mi>j</mi></mstyle><mrow><mo>*</mo><mo>+</mo></mrow></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>=</mo><mstyle mathvariant="bold-italic"><mi>B</mi></mstyle><mo stretchy="false">(</mo><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo><msup><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle><mo>*</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>+</mo><mi>R</mi><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><msup><mstyle mathvariant="bold-italic"><mi>j</mi></mstyle><mrow><mo>*</mo><mo>−</mo></mrow></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where the response matrices are identical to those in Eqns. 26 and 28. The vector of the initial adjoint source <math id="M93"><mrow><msup><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle><mo>*</mo></msup><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></math> is <math display="block" id="M94"><mrow><msup><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle><mo>*</mo></msup><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mo>=</mo><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mi>r</mi><msub><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>z</mi></msub><mrow><mo>(</mo><mi>z</mi><mo>)</mo></mrow><mo>⊗</mo><mstyle mathvariant="bold-italic"><mi>g</mi></mstyle><mrow><mo>(</mo><mrow><mi>x</mi><mn>,</mn><mi>y</mi></mrow><mo>)</mo></mrow><msup><mstyle mathvariant="bold-italic"><mi>q</mi></mstyle><mo>*</mo></msup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M95"><mtable><mtr><mtd><msubsup><mi>q</mi><mi>g</mi><mo>*</mo></msubsup><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><msup><mi>g</mi><mo>'</mo></msup><mo>≠</mo><mi>g</mi></mrow></munder><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mi>s</mi><msup><mi>g</mi><mo>′</mo></msup><mi>g</mi></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo></mrow></mstyle></mrow></mstyle><msubsup><mi>ϕ</mi><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mo>*</mo></msubsup><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mo>+</mo><msubsup><mi>k</mi><mrow><mtext>eff</mtext></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msubsup><mi>ν</mi><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><mi>g</mi></mrow></msub><mrow><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo></mrow></mstyle><mstyle displaystyle="true"><munder><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup></munder><mrow><msubsup><mi>χ</mi><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mo>*</mo></msubsup></mrow></mstyle><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><msubsup><mi>ϕ</mi><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mo>*</mo></msubsup><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>,</mo></mtd></mtr></mtable></math> <math display="block" id="M96"><mrow><msubsup><mi>χ</mi><mi>g</mi><mo>*</mo></msubsup><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>=</mo><msub><mi>χ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo stretchy="false">)</mo><mrow><mo>[</mo><mrow><mn>1</mn><mo>−</mo><mi>β</mi><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover></mstyle><msub><mi>χ</mi><mrow><mi>i</mi><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>β</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow><mo>.</mo></mrow></math> Note that eigenvalue k eff in the adjoint equation should be identical to that in the forward equation. EPKE 3 The EPKE is derived from the transient transport equation[ 28 28 ]. Similarly, the neutron flux is factorized into the shape function and amplitude function, as follows: <math display="block" id="M97"><mrow><msub><mi>ψ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><mi>T</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M98"><mrow><mi>T</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> is the amplitude function, and <math id="M99"><mrow><msub><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></math> is the shape function. Similarly, to make the decomposition unique, the following constraint is introduced: <math display="block" id="M100"><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub></mrow></mstyle></mrow></mstyle><mtext>d</mtext><mi>r</mi><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mtext>Ω</mtext><msubsup><mi>ψ</mi><mi>g</mi><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mfrac><mn>1</mn><mrow><msub><mi>v</mi><mi>g</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></mfrac><msub><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><mi>T</mi><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M101"><mrow><msubsup><mi>ψ</mi><mi>g</mi><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></math> is the adjoint angular flux in the steady state. It indicates that the integral of the shape function with the initial adjoint flux does not change with time. By inserting Eq. 44 into Eqns. 10 and 12, and multiplying by the adjoint flux on both sides, then integrating over all variables, the EPKE is obtained, as follows: <math display="block" id="M102"><mrow><mfrac><mtext>d</mtext><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mi>T</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mi>ρ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><mover accent="true"><mi>β</mi><mo stretchy="true">¯</mo></mover><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow><mrow><mi>Λ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mfrac><mi>T</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover><mrow><msub><mrow><mover accent="true"><mi>λ</mi><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub></mrow></mstyle><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mi>ζ</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M103"><mrow><mfrac><mtext>d</mtext><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><msub><mi>ζ</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><msub><mrow><mover accent="true"><mi>β</mi><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow><mrow><mi>Λ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mfrac><mi>T</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><msub><mrow><mover accent="true"><mi>λ</mi><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mi>ζ</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M104"><mrow><mi>ρ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> is the reactivity, <math id="M105"><mrow><mover accent="true"><mi>β</mi><mo stretchy="true">¯</mo></mover><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> and <math id="M106"><mrow><msub><mrow><mover accent="true"><mi>β</mi><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> are effective fractions of delayed neutrons, <math id="M107"><mrow><msub><mrow><mover accent="true"><mi>λ</mi><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> is the effective decay constant, <math id="M108"><mrow><mi>Λ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> is the prompt neutron lifetime, and <math id="M109"><mrow><msub><mi>ζ</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> is the reduced precursor concentration. The detailed formulation of the EPKE is available in [ 17 17 ]. Eqns. 46 and 47 can be solved using ordinary differential equation solvers. In VITAS, ordinary differential equations are solved using the SCM for EPKE [ 19 19 ]. As shown in Sect. 2.2, the even-parity integral method implemented in VITAS leads only to a solution of the scalar flux. Reconstructing the angular flux requires extra storage of temporary matrices for discretized operators[ 11 11 ] and extra algebraic operations of these matrices to the scalar flux. In addition, storing the angular flux significantly increases memory usage. Therefore, the scalar flux is used as the weighting function to calculate the KPs in VITAS. Note that reactivity is not included in the KPs in this section. In detail, the KPs calculated using the scalar flux <math id="M110"><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover></math> and ϕ † are <math display="block" id="M111"><mtable><mtr><mtd><mi>F</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mi>r</mi></mrow></mstyle><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mfrac><mrow><msub><mi>χ</mi><mi>g</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><msubsup><mi>ϕ</mi><mi>g</mi><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>×</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mi>G</mi></msubsup><mi>ν</mi></mstyle><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><msup><mi>g</mi><mo>'</mo></msup></mrow></msub><mrow><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mstyle><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow><mo>}</mo></mrow></mtd></mtr><mtr><mtd><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mi>ν</mi></munder><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mi>k</mi></munder><mrow><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mfrac><mrow><msub><mi>χ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><msub><mi>A</mi><mrow><mi>ν</mi><mi>k</mi></mrow></msub><msubsup><mi>ϕ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>×</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><msup><mi>g</mi><mo>′</mo></msup><mi>G</mi></msubsup><mi>ν</mi></mstyle><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><msup><mi>g</mi><mo>'</mo></msup></mrow></msub><mrow><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mstyle><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><msup><mi>g</mi><mo>'</mo></msup></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow><mo>}</mo></mrow></mrow></mstyle></mrow></mstyle></mrow></mstyle><mo>,</mo></mtd></mtr></mtable></math> <math display="block" id="M112"><mtable><mtr><mtd><msub><mover accent="true"><mi>β</mi><mo stretchy="true">¯</mo></mover><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mi>r</mi></mrow></mstyle><mrow><mo>{</mo> <mrow><mtable><mtr><mtd><mrow><mfrac><mrow><msub><mi>χ</mi><mrow><mi>i</mi><mi>g</mi></mrow></msub><mrow><mo>(</mo><mtext>r</mtext><mo>)</mo></mrow></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><msubsup><mi>ϕ</mi><mi>g</mi><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><msub><mi>β</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>×</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mi>G</mi></msubsup><mi>ν</mi></mstyle><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><msup><mi>g</mi><mo>'</mo></msup></mrow></msub><mrow><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mstyle><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow> <mo>}</mo></mrow></mrow><mrow><mi>F</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mfrac></mtd></mtr><mtr><mtd><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></msubsup><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mi>ν</mi></msub><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mi>k</mi></msub><mrow><mrow><mo>{</mo> <mrow><mtable><mtr><mtd><mrow><mfrac><mrow><msub><mi>χ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi><mi>g</mi></mrow></msub></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><msub><mi>A</mi><mrow><mi>ν</mi><mi>k</mi></mrow></msub><msubsup><mi>ϕ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><msub><mi>β</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><mo>×</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mi>G</mi></msubsup><mi>ν</mi></mstyle><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mstyle mathvariant="bold-italic"><mi>f</mi></mstyle><msup><mi>g</mi><mo>'</mo></msup></mrow></msub><mrow><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mstyle><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><msup><mi>g</mi><mo>'</mo></msup></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow> <mo>}</mo></mrow></mrow></mstyle></mrow></mstyle></mrow></mstyle></mrow><mrow><mi>F</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mfrac><mo>,</mo></mtd></mtr></mtable></math> <math display="block" id="M113"><mrow><mi>Λ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></msubsup><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mi>r</mi></mrow></mstyle><mfrac><mrow><msubsup><mi>ϕ</mi><mi>g</mi><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>v</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mfrac></mrow><mrow><mi>F</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mfrac><mo>=</mo><mfrac><mrow><mi>T</mi><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>F</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mfrac><mo>.</mo></mrow></math> <math display="block" id="M114"><mtable><mtr><mtd><msub><mover accent="true"><mi>λ</mi><mo stretchy="true">¯</mo></mover><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></msubsup><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mi>r</mi></mrow></mstyle><msubsup><mi>ϕ</mi><mi>g</mi><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><msub><mi>χ</mi><mrow><mi>i</mi><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>λ</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>C</mi><mi>i</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></msubsup><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mi>r</mi></mrow></mstyle><msubsup><mi>ϕ</mi><mi>g</mi><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><msub><mi>χ</mi><mrow><mi>i</mi><mi>g</mi></mrow></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow><msub><mi>C</mi><mi>i</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mfrac></mtd></mtr><mtr><mtd><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></msubsup><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mi>ν</mi></msub><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mi>k</mi></msub><mrow><msub><mi>A</mi><mrow><mi>ν</mi><mi>k</mi></mrow></msub><msubsup><mi>ϕ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow><mo>*</mo></msubsup></mrow></mstyle></mrow></mstyle></mrow></mstyle><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><msub><mi>χ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi><mi>g</mi></mrow></msub><msub><mi>λ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><msub><mi>C</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></msubsup><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mi>ν</mi></msub><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mi>k</mi></msub><mrow><msub><mi>A</mi><mrow><mi>ν</mi><mi>k</mi></mrow></msub><msubsup><mi>ϕ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow><mo>*</mo></msubsup></mrow></mstyle></mrow></mstyle></mrow></mstyle><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><msub><mi>χ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi><mi>g</mi></mrow></msub><msub><mi>C</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>i</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mfrac><mo>.</mo></mtd></mtr></mtable></math> Previous work showed that the KPs calculated using the scalar flux agreed well with those calculated using the angular flux [ 5 5 ]. However, using the scalar flux as the weighting function could result in significant reactivity deviations [ 5 5 , 29 29 ]. For example, in the moderation density reduction problems of the C5G7-TD benchmark, the scalar flux weighting function underestimated the reactivity by up to 15.6%[ 5 5 ]. The different reactivities resulting from the scalar flux were mainly owing to the strong anisotropy of the leakage operator <math id="M115"><mrow><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>⋅</mo><mo>∇</mo></mrow></math> . Therefore, adopting the isotropic approximation for evaluating the perturbation of the leakage operator could cause significant reactivity deviations. To obtain accurate reactivity values for the EPKE calculation, the reactivity in the present work is evaluated [ 30 30 ] as <math display="block" id="M116"><mrow><mi>ρ</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><mn>1</mn></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mfrac><mo>,</mo></mrow></math> where <math id="M117"><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> is the dynamic eigenvalue, which is defined as the eigenvalue of Eq. 13, obtained by setting frequencies μ and ω to zeros. Since <math id="M118"><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> is obtained by solving the transport EVP, Eq. 52 is equivalent to calculating the reactivity using the angular flux. It eliminates the necessity of reconstructing and storing the angular flux. In addition, solving for the dynamic eigenvalue does not increase computational cost, because the SCM requires solving EVPs iteratively, as shown in Fig. 1 Fig. 1 . Denoting the reactivity evaluation approach using Eq. 52 as Method 1 and using the weight of the scalar flux as Method 2, the comparison of the two methods is shown in Fig. 6 Fig. 6 . The results indicate that the reactivity evaluation using the scalar flux underestimates the reactivity by -15.0% and -10.0% in TD3-4 and TD5-1 of the C5G7-TD problem respectively, which is consistent with the results reported by [ 5 5 ]. Overall PCQ-SCM scheme 3 In the PCQ-SCM solver of VITAS, the initial adjoint scalar flux <math id="M119"><mrow><msubsup><mi>ϕ</mi><mi>g</mi><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><mtext>r</mtext><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></math> and the forward scalar flux <math id="M120"><mrow><msub><mi>ϕ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mtext>r</mtext><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></math> are obtained by solving the initial EVPs. The initial amplitude function T ( t 0 ) is determined using P ( t 0 ). Then, in each coarse time step [ t n -1 , t n ], the TEVP and EPKE are solved as follows: 1. Solve the predicted scalar flux <math id="M121"><mrow><msub><mi>ϕ</mi><mrow><mi>p</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mtext>r</mtext><mn>,</mn><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow></math> , flux amplitude P p ( t n ), and precursor concentrations C p,i ( 𝒓,t ) using the TEVP in Eq. 13 with k-ω iteration at t n . In the first iteration, the frequencies μ and ω are set to zero. Then, the eigenvalue obtained in the first iteration is the dynamic eigenvalue <math id="M122"><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow></math> required in Eq. 52. 2. Determine the reactivity ρ ( t n ) using Eq. 52. 3. Calculate the predicted shape function <math id="M123"><mrow><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mrow><mi>p</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mtext>r</mtext><mn>,</mn><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow></math> with the constraint in Eq. 45 as <math display="block" id="M124"><mtable><mtr><mtd><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mrow><mi>p</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><msub><mi>ϕ</mi><mrow><mi>p</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>×</mo><mfrac><mrow><mi>T</mi><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></msubsup><mrow><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub></mrow></mstyle><mtext>d</mtext><mi>r</mi></mrow></mstyle><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mtext>d</mtext><mtext>Ω</mtext><msubsup><mi>ϕ</mi><mi>g</mi><mo>*</mo></msubsup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mfrac><mn>1</mn><mrow><msub><mi>v</mi><mi>g</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mo>)</mo></mrow></mrow></mfrac><msub><mi>ϕ</mi><mrow><mi>p</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow></mfrac><mo>.</mo></mtd></mtr></mtable></math> 4. Calculate the KPs at t n with Eqns. 48-51, and solve the EPKE to obtain the amplitude function T ( t n ) at the fine time step using the interpolated KPs from the values at t n -1 and t n . 5. Determine the corrected flux amplitude P c ( t n ) and the scalar flux <math id="M125"><mrow><msub><mi>ϕ</mi><mrow><mi>c</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow></math> , as follows: <math display="block" id="M126"><mrow><msub><mi>P</mi><mi>c</mi></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mi>T</mi><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mi>P</mi><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>T</mi><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></mfrac><mo>,</mo></mrow></math> <math display="block" id="M127"><mtable><mtr><mtd><msub><mi>ϕ</mi><mrow><mi>c</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><msub><mi>P</mi><mi>c</mi></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>P</mi><mi>p</mi></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow></mfrac><msub><mi>ϕ</mi><mrow><mi>p</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>=</mo><msub><mi>P</mi><mi>c</mi></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><msub><mover accent="true"><mi>ϕ</mi><mo>^</mo></mover><mrow><mi>p</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold-italic"><mi>r</mi></mstyle><mn>,</mn><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>.</mo></mtd></mtr></mtable></math> 6. Re-calculate the precursor concentrations with the corrected scalar flux using Eq. 35. Then, iterate Steps 1 to 6 in the next coarse time step until the end of the transient event. The overall flowchart of the PCQ-SCM is shown in Fig. 7 Fig. 7 . In addition, in the PCQ-SCM solver, the amplitude frequency ω T is approximated as a constant in [ t n -1 , t n ], to improve the numerical stability in the case of large time steps. Although the linear approximation of Eq. 33 shows higher accuracy to the constant approximation, the amplitude is corrected with the solution of the EPKE. Therefore, the constant approximation offers better stability in large time-step cases, without affecting numerical accuracy. In the PCQ-SCM solver, Eq. 34 is replaced by <math display="block" id="M128"><mtable><mtr><mtd><msub><mi>ϕ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><msub><mi>ϕ</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow><msup><mi>e</mi><mrow><msub><mi>ω</mi><mi>T</mi></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mi>Δ</mi><msub><mi>t</mi><mi>n</mi></msub><mo>+</mo><msub><mi>ω</mi><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mtext>S</mtext><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mi>Δ</mi><msub><mi>t</mi><mi>n</mi></msub></mrow></msup></mtd></mtr><mtr><mtd><mo>=</mo><msub><mover accent="true"><mi>ϕ</mi><mo>^</mo></mover><mrow><mi>ν</mi><mi>k</mi><mn>,</mn><mi>g</mi></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><msup><mi>e</mi><mrow><msub><mi>ω</mi><mi>T</mi></msub><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mi>Δ</mi><msub><mi>t</mi><mi>n</mi></msub></mrow></msup><mo>.</mo></mtd></mtr></mtable></math> C5G7-TD benchmark problem 1 The C5G7-TD transient benchmark is derived from the C5G7 benchmark model [ 31 31 , 32 32 ]. It includes both 2D and 3D models that are used to verify the transient capability for modeling a heterogeneous light-water reactor without thermal feedback. In the following sections, the geometry, materials, and transient exercises in C5G7-TD are described. Geometry description 2 The benchmark problem simulates a light-water reactor with eight uranium oxide (UO2) assemblies and eight mixed oxide (MOX) assemblies surrounded by a water reflector. The model is simplified into a quarter core using reflective boundaries on the north and west boundaries and vacuum boundaries on the south and east boundaries, as shown in Fig. 8-a Fig. 8-a . The four assemblies are numbered 1-4, respectively. The dimensions of the 2D model are 64.26 cm × 64.26 cm. The 3D model has the same planar layout as the 2D model; however, in the axial direction, two identical water reflectors are added above and below the active region of the core, as shown in Fig. 8-a Fig. 8-a . The height of the 3D configuration is 171.36 cm while that of the active region is 128.52 cm. The UO2 assembly and MOX assembly have identical geometric configurations but different fuel-pin compositions. Their dimensions are 21.42 cm × 21.42 cm, with a 17 × 17 pin layout. Each assembly consists of 289 pin cells, including 264 fuel pins, 24 guide tubes, and one fission chamber, as shown in Fig. 8-b Fig. 8-b . The UO2 assembly contains one type of fuel pin, whereas the MOX assembly contains three types of fuel pins with enrichments of 4.3%, 7.0%, and 8.7%, respectively. Each pin cell has a pitch of 1.26 cm and is simplified into two zones, as shown in Fig. 8-c Fig. 8-c . Zone 1 is the circular area with the radius of 0.54 cm. It is filled with fuel, control rod, or fission chamber. Zone 2 is located outside the circular zone and is filled with a moderator. Material description 2 The C5G7-TD benchmark adopts multi-group macroscopic cross-sections and KPs for eight materials: UO 2 fuel, 7.0% MOX fuel, 4.3% MOX fuel, 8.7% MOX fuel, guide tube, fission chamber, moderator, and control rod[ 32 32 ]. The cross-sections are provided in a 7-group structure including transport-corrected total cross-sections ∑ tg , absorption cross-sections ∑ ag , fission cross-sections ∑ fg , fission spectra χ g , fission neutron yields ν , and scattering cross-sections <math id="M129"><mrow><mstyle displaystyle="true"><mo>∑</mo><mrow><msub><mrow/><mrow><mtext>s</mtext><mi>g</mi><msup><mi>g</mi><mo>′</mo></msup></mrow></msub></mrow></mstyle></mrow></math> . The KPs are provided in an 8-delayed-group structure including neutron velocities v g , delayed neutron fission spectra χ ig , delayed neutron fractions β i , and decay constants λ i . For 3D control-rod movement problems, such as TD4 problems, to reduce the rod-cusping effects, the effective volume fraction proposed by[ 6 6 ] is applied to homogenize the cross-sections of the guide tube and control rod in the partially rodded node. The effective volume fraction ( y ) of the rod is a polynomial function of the original volume fraction ( x ), which is the ratio of the inserted rod length to the node height. The polynomial function is <math display="block" id="M130"><mtable><mtr><mtd><mi>y</mi><mo>=</mo><mn>0.3867848</mn><mi>x</mi><mo>+</mo><mn>0.1707878</mn><msup><mi>x</mi><mn>2</mn></msup><mo>+</mo><mn>0.9887881</mn><msup><mi>x</mi><mn>3</mn></msup></mtd></mtr><mtr><mtd><mo>−</mo><mn>5.9535775</mn><msup><mi>x</mi><mn>4</mn></msup><mo>+</mo><mn>25.1805898</mn><msup><mi>x</mi><mn>5</mn></msup><mo>−</mo><mn>63.1841252</mn><msup><mi>x</mi><mn>6</mn></msup></mtd></mtr><mtr><mtd><mo>+</mo><mn>98.3898802</mn><msup><mi>x</mi><mn>7</mn></msup><mo>−</mo><mn>92.4982596</mn><msup><mi>x</mi><mn>8</mn></msup><mo>+</mo><mn>48.0373368</mn><msup><mi>x</mi><mn>9</mn></msup></mtd></mtr><mtr><mtd><mo>−</mo><mn>10.5182051</mn><msup><mi>x</mi><mrow><mn>10</mn></mrow></msup></mtd></mtr></mtable></math> Figure 9 Figure 9 shows the relationships between the original and effective volume fractions, for Eq. 57 and linear assumption, respectively. Transient description 2 The C5G7-TD benchmark problem consists of six exercises, labelled TD0-TD5. Each exercise includes a number of cases that simulate the movement of the control-rod banks or changes in the moderator density. To model the control-rod movements, the control rods in an assembly are referred to as rod banks. Control rods in the same bank are simultaneously inserted or withdrawn. As illustrated in Table 1 Table 1 , TD0, TD1, and TD2 are 2D exercises that require the adjustment of different control-rod banks in various ways. TD0 contains five test problems. In TD0, the control rods are abruptly inserted into the active region (10% of the core height) at the initial time and remain there for 1 s. Then, the control rods are partially extracted (5% of the core height) and remain so for 1 s. After 2 s, the control rods are completely withdrawn from the active region. 1 1 Test problem 1 1 TD0 1 1 TD1 1 1 TD2 1 1 1 1 1 Bank 1 1 1 Bank 1 1 1 Bank 1 1 1 2 1 1 Bank 3 1 1 Bank 3 1 1 Bank 3 1 1 3 1 1 Bank 4 1 1 Bank 4 1 1 N/A 1 1 4 1 1 Bank 1, 3 and 4 1 1 Bank 1, 3 and 4 1 1 N/A 1 1 5 1 1 Bank 1-4 1 1 Bank 1-4 1 1 N/A TD1 and TD2 each contain five and three test problems, respectively. In TD1 and TD2, control rods move linearly at a constant speed. From the starting point at t =0 s, the control rods are inserted to the maximal depth at t =1 s and then returned to the starting point at t =2 s. The maximal insertion depths for TD1 and TD2 are 1% and 10% of the core heights, respectively. In the 2D model, the control-rod movement is simulated as a linear substitution of the moderator-filled guide tube material with the control-rod material in Zone 1 of the specified pin cells. TD3 is a 2D exercise simulating the reactivity insertion caused by the moderator density change in the reactor core. It has four test problems. During the first second, the moderator density in all assemblies declines simultaneously and at the same rate, from the nominal value to the minimal value. After t =1 s, the density increases to a nominal value within the following second: These four test problems differ in the ratio of the minimum to the nominal value, as shown in Fig. 10 Fig. 10 . TD4 is a 3D exercise that simulates the process of insertion and withdrawal of the control rods. It contains five problems. In the initial phase, all the control rods are placed in the top water reflector outside the active region. These five problems involve the continual insertion and removal of various combinations of control-rod banks. Fig. 11 Fig. 11 depicts the control-rod banks and the relative insertion length to the active-core length, for each problem. Similar to TD3, TD5 is a 3D transient event initiated by the moderator density change with all control rods fully withdrawn. The moderator density in the same assembly changes simultaneously. For each of the four test problems, the density of the moderator varies for each of the four assemblies, as shown in 10. Numerical results 1 Unless otherwise noted, all benchmark problems were simulated using VITAS, on a workstation, with 16 processors for the 2D problems and 48 processors for the 3D problems. Each fuel pin was meshed with three radial rings for the fuel zone, one radial ring for the moderator zone, and eight azimuthal sectors as determined by a preliminary h - p sensitivity analysis[ 11 11 ]. In total, the mesh comprised 32 elements, as shown in Fig. 2 Fig. 2 . In the 3D calculation, the axial direction was discretized into 32 layers, and the thickness of each layer was 5.355 cm for the TD5 exercise. For the TD4 exercise, the axial direction was discretized into 48 layers to better describe the control-rod movement. For the spatial expansion within the node, the x - y plane was expanded with continuous FE basis functions of the second order, and the axial direction was expanded with orthonormal polynomials of the second order. In addition, the interface was approximated using orthonormal polynomials of the second order. For the angular expansion, a 25 × 25 square Legendre-Chebyshev cubature was used to evaluate the angular integrals within the node, and the spherical harmonic functions of P 23 P 3 and P 3 were implemented on the lateral and axial interfaces, respectively. Table 2 Table 2 summarizes these parameters. To reduce the computational effort, the 2D/1D approach (for the 3D problems) and a quasi-reflected interface condition of up to P3 were applied. 1 1 Model parameters 1 1 Value 1 1 Fuel-pin mesh (fuel-zone rings/moderator rings/azimuthal sectors) 1 1 3/1/8 1 1 Volume spatial expansion order 1 1 2 1 1 Surface spatial expansion order 1 1 2 1 1 Volume angular integrals 1 1 25×25 Square Legendre-Chebyshev cubature 1 1 Surface P N order (lateral interfaces/axial interfaces) 1 1 P 23 P 3 / P 3 1 1 Number of axial meshes in 3D 1 1 48 (TD4), 32 (TD5) Steady-state results 2 For steady-state calculations, Table 3 Table 3 compares the eigenvalue results for the 2D and 3D initial eigenvalue problems calculated using different solvers[ 5 5 , 6 6 ]. VITAS results were consistent with those obtained using the other solvers. For the 2D case, the difference between the eigenvalues of VITAS and MCNP5 was 29 pcm, less than the Monte Carlo uncertainty of 0.07%, whereas for the 3D case, the difference between VITAS and PROTEUS-MOC was 25 pcm. Minor deviations indicate that the model was accurate with respect to the initial conditions. In addition, the eigenvalues of the forward and adjoint calculations were consistent, with a difference of 1 pcm. 1 1 Code 1 1 Method 1 1 2D eigenvalue 1 1 3D eigenvalue 1 1 MCNP5 1 1 MC 1 1 1.18646±0.07% 1 1 PROTEUS-MOC 1 1 3D MOC 1 1 1.18651 1 1 1.16469 1 1 MPACT 1 1 2D/1D MOC 1 1 1.18666 1 1 1.16359 1 1 NECP-X 1 1 2D/1D MOC 1 1 1.18695 1 1 PANDAS-MOC 1 1 2D/1D MOC 1 1 1.18631 1 1 1.16512 1 1 VITAS (Forward) 1 1 VNM 1 1 1.18675 1 1 1.16494 1 1 VITAS (Adjoint) 1 1 VNM 1 1 1.18674 1 1 1.16493 The normalized pin-wise fission rate distributions for the 2D and 3D problems are shown in Fig. 13 Fig. 13 and Fig. 14 Fig. 14 , respectively. The normalization was implemented by scaling the initial total fission rate to unity. Verification of the direct-SCM solver 2 To verify the direct-SCM solver in VITAS, the VITAS results were compared with the results of MPACT[ 4 4 ] and PANDAS-MOC[ 6 6 ]. MPACT uses a novel transient multilevel (TML) approach that involves three calculation levels: 3D-transport, 3D-coarse mesh finite difference (CMFD), and EPKE, to enhance computational efficiency. The coarse time-step solution of the 3D-transport calculation was corrected with the fine time-step solutions of the CMFD and EPKE calculations. In the PANDAS-MOC approach, the transient solution is solved by coupling the 2D MOC and 1D nodal methods accelerated by multi-level CMFD methods. Owing to the large reactivity insertion, the core power in the TD3-4 problem fluctuates rapidly. Consequently, the TD3-4 problem was chosen as a representative case for the time-step size sensitivity analysis. The fractional core fission rates obtained with time-step sizes ranging from 10 ms to 200 ms are shown in Fig. 15 Fig. 15 . Fig. 15 Fig. 15 illustrates the reference solution and the percentage differences between the reference solution and solutions with larger time steps. The reference solution was obtained with a time step of 1 ms using the direct-SCM solver. The percentage difference ε was calculated as <math display="block" id="M131"><mrow><mi>ε</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mi>P</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><msup><mi>P</mi><mo>⋆</mo></msup><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow><mrow><mi>P</mi><mrow><mo>(</mo><mrow><msub><mi>t</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow></mfrac><mo>×</mo><mn>100</mn><mi>%</mi><mo>,</mo></mrow></math> where the superscript * denotes the reference solution. Fig. 15 Fig. 15 indicates that as the time step increases to Δ t =100 ms or 200 ms, fission-rate oscillations appear in the asymptotic regime ( t =2-3 s). Similarly, the oscillation of the percentage difference is shown in Fig. 15 Fig. 15 . The numerical oscillation is attributed to the linear approximation of the amplitude frequency for large time steps. In the asymptotic regime, control rods are completely withdrawn from the active region, and the fission rate increases gradually owing to the delayed neutrons from precursors. Owing to the slower growth rate of the fission rate, the numerical oscillation becomes dominant in the asymptotic regime. This also explains why there are no observable oscillations when the fission rate changes rapidly during t =0-2 s. When the time step was reduced to 20 ms, the oscillation decayed rapidly. To further quantify the effect of the time-step size, Table 4 Table 4 summarizes the highest percentage differences for different time-step sizes. The root mean square (RMS) error was calculated using Eq. 59. According to Table 4 Table 4 , for the time-step size of 20 ms, the solution agreed well with the reference solution, with the RMS error under 0.1%. In addition, the benchmark suggested that the time-step size during the transient should not exceed 25 ms[ 32 32 ]. Therefore, 20 ms may be an effective time-step size for the direct-SCM solver to reduce the computational cost without losing numerical resolution. If not otherwise specified, the time-step size used by the direct-SCM solver in the following section was Δ t =20 ms. <math display="block" id="M132"><mrow><mtext>RMS</mtext><mo>=</mo><msqrt><mrow><mfrac><mrow><msup><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>0</mn></mrow><mi>N</mi></msubsup><mrow><mrow><mo>[</mo><mrow><mi>P</mi><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mo>−</mo><msup><mi>P</mi><mo>⋆</mo></msup><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow></mrow></mstyle></mrow><mn>2</mn></msup></mrow><mi>N</mi></mfrac></mrow></msqrt></mrow></math> 1 1 Step size (ms) 1 1 Max. Diff. (%) 1 1 RMS error (%) 1 1 10 1 1 0.463% 1 1 0.035% 1 1 20 1 1 0.826% 1 1 0.099% 1 1 100 1 1 1.370% 1 1 0.499% 1 1 200 1 1 2.889% 1 1 1.288% Results of 2D studies 3 In this section, the results of the 2D exercises with ramp perturbations are discussed, including TD-1, TD-2, and TD-3. To assess the accuracy of the simulations, the VITAS results were compared with those of MPACT[ 4 4 ] and PANDAS-MOC[ 6 6 ]. Rela tive differences were evaluated using MPACT and PANDAS-MOC as references. The time-step sizes employed in MPACT and PANDAS-MOC were 10 ms and 20 ms, respectively. It should also be noted that the PANDAS-MOC method uses longer time steps of 100 ms and 50 ms in the asymptotic regime, for the TD1/2 and TD3 exercises. The fractional fission rates for the TD1 and TD2 exercises are shown in Fig. 16 Fig. 16 . For these problems, continuous changes in fractional fission rates were observed owing to the constant-speed movement of control rods. Despite the fact that the speeds of insertion and withdrawal were identical, the increase rate of the fission rate after 1 s was smaller than the decrease rate before 1 s, owing to the reduction in the precursor concentration. The fission rate then asymptotically increased to a value smaller than the initial value. Owing to the deeper insertion of the control rods in the TD2 case, the fission rate changed more rapidly in that case. Fig. 17 Fig. 17 and Fig. 18 Fig. 18 compare the direct-SCM solver’s results with those obtained using MPACT and PANDAS-MOC. Comparisons reveal that the results obtained using VITAS with the direct-SCM agree well with those obtained using the other codes, for the 2D control-rod movement problems. Figure 19 Figure 19 shows the fractional fission rate histories, for the TD3 case. Owing to the similar patterns of perturbations caused by the moderator density changes, the fission rates exhibit similar tendencies. As seen in Fig. 20 Fig. 20 , the direct-SCM solver’s results for the TD3 case were also consistent with those obtained using the other codes. Table 5 Table 5 summarizes the RMS errors and maximal percentage differences, for the TD1, TD2, and TD3 cases. Compared with the other cases, the solutions for the TD2-1 and TD3-4 cases exhibit greater deviations, with a maximal percentage difference of 1.0%. As shown in Fig. 18 Fig. 18 and Fig. 20 Fig. 20 , the maximal differences appeared after a few time steps since the introduction of the reactivity. Because the reactivity insertion was substantial and the fission rate changed rapidly in the TD2-1 and TD3-4 cases, different temporal discretization schemes could lead to larger discrepancies. However, these discrepancies would be under 1.0%, implying that the direct-SCM solver of VITAS can perform accurate calculations for 2D ramp-perturbation problems. 1 1 Problem 1 1 v.s. MPACT Max. Diff. (%) 1 1 v.s. MPACT RMS Err. (%) 1 1 v.s. PANDAS Max. Diff. (%) 1 1 v.s. PANDAS RMS Err. (%) 1 1 1-1 1 1 -0.13 1 1 0.08 1 1 0.14 1 1 0.06 1 1 1-2 1 1 -0.16 1 1 0.08 1 1 -0.16 1 1 0.08 1 1 1-3 1 1 -0.14 1 1 0.07 1 1 -0.14 1 1 0.06 1 1 1-4 1 1 -0.17 1 1 0.07 1 1 0.17 1 1 0.06 1 1 1-5 1 1 -0.17 1 1 0.08 1 1 0.27 1 1 0.08 1 1 2-1 1 1 0.91 1 1 0.10 1 1 0.91 1 1 0.10 1 1 2-2 1 1 -0.21 1 1 0.11 1 1 -0.23 1 1 0.15 1 1 2-3 1 1 -0.12 1 1 0.05 1 1 -0.11 1 1 0.04 1 1 3-1 1 1 0.25 1 1 0.11 1 1 0.29 1 1 0.16 1 1 3-2 1 1 0.45 1 1 0.11 1 1 0.56 1 1 0.16 1 1 3-3 1 1 0.72 1 1 0.11 1 1 0.72 1 1 0.15 1 1 3-4 1 1 0.97 1 1 0.11 1 1 0.97 1 1 0.14 Results of 3D studies 3 In this section, the results of the 3D TD4 and TD5 exercises are discussed. For the TD4 and TD5 exercises, VITAS employed time-step sizes of 25 ms and 20 ms, respectively. MPACT and PANDAS-MOC employed time-step sizes of 25 ms and 20 ms, respectively. In the asymptotic regime, the PANDAS-MOC time-step size was increased to 100 ms. The fractional fission-rate history for the TD4 exercise is shown in Fig. 21 Fig. 21 . For the TD4-4 and TD4-5 problems, the fission rates are more complicated because multiple control-rod banks were inserted and withdrawn simultaneously. Fig. 22 Fig. 22 compares the VITAS results with those obtained using MPACT and PANDAS-MOC. The MPACT results agree well with the results obtained using the other codes, with the VITAS results being the most comparable to the MPACT results. Oscillations in the percentage difference over a certain period can be observed in the stage of the control-rod movement. This was attributed to the different de-cusping techniques implemented in the different codes. As shown in Table 6 Table 6 , the maximal difference between the MPACT and VITAS results was -0.97%, whereas the maximal difference between the PANDAS-MOC and VITAS results was 1.57%. The discrepancies between these codes in the TD4 case were caused not only by the temporal or spatial discretization schemes but also by the different de-cusping techniques. 1 1 Problem 1 1 v.s. MPACT Max. Diff. (%) 1 1 v.s. MPACT RMS Err. (%) 1 1 v.s. PANDAS Max. Diff. (%) 1 1 v.s. PANDAS RMS Err. (%) 1 1 4-1 1 1 -0.97 1 1 0.28 1 1 1.57 1 1 0.69 1 1 4-2 1 1 -0.77 1 1 0.28 1 1 0.68 1 1 0.28 1 1 4-3 1 1 -0.77 1 1 0.29 1 1 0.68 1 1 0.31 1 1 4-4 1 1 -0.75 1 1 0.34 1 1 0.74 1 1 0.25 1 1 4-5 1 1 -0.97 1 1 0.39 1 1 1.57 1 1 0.62 Figure 23 Figure 23 presents the evolution of the relative fission rates for the TD5 exercise, in which the moderator density varied linearly in different fuel assemblies, as shown in Fig. 12 Fig. 12 . During the first two seconds, the fission rates decreased to minimal values as the moderator densities decreased. In the subsequent two seconds, the moderator densities were restored to nominal values, and the fission rate increased gradually. Fig. 24 Fig. 24 compares the results obtained using several different codes. The comparison demonstrates that the solutions obtained using the direct-SCM agreed well with the MPACT and PANDAS-MOC results. As shown in Table 7 Table 7 , compared with the corresponding MPACT and PANDAS-MOC solutions for the TD5-2 problem, the results obtained using the direct-SCM exhibited maximal differences of 0.87% and 0.85%, respectively. 1 1 Problem 1 1 v.s. MPACT Max. Diff.(%) 1 1 v.s. MPACT RMS Err.(%) 1 1 v.s. PANDAS Max. Diff.(%) 1 1 v.s. PANDAS RMS Err.(%) 1 1 5-1 1 1 0.69 1 1 0.43 1 1 0.74 1 1 0.52 1 1 5-2 1 1 0.87 1 1 0.45 1 1 0.85 1 1 0.54 1 1 5-3 1 1 0.85 1 1 0.38 1 1 0.78 1 1 0.47 1 1 5-4 1 1 0.33 1 1 0.10 1 1 0.34 1 1 0.12 Step-perturbation results 3 The TD0 results obtained using VITAS and MPACT were compared to verify the direct-SCM solver for step-perturbation problems. Typically, a step-perturbation problem requires finer time steps to capture the rapid changes induced by the step reactivity[ 7 7 ]. For the PCQM solver, the time steps are further discretized into several micro time steps for the EPKE solution. Therefore, the PCQM solver can simulate the step reactivity with rather large time steps, without any significant loss of accuracy. However, the direct transient solver requires finer time steps. To compare the solutions obtained using the direct-SCM solver with those obtained using MPACT in the TD0 exercise, variable time steps were employed in VITAS. In the TD0 exercise, reactivity was introduced abruptly at t =0,1,2 s, while the cross-sections remained constant for the remaining transient activity. Consequently, smaller time steps could be employed immediately after the step reactivity change; otherwise, larger time steps could be used. The time-step size used in MPACT was 10 ms. Fig. 25 Fig. 25 illustrates the variable time-step size employed in VITAS, which was 1 ms during the first 20 ms following a step reactivity change, and was set to 20 ms for the remainder of the transient process. The fractional core fission rates for the TD0 exercise are shown in Fig. 26 Fig. 26 . The fission rates decreased at 0 s and 1 s as a result of the abrupt insertion of control rods, and then increased promptly at 2 s with the withdrawal of the control rods. After 2 s, all the control rods were completely out of the active region, and the fission rates increased asymptotically to a value smaller than the initial value. Fig. 27 Fig. 27 compares the VITAS and MPACT solutions, while Table 8 Table 8 summarizes the error. The maximal relative difference across all of the step-perturbation problems was approximately 0.21% for the TD0-2 problem, whereas the differences for the other problems were within 0.12%. These results demonstrate that the solutions obtained using VITAS and MPACT were highly consistent, indicating that the direct-SCM solver can generate accurate solutions for step-perturbation problems. 1 1 Problem 1 1 v.s. MPACT Max. Diff. (%) 1 1 v.s. MPACT RMS Err. (%) 1 1 0-1 1 1 -0.12 1 1 0.05 1 1 0-2 1 1 -0.20 1 1 0.11 1 1 0-3 1 1 -0.10 1 1 0.05 1 1 0-4 1 1 -0.11 1 1 0.05 1 1 0-5 1 1 -0.11 1 1 0.05 Performance of the PCQ-SCM solver 2 The PCQ-SCM solver was compared with the direct-SCM solver using three problems: TD0-5, TD3-4, and TD5-1. The results of this analysis are presented below. The numerical accuracy and computational cost were used as metrics for evaluating the methods’ performance. Figure 28 Figure 28 shows that the reactivity evaluated using VITAS was consistent with that evaluated using MPACT. Figure 29 Figure 29 compares the results obtained using the direct-SCM and PCQ-SCM solvers, for the TD0-5 problem. For the direct-SCM solver the solution was obtained with a constant time step of 1 ms, while for the PCQ-SCM solver the solution was obtained using a variable time step. To capture the step reactivity, the step size was 1 ms during the first 1 ms following the step reactivity change, whereas step sizes of 100 or 200 ms were used in the remainder. As shown in Fig. 29 Fig. 29 , the solutions obtained using the PCQ-SCM solver were consistent with those obtained using the direct-SCM solver, with the exception of the first 1 ms following the step reactivity change. As shown in Table 9 Table 9 , the percentage differences were within 0.08% for the both time-step sizes. Note that when calculating the RMS error and maximal percentage difference shown in Table 9 Table 9 , the results for the first 1 ms following the reactivity change were not included. The significant discrepancies following the reactivity change were caused by errors associated with the direct-SCM solution. The step size of 1 ms in the direct solver was insufficient for tracing the step response of the fission rate; however, for the PCQ-SCM solver, the first 1 ms was further subdivided into 100 micro-steps. Using micro-steps (0.01 ms) it was possible to simulate the step response accurately. 1 1 Method 1 1 Δ t (ms) 1 1 Number of TEVPS 1 1 Max. Diff. (%) 1 1 RMS Err. (%) 1 1 Computing time (min) 1 1 Speedup factor 1 1 Direct 1 1 1 1 1 5000 1 1 N/A 1 1 N/A 1 1 14 1 1 N/A 1 1 PCQ 1 1 100 1 1 53 1 1 0.08 1 1 0.04 1 1 7 1 1 2.0 1 1 PCQ 1 1 200 1 1 28 1 1 0.08 1 1 0.04 1 1 6 1 1 2.4 In addition, Table 9 Table 9 lists the compute times for the different solvers. Compared with the direct solution, the number of TEVPs in the PCQ-SCM solution with the step size of 100 ms decreased from 5000 to 53, and the compute time decreased from 14 min to 7 min. The incommensurate decrease in the compute time was attributed to the fact that the compute time was dominated by the response matrix formation in 2D calculations. In the transient solvers of VITAS, the response matrices were updated only when the cross-sections were altered. For the TD0-5 problem, increasing the time-step size to 100 ms accelerated the iteration for flux solving, but not the response matrix formation. Therefore, the reduction in the computational cost owing to the PCQ-SCM solver was insignificant for the TD0 exercises. Fig. 30 Fig. 30 illustrates the reactivity history of the TD3-4 problem, for VITAS and MPACT. Although the cross-sections varied linearly in this problem, the reactivity curve versus time was nonlinear. This behavior was attributed to the spatial effects caused by variations in the flux shape. Fig. 31 Fig. 31 compares the solutions of the direct-SCM and PCQ-SCM solvers, for the TD3-4 problem. The solutions of the direct-SCM solver were obtained with fixed time-step sizes of 1 and 100 ms, and the PCQ-SCM calculations were performed with the time-step sizes of 100, 500, and 1000 ms, respectively. Fig. 31 Fig. 31 shows that as the time step of the PCQ-SCM decreased, the PCQ-SCM solution gradually converged to the direct solution. The maximal percentage difference for the solution with the step size of 100 ms was -0.23%, as shown in Table 10 Table 10 . This was less than the direct solution of 10 ms, as shown in Table 4 Table 4 . In addition, Fig. 31 Fig. 31 indicates that in the direct solution with 100 ms, there were significant fluctuations in the asymptotic regime, leading to the maximal difference exceeding 1.0%. This indicates that the PCQ-SCM solver can accurately predict the fission rate and eliminate numerical fluctuations with large time steps. 1 1 Method 1 1 Δ t (ms) 1 1 Number of TEVPS 1 1 Max. Diff. (%) 1 1 RMS Err. (%) 1 1 Computing time (min) 1 1 Speedup factor 1 1 Direct 1 1 1 1 1 10000 1 1 N/A 1 1 N/A 1 1 6573 1 1 N/A 1 1 PCQ 1 1 100 1 1 100 1 1 -0.23 1 1 0.12 1 1 83 1 1 78.9 1 1 PCQ 1 1 500 1 1 20 1 1 -0.42 1 1 0.21 1 1 16 1 1 405.0 1 1 PCQ 1 1 1000 1 1 10 1 1 -1.21 1 1 0.65 1 1 11 1 1 598.6 Meanwhile, the number of the TEVP calculations for the PCQ-SCM calculation with the 100 ms step decreased from 10000 to 100, where the linear speedup was 100, and the compute time decreased by a factor of 78.9. The speedup of the TEVP calculations was not proportional to the speedup of the total compute time. This was because larger time steps made it more difficult for the amplitude frequency, fission source, and flux to converge in the k-ω iteration. However, the speedup of the overall compute time was close to the linear speedup of the TEVP calculations. This implies that the PCQ-SCM solver efficiently reduced the compute time without a significant loss of numerical accuracy, for the TD3-4 problem. Fig. 32 Fig. 32 shows the reactivity history for TD5-1, as evaluated by VITAS, MPACT, and PANDAS-MOC. Notably, the reactivity in Fig. 32 Fig. 32 was defined as <math id="M133"><mrow><mi>ρ</mi><mo>/</mo><mover accent="true"><mi>β</mi><mo stretchy="true">¯</mo></mover></mrow></math> with the units of $, whereas it was defined as ρ with the units of pcm in Fig. 32 Fig. 32 . Fig. 32 Fig. 32 reveals that, compared with the reactivity of MPACT and PANDAS-MOC, VITAS underestimated <math id="M134"><mrow><mi>ρ</mi><mo>/</mo><mover accent="true"><mi>β</mi><mo stretchy="true">¯</mo></mover></mrow></math> . However, in Fig. 32 Fig. 32 , the reactivity ρ evaluated by VITAS was consistent with that estimated by PANDAS-MOC. Consequently, the deviations in Fig. 32 Fig. 32 were attributed to the discrepancies of the effective delayed neutron fractions <math id="M135"><mrow><mover accent="true"><mi>β</mi><mo stretchy="true">¯</mo></mover></mrow></math> computed by different codes. Figure 33 Figure 33 compares the results obtained using the direct-SCM and PCQ-SCM solvers, for the TD5-1 problem. Direct solutions were obtained with the time steps of 1 and 100 ms. The macro time steps used in the PCQ-SCM calculation were 100 ms and 500 ms. As shown in Fig. 33 Fig. 33 , during the first four seconds, the PCQ-SCM calculations underestimated the fission rate with respect to the direct solution. However, when the cross-sections reverted to the nominal value in the asymptotic regime, the fission rate values became closer to those of the reference solution. These discrepancies are listed in Table 11 Table 11 . The maximal difference between the direct solution and PCQ-SCM solution with the step size of 500 ms was -1.01%, and the RMS error was 0.66%. When the macro time step was set to 100 ms, the solution and error showed no significant changes. In addition, as shown in Fig. 33 Fig. 33 , compared with the PCQ solutions, the direct solution with the time step of 100 ms exhibited small discrepancies, below 0.5%. The discrepancies between the direct and PCQ-SCM solutions could be attributed to the inconsistent formulation of KPs in VITAS transient models [ 5 5 ]. 1 1 Method 1 1 Δ t (ms) 1 1 Number of TEVPS 1 1 Max. Diff. (%) 1 1 RMS Err. (%) 1 1 Computing time (min) 1 1 Speedup factor 1 1 Direct 1 1 1 1 1 500 1 1 N/A 1 1 N/A 1 1 6820 1 1 N/A 1 1 PCQ 1 1 100 1 1 50 1 1 -1.00 1 1 0.66 1 1 271 1 1 25.2 1 1 PCQ 1 1 500 1 1 10 1 1 -1.01 1 1 0.65 1 1 99 1 1 68.9 Table 11 Table 11 further summarizes the speedup for the PCQ-SCM solutions. This shows that the speedup factors for the TD5-1 problem were significantly lower than those for the TD3-4 problem. The speedup discrepancy was owing to the different fractions of compute time for 2D and 3D problems. In VITAS, 99% of the overall runtime was consumed by obtaining the response matrices and solving the transport equation. Therefore, the following runtime analysis focused on matrix formation and equation solving. Fig. 34 Fig. 34 shows the fractional runtime of the matrix formation and equation solving for the TD5-1 problem. When the time step was 1 ms, matrix formation and equation solving accounted for 57.16% and 42.80% of the runtime, respectively. However, when using the time steps of 100 ms or 500 ms for the PCQ-SCM solutions, equation solving became dominant over the total runtime. Fig. 34 Fig. 34 shows the speedups for the matrix formation, equation solving, and total calculation. The speedup for the matrix formation was close to linear, whereas the speedup for equation solving was significantly sublinear. As previously noted, the fission source and flux required more iterations to converge for larger time steps. Therefore, the equation-solving speedup was not linear. As the time step increased, the run time became dominated by equation solving, and the speedup of the total calculation was primarily determined by the speedup of equation solving. Conclusion 1 This work discussed the direct-SCM and PCQ-SCM transient models of the VITAS code, as well as the verification of these models using the C5G7-TD benchmark. In the formulation of the direct-SCM, the SCM transformed the transient transport equation into a TEVP, and then the TEVP was solved by a heterogeneous VNM. In addition, the PCQ-SCM was proposed and developed to reduce the computational cost. The predicted neutron flux in the PCQ-SCM was obtained by solving the TVEP, and was then corrected by the EPKE solution. For the formulation of the EPKE, methods for computing the reactivity and adjoint flux within the heterogeneous VNM were discussed. Owing to the second-order even-parity form of the transport equation in the VNM, the adjoint equation was easily solved by creating an adjoint source without generating additional response matrices. To calculate the reactivity, the exact reactivity was determined using the dynamic eigenvalue acquired from the TEVP, without the need to reconstruct the angular flux or apply perturbation theory methods. To verify the direct-SCM solver, four C5G7-TD 2D and two 3D benchmark exercises were considered. The direct-SCM solver’s transient results agreed well with those previously obtained using MPACT and PANDAS-MOC. For the 2D problems, the percentage discrepancy from the other reference solutions was within 0.97%, whereas for the 3D problems, the discrepancy was within 0.97% and 1.57%, compared with the results obtained using MPACT and PANDAS-MOC, respectively. This verified the direct-SCM solver implemented in VITAS. In addition, the performance of the PCQ-SCM solver was compared with that of the direct-SCM solution for the TD0-5, TD3-4, and TD5-1 problems. Owing to the response matrix update strategy implemented in VITAS, the PCQ-SCM minimized the compute time for the step-perturbation problem TD0-5 by approximately twofold, as it decreased the number of iterations for larger time steps. For ramp-perturbation problems, the PCQ-SCM with larger time steps drastically lowered the number of response matrices and power iterations. The compute time for various problems, such as TD3-4, could be reduced by one to two orders of magnitude without sacrificing accuracy. Meanwhile, the analysis of the compute time revealed that the speedup efficiency of the PCQ-SCM was determined by the runtime fractions of matrix formation and equation solving. For the 2D ramp-perturbation problem, where the matrix formation typically dominates the runtime, the PCQ-SCM exhibited a high speedup efficiency, whereas for the 3D problems it exhibited a lower speedup efficiency, because equation solving dominated the runtime. Future efforts will focus on improving the parallel computation efficiency of the code, allowing to deploy it on computer clusters, which in turn will enable solving increasingly complex problems. Moreover, the high-fidelity transient multiphysics coupling of the neutron transport method with thermal-hydraulic simulations and uncertainty analysis [ 33 33 ] of transient multi-physics calculations will be performed. References G. Zhong , K. Xu , Y. Lu , et al. , Study on application of DAG-OPENMC in fusion neutronics analysis . Nuclear Techniques 45 ( 5 ), 050602 ( 2022 ). doi: 10.11889/j.0253-3219.2022.hjs.45.050602 http://doi.org/10.11889/j.0253-3219.2022.hjs.45.050602 ( in Chinese ) Study on application of DAG-OPENMC in fusion neutronics analysis Z. Zeng , S. Chen , Y. Zhang , et al. , Neutronics experiments of dual functional lithium-lead blanket based on D-T fusion neutron source . Nuclear Techniques 45 ( 4 ), 040601 ( 2022 ). doi: 10.11889/j.0253-3219.2022.hjs.45.040601 http://doi.org/10.11889/j.0253-3219.2022.hjs.45.040601 ( in Chinese ) Neutronics experiments of dual functional lithium-lead blanket based on D-T fusion neutron source X. Du , Y. Wang , Y. Zheng , et al. , The steady-state neutronic analysis and transient simulation of ADANES reactor design based on deterministic method . Nuclear Techniques 45 ( 10 ), 100601 ( 2022 ). doi: 10.11889/j.0253-3219.2022.hjs.45.100601 http://doi.org/10.11889/j.0253-3219.2022.hjs.45.100601 The steady-state neutronic analysis and transient simulation of ADANES reactor design based on deterministic method Q. Shen , Y. Wang , D. Jabaay , et al. , Transient analysis of C5G7-TD benchmark with MPACT . Ann. Nucl. Energy 125 , 107 - 120 ( 2019 ). doi: 10.1016/j.anucene.2018.10.049 http://doi.org/10.1016/j.anucene.2018.10.049 Transient analysis of C5G7-TD benchmark with MPACT A. Hsieh , G. Zhang , W. S. Yang , Consistent Transport Transient Solvers of the High-Fidelity Transport Code PROTEUS-MOC . Nucl. Sci. Eng. 194 ( 7 ), 508 - 540 ( 2020 ). doi: 10.1080/00295639.2020.1746619 http://doi.org/10.1080/00295639.2020.1746619 Consistent Transport Transient Solvers of the High-Fidelity Transport Code PROTEUS-MOC S. Tao , Y. Xu , Neutron transport analysis of C5G7-TD benchmark with PANDAS-MOC . Ann. Nucl. Energy. 169 , 108966 ( 2022 ). doi: 10.1016/j.anucene.2022.108966 http://doi.org/10.1016/j.anucene.2022.108966 Neutron transport analysis of C5G7-TD benchmark with PANDAS-MOC B. Wang , Z. Liu , J. Chen , et al. , A modified predictor-corrector quasi-static method in NECP-X for reactor transient analysis based on the 2D/1D transport method . Prog. Nucl. Energy. 108 , 122 - 135 ( 2018 ). doi: 10.1016/j.pnucene.2018.05.014 http://doi.org/10.1016/j.pnucene.2018.05.014 A modified predictor-corrector quasi-static method in NECP-X for reactor transient analysis based on the 2D/1D transport method Z. Shen , Q. Sun , D. He , et al. , Comparison and verification of NECP-X and OpenMC using high-fidelity BEAVRS benchmark models . Nuclear Techniques. 45 ( 01 ), 010602 ( 2022 ). doi: 10.11889/j.0253-3219.2022.hjs.45.010602 http://doi.org/10.11889/j.0253-3219.2022.hjs.45.010602 ( in Chinese ) Comparison and verification of NECP-X and OpenMC using high-fidelity BEAVRS benchmark models Y. Wang , S. Schunert , J. Ortensi , et al. , Rattlesnake: A MOOSE-based multiphysics multischeme radiation transport application . Nucl. Technol. 207 ( 7 ), 1047 - 1072 ( 2021 ). doi: 10.1080/00295450.2020.1843348 http://doi.org/10.1080/00295450.2020.1843348 Rattlesnake: A MOOSE-based multiphysics multischeme radiation transport application T. Zhang , Z. Li , Variational nodal methods for neutron transport: 40 years in review . Nucl. Eng. Technol. 54 ( 9 ), 3181 - 3204 ( 2022 ). doi: 10.1016/j.net.2022.04.012 http://doi.org/10.1016/j.net.2022.04.012 Variational nodal methods for neutron transport: 40 years in review T. Zhang , Y. Wang , E. E. Lewis , et al. , A three-dimensional variational nodal method for pin-resolved neutron transport analysis of pressurized water reactors . Nucl. Sci. Eng. 188 ( 2 ), 160 - 174 ( 2017 ). doi: 10.1080/00295639.2017.1350002 http://doi.org/10.1080/00295639.2017.1350002 A three-dimensional variational nodal method for pin-resolved neutron transport analysis of pressurized water reactors T. Zhang , E.E. Lewis , M.A. Smith , et al. , A variational nodal approach to 2D/1D pin resolved neutron transport for pressurized water reactors . Nucl. Sci. Eng. 186 ( 2 ), 120 – 133 ( 2017 ). doi: 10.1080/00295639.2016.1273023 http://doi.org/10.1080/00295639.2016.1273023 A variational nodal approach to 2D/1D pin resolved neutron transport for pressurized water reactors M.A. Smith , E.E. Lewis , E.R. Shemon , DIF3D-VARIANT 11.0: A Decade of Updates ( 2017 ). doi: 10.2172/1127298 http://doi.org/10.2172/1127298 DIF3D-VARIANT 11.0: A Decade of Updates Y. Wang , H. Wu , Y. Li , Comparison of two three-dimensional heterogeneous variational nodal methods for PWR control rod cusping effect and pin-by-pin calculation . Prog. Nucl. Energy. 101 , 370 – 380 ( 2017 ). doi: 10.1016/j.pnucene.2017.06.002 http://doi.org/10.1016/j.pnucene.2017.06.002 Comparison of two three-dimensional heterogeneous variational nodal methods for PWR control rod cusping effect and pin-by-pin calculation T. Zhang , W. Xiao , H. Yin et al. , VITAS: A multi-purpose simulation code for the solution of neutron transport problems based on variational nodal methods . Ann. Nucl. Energy 178 , 109335 ( 2022 ). doi: 10.1016/j.anucene.2022.109335 http://doi.org/10.1016/j.anucene.2022.109335 VITAS: A multi-purpose simulation code for the solution of neutron transport problems based on variational nodal methods F. He , X. Cai , W. Guo et al. , The transient analysis of molten salt reactor reactivity insertion based on RELAP5/FLUENT coupled program . Nuclear Techniques. 44 ( 3 ), 030601 ( 2021 ). doi: 10.11889/j.0253-3219.2021.hjs.44.030601 http://doi.org/10.11889/j.0253-3219.2021.hjs.44.030601 ( in Chinese ) The transient analysis of molten salt reactor reactivity insertion based on RELAP5/FLUENT coupled program S. Dulla , E.H. Mund , P. Ravetto , The quasi-static method revisited . Prog. Nucl. Energy. 50 ( 8 ), 908 – 920 ( 2008 ). doi: 10.1016/j.pnucene.2008.04.009 http://doi.org/10.1016/j.pnucene.2008.04.009 The quasi-static method revisited D. Caron , S. Dulla , P. Ravetto , New aspects in the implementation of the quasi-static method for the solution of neutron diffusion problems in the framework of a nodal method . Ann. Nucl. Energy. 87 , 34 - 48 ( 2016 ). doi: 10.1016/j.anucene.2015.02.035 http://doi.org/10.1016/j.anucene.2015.02.035 New aspects in the implementation of the quasi-static method for the solution of neutron diffusion problems in the framework of a nodal method Y.-A. Chao , A. Attard , A resolution of the stiffness problem of reactor kinetics . Nucl. Sci. Eng. 90 ( 1 ), 40 – 46 ( 1985 ). doi: 10.13182/NSE85-A17429 http://doi.org/10.13182/NSE85-A17429 A resolution of the stiffness problem of reactor kinetics S. Aoki , T. Suemura , J. Ogawa , et al. , The verification of 3 dimensional nodal kinetics code ANCK using transient benchmark problems . J. Nucl. Sci. Technol. 44 ( 6 ), 862 - 868 ( 2007 ). doi: 10.1080/18811248.2007.9711323 http://doi.org/10.1080/18811248.2007.9711323 The verification of 3 dimensional nodal kinetics code ANCK using transient benchmark problems B. W. Park , H. G. Joo , Improved stiffness confinement method within the coarse mesh finite difference framework for efficient spatial kinetics calculation . Ann. Nucl. Energy 76 , 200 – 208 ( 2015 ). doi: 10.1016/j.anucene.2014.09.029 http://doi.org/10.1016/j.anucene.2014.09.029 Improved stiffness confinement method within the coarse mesh finite difference framework for efficient spatial kinetics calculation C. Tang , Application of stiffness confinement method in transient transport calculation . J. Nucl. Eng. Radiat. Sc. 6 ( 1 )( 2019 ). doi: 10.1115/1.4044749 http://doi.org/10.1115/1.4044749 Application of stiffness confinement method in transient transport calculation W. Xiao , Q. Sun , X. Liu et al. , Application of stiffness confinement method within variational nodal method for solving time-dependent neutron transport equation . Comput. Phys. Commun. 279 , 108450 ( 2022 ). doi: 10.1016/j.cpc.2022.108450 http://doi.org/10.1016/j.cpc.2022.108450 Application of stiffness confinement method within variational nodal method for solving time-dependent neutron transport equation E.E. Lewis , W.F. Miller , Computational Methods of Neutron Transport . ( Wiley-Interscience , 1993 ) Computational Methods of Neutron Transport K. O. Ott , D. A. Meneley , Accuracy of the quasistatic treatment of spatial reactor kinetics . Nucl. Sci. Eng. 36 ( 3 ), 402 - 411 ( 1969 ). doi: 10.13182/NSE36-402 http://doi.org/10.13182/NSE36-402 Accuracy of the quasistatic treatment of spatial reactor kinetics M. Kheradmand Saadi , A. Abbaspour , Effective point kinetic parameters calculation in Tehran research reactor using deterministic and probabilistic methods . Nucl. Sci. Tech. 28 ( 12 ), 171 ( 2017 ). doi: 10.1007/s41365-017-0323-7 http://doi.org/10.1007/s41365-017-0323-7 Effective point kinetic parameters calculation in Tehran research reactor using deterministic and probabilistic methods K. F. Laurin-Kovitz , E. E. Lewis , Variational nodal transport perturbation theory . Nucl. Sci. Eng. 123 ( 3 ), 369 - 380 ( 1996 ). doi: 10.13182/NSE96-A24200 http://doi.org/10.13182/NSE96-A24200 Variational nodal transport perturbation theory A. Zhu , Y. Xu , T. Downar , A multilevel quasi-static kinetics method for pin-resolved transport transient reactor analysis . Nucl. Sci. Eng. 182 ( 4 ), 435 – 451 ( 2016 ). doi: 10.13182/NSE15-39 http://doi.org/10.13182/NSE15-39 A multilevel quasi-static kinetics method for pin-resolved transport transient reactor analysis F. Alcaro , S. Dulla , P. Ravetto , Implementation of the quasi-static method for neutron transport . International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2011) , Rio de Janeiro, RJ, Brazil , May 8-12, 2011 . Implementation of the quasi-static method for neutron transport W. M. Stacey , Nuclear reactor physics , 2 nd edn. ( WILEY-VCH ; Weinheim 2007 ), pp. 603 - 604 . Nuclear reactor physics J. Hou , K. Ivanov , V. Boyarinov , et al. , C5G7-TD Benchmark for Time-Dependent Heterogeneous Neutron Transport Calculations . International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2011) , Jeju, Korea , 2017 . C5G7-TD Benchmark for Time-Dependent Heterogeneous Neutron Transport Calculations J. Hou , K. N. Ivanov , V. F. Boyarinov , et al. , OECD/NEA benchmark for time-dependent neutron transport calculations without spatial homogenization . Nucl. Eng. Des. 317 , 177 - 189 ( 2017 ). doi: 10.1016/j.nucengdes.2017.02.008 http://doi.org/10.1016/j.nucengdes.2017.02.008 OECD/NEA benchmark for time-dependent neutron transport calculations without spatial homogenization C. Hao , J. Ma , X. Zhou , et al. , Development, verification and application of the uncertainty analysis platform CUSA . Ocean Engineering 261 , 112160 ( 2022 ). doi: 10.1016/j.oceaneng.2022.112160 http://doi.org/10.1016/j.oceaneng.2022.112160 Development, verification and application of the uncertainty analysis platform CUSA All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Wei Xiao. The first draft of the manuscript was written by Wei Xiao and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. We would like to thank Professor Yun-Lin Xu of Purdue University for providing the detailed C5G7-TD benchmark results data of the PANDAS-MOC. 10.1007/s41365-023-01170-x This research was supported by the National Natural Science Foundation of China (Nos. 12175138, U20B2011) and the Young Talent Project of the China National Nuclear Corporation. 2022-08-29 2022-12-06 2022-12-06 2023-02-20 2023-03-20 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2023 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-023-01170-x http://dx.doi.org/10.1007/s41365-023-01170-x . Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law. 2023 Nuclear Science and Techniques 2 34