Frequency-transformed equation

The dynamic frequency of the scalar flux is defined as ${\omega }_g\left(r\mathrm{,}t\right)\equiv \frac{1}{{\varphi }_g\left(r\mathrm{,}t\right)}\frac{\partial }{\partial t}{\varphi }_g\left(r\mathrm{,}t\right),$ (1) where the scalar flux ${\varphi }_g\left(r\mathrm{,}t\right)$ of group g is given by ${\varphi }_g(r\mathrm{,}t)={\displaystyle \int }\text{d}\text{Ω}{\psi }_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right),$ (2) where ${\psi }_g\left(\text{r}\mathrm{,}\text{Ω}\mathrm{,}t\right)$ is the angular flux at position 𝒓 in direction $\text{Ω}$. Note that to remove factors of π from these equations, dΩ is normalized such that ${\displaystyle \int }\text{d}\text{Ω}=1$. Generally, it is assumed that the dynamic frequency of the angular flux is isotropic, such that ${\omega }_g\left(r\mathrm{,}t\right)=\frac{1}{{\psi }_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right)}\frac{\partial }{\partial t}{\psi }_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right).$ (3) Then, the time-variant angular flux can be expressed as ${\psi }_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right)={\psi }_g\left(r\mathrm{,}\text{Ω}\mathrm{,}{t}_0\right)e^{{\displaystyle {\int }_{t_0}^t}d{t}^{\text{'}}{\omega }_g\left(\text{r}\mathrm{,}{t}^{\text{'}}\right)}.$ (4) The dynamic frequency can be further decomposed into ${\omega }_g(r\mathrm{,}t)={\omega }_{S\mathrm{,}g}(r\mathrm{,}t)+{\omega }_{\text{T}}(t)$[21]. The flux amplitude frequency ω_T(t) represents a global quantity and depends only on time, and the flux shape frequency ${\omega }_{\text{S}\mathrm{,}g}(r\mathrm{,}t)$ depends on space, time, and energy. With this approximation, the angular flux can be decomposed using the amplitude and shape terms as ${\psi }_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right)=P\left(t\right){\widehat{\psi }}_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right),$ (5) where P(t) is the flux amplitude and ${\widehat{\psi }}_g\left(r\mathrm{,}\text{Ω}\mathrm{,}{t}_0\right)$ is the flux shape. The time-variant flux amplitude and flux shape are evaluated as $P(t)=P(t_0)e^{{\displaystyle {\int }_{t^0}^t}\text{d}{t}^{\text{'}}{\omega }_{\text{T}}(t^{\text{'}})},$ (6) ${\widehat{\psi }}_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right)={\widehat{\psi }}_g\left(r\mathrm{,}\text{Ω}\mathrm{,}{t}_0\right)e^{{\displaystyle {\int }_{t_0}^t}d{t}^{\text{'}}{\omega }_{S\mathrm{,}g}\left(\text{r}\mathrm{,}{t}^{\text{'}}\right)}.$ (7) To make the decomposition unique[21], a physics-based constraint on the shape frequency is introduced, as follows: ${\displaystyle \sum _{g=1}^G{\displaystyle {\int }_V}\text{d}r}\nu {\Sigma }_{f\mathrm{,}g}(r\mathrm{,}t){\displaystyle \int }d\text{Ω}{\widehat{\psi }}_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right)=P\left(t_0\right),$ (8) where ν is the number of released neutrons per fission reaction and ∑_fg(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 ${\mu }_i(r\mathrm{,}t)\equiv \frac{1}{C_i(r\mathrm{,}t)}\frac{\partial }{\partial t}{C}_i(r\mathrm{,}t),$ (9) where Ci(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] with isotropic scattering, the following frequency-transformed neutron transport equations are obtained: $\begin{array}{l}\frac{{\omega }_{\text{T}}\left(t\right)+{\omega }_{\text{S}g}\left(r\mathrm{,}t\right)}{v_g\left(r\right)}{\psi }_g(r\mathrm{,}\text{Ω}\mathrm{,}t)+\text{Ω}\cdot \nabla {\psi }_g(r\mathrm{,}\text{Ω}\mathrm{,}t)\\ +{\displaystyle {\sum }_{\text{t}g}(r\mathrm{,}t)}{\psi }_g(r\mathrm{,}\text{Ω}\mathrm{,}t)\\ ={\displaystyle {\sum }_{\text{s}g}(r\mathrm{,}t)}{\varphi }_g(r\mathrm{,}t)+q_g(r\mathrm{,}t)g=\mathrm{1,}\cdots \mathrm{,}G\end{array}$ (10) where ${\displaystyle {\sum }_{\text{t}g}(r\mathrm{,}t)}$ and ${\displaystyle {\sum }_{\text{s}g}(r\mathrm{,}t)}$ are the macroscopic total and scattering cross sections, respectively. The group source qg(r,t) is the contribution from scattering, prompt, and delayed neutrons: $\begin{array}{c}{q}_g\left(r\mathrm{,}t\right)={\displaystyle \sum _{g^{\text{'}}\ne g}{\displaystyle {\sum }_{\text{s}g{g}^{\prime }}(r\mathrm{,}t){\varphi }_{g^{\prime }}(r\mathrm{,}t)}}\\ +\left\{\begin{array}{c}{\chi }_g\left(r\right)\left[1-\beta \left(r\right)\right]\\ +{\displaystyle {\sum }_{i=1}^I\frac{{\chi }_{ig}\left(r\right){\beta }_i\left(r\right){\lambda }_i\left(r\right)}{{\mu }_i\left(r\mathrm{,}t\right)+{\lambda }_i\left(r\right)}}\end{array}\right\}\\ \times {\displaystyle \sum _{g^{\prime }}\nu {\displaystyle {\sum }_{f{g}^{\prime }}(r\mathrm{,}t){\varphi }_{g^{\prime }}(r\mathrm{,}t)}}\end{array}$ (11) where ${\displaystyle {\sum }_{\text{s}g{g}^{\prime }}\left(r\mathrm{,}t\right)}$ is the macroscopic scattering cross-section, ${\chi }_g\left(r\right)$ and ${\chi }_{ig}\left(r\right)$ are the fission spectra of prompt and delayed neutrons, ${\beta }_i\left(r\right)$ and $\beta \left(r\right)$ are the fractions of delayed neutrons, and ${\lambda }_i\left(r\right)$ is the delayed constant of the precursor. The frequency-transformed precursor equations are given as $\begin{array}{l}{\mu }_i(r\mathrm{,}t)C_i(r\mathrm{,}t)={\beta }_i\left(r\right){\displaystyle \sum _{g^{\prime }}\nu }{\displaystyle {\sum }_{f{g}^{\prime }}(r\mathrm{,}t)}{\varphi }_{g^{\prime }}(r\mathrm{,}t)\\ -{\lambda }_i\left(r\right)C_i(r\mathrm{,}t)i=\mathrm{1,...,}I\end{array}$ (12) 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 $\begin{array}{l}\text{Ω}\cdot \nabla {\psi }_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right)+{\displaystyle {\sum }_{\text{t}g}(r\mathrm{,}t)}{\psi }_g(r\mathrm{,}\text{Ω}\mathrm{,}t)\\ ={\displaystyle {\sum }_{\text{s}g}(r\mathrm{,}t)}{\varphi }_g(r\mathrm{,}t)+q_g^{\text{'}}(r\mathrm{,}t),\end{array}$ (13) $\begin{array}{l}{q^{\prime }}_g(r\mathrm{,}t)={\displaystyle \sum _{g^{\text{'}}\ne g}{\displaystyle {\sum }_{\text{s}g{g}^{\prime }}(r\mathrm{,}t)}}{\varphi }_{g^{\text{'}}}(r\mathrm{,}t)\\ -\frac{{\omega }_{\text{T}}\left(t\right)+{\omega }_{\text{S}g}\left(r\mathrm{,}t\right)}{v_g\left(r\right)}{\varphi }_g(r\mathrm{,}t)\\ +k_{\text{D}}^{-1}{{\chi }^{\prime }}_g\left(r\mathrm{,}t\right){\displaystyle \sum _{g^{\prime }}\nu }{\displaystyle {\sum }_{f{g}^{\prime }}(r\mathrm{,}t)}{\varphi }_{g^{\prime }}(r\mathrm{,}t),\end{array}$ (14) where ${\chi }_g^{\text{'}}\left(r\mathrm{,}t\right)$ is the dynamic fission spectrum, defined as ${\chi }_g^{\text{'}}(r\mathrm{,}t)\equiv {\chi }_g(r)\left[1-\beta \left(r\right)\right]+{\displaystyle \sum _{i=1}^I\frac{{\chi }_{ig}\left(r\right){\beta }_i\left(r\right){\lambda }_i\left(r\right)}{{\mu }_i(r\mathrm{,}t)+{\lambda }_i\left(r\right)}}.$ (15) 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], and is illustrated in Fig. 1.

Fig. 1Full-screen View

Flowchart of the k-ω algorithm

Variational nodal formulation

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.

2.2.1

The discretized functional

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] $\begin{array}{c}{F}_{\nu }\left[{\psi }^{+}\mathrm{,}{\psi }^{-}\right]={\displaystyle {\int }_{\nu }}\text{d}r\left\{\begin{array}{c}{\displaystyle \int }\text{d}\text{Ω}\left[{\displaystyle {\sum }_{\text{t}}^{\text{'}-1}{\left(\text{Ω}\cdot \nabla {\psi }^{+}\right)}^2}+{\displaystyle {\sum }_{\text{t}}^{\text{'}}{\psi }^{+2}}\right]\\ -{\Sigma }_{\text{s}}{\varphi }^2-2\varphi q^{\prime }\end{array}\right\}\\ +2{\displaystyle \int }\text{d}z{\displaystyle {\int }_{\Gamma }}\text{d}\Gamma {\displaystyle \int }\text{d}\text{Ω}{n}_p\cdot \text{Ω}{\psi }^{+}{\psi }^{-}\\ +2{\displaystyle {\int }_A}\text{d}A{\displaystyle \int }\text{d}\text{Ω}\left(\begin{array}{c}{n}_{z^{+}}\cdot \text{Ω}{\psi }^{+}{\psi }^{-}{|}_{z^{+}}\\ +n_{z^{-}}\cdot \text{Ω}{\psi }^{+}{\psi }^{-}{|}_{z^{-}}\end{array}\right).\end{array}$ (16) 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 $n_{z^{+}}$ and $n_{z^{-}}$ 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 $F\left[{\psi }^{+}\mathrm{,}{\psi }^{-}\right]={\displaystyle {\sum }_{\nu }{F}_{\nu }}\left[{\psi }^{+}\mathrm{,}{\psi }^{-}\right]$.

Within the node, the spatial distribution of the even-parity flux is approximated by ${\psi }^{+}(r\mathrm{,}\text{Ω}\mathrm{,}t)\approx f_z^{\text{T}}(z)\otimes {\text{g}}^{\text{T}}\left(x\mathrm{,}y\right)\psi (\text{Ω}\mathrm{,}t)\mathrm{.}$ (17) The scalar flux can be expanded as $\varphi (r\mathrm{,}t)=f_z^T(z)\otimes {\text{g}}^T\left(x\mathrm{,}y\right)\varphi \left(t\right)\mathrm{.}$ (18) where the time-dependent scalar flux moments are defined as $\varphi \left(t\right)={\displaystyle \int }\text{d}\text{Ω}\psi (\text{Ω}\mathrm{,}t)$. $f_z\left(z\right)$ is a vector of orthogonal polynomials defined at node ν. The polynomials are governed by ${\displaystyle \int }\text{d}z{f}_z^{\text{T}}\left(z\right)f_z\left(z\right)=\Delta zI,$ (19) where 𝑰 denotes the identity matrix. $g\left(x\mathrm{,}y\right)$ is a vector of continuous FE basis functions. Triangular and quadrilateral iso-parametric quadratic FEs [11] were employed to map the curved interfaces between materials. Figure 2 shows the FE mesh with 32 elements, used for describing a fuel pin cell in this paper.

Fig. 2Full-screen View

(Color online) An FE mesh for a fuel pin. The orange and blue colors represent the fuel region and the moderator region, respectively

The odd-parity flux at the axial interface is approximated as $\begin{array}{c}{\psi }^{-}(r\mathrm{,}\text{Ω}\mathrm{,}t)\approx f_z^{\text{T}}\left(\pm \Delta z/2\right){\text{y}}_z^{\text{T}}\left(\text{Ω}\right)\otimes h^{\text{T}}\left(x\mathrm{,}y\right){\chi }_z\left(t\right)\\ +f_z^{\text{T}}\left(\pm \Delta z/2\right){\text{y}}_{zo}^{\text{T}}\left(\text{Ω}\right)\otimes h^{\text{T}}\left(x\mathrm{,}y\right){\chi }_{zo}\left(t\right),\end{array}$ (20) where $h\left(x\mathrm{,}y\right)$ denotes a vector of piecewise constants. $y_z^{\text{T}}\left(\text{Ω}\right)$ and $y_{zo}^{\text{T}}\left(\text{Ω}\right)$ are odd-parity spherical harmonic vectors comprising both sine and cosine functions[24]. The vectors $y_z^{\text{T}}\left(\text{Ω}\right)$ are low-order Pn approximations, whereas $y_{zo}^{\text{T}}\left(\text{Ω}\right)$ and $y_{\gamma o}^{\text{T}}\left(\text{Ω}\right)$ contain higher-order terms ranging from n+2 to a larger number N. In this study, we refer to these as PNPn approximations. ${\chi }_z\left(t\right)$ and ${\chi }_{zo}\left(t\right)$ 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 $\begin{array}{c}{\psi }^{-}\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right)={\text{y}}_{\gamma }^{\text{T}}\left(\text{Ω}\right)\otimes f_{\gamma }^{\text{T}}\left(\zeta \right){\chi }_{\gamma }\left(t\right)\\ +y_{\gamma o}^{\text{T}}\left(\text{Ω}\right)\otimes f_{\gamma o}^{\text{T}}\left(\zeta \right){\chi }_{\gamma o}\left(t\right)\\ \gamma =x\mathrm{,}y\zeta =y\mathrm{,}x\end{array}$ (21) where γ =1,2,3,4 corresponds to the γ =x^-,x⁺,y^-,y⁺ lateral interface. In VITAS, $f_{\gamma }^{\text{T}}\left(\zeta \right)$ and $f_{\gamma o}^{\text{T}}\left(\zeta \right)$ are chosen as fourth-order polynomials. Similarly, $y_{\gamma }^{\text{T}}\left(\text{Ω}\right)$ and $y_{\gamma o}^{\text{T}}\left(\text{Ω}\right)$ are low- and high-order Pn approximations, while ${\chi }_{\gamma }\left(t\right)$ and ${\chi }_{\gamma o}\left(t\right)$ 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: $\begin{array}{c}F\left[\psi \mathrm{,}{\chi }_{\gamma }\mathrm{,}{\chi }_z\right]={\displaystyle \int }\text{d}\text{Ω}{\psi }^{\text{T}}(\text{Ω}\mathrm{,}t)\text{A}(\text{Ω}\mathrm{,}t)\psi (\text{Ω}\mathrm{,}t)\\ -{\varphi }^{\text{T}}\left(t\right){\text{I}}_z\otimes {\text{F}}_s\left(t\right)\varphi \left(t\right)-2\varphi {\left(t\right)}^{\text{T}}{q}^{\prime }\left(t\right)\\ +2{\displaystyle \sum _{\gamma }{\displaystyle \int }\text{d}\text{Ω}{\psi }^{\text{T}}}(\text{Ω}\mathrm{,}t)E_{\gamma }(\text{Ω}){\chi }_{\gamma }\left(t\right)\\ +2{\displaystyle \sum _{\gamma }{\displaystyle \int }\text{d}\text{Ω}{\psi }^{\text{T}}}(\text{Ω}\mathrm{,}t)E_z(\text{Ω}){\chi }_z\left(t\right)\\ +2{\displaystyle \sum _{\gamma }{\displaystyle \int }\text{d}\text{Ω}{\psi }^{\text{T}}}(\text{Ω}\mathrm{,}t)E_{\gamma o}(\text{Ω}){\chi }_{\gamma o}\left(t\right)\\ +2{\displaystyle \sum _{\gamma }{\displaystyle \int }\text{d}\text{Ω}{\psi }^{\text{T}}}(\text{Ω}\mathrm{,}t)E_{zo}(\text{Ω}){\chi }_{zo}\left(t\right),\end{array}$ (22) $q^{\prime }\left(t\right)={\displaystyle \int }\text{d}r{f}_z\left(z\right)\otimes g\left(x\mathrm{,}y\right)q^{\text{'}}\left(r\mathrm{,}t\right),$ (23) 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, 12]. Requiring the discretized functional in Eq. 22 to be stationary with respect to variations in $\psi (\text{Ω}\mathrm{,}t)$ yields: $\begin{array}{l}A(\text{Ω}\mathrm{,}t)\psi (\text{Ω}\mathrm{,}t)-I_z\otimes F_s\varphi \left(t\right)\\ =q^{\text{'}}\left(t\right)-E_l(\text{Ω}){\chi }_l\left(t\right)-E_o(\text{Ω}){\chi }_o\left(t\right).\end{array}$ (24) Note that the terms in Eq. 24 are re-grouped by low-order (with subscript l) and high-order (with subscript o) terms.

2.2.2

Nodal response matrix equations

The odd-parity moments ${\chi }_l$ and ${\chi }_o$ are defined as continuous across interfaces. Also, requiring Eq. 22 to be stationary with respect to variations in ${\chi }_l$ and ${\chi }_o$ implies that the even-parity flux is continuous across interfaces: ${\phi }_l\left(t\right)={\text{Λ}}^{-1}{\displaystyle \int }\text{d}\text{Ω}{\text{E}}_l^T\left(\text{Ω}\right)\psi \left(\text{Ω}\mathrm{,}t\right).$ (25) By integrating the angle and matrix operations[11] into Eq. 24, the scalar flux can be determined in terms of the response matrix equations, given as $\varphi (t)=V\left(t\right)q^{\text{'}}\left(t\right)-C(t)\left[j^{+}(t)-j^{-}(t)\right].$ (26) Coefficient matrices can also be found in [11, 12]. $j^{+}$ and $j^{-}$ correspond to the interfacial expansion moments of the outgoing and incoming neutron currents, respectively. They are defined as $j^{\pm }\left(t\right)=\frac{1}{4}{\phi }_l(t)\pm \frac{1}{2}{\chi }_l(t)\mathrm{.}$ (27) The outgoing current moments $j^{+}$ satisfy $j^{+}(t)=B(t)q^{\text{'}}(t)+R\left(t\right)j^{-}(t)\mathrm{.}$ (28) In summary, the flowchart of the VNM method is shown in Fig. 3.

Fig. 3Full-screen View

Flowchart of the VNM method

Direct-SCM in the VNM

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 $\overline{\varphi }\left(t\right)={\text{Ξ}}^{-1}{\displaystyle \int }\text{d}x\text{d}yh\left(x\mathrm{,}y\right)g^T\left(x\mathrm{,}y\right){\varphi }_0(t),$ (29) where ${\varphi }_0(t)$ is the segment of $\varphi (t)$ corresponding to the components of the vector of orthogonal polynomials $f_z\left(z\right)$ of order 0, and $\Xi$ is a diagonal matrix composed of the FE areas. $h\left(x\mathrm{,}y\right)$ 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. 4Full-screen View

Flowchart of the transient calculations for the direct-SCM with the VNM

In the k-ω iteration, the secant method [21] can be applied to update the amplitude frequency. The updated formula of the amplitude frequency based on the secant method is $\begin{array}{c}{\omega }_{\text{T}}^{(m+1)}\left(t_n\right)={\omega }_{\text{T}}^{(m)}\left(t_n\right)\\ +\left[{\omega }_{\text{T}}^{(m-1)}\left(t_n\right)-{\omega }_{\text{T}}^{(m)}\left(t_n\right)\right]\frac{1-k_{\text{D}}^{(m)}}{k_{\text{D}}^{(m-1)}-k_{\text{D}}^{(m)}},\end{array}$ (30) where m denotes the iteration step of the k-ω iteration. The normalized element-wise average scalar flux ${\widehat{\varphi }}_{\nu k\mathrm{,}g}\left(t_n\right)$ is $\begin{array}{c}{\widehat{\varphi }}_{\nu k\mathrm{,}g}\left(t_n\right)={\tilde{\varphi }}_{\nu k\mathrm{,}g}\left(t_n\right)\\ \times \frac{P\left(t_{n-1}\right)}{{\displaystyle {\sum }_{\nu }{\displaystyle {\sum }_k{\displaystyle {\sum }_{g=1}^G\nu }}}{\displaystyle {\sum }_{\nu k\mathrm{,}fg}\left(t_n\right)}{A}_{\nu k}{\tilde{\varphi }}_{\nu k\mathrm{,}g}\left(t_n\right)},\end{array}$ (31) in which Aνk is the element area and ${\tilde{\varphi }}_{\nu k\mathrm{,}g}\left(t_n\right)$ 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 ${\omega }_{\nu k\mathrm{,}S\mathrm{,}g}\left(t_n\right)=\frac{1}{\Delta t_n}\ln \left[\frac{{\widehat{\varphi }}_{\nu k\mathrm{,}g}\left(t_n\right)}{{\varphi }_{\nu k\mathrm{,}g}\left(t_{n-1}\right)}\right],$ (32) where Δtn=tn-t_n-1, and ων k,S,g(tn) is the element-wise shape frequency within [t,_n-1,tn].

It is assumed that the amplitude frequency varies linearly within [t_n-1,tn]. ${\omega }_{\text{T}}\left(t\right)={\omega }_{\text{T}}(t_{n-1})+\frac{{\omega }_{\text{T}}(t_n)-{\omega }_{\text{T}}(t_{n-1})}{\Delta t}\left(t-t_{n-1}\right)$ (33) Then, according to Eq. 5, the actual element-wise average scalar flux is $\begin{array}{c}{\varphi }_{\nu k\mathrm{,}g}\left(t_n\right)={\varphi }_{\nu k\mathrm{,}g}\left(t_{n-1}\right)e^{\frac{{\omega }_{\text{T}}\left(t_n\right)+{\omega }_{\text{T}}\left(t_{n-1}\right)}{2}\Delta t_n+{\omega }_{\nu k\mathrm{,}\text{S}\mathrm{,}g}\left(t_n\right)\Delta t_n}\\ ={\widehat{\varphi }}_{\nu k\mathrm{,}g}\left(t_n\right)e^{\frac{{\omega }_{\text{T}}\left(t_n\right)+{\omega }_T\left(t_{n-1}\right)}{2}\Delta t_n}.\end{array}$ (34) 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,tn]. The analytical solution for the element-wise precursor concentration is $\begin{array}{c}{C}_{\nu k\mathrm{,}i}\left(t_n\right)=C_{\nu k\mathrm{,}i}\left(t_{n-1}\right)e^{-{\lambda }_{\nu k\mathrm{,}i}\Delta t_n}\\ +{\beta }_{\nu k\mathrm{,}i}{e}^{-{\lambda }_{\nu k\mathrm{,}i}\Delta t_n}{\displaystyle {\int }_{t_{n-1}}^{t_n}}{Q}_{\nu k}\left(t\right)e^{{\lambda }_{\nu k\mathrm{,}i}t}\text{d}t,\end{array}$ (35) in which the element-wise fission source $Q_{\nu k}\left(t\right)$ is expressed as $Q_{\nu k}\left(t\right)=Q_{\nu k}(t_{n-1})+\frac{Q_{\nu k}(t_n)-Q_{\nu k}(t_{n-1})}{\Delta t_n}\left(t-t_{n-1}\right),$ (36) $Q_{\nu k}\left(t_n\right)={\displaystyle \sum _g\kappa }{\displaystyle {\sum }_{\nu k\mathrm{,}fg}\left(t_n\right)}{\varphi }_{\nu k\mathrm{,}g}\left(t_n\right)\mathrm{,}$ (37) The updated formula for the element-wise precursor frequency is ${\mu }_{\nu k\mathrm{,}i}(t_n)=\{\begin{array}{cc}{\beta }_{\nu k\mathrm{,}i}\frac{Q_{\nu }(t_n)}{C_{\nu k\mathrm{,}i}(t_n)}-{\lambda }_{{}_{\nu k\mathrm{,}i}}& C_{\nu k\mathrm{,}i}(t_n)\ne 0\\ 0& C_{\nu k\mathrm{,}i}(t_n)=0\end{array}$ (38)

PCQ-SCM in the VNM

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. The error in the neutron flux resulting from the coarse time step is then corrected using the amplitude function.

Fig. 5Full-screen View

Illustration of the coarse and fine time steps in the PCQ-SCM

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.

2.4.2

EPKE

The EPKE is derived from the transient transport equation[28]. Similarly, the neutron flux is factorized into the shape function and amplitude function, as follows: ${\psi }_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right)=T\left(t\right){\tilde{\psi }}_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right),$ (44) where $T\left(t\right)$ is the amplitude function, and ${\tilde{\psi }}_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right)$ is the shape function. Similarly, to make the decomposition unique, the following constraint is introduced: ${\displaystyle \sum _{g=1}^G{\displaystyle {\int }_V}}\text{d}r{\displaystyle \int }\text{d}\text{Ω}{\psi }_g^{*}\left(r\mathrm{,}\text{Ω}\mathrm{,}{t}_0\right)\frac{1}{v_g\left(r\right)}{\tilde{\psi }}_g\left(r\mathrm{,}\text{Ω}\mathrm{,}t\right)=T\left(t_0\right),$ (45) where ${\psi }_g^{*}\left(r\mathrm{,}\text{Ω}\mathrm{,}{t}_0\right)$ 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: $\frac{\text{d}}{\text{d}t}T\left(t\right)=\frac{\rho \left(t\right)-\overline{\beta }\left(t\right)}{\Lambda \left(t\right)}T\left(t\right)+{\displaystyle \sum _{i=1}^I{\overline{\lambda }}_i}\left(t\right){\zeta }_i\left(t\right),$ (46) $\frac{\text{d}}{\text{d}t}{\zeta }_i\left(t\right)=\frac{{\overline{\beta }}_i\left(t\right)}{\Lambda \left(t\right)}T\left(t\right)-{\overline{\lambda }}_i\left(t\right){\zeta }_i\left(t\right),$ (47) where $\rho \left(t\right)$ is the reactivity, $\overline{\beta }\left(t\right)$ and ${\overline{\beta }}_i\left(t\right)$ are effective fractions of delayed neutrons, ${\overline{\lambda }}_i\left(t\right)$ is the effective decay constant, $\Lambda \left(t\right)$ is the prompt neutron lifetime, and ${\zeta }_i\left(t\right)$ is the reduced precursor concentration. The detailed formulation of the EPKE is available in [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].

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] 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 $\tilde{\varphi }$ and ϕ^† are $\begin{array}{c}F\left(t\right)={\displaystyle \sum _{g=1}^G{\displaystyle \int }\text{d}r}\left\{\begin{array}{c}\frac{{\chi }_g\left(r\right)}{k_{\text{eff}}}{\varphi }_g^{*}\left(r\mathrm{,}{t}_0\right)\\ \times {\displaystyle {\sum }_{g^{\text{'}}}^G\nu }{\displaystyle {\sum }_{f{g}^{\text{'}}}\left(r\mathrm{,}t\right)}{\tilde{\varphi }}_{g^{\text{'}}}\left(r\mathrm{,}t\right)\end{array}\right\}\\ ={\displaystyle \sum _{g=1}^G{\displaystyle \sum _{\nu }{\displaystyle \sum _k\left\{\begin{array}{c}\frac{{\chi }_{\nu k\mathrm{,}g}}{k_{\text{eff}}}{A}_{\nu k}{\varphi }_{\nu k\mathrm{,}g}^{*}\left(t_0\right)\\ \times {\displaystyle {\sum }_{g^{\prime }}^G\nu }{\displaystyle {\sum }_{\nu k\mathrm{,}f{g}^{\text{'}}}\left(t\right)}{\tilde{\varphi }}_{\nu k\mathrm{,}{g}^{\text{'}}}\left(t\right)\end{array}\right\}}}},\end{array}$ (48) $\begin{array}{c}{\overline{\beta }}_i\left(t\right)=\frac{{\displaystyle \sum _{g=1}^G{\displaystyle \int }\text{d}r}\left\{\begin{array}{c}\frac{{\chi }_{ig}\left(\text{r}\right)}{k_{\text{eff}}}{\varphi }_g^{*}\left(r\mathrm{,}{t}_0\right){\beta }_i\left(r\right)\\ \times {\displaystyle {\sum }_{g^{\text{'}}}^G\nu }{\displaystyle {\sum }_{f{g}^{\text{'}}}\left(r\mathrm{,}t\right)}{\tilde{\varphi }}_{g^{\text{'}}}\left(r\mathrm{,}t\right)\end{array}\right\}}{F\left(t\right)}\\ =\frac{{\displaystyle {\sum }_{g=1}^G{\displaystyle {\sum }_{\nu }{\displaystyle {\sum }_k\left\{\begin{array}{c}\frac{{\chi }_{\nu k\mathrm{,}ig}}{k_{\text{eff}}}{A}_{\nu k}{\varphi }_{\nu k\mathrm{,}g}^{*}\left(t_0\right){\beta }_{\nu k\mathrm{,}i}\\ \times {\displaystyle {\sum }_{g^{\text{'}}}^G\nu }{\displaystyle {\sum }_{\nu k\mathrm{,}f{g}^{\text{'}}}\left(t\right)}{\tilde{\varphi }}_{\nu k\mathrm{,}{g}^{\text{'}}}\left(t\right)\end{array}\right\}}}}}{F\left(t\right)},\end{array}$ (49) $\Lambda \left(t\right)=\frac{{\displaystyle {\sum }_{g=1}^G{\displaystyle \int }\text{d}r}\frac{{\varphi }_g^{*}\left(r\mathrm{,}{t}_0\right){\tilde{\varphi }}_g\left(r\mathrm{,}t\right)}{v_g\left(r\mathrm{,}t\right)}}{F\left(t\right)}=\frac{T\left(t_0\right)}{F\left(t\right)}.$ (50) $\begin{array}{c}{\overline{\lambda }}_i\left(t\right)=\frac{{\displaystyle {\sum }_{g=1}^G{\displaystyle \int }\text{d}r}{\varphi }_g^{*}\left(r\mathrm{,}{t}_0\right){\chi }_{ig}\left(r\right){\lambda }_i\left(r\right)C_i\left(r\mathrm{,}t\right)}{{\displaystyle {\sum }_{g=1}^G{\displaystyle \int }\text{d}r}{\varphi }_g^{*}\left(r\mathrm{,}{t}_0\right){\chi }_{ig}\left(r\right)C_i\left(r\mathrm{,}t\right)}\\ =\frac{{\displaystyle {\sum }_{g=1}^G{\displaystyle {\sum }_{\nu }{\displaystyle {\sum }_k{A}_{\nu k}{\varphi }_{\nu k\mathrm{,}g}^{*}}}}\left(t_0\right){\chi }_{\nu k\mathrm{,}ig}{\lambda }_{\nu k\mathrm{,}i}{C}_{\nu k\mathrm{,}i}\left(t\right)}{{\displaystyle {\sum }_{g=1}^G{\displaystyle {\sum }_{\nu }{\displaystyle {\sum }_k{A}_{\nu k}{\varphi }_{\nu k\mathrm{,}g}^{*}}}}\left(t_0\right){\chi }_{\nu k\mathrm{,}ig}{C}_{\nu k\mathrm{,}i}\left(t\right)}.\end{array}$ (51) Previous work showed that the KPs calculated using the scalar flux agreed well with those calculated using the angular flux [5]. However, using the scalar flux as the weighting function could result in significant reactivity deviations [5, 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]. The different reactivities resulting from the scalar flux were mainly owing to the strong anisotropy of the leakage operator $\text{Ω}\cdot \nabla$. 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] as $\rho \left(t\right)=\frac{k_{\text{eff}}\left(t\right)-1}{k_{\text{eff}}\left(t\right)},$ (52) where $k_{\text{eff}}\left(t\right)$ is the dynamic eigenvalue, which is defined as the eigenvalue of Eq. 13, obtained by setting frequencies μ and ω to zeros. Since $k_{\text{eff}}\left(t\right)$ 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. 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. 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].

Fig. 6Full-screen View

Reactivity evaluation using different methods, for the C5G7-TD problems

2.4.3

Overall PCQ-SCM scheme

In the PCQ-SCM solver of VITAS, the initial adjoint scalar flux ${\varphi }_g^{*}\left(\text{r}\mathrm{,}{t}_0\right)$ and the forward scalar flux ${\varphi }_g\left(\text{r}\mathrm{,}{t}_0\right)$ are obtained by solving the initial EVPs. The initial amplitude function T(t₀) is determined using P(t₀). Then, in each coarse time step [t_n-1,tn], the TEVP and EPKE are solved as follows:

1. Solve the predicted scalar flux ${\varphi }_{p\mathrm{,}g}\left(\text{r}\mathrm{,}{t}_n\right)$, flux amplitude Pp(tn), and precursor concentrations Cp,i(𝒓,t) using the TEVP in Eq. 13 with k-ω iteration at tn. In the first iteration, the frequencies μ and ω are set to zero. Then, the eigenvalue obtained in the first iteration is the dynamic eigenvalue $k_{\text{eff}}\left(t_n\right)$ required in Eq. 52.

2. Determine the reactivity ρ(tn) using Eq. 52.

3. Calculate the predicted shape function ${\tilde{\varphi }}_{p\mathrm{,}g}\left(\text{r}\mathrm{,}{t}_n\right)$ with the constraint in Eq. 45 as $\begin{array}{c}{\tilde{\varphi }}_{p\mathrm{,}g}\left(r\mathrm{,}{t}_n\right)={\varphi }_{p\mathrm{,}g}\left(r\mathrm{,}{t}_n\right)\\ \times \frac{T\left(t_0\right)}{{\displaystyle {\sum }_{g=1}^G{\displaystyle {\int }_V}\text{d}r}{\displaystyle \int }\text{d}\text{Ω}{\varphi }_g^{*}\left(r\mathrm{,}{t}_0\right)\frac{1}{v_g\left(r\right)}{\varphi }_{p\mathrm{,}g}\left(r\mathrm{,}{t}_n\right)}.\end{array}$ (53) 4. Calculate the KPs at tn with Eqns. 48-51, and solve the EPKE to obtain the amplitude function T(tn) at the fine time step using the interpolated KPs from the values at t_n-1 and tn.

5. Determine the corrected flux amplitude Pc(tn) and the scalar flux ${\varphi }_{c\mathrm{,}g}\left(r\mathrm{,}{t}_n\right)$, as follows: $P_c\left(t_n\right)=\frac{T\left(t_n\right)P\left(t_0\right)}{T\left(t_0\right)},$ (54) $\begin{array}{c}{\varphi }_{c\mathrm{,}g}\left(r\mathrm{,}{t}_n\right)=\frac{P_c\left(t_n\right)}{P_p\left(t_n\right)}{\varphi }_{p\mathrm{,}g}\left(r\mathrm{,}{t}_n\right)\\ =P_c\left(t_n\right){\widehat{\varphi }}_{p\mathrm{,}g}\left(r\mathrm{,}{t}_n\right).\end{array}$ (55) 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. 7Full-screen View

Flowchart of the transient calculations with the PCQ-SCM for the VNM

In addition, in the PCQ-SCM solver, the amplitude frequency ωT is approximated as a constant in [t_n-1,tn], 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 $\begin{array}{c}{\varphi }_{\nu k\mathrm{,}g}\left(t_n\right)={\varphi }_{\nu k\mathrm{,}g}\left(t_{n-1}\right)e^{{\omega }_T\left(t_n\right)\Delta t_n+{\omega }_{\nu k\mathrm{,}\text{S}\mathrm{,}g}\left(t_n\right)\Delta t_n}\\ ={\widehat{\varphi }}_{\nu k\mathrm{,}g}\left(t_n\right)e^{{\omega }_T\left(t_n\right)\Delta t_n}.\end{array}$ (56)

Geometry description

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. 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. The height of the 3D configuration is 171.36 cm while that of the active region is 128.52 cm.

Fig. 8Full-screen View

Geometry and composition of the C5G7-TD benchmark

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. 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. 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

The C5G7-TD benchmark adopts multi-group macroscopic cross-sections and KPs for eight materials: UO₂ fuel, 7.0% MOX fuel, 4.3% MOX fuel, 8.7% MOX fuel, guide tube, fission chamber, moderator, and control rod[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 ${\displaystyle \sum {}_{\text{s}g{g}^{\prime }}}$. The KPs are provided in an 8-delayed-group structure including neutron velocities vg, 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] 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 $\begin{array}{c}y=0.3867848x+0.1707878x^2+0.9887881x^3\\ -5.9535775x^4+25.1805898x^5-63.1841252x^6\\ +98.3898802x^7-92.4982596x^8+48.0373368x^9\\ -10.5182051x^{10}\end{array}$ (57) Figure 9 shows the relationships between the original and effective volume fractions, for Eq. 57 and linear assumption, respectively.

Fig. 9Full-screen View

Relationship between original and effective rodded volume fractions

Transient description

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, 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.

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. 10Full-screen View

Fractional moderator density change in the TD3 exercise

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 depicts the control-rod banks and the relative insertion length to the active-core length, for each problem.

Fig. 11Full-screen View

Relative depths of control-bank movements, in the TD4 exercise

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.

Steady-state results

For steady-state calculations, Table 3 compares the eigenvalue results for the 2D and 3D initial eigenvalue problems calculated using different solvers[5, 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.

The normalized pin-wise fission rate distributions for the 2D and 3D problems are shown in Fig. 13 and Fig. 14, respectively. The normalization was implemented by scaling the initial total fission rate to unity.

Fig. 13Full-screen View

(Color online) Pin-wise fission rate distribution in the steady state for the 2D cases

Fig. 14Full-screen View

(Color online) Pin-wise fission rate distribution in the steady state for the 3D cases

Verification of the direct-SCM solver

To verify the direct-SCM solver in VITAS, the VITAS results were compared with the results of MPACT[4] and PANDAS-MOC[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 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 $\epsilon \left(t\right)=\frac{P\left(t\right)-P^{\star }\left(t\right)}{P\left(t_0\right)}\times 100\%,$ (58) where the superscript ^* denotes the reference solution. 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. 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.

Fig. 15Full-screen View

Fractional fission rate and percentage difference, for different time-step sizes in the TD3-4 case

To further quantify the effect of the time-step size, 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, 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]. 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. $\text{RMS}=\sqrt{\frac{{{\displaystyle {\sum }_{n=0}^N\left[P\left(t_n\right)-P^{\star }\left(t_n\right)\right]}}^2}{N}}$ (59)

4.2.1

Results of 2D studies

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] and PANDAS-MOC[6]. Relative 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. 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 and 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 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, the direct-SCM solver’s results for the TD3 case were also consistent with those obtained using the other codes.

Fig. 16Full-screen View

Fractional fission rates, for the TD1 and TD2 exercises

Fig. 17Full-screen View

Comparison of fractional fission rate changes of VITAS, MPACT and PANDAS-MOC for the TD1 exercise

Fig. 18Full-screen View

Comparison of fractional fission rate changes of VITAS, MPACT and PANDAS-MOC for the TD2 exercise

Fig. 19Full-screen View

Fractional fission rate for the TD3 exercise

Fig. 20Full-screen View

Comparison of fractional fission rate changes of VITAS, MPACT and PANDAS-MOC for the TD3 exercise

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 and 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.

Table 5Full-screen View

Maximal difference and RMS error, for the 2D results

Problem	v.s. MPACT Max. Diff. (%)	v.s. MPACT RMS Err. (%)	v.s. PANDAS Max. Diff. (%)	v.s. PANDAS RMS Err. (%)
1-1	-0.13	0.08	0.14	0.06
1-2	-0.16	0.08	-0.16	0.08
1-3	-0.14	0.07	-0.14	0.06
1-4	-0.17	0.07	0.17	0.06
1-5	-0.17	0.08	0.27	0.08
2-1	0.91	0.10	0.91	0.10
2-2	-0.21	0.11	-0.23	0.15
2-3	-0.12	0.05	-0.11	0.04
3-1	0.25	0.11	0.29	0.16
3-2	0.45	0.11	0.56	0.16
3-3	0.72	0.11	0.72	0.15
3-4	0.97	0.11	0.97	0.14

Show more

4.2.2

Results of 3D studies

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. 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 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, 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.

Table 6Full-screen View

Maximal difference and RMS error, for the TD4 case

Problem	v.s. MPACT Max. Diff. (%)	v.s. MPACT RMS Err. (%)	v.s. PANDAS Max. Diff. (%)	v.s. PANDAS RMS Err. (%)
4-1	-0.97	0.28	1.57	0.69
4-2	-0.77	0.28	0.68	0.28
4-3	-0.77	0.29	0.68	0.31
4-4	-0.75	0.34	0.74	0.25
4-5	-0.97	0.39	1.57	0.62

Show more

Fig. 21Full-screen View

Fractional fission rate, for the TD4 exercise

Fig. 22Full-screen View

Comparison of fractional fission rate changes of VITAS, MPACT and PANDAS-MOC for the TD4 exercise

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. 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 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, 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.

Table 7Full-screen View

Maximal difference and RMS error, for the TD5 case

Problem	v.s. MPACT Max. Diff.(%)	v.s. MPACT RMS Err.(%)	v.s. PANDAS Max. Diff.(%)	v.s. PANDAS RMS Err.(%)
5-1	0.69	0.43	0.74	0.52
5-2	0.87	0.45	0.85	0.54
5-3	0.85	0.38	0.78	0.47
5-4	0.33	0.10	0.34	0.12

Show more

Fig. 23Full-screen View

Fractional fission rate, for the TD5 exercise

Fig. 12Full-screen View

Fractional moderator density change in the TD5 exercise

Fig. 24Full-screen View

Comparison of fractional fission rate changes of VITAS, MPACT and PANDAS-MOC for the TD5 exercise

4.2.3

Step-perturbation results

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]. 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 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.

Fig. 25Full-screen View

Variable time-step sizes of VITAS, for the TD0 exercise

The fractional core fission rates for the TD0 exercise are shown in 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 compares the VITAS and MPACT solutions, while 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.

Table 8Full-screen View

Maximal difference and RMS error, for the TD0 case

Problem	v.s. MPACT Max. Diff. (%)	v.s. MPACT RMS Err. (%)
0-1	-0.12	0.05
0-2	-0.20	0.11
0-3	-0.10	0.05
0-4	-0.11	0.05
0-5	-0.11	0.05

Show more

Fig. 26Full-screen View

Fractional fission rate, for the TD0 exercise

Fig. 27Full-screen View

Comparison of fractional fission rate changes of VITAS and MPACT for the TD0 exercise

Performance of the PCQ-SCM solver

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 shows that the reactivity evaluated using VITAS was consistent with that evaluated using MPACT. 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, 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, 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, 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.

Fig. 28Full-screen View

Comparison of the reactivity history for the TD0-5 problem

Fig. 29Full-screen View

Comparison of the direct-SCM and PCQ-SCM solutions for the TD0-5 problem

In addition, 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 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 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 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. This was less than the direct solution of 10 ms, as shown in Table 4. In addition, 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.

Fig. 30Full-screen View

Comparison of the reactivity history for the TD3-4 problem

Fig. 31Full-screen View

Comparison of the direct-SCM and PCQ-SCM solutions for the TD3-4 problem

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 shows the reactivity history for TD5-1, as evaluated by VITAS, MPACT, and PANDAS-MOC. Notably, the reactivity in Fig. 32 was defined as $\rho /\overline{\beta }$ with the units of $, whereas it was defined as ρ with the units of pcm in Fig. 32. Fig. 32 reveals that, compared with the reactivity of MPACT and PANDAS-MOC, VITAS underestimated $\rho /\overline{\beta }$. However, in Fig. 32, the reactivity ρ evaluated by VITAS was consistent with that estimated by PANDAS-MOC. Consequently, the deviations in Fig. 32 were attributed to the discrepancies of the effective delayed neutron fractions $\overline{\beta }$ computed by different codes.

Fig. 32Full-screen View

Comparison of the reactivity history for the TD5-1 problem

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, 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. 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, 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].

Fig. 33Full-screen View

Comparison of the direct-SCM and PCQ-SCM solutions for the TD5-1 problem

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 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 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.

Fig. 34Full-screen View

(Color online) Runtime analysis of the PCQ-SCM solutions for the TD5-1 problem