An improved semi-implicit direct kinetics method for transient analysis of nuclear reactors

An improved semi-implicit direct kinetics method for transient analysis of nuclear reactors 1 corresp first-author Corresponding author, roozbehvadi@yahoo.com Corresponding author, roozbehvadi@yahoo.com roozbehvadi@yahoo.com 1 School of Nuclear Reactors and Safety, Nuclear Science and Technology Research Institute, End of North Karegar Ave., 14359-836, Tehran, Iran School of Nuclear Reactors and Safety, Nuclear Science and Technology Research Institute, End of North Karegar Ave., 14359-836, Tehran, Iran Semi-implicit direct kinetics (SIDK) is an innovative method for the temporal discretization of neutronic equations proposed by J. Banfield. The key approximation of the SIDK method is to substitute a time-averaged quantity for the fission source term in the delayed neutron differential equations. Hence, these equations are decoupled from prompt neutron equations and an explicit analytical representation of precursor groups is obtained, which leads to significant reduction in computational cost. As the fission source is not known in a time-step, the original study suggested using a constant quantity pertaining to the previous time-step for this purpose, and a reduction in the size of the time-step was proposed to lessen the imposed errors. However, this remedy notably diminishes the main advantage of the SIDK method. We discerned that if the original method is properly introduced into the algorithm of the point-implicit solver along with some modifications, the mentioned drawbacks will be mitigated adequately. To test this idea, a novel multi-group, multi-dimensional diffusion code using the finite-volume method and a point-implicit solver is developed which works in both transient and steady states. In addition to the SIDK, two other kinetic methods, i.e. direct kinetics and higher-order backward discretization (HOBD), are programmed into the diffusion code for comparison with the proposed model. The final code is tested at different conditions of two well-known transient benchmark problems. Results indicate that while the accuracy of the improved SIDK is comparable with the best available kinetic methods, it reduces the total time required for computation by up to 24%. Nuclear kinetics Semi-implicit direct kinetics Higher-order backward discretization Finite volume Point-implicit solver 1. Introduction 1 Insatiable growth in energy consumption has increased the demands on the rate and efficiency of nuclear power production systems [ 1 1 ] . It means higher discharge burnup, longer operation cycles, and more intense transient regimes, which lead to an increase in core loads on the fuel and core internals. In order to fulfill these demanding criteria, the technology incorporated into the design of these systems has to improve constantly, especially with regard to geometry optimization, and therefore the conventional one-dimensional numerical methodologies have been losing their function as reliable and accurate tools for both safety assessment and design purposes [ 2 2 ] . As a result, national and international programs such as Nuclear Energy Advanced Modeling and Simulation (NEAMS) [ 3 3 ] and CASL [ 4 4 ] in the US and Basic Research for Innovative Fuels Design for GEN IV systems (F-BRIDGE) in the European Union [ 5 5 ] , and similar projects in the other countries [ 6 6 - 8 8 6-8 ] , have been launched to develop an advanced generation of integrated modeling and safety codes. This study could be considered part of an independent academic effort that pursues the same objective [ 9 9 - 11 11 9-11 ] . The Fukushima Daiichi accident has particularly demonstrated the need to enhance transient modeling capabilities. Accordingly, the International Atomic Energy Agency (IAEA) launched a four-year Coordinated Research Project (CRP) named Fuel Modeling in Accident Conditions (FUMAC) [ 12 12 ] in 2014, and now, a new CRP in this regard is planned to be initiated in May 2019 [ 13 13 ] with a focus on accident-tolerant advanced-technology fuels and providing reliable experimental benchmarks for validation of the developed codes under extreme conditions. Transient analysis using these sorts of advanced codes has an intrinsically high computational cost, which increases with complexity and details of the target problem [ 1 1 ] . While the available numerical resources for running these codes have also improved continually, nevertheless, the necessity of optimizing the underlying numerical procedures is undeniable. One of the approaches to achieve this goal is to modify the reactor kinetic methods. SIDK is an innovative kinetic method developed by Banfield to reach the very same objective [ 14 14 ] . Banfield started the derivation of SIDK theory from factorization of the Boltzmann neutron transport equation, and then by neglecting the time-dependency of angular flux, which appears in the fission source term of the delayed neutron differential equations, this set of equations is decoupled from the set of prompt neutron equations and an explicit analytical solution to delayed neutron equations is achieved that is dependent upon a time-averaged fission power, while the prompt neutron equations are still solved in a fully implicit manner, similar to the original direct kinetic method [ 1 1 ] . However, decoupling the prompt and delayed equations can be traced back to an earlier study conducted by Ginestar et al. [ 15 15 ] in which, first, time-integration is applied to the transient form of the neutron diffusion equations, and then, by substituting a time-dependent extrapolation correlation for the fission source term in the delayed neutron equations, coefficients of the correlation are appeared as unknown parameters in the final analytical relation, calculated using the quantities of the fission sources obtained in the former time-steps. It is obvious that replacing the solution of up to six differential equations (pertaining to the number of delayed neutron precursor groups) with computation of the same number of simple analytical relations would result in a significant reduction in computational cost. However, the problem is that the fission source is neither known nor constant in a time-step, and the estimation techniques described above impose an error which drastically increases with the magnitude of the induced perturbation (as discussed and demonstrated in Sec. 2.3 and 4.2). A reduction in the size of the time-step was proposed to lessen the mentioned errors [ 1 1 ] . However, this remedy increases the computational cost, hence it notably diminishes the main advantage of the SIDK method. In this study, the key feature of the point-implicit solver [ 16 16 ] , which is the separated solution of the governing equations (as described in Sect. 2.2), is employed to mitigate the mentioned drawbacks. Accordingly, the SIDK’s analytical relations are inserted into the numerical iteration cycle of the point-implicit solver as independent explicit steps, while the quantity of fission source is not constant but is updated, according to a modified formulation, using the most recent values of neutron group fluxes obtained from solving the prompt neutron equations. This remedy should result in attaining a much more accurate estimation of the time-averaged fission source term in the current time-step, while retaining a significant portion of the computational cost reduction. In order to test this idea, a novel multi-group, multi-dimensional diffusion code using the finite-volume method (FVM) for spatial discretization and point-implicit solver is developed, which works in both transient and steady states. It is worth mentioning that although the FVM in its pure cell-based and geometric form (as used in this research) has many intrinsic advantages for solving neutronic equations, it has rarely been employed for this purpose, and hence, it could be considered as an innovation of the current study. Nonetheless, the specific matter of concern in this regard is that the FVM is highly compatible with the point-implicit solver [ 16 16 ] that is so crucial to this study, and it is extremely flexible regarding the geometric details of the solution domain and the boundary conditions applied to it [ 17 17 ] . The method of applying the FVM to the governing equations is delineated in Sect. 2.1, whereas an introduction to the point-implicit solver along with the resultant solution algorithms and a brief review of the secondary numerical methods are presented in Sect. 2.2. In addition to the SIDK method in both its original form and with the modifications applied in this study, two other kinetic methods, i.e. direct kinetics [ 18 18 ] and HOBD [ 15 15 ] , are programmed into the diffusion code in order to evaluate the proficiency of the proposed model in comparison with the best available and well-known kinetic methods. The descriptions in these regards are presented in Sect. 2.3. It should be noted that while the fully implicit solver has been employed for the original direct kinetic and HOBD methods [ 1 1 ] , the point-implicit solver is used for all kinetic methods developed in this research, in order to solely focus on the advantages provided by the modifications applied in this study to the SIDK method without involving the type of solver. However, note that application of the point-implicit solver could single-handedly reduce the computational cost by less than half compared to the fully implicit solver [ 19 19 ] . The final code is tested at different conditions of two well-known transient benchmark problems, i.e. TWIGL and LRA, which have been used in many studies for validation purposes. Therefore, the available results for comparison are abundant. An introduction to these benchmark problems is presented in Sec. 3, and the results and discussion are presented in Sect. 4. Accordingly, the main contributions of this study could be summarized by the following two points: • The original SIDK method leads to a notable reduction in the computational cost of transient analysis at the expense of an imposed inaccuracy, which escalates with intensity of the transient and details of the target model. Our improvements to the SIDK method result in resolving the accuracy issue while a significant portion of the computational cost reduction is retained (see Sect. 4). • In order to apply the improvements and test the new ideas, a novel multi-group, multi-dimensional diffusion code has been developed which works in both transient and steady states. The results of testing this code, in which the both improved and conventional transient methods have been embedded, prove the proficiency and accuracy of the novel code and the improved SIDK methods presented in this study. 2. Methodology 1 2.1. Spatial discretization of the governing equations 2 In this study, the FVM [ 17 17 ] has been employed for special discretization of the governing set of differential equations. For this purpose, the following three steps should be taken. 1. The solution domain should be divided into a proper number of control volumes (or numerical cells) using numerical grid generation applications. 2. All governing differential equations should be integrated about each control volume, yielding discrete equations which express and ensure the conservation laws on that control volume. Accordingly, the integral form of multi-group diffusion equations in the transient state with M groups of neutron precursors for an arbitrary control volume (as a part of the numerical grid) could be written as follows: <math display="block" id="M1"><mrow><mrow><mo>{</mo><mtable columnalign="left"><mtr><mtd><mfrac><mn>1</mn><mrow><msub><mi>v</mi><mi>g</mi></msub></mrow></mfrac><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ϕ</mi><mi>g</mi></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac></mrow></mrow></mstyle><mtext>d</mtext><mi>V</mi><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><msub><mi>A</mi><mtext>f</mtext></msub></mrow></msub><mrow><msub><mi>D</mi><mi>g</mi></msub><msub><mrow><mover accent="true"><mrow><mo>∇</mo><mi>ϕ</mi></mrow><mo stretchy="true">→</mo></mover></mrow><mi>g</mi></msub><mo>.</mo><mtext>d</mtext><mover accent="true"><mi>A</mi><mo>→</mo></mover></mrow></mrow></mstyle><mo>−</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><mo stretchy="false">(</mo><msub><mi>Σ</mi><mrow><mi>R</mi><mi>g</mi></mrow></msub><msub><mi>ϕ</mi><mi>g</mi></msub><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mtable columnalign="left"><mtr><mtd><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow></mtd></mtr><mtr><mtd><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>≠</mo><mi>g</mi></mrow></mtd></mtr></mtable><mi>G</mi></munderover><mrow><msub><mi>Σ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow></msub><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub></mrow></mstyle><mo>−</mo><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><msub><mi>χ</mi><mi>g</mi></msub><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mi>ν</mi><msub><mi>Σ</mi><mrow><mi>f</mi><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub></mrow></mstyle><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>m</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover><mrow><msub><mi>λ</mi><mi>m</mi></msub><msub><mi>C</mi><mi>m</mi></msub></mrow></mstyle><mo stretchy="false">)</mo><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle></mtd></mtr><mtr><mtd><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><mfrac><mrow><mo>∂</mo><msub><mi>C</mi><mi>m</mi></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><mo stretchy="false">(</mo><msub><mi>β</mi><mi>m</mi></msub><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mi>ν</mi><msub><mi>Σ</mi><mrow><mi>f</mi><mi>g</mi></mrow></msub><msub><mi>ϕ</mi><mi>g</mi></msub></mrow></mstyle><mo>−</mo><msub><mi>λ</mi><mi>m</mi></msub><msub><mi>C</mi><mi>m</mi></msub></mrow></mrow></mstyle><mo stretchy="false">)</mo><mtext>d</mtext><mi>V</mi><mo>;</mo><mtext>     </mtext><mi>g</mi><mtext>=1,2,</mtext><mn>...</mn><mtext>,</mtext><mi>G</mi><mtext>     ,     </mtext><mi>m</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>M</mi></mtd></mtr></mtable></mrow></mrow></math> in which ϕ , C , G , and V are neutron flux, concentration of delayed neutron precursor, number of neutron energy groups, and volume of the arbitrary control volume, respectively. The two-dimensional cell ( c 0 ) shown in Fig. 1 Fig. 1 is an example of such a control volume on which the geometric notations used in this section have been illustrated. Accordingly, the spatial discretization of Eq. (1) on the given cell ( c 0 ) yields: <math display="block" id="M2"><mrow><mrow><mo>{</mo><mtable columnalign="left"><mtr><mtd><mfrac><mn>1</mn><mrow><msub><mi>v</mi><mi>g</mi></msub></mrow></mfrac><mfrac><mrow><mo>∂</mo><msub><mi>ϕ</mi><mi>g</mi></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mi>V</mi><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>f</mi><mrow><msub><mi>N</mi><mtext>f</mtext></msub></mrow></munderover><mrow><msub><mi>D</mi><mi>g</mi></msub><msubsup><mrow><mover accent="true"><mrow><mo>∇</mo><mi>ϕ</mi></mrow><mo stretchy="true">→</mo></mover></mrow><mi>g</mi><mi>f</mi></msubsup><mo>.</mo><msub><mover accent="true"><mi>A</mi><mo>→</mo></mover><mtext>f</mtext></msub></mrow></mstyle><mo>−</mo><mo stretchy="false">(</mo><msub><mi>Σ</mi><mrow><mi>R</mi><mi>g</mi></mrow></msub><msub><mi>ϕ</mi><mi>g</mi></msub><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mtable columnalign="left"><mtr><mtd><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow></mtd></mtr><mtr><mtd><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>≠</mo><mi>g</mi></mrow></mtd></mtr></mtable><mi>G</mi></munderover><mrow><msub><mi>Σ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow></msub><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub></mrow></mstyle><mo>−</mo><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><msub><mi>χ</mi><mi>g</mi></msub><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mi>ν</mi><msub><mi>Σ</mi><mrow><mi>f</mi><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub></mrow></mstyle><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>m</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover><mrow><msub><mi>λ</mi><mi>m</mi></msub><msub><mi>C</mi><mi>m</mi></msub></mrow></mstyle><mo stretchy="false">)</mo><mi>V</mi></mtd></mtr><mtr><mtd><mfrac><mrow><mo>∂</mo><msub><mi>C</mi><mi>m</mi></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mi>V</mi><mo>=</mo><mo stretchy="false">(</mo><msub><mi>β</mi><mi>m</mi></msub><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mi>ν</mi><msub><mi>Σ</mi><mrow><mi>f</mi><mi>g</mi></mrow></msub><msub><mi>ϕ</mi><mi>g</mi></msub></mrow></mstyle><mo>−</mo><msub><mi>λ</mi><mi>m</mi></msub><msub><mi>C</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mi>V</mi><mo>;</mo><mtext>     </mtext><mi>g</mi><mtext> =1,2,</mtext><mn>...</mn><mtext>,</mtext><mi>G</mi><mtext>,     </mtext><mi>m</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>M</mi></mtd></mtr></mtable></mrow></mrow></math> where A f is surface area vector and ∇ ϕ is gradient of neutron flux (see Fig. 1 Fig. 1 ). Discretization of temporal terms will be discussed in Sect. 2.3. 3. In order to linearize and solve Eq. (2), along with similar equations associated with the other cells which constitute the solution domain as a set, the gradient of ϕ should be evaluated at the faces enclosing each cell. In this study, an improved form of the least-square cell-based scheme that we introduced in Ref. [ 9 9 ] is employed for this purpose, which is essentially based on second-order linear variation of ϕ at cell centers and can be stated as: <math display="block" id="M3"><mrow><mo stretchy="false">[</mo><mi>J</mi><mo stretchy="false">]</mo><msubsup><mrow><mrow><mo>(</mo><mrow><mo>∇</mo><mi>ϕ</mi></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>c</mi><mtext>o</mtext></msub></mrow><mtext>f</mtext></msubsup><mo>=</mo><mi>Δ</mi><mi>ϕ</mi><mo>,</mo></mrow></math> in which [ J ] is a geometric-based coefficient matrix and ∆ ϕ is a matrix with components of ϕ ci - ϕ co , where " i " is the counter of the faces (or the cells) enclosing each target cell [ 9 9 ] . In steady-state conditions, temporal derivatives are eliminated from the left side of Eq. (1) and the remaining equations form a set of eigenvalue differential equations with a trivial solution of all zeros for spatial distribution of neutron group flux. In order to find the target non-trivial solution with positive quantities for neutron flux, the equation-based power iteration method [ 19 19 ] is employed in this study, in which fission cross-sections are divided by the effective multiplication factor (denoted by k eff ), which is the first and largest eigenvalue of neuron diffusion equation in steady-state conditions: <math display="block" id="M4"><mrow><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><msub><mi>A</mi><mtext>f</mtext></msub></mrow></msub><mrow><msub><mi>D</mi><mi>g</mi></msub><msub><mrow><mover accent="true"><mrow><mo>∇</mo><mi>ϕ</mi></mrow><mo stretchy="true">→</mo></mover></mrow><mi>g</mi></msub><mo>.</mo><mtext>d</mtext><mover accent="true"><mi>A</mi><mo>→</mo></mover></mrow></mrow></mstyle><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><mo stretchy="false">(</mo><msub><mi>Σ</mi><mrow><mi>R</mi><mi>g</mi></mrow></msub><msub><mi>ϕ</mi><mi>g</mi></msub><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mtable columnalign="left"><mtr><mtd><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow></mtd></mtr><mtr><mtd><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>≠</mo><mi>g</mi></mrow></mtd></mtr></mtable><mi>G</mi></munderover><mrow><msub><mi>Σ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow></msub><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub></mrow></mstyle><mo>−</mo><msub><mi>χ</mi><mi>g</mi></msub><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mi>ν</mi><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>f</mi><msup><mi>g</mi><mo>′</mo></msup></mrow></msub></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub></mrow></mstyle><mo stretchy="false">)</mo><mtext>d</mtext><mi>V</mi></mrow></mrow></mstyle><mo>.</mo></mrow></math> 2.2. Solution algorithm and numerical methods 2 After applying the discretization and linearization procedures described in the preceding section to the governing equations at each numerical cell, a set of highly coupled algebraic equations is obtained for the whole solution domain that should be solved. The point-implicit solver is employed for this purpose, in which the equations are solved separately and sequentially, according to the algorithm represented in Fig. 2(a) Fig. 2(a) . In other words, each equation is solved separately in an implicit manner for the specific variable associated with that equation (for example, the diffusion equation pertaining to the first group of neutron energy is solved for ϕ 1 ), while all other variables that appeared in that equation are replaced with the most updated quantities available for those variables. Then, the new quantities obtained for the target variable ( ϕ 1 ) are substituted for the older ones and used in the solution to the next equation. This process continues until a proper level of convergence is achieved for all variables, including the group neutron flux, according to the following relations: <math display="block" id="M5"><mrow><mfrac><mrow><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><msubsup><mi>ϕ</mi><mi>g</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup></mrow></mrow></mstyle><mtext>d</mtext><mi>V</mi><mo>−</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><msubsup><mi>ϕ</mi><mi>g</mi><mi>n</mi></msubsup></mrow></mrow></mstyle><mtext>d</mtext><mi>V</mi></mrow><mrow><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>V</mi></msub><mrow><msubsup><mi>ϕ</mi><mi>g</mi><mi>n</mi></msubsup></mrow></mrow></mstyle><mtext>d</mtext><mi>V</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>C</mi><mi>e</mi></mrow></munderover><mrow><msubsup><mi>ϕ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><msub><mi>V</mi><mi>i</mi></msub></mrow></mstyle><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>C</mi><mi>e</mi></mrow></munderover><mrow><msubsup><mi>ϕ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mi>n</mi></msubsup><msub><mi>V</mi><mi>i</mi></msub></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>C</mi><mi>e</mi></mrow></munderover><mrow><msubsup><mi>ϕ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow><mi>n</mi></msubsup><msub><mi>V</mi><mi>i</mi></msub></mrow></mstyle></mrow></mfrac><mo><</mo><msub><mi>ε</mi><mi>g</mi></msub><mtext>         </mtext><mi>g</mi><mo>=</mo><mn>1</mn><mo>,</mo><mtext> 2, </mtext><mn>.....</mn><mtext>, </mtext><mi>G</mi></mrow></math> where n is the counter of the solver’s iterations and Ce is the total number of numerical cells which constitute the solution domain. A very similar relation is applied for convergence of the delayed neutron precursor equations as well. Afterwards, the calculations pertaining to the next time-step start and the whole procedure is repeated until the end of the time-steps. In steady-state conditions, the set of equations presented in Eq. 4 are solved with the same solution approach but based on the equation-based power iteration method [ 20 20 ] , according to the algorithm shown in Fig. 2(b) Fig. 2(b) . Table 1 Table 1 presents a brief survey of the numerical models and convergence criterion employed in this study, along with the references in which further information regarding each method can be found (Ref. [ 14 14 - 18 18 14-18 , 20 20 - 24 24 20-24 ] ). 2 1 Title 1 1 Applied method 2 1 Spatial discretization 1 1 Finite volume [ 17 17 ] 2 1 Temporal discretization 1 1 Original SIDK [ 14 14 ] 2 1 1 1 Improved SIDK of first and second order 2 1 1 1 Direct kinetics [ 18 18 ] 2 1 1 1 HOBD [ 15 15 ] 2 1 Eigenvalue evaluation (steady-state solution) 1 1 Equation-based power iteration method [ 20 20 ] 2 1 Solver approach 1 1 Point-implicit solver [ 16 16 ] 2 1 Linearization shape function 1 1 Second-order linear function [ 17 17 ] 2 1 Solution method of linear algebraic equations (LAEs) set 1 1 LU decomposition with partial implicit pivoting [ 21 21 ] 2 1 Gradients evaluation method 1 1 Least-squares cell-based [ 22 22 ] 2 1 Gradient (slope) limiter 1 1 Differential type [ 23 23 ] 2 1 Under-relaxation method 1 1 Variables (explicit) relaxation scheme [ 17 17 ] 2 1 Residual convergence criterion 1 1 1 1 Neutron flux 1 1 10 -6 (for all groups) 1 1 1 1 Precursor concentration 1 1 10 -6 (for all groups) 2 1 Solution accelerator 1 1 Algebraic Multi-Grid (AMG) with V-cycle recursive procedure [ 24 24 ] 2.3. Temporal discretization and the modifications applied to improve the SIDK method 2 According to the HOBD (higher-order backward discretization) method [ 15 15 , 25 25 ] , which belongs to a category of differentiation methods known as the multi-step method in numerical calculation science [ 26 26 ] , the discretized form of the temporal term for group flux equations can be written as follows [ 1 1 , 15 15 ] : <math display="block" id="M6"><mrow><msubsup><mi>ϕ</mi><mi>g</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>+</mo><msub><mi>a</mi><mn>1</mn></msub><msubsup><mi>ϕ</mi><mi>g</mi><mi>k</mi></msubsup><mo>+</mo><msub><mi>a</mi><mn>2</mn></msub><msubsup><mi>ϕ</mi><mi>g</mi><mrow><mi>k</mi><mo>−</mo><mn>1</mn></mrow></msubsup><mo>+</mo><mo>⋯</mo><mo>+</mo><msub><mi>a</mi><mi>m</mi></msub><msubsup><mi>ϕ</mi><mi>g</mi><mrow><mi>k</mi><mo>−</mo><mi>m</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>=</mo><mi>Δ</mi><mi>t</mi><mo>×</mo><mi>b</mi><mo>×</mo><msubsup><mi>F</mi><mi>g</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo stretchy="false">(</mo><mi>ϕ</mi><mo>,</mo><mi>C</mi><mo stretchy="false">)</mo></mrow></math> in which " k " is the counter for time-step and " m " is the order of the HOBD method that denotes the number of backward time-steps at which quantities of all solution variables should be stored and incorporated into the calculations associated with the current time-step. " a " and " b " are constants calculated for each specific order ( m ) based on the simple method explained in Ref. [ 15 15 ] . The HOBD of first order ( m  = 1) leads to the well-known direct kinetics method [ 18 18 ] , which has been frequently employed for transient analysis of nuclear systems [ 1 1 ] . <math display="block" id="M7"><mrow><mfrac><mrow><msubsup><mi>ϕ</mi><mi>g</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mi>ϕ</mi><mi>g</mi><mi>k</mi></msubsup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mo>=</mo><msubsup><mi>F</mi><mi>g</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo stretchy="false">(</mo><mi>ϕ</mi><mo>,</mo><mi>C</mi><mo stretchy="false">)</mo></mrow></math> The second-order HOBD method ( m  = 2) takes the following form: <math display="block" id="M8"><mrow><mfrac><mrow><mn>3</mn><msubsup><mi>ϕ</mi><mi>g</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><mn>4</mn><msubsup><mi>ϕ</mi><mi>g</mi><mi>k</mi></msubsup><mo>+</mo><msubsup><mi>ϕ</mi><mi>g</mi><mrow><mi>k</mi><mo>−</mo><mn>1</mn></mrow></msubsup></mrow><mrow><mn>2</mn><mi>Δ</mi><mi>t</mi></mrow></mfrac><mo>=</mo><msubsup><mi>F</mi><mi>g</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo stretchy="false">(</mo><mi>ϕ</mi><mo>,</mo><mi>C</mi><mo stretchy="false">)</mo><mo>.</mo></mrow></math> It is obvious that the accuracy of the solution increases with " m ". However, at the same time, it would also lead to exponential escalation in computational cost due to the necessity of storing all variables in the whole calculation domain for every backward time-step involved in the solution. As a result, there should be a limit on the growth in " m " that it is decided to be the second order (i.e. Eq. (8)) in this study. Please note that while the HOBD method requires the quantities of variables at two former time-steps, only one set of quantities, i.e. the initial values, is available at the first time-step. In that case, Eq. (7) (the lower-order method) is used in the first time-step and application of Eq. (8) is started from the second time-step. As for the SIDK method, if the neutron fission source term in the delayed neutron precursors equations is assumed to be a constant time-averaged quantity (denoted by " S ") at each time-step, the second set of equations presented in Eq. 1 will be modified as follows: <math display="block" id="M9"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>C</mi><mi>m</mi></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><msub><mi>β</mi><mi>m</mi></msub><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mi>ν</mi><msub><mi>Σ</mi><mrow><mi>f</mi><mi>g</mi></mrow></msub><msub><mi>ϕ</mi><mi>g</mi></msub></mrow></mstyle><mo>−</mo><msub><mi>λ</mi><mi>m</mi></msub><msub><mi>C</mi><mi>m</mi></msub><mo>=</mo><msub><mi>β</mi><mi>m</mi></msub><mi>S</mi><mo>−</mo><msub><mi>λ</mi><mi>m</mi></msub><msub><mi>C</mi><mi>m</mi></msub><mo>;</mo><mtext>     </mtext><mi>m</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>M</mi></mrow></math> Now, if the integral factor, <math id="M10"><mrow><msup><mi>e</mi><mrow><mo>−</mo><msub><mi>λ</mi><mi>m</mi></msub><mi>Δ</mi><mi>t</mi></mrow></msup></mrow></math> is multiplied by both sides of Eq. 9 and then both sides are integrated about 0 to t , after simple reordering, the following analytical relations will be obtained. <math display="block" id="M11"><mrow><msubsup><mi>C</mi><mi>m</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>=</mo><msubsup><mi>C</mi><mi>m</mi><mi>k</mi></msubsup><mo>×</mo><msup><mi>e</mi><mrow><msub><mi>λ</mi><mi>m</mi></msub><mi>Δ</mi><mi>t</mi></mrow></msup><mo>+</mo><mfrac><mrow><msub><mi>β</mi><mi>m</mi></msub></mrow><mrow><msub><mi>λ</mi><mi>m</mi></msub></mrow></mfrac><mo>×</mo><mi>S</mi><mo>×</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msup><mi>e</mi><mrow><msub><mi>λ</mi><mi>m</mi></msub><mi>Δ</mi><mi>t</mi></mrow></msup></mrow><mo>)</mo></mrow><mo>;</mo><mtext>     </mtext><mi>m</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>M</mi></mrow></math> Accordingly, the solution of M differential equations of delayed neutron precursors at each time-step is replaced by a simple calculation of the above analytical relations, whereby computational cost is expected to reduce notably. However, this method has two drawbacks. The fission source depends on the summation of group fluxes and this is unknown at each time-step, so it is suggested in the original study [ 14 14 ] to use the quantities associated with the previous time-step instead, as constant values. Furthermore, fission source term changes during the time-step and replacing constant values would impose deviation. A reduction in the size of time-steps was suggested in the original research [ 12 12 ] as a solution to the mentioned issues. However, this remedy increases the computational cost and, hence, it notably reduces the main advantage of the original SIDK method. In this study, the first idea for dealing with the mentioned drawbacks was to use the quantities of fission sources calculated in the previous time-steps in order to extrapolate the fission source at the current time-step, similar to the methodology adopted by Ginestar et al. [ 15 15 ] . However, this approach also has two problems. • Even for a linear extrapolation which is not so accurate, the neutron flux associated with all groups of neutron energy should be stored for at least two previous time-steps for the whole simulation domain. Naturally, if higher accuracy is required, much more data needs to be stored. The imposed computational cost reduces the advantage of the SIDK method. • Transient analysis of large and fast temporal perturbations (as we will see regarding the LRA problem) usually encounter sharp peaks, especially for the quantities of power and neutron flux. Accurate estimation of these peaks and the time at which they occur are of the utmost importance. It is obvious that an extrapolation technique based on the quantities at previous time-steps could not accurately estimate such abrupt variations and would cause significant deviation. Then, we realized that if the analytical relations (presented in Eq. (10)) are inserted into the numerical iteration cycle of the point-implicit solver instead of the equivalent steps pertaining to neutron precursor equations (see Fig. 2 (a) Fig. 2 (a) ), while the quantity of " S " is being updated using the most recent values of neutron group fluxes obtained from solving the prompt neutron equations at each iteration, a much more accurate estimation of the fission source term at the current time-step will be incorporated into the calculation of the precursor concentrations. However, direct replacement of the analytical relations would mean that the quantity of fission source at the end of each time-step should be used, whereas the correct quantity of fission source is the average of this term during the time-step. Accordingly, the final analytical relations employed in this study have been modified as follows: <math display="block" id="M12"><mrow><msubsup><mi>C</mi><mi>m</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>=</mo><msubsup><mi>C</mi><mi>m</mi><mi>k</mi></msubsup><mo>×</mo><msup><mi>e</mi><mrow><msub><mi>λ</mi><mi>m</mi></msub><mi>Δ</mi><mi>t</mi></mrow></msup><mo>+</mo><mfrac><mrow><msub><mi>β</mi><mi>m</mi></msub></mrow><mrow><msub><mi>λ</mi><mi>m</mi></msub></mrow></mfrac><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><mi>ν</mi><msub><mi>Σ</mi><mrow><mi>f</mi><mi>g</mi></mrow></msub><msubsup><mi>ϕ</mi><mi>g</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup></mrow></mstyle><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mi>ν</mi><msub><mi>Σ</mi><mrow><mi>f</mi><mi>g</mi></mrow></msub><msubsup><mi>ϕ</mi><mi>g</mi><mi>k</mi></msubsup></mrow></mstyle></mrow><mn>2</mn></mfrac><mo>×</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msup><mi>e</mi><mrow><msub><mi>λ</mi><mi>m</mi></msub><mi>Δ</mi><mi>t</mi></mrow></msup></mrow><mo>)</mo></mrow><mo>;</mo><mtext>     </mtext><mi>m</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>M</mi></mrow></math> These relations state that the quantity of fission source term at each time-step is the arithmetic average of the respective initial and final values. The last point regarding the SIDK approach is selection of the method that should be used for temporal discretization of the prompt neutron equations (the first set in Eq. (1)) in combination with Eq. (11). The direct kinetics method (Eq. (7)) was used in the original research [ 14 14 ] for this purpose; however, the HOBD method (Eq. (8)) could also be employed. Therefore, two orders of the SIDK method are defined in this study based on whether the direct kinetics method (first order) or HOBD method (second order) is used for temporal discretization of the prompt neutron equations. 3. Numerical solution setup 1 3.1. TWIGL problem 2 The TWIGL benchmark problem is a 2D model of a seed-blanket reactor with two neutron energy groups and one delayed neutron precursor group [ 27 27 ] . The geometry of the solution domain is depicted in Fig. 3 Fig. 3 along with dimensions, boundary conditions and region assignments. Note that due to symmetry, a quadrant of the reactor is modeled. A uniform quadrilateral mesh with a size of 1 cm is applied to the geometry. The group constants for the steady-state calculations are presented in Table 2 Table 2 . Delayed neutron data and other constants required for solving the TWIGL problem are tabulated in Table 3 Table 3 . The transient state consists of two separate problems, which are both models of subcritical reactivity insertion accidents initiated from the steady state. In both problems, the absorption cross-sections of the thermal (second) group are increased but in two different forms of step and ramp perturbations, defined respectively as follows [ 27 27 ] : 1 1 Region 1 1 Group ( g ) 1 1 D g [cm] 1 1 <math id="M13"><mrow><msubsup><mi>Σ</mi><mi>g</mi><mi>a</mi></msubsup></mrow></math> [cm -1 ] 1 1 <math id="M14"><mrow><mi>ν</mi><msubsup><mi>Σ</mi><mi>g</mi><mi>f</mi></msubsup></mrow></math> [cm -1 ] 1 1 Σ 12 [cm -1 ] 1 1 1 1 1 1 1 1 1.4 1 1 0.01 1 1 0.007 1 1 0.01 1 1 1 1 2 1 1 0.4 1 1 0.15 1 1 0.2 1 1 ------ 1 1 2 1 1 1 1 1 1.4 1 1 0.01 1 1 0.007 1 1 0.01 1 1 1 1 2 1 1 0.4 1 1 0.15 1 1 0.5 1 1 ------ 1 1 3 1 1 1 1 1 1.3 1 1 0.008 1 1 0.003 1 1 0.01 1 1 1 1 2 1 1 0.5 1 1 0.05 1 1 0.06 1 1 ------ 1 1 Parameter 1 1 β 1 1 λ [1/s] 1 1 ν 1 1 ν 1 [m/s] 1 1 ν 2 [m/s] 1 1 χ 1 1 1 χ 2 1 1 Value 1 1 0.0075 1 1 0.08 1 1 2.43 1 1 1.0×10 5 1 1 2.0×10 3 1 1 1.0 1 1 0.0 <math display="block" id="M15"><mrow><mi>Δ</mi><msub><mi>Σ</mi><mrow><mi>a</mi><mn>2</mn></mrow></msub><mo>=</mo><mo>−</mo><mn>0.0035</mn><msup><mrow><mtext>  [cm</mtext></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup><mtext>]          </mtext><mi>t</mi><mo>=</mo><mn>0.0</mn></mrow></math> <math display="block" id="M16"><mrow><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>a</mi><mn>2</mn></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>Σ</mi><mrow><mi>a</mi><mn>2</mn></mrow></msub><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo></mrow></mfrac><mo>=</mo><mrow><mo>{</mo><mtable columnalign="left"><mtr><mtd><mn>1</mn><mo>−</mo><mn>0.11667</mn><mo>×</mo><mi>t</mi><mtext>        </mtext><mi>t</mi><mo>≤</mo><mn>0.2</mn></mtd></mtr><mtr><mtd><mn>0.97666</mn><mtext>                 </mtext><mi>t</mi><mo>></mo><mn>0.2</mn><mtext>   </mtext></mtd></mtr></mtable></mrow></mrow></math> The total time of computation is 0.5 s [ 27 27 ] . The time-step is selected to be 0.001 s. According to Ref. [ 27 27 ], the initial average power density of the core (based on fission production) should be assumed as follows, and the quantities calculated at the following time-steps for this parameter should be normalized and presented accordingly. <math display="block" id="M17"><mrow><mover accent="true"><mi>P</mi><mo>¯</mo></mover><mo stretchy="false">(</mo><mi>t</mi><mo>=</mo><mn>0</mn><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>V</mi><mrow><mtext>core</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><msub><mi>V</mi><mrow><mtext>core</mtext></mrow></msub></mrow></msub><mrow><mrow><mo>(</mo><mrow><msub><mi>Σ</mi><mrow><mi>f</mi><mn>1</mn></mrow></msub><msub><mi>ϕ</mi><mn>1</mn></msub><mo>+</mo><msub><mi>Σ</mi><mrow><mi>f</mi><mn>2</mn></mrow></msub><msub><mi>ϕ</mi><mn>2</mn></msub></mrow><mo>)</mo></mrow></mrow></mrow></mstyle><mtext>d</mtext><mi>V</mi><mo>=</mo><mn>1.0</mn><mtext> </mtext><mrow><mo>[</mo><mrow><mfrac><mo>#</mo><mrow><mtext>sec</mtext><mo>⋅</mo><msup><mrow><mtext>cm</mtext></mrow><mn>3</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow></mrow></math> 3.2. LRA problem 2 The LRA benchmark problem is a simplified multi-dimensional model of a boiling water reactor (BWR) with two neutron energy groups, two delayed neutron precursor groups, an adiabatic heat-up, and Doppler thermal feedback models [ 28 28 ] . A schematic view of the reactor horizontal cross-section is depicted in Fig. 4 Fig. 4 , in which dimensions, boundary conditions and region assignments are denoted. Note that due to symmetry, a quadrant of the geometry is modeled. A uniform hexahedral mesh with a size of 1 cm is applied to the solution domain. Initial two-group constants and delayed neutron data are tabulated in Tables 4 Tables 4 and 5 5 , respectively. Note that according to the instructions from Ref. [ 28 28 ] , k eff should be calculated and reported for two separate steady-state configurations of the LRA problem’s reactor core based on the position of the control-rod bundle. The first configuration is the case in which the control-rod bundle is fully inserted into the core (or the "initial" critical state), and the second one is the case in which the bundle is fully withdrawn from the core (or the "final" state) for which the thermal feedback effects should be neglected or, in other words, it should be assumed that the reactor is in the cold shut-down state. Since the control-rod bundle is located in region 6 (denoted in Fig. 4 Fig. 4 ), the group constants of this region are affected by repositioning of the control elements, as indicated by the words "Initial" and "Final" in Table 4 Table 4 . 1 1 Region 1 1 Group ( g ) 1 1 <math id="M18"><mrow><msubsup><mi>Σ</mi><mi>g</mi><mi>a</mi></msubsup></mrow></math> [cm -1 ] 1 1 <math id="M19"><mrow><mi>ν</mi><msubsup><mi>Σ</mi><mi>g</mi><mi>f</mi></msubsup></mrow></math> [cm -1 ] 1 1 Σ 12 [cm -1 ] 1 1 D g [cm] 1 1 1 1 1 1 1 1 0.008252 1 1 0.004602 1 1 0.02533 1 1 1.255 1 1 1 1 2 1 1 0.1003 1 1 0.1091 1 1 ------ 1 1 0.211 1 1 2 1 1 1 1 1 0.007181 1 1 0.004609 1 1 0.02767 1 1 1.268 1 1 1 1 2 1 1 0.07047 1 1 0.08675 1 1 ------ 1 1 0.1902 1 1 3 1 1 1 1 1 0.008002 1 1 0.004643 1 1 0.02617 1 1 1.259 1 1 1 1 2 1 1 0.08344 1 1 0.1021 1 1 ------ 1 1 0.2091 1 1 4 1 1 1 1 1 0.008002 1 1 0.004643 1 1 0.02617 1 1 1.259 1 1 1 1 2 1 1 0.073324 1 1 0.1021 1 1 ------ 1 1 0.2091 1 1 5 1 1 1 1 1 0.0006034 1 1 0.0 1 1 0.04754 1 1 1.257 1 1 1 1 2 1 1 0.01911 1 1 0.0 1 1 ------ 1 1 0.1592 1 1 6 1 1 1 1 1 0.008002 1 1 0.004643 1 1 0.02617 1 1 1.259 1 1 1 1 2 1 1 1 1 Initial 1 1 0.08344 1 1 0.1021 1 1 ------ 1 1 0.2091 1 1 1 1 Final 1 1 0.07332 1 1 1 1 1 1 1 1 Parameter 1 1 β 1 1 1 β 2 1 1 λ 1 [1/s] 1 1 λ 2 [1/s] 1 1 α [K.cm 3 ] 1 1 γ [K 1/2 ] 1 1 ν 1 1 ν 1 [m/s] 1 1 ν 2 [m/s] 1 1 χ 1 1 1 χ 2 1 1 Value 1 1 0.0054 1 1 0.001087 1 1 0.00654 1 1 1.35 1 1 3.83×10 -11 1 1 2.043×10 -3 1 1 2.43 1 1 3.0×10 5 1 1 3.0×10 3 1 1 1.0 1 1 0.0 The transient state associated with the LRA problem is a model of a fast, super-prompt critical reactivity insertion accident initiated from the initial critical configuration by rapid simultaneous withdrawal of the four control elements form the reactor core. This perturbation is modeled by a reduction in the absorption cross-section of the second neutron energy group in region 6 (denoted in Fig. 4 Fig. 4 ), according to the following relation: <math display="block" id="M20"><mrow><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>a</mi><mn>2</mn></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>Σ</mi><mrow><mi>a</mi><mn>2</mn></mrow></msub><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo></mrow></mfrac><mo>=</mo><mrow><mo>{</mo><mtable columnalign="left"><mtr><mtd><mn>1</mn><mo>−</mo><mn>0.0606184</mn><mo>×</mo><mi>t</mi><mtext>     </mtext><mi>t</mi><mo><</mo><mn>2.0</mn></mtd></mtr><mtr><mtd><mn>0.8787631</mn><mtext>               </mtext><mi>t</mi><mo>≥</mo><mn>2.0</mn><mtext>   </mtext></mtd></mtr></mtable></mrow></mrow></math> The time-step and total time for the modeled transient are 0.001 and 3 s, respectively. The adiabatic heat-up and Doppler thermal feedback models are stated by the following equations, respectively. <math display="block" id="M21"><mrow><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mi>α</mi><mrow><mo>(</mo><mrow><msub><mi>Σ</mi><mrow><mi>f</mi><mn>1</mn></mrow></msub><msub><mi>ϕ</mi><mn>1</mn></msub><mo>+</mo><msub><mi>Σ</mi><mrow><mi>f</mi><mn>2</mn></mrow></msub><msub><mi>ϕ</mi><mn>2</mn></msub></mrow><mo>)</mo></mrow></mrow></math> <math display="block" id="M22"><mrow><msub><mi>Σ</mi><mrow><mi>a</mi><mn>1</mn></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><msub><mi>Σ</mi><mrow><mi>a</mi><mn>1</mn></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo>=</mo><mn>0</mn><mo stretchy="false">)</mo><mo>×</mo><mrow><mo>[</mo><mrow><mn>1</mn><mo>+</mo><mi>γ</mi><mrow><mo>(</mo><mrow><mi>T</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mi>T</mi><mi>o</mi></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow></mrow></math> The constants mentioned in the above equations have been given in Table 5 Table 5 along with the delayed neutron data. 4. Results and discussion 1 4.1. TWIGL problem 2 The effective multiplication factor ( k eff ) calculated for the steady-state conditions of the TWIGL problem is presented in Table 6 Table 6 , along with analogous data from various studies [ 1 1 , 29 29 , 30 30 ] . The results indicate less than 3 pcm (nearly 0.003%) relative difference between k eff calculated in the current study and the most accurate data presented in Ref. [ 29 29 , 30 30 ] , in which 2D and 3D nodal analytical methods were employed, respectively. 1 1 Ref. 1 1 Current study 1 1 CONQUEST [29] 1 1 Banfield [1] 1 1 QUANDRY [30] 1 1 k eff 1 1 0.913176 1 1 0.91320 1 1 0.91385 1 1 0.91321 As for the transient state, variations in the relative power density of the core versus time for the step and ramp perturbations are demonstrated in Fig. 5 Fig. 5 , along with the most accurate analogous data from two other studies (QUANDRY [ 30 30 ] and TDTort [ 31 31 ] ). Temporal variations in the relative errors of the original and improved SIDK methods with respect to the HOBD method for the same outputs are depicted in Fig. 6 Fig. 6 . To avoid excessive indistinguishable graphs, due to the comparable accuracy of the kinetic methods developed for this study from a graphical perspective, the graphs associated with the direct kinetics and first-order SIDK methods have not been drawn in Figs. 5 Figs. 5 and 6 6 ; however, the results of all kinetic methods, along with equivalent data from various studies [ 29 29 - 31 31 29-31 ] , are directly compared in Tables 7 Tables 7 and 8 8 at every 0.1 s interval (from the start to the end of transients) for the step and ramp perturbations, respectively. 1 1 Time [s] 1 1 0.0 1 1 0.1 1 1 0.2 1 1 0.3 1 1 0.4 1 1 0.5 1 1 HOBD 1 1 1.000 1 1 2.06129 1 1 2.07859 1 1 2.09597 1 1 2.11350 1 1 2.13117 1 1 Direct kinetics 1 1 1.000 1 1 2.06135 1 1 2.07868 1 1 2.09606 1 1 2.11359 1 1 2.13126 1 1 Original SIDK 1 1 1.000 1 1 2.02818 1 1 2.04591 1 1 2.06791 1 1 2.08551 1 1 2.10132 1 1 Improved SIDK of first order 1 1 1.000 1 1 2.06132 1 1 2.07865 1 1 2.09604 1 1 2.11358 1 1 2.13125 1 1 Improved SIDK of second order 1 1 1.000 1 1 2.06130 1 1 2.07860 1 1 2.09598 1 1 2.11351 1 1 2.13117 1 1 QUANDRY [ 30 30 ] 1 1 1.000 1 1 2.061 1 1 2.078 1 1 2.095 1 1 2.113 1 1 2.131 1 1 CONQUEST [ 29 29 ] 1 1 1.000 1 1 2.060 1 1 2.078 1 1 2.095 1 1 2.112 1 1 2.130 1 1 TDTort [ 31 31 ] 1 1 1.000 1 1 2.061 1 1 2.079 1 1 2.096 1 1 2.114 1 1 2.131 1 1 Time [s] 1 1 0.0 1 1 0.1 1 1 0.2 1 1 0.3 1 1 0.4 1 1 0.5 1 1 HOBD 1 1 1.000 1 1 1.30849 1 1 1.95955 1 1 2.07549 1 1 2.09286 1 1 2.11037 1 1 Direct kinetics 1 1 1.000 1 1 1.30859 1 1 1.96004 1 1 2.07557 1 1 2.09295 1 1 2.11046 1 1 Original SIDK 1 1 1.000 1 1 1.30861 1 1 1.95858 1 1 2.06307 1 1 2.06885 1 1 2.07457 1 1 Improved SIDK of first order 1 1 1.000 1 1 1.30855 1 1 1.95975 1 1 2.07553 1 1 2.09293 1 1 2.11044 1 1 Improved SIDK of second order 1 1 1.000 1 1 1.30850 1 1 1.95960 1 1 2.07551 1 1 2.09287 1 1 2.11038 1 1 QUANDRY [ 30 30 ] 1 1 1.000 1 1 1.307 1 1 1.957 1 1 2.074 1 1 2.092 1 1 2.109 1 1 CONQUEST [ 29 29 ] 1 1 1.000 1 1 1.311 1 1 1.966 1 1 2.074 1 1 2.091 1 1 2.109 1 1 TDTort [ 31 31 ] 1 1 1.000 1 1 1.305 1 1 1.951 1 1 2.075 1 1 2.093 1 1 2.110 From Figs. 5 Figs. 5 and 6 6 , the higher accuracy of the improved SIDK methods developed in the current study compared to the original SIDK method is quite clear. The results also indicate that deviation of the original SIDK method is highly dependent on the type of perturbation. In the case of ramp perturbation, the results of the original SIDK method are nearly same as the other methods up to 0.169 s, but afterwards the deviation starts to grow in an almost linear (ramp) trend until it reaches a maximum error of 1.7% at 3 s; whereas the maximum deviation of the improved SIDK methods is about 0.06% in similar conditions. In the case of step perturbation, deviation of the results of the original SIDK method initiate from the start of computation in a rapid manner, while very low-magnitude oscillations are observable in the data pertaining to the initial jump of the graph and maximum deviation is more than 1.4%; whereas the maximum deviation of the improved SIDK methods is 0.012% in similar conditions. From Tables 7 Tables 7 and 8 8 , the HOBD and second-order SIDK methods demonstrate higher accuracy than the direct kinetics and first-order SIDK methods, respectively, which complies with expectations. However, the first-order SIDK method appears to be more accurate than the direct kinetics method, while the second-order SIDK method indicates lower precision compared to the HOBD method. These results imply that using the analytical solution with an average fission source term is more accurate than the first-order backward differentiation for the temporal term, but the reverse is true regarding the second-order backward differentiation. The total time required for conducting computations for both step and ramp perturbations using all the kinetic methods developed and employed in this study, along with the time-reduction ratios relative to the HOBD as the most accurate method and the portion of the time reduction of the original SIDK method retained by the improved SIDK methods, are tabulated in Table 9 Table 9 . 2 1 Method 1 1 HOBD 1 1 Direct kinetics 1 1 Improved SIDK of second order 1 1 Improved SIDK of first order 1 1 Original SIDK 7 1 Step perturbation 1 1 1 1 Total run-time [s] 1 1 129.7 1 1 125.8 1 1 122.6 1 1 116.9 1 1 111.7 1 1 1 1 Time-reduction ratio relative to the HOBD method % 1 1 NA 1 1 3.0 1 1 5.5 1 1 9.8 1 1 13.9 1 1 1 1 Ratio of run-time of improved SIDK methods relative to original SIDK method % 1 1 NA 1 1 NA 1 1 39.5 1 1 70.5 1 1 NA 7 1 Ramp perturbation 1 1 1 1 Total run-time [s] 1 1 141.1 1 1 136.9 1 1 131.2 1 1 127.0 1 1 121.3 1 1 1 1 Time-reduction ratio relative to the HOBD method % 1 1 NA 1 1 3.0 1 1 5.6 1 1 10.0 1 1 14.0 1 1 1 1 Ratio of run-time of improved SIDK methods relative to original SIDK method % 1 1 NA 1 1 NA 1 1 40.0 1 1 71.4 1 1 NA It should be confirmed that in the original SIDK method, the analytical solutions to the precursor equations are solved once per time-step with the quantity of fission source pertaining to the previous time-step; meanwhile, in the improved SIDK methods, the fission source term is updated in the inner iteration cycle of the point-implicit solver until the exact average quantity associated with the current time-step is achieved, and these facts justify the trends observed in results. While Fig. 6 Fig. 6 clearly indicates that the accuracy issue of the original SIDK method has been solved by the improvements suggested in this study, Table 9 Table 9 demonstrates that a significant portion of the computational time reduction of the original method (which is more than 71% for this specific benchmark problem) has also been retained. Temporal sequences of normalized contours of the thermal neutron flux and power density distribution during an accident are illustrated for step and ramp perturbations in Figs. 7 Figs. 7 and 8 8 , respectively. It is observed that since the TWIGL problem is a model of two subcritical reactivity insertion accidents, both perturbations predominantly affect the magnitudes of the variables and not their generic shapes of distribution. 4.2. LRA problem 2 Effective multiplication factors calculated for both fully inserted and fully withdrawn control-rod configurations associated with the LRA problem are presented in Table 10 Table 10 , along with the equivalent data from various studies [ 1 1 , 28 28 , 29 29 , 32 32 ] . Note that many of the available studies on the LRA problem did not report k eff for the second configuration. Nevertheless, the results indicate a maximum relative difference of 4 pcm (nearly 0.004%) between k eff calculated in the current study and the most accurate data reported in Ref. [ 28 28 , 29 29 , 32 32 ] for the initial configuration, and reported in Ref. [ 28 28 ] for the final configuration. 2 1 Reference 1 1 Current study 1 1 CUBBOX [28] 1 1 ISQBOX [28] 1 1 Shober [32] 1 1 CONQUEST [29] 1 1 Banfield [1] 8 1 k eff 1 1 1 1 For configuration of inserted rod (initial) 1 1 0.996325 1 1 0.99633 1 1 0.99631 1 1 0.99636 1 1 0.996329 1 1 0.99664 1 1 1 1 For configuration of withdrawn rod (final) 1 1 1.015453 1 1 1.01546 1 1 1.01541 1 1 ------- 1 1 ------- 1 1 -------- As for the transient state, variations in the average power density and temperature of the core versus time are demonstrated in Fig. 9 Fig. 9 along with the most accurate equivalent data from two other studies (CONQUEST [ 29 29 ] and Shober [ 32 32 ] ). Temporal variations in the relative errors of the original and improved SIDK methods with respect to the HOBD method for the same outputs are depicted in Fig. 10 Fig. 10 . Due to the comparable accuracy of the kinetic methods developed for this study from a graphical perspective, the graphs associated with the direct kinetics and first-order SIDK methods have not been drawn in Fig. 8 Fig. 8 in order to avoid excessive indistinguishable data; however, the prime results of accident analysis pertaining to all kinetic methods, along with equivalent data from two other studies [ 29 29 , 32 32 ] , are directly compared in Table 11 Table 11 . 1 1 Method 1 1 HOBD 1 1 Direct kinetics 1 1 Original SIDK 1 1 Improved SIDK of first order 1 1 Improved SIDK of second order 1 1 Shober [32] 1 1 CONQUEST [29] 1 1 Time to first peak [s] 1 1 1.435 1 1 1.431 1 1 1.428 1 1 1.433 1 1 1.434 1 1 1.436 1 1 1.437 1 1 Mean power at first peak [W/cm 3 ] 1 1 5391.8 1 1 5359.94 1 1 5052.88 1 1 5385.1 1 1 5390.01 1 1 5411 1 1 5505 1 1 Time to second peak [s] 1 1 1.999 1 1 1.995 1 1 1.983 1 1 1.997 1 1 1.999 1 1 2.0 1 1 2.0 1 1 Mean power at second peak [W/cm 3 ] 1 1 781.08 1 1 775.4 1 1 697.58 1 1 779.95 1 1 780.76 1 1 784.11 1 1 798 1 1 Maximum temperature at t = 3.0 s [K] 1 1 2867.2 1 1 2855.04 1 1 2709.68 1 1 2860.56 1 1 2866.11 1 1 2948 1 1 3023 1 1 Average temperature at t = 3.0 s [K] 1 1 1069.69 1 1 1066.55 1 1 1011.93 1 1 1067.52 1 1 1068.01 1 1 1087.32 1 1 1107 1 1 Mean power at t = 3.0 s [W/cm 3 ] 1 1 94.61 1 1 92.56 1 1 79.46 1 1 93.61 1 1 94.27 1 1 96.18 1 1 99.1 Please note that according to the instruction of Ref. [ 28 28 ], a logarithmic scale should be used for the ordinate or y -axis of the graph presented in Fig. 9 (a) Fig. 9 (a) , due to the notable intensity of the accident and the consequently large range of variation in the core power density associated with the LRA problem, which is clearly observable in this figure. As a result, discrepancies among the different sets of results in Fig. 9 (a) Fig. 9 (a) are more distinguishable for the lower quantities than the higher ones; however, the usual linear scale has been used for the other graphs represented in this study, including Fig. 9 (b) Fig. 9 (b) and Fig. 10 Fig. 10 , which clearly demonstrate how the accuracy issue of the original SIDK method has been resolved due to improvements suggested in this study for this specific benchmark problem. From Figs. 9 Figs. 9 and 10 10 and Table 11 Table 11 , the kinetic methods employed in this study could be arranged as HOBD, second-order SIDK, first-order SIDK, direct kinetics, and original SIDK in descending order of accuracy relative to the reference data [ 29 29 , 32 32 ] , which is in agreement with the order obtained for the TWIGL problem; hence, the same justification mentioned in Sec. 4.1 applies to this case as well. Moreover, similar to the case of ramp perturbation pertaining to the TWIGL problem, the results of the original SIDK method are very close to those of the other methods in the first part of the graphs ( Fig. 9 Fig. 9 ) associated with the initial jump (which terminates at about 1.44 s in the accident modeling time), but afterwards the deviation increases in an almost linear (ramp) trend until it reaches maximum errors of 17.4% and 6.9% at the end of the simulation time (3 s) for average power density and temperature of the core, respectively; whereas the maximum deviations of the improved SIDK methods are about 2.67% and 1.82% in similar conditions (based on Table 11 Table 11 ). Note that the quantity of errors of the original SIDK method has drastically increased in comparison with their equivalents in the TWIGL problem (Sec. 4.1), mainly due to the fact that the magnitude of the LRA accident is considerably greater than that of the TWIGL accident (a super-prompt critical reactivity insertion accident as compared to subcritical ones). The total time required for conducting computations for both average power density and temperature of the core using all kinetic methods developed and employed in this study, along with the time-reduction ratios relative to the HOBD as the most accurate method and the portion of the time-reduction of the original SIDK method retained by the improved SIDK methods, are tabulated in Table 12 Table 12 . 1 1 Method 1 1 HOBD 1 1 Direct kinetics 1 1 Improved SIDK of second order 1 1 Improved SIDK of first order 1 1 Original SIDK 1 1 Total run-time [s] 1 1 817.2 1 1 769.8 1 1 726.5 1 1 613.7 1 1 562.2 1 1 Time-reduction ratio relative to the HOBD method % 1 1 NA 1 1 5.8 1 1 11.1 1 1 24.9 1 1 31.2 1 1 Ratio of improved SIDK methods run-time relative to original SIDK % 1 1 NA 1 1 NA 1 1 35.6 1 1 79.8 1 1 NA It is observed that the order of time reductions of the kinetic methods is also in accordance with the order obtained for the TWIGL problem; however, the quantities of reduction ratios are notably larger than their analogues in Sec. 4.1 for two major reasons. First, the magnitude of the LRA accident is considerably higher than that of the TWIGL accident, and second, the LRA problem has two groups of delayed neutron precursors while the TWIGL problem has only one group for which the analytical solution obtained in Sect. 2.3 (Eq. (11)) should be employed. In other words, the advantages of the SIDK methods become more apparent as the magnitude of the target accident and number of precursor groups increase, but at the same time, the accuracy of the original SIDK method also reduces drastically, whereas that is not the case with regard to the improved SIDK methods developed in this study. While Fig. 10 Fig. 10 clearly indicates that the accuracy issue of the original SIDK method has been solved by the improvements suggested in this study, Table 12 Table 12 demonstrates that a significant portion of the computational time reduction of the original method (which is nearly 80% for this specific benchmark problem) has also been retained. The temporal sequence of normalized contours of power density distribution during the LRA accident is illustrated in Fig. 11 Fig. 11 . It is observed that, in contrast to the TWIGL problem (Sect. 4.1), the LRA accident considerably affects both the magnitudes of variables and their generic shapes of distribution. The reason also lies in the difference in magnitude of the accidents associated with these two problems. Adequate agreement between the outputs and reference data proves the fidelity and validity of the kinetics methods and the code developed in the current research for these sorts of severe conditions. The temporal sequence of contours of temperature distribution during the accident is illustrated in Fig. 12 Fig. 12 . With such detailed and comprehensive outputs, not only the critical times of the accident but also the most vulnerable positions of the reactor to effects of the imposed perturbation (which is region 4 in the case of the LRA accident (see Fig. 4 Fig. 4 for the region assignment)) could be recognized in the clearest and most accurate manner. It should be noted that in almost all available studies that have addressed the LRA accident, including the ones [ 29 29 , 32 32 ] used as the frame of reference in the current study, the nodal method has been employed for spatial discretization of the neutronic equations, which demands a low-density numerical gird; however, the finite-difference method with the same low number of meshes has been used for discretization of the adiabatic heat-up equation (Eq. (16)). As the FVM with a high-density numerical grid has been employed in this study (see Sect. 3.2) for the same purpose, it is expected that the results of temperature distribution (and hence the whole output of this study concerning the LRA problem due to the effect of feedbacks via Eq. (17)) should be more accurate than in the references. This matter justifies the higher deviations mentioned in the second paragraph of this section for the improved SIDK methods developed in this study as compared to their analogues in Sect. 4.1, because the TWIGL problem does not involve computation of feedback effects or temperature distribution. Nevertheless, the maximum error of less than 3% should be acceptable for most practical applications. According to the results obtained at different conditions of both the TWIGL and LRA problems, and in view of both accuracy and computation time, the first-order SIDK has emerged as the best method among the five kinetic methods developed and evaluated in the current study. While the accuracy of this method is quite comparable with the HOBD method (as the most accurate method), it has reduced the computational time by up to 24.9% for the selected benchmark problems, while only two groups of delayed neutron precursor groups were involved in these problems at most. Accordingly, this reduction ratio is expected to increase notably with severity and details of the target accident, and number of delayed neutron precursor groups involved in the simulation. 5. Conclusion 1 SIDK is a reactor kinetic method developed by Banfield to optimize temporal discretization of neutronic equations. The main assumption of the SIDK method is to substitute a time-averaged quantity for a fission source term which appears in the delayed neutron equations. As a result, these equations are decoupled from prompt neutron equations, and analytical solutions to the delayed equations are obtained which leads to significant reduction in computational cost. The problem is that the fission source averaged over a time-step is not known, so the original study [ 14 14 ] suggested using the constant quantity pertaining to the previous time-step for this purpose, and a reduction in the size of the time-step was proposed to lessen the imposed errors, which drastically escalate with magnitude of the induced perturbation. However, this remedy leads to an increase in computational cost, which nullifies the main advantage of the SIDK method. In this study, we developed a novel multi-group, multi-dimensional diffusion code using FVM and a point-implicit solver, and inserted the SIDK’s analytical relations into the numerical iteration cycle of the point-implicit solver as explicit steps, while the quantity of the fission source was updated, according to a modified formulation, using the most recent values of group fluxes obtained from solving the prompt neutron equations. In addition to SIDK, two other kinetic methods, i.e. direct kinetics and HOBD, were programmed and the resultant code was tested at different conditions of two well-known benchmark problems, modeling subcritical to super-prompt critical reactivity insertion accidents with one and two groups of delayed neutron precursors. The results indicated that since the method proposed in this article resulted in attaining a much more accurate estimation of the time-averaged fission source in a time-step, the accuracy issue of the original SIDK method was adequately resolved, while a significant portion of the computational cost reduction of this method was saved, which was evaluated to be up to 80% (or more than 24% of the total computational time) for the cases tested in this study. However, according to the comparison of results obtained at different sets of conditions, the time-reduction ratio is expected to increase notably with severity of the target accident and number of the delayed neutron precursor groups involved in the simulation. References J. E. Banfield , "Semi-implicit direct kinetics methodology for deterministic, time-dependent, and fine-energy neutron transport solutions", Ph.D. diss . University of Tennessee , 2013 , http://trace.tennessee.edu/utk.graddiss/1694 http://trace.tennessee.edu/utk.graddiss/1694 . "Semi-implicit direct kinetics methodology for deterministic, time-dependent, and fine-energy neutron transport solutions", Ph.D. diss AMSO: Advanced Modeling and Simulation Office. (the USA DOE/NE) , http://energy.gov/ne/nuclear-reactor-technologies/advanced-modeling-simulation http://energy.gov/ne/nuclear-reactor-technologies/advanced-modeling-simulation , Accessed 17 Dec 2014 . AMSO: Advanced Modeling and Simulation Office. (the USA DOE/NE) NEAMS: Nuclear Energy Advanced Modeling and Simulation program. (the USA DOE/NE) , http://neams.ne.anl.gov http://neams.ne.anl.gov , Accessed 17 Dec 2014 . NEAMS: Nuclear Energy Advanced Modeling and Simulation program. (the USA DOE/NE) CASL: The Consortium for Advanced Simulation of Light Water Reactors. (the USA DOE/NE) , http://www.casl.gov http://www.casl.gov , Accessed 17 Dec 2014 . CASL: The Consortium for Advanced Simulation of Light Water Reactors. (the USA DOE/NE) F-BRIDGE: Framework of Basic Research for Innovative Fuel Design for GEn IV systems project. (the European Commission INSC project) , http://www.f-bridge.eu http://www.f-bridge.eu , Accessed 17 Dec 2014 . F-BRIDGE: Framework of Basic Research for Innovative Fuel Design for GEn IV systems project. (the European Commission INSC project) 3D flow, ERCOSAM-SAMRA and RSQUD projects. (Russian Federation, Nuclear Safety Institute of the Russian Academy of Science (IBRAE)) , http://www.ibrae.ac.ru http://www.ibrae.ac.ru , Accessed 17 Dec 2014 . 3D flow, ERCOSAM-SAMRA and RSQUD projects. (Russian Federation, Nuclear Safety Institute of the Russian Academy of Science (IBRAE)) SMART project: Integrated Modular Advanced Reactor project, (Korea Atomic Energy Research Institute) , https://www.kaeri.re.kr https://www.kaeri.re.kr , Accessed 17 Dec 2014 . SMART project: Integrated Modular Advanced Reactor project, (Korea Atomic Energy Research Institute) IVNET: Innovative and Viable Nuclear Energy Technology Development project. (Japan Atomic Energy Agency, Research Group for Thermal and Fuel Energy) , https://nsec.jaea.go.jp https://nsec.jaea.go.jp , Accessed 17 Dec 2014 . IVNET: Innovative and Viable Nuclear Energy Technology Development project. (Japan Atomic Energy Agency, Research Group for Thermal and Fuel Energy) R. Vadi , K. Sepanloo . Investigation of a LOCA in a typical MTR by a novel best-estimate code . Progress in Nuclear Energy 86 , 141 - 161 ( 2016 ). doi: 10.1016/j.pnucene.2015.10.013 http://doi.org/10.1016/j.pnucene.2015.10.013 Investigation of a LOCA in a typical MTR by a novel best-estimate code R. Vadi , K. Sepanloo , An improved porous media model for nuclear reactor analysis . Nuclear Science and Techniques 27 , 1 - 24 ( 2016 ). doi: 10.1007/s41365-016-0016-7 http://doi.org/10.1007/s41365-016-0016-7 An improved porous media model for nuclear reactor analysis R. Vadi , K. Sepanloo , Reassessment of the generic assumptions applied to the conventional analysis of the reactivity insertion accident in the MTRs using a novel coupled code . Progress in Nuclear Energy 93 , 96 - 115 ( 2016 ). doi: 10.1016/j.pnucene.2016.08.008 http://doi.org/10.1016/j.pnucene.2016.08.008 Reassessment of the generic assumptions applied to the conventional analysis of the reactivity insertion accident in the MTRs using a novel coupled code IAEA-CRP (Coordinated Research Project), FUMAC: Fuel Modeling in Accident Condition, 2014 to 2018 . https://www.iaea.org/projects/crp/t12028 https://www.iaea.org/projects/crp/t12028 . IAEA-CRP (Coordinated Research Project), FUMAC: Fuel Modeling in Accident Condition, 2014 to 2018 IAEA technical meeting on modelling of fuel behavior in design basis accidents and design extension conditions (EVT1803345), Shenzhen, China, 13-16 May 2019 . https://www.iaea.org/events/evt1803345 https://www.iaea.org/events/evt1803345 . IAEA technical meeting on modelling of fuel behavior in design basis accidents and design extension conditions (EVT1803345), Shenzhen, China, 13-16 May 2019 J. E. Banfield , et al ., A new semi-implicit direct kinetics method with analytical representation of delayed neutrons . In proceeding of 2012 ANS Annual Winter Meeting , San Diego, CA , November 11-15, 2012 , vol. 107 , pp. 1111 - 1114 . A new semi-implicit direct kinetics method with analytical representation of delayed neutrons D. Ginestar , High-order backward discretization of the neutron diffusion equation . Ann. Nuclear Energy 25 ( 3 ), 47 - 644 ( 1998 ). doi: 10.1016/S0306-4549(97)00046-7 http://doi.org/10.1016/S0306-4549(97)00046-7 High-order backward discretization of the neutron diffusion equation A. J. Chorin , Numerical solution of Navier–Stokes equations . Mathematics of Computation 22 , 745 - 762 ( 1968 ). doi: 10.1090/S0025-5718-1968-0242392-2 http://doi.org/10.1090/S0025-5718-1968-0242392-2 Numerical solution of Navier–Stokes equations R. J. Leveque . " Finite Volume Methods for Hyperbolic Problems ", ( Cambridge University Press , UK , 2004 ), pp. 1 - 558 . Finite Volume Methods for Hyperbolic Problems D. L. Hetric , " Dynamics of nuclear reactors ", ( University of Chicago Press , Chicago , 1971 ), pp. 87 - 104 . Dynamics of nuclear reactors ANSYS, Inc , " ANSYS Fluent Theory Guide, Release 19.2, August 2018, Chapter 21 , pp. 687 ", Canonsburg, Pennsylvania, USA , http://www.ansys.com http://www.ansys.com . ANSYS Fluent Theory Guide, Release 19.2, August 2018, Chapter 21 Y. Saad , " Numerical Methods for Large Eigenvalue Problems ", ( Halstead Press , New York , 1992 ), pp. 62 - 77 . Numerical Methods for Large Eigenvalue Problems W. H. Press , S. A. Teukolsky , W. T. Vetterling , B. P. Flannery , " Numerical Recipes in C ", 2 nd edition, ( Cambridge University Press , UK , 1992 ), pp. 312 - 347 . Numerical Recipes in C W. Anderson , D. L. Bonhus , An Implicit Upwind Algorithm for Computing Turbulent Flows on Unstructured Grids , Computers Fluids 23 ( 1 ), 1 - 21 ( 1994 ). doi: 10.1016/0045-7930(94)90023-X http://doi.org/10.1016/0045-7930(94)90023-X An Implicit Upwind Algorithm for Computing Turbulent Flows on Unstructured Grids Z. J. Wang , A Fast Nested Multi-grid Viscous Flow Solver for Adaptive Cartesian/Quad Grids . International Journal for Numerical Methods in Fluids 33 , 657 - 680 ( 2000 ). doi: 10.1002/1097-0363(20000715)33:5<657::AID-FLD24>3.0.CO;2-G http://doi.org/10.1002/1097-0363(20000715)33:5<657::AID-FLD24>3.0.CO;2-G A Fast Nested Multi-grid Viscous Flow Solver for Adaptive Cartesian/Quad Grids B. R. Hutchinson , G. D. Raithby , A Multi-grid Method Based on the Additive Correction Strategy . Numerical Heat Transfer 9 , 511 - 537 ( 1986 ). doi: 10.1080/10407788608913491 http://doi.org/10.1080/10407788608913491 A Multi-grid Method Based on the Additive Correction Strategy T. J. Trahan . "An asymptotic, homogenized, anisotropic, multi-group, diffusion approximation to the neutron transport equation". Ph.D. diss . University of Michigan , Ann Arbor, MI, USA , 2014 . "An asymptotic, homogenized, anisotropic, multi-group, diffusion approximation to the neutron transport equation". Ph.D. diss R. L. Burden , J. D. Faires , " Numerical Analysis ", Ninth Edition ( Brooks/Cole Press , Florence , 2010 ), pp. 101 - 146 . Numerical Analysis L. A. Hageman , J. B. Yasinsky . Comparison of alternating-direction time-differencing methods and other implicit methods for the solution of the neutron group diffusion equations . Nuclear Sci. Eng . 38 , 8 ( 1969 ). doi: 10.13182/NSE38-8 http://doi.org/10.13182/NSE38-8 Comparison of alternating-direction time-differencing methods and other implicit methods for the solution of the neutron group diffusion equations Computational benchmark committee of the mathematics and computation division of the American nuclear society , "National energy software centre: benchmark problem book", ANL Technical Report, ANL-7416, supplement 2 and 3 , Argonne, Illinois, USA , 1985 . "National energy software centre: benchmark problem book", ANL Technical Report, ANL-7416, supplement 2 and 3 J. C. Gehin . "A quasi-static polynomial nodal method for nuclear reactor analysis", Ph.D. diss . MIT University , Massachusetts, USA , 1992 . "A quasi-static polynomial nodal method for nuclear reactor analysis", Ph.D. diss K. S. Smith . "An analytical nodal method for two-group multi-dimensional static and transient neutron diffusion equations", M.S. thesis . MIT University , Massachusetts, USA , 1979 . "An analytical nodal method for two-group multi-dimensional static and transient neutron diffusion equations", M.S. thesis S. Goluoglu , H. L. Dodds , A time-dependent, three-dimensional neutron transport methodology . Nuclear Science and Engineering 139 , 248 - 261 ( 2001 ). doi: 10.13182/NSE01-A2235 http://doi.org/10.13182/NSE01-A2235 A time-dependent, three-dimensional neutron transport methodology R. A. Shober . A nodal method for solving transient fewgroup neutron diffusion equations. ANL Technical Report. ANL-78-51 , Argonne, Illinois, USA , 1978 . A nodal method for solving transient fewgroup neutron diffusion equations. ANL Technical Report. ANL-78-51 10.1007/s41365-019-0690-3 2019-03-05 2019-05-31 2019-06-26 2019-10-21 2019-11-01 2021-04-11 © 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. 2019 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-019-0690-3 http://dx.doi.org/10.1007/s41365-019-0690-3 . 2019 Nuclear Science and Techniques 11 30