General solution of reactor kinetic equations for one group of delayed neutrons and time dependent reactivities

Journal of NuclearEnergy, Vol.26,pp.213to230.Persamon Press1972.Printed in Northern Ireland

GENERAL SOLUTION OF REACTOR KINETIC EQUATIONS FOR ONE GROUP OF DELAYED NEUTRONS AND TIME DEPENDENT REACTIVITIES G.

M. SANDQUIST

University of Utah, Mechanical Engineering Department, U.S.A. (Received

16 December

Salt Lake City, Utah,

1971)

Abstract-A new method of obtaining analytic solutions of the space independent reactor kinetic equations for arbitrary time dependent reactivities is presented. The method is based upon transforming the kinetic equations with one group of delayed neutrons into a single Ricatti equation for the system. This single equation permits a variety of useful approximate solutions which are generally considerably more accurate than approximate solutions obtained for the standard kinetic equations. General perturbation solutions of the single equation for the system are found and are shown to be valid even for very large changes of the dependent variables. Specific solutions for step, impulse, ramp sinusoidal and generalized reactivities are found. The solution for sinusoidal reactivities confirms that unstable exponential flux growth occurs due to the presence of delayed neutrons, while without delayed neutrons the solution is oscillatory but bounded. Finally the exact solution is found for the system when the reactivity is an analytic function. INTRODUCTION

independent reactor kinetic equations which are the basis of most analytic kinetic studies of reactor systems have been studied at length and subjected to numerous mathematical and simulation methods and techniques. However, in fact, the repertory of useful and satisfactory analytic solutions for the general class of reactivity functions that are of interest to kineticists, reactor operators and those responsible for safety and hazard analysis is rather limited and incomplete even for the linear system. Because of the analytical difficulties computer simulation (KEEPIN and Cox, 1960) of the system equations is generally employed for system studies. However, there is a real and significant loss and masking of information and insight into such areas as general characteristics, singularities, stability, parameter dependence and sensitivity, and optimization possibilities when computer methods are employed. For a complete study of the kinetic behavior of the system, all groups of delayed neutrons and reactivity feedback terms as well as spatial effects should be considered. However, we shall consider only time-dependent reactivity functions with no dependent variable feedback. Also for simple analysis and with generally good success, one group of delayed neutrons with appropriate spatial averaging usually suffices to adequately describe the time behavior of the system. By appropriately selecting and weighing this one effective group of delayed neutrons (KIRCHENMAYER, 1960; SKINNER and COHEN, 1960), a broad range of time-dependent reactivity functions can be treated. A complete survey of previous work in the area of space independent kinetics would be voluminous and beyond the scope of this study. Of particular interest, however, is the work of SMETS (1958) and AKCASU (1958). Also KEEPIN (1965) and others have provided calculations and measurements of the parameter values associated with the kinetic equations. The primary purpose of this study is to present a scheme for reducing the system equations for the neutron density and the delayed neutron precursor density for one THE SPACE

1

213

214

G. M. SANDQUIST

delayed group into a single first order non-linear differential equation by a transformation involving the ratio of the neutron and precursor density. This variable transformation is bounded for all time for any bounded reactivity input and approximate solution methods such as perturbation techniques are generally valid even for very large changes of the neutron and precursor density. SYSTEM

DESCRIPTION

The space independent nuclear reactor kinetic equations (with no nonfission neutron sources) which can be derived from the neutron transport equation under reasonable assumptions (HENRY, 1958 ; HANSENand MAIER, 1960) are

C

(~-l)N+Cb~c~

;=A

2=

1,(N - C,),

i 1,

(l-1) (1*a

(i=l,...,M),

where N and Ci are the normalized neutron and i’th delayed neutron precursor densities, p(t) is the reactivity in dollars, &(sec-l) and A[= fi/l(sec-l)] are the i’th delayed neutron decay constant and prompt neutron time constant and b,(= &/,!I)is the effective fractional fission yield for the i’th delayed neutron precursor. Equilibrium values for equations (1) are N = Ci > 0 for all i and p = 0 and also for any value of p when N = C, = 0 for all i. The first set of equilibrium values (p = 0) corresponds to equilibrium operation at any arbitrary power level while the second set of equilibrium values (N = C, = 0) corresponds to reactor shutdown. It is possible to reduce this system of kinetic equations of order M + 1 to an M’th order system by the introduction of a new set of dependent variables given by (for N # 0) (i = 1, . . . , M). 5i = CJN, (2) So then equations (1) become i

fi =

a,(1- EJ - A

(p- 1 + 2j bjCj)

Ei

fori,j= 1, . . . ,M. The result of this change of variable is to reduce the order of the system but the resulting equations form a non-linear system of Ricatti type equations which are difficult to solve analytically although certain solutions are known. However, if we assume that one group of delayed neutrons can provide an acceptable description of the behavior of the system, then equations (3) reduce to the following first order non-linear equation of the Ricatti type (for M = 1). $5

= A(1 - 8 - A(p - 1 + 05

(4.1)

As an alternative formulation the variable definition 17= l/5 = N/C (for M = 1 and C # 0) yields as an alternative system description

z ?,j =A(1 -

r])?j + A[@ - l)r + 11,

(4.2)

General solutions of reactor kinetic equations for one group of delayed neutrons

215

where from physical requirements we have that 7 2 0 and 6 2 0 since N 2 0 and all Ci 2 0. It is possible to obtain additional equations of the Ricatti type for the one delayed neutron group system by eliminating through differentiation one of the dependent variables to obtain a single second order equation. Then with a well-known variable transformation (INCE, 1956), the second-order equation can be reduced to a Ricatti type equation. SMETS(1958) has done this and obtained solutions for the system for a certain limited class of reactivities under the assumption that the neutron density time derivative and reactivity are differentiable. However, of the many Ricatti type equation forms that could be formulated from equation (l), the particular variable transformation employed to obtain equation (4) (7 = l/E = N/C) is of particular significance and value since we shall find that 7 and [ are bounded for all time for any bounded reactivity function and in general tend to vary over a smaller range of values than do the corresponding variables N and C individually for a given reactivity function. Furthermore, no assumption of differentiability of the reactivity is made here which is necessary for certain classes of reactivity functions. Either equations (4.1) or (4.2) may be used as the defining equation for a system model with one effective group of delayed neutrons where the system reactivity is a function of time. We shall find that although equations (4) are non-linear, they admit many useful and accurate approximate solutions and even exact solutions under certain assumptions. Furthermore, we should observe that when we have any solution for equations (4) we may obtain the associated solution for N and C by noting from equation (1.2) that the precursor density for one group of delayed neutrons can be expressed as i In C = A(7j - 1) = 1(1/l - 1) which has as a solution given that C(0) = 1 C(t) = exp

ii (7 - 1) dt [s,’

I- [s, - exp

1 ‘(l/l - 1) dt

1

The neutron density is then found by recalling from previous definitions that N = yC = C/t so with N(0) = 1 then N(t)=

qexp kjO’qdt

- it)

(6.1)

and N(t) = 1 exp ( A]o’;dt-,t).

(6.2)

If exact solutions are found for 7 and 5 or approximate solutions for which the approximations are equivalent in equations (4) then the resulting solutions given by equations (6) will be identical. However, in general if nonequivalent approximations are made the resulting solutions for the neutron density given by equations (6) will not be identical and the more accurate solution for a given reactivity function should be used. Note that the first time derivative of the neutron density at zero time is identical for

G. M. SANDQUIST

216

both equations (6) and is given by dN dt rzO= AP(O) Also, observe that the solution for the system without delayed neutrons is obtained by setting L = 0 and p = 0 so that C = 1 for all time and N(f) = r(t) = l/[(t). With ;Z= 0, equation (4.1) becomes a Bernoulli equation and is easily solved while equation (4.2) is linear and also easily solved. For a step reactivity input p,,, equations (4) are separable and the resulting solution with r](O) = t(O) = 1 is

770)= l/W = 1 +

m2

(1 - exp [Cm2- m&l

--j-m

2

-

exp

-

(7)

[(m2- m&l

ml

where m2 Z -

m1

4 31 -&2/(a” +

4p&),

a = A(1 - po) + A We find that for all real values of p,, that m2 - ml is positive real so we have as the asymptotic values for 7 and 5 from equation (7) as t --f co for all real values of p,, that qm = l/Em = & ; (PO or

1) +

1 +

2 (PO -

II2 + 2;

(PO +

1) +

11

equivalently 2PO

qm = l/E, = 1 + 1-po+x+

;I

2

(PO -

II2 + 211 (PO +

\i[

1) +

1

( II ;i

a particularly useful form when evaluating the asymptotic values for p,, < 0 andA >A. Figure 1 displays the asymptotic values qa, and 6, as functions of the time constant ratio A/A for various step insertions of reactivity in dollars. It is significant to observe that r(t) and t(t) are bounded for all time for all finite valued functions of reactivity even when the corresponding values of neutron and precursor density become unbounded. Thus, in general, approximate solutions of equations (4) will exhibit greater accuracy and parameter range capacity than approximate solutions obtained for equations (1) and alternate formulations for one group of delayed neutrons such as the small reactivity signal assumption used to obtain the usual system transfer function (SCHULTZ, 1961). which is

Also it is evident from Fig. 1 that

0 I & 5 1 for 0 5~~

I1

for

p. 2 p.

0, 20

which is useful information in establishing limiting range values, particularly when numerical or perturbation methods are employed in solving equations (4). Furthermore, these ranges indicate in general that solutions for E are preferred for reactivity

: Scale forOSnmSl is linear // p,- IO.00 2.00

01

-10.00

FIG. I.-Asymptotic

I

-

A/X

IO’ c

IO *

values of q, and E, for various step reactivities

IO7

&&.250

Note

in dollars.

Id

m

IO'

’

_

218

G. M.

SANDQUIST

functions such that p(t) > 0 and q is preferred when p(t) < 0. However, we shall find that the range of parameter values associated with physical reactor systems will modify this observation somewhat. PERTURBATION

SOLUTIONS

OF THE

SYSTEM

First order perturbation solutions for equations (4) are found by letting 1;1= 1 + 67 and E = 1 + SE and assuming that the following terms of second degree in 67 and St are negligible for all time. That is

(8.1) and

The value range for the parameter ratio A/L for physical reactor systems is approximately (KEEPIN, 1965) 1 1 (typically R/I M 103 and expression (8.1) is more easily satisfied than expression (8.2) particularly as the perturbations 67 and SE become large. Note, however, the special case when 67 = p/(1 - p), then expression (8.1) implies that 1 < 11 - l/p1 independent of A/3, and expression (8.2) implies that 1 <<11 + n/(&l which are both satisfied if p is sufficiently small. In general, for monotonic reactivities such that p(t) > 0 for all time, both expressions are potentially satisfied while for monotonic reactivities such that p(t) < 0 for all time expression (8.1) is better satisfied. For nonmonotonic reactivities (particularly periodic functions with a zero mean) again both expressions are potentially satisfied. Since equation (4.2) for 7 will be the basis for most of the solutions developed in this study because of its broad potential application and simpler integration requirements, it is particularly important to observe the range capacity for which expression (8.1) is satisfied. In terms of the neutron and precursor density expression (8.1) becomes (for A/3, > 1) N

A h 1++11);

JI

t1. (8.3) I) For example for prompt critical operation (p = 1) and conservatively assuming that C does not increase during the rapid increase of N (i.e. C = l), then so long as N Q d(A/l) perturbation solutions of equation (4.2) are valid. If A/At w 104, then we require that N < 100. This is in sharp contrast to the usual perturbation assumptions made for equations (1) which require that the change in the neutron density be small compared to the equilibrium density. Note that for fast reactor systems A/A --t lo6 to lo* and expression (8.3) is particularly well satisfied. Assuming the validity of the approximations given by expressions (8) the C<

* Vertical line enclosures denote the absolute value of the enclosed expression. t For prompt critical operation and similar rapid reactivity changes, the effective il which should be used is essentially determined by the first one or two fastest decaying delayed neutron groups because of the rapid time behavior of the system.

General solutions of reactor kinetic equations for one group of delayed neutrons

219

perturbation solutions for equations (4) with initial conditions that 67(O) = @(O) = 0 (9.1)

Q(O = [ 1 + y Sdexp(yt-ASdpdr)dl]exp(nSbpdr-yr)-1, secr)=[I+rS:exp(yf+A~pdr)dr]exp(-A~~df-yl)

-1

are where y = A + 1. We see that the perturbation approximations expressions (8) are not equivalent. If they were, we would have since

(9.2)

assumed in

YE = (1 + SV)(l + 65) = 1 that

but equations (9.1) and (9.2) do not identically satisfy this equation, From equations (6) we have as the general perturbation solution for the neutron density when N(0) = 1 that N(t)=

(10.1)

(1 + 6q)exp [Ii’hdt],

(10.2)

where 6r] and SE are given by equations (9.1) and (9.2) respectively. Of course, equations (10.1) and (10.2) are not identical and give different solutions for a given reactivity function. For example, if p is a constant step input p,,, we find from equation (10.1) that N(t) = [ 1 + p0 c (1 - eYa’)] exp [pO $

(crt + e-at - l)]

(11.1)

and from equation (9.2) that N(t) =

I

1-

p. a +:,,

(1 - exp {--(a

+ 2wW~)]~A’aexp

[,Yt]

(11.2)

where a = A + A(1 - p,,) and y = A + 1. However, for small p0 and A > 1 equations (11.1) and (11.2) are essentially equivalent. For then both equations are approximately given by N(t)=

[I + p,i(l-e-at)]

exp [pOGt]

(11.3)

where the first term in brackets describes the initial prompt change and the second term, the exponential term, describes the response after prompt effects have disappeared. A generalized plot of equation (11.3) given in Fig. 2 reveals the general time behavior anticipated for the neutron density given a small step change in the reactivity. Furthermore equation (11.3) may be rearranged to give

General solutions of reactor kinetic equations for one group of delayed neutrons

221

and as t --+ cc we have that N( co) = exp [A/v(e”* - l)]. If aR < 1 and A > 1,, then equation (12.2) becomes N(t) = (1 + afl e-?“) exp [aL(l - e-‘?]

(12.3)

and N(co) = exp (an). Thus, if a given physical reactor system can be safely pulsed, the insertion of a known impulse of reactivity could provide an experimental measure of the effective delayed neutron decay constant to be used for reactivity functions imposed upon the system which change rapidly in time. The exact solution for the neutron density from equations (1) for an initial impulse reactivity input is N(f)

= 1 +

anll+ an2

e-Yt

Y

Y

and the required conditions for equivalence of this equation and equations (12.2) and (12.3) (viz. that aA < 1) are found to be consistent with the assumption contained in expression (8.1). 2. Solution for ramp reactivities If the reactivity input is a ramp (linear) function in time (AKCASU, 1958; WALLACH, by p = 2at, then we find from equation (9.1) that

1950) given

t q(t) = 1 +6q

=

exp (yt - m2t2) dt

1+ y

J0

[

1

exp (m2t2 - yt).

Performing the required integration we have that r(t) = [1 + 2/(n) pep2 {erf(p) - erf(p - mt)}] exp (m2t2 - yt)

(13.1)

where ~=A+13 erf(z) = *)

m = 2/(aA), 2 s

p = y/2m,

z o exp (-t2) dt

(error function).

The response of the neutron density for a ramp input of reactivity is then given by N(t) = rj exp [Ai{erf(--ip) - erf(imt - ip)} - B(1 erf (p) - I erf (p - mt)} - At] (13.2) where

A=

z

B = pL/m,

[2/(r) p

erf(p)

+ exp (-p2)],

i = 1/(-l),

I erf (z) = d\/(r) lo2 exp (t2) erf (t) dt = $1

.

3 5 ‘;;:I . ...

l)(n + 1).

It is important to note that the functions Ai[erf (-ip)

- erf (imt - ip)] and

BII erf (p) - I erf (p - mt)]

222

G. M. SANDQUIST

are real valued functions of time for all values of a. Also it can be shown that the function I erf (z) converges uniformly for all finite values of z. The Handbook of Mathematical Functions (1964) is a useful reference in deriving and evaluating these results. In the absence of delayed neutrons equation (13.2) reduces to the standard solution for ramp reactivities found from equation (1) [viz. N = exp (aAP)]. For sufficiently small values of time, the function I erf (z) may be adequately represented by a truncated Taylor series in time around p/m so that we obtain for the neutron density N(t) = [l + Z/“rpePa {erf (p) - erf (p - mt)}]

x exp g

(3 + 2p2)t5 - m$

t4 + i m2t3 + m2t2 - yt

[

1

3. Solution for sinusoidal reactivities If the reactivity input is a sinusoidal function (AKCASU, 1958; SMETS,1958) given by p = a sin cot, then from equation (9.1) 1+ y

q(t) =

’exp

(yt - q(l - cos wt)} dt] exp [q(l - cos

cot) -

p]

1 0

c

1 (15.1)

where q = an/w. To perform the required integration we express the following term in the integrand as exp [i(2m - n)wt] m 4” exp(qcosW= J ;cosnmt=$ i 4” 2~m,(n_m), . (16) n=Om=o Substituting equation (16) into equation (15.1) and performing the integration we have that q(t) = exp [q(l - cos ot) - p] + exp (-q

cos ot - yt)L(n, m)[exp {yt + i(2m - n)ot} - l]

(15.2)

where L(n, m) is a summing operator to be defined presently. Since the imaginary terms cancel from equation (15.2) then 7 can be expressed also as q(t) = exp [q(l - cos cot) + exp (-q

yt]

cos ut)M(n, m)[y cos ((2m - n)ot) -ye+ + w(2m - n) sin ((2m - n)ot)]

(15.3)

where the summing operators M(n, m) and L&z, m) are given by

Wn, m) = y

2 5

4”

n=O m=o 2”m!(n - m)![y2 + 02(2m - n)2] ’

L(n, m) = M(n, m)[y - ico(2m - n)]. Now substituting equations (15.2) for q into equation (6.1) and using the expansion given by equation (16) with q replaced by -4, the resulting expression may be integrated to give as the solution for the neutron density with a sinusoidal reactivity

N(t) = 7(t) expbM> + +2(f)+

+3(t)l

(17)

General solutions of reactor kinetic equations

for one group of delayed neutrons

223

where q is given by equation (15.3) and for the $‘s we have

+l(t>=

w, m>[r - ye-Yt cos ((2~

+

- n)wt} + w(2m - n)e-yt sin ((2m - nW}l

cos [(2m + 2j (2m - n)” -2mf2j-n-k

n - k)ot]}

17

$3(t) = ilR(n, m, k, j){ [y2 + 02(2j - k)(2m - n)][e-‘” cos ((2j - k)wt} - l] - uy(2j - k + n - 2m)edYt sin ((2j - k)ot}} where the summing operators P(n, m), Q( n, m, k, j) and R(n, m, k,j) are given by (- l)“q” P(n, m) = e* 2 f$ n=~ l,l=o2”m!(n - m)![y” -t (2m - n)2~2] ’ ( - I )kq”+k

Qh myk,.d = ~~~on~ok~oj~02 “+gm!(n - m)!j!(k - j)![y” + d(2m - n)“] 1 Rh

m, k,

j) = Qh

m, k, j)

[y2 + w2(2j - k)2] ’

We see in equation (17) for the neutron density that the presence of delayed neutrons excites all harmonics of the fundamental frequency of oscillation. For a system model with no delayed neutrons (i.e. A = 0), the neutron density is bounded for all time and is given by N = q where 7 is given by equation (15.3). It is also important to observe that for sinusoidal reactivities of any amplitude or frequency, the presence of delayed neutrons makes the system unstable. Furthermore, perturbation approximations of equations (4) do not eliminate this real system instability as do the standard perturbation approximations of equations (1). This unstable response can be demonstrated by noting that although q, & and 43 in equation (17) for N are bounded in time, the term $2 for the value of the sum when 2m + 2j = n + k becomes +,(t, n + k = 2m + 2j) = -At + itS(p, k, m, j) where S(p, k, m, j) = 2

2 “i’ 2

p=Ok=O

?n=O~=O

(- l)kq2”y2 ’ 22pm!(2p - k - m)!j!(k - j)![y2 + w2(2m - 2p + k)2] *

Now we observe that for p = 0, then S(0, k, m, j) = 1 and q$(t, n + k = 2m + 2j) = LtS(p 2 1, k, m, j) = s2t,

since the time coefficient S is always postive for (u2 > 0. Thus, we have that the neutron density behaves asymptotically as N + exp (c2t). So the neutron density with a sinusoidal reactivity and delayed neutrons exhibits exponential growth with a period given by 1/c2. If q <<1 (i.e. for sufficiently small reactivity amplitudes and oscillation periods), then &2is well approximated by the p = 1 term for S and c2 m

S(p = 1, k, m, j) = 4

where q = aR/co and y = A + A.

1 4 - + y2 + c.2 ^J2+ 402

(18)

G. M. SANDQLIIST

224

Figure 3 presents the time behavior of the slope of the maximum envelope bounding the oscillating neutron flux for the indicated sinusoidal reactivity input as determined by computer simulation of equations (1). Also plotted as dashed lines are the associated values of s2 as given by the approximation, equation (18). Note that asymptotic agreement is apparent for all reactivities considered. If .s2is sufficiently small so that the growth of the system response is slow compared to prompt neutron transients, then the response of the neutron density after prompt neutron transients have become negligible (i.e. for t > 3/y) is given by equations (17) and (15.3) with the term exp (-yt) set to zero. 4. Solution for generalized reactivities An analytical solution for the neutron density can be readily found from equations (9.1) and (9.2) for example for any choice of reactivity function such that the following expressions are analytically integrable.

s c

s

eYtf(O dt

dt,

f(t)

s

e-“ - W f(t)

= g(t),

Computer ----

Simulation

Asymptotic

&

dt

slope

given

by f

2

::‘I:: ,----p=o IO S,” 5 f

--

p=o-IO SI” IO f

---

4

s1n20 1 -----

p=o,20

\

-
05s1n5t

----I

I

icoo

I

Time,

FIG. 3.-Time

I 2000

I

I

3000

set

behavior of the slope of the maximum envelope bounding the neutron flux for various values of reactivity.

General solutions of reactor kinetic equations for one group of delayed neutrons

225

where f(t) = exp (--ASP dt). Usually the most fruitful procedure is to assume reasonable functional form for J(f) and determine the resulting reactivity function from the relation p(l) = - i i In f(t) As an example a useful reactivity model is found by letting f(t) = eMut- xebt where a, b and x are arbitrary parameters. It can then be shown that the reactivity function is given by 1 - exp {(x - 1)AApt) (19.1) ~(0 = PO + APX 1 - x exp ((x - 1)RApt) where a = A(p, + xAp) and b = -A(po + Ap). Note that p(0) = p. and p(co) = p. + Ap if (x - 1)AAp > 0. To permit a universal plot of equation (19.1) for the reactivity model for all parameter values, we define new variables such that e(s-l)T

1 -

8(x,

T) = y

= x

1 _

x

e(“-l)r

(19.2)

where T = AApt. Note that although t 2 0, the range of T is ---co I T I co depending upon the sign of Ap. A plot of equation (19.2) given in Fig. 4 for various

FIG. 4.-Generalized

reactivity function curves.

226

G. M. SANDQUIST

values of x displays the range of reactivity forms that can be simulated by the model. For example, a negative ramp insertion with termination is approximated in the region where T 2 0 and x < 0 (particularly for x < -4). Since p,, is arbitrary, this region also could be used to simulate an impulse input with variable die-away. The solution for 17from equations (19.1) and (9.1) is found to be KlemYt+ KZeeot + K,ebt 7(t) =

(20)

pzt - xebt

where

K,=l--+_t, Y--a

K,=L,

Y+b

KS=--$.

y=Ai-1.

y--a The solution for the neutron density for the generalized reactivity function given by equation (19.1) is then

IK, mzO$

(eDnt - 1);

for (a + b)t > -4 In x2

n

F(r) =

AK, 2 - ’ (e--Bet _ 1); 12=0Enx”+l

for (a + b)t < -4 In x2

and D, = n(a + b) + a - y, En = n(a + b) + b + y. As a specific example of the usage of this reactivity model suppose we wish to determine the system response for a negative ramp reactivity of slope (-$l.OO/sec) inserted for one second and then terminated. Essentially we want p = --t for 0 I t I 1 and p = - 1 for t > 1 with p smooth and differentiable at t = 1. Assuming that A = 5 set-l and using the information that P(0) = 0,

p(1) = -1,

and

7(t = 1) = 1,

we find that Ap = i and x = -5. Using these values for po, Ap and x we then find that a = -5 and b = -1. With the parameters a and b known, the neutron density is given by equation (21) if ;I is known. Of course, more ambitious reactivity models could be determined. For example, the selection of f(t) = e-=* - xebt + y where y is an additional parameter includes the model presented here as a special case when y = 0. The mathematical details for obtaining solutions are rather extensive, however. APPROXIMATE

NON-LINEAR

SOLUTIONS

OF THE

SYSTEM

Equation (4) has the general form dx z + f(p)x + x2 + K = 0.

(22)

where K is a constant and f[p(~)] is a linear function of the reactivity p. An obvious approximation to consider is that of neglecting the constant K if possible and then

General solutions of reactor kinetic equations for one group of delayed neutrons

227

solving the resulting Bernoulli equation to obtain as the general solution

exp(-SfW

X(T) = x0

-I- .f exp (-JfdT)

d7

where x0 is the arbitrary integration constant. The assumption term I is negligible in equation (4.1) implies that

that the constant

(23.1) and in equation (4.2) the constant term A is negligible if A ;I

A <

7

(P2- 1) + 1 -

(23.2)

1;1 11.

The terms p1 and pz are to be distinguished from the formal definition of reactivity since if the expressions are satisfied p1 = pZ = 0 does not correspond to equilibrium operation as required in equation (1). In general since A > 1 for reactor systems we find that expression (23.1) is usually more easily satisfied than expression (23.2.) Furthermore neglecting iz in equation (4.1) implies for critical operation that p1 = -1/A which is small and generally an acceptable approximation. However, neglecting A in equation (4.2) implies that at critical operation ps = 1 which corresponds to a prompt critical system. Of course the functions p1 and pZ should be redefined so as to render the resulting approximate system equations with equilibrium operation occurring at zero reactivity. Thus the system equations utilizing the approximations given by expression (23) are

$ 5 + 4p $7 where p = p1 + J/A q(O) = 5(O) = 1,

CzI

- 1)E+ AiT2= 0,

(24.1)

- (A + Ap)r + Jr2 = 0

(24.2)

p2 - 1. Solutions for 17 and E are the following where

t(t) = epcc) ( 1 + A /oe”‘Ldr)W1

(25.1)

and

where F(t) = At -

ts

A o P dt,

(25.2)

G(t) = /It + A

s0

‘p dt.

The resulting solutions for the neutron density recalling equation (10) are N(t) = (1 + AJoteP dt) exp p /o’e-P (1 - A leF

dt) dt - F - At]

(26.1)

228

G. M. SANDQUIST

and N(t)=

(1 f

%~ote”dt)-lexp

[flSde@

(1 +

i,JoteGdr)-‘dt

+

G-it].

(26.2)

For the case where the reactivity input is a constant step input p0 then F(t) =

A(1

-

p,Jt =

F,t

and

G(t) = (3, + Ap,Jt

and the solutions for the neutron density are given by APO

N(t) = 1 + F

(1 -

e-Fat) exp

(W

+ e-Fat- I)]

(27.1)

0

associated with system equation (26.1) and (27.2)

N(t) = exp (Apot)

which is the solution associated with equation (26.2). Equation (27.2) does not give even the correct asymptotic period for the system response and is an unacceptable approximate solution. Thus expression (24.2) is not valid and equations (25.2) and (26.2) are unacceptable approximations of the system behavior. However equations (25.1) and (26.1) are acceptable approximations for the system behavior and are useful companions for the results found from perturbation approximation. It is interesting to observe that for small p. and A > 3, then F, + a and equation (27.1) reduces to equation (11.1) which provides corroboration for both approximations. EXACT

SOLUTIONS

FOR

THE

SYSTEM

It is possible to determine exact solutions for equations (4) and even equation (3) when the reactivity is a known analytic function expressible as a power series in time. If we consider the general form of equation (4) as given by equation (22) with f(p) expressed as a power series expansion of the independent variable, since f(p) is a linear function of p(t) which is analytic, then dx z + x2 + x # b,r” + K = 0. We assume an analytic solution of the form X(T) =

2 umTTn.

(29)

0

Substituting the series expansion for x into equation (28) we find that

$ ma,r-1

+

[za,~“]~+ f, $

b,a,P+

+K =

0.

The second term can be expressed as

2= and so we have that akaf_krj t 2 f b,a,r”‘+” +K=O. 0 k=o

By

0

0

redefining the dummy variables of summation as m = p +

1,

j = p and m + n = p

General solutions of reactor kinetic equations for one group of delayed neutrons

229

in the first three terms respectively and grouping on equal powers of 7 we have that ao2 + a, + K + $

Since we

[

(p + l)~,,

+ &La,,-,

+ b,-,)

19=0.

want this expression to be satisfied for all values of T we require that the coefficients of each different power of r vanish so we require that a, = -K

- ao2

and for p 2 1 that

The constant a,, is the arbitrary solution (or integration) constant which is determined by the choice of an initial condition for x. All the other constant coefficients a, (for n 2 1) are defined in terms of a, and other known quantities and thus we have found a solution of equation (28) in the analytic form given by equation (29). It can be shown that the series expansion give by equation (29) is the complete solution (HILLE, 1969) and converges uniformly (TITCHMARSH, 1939) within the common radius of convergence of equation (29) and the analytic reactivity function p(t). If we express the general solution of equation (28) in terms of 17as given by equation (4.2) we have that (30) where N, = N(0) and C,, = C(0) are the initial conditions for the neutron density and delayed neutron precursor density respectively. Note that ~(0) = N(O)/C(O). The neutron density for arbitrary initial conditions and any analytic function of reactivity is given by the following equation. N(t) = No

[

1+

z

A,(h) n] exp [@-

1) il+$$&g]

(31)

where

A Co No

A’=hN,-C,’

A

A, = -

2 0 (A, + B,),

2A, + B, + B, + ‘$I&$,_,

9+1 =

-

1 p+lCo

3

1

+ B&

1

for p 2 2. The arbitrary analytic reactivity function is determined by specifying the B, coeffcients for n 2 0 in the following power series expansion. (32) By a similar procedure it is possible to determine a general solution for the neutron density for M groups of delayed neutrons when the reactivity is a known analytic function of time. However a reversion of the series expansion for ti is required to perform the integration step indicated by equation (6.2).

230

G. M. SANDQUIST

CLOSURE

By a transformation involving the ratio of the dependent variables we have shown that the reactor kinetic equations for one group of delayed neutrons can be reduced to a single first order equation. Although this single system equation is non-linear exact series solutions have been found for analytic reactivity functions. Of more practical value, however, perturbation solutions of this equation have been shown to be accurate and useful for large changes of the system variables and therefore can provide practical analytic solutions for the kinetic equations. Acknowledgments-The author wishes to express appreciation to Professor E. P. Gyftopoulos for valuable discussion and criticism of this study, Forster Betts, and to the National Science Foundation for providing support through its Science Faculty Fellowship Program for postdoctoral study at the Massachusetts Institute of Technology. REFERENCES ABRAMOWITZW. and ST~GUN I. A. Ed. (1964) Handbook of Mathematical Functions, Chap. 7, National Bureau of Standards. AKCA~UA. Z. (1958) Nucl. Sci. Engng 3,456. HANSENG. E. and MAIER C. (1960) Nucl. Sci. Engng 8, 532. HENRY A. F. (1958) Nucl. Sci. Engng 3, 52. HILLE E. (1969) Lectures on Ordinary D@rential Equations, pp. 48-56. Addison-Wesley, Massachussetts. INCB E. L. (1956) Ordinary D@rential Equations, pp. 22,23,295. Dover, New York. KEEPIN G. R. (1965) Physics of Nuclear Kinetics, Addison-Wesley, Massachussetts. KEEPIN G. R. and Cox C. W. (1960) Nucl. Sci. Engng 8,670. KIRCHENMAYER A. (1960) Nucl. Sci. Engng 8,720. SCHULTZ M. C. (1961) Control of Nuclear Reactors and Pobver Plants, 2nd edn. McGraw-Hill, New York. SKINNERR. E. and COHENE. R. (1959) Nucl. Sci. Engng 5,291. SMETSH. B. (1958), Low and high power nuclear reactor kinetics, Sc.D. Thesis, MIT. SMETSH. B. (1958), Proc. U.N. Intern. Co@ Peaceful Uses Atom. Energy, 2nd, Geneva, 11,237. TITCHMARSHE. C. (1939), The Theory of Functions, 2nd edn., Chap. I and VII, Oxford University Press, Oxford. WALLACHS. (1950), Solutions of the pile kinetic equations when the reactivity is a linear function of time. WAPD-13, Westinghouse Electric Corp.