State change modal method for numerical simulation of dynamic processes in a nuclear reactor

TL;DR

This paper introduces a state change modal method for simulating dynamic processes in nuclear reactors, utilizing eigenfunction expansions and spectral analysis to improve modeling accuracy, demonstrated through simulations of a VVER-1000 reactor.

Contribution

The paper presents a novel numerical-analytical modal method based on dominant eigenvalues for dynamic reactor modeling, accounting for complex eigenvalues in a two-group diffusion approximation.

Findings

Effective simulation of reactor dynamics using the modal method.

Inclusion of complex eigenvalues improves model accuracy.

Application to VVER-1000 reactor demonstrates method viability.

Abstract

Modeling of dynamic processes in nuclear reactors is carried out, mainly, on the basis of the multigroup diffusion approximation for the neutron flux. The basic model includes a multidimensional set of coupled parabolic equations and ordinary differential equations. Dynamic processes are modeled by a successive change of the reactor states, which are characterized by given coefficients of the equations. In the modal method, the approximate solution is represented as an expansion on the first eigenfunctions of some spectral problem. The numerical-analytical method is based on the use of dominant time-eigenvalues of a multigroup diffusion model taking into account delayed neutrons. Numerical simulations of the dynamic process were performed in the framework of the two-group approximation for the VVER-1000 reactor test model. The last is characterized by the fact that some eigenvalues are complex.

Tables

Table 1. Table 1: Diffusion neutronics constants for VVER-1000

Material	1	2	3	4	5
$D_1$	1.38320e-0	1.38299e-0	1.39522e-0	1.39446e-0	1.39506e-0
$D_2$	3.86277e-1	3.89403e-1	3.86225e-1	3.87723e-1	3.84492e-1
${\Sigma }_1+{\Sigma }_{s,1\to 2}$	2.48836e-2	2.62865e-2	2.45662e-2	2.60117e-2	2.46141e-2
${\Sigma }_2$	6.73049e-2	8.10328e-2	8.44801e-1	9.89671e-2	8.93878e-2
${\Sigma }_{s,1\to 2}$	1.64977e-2	1.47315e-2	1.56219e-2	1.40185e-2	1.54981e-2
$\nu {\Sigma }_{f1}$	4.81619e-3	4.66953e-3	6.04889e-3	5.91507e-3	6.40256e-3
$\nu {\Sigma }_{f2}$	8.46154e-2	8.52264e-2	1.19428e-1	1.20497e-1	1.29281e-1

Table 2. Table 2: Eigenvalues α n = λ n ( α ) , n = 1 , 2 , … , 5 formulae-sequence subscript 𝛼 𝑛 superscript subscript 𝜆 𝑛 𝛼 𝑛 1 2 … 5 \alpha_{n}=\lambda_{n}^{(\alpha)},\ n=1,2,...,5

$p$	$\kappa$	${\alpha }_1$	${\alpha }_2,{\alpha }_3$	${\alpha }_4,{\alpha }_5$
	6	-0.22557	0.04241 $\mp$ 3.08808e-06$i$	0.06588 $\mp$ 4.80449e-07$i$
1	24	-0.82690	0.03777 $\mp$ 5.37884e-06$i$	0.06489 $\mp$ 1.37315e-06$i$
	96	-1.74998	0.03619 $\mp$ 5.69002e-06$i$	0.06456 $\mp$ 1.40299e-06$i$
	6	-2.10154	0.03592 $\mp$ 4.96474e-06$i$	0.06452 $\mp$ 1.21320e-06$i$
2	24	-2.46601	0.03562 $\mp$ 5.78277e-06$i$	0.06445 $\mp$ 1.40897e-06$i$
	96	-2.50375	0.03559 $\mp$ 5.80693e-06$i$	0.06444 $\mp$ 1.41324e-06$i$
	6	-2.47975	0.03561 $\mp$ 5.83718e-06$i$	0.06445 $\mp$ 1.41869e-06$i$
3	24	-2.50294	0.03559 $\mp$ 5.80783e-06$i$	0.06444 $\mp$ 1.41341e-06$i$
	96	-2.51280	0.03558 $\mp$ 5.80954e-06$i$	0.06444 $\mp$ 1.41362e-06$i$

Table 3. Table 3: Eigenvalues α n = λ n ( α ) , n = 6 , 7 , … , 10 formulae-sequence subscript 𝛼 𝑛 superscript subscript 𝜆 𝑛 𝛼 𝑛 6 7 … 10 \alpha_{n}=\lambda_{n}^{(\alpha)},\ n=6,7,...,10

$p$	$\kappa$	${\alpha }_6$	${\alpha }_7$	${\alpha }_8$	${\alpha }_9,{\alpha }_{10}$
	6	0.07107	0.07214	0.07323	0.07397 $\mp$ 2.04990e-08$i$
1	24	0.07050	0.07167	0.07283	0.07362 $\mp$ 3.65907e-08$i$
	96	0.07033	0.07152	0.07269	0.07351 $\mp$ 3.91936e-08$i$
	6	0.07030	0.07151	0.07268	0.07349 $\mp$ 3.69824e-08$i$
2	24	0.07027	0.07147	0.07265	0.07347 $\mp$ 4.03121e-08$i$
	96	0.07026	0.07147	0.07265	0.07347 $\mp$ 4.02324e-08$i$
	6	0.07027	0.07147	0.07265	0.07347 $\mp$ 4.02573e-08$i$
3	24	0.07026	0.07147	0.07265	0.07347 $\mp$ 4.02248e-08$i$
	96	0.07026	0.07147	0.07265	0.07347 $\mp$ 4.02332e-08$i$

Table 4. Table 4: Eigenvalues α n = λ n ( α ) , n = 1 , 2 , … , 10 formulae-sequence subscript 𝛼 𝑛 superscript subscript 𝜆 𝑛 𝛼 𝑛 1 2 … 10 \alpha_{n}=\lambda_{n}^{(\alpha)},\ n=1,2,...,10 for direct and adjoint problems

$n$	${\alpha }_n$ for problem (14)	${\alpha }_n$ for problem (16)
1	-2.51280117966	-2.51280117972
2,3	0.0355815000364 $\mp$ 5.80954455861e-06	0.0355815000365 $\mp$ 5.80954421646e-06
4,5	0.0644427013767 $\mp$ 1.41362187449e-06	0.0644427013767 $\mp$ 1.41362190730e-06
6	0.0702618501639	0.0702618501639
7	0.0714652882224	0.0714652882164
8	0.0726456060606	0.0726456060606
9,10	0.0734708921578 $\mp$ 4.02332269037e-08	0.0734708921578 $\mp$ 4.02332146248e-08

Table 5. Table 5: Scalar product ( ϕ 1 ( n ) , ϕ 1 ( m ) ) , n , m = 1 , 2 , … , 10 formulae-sequence superscript subscript italic-ϕ 1 𝑛 superscript subscript italic-ϕ 1 𝑚 𝑛 𝑚 1 2 … 10 (\phi_{1}^{(n)},\phi_{1}^{(m)}),\ n,m=1,2,...,10

$n$\$m$	1	2	3	4	5	6	7	8	9	10
1	1.0e-00	1.3e-08	2.2e-08	-3.8e-08	9.8e-09	-1.8e-09	1.0e-02	-3.2e-09	-2.2e-08	1.6e-09
2	1.3e-08	1.0e-00	-1.6e-08	-1.6e-08	1.4e-08	4.1e-08	1.2e-09	-2.0e-07	-3.1e-03	7.5e-03
3	2.2e-08	-1.6e-08	1.0e-00	-9.8e-09	-1.1e-08	-1.8e-08	1.1e-08	-3.3e-08	-7.5e-03	-3.1e-03
4	-3.8e-08	-1.6e-08	-9.8e-09	1.0e-00	-3.9e-10	-1.1e-08	1.4e-08	4.0e-09	3.0e-09	-1.1e-08
5	9.8e-09	1.4e-08	-1.1e-08	-3.9e-10	1.0e-00	2.9e-09	-1.6e-08	-1.9e-08	6.3e-09	6.3e-09
6	-1.8e-09	4.1e-08	-1.8e-08	-1.1e-08	2.9e-09	1.0e-00	-4.2e-09	-5.6e-03	4.1e-08	-1.2e-07
7	1.0e-02	1.2e-09	1.1e-08	1.4e-08	-1.6e-08	-4.2e-09	1.0e-00	-2.1e-09	-1.8e-08	8.0e-09
8	-3.2e-09	-2.0e-07	-3.3e-08	4.0e-09	-1.9e-08	-5.6e-03	-2.1e-09	1.0e-00	-5.2e-08	2.3e-07
9	-2.2e-08	-3.1e-03	-7.5e-03	3.0e-09	6.3e-09	4.1e-08	-1.8e-08	-5.2e-08	1.0e-00	-5.5e-07
10	1.6e-09	7.5e-03	-3.1e-03	-1.1e-08	6.3e-09	-1.2e-07	8.0e-09	2.3e-07	-5.5e-07	1.0e-00

Table 6. Table 6: Subcritical state: α n = λ n ( α ) , n = 1 , 2 , … , 5 formulae-sequence subscript 𝛼 𝑛 superscript subscript 𝜆 𝑛 𝛼 𝑛 1 2 … 5 \alpha_{n}=\lambda_{n}^{(\alpha)},\ n=1,2,...,5

$p$	$\kappa$	${\alpha }_1$	${\alpha }_2,{\alpha }_3$	${\alpha }_4,{\alpha }_5$
	6	0.03602	0.05760 $\mp$ 1.49652e-06$i$	0.06890 $\mp$ 4.92606e-07$i$
1	24	0.02656	0.05502 $\mp$ 2.06007e-06$i$	0.06804 $\mp$ 1.01253e-06$i$
	96	0.02276	0.05411 $\mp$ 2.16813e-06$i$	0.06774 $\mp$ 1.03843e-06$i$
	6	0.02250	0.05404 $\mp$ 1.81823e-06$i$	0.06772 $\mp$ 8.73562e-07$i$
2	24	0.02144	0.05380 $\mp$ 2.19400e-06$i$	0.06765 $\mp$ 1.04253e-06$i$
	96	0.02125	0.05376 $\mp$ 2.20812e-06$i$	0.06763 $\mp$ 1.04715e-06$i$
	6	0.02139	0.05379 $\mp$ 2.22579e-06$i$	0.06764 $\mp$ 1.05369e-06$i$
3	24	0.02124	0.05376 $\mp$ 2.20883e-06$i$	0.06763 $\mp$ 1.04736e-06$i$
	96	0.02122	0.05376 $\mp$ 2.20951e-06$i$	0.06763 $\mp$ 1.04756e-06$i$

Table 7. Table 7: Subcritical state: α n = λ n ( α ) , n = 6 , 7 , … , 10 formulae-sequence subscript 𝛼 𝑛 superscript subscript 𝜆 𝑛 𝛼 𝑛 6 7 … 10 \alpha_{n}=\lambda_{n}^{(\alpha)},\ n=6,7,...,10

$p$	$\kappa$	${\alpha }_6$	${\alpha }_7$	${\alpha }_8$	${\alpha }_9,{\alpha }_{10}$
	6	0.07276	0.07363	0.07369	0.07466 $\mp$ 2.47162e-08$i$
1	24	0.07222	0.07316	0.07329	0.07429 $\mp$ 1.08814e-08$i$
	96	0.07204	0.07301	0.07316	0.07417 $\mp$ 1.38093e-08$i$
	6	0.07203	0.07300	0.07315	0.07416 $\mp$ 1.26708e-08$i$
2	24	0.07199	0.07296	0.07312	0.07413 $\mp$ 1.50527e-08$i$
	96	0.07198	0.07295	0.07312	0.07413 $\mp$ 1.49196e-08$i$
	6	0.07198	0.07295	0.07312	0.07413 $\mp$ 1.52256e-08$i$
3	24	0.07198	0.07295	0.07312	0.07413 $\mp$ 1.49141e-08$i$
	96	0.07198	0.07295	0.07311	0.07413 $\mp$ 1.49178e-08$i$

Table 8. Table 8: Subcritical asymmetric state: α n = λ n ( α ) , n = 1 , 2 , … , 5 formulae-sequence subscript 𝛼 𝑛 superscript subscript 𝜆 𝑛 𝛼 𝑛 1 2 … 5 \alpha_{n}=\lambda_{n}^{(\alpha)},\ n=1,2,...,5

$p$	$\kappa$	${\alpha }_1$	${\alpha }_2$	${\alpha }_3$	${\alpha }_4$	${\alpha }_5$
	6	0.03347	0.05728	0.05788	0.06884	0.06889
1	24	0.02333	0.05467	0.05528	0.06797	0.06802
	96	0.01925	0.05374	0.05436	0.06768	0.06772
	6	0.01894	0.05367	0.05429	0.06765	0.06770
2	24	0.01782	0.05343	0.05405	0.06758	0.06762
	96	0.01763	0.05339	0.05401	0.06757	0.06761
	6	0.01777	0.05342	0.05404	0.06758	0.06762
3	24	0.01762	0.05339	0.05400	0.06757	0.06761
	96	0.01760	0.05338	0.05400	0.06757	0.06761

Table 9. Table 9: Subcritical asymmetric state: α n = λ n ( α ) , n = 6 , 7 , … , 10 formulae-sequence subscript 𝛼 𝑛 superscript subscript 𝜆 𝑛 𝛼 𝑛 6 7 … 10 \alpha_{n}=\lambda_{n}^{(\alpha)},\ n=6,7,...,10

$p$	$\kappa$	${\alpha }_6$	${\alpha }_7$	${\alpha }_8$	${\alpha }_9$	${\alpha }_{10}$
	6	0.07274	0.07355	0.07369	0.07464	0.07468
1	24	0.07220	0.07309	0.07329	0.07427	0.07430
	96	0.07202	0.07294	0.07316	0.07415	0.07419
	6	0.07201	0.07293	0.07314	0.07414	0.07417
2	24	0.07197	0.07289	0.07311	0.07412	0.07415
	96	0.07196	0.07288	0.07311	0.07411	0.07414
	6	0.07196	0.07289	0.07311	0.07411	0.07415
3	24	0.07196	0.07288	0.07311	0.07411	0.07414
	96	0.07196	0.07288	0.07311	0.07411	0.07414

Equations

$$\begin{array}{rl}\frac{1}{v_g}\frac{\partial {\varphi }_g}{\partial t}-& \nabla \cdot D_g\nabla {\varphi }_g+{\Sigma }_{rg}{\varphi }_g-\sum _{g\ne g^{\prime }=1}^G{\Sigma }_{s,g^{\prime }\to g}{\varphi }_{g^{\prime }}\\ =& (1-\beta ){\chi }_g\sum _{g^{\prime }=1}^G\nu {\Sigma }_{f{g}^{\prime }}{\varphi }_{g^{\prime }}+{\tilde{\chi }}_g\sum _{m=1}^M{\lambda }_m{c}_m,g=1,2,...,G.\end{array}$$

$$\frac{\partial c_m}{\partial t}+{\lambda }_m{c}_m={\beta }_m\sum _{g=1}^G\nu {\Sigma }_{fg}{\varphi }_g,m=1,2,...,M,$$

$$\beta =\sum _{m=1}^M{\beta }_m.$$

$$D_g\frac{\partial {\varphi }_g}{\partial n}+{\gamma }_g{\varphi }_g=0,g=1,2,...,G,$$

$${\varphi }_g(\mathbf{x},0)={\varphi }_g^0(\mathbf{x}),c_m(\mathbf{x},0)=c_m^0(\mathbf{x}).$$

$$V$$

$$D$$

$$S$$

$$R$$

$$B$$

$$\Lambda$$

$$Q$$

$$g,g^{\prime }$$

$${\delta }_{g{g}^{\prime }}=\{\begin{array}{cc}1,& g=g^{\prime },\\ 0,& g\ne g^{\prime },\end{array}$$

$$\begin{array}{rl}V(t)\frac{d\mathbf{\varphi }}{dt}+(D(t)+S(t))\mathbf{\varphi }& =R(t)\mathbf{\varphi }+B(t)\mathbf{c},\\ \frac{d\mathbf{c}}{dt}+\Lambda (t)\mathbf{c}& =Q(t)\mathbf{\varphi }.\end{array}$$

$$\mathbf{\varphi }(0)={\mathbf{\varphi }}^0,\mathbf{c}(0)={\mathbf{c}}^0,$$

$$V(t)=V(t_s),D(t)=D(t_s),S(t)=S(t_s),R(t)=R(t_s),B(t)=B(t_s),$$

$$\Lambda (t)=\Lambda (t_s),Q(t)=Q(t_s),t_{s-1}<t\le t_s,s=1,2,....$$

$$\begin{array}{rl}V(t_s)\frac{d\mathbf{\varphi }}{dt}+(D(t_s)+S(t_s))\mathbf{\varphi }& =R(t_s)\mathbf{\varphi }+B(t_s)\mathbf{c},\\ \frac{d\mathbf{c}}{dt}+\Lambda (t_s)\mathbf{c}& =Q(t_s)\mathbf{\varphi },t_{s-1}<t\le t_s,\end{array}$$

$$\mathbf{\varphi }(t_{s-1})={\mathbf{\varphi }}^s,\mathbf{c}(t_{s-1})={\mathbf{c}}^s.$$

$$\mathbf{B}\frac{d\mathbf{u}}{dt}+\mathbf{A}\mathbf{u}=0,t_{s-1}<t\le t_s,$$

$$\mathbf{A}=\left(\begin{array}{cc}D(t_s)+S(t_s)-R(t_s)& -B(t_s)\\ -Q(t_s)& \Lambda (t_s)\end{array}\right),\mathbf{B}=\left(\begin{array}{cc}V(t_s)& 0\\ 0& I\end{array}\right),$$

$$\mathbf{u}(t_{s-1})={\mathbf{u}}^s.$$

$${\mathbf{B}}_h\frac{d{\mathbf{u}}_h}{dt}+{\mathbf{A}}_h{\mathbf{u}}_h=0,t_{s-1}<t\le t_s,$$

$${\mathbf{u}}_h(t_{s-1})={\mathbf{u}}_h^s,$$

$${\mathbf{u}}_N(\mathbf{x},t)=\sum _{n=1}^N{a}_n(t){\mathbf{w}}_n(\mathbf{x}),$$

$${\mathbf{A}}_h\mathbf{v}=\lambda {\mathbf{B}}_h\mathbf{v}.$$

$${\lambda }_1\le {\lambda }_2\le ...\le {\lambda }_{N_h}.$$

$${\mathbf{u}}_h(\mathbf{x},t)=\sum _{n=1}^{N_h}{b}_n\exp (-{\lambda }_n(t-t_{s-1})){\mathbf{v}}_n(\mathbf{x}),$$

$$a_n(t)=b_n\exp (-{\lambda }_n(t-t_{s-1})),{\mathbf{w}}_n(\mathbf{x})={\mathbf{v}}_n(\mathbf{x}),n=1,2,...,N.$$

$${\lambda }_{n+1}=\mathrm{Re}{\lambda }_n-i\mathrm{Im}{\lambda }_n.$$

$$\begin{split}a_{n}(t)\bm{w}_{n}(\bm{x})&=b_{n}\mathrm{Re}\big{(}\exp(-\lambda_{n}(t-t_{s-1}))\bm{v}_{n}(\bm{x})\big{)},\\ a_{n+1}(t)\bm{w}_{n+1}(\bm{x})&=b_{n+1}\mathrm{Im}\big{(}\exp(-\lambda_{n}(t-t_{s-1}))\bm{v}_{n}(\bm{x})\big{)}.\end{split}$$

$${\mathbf{u}}_h^s(\mathbf{x})=\sum _{n=1}^{N_h}{b}_n{\mathbf{v}}_n(\mathbf{x}).$$

$${\mathbf{A}}_h^T\tilde{\mathbf{v}}=\lambda {\mathbf{B}}_h^T\tilde{\mathbf{v}}.$$

Full text

††journal: arXiv.org

1 Introduction

Among the basic physical processes occurring in a nuclear reactor (Duderstadt and Hamilton, 1976), the focus is on neutron transport. To describe that process, a rather complex integro-differential equation is involved, in which the neutron flux distribution depends on time, energy, spatial and angular variables (Hetrick, 1971; Stacey, 2007). In practical calculations of nuclear reactors, as a rule, simpler systems of equations in the multigroup diffusion approximation are used (Marchuk and Lebedev, 1986; Lewis and Miller, 1993; Sutton and Aviles, 1996; Cho, 2005). Currently diffusion models are derived and applied using sophisticated homogenization methodologies (see Sanchez (2009); Dugan et al. (2016)) which define parameters of the multigroup diffusion equations that enable one to take into account transport effects.

The computational model is a system of coupled second-order parabolic equations, which is supplemented by a system of ordinary differential equations to take into account delayed neutrons. Engineering neutronics codes designed to simulate neutron transport in the few-group diffusion approximation are based, mainly, on finite-difference spatial approximations. To improve the calculation accuracy, nodal methods are widely used (see, for example, Smith (1979); Lawrence (1986)). Applying these methods, a rather coarse grid can be used: several points per assembly in plane and several dozen layers in height. Nodal methods are based on the representation of the neutron flux within the computational element as a small degree polynomial or a set of functions along one of the coordinates (or in plane). In some cases, nodal methods can be related (Grossman and Hennart, 2007) with special variants of finite-element approximation. It should be noted that to increase the finite-element accuracy of approximate solution of boundary value problems, it is more appropriate to use standard procedures with condensed computational grids finite elements of higher degree. Such technology is used in Vidal-Ferrandiz et al. (2014); Avvakumov et al. (2017) when considering spectral problems for the multigroup diffusion equations. Modern computational algorithms for solving spectral problems in neutron theory are discussed by Gill and Azmy (2011); Willert et al. (2014).

When modeling the dynamics of neutron-physical processes, standard methods for the approximate solution of non-stationary problems are used (Sutton and Aviles, 1996; Cho, 2005; Stacey, 2007). The most attention is paid to two-layer schemes with weights ($\theta$-method) (Ascher, 2008; LeVeque, 2007; Hundsdorfer and Verwer, 2003), also Runge-Kutta and Rosenbroek schemes (Butcher, 2008; Hairer and Wanner, 2010) are used. To characterize the reactor dynamics, the spectral parameter $\alpha$ is used. It is defined as the main eigenvalue of the time-eigenvalue ($\alpha$-eigenvalue) problem, which is associated with the nonstationary neutron diffusion equations (Bell and Glasstone, 1970; Modak and Gupta, 2007; Verdu et al., 2010). By analogy with the usual problems of heat conduction (see, for example, Luikov (1968); Samarskii and Vabishchevich (1996)) we can discriminate the regular regime of the reactor. For large times, the neutron flux behavior has an asymptotic character when one can speak of the space-time factorization of a solution whose amplitude is $\exp(-\alpha t)$, the form-function is the spectral problem eigenfunction. To model the reactor regular regime, it is necessary to focus on the use of completely implicit schemes, while the Crank-Nicholson scheme is unsuitable (Avvakumov et al., 2016).

Real three-dimensional dynamic neutronic calculations require the use of large computational grids for large characteristic times which in turn determine the use of modern multiprocessor computer systems. Parallel computational algorithms are based on a sequence of solving more simple problems for individual processes. Advance is achieved through the use of decoupling technology: splitting by physical processes (Vabishchevich, 2014)), decomposition of the computational domain into subdomains (Toselli and Widlund, 2005), iterative methods for solving systems of algebraic equations (Saad, 2003). In the case of spectral problems for the neutron diffusion, domain decomposition methods are used, for example, in Guérin et al. (2010). Features of solving non-stationary problems on parallel computers are taken into account by constructing special iterative methods such as the parareal in time algorithm (Maday and Turinici, 2005). This approach was implemented Baudron et al. (2014) for numerical calculation of transient multigroup neutron kinetics equations involving a time delayed contributions.

In the theory and practice of neutronics calculations, fast methods to obtain approximate solutions are given great attention. In this connection, a class of methods for modeling the nonstationary group neutron diffusion should be noted. It is associated with the multiplicative representation of the space-time factorization methods and the quasistatic method (Dodds Jr, 1976; Chou et al., 1990; Goluoglu and Dodds, 2001; Dulla et al., 2008; Dahmani et al., 2001). In this case, an approximate solution is represented as a product of two functions: one of which depends on the time and is related to the amplitude, the second (shape function) describes the spatial distribution. The shape function is often associated with the fundamental eigenfunction of certain eigenvalue problems for neutron diffusion equations.

When using the quasistatic method, the problem is significantly simplified, thus, it is doubtful to obtain good accuracy for an approximate solution, in particular, in dynamic regimes with a complex redistribution of the neutron flux. For this reason, the more general approach of the modal method has been successfully developed (Stacey, 1967, 1969; Sutton and Aviles, 1996). In this case, the solution is represented in the form of a sum of several dominant eigenvalues with time-dependent coefficients.

To characterize the reactor dynamic processes some spectral problems (Bell and Glasstone, 1970; Hetrick, 1971; Stewart, 1976; Stacey, 2007) were considered. The processes occurring in a nuclear reactor are essentially non-stationary. The stationary state of neutron flux is characterised by local balancing of neutron absorption and generation and is usually described by solution of a spectral problem (Lambda modes problem, $\lambda$-eigenvalue problem). The fundamental eigenvalue (maximal eigenvalue) is called $k$-effective of the reactor core. The nodal method based on the use of the $\lambda$-eigenvalue problem is discussed, for example, Verdú et al. (1998); Miró et al. (2002); González-Pintor et al. (2009). In particular, it should be paid attention to issues related with the calculation of the related system of equations for the nonstationary expansion coefficients.

Nonstationary processes can naturally be described on the basis of the approximate solution expansion in time-eigenvalue of $\alpha$-eigenvalue problem (Ginestar et al., 2002; Verdu et al., 2010; Verdú and Ginestar, 2014). In a simpler model without delayed neutrons, modal methods were used in Modak and Gupta (2007). The principal point is connected with the fact that when using this approach, we deal with an very big system of equations for the coefficients. It should also be noted that the eigenvalues are complex for both $\lambda$- and $\alpha$-eigenvalue problem. To set the initial state, this leads to the need to solve the appropriate adjoint spectral problems.

In this paper, we formulate a general strategy for the approximate solution of nonstationary problems of neutron transport in nuclear reactors, which is oriented to fast real-time calculations using the State Change Modal (SCM) method. The dynamic behavior of a nuclear reactor is considered as a sequence of its states characterized by its set of constant coefficients of the multigroup diffusion equations. It is considered that the transition from one state to another occurs instantaneously. For a separate state, the neutron flux is calculated using the modal method to represent the problem’s solution in the form of expansion in the dominant eigenfunctions of the $\alpha$-eigenvalue problem taking into account delayed neutrons (the full Alpha method). Real-time calculations are provided by the fact that the desired set of eigenvalues and eigenfunctions of the reactor deterministic state are calculated in advance.

The paper is organized as follows. The dynamic model of a nuclear reactor based on the multigroup diffusion equations is given in Section 2. The general strategy of numerical modeling of nonstationary processes based on the SCM method is described in Section 3. The key computational aspects of the SCM technology are discussed in Section 4. Two-dimensional test problem for VVER-1000 reactor is discussed in Section 5. Modeling of the dynamic process, which corresponds to two reactor states (supercritical and subcritical): a transition from the regular regime of the critical state to the subcritical state of the reactor. The results of the work are summarized in Section 6.

2 Problem statement

The neutron flux is modelled in multigroup diffusion approximation. The neutron dynamics is considered in the bounded convex two-dimensional or three-dimensional area $\Omega$ ($\bm{x}=\{x_{1},...,x_{d}\}\in\Omega,\ d=2,3$) with boundary $\partial\Omega$. The neutron transport is described by the system of equations:

[TABLE]

Here $\phi_{g}(\bm{x},t)$ — neutron flux of $g$ group at point $\bm{x}$ and time $t$, $G$ — number of energy groups, $v_{g}$ — effective velocity of neutrons in the group $g$, $D_{g}(\bm{x})$ — diffusion coefficient, $\Sigma_{rg}(\bm{x},t)$ — removal cross-section, $\Sigma_{s,g^{\prime}\rightarrow g}(\bm{x},t)$ — scattering cross-section from group $g^{\prime}$ to group $g$, $\beta$ — effective fraction of delayed neutrons, $\chi_{g}$, $\widetilde{\chi}_{g}$ — spectra of prompt and delayed neutrons, $\nu\Sigma_{fg}(\bm{x},t)$ — generation cross-section of group $g$, $c_{m}$ — density of sources of delayed neutrons of $m$-type, $\lambda_{m}$ — decay constant of sources of delayed neutrons, $M$ — number of types of delayed neutrons. The density of sources of delayed neutrons is described by the equations:

[TABLE]

where $\beta_{m}$ is a fraction of delayed neutrons of $m$-type, and

[TABLE]

System of equations (1), (2) is supplemented with corresponding initial and boundary conditions.

The albedo-type conditions are set at the boundary $\partial\Omega$ of the area $\Omega$:

[TABLE]

where $n$ — outer normal to the boundary $\partial\Omega$. Initial state is defined in the following manner:

[TABLE]

Let’s write the boundary problem (1)–(4) in operator form. The vectors $\bm{\phi}=\{\phi_{1},\phi_{2},...,\phi_{G}\}$, $\bm{c}=\{c_{1},c_{2},...,c_{M}\}$ and matrices are defined as follows:

[TABLE]

where

[TABLE]

is the Kronecker symbol. We shall use the set of vectors $\bm{\phi}$, whose components satisfy the boundary conditions (3). Using the set definitions, the system of equations (1), (2) can be written in the form of first-order equation of evolution:

[TABLE]

The Cauchy problem is formulated for equations (5) when

[TABLE]

where taken into account (4) $\bm{\phi}^{0}=\{\phi_{1}^{0},\phi_{2}^{0},...,\phi_{G}^{0}\}$, $\bm{c}^{0}=\{c_{1}^{0},c_{2}^{0},...,c_{M}^{0}\}$.

For an approximate solution of the Cauchy problem (5), (6) two basic approaches are used. The first of these is connected with the use of standard two-layer or three-layer schemes, which are widely used in the numerical solution of parabolic problems (Samarskii, 2001).

The operator matrices $V,D$ are diagonal, and $S$ is the lower triangular matrix. The essential binding of the equations is due only to the neutron generation operator $R$. The problem of choosing design schemes amongst stable difference schemes, which is optimal for some additional criteria, is topical. In the theory of difference schemes, a class of asymptotically stable difference schemes is distinguished, which (Samarskii and Vabishchevich, 1996) provide the correct behavior of the approximate solution for large times. In the theory of numerical methods for solving systems of ordinary differential equations (Butcher, 2008; Gear, 1971) the concept of $L$-stable methods is introduced, in which, from several other positions, the asymptotic behavior of the approximate solution is also monitored for large times. A completely implicit scheme has better asymptotic properties than a symmetric scheme (the Crank-Nicholson scheme) (see Vabishchevich (2012)), which is important in the study of the regular regime of a nuclear reactor (Avvakumov et al., 2016).

The second class is the numerical-analytical methods for solving problem (5), (6), the most striking example of which are the modal methods noted above (Stacey, 1967, 1969; Sutton and Aviles, 1996). They take into account the linear nature of the problem under consideration, the independence of the coefficients of the system of equations as a function of time. This allows us to construct an approximate solution by the method of separation of variables in the numerical determination of the dependence of the solution on the spatial variables.

3 State change modal method

The nuclear reactor is always non-stationary. The limiting case of reaching a steady state (critical reactor) is observed only for certain coefficients of the equations system (5). We will use the following simplified description of the dynamic processes in a nuclear reactor.

In selected time interval, the non-stationary neutron flux is determined by the nuclear reactor state. The state of the reactor is characterized by the constant coefficients of the system of multigroup diffusion equations (1), (2).

Dynamic processes in a nuclear reactor can be considered as a change of states (see Fig.1). At a certain time $t=t_{s},\ s=1,2,...$ an instantaneous change of state occurs. The state $s$ is defined by the parameters in equations (5):

[TABLE]

[TABLE]

Simulation of the dynamic behavior of the reactor consists in solving the sequence of subtasks for the individual states of the reactor. The initial condition for the state $s$ (at $t=t_{s-1}$) is the final state of the reactor for the state $s-1$.

An approximate description of the non-stationary process at a separate stage is based on modal approximation. An approximate solution is sought in the form of decomposition in eigenfunctions of time and $\alpha$-eigenvalue problem. Finite-element approximation in space is used.

In a separate stage $s$ the following system of equations is considered

[TABLE]

supplemented by the corresponding initial conditions

[TABLE]

Let’s $\bm{u}=\{\bm{\phi},\bm{c}\}$. Rewrite the system of equations (7) as

[TABLE]

with constants

[TABLE]

where $I$ is the identity matrix. From (8) one can obtain

[TABLE]

After approximating over the space by the finite volume method or by the finite element method from (9), (10) we turn to the Cauchy problem for a linear system of ordinary differential equations with constant coefficients:

[TABLE]

[TABLE]

where $h$ is the discretization parameter. The main feature of the problems we are considering is that the matrices $\bm{A}_{h}$ and $\bm{B}_{h}$ are real and asymmetric.

The modal approximation corresponds to the representation of the approximate solution ($\bm{u}_{h}\approx\bm{u}_{N}$) of problem (11), (12) in the following form

[TABLE]

where $N$ is the number of dominant eigenvalues of the spectral problem, $\bm{w}_{n}(\bm{x})$ — corresponding eigenfunctions.

Let us define eigenfunctions and eigenvalues as the solution of the $\alpha$-eigenvalue problem:

[TABLE]

In the simplest case all eigenvalues of the spectral problem (14) are real:

[TABLE]

Under these conditions (Laub, 2005; Ortega, 1987) the general solution of equation (11) is

[TABLE]

that is in (13)

[TABLE]

In the general case, eigenfunctions and eigenvalues of the spectral problem (14) are complex. Taking into account the reality of the matrix coefficients $\bm{A}_{h},\ \bm{B}_{h}$ complex eigenvalues appear as pairs of complex conjugate numbers. For example, we have a pair of $n,n+1$:

[TABLE]

Then in the representation (13) we obtain

[TABLE]

A special attention should be paid to define the coefficients $a_{n}(t_{s-1})=b_{n},\ n=1,2,...,N$. For this, the initial condition (12) is involved. For example, in the case of real eigenvalues, we have

[TABLE]

This representation is not very suitable for practical use with modal approximation, when we work only with dominant eigenfunctions.

The initial condition includes two components $\bm{u}_{h}^{s}(\bm{x})=(\bm{\phi}_{h}^{s}(\bm{x}),\bm{c}_{h}^{s}(\bm{x}))$. Dynamic behaviour of these components is due to different time-scale processes. Delayed neutrons source determines slow processes, when $\bm{c}(\bm{x},t)$ changes slightly with the reactor state change. In contrast, neutron flux $\bm{\phi}(\bm{x},t)$ determines fast processes when the reactor state changes. By virtue of this separation of dynamic processes, we model the slow phase of the dynamics of the reactor with modal approximation and orientate ourselves on the approximate prediction of the initial state for delayed neutrons, only the function $\bm{c}_{h}^{s}(\bm{x})$ is approximated. The approximation $\bm{\phi}_{h}^{s}(\bm{x})$ is not of interest to us; we do not model a fast phase of the state change.

The standard approach for the decomposition of the function $\bm{u}_{h}^{s}(\bm{x})$ over the system of non-orthogonal functions $\bm{v}_{n}(\bm{x}),\ n=1,2,...,N_{h}$ consists in using the biorthogonal system of functions (Henry, 1975; Brezinski, 1991). Consider the spectral problem adjoint to (14)

[TABLE]

The eigenfunctions of problems (14) and (16) are orthogonal (Laub, 2005; Ortega, 1987) in the sense of the equality

[TABLE]

where $(\cdot,\cdot)$ means corresponding scalar product. In view of this, one can obtain

[TABLE]

With known solutions of the spectral problems (14), (16) the solution is represented in the form (15), (17).

In the approximate solution of problem (11), (12) only the first $N$ coefficients $b_{n}$ in (17) are used (see (13)):

[TABLE]

where $\bm{v}_{n}(\bm{x})=(\bm{\phi}_{n}(\bm{x}),\bm{c}_{n}(\bm{x}))$. In this case, the spectral problems (14), (16) are solved for $N$ dominant eigenvalues.

The solution of the adjoint spectral problem is involved only for calculating the initial condition coefficients. This complication of the problem is not always justified. Therefore, it is worth to use simpler algorithms for obtaining the coefficients $b_{n},\ n=1,2,...,N$ in (18). We can define them, for example, based on linear least squares (Björck, 1996; Verdú and Ginestar, 2014). In this case, one can obtain

[TABLE]

To find the coefficients, a system of linear equations is solved.

The state change modal method is based on the following calculating scheme.

Off-line calculation.

Calculation of the coefficients of the mathematical model of the multigroup diffusion approximation for the isolated reactor states, which is performed in advance. The status passport also includes calculated dominant eigenvalues and eigenfunctions of the $\alpha$-eigenvalue problem (14). These data can be supplemented by dominant eigenvalues and eigenvalues of the conjugate eigenvalue problem (16).

On-line calculation.

Real-time modeling is carried out on the basis of the modal solution of the problem (11), (12). The coefficients in the representation (18) are calculated from the initial condition using (17) or (19). The solution for other time intervals is determined according to (15).

4 Computational aspects of the state change modal method

In practical implementation of this approach for approximate solution of reactor dynamic problems, the key aspects are spatial approximation and numerical solution of spectral problems. We will not discuss the available opportunities in this direction, but will give only a brief description of how these issues are solved for the below calculations.

To approximate in space variables, we use the finite element method (Brenner and Scott, 2008; Quarteroni and Valli, 2008). Let $H^{1}(\Omega)$ is the Sobolev space consisting of scalar functions $v$ such that $v^{2}$ and $|\nabla v|^{2}$ have a finite integral in $\Omega$. For vector functions $\bm{v}=\{v_{1},v_{2},...,v_{d}\}$ we define similarly $V^{d}=[H^{1}(\Omega)]^{d}$. For test functions we use the notation $\bm{\chi}=\{\chi_{1},\chi_{2},...,\chi_{G}\}$, $\bm{s}=\{s_{1},s_{2},...,s_{M}\}$. In the variational form, problem (5) has the following form: find $\bm{\phi}\in V^{D},\ \bm{c}\in V^{M}$, for which

[TABLE]

for all $\bm{\chi}\in V^{D},\ \bm{s}\in V^{M}$.

Further, we must pass from the continuous variational problem (20) to the discrete problem. We introduce finite-dimensional finite element spaces $V_{h}^{D}\subset V^{D}$, $V_{h}^{M}\subset V^{M}$. The discrete variational problem is formulated as follows: find $\bm{\phi}^{h}\in V_{h}^{D},\ \bm{c}^{h}\in V_{h}^{M}$, such that

[TABLE]

for all $\bm{\chi}^{h}\in V_{h}^{D},\ \bm{s}^{h}\in V_{h}^{M}$. For two-dimensional problems, scalar functions (components of vector functions) are approximated on a triangular grid using Lagrangian finite elements with polynomials of degree 1, 2 and 3. The calculation of the first eigenvalues and the corresponding eigenfunctions is a standard problem in computational mathematics (Saad, 2011). It is necessary to note some fundamental features of the spectral problems (14) и (16). The first feature is related to the fact that the problems under consideration are large-scale problems: two-dimensional or three-dimensional in space and many unknown quantities (a system of equations). This means that in applied problems, we have to focus on the use of computers of parallel architecture. The second is due to the asymmetry of the matrices. This leads to the appearance of complex roots.

In our study (see Avvakumov et al. (2017, 2016)) we focus on the use of well-designed algorithms and relevant free software. To solve spectral problems with non-symmetrical matrices we use the SLEPc (Scalable Library for Eigenvalue Problem Computations, http://slepc.upv.es/). This library has traditionally been widely used (see, for instance, Hernandez et al. (2003, 2005)) for numerical solution of the spectral problems in nuclear reactor calculations. We use a Krylov-Schur algorithm, a variation of Arnoldi method, proposed by (Stewart, 2001).

5 The test: the dynamics of the VVER-1000 reactor during the transition from the supercritical mode to the subcritical mode

A test problem for a VVER-1000 reactor without a reflector (Chao and Shatilla, 1995) is considered in the two-dimensional approximation ($\Omega$ is a reactor core cross-section).

5.1 General description

The geometric model of the VVER-1000 core consists of a set of hexagonal-shaped cassettes and is shown in Fig.2, where fuel assemblies of various types are shown. The assembly wrench size is 23.6 cm.

For an approximate solution of the problem, regular triangular grids are used. The number of triangles per cassette $\kappa$ varies from 6 to 96 (Fig.3). A two-group approximation is used taken into account delayed neutrons:

[TABLE]

The supercritical state of the reactor is characterized by a set of coefficients, which are given in Table 1. The following boundary conditions (3) are used: $\gamma_{g}=0.5,\ g=1,2$. The following delayed neutrons parameters are used: one group of delayed neutrons with effective fraction $\beta_{1}=6.5\cdot 10^{-3}$ and decay constant $\lambda_{1}=0.08$ s*-1*. Neutron velocity $v_{1}=1.25\cdot 10^{7}$ cm/s and $v_{2}=2.5\cdot 10^{5}$ cm/s.

5.2 Supercritical state: $\alpha$-eigenvalue problem

Below the results of a numerical solution of the $\alpha$-eigenvalue problem (14) are presented. In the framework of the two-group approximation and taking into account the delayed neutrons, one can write

[TABLE]

The dominant eigenvalues $\alpha_{n}=\lambda_{n}^{(\alpha)},\ n=1,2,...,N$ are searched for at

[TABLE]

Similar calculations of the eigenvalues for the VVER-1000 test problem without delayed neutrons can be found in Avvakumov et al. (2017).

The results of solving the spectral problem (23) for the first eigenvalues $\alpha_{n}=\lambda_{n}^{(\alpha)},\ n=1,2,...,N$, $N=10$ on different computational grids using different finite element approximations are shown in Table 2, 3. The eigenvalues $\alpha_{2},\alpha_{3}$, $\alpha_{4},\alpha_{5}$, $\alpha_{9},\alpha_{10}$ of the spectral problem (23) are complex with small imaginary parts, the eigenvalues $\alpha_{1},\alpha_{6}$, $\alpha_{7},\alpha_{8}$ are real.

In our example, the main eigenvalue is negative and therefore the major harmonic will increase, and all others will fade. This demonstrates the regular mode of the reactor operation. The value $\alpha=\lambda_{1}^{(\alpha)}$ itself determines the neutron flux amplitude and is directly related to the reactor period in the regular regime.

Dominant eigenfunctions of the spectral problem (23) for the first group are shown in Figs.4–8. The calculations are performed on a grid with $\kappa=96$ using finite elements of degree $p=3$. The main eigenfunctions for the second group of $\varphi^{(1)}_{2}$ and delayed neutrons source $s^{(1)}$ are shown in Fig.9.

5.3 Adjoint spectral problem

Analogous data were obtained for approximate solution of the adjoint spectral problem (16). The eigenvalues of the spectral problems (14) and (16) coincide. Their difference from each other is an indirect measure of the accuracy of the numerical solution. Data on the dominant eigenvalues, which are given in Table 4 ($k=96,\ p=3$), show that the eigenvalues of the main and adjoint spectral problems are close to each other with good accuracy.

The spectral problems under consideration are characterized by small imaginary parts of the eigenvalues. Therefore, we can expect that the eigenfunctions of problem (14) are close to orthogonal. As illustration, Table 5 contains data for the scalar products $(\phi_{1}^{(n)},\phi_{1}^{(m)})$ for the first 10 eigenfunctions. For convenience of comparison, the eigenfunctions are normalized in $L_{2}(\Omega)$:

[TABLE]

The maximum non-orthogonality (for $(\phi_{1}^{(1)},\phi_{1}^{(7)})$) does not exceed 10 %. The biorthogonality condition of the eigenfunctions of the fundamental functions (see (14)) and the adjoint spectral problems (see (16)) is valid with approximately the same accuracy. This observed error can be related to an approximate calculation of eigenvalues and eigenfunctions.

Within the modal method, we can not rely on high accuracy when considering a relatively small number of dominant eigenvalues. Therefore, in the example under consideration, we can assume that the eigenvalues are real, and the corresponding eigenfunctions are orthogonal. Instead of (17) coefficients are used

[TABLE]

to approximate the initial condition.

5.4 Subcritical state

In the supercritical mode, due to the sufficiently large magnitude of the main eigenvalue, the regular regime of the reactor is rapidly developing, where

[TABLE]

Here $\bm{v}_{1}^{0}(\bm{x})$ is the first mode of the supercritical state. We consider the problem with the transition from this supercritical state at $t_{0}=0$ to the subcritical state.

The subcritical stage is characterized by a 15% increase in the coefficient $\Sigma_{2}$ for material 4 in the VVER-1000 test diffusion constants (see Table 1). Thus, the dynamics of the reactor is as follows:

[TABLE]

The initial state is characterized by specifying the initial conditions at $t_{0}=0$ as

[TABLE]

The calculational results of the dominant eigenvalues for the subcritical state are presented in Tables 6, 7. In this case even the first eigenvalues not significantly differ from each other.

For an approximate solution, we use the following formulation:

[TABLE]

where the coefficients $b_{n},\ n=1,2,...,N$ are calculated according to the given initial condition (24). These coefficients for $N=50$ are shown in Fig.10. As we see, an approximate solution can be described by first mode only.

We distinguish two phases of the dynamic process: fast and slow. In the fast phase, the initial condition (25) is rearranged to the initial condition, which corresponds to (26): from the function $\bm{u}(\bm{x},0)$ to the function $\bm{u}_{N}(\bm{x},0)$. The slow phase is associated with the evolution of the solution according to (26). Within the state change modal technology, the fast phase is not modeled at all.

The beginning and the end of the fast phase are illustrated through the calculational data shown in Fig. 11. The results were obtained with $N=10$. We can note very small changes in the topology of the initial and reconstructed initial conditions. Let us pay attention to the substantial restructuring of the solution, which is illustrated by large changes in the neutron flux amplitudes for the first and second groups.

5.5 The asymmetric perturbation

Consider a more complex transition to a subcritical state. The subcritical stage will be characterized by a different increase in the coefficient $\Sigma_{2}$ for material 4 in the diffusion constants in the upper and lower half of the reactor cross-section (see Fig.2). Now let the reactor dynamics corresponds to the following transformation

[TABLE]

The calculational results for the dominant eigenvalues of the reactor subcritical asymmetric state are presented in Tables 8, 9. All eigenvalues for this reactor state are real.

The coefficients $b_{n},\ n=1,2,...,N$, $N=50$ of the approximate solution (26) with the initial condition (25) are shown in Fig.12. In the case under consideration, the approximate solution contains several modes and can not be described only by the first mode. The fast transition phase is illustrated in Fig. 13. Calculations of the end of the fast phase (the functions $\bm{u}_{N}(\bm{x},0)$) are performed at $N=10$.

5.6 Comparison with the nonstationary problem solution

An approximate solution, obtained using modal approximation, can be compared with the dynamic problem solution. The boundary value problem for the system of equations (22) is solved. Fully implicit scheme on a uniform grid in time with a sufficiently small step $\tau=0.0025$ is used (see details in Avvakumov et al. (2016)). The dynamics of the neutron power of the nuclear reactor $P$ and the delayed neutrons source $C$ at the initial stage during the transition from the critical state to the subcritical in the case of a symmetric perturbation is shown in Fig.14. Here

[TABLE]

There is a rapid change in neutron power over a short period of time, while the delayed neutrons source changes rather slow. The dynamics of the slow phase is illustrated in Figs.15,16. Here the solution of full equations (dynamic in Figs.15,16), and the modal approximation solution (modal) are presented. Similar data for the asymmetric perturbation of the reactor state are shown in Figs.17,18. We see that the integral characteristics of the reactor dynamics at slow stage are calculated with good accuracy.

6 Conclusions

The problem of simulation of reactor dynamic processes is considered on the basis of multigroup neutron diffusion equations accounting for delayed neutrons. The modal approximation is used: an approximate solution is represented as an expansion on limited number of dominant eigenfunctions of the $\alpha$-eigenvalue spectral problem.

Numerical simulation of reactor non-stationary processes is carried out on the basis of a successive change in the states of the reactor. These states are characterized by a set of constant parameters to describe the multigroup neutron flux behavior. The state change modal method was developed. The phase, which described fast transition to the approximate solution, is selected as a set of dominant modes. At a slow phase of the reactor dynamics, the solution is based on the evolution of dominant modes.

The computational implementation of the state change modal method is based on the previously calculated (of-line calculation) eigenfunctions and eigenvalues of the $\alpha$-eigenvalue spectral problem. Fast determination of dominant modes and calculation of the reactor neutron flux at selected times are based on on-line calculation. The classical Lagrange finite elements $p=1,2,3$ are used for the spatial approximation. Accuracy control is performed using condensed grids. Spectral problems are solved numerically using well-developed free software SLEPc.

Test calculations are made in two-dimensional analysis using a two-group diffusion approximation. Calculations of dominant modes for a reactor supercritical state are performed. The major mode solution, which determines the reactor regular regime, is used as the initial condition for transition to the subcritical state. The modeling of the reactor state change as a transfer from one state to another state is carried out for two cases. The first of them (symmetric perturbation) is due to the uniform change in the absorbing material properties. The second case (asymmetric perturbation) deals with a non-uniform change in the absorbing material properties (in two halves over the reactor cross-section).

Comparison of the calculational results obtained by using two methods (one based on modal approximation and another based on the full dynamics calculation) shows the acceptable accuracy in calculation of neutron power and delayed neutrons source for the VVER-1000 test problem.

Acknowledgements

This work are supported by the Russian Foundation for Basic Research (# 16-08-01215) and by the grant of the Russian Federation Government

(# 14.Y26.31.0013).