A numerical method for the solution of the kinetic equation of a reactor

Abstract

The phase density N( r, v, t) of neutrons in a reactor obeys Boltzmann's linearized equation. The solution of the equation is complicated and in many applied problems we can be satisfied by knowing how the total number n( t) of neutrons in the reactor depends on time. Under certain conditions, the phase density of the neutrons can be found in the form N = N 0( r, v) n( t). We then obtain the following integro-differential equation for n( t) (see [1], [2]) from the Boltzmann equation: dn dt =(K (t)− β l) n(t)+ ∑ j β jλ j l ∫ −∞ t n(t′)e −λ j(t−t′) dt′ . where K(t) = (k eff − 1) lk eff , l is the average life of a neutron in the reactor, λ j is the decay constant of predecessors of the jth group of delayed neutrons, β j is the effective fraction of delayed neutrons in the jth group relative to the total number of fission neutrons, and β = ∑ j β j . The quantities k eff , β j and l depend, in general, on time. However, in many cases the fact that β j and l are not constant can be ignored [1]. The work which follows is concerned with a numerical method for the solution of equation (1) where β j = const, l = const and K( t) is an arbitrary function from a sufficiently wide class. Since (1) is equivalent to a system of ordinary linear differential equations no difficulties will arise as far as the existence and uniqueness of its solution are concerned. In finding a numerical method to solve equation (1) we have been principally interested in the possibilities of controlling the accuracy of the computation internally, and of choosing the optimal size of the time step. We have chosen the latter in the following way. Suppose that A = A( h) is some parameter to be worked out in the approximate solution of the problem on a step of length h. We stipulate that h shall satisfy the conditions ε − ⩽ A ( h) ⩽ ε +, where ge − and ge su+ are pre-set constants. Since A( h) depends on the behaviour of the reactivity K( t) and the solution n( t) within the step, h is, generally speaking, a variable function of the step number. If A( h) is related in some way to the accuracy of the computation in the given step, then any change in the step length will automatically keep the error in the approximate method introduced during each step within constant limits, depending on the values of the constants ge + and ge −. It might happen for some values of ge + and ge − that the step h is near its optimal value h 0, for which the calculation error of the unknown function at the end of the step will be least (or will vanish altogether). In the latter case a small change in the step length can lead to a change in the sign of the error in this step, thus preventing the build-up of errors from step to step. This procedure for choosing the step length is easily programmed. Actual calculations have confirmed the considerations discussed above to be valid, and have shown that a very high degree of accuracy in the solution of equation (1) with variable reactivity can be obtained with a moderate increase in machine time.