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

by   Alexander V. Avvakumov, et al.

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.


page 12

page 14

page 15

page 16

page 17


Change point inference in ergodic diffusion processes based on high frequency data

We deal with the change point problem in ergodic diffusion processes bas...

Numerical modeling of neutron transport in SP3 approximation by finite element method

The SP3 approximation of the neutron transport equation allows improving...

Modeling and simulation of nuclear architecture reorganization process using a phase field approach

We develop a special phase field/diffusive interface method to model the...

Mathematical modelling of nuclear medicine data

Positron Emission Tomography using 2-[18F]-2deoxy-D-glucose as radiotrac...

Global existence of classical solutions and numerical simulations of a cancer invasion model

In this paper, we study a cancer invasion model both theoretically and n...

Efficient Low-Order Approximation of First-Passage Time Distributions

We consider the problem of computing first-passage time distributions fo...

Neural ODE and DAE Modules for Power System Dynamic Modeling

The time-domain simulation is the fundamental tool for power system tran...

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 (-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 is used. It is defined as the main eigenvalue of the time-eigenvalue (-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 , 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, -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á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 -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 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 -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 () with boundary . The neutron transport is described by the system of equations:


Here — neutron flux of group at point and time , — number of energy groups, — effective velocity of neutrons in the group , — diffusion coefficient, — removal cross-section, — scattering cross-section from group to group , — effective fraction of delayed neutrons, , — spectra of prompt and delayed neutrons, — generation cross-section 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:


where is a fraction of delayed neutrons of -type, and

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

The albedo-type conditions are set at the boundary of the area :


where — outer normal to the boundary . Initial state is defined in the following manner:


Let’s write the boundary problem (1)–(4

) in operator form. The vectors

, and matrices are defined as follows:


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 first-order equation of evolution:


The Cauchy problem is formulated for equations (5) when


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 two-layer or three-layer 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 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).




Figure 1: State change scheme.

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 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 -eigenvalue problem. Finite-element approximation in space is used.

In a separate stage the following system of equations is considered


supplemented by the corresponding initial conditions


Let’s . Rewrite the system of equations (7) as


with constants


is the identity matrix. From (

8) one can obtain


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:


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


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:


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


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 time-scale 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 non-orthogonal functions consists in using the biorthogonal system of functions (Henry, 1975; Brezinski, 1991). Consider the spectral problem adjoint to (14)


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


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)):


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


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


for all .

Further, we must pass from the continuous variational problem (20) to the discrete problem. We introduce finite-dimensional finite element spaces , . The discrete variational problem is formulated as follows: find , such that


for all . 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, 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 ( 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.

Figure 2: Geometrcial model of the VVER-1000 reactor core.
Figure 3: Discretization of assembly into 6, 24 and 96 finite elements.

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 two-group approximation is used taken into account delayed neutrons:

Material 1 2 3 4 5
1.38320e-0 1.38299e-0 1.39522e-0 1.39446e-0 1.39506e-0
3.86277e-1 3.89403e-1 3.86225e-1 3.87723e-1 3.84492e-1
2.48836e-2 2.62865e-2 2.45662e-2 2.60117e-2 2.46141e-2
6.73049e-2 8.10328e-2 8.44801e-1 9.89671e-2 8.93878e-2
1.64977e-2 1.47315e-2 1.56219e-2 1.40185e-2 1.54981e-2
4.81619e-3 4.66953e-3 6.04889e-3 5.91507e-3 6.40256e-3
8.46154e-2 8.52264e-2 1.19428e-1 1.20497e-1 1.29281e-1
Table 1: Diffusion neutronics constants for VVER-1000

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 two-group approximation and taking into account the delayed neutrons, one can write


The dominant eigenvalues are searched for at

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

6 -0.22557 0.04241 3.08808e-06 0.06588 4.80449e-07
1 24 -0.82690 0.03777 5.37884e-06 0.06489 1.37315e-06
96 -1.74998 0.03619 5.69002e-06 0.06456 1.40299e-06
6 -2.10154 0.03592 4.96474e-06 0.06452 1.21320e-06
2 24 -2.46601 0.03562 5.78277e-06 0.06445 1.40897e-06
96 -2.50375 0.03559 5.80693e-06 0.06444 1.41324e-06
6 -2.47975 0.03561 5.83718e-06 0.06445 1.41869e-06
3 24 -2.50294 0.03559 5.80783e-06 0.06444 1.41341e-06
96 -2.51280 0.03558 5.80954e-06 0.06444 1.41362e-06
Table 2: Eigenvalues
6 0.07107 0.07214 0.07323 0.07397 2.04990e-08
1 24 0.07050 0.07167 0.07283 0.07362 3.65907e-08
96 0.07033 0.07152 0.07269 0.07351 3.91936e-08
6 0.07030 0.07151 0.07268 0.07349 3.69824e-08
2 24 0.07027 0.07147 0.07265 0.07347 4.03121e-08
96 0.07026 0.07147 0.07265 0.07347 4.02324e-08
6 0.07027 0.07147 0.07265 0.07347 4.02573e-08
3 24 0.07026 0.07147 0.07265 0.07347 4.02248e-08
96 0.07026 0.07147 0.07265 0.07347 4.02332e-08
Table 3: Eigenvalues

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.

Figure 4: The eigenfunction (left) and real part of eigenfunctions (right).
Figure 5: Imaginary part of eigenfunctions (left) and real part of eigenfunctions (right).
Figure 6: Imaginary part of eigenfunctions (left) and eigenfunction (right).
Figure 7: The eigenfunction (left) and eigenfunction (right).
Figure 8: Real part of eigenfunctions (left) and imaginary part of eigenfunctions (right).

Dominant eigenfunctions of the spectral problem (23) for the first group are shown in Figs.48. 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.

Figure 9: The eigenfunction (left) and the eigenfunctions (right).
for problem (14) for problem (16)
1 -2.51280117966 -2.51280117972
2,3 0.0355815000364 5.80954455861e-06 0.0355815000365 5.80954421646e-06
4,5 0.0644427013767 1.41362187449e-06 0.0644427013767 1.41362190730e-06
6 0.0702618501639 0.0702618501639
7 0.0714652882224 0.0714652882164
8 0.0726456060606 0.0726456060606
9,10 0.0734708921578 4.02332269037e-08 0.0734708921578 4.02332146248e-08
Table 4: Eigenvalues for direct and adjoint problems

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 non-orthogonality (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.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 5: Scalar product

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


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 VVER-1000 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


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.49652e-06 0.06890 4.92606e-07
1 24 0.02656 0.05502 2.06007e-06 0.06804 1.01253e-06
96 0.02276 0.05411 2.16813e-06 0.06774 1.03843e-06
6 0.02250 0.05404 1.81823e-06 0.06772 8.73562e-07
2 24 0.02144 0.05380 2.19400e-06 0.06765 1.04253e-06
96 0.02125 0.05376 2.20812e-06 0.06763 1.04715e-06
6 0.02139 0.05379 2.22579e-06 0.06764 1.05369e-06
3 24 0.02124 0.05376 2.20883e-06 0.06763 1.04736e-06
96 0.02122 0.05376 2.20951e-06 0.06763 1.04756e-06
Table 6: Subcritical state:
6 0.07276 0.07363 0.07369 0.07466 2.47162e-08
1 24 0.07222 0.07316 0.07329 0.07429 1.08814e-08
96 0.07204 0.07301 0.07316 0.07417 1.38093e-08
6 0.07203 0.07300 0.07315 0.07416 1.26708e-08
2 24 0.07199 0.07296 0.07312 0.07413 1.50527e-08
96 0.07198 0.07295 0.07312 0.07413 1.49196e-08
6 0.07198 0.07295 0.07312 0.07413 1.52256e-08
3 24 0.07198 0.07295 0.07312 0.07413 1.49141e-08
96 0.07198 0.07295 0.07311 0.07413 1.49178e-08
Table 7: Subcritical state:

For an approximate solution, we use the following formulation:


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.

Figure 10: Approximate solution coefficients (26).

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.

1 2 a b c
Figure 11: Function (line 1) and function (line 2): a — neutron flux of group 1, b — neutron flux of group 2, c — delayed neutrons source.

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 cross-section (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
Table 8: Subcritical asymmetric state:
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
Table 9: Subcritical asymmetric state:

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 .

Figure 12: Approximate solution coefficients (26) — asymmetric perturbation.
1 2 a b c
Figure 13: Function (line 1) and function (line 2) for asymmetric perturbation: a — neutron flux of group 1, b — neutron flux of group 2, c — delayed neutrons source.

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.

Figure 14: Fast stage of reactor state: neutronic power.
Figure 15: Slow stage of reactor state: neutronic power.
Figure 16: Slow stage of reactor state: delayed neutrons sourse.
Figure 17: Slow stage of reactor state for the asymmetric perturbation: neutronic power.
Figure 18: Slow stage of reactor state for the asymmetric perturbation: delayed neutrons source.

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


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


  • 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 methods-II: Implementation in the ANC-H 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 three-dimensional space-time 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 semi-implicit 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 two-dimensional 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. Cross-section homogenization for reactivity-induced 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 quasi-static 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 time-dependent, three-dimensional neutron transport methodology. Nuclear Science and Engineering 139 (3), 248–261.
  • González-Pintor et al. (2009) González-Pintor, 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 space-time 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 Differential-Algebraic Problems. Springer Verlag.
  • Henry (1975) Henry, A. F., 1975. Nuclear-Reactor 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 Time-Dependent Advection-Diffusion-Reaction 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: Steady-State and Time-Dependent 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 time-eigenvalues 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 two-group, 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. Space-Time 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. SM-stability of operator-difference schemes. Computational Mathematics and Mathematical Physics 52 (6), 887–894.
  • Vabishchevich (2014) Vabishchevich, P. N., 2014. Additive Operator-Difference Schemes: Splitting Schemes. de Gruyter.
  • Verdú and Ginestar (2014) Verdú, G., Ginestar, D., 2014. Modal decomposition method for BWR stability analysis using Alpha-modes. 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.
  • Vidal-Ferrandiz et al. (2014) Vidal-Ferrandiz, 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 k-eigenvalue problem. Journal of Computational Physics 274, 681–694.