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 integrodifferential 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 secondorder 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 fewgroup diffusion approximation are based, mainly, on finitedifference 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 finiteelement approximation. It should be noted that to increase the finiteelement 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 VidalFerrandiz 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 neutronphysical processes, standard methods for the approximate solution of nonstationary problems are used (Sutton and Aviles, 1996; Cho, 2005; Stacey, 2007). The most attention is paid to twolayer schemes with weights (method) (Ascher, 2008; LeVeque, 2007; Hundsdorfer and Verwer, 2003), also RungeKutta and Rosenbroek schemes (Butcher, 2008; Hairer and Wanner, 2010) are used. To characterize the reactor dynamics, the spectral parameter is used. It is defined as the main eigenvalue of the timeeigenvalue (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 spacetime factorization of a solution whose amplitude is , the formfunction 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 CrankNicholson scheme is unsuitable (Avvakumov et al., 2016).
Real threedimensional 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 nonstationary 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 spacetime 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 timedependent 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 nonstationary. 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, eigenvalue problem). The fundamental eigenvalue (maximal eigenvalue) is called effective of the reactor core. The nodal method based on the use of the eigenvalue problem is discussed, for example, Verdú et al. (1998); Miró et al. (2002); GonzálezPintor 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 timeeigenvalue of 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  and 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 realtime 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 eigenvalue problem taking into account delayed neutrons (the full Alpha method). Realtime 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. Twodimensional test problem for VVER1000 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 twodimensional or threedimensional area () with boundary . The neutron transport is described by the system of equations:
(1) 
Here — neutron flux of group at point and time , — number of energy groups, — effective velocity of neutrons in the group , — diffusion coefficient, — removal crosssection, — scattering crosssection from group to group , — effective fraction of delayed neutrons, , — spectra of prompt and delayed neutrons, — generation crosssection of group , — density of sources of delayed neutrons of type, — decay constant of sources of delayed neutrons, — number of types of delayed neutrons. The density of sources of delayed neutrons is described by the equations:
(2) 
where is a fraction of delayed neutrons of type, and
System of equations (1), (2) is supplemented with corresponding initial and boundary conditions.
The albedotype conditions are set at the boundary of the area :
(3) 
where — outer normal to the boundary . Initial state is defined in the following manner:
(4) 
Let’s write the boundary problem (1)–(4
) in operator form. The vectors
, and matrices are defined as follows:where
is the Kronecker symbol. We shall use the set of vectors , whose components satisfy the boundary conditions (3). Using the set definitions, the system of equations (1), (2) can be written in the form of firstorder equation of evolution:
(5) 
The Cauchy problem is formulated for equations (5) when
(6) 
where taken into account (4) , .
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 twolayer or threelayer schemes, which are widely used in the numerical solution of parabolic problems (Samarskii, 2001).
The operator matrices are diagonal, and is the lower triangular matrix. The essential binding of the equations is due only to the neutron generation operator . 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 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 CrankNicholson 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 numericalanalytical 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 nonstationary. 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 nonstationary 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 an instantaneous change of state occurs. The state is defined by the parameters in equations (5):
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 (at ) is the final state of the reactor for the state .
An approximate description of the nonstationary 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 eigenvalue problem. Finiteelement approximation in space is used.
In a separate stage the following system of equations is considered
(7) 
supplemented by the corresponding initial conditions
(8) 
Let’s . Rewrite the system of equations (7) as
(9) 
with constants
where
is the identity matrix. From (
8) one can obtain(10) 
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:
(11) 
(12) 
where is the discretization parameter. The main feature of the problems we are considering is that the matrices and are real and asymmetric.
The modal approximation corresponds to the representation of the approximate solution () of problem (11), (12) in the following form
(13) 
where is the number of dominant eigenvalues of the spectral problem, — corresponding eigenfunctions.
Let us define eigenfunctions and eigenvalues as the solution of the eigenvalue problem:
(14) 
In the simplest case all eigenvalues of the spectral problem (14) are real:
Under these conditions (Laub, 2005; Ortega, 1987) the general solution of equation (11) is
(15) 
that is in (13)
In the general case, eigenfunctions and eigenvalues of the spectral problem (14) are complex. Taking into account the reality of the matrix coefficients complex eigenvalues appear as pairs of complex conjugate numbers. For example, we have a pair of :
Then in the representation (13) we obtain
A special attention should be paid to define the coefficients . For this, the initial condition (12) is involved. For example, in the case of real eigenvalues, we have
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 . Dynamic behaviour of these components is due to different timescale processes. Delayed neutrons source determines slow processes, when changes slightly with the reactor state change. In contrast, neutron flux 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 is approximated. The approximation 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 over the system of nonorthogonal functions consists in using the biorthogonal system of functions (Henry, 1975; Brezinski, 1991). Consider the spectral problem adjoint to (14)
(16) 
The eigenfunctions of problems (14) and (16) are orthogonal (Laub, 2005; Ortega, 1987) in the sense of the equality
where means corresponding scalar product. In view of this, one can obtain
(17) 
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 coefficients in (17) are used (see (13)):
(18) 
where . In this case, the spectral problems (14), (16) are solved for 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 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
(19) 
To find the coefficients, a system of linear equations is solved.
The state change modal method is based on the following calculating scheme.
 Offline 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 eigenvalue problem (14). These data can be supplemented by dominant eigenvalues and eigenvalues of the conjugate eigenvalue problem (16).
 Online calculation.
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 is the Sobolev space consisting of scalar functions such that and have a finite integral in . For vector functions we define similarly . For test functions we use the notation , . In the variational form, problem (5) has the following form: find , for which
(20) 
for all .
Further, we must pass from the continuous variational problem (20) to the discrete problem. We introduce finitedimensional finite element spaces , . The discrete variational problem is formulated as follows: find , such that
(21) 
for all . For twodimensional 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 largescale problems: twodimensional or threedimensional 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 welldesigned algorithms and relevant free software. To solve spectral problems with nonsymmetrical 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 KrylovSchur algorithm, a variation of Arnoldi method, proposed by (Stewart, 2001).
5 The test: the dynamics of the VVER1000 reactor during the transition from the supercritical mode to the subcritical mode
A test problem for a VVER1000 reactor without a reflector (Chao and Shatilla, 1995) is considered in the twodimensional approximation ( is a reactor core crosssection).
5.1 General description
The geometric model of the VVER1000 core consists of a set of hexagonalshaped 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 varies from 6 to 96 (Fig.3). A twogroup approximation is used taken into account delayed neutrons:
(22) 
Material  1  2  3  4  5 

1.38320e0  1.38299e0  1.39522e0  1.39446e0  1.39506e0  
3.86277e1  3.89403e1  3.86225e1  3.87723e1  3.84492e1  
2.48836e2  2.62865e2  2.45662e2  2.60117e2  2.46141e2  
6.73049e2  8.10328e2  8.44801e1  9.89671e2  8.93878e2  
1.64977e2  1.47315e2  1.56219e2  1.40185e2  1.54981e2  
4.81619e3  4.66953e3  6.04889e3  5.91507e3  6.40256e3  
8.46154e2  8.52264e2  1.19428e1  1.20497e1  1.29281e1 
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: . The following delayed neutrons parameters are used: one group of delayed neutrons with effective fraction and decay constant s. Neutron velocity cm/s and cm/s.
5.2 Supercritical state: eigenvalue problem
Below the results of a numerical solution of the eigenvalue problem (14) are presented. In the framework of the twogroup approximation and taking into account the delayed neutrons, one can write
(23) 
The dominant eigenvalues are searched for at
Similar calculations of the eigenvalues for the VVER1000 test problem without delayed neutrons can be found in Avvakumov et al. (2017).
6  0.22557  0.04241 3.08808e06  0.06588 4.80449e07  

1  24  0.82690  0.03777 5.37884e06  0.06489 1.37315e06 
96  1.74998  0.03619 5.69002e06  0.06456 1.40299e06  
6  2.10154  0.03592 4.96474e06  0.06452 1.21320e06  
2  24  2.46601  0.03562 5.78277e06  0.06445 1.40897e06 
96  2.50375  0.03559 5.80693e06  0.06444 1.41324e06  
6  2.47975  0.03561 5.83718e06  0.06445 1.41869e06  
3  24  2.50294  0.03559 5.80783e06  0.06444 1.41341e06 
96  2.51280  0.03558 5.80954e06  0.06444 1.41362e06 
6  0.07107  0.07214  0.07323  0.07397 2.04990e08  
1  24  0.07050  0.07167  0.07283  0.07362 3.65907e08 
96  0.07033  0.07152  0.07269  0.07351 3.91936e08  
6  0.07030  0.07151  0.07268  0.07349 3.69824e08  
2  24  0.07027  0.07147  0.07265  0.07347 4.03121e08 
96  0.07026  0.07147  0.07265  0.07347 4.02324e08  
6  0.07027  0.07147  0.07265  0.07347 4.02573e08  
3  24  0.07026  0.07147  0.07265  0.07347 4.02248e08 
96  0.07026  0.07147  0.07265  0.07347 4.02332e08 
The results of solving the spectral problem (23) for the first eigenvalues , on different computational grids using different finite element approximations are shown in Table 2, 3. The eigenvalues , , of the spectral problem (23) are complex with small imaginary parts, the eigenvalues , 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 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 using finite elements of degree . The main eigenfunctions for the second group of and delayed neutrons source are shown in Fig.9.
for problem (14)  for problem (16)  

1  2.51280117966  2.51280117972 
2,3  0.0355815000364 5.80954455861e06  0.0355815000365 5.80954421646e06 
4,5  0.0644427013767 1.41362187449e06  0.0644427013767 1.41362190730e06 
6  0.0702618501639  0.0702618501639 
7  0.0714652882224  0.0714652882164 
8  0.0726456060606  0.0726456060606 
9,10  0.0734708921578 4.02332269037e08  0.0734708921578 4.02332146248e08 
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 (), 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 for the first 10 eigenfunctions. For convenience of comparison, the eigenfunctions are normalized in :
The maximum nonorthogonality (for ) 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.
\  1  2  3  4  5  6  7  8  9  10 

1  1.0e00  1.3e08  2.2e08  3.8e08  9.8e09  1.8e09  1.0e02  3.2e09  2.2e08  1.6e09 
2  1.3e08  1.0e00  1.6e08  1.6e08  1.4e08  4.1e08  1.2e09  2.0e07  3.1e03  7.5e03 
3  2.2e08  1.6e08  1.0e00  9.8e09  1.1e08  1.8e08  1.1e08  3.3e08  7.5e03  3.1e03 
4  3.8e08  1.6e08  9.8e09  1.0e00  3.9e10  1.1e08  1.4e08  4.0e09  3.0e09  1.1e08 
5  9.8e09  1.4e08  1.1e08  3.9e10  1.0e00  2.9e09  1.6e08  1.9e08  6.3e09  6.3e09 
6  1.8e09  4.1e08  1.8e08  1.1e08  2.9e09  1.0e00  4.2e09  5.6e03  4.1e08  1.2e07 
7  1.0e02  1.2e09  1.1e08  1.4e08  1.6e08  4.2e09  1.0e00  2.1e09  1.8e08  8.0e09 
8  3.2e09  2.0e07  3.3e08  4.0e09  1.9e08  5.6e03  2.1e09  1.0e00  5.2e08  2.3e07 
9  2.2e08  3.1e03  7.5e03  3.0e09  6.3e09  4.1e08  1.8e08  5.2e08  1.0e00  5.5e07 
10  1.6e09  7.5e03  3.1e03  1.1e08  6.3e09  1.2e07  8.0e09  2.3e07  5.5e07  1.0e00 
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
(24) 
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
Here is the first mode of the supercritical state. We consider the problem with the transition from this supercritical state at to the subcritical state.
The subcritical stage is characterized by a 15% increase in the coefficient for material 4 in the VVER1000 test diffusion constants (see Table 1). Thus, the dynamics of the reactor is as follows:
The initial state is characterized by specifying the initial conditions at as
(25) 
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.
6  0.03602  0.05760 1.49652e06  0.06890 4.92606e07  

1  24  0.02656  0.05502 2.06007e06  0.06804 1.01253e06 
96  0.02276  0.05411 2.16813e06  0.06774 1.03843e06  
6  0.02250  0.05404 1.81823e06  0.06772 8.73562e07  
2  24  0.02144  0.05380 2.19400e06  0.06765 1.04253e06 
96  0.02125  0.05376 2.20812e06  0.06763 1.04715e06  
6  0.02139  0.05379 2.22579e06  0.06764 1.05369e06  
3  24  0.02124  0.05376 2.20883e06  0.06763 1.04736e06 
96  0.02122  0.05376 2.20951e06  0.06763 1.04756e06 
6  0.07276  0.07363  0.07369  0.07466 2.47162e08  
1  24  0.07222  0.07316  0.07329  0.07429 1.08814e08 
96  0.07204  0.07301  0.07316  0.07417 1.38093e08  
6  0.07203  0.07300  0.07315  0.07416 1.26708e08  
2  24  0.07199  0.07296  0.07312  0.07413 1.50527e08 
96  0.07198  0.07295  0.07312  0.07413 1.49196e08  
6  0.07198  0.07295  0.07312  0.07413 1.52256e08  
3  24  0.07198  0.07295  0.07312  0.07413 1.49141e08 
96  0.07198  0.07295  0.07311  0.07413 1.49178e08 
For an approximate solution, we use the following formulation:
(26) 
where the coefficients are calculated according to the given initial condition (24). These coefficients for 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 to the function . 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 . 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 for material 4 in the diffusion constants in the upper and lower half of the reactor crosssection (see Fig.2). Now let the reactor dynamics corresponds to the following transformation
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.
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 
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 
The coefficients , 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 ) are performed at .
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 is used (see details in Avvakumov et al. (2016)). The dynamics of the neutron power of the nuclear reactor and the delayed neutrons source 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
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 eigenvalue spectral problem.
Numerical simulation of reactor nonstationary 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 (ofline calculation) eigenfunctions and eigenvalues of the eigenvalue spectral problem. Fast determination of dominant modes and calculation of the reactor neutron flux at selected times are based on online calculation. The classical Lagrange finite elements are used for the spatial approximation. Accuracy control is performed using condensed grids. Spectral problems are solved numerically using welldeveloped free software SLEPc.
Test calculations are made in twodimensional analysis using a twogroup 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 nonuniform change in the absorbing material properties (in two halves over the reactor crosssection).
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 VVER1000 test problem.
Acknowledgements
This work are supported by the Russian Foundation for Basic Research (# 160801215)
and by the grant of the Russian Federation Government
(# 14.Y26.31.0013).
References
 Ascher (2008) Ascher, U. M., 2008. Numerical Methods for Evolutionary Differential Equations. Society for Industrial Mathematics.
 Avvakumov et al. (2017) Avvakumov, A. V., Strizhov, V. F., Vabishchevich, P. N., Vasilev, A. O., 2017. Spectral properties of dynamic processes in a nuclear reactor. Annals of Nuclear Energy 99, 68–79.
 Avvakumov et al. (2016) Avvakumov, A. V., Vabishchevich, P. N., Vasilev, A. O., Strizhov, V. F., 2016. Numerical modelling neutron diffusion unsteady problems. Mathematical Models and Computer Simulations (submitted).
 Baudron et al. (2014) Baudron, A.M., Lautard, J.J., Maday, Y., Riahi, M. K., Salomon, J., 2014. Parareal in time 3D numerical solver for the LWR benchmark neutron diffusion transient model. Journal of Computational Physics 279, 67–79.
 Bell and Glasstone (1970) Bell, G. I., Glasstone, S., 1970. Nuclear Reactor Theory. Van Nostrand Reinhold Company.
 Björck (1996) Björck, A., 1996. Numerical Methods for Least Squares Problems. Society for Industrial and Applied Mathematics.
 Brenner and Scott (2008) Brenner, S. C., Scott, R., 2008. The Mathematical Theory of Finite Element Methods. Springer.
 Brezinski (1991) Brezinski, C., 1991. Biorthogonality and Its Applications to Numerical Analysis. CRC Press.
 Butcher (2008) Butcher, J. C., 2008. Numerical Methods for Ordinary Differential Equations. Wiley.
 Chao and Shatilla (1995) Chao, Y. A., Shatilla, Y. A., 1995. Conformal mapping and hexagonal nodal methodsII: Implementation in the ANCH Code. Nuclear Science and Engineering 121, 210–225.
 Cho (2005) Cho, N. Z., 2005. Fundamentals and recent developments of reactor physics methods. Nuclear Engineering and Technology 37 (1), 25–78.
 Chou et al. (1990) Chou, H. P., Lu, J. R., Chang, M. B., 1990. A threedimensional spacetime model and its use in pressurized water reactor rod ejection analyses. Nuclear Technology 90 (2), 142–154.
 Dahmani et al. (2001) Dahmani, M., Baudron, A. M., Lautard, J. J., Erradi, L., 2001. A 3D nodal mixed dual method for nuclear reactor kinetics with improved quasistatic model and a semiimplicit scheme to solve the precursor equations. Annals of Nuclear Energy 28 (8), 805–824.
 Dodds Jr (1976) Dodds Jr, H. L., 1976. Accuracy of the quasistatic method for twodimensional thermal reactor transients with feedback. Nuclear Science and Engineering 59 (3), 271–276.
 Duderstadt and Hamilton (1976) Duderstadt, J. J., Hamilton, L. J., 1976. Nuclear Reactor Analysis. Wiley.
 Dugan et al. (2016) Dugan, K., Zmijarevic, I., Sanchez, R., 2016. Crosssection homogenization for reactivityinduced transient calculations. Journal of Computational and Theoretical Transport 45 (6), 425–441.
 Dulla et al. (2008) Dulla, S., Mund, E. H., Ravetto, P., 2008. The quasistatic method revisited. Progress in Nuclear Energy 50 (8), 908–920.
 Gear (1971) Gear, C. W., 1971. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice Hall, NJ.
 Gill and Azmy (2011) Gill, D. F., Azmy, Y. Y., 2011. Newton’s method for solving eigenvalue problems in neutron diffusion theory. Nuclear Science and Engineering 167 (2), 141–153.
 Ginestar et al. (2002) Ginestar, D., Miro, R., Verdu, G., Hennig, D., 2002. A transient modal analysis of a BWR instability event. Journal of Nuclear Science and Technology 39 (5), 554–563.
 Goluoglu and Dodds (2001) Goluoglu, S., Dodds, H. L., 2001. A timedependent, threedimensional neutron transport methodology. Nuclear Science and Engineering 139 (3), 248–261.
 GonzálezPintor et al. (2009) GonzálezPintor, S., Ginestar, D., Verdú, G., 2009. High order finite element method for the lambda modes problem on hexagonal geometry. Annals of Nuclear Energy 36 (9), 1450–1462.
 Grossman and Hennart (2007) Grossman, L. M., Hennart, J.P., 2007. Nodal diffusion methods for spacetime neutron kinetics. Progress in Nuclear Energy 49 (3), 181–216.
 Guérin et al. (2010) Guérin, P., Baudron, A.M., Lautard, J.J., 2010. Domain decomposition methods for the neutron diffusion problem. Mathematics and Computers in Simulation 80 (11), 2159–2167.
 Hairer and Wanner (2010) Hairer, E., Wanner, G., 2010. Solving Ordinary Differential Equations II: Stiff and DifferentialAlgebraic Problems. Springer Verlag.
 Henry (1975) Henry, A. F., 1975. NuclearReactor Analysis. MIT Press.
 Hernandez et al. (2005) Hernandez, V., Roman, J. E., Vidal, V., 2005. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems. ACM Transactions on Mathematical Software (TOMS) 31 (3), 351–362.
 Hernandez et al. (2003) Hernandez, V., Roman, J. E., Vidal, V., Verdu, G., Ginestar, D., 2003. Resolution of the neutron diffusion equation with SLEPc, the Scalable Library for Eigenvalue Problem Computations. In: Nuclear Mathematical and Computational Sciences: A Century in Review, A Century Anew Gatlinburg. American Nuclear Society, pp. 1–10.
 Hetrick (1971) Hetrick, D. L., 1971. Dynamics of Nuclear Reactors. University of Chicago Press.
 Hundsdorfer and Verwer (2003) Hundsdorfer, W. H., Verwer, J. G., 2003. Numerical Solution of TimeDependent AdvectionDiffusionReaction Equations. Springer Verlag.
 Laub (2005) Laub, A. J., 2005. Matrix Analysis for Scientists and Engineers. Society for Industrial and Applied Mathematics.
 Lawrence (1986) Lawrence, R. D., 1986. Progress in nodal methods for the solution of the neutron diffusion and transport equations. Progress in Nuclear Energy 17 (3), 271–301.

LeVeque (2007)
LeVeque, R. J., 2007. Finite Difference Methods for Ordinary and Partial Differential Equations: SteadyState and TimeDependent Problems. Society for Industrial Mathematics.
 Lewis and Miller (1993) Lewis, E. E., Miller, W. F., 1993. Computational Methods of Neutron Transport. American Nuclear Society.
 Luikov (1968) Luikov, A., 1968. Analytical Heat Diffusion Theory. Academic Press.
 Maday and Turinici (2005) Maday, Y., Turinici, G., 2005. The parareal in time iterative solver: a further direction to parallel implementation. In: Domain Decomposition Methods in Science and Engineering. Springer, pp. 441–448.
 Marchuk and Lebedev (1986) Marchuk, G. I., Lebedev, V. I., 1986. Numerical Methods in the Theory of Neutron Transport. Harwood Academic Pub.
 Miró et al. (2002) Miró, R., Ginestar, D., Verdú, G., Hennig, D., 2002. A nodal modal method for the neutron diffusion equation. application to BWR instabilities analysis. Annals of Nuclear Energy 29 (10), 1171–1194.
 Modak and Gupta (2007) Modak, R., Gupta, A., 2007. A scheme for the evaluation of dominant timeeigenvalues of a nuclear reactor. Annals of Nuclear Energy 34 (3), 213–221.
 Ortega (1987) Ortega, J. M., 1987. Matrix Theory: A Second Course. Springer.
 Quarteroni and Valli (2008) Quarteroni, A., Valli, A., 2008. Numerical Approximation of Partial Differential Equations. Springer.
 Saad (2003) Saad, Y., 2003. Iterative Methods for Sparse Linear Systems. Society for Industrial Mathematics.
 Saad (2011) Saad, Y., 2011. Numerical Methods for Large Eigenvalue Problems. SIAM.
 Samarskii (2001) Samarskii, A. A., 2001. The Theory of Difference Schemes. Marcel Dekker, New York.
 Samarskii and Vabishchevich (1996) Samarskii, A. A., Vabishchevich, P. N., 1996. Computational Heat Transfer. Wiley.
 Sanchez (2009) Sanchez, R., 2009. Assembly homogenization techniques for core calculations. Progress in Nuclear Energy 51 (1), 14–31.
 Smith (1979) Smith, K. S., 1979. An analytic nodal method for solving the twogroup, multidimensional, static and transient neutron diffusion equations. Ph.D. thesis, Massachusetts Institute of Technology.
 Stacey (1967) Stacey, W. M., 1967. Modal Approximations: Theory and an Application to Reactor Physics. The MIT Press.
 Stacey (1969) Stacey, W. M., 1969. SpaceTime Nuclear Reactor Kinetics. Academic Press.
 Stacey (2007) Stacey, W. M., 2007. Nuclear Reactor Physics. Wiley.
 Stewart (2001) Stewart, G. W., 2001. A Krylov–Schur algorithm for large eigenproblems. SIAM Journal on Matrix Analysis and Applications 23 (3), 601–614.
 Stewart (1976) Stewart, H. B., 1976. Spectral theory of heterogeneous diffusion systems. Journal of Mathematical Analysis and Applications 54 (1), 59–78.
 Sutton and Aviles (1996) Sutton, T. M., Aviles, B. N., 1996. Diffusion theory methods for spatial kinetics calculations. Progress in Nuclear Energy 30 (2), 119–182.
 Toselli and Widlund (2005) Toselli, A., Widlund, O., 2005. Domain Decomposition Methods – Algorithms and Theory. Springer.
 Vabishchevich (2012) Vabishchevich, P. N., 2012. SMstability of operatordifference schemes. Computational Mathematics and Mathematical Physics 52 (6), 887–894.
 Vabishchevich (2014) Vabishchevich, P. N., 2014. Additive OperatorDifference Schemes: Splitting Schemes. de Gruyter.
 Verdú and Ginestar (2014) Verdú, G., Ginestar, D., 2014. Modal decomposition method for BWR stability analysis using Alphamodes. Annals of Nuclear Energy 67, 31–40.
 Verdu et al. (2010) Verdu, G., Ginestar, D., Roman, J., Vidal, V., 2010. 3D alpha modes of a nuclear power reactor. Journal of Nuclear Science and Technology 47 (5), 501–514.
 Verdú et al. (1998) Verdú, G., Ginestar, D., Vidal, V., Miró, R., 1998. Modal decomposition method for BWR stability analysis. Journal of Nuclear Science and Technology 35 (8), 538–546.
 VidalFerrandiz et al. (2014) VidalFerrandiz, A., Fayez, R., Ginestar, D., Verdú, G., 2014. Solution of the lambda modes problem of a nuclear power reactor using an h–p finite element method. Annals of Nuclear Energy 72, 338–349.
 Willert et al. (2014) Willert, J., Park, H., Knoll, D. A., 2014. A comparison of acceleration methods for solving the neutron transport keigenvalue problem. Journal of Computational Physics 274, 681–694.