Solution of two-dimensional space–time multigroup reactor kinetics equations by generalized Padé and cut-product approximations

Abstract

New methods based on the class of Padé and cut-product approximations are applied to the solution of the multigroup diffusion theory reactor kinetics equations in two space dimensions. The methods are shown to be consistent and numerically stable. The stability of the developed approximations is studied and it is in general A(α)-stable, where 0 ⩽ α ⩽ π/2 and therefore they are quite efficient in the presence of extreme stiffness, as shown by numerical results for some typical test cases. The system of stiff coupled differential equations is represented in a matrix form and due to the particular structure of this matrix, new algorithm based on the back-substitution is proposed for analytical inversion of the “A” matrix. The method has been tested for different benchmark problems. The results indicated that this method can allow much larger time steps and thus save much computational time as compared to the other conventional methods.

Introduction

An approximate solution of the space-dependent multigroup diffusion theory kinetics equations is particularly difficult to obtain by standard finite difference techniques, especially in more than one space dimension. For reactivities less than prompt critical, the reactor kinetics equations are “stiff” equations from a mathematical point of view. Stiffness in the above connection will mean that the equations have time constants spanning a wide range of values. For the kinetics equations this is due to the physical fact that the decay constants of the precursors are much smaller than the time constants characteristic of either the neutron slowing down or the neutron diffusion process. These precursors retard the response of the reactor to a perturbation from a critical state, a well-known fact which facilitates the control of a nuclear reactor. The prediction of this response through numerical solutions of the kinetics equations is made difficult by the time constants of the neutron diffusion, but the response of the precursors greatly lengthens the time interval during which the neutron flux adjusts to a new shape. An unacceptably large number of time steps are necessary to predict this most interesting region of the transient response.

It has been known for some time that space-dependent calculations of reactor transient are required in certain situations in order to accurately predict the power behavior and the spatial distribution of neutron flux. For example, one of the most important reactor analysis safety issues involves the postulated ejection of the highest reactivity worth control rod from the core of a pressurized water reactor. Such a strongly localized perturbation causes the time-dependent neutron flux to become a function that can not be separated into the product of a time-dependent amplitude function. In this situation the point kinetics model is invalid and incapable of predicting the detailed behavior of the transient.

The motivation behind the development of the methods for solving the energy-, space-, and time-dependent kinetic equations is not only the challenge of developing a method for solving a large set of coupled partial differential equations, but also a real need to predict the performance and assess the safety of large commercial reactors, both those presently operating and those being designed for the future.

The point kinetic approximation, which assumes complete separability in time and space, is valid for the tightly coupled small cores where the spatial distribution is not sensitive to local changes in the reactor properties. Changes in the nuclear properties which occur during a transient do, of course, cause changes in the space- and energy-distribution of neutrons within a reactor. The fundamental source of error in the basic point kinetic model is the failure to account for these changes in neutron flux distribution. Numerical experiments by Yasinsky and Henery (1965) have clearly demonstrated the inadequacy of the point kinetic solution and the need for more sophisticated methods. Since then many comparative analysis between the point kinetics and various spatial approximations were performed using different approximations for the full space solutions, resulting in much better knowledge on the range of applicability of the various approximations.

Various methods have been developed by many investigators to find numerical solution to the multigroup space–time dependent neutron diffusion equations. These methods are usually classified as finite difference methods (Chou et al., 1990, Urku and Christenson, 1994, Jagannathan, 1985a, Jagannathan, 1985b), finite element methods (Koclas, 1998), synthesis or nodal expansion methods (Kaplan, 1962, Dougherty and Shen, 1962, Kaplan et al., 1964, Yasinsky, 1967, Yasinsky and Kaplan, 1967, Stacey, 1968, Jagannathan, 1985a, Jagannathan, 1985b, Crouzet, 1996), nodal or coarse mesh methods (Ott, 1969, Dodds, 1976, Camiciola et al., 1986, Favorite and Stacey, 1997a, Favorite and Stacey, 1997b, Koclas et al., 1997, Dahmani et al., 2001) which later evolved into response matrix methods and quasistatic and adiabatic methods.

A broad class of methods, called alternating semi-implicit method, for the solution of reactor kinetics equations are examined by Reed and Hansen (1970). This class of methods includes such well known schemes as the alternating direction implicit, or ADI, method and the alternating direction explicit, or ADE, method, Varga (1962). Weight et al. (1971) have introduced an approximate solution of the multigroup neutron diffusion kinetics equations with delayed neutrons in two-dimensional geometry by matrix splitting methods based on an alternating-direction implicit (ADI) scheme. However, this conventional ADI method is unstable for heterogeneous problems unless extremely small time steps are used. Recently, the suitability of the ADI method for parallel computation has been developed by Chen et al. (1992). This method becomes extremely attractive for parallel computer applications. Since, it is based on the solution of a system of independent block-tridiagonal matrix equations that can be solved in parallel by improving the stability of the ADI method. However, it still suffered from the computational time for both homogeneous and heterogeneous two-group problems.

For many reasons, it is sometimes convenient to approximate the exponential matrix by a rational approximation. Two techniques are considered in this work for accomplishing this goal for the ‘approximation problem’ of network synthesis and automatic control.

The present work is concerned with the development of an approximate solution of the multigroup neutron diffusion kinetics equations with delayed neutrons in two-dimension geometry based on the class of Padé and cut-product approximations. These methods are applicable to a general class of problems and the amount of computer time, the stability of approximations, and the accuracy in between and/or with the conventional methods are discussed.

In Section 2 the multigroup neutron diffusion kinetics equations with delayed neutrons are investigated and presented in their matrix form. In Section 3, the basic method is introduced and analyzed from a mathematical point of view; stability of the different approximations of the method are presented and discussed. Finally, numerical solutions for several sample problems have been obtained using these methods, as well as using other conventional methods. These solutions are compared in Section 4.