SELECTED TOPICS in REACTOR PHYSICS II.

By Professor Zoltán SZATMÁRY

From Budapest University of Technology and Economics

 

Table of content

The adjoint function and its applications. 1

1. The adjoint equation. 2

Neutron field and detector field. 2

Neutron importance. 5

The adjoint transport equation. 7

The adjoint diffusion equation. 8

2. Time dependence of the adjoint function. 10

Source importance. 12

Static eigenfunction. 13

3. Variational principles. 14

The Lagrangian of transport theory. 14

Stationary formulae. 15

4. Perturbation theory. 17

Perturbation of the reactivity. 17

The characteristic of a control rod. 19

5. The point kinetic equation and the effective delayed neutron fraction. 20


 

The adjoint function and its applications

 

            Neither the transport equation nor the diffusion equation is self-adjoint – in contrast, for example, to the Schrödinger equation. It would not be difficult to formulate the equations adjoint to the transport or diffusion equations (at least formally). However, we are interested in a more fundamental question: do the solutions of the adjoint equations have some physical meaning? We shall see that the answer is: yes! As a matter of fact, a proper understanding of the transport equation is not possible without understanding the adjoint equation. The physical meaning can be best understood in the framework of the simplest model e.g. the one-group diffusion equation. This is the only approximation of the transport equation which self-adjoint.

            Consider a bare critical reactor in which the neutron flux is constant in time. This flux is proportional to F1(r) that is the eigenfunction of the Helmholtz equation:

(1)

belonging to the smallest eigenvalue . Assume that a neutron appears at the point r0 in the moment t = 0. After this, the flux will become time dependent. The initial condition for it can be written in the form

(2)

where c is some constant. For later moments t > 0, the time dependence of the flux is determined by the time-dependent one-group diffusion equation:

.

(3)

We write its solution in the form of the series

(4)

where the functions Fn(r) are the eigenfunctions of Equation 1 (n = 1, 2, ...). The coefficients fn are the expansion coefficients of the initial distribution given by Equation 2:

.

 

For a critical reactor, w1 = 0 and wn < 0 for n > 1. Thus, for sufficiently large times t, the flux goes over in a constant value again:

.

(5)

We have found that the appearance of a single neutron increases the steady-state flux by a value proportional to F1(r0). We can put it also in an other way: for the chain reaction, it is the worth of a neutron injected at point r0. This is the idea of neutron importance. Equation 5 shows that, in one-group diffusion theory, the importance is equal to the flux.

            In the next section, we shall derive the equation satisfied by the neutron importance. We shall find that this equation is the adjoint of the transport equation. In the present section, our starting point was the one-group diffusion equation (see Equation 3). All operators in the right hand side are self-adjoint. Thus, it was not by mere chance that the neutron importance was found to be equal to the neutron flux.

            There are problems in transport theory which, principally, could be solved within the framework of transport theory, but, practically, they cannot be solved without the adjoint equation. An important set of such problems is the treatment of small perturbations: we change some characteristic of the reactor by some small amount (for example, we move a control rod a little), and we want to determine the corresponding changes of the reactivity and the neutron flux. Since we want to predict the response of the reactor to small changes, such problems can not be treated by the usual numerical methods. Let us consider the example of the reactivity change. Denote by  and  the values of the multiplication factor corresponding to the initial and final positions of the control rod, respectively. Then the reactivity changes by

.

 

The usual order of magnitude of Dr is 10–5 and we want to get it as the difference of two numbers which are both close to 1. Any numerical solution of the transport equation is only an approximation. We would get a rather poor accuracy by calculating the small difference of two large quantities. In such cases, it is always better to look for such methods which, although approximate as well, give the difference directly. This is the basic idea of perturbation theory. It was first elaborated for the Schrödinger equation. Later, it was developed for the transport equation as well. However, there is a difference with respect to the Schrödinger equation: it is necessary to solve both the transport equation and its adjoint.

            Although the flux F is calculated for most practical applications, the solutions of the adjoint equation can be interpreted consequently only if we start out from the transport equation formulated in terms of the neutron density n. One-group theory is an exception: the neutron speed v is not a variable of the equations, thus, the flux and density are proportional to each other.

 

1. The adjoint equation

Neutron field and detector field

            Since – as far as possible – we do not wish to repeat our considerations for all approximations of the transport equation, we put our starting equations in an operator form:

(1.1)

where

.

(1.2a)

The operators are defined for the transport equation as follows:

 

(1.2b)

is the destruction operator describing the disappearance of the neutrons while

(1.2c)

is the production operator describing the neutrons produced by fission. In Equation 1.2b, the kernel S(r,E¢®E,) is the sum of elastic and inelastic scattering kernels. These operators have the following form in diffusion theory:

 

(1.3a)

and

(1.3b)

where – as before –  the kernel S(r,E¢®E) is the sum of elastic and inelastic scattering kernels. The solution of Equation 1.1 has to satisfy the boundary conditions belonging to the model chosen. The initial condition is written in the form

(1.4)

where the function ni has to satisfy the same boundary conditions as the neutron density n itself.

            There is no instrument which would allow the measurement either of the neutron density or of the neutron flux. All what can be done consists in putting some neutron detector in the reactor and in counting the nuclear reactions induced in it by the neutrons. The result of the counting will be called meter reading in the following. If the cross section of the detector is Sd(r,E) for neutrons of energy E at point r of the reactor, the number of nuclear reactions in unit of time, in volume element dV at r is given by the integral

.

 

The total meter reading at time t is then

 

(1.5a)

where V is the volume of the reactor. In diffusion theory, the same thing can be written as

.

(1.5b)

We cite two examples for explaining what we mean by “detector”. If it is a BF3 counter, Sd is proportional to the absorption cross section of 10B, and the integrals for r in Equations 1.5 are extended only to the detector’s volume since Sd vanishes outside the detector’s volume. Each neutron absorbed in the detector results in an electric pulse, thus, the meter reading is obtained by counting the pulses. Another example is the reactor power: Sd is then proportional to the fission cross section of the reactor fuel, thus, the integrals in Equations 1.5 are extended to the whole volume of the reactor core.

            It follows from these considerations that the observable quantities cannot be described by the neutron field alone but we need another field, too, that is the detector field. Any function belonging the detector field will be identified by a superscript “+” such as y+. The variables of the detector functions are the same as those of the model actually used. This means that we have y+(r,E,,t) and y+(r,E,t) in transport theory and diffusion theory, respectively. The physical meaning of the former is the contribution of a neutron characterized by (r,E,) at time t to the meter reading measured at time t.[1] The meaning of y+(r,E,t) is analogous. According to this definition, we have in Equations 1.5

.

 

The integrals given in Equations 1.5 are inner products of the detector function and the neutron density. We use the same notation for them as in functional analysis:

.

(1.6)

Variables r, E, and  are omitted in formulae like this but the time variable t is given when, in the actual context, it is necessary to specify for which moment of time the inner product is defined.

            The detector functions y+ may be arbitrary in principle but we suppose them to satisfy some boundary conditions which follow from their physical meaning. Since a neutron flying outwards at the outer boundary of the reactor gives no contribution to any meter reading, we require:

(1.7a)

where Ns is the normal vector of the boundary pointing outwards (Figure 1.1). Note that this is just complementary to the condition set for the neutron density:

.

(1.7b)

In diffusion theory, for the same reason, we require that y+(r,E,t) vanished at the extrapolated reactor boundary. Note that the same boundary condition is set for the neutron density.

Figure 1.1. Illustration of the boundary condition set for the detector functions in transport theory

Neutron importance

            Assume that a measurement will be performed at t = tf by using a detector with detector function .[2] The meter reading will be by definition

.

(1.8)

In order to calculate n(tf), we need to know the initial neutron density ni for some time ti < tf. Equation 1.1 determines how the neutron density develops between the two moments of time. Now we formulate two definitions:

·     For t < t¢, consider a neutron in both of these moments of time. The latter is called a descendant of the former if its birth can be traced back to the former through translation, scattering and fission.

·     We call n+(r,E,,t) neutron importance if it gives the total contribution of all descendants of a neutron characterized by (r,E,) at moment t to the meter reading measured in the moment tf (t < tf).

It follows from the definition that the importance has to satisfy boundary condition 1.7 as well as the following “final condition”:

.

(1.9)

The neutron importance n+(r,E,t) of diffusion theory can be defined in an analogous way.

            The neutron importance defined in the introduction is a special case of this general definition. As the relation between the two definitions is not trivial, we give the following explanation. Firstly, we used one-group theory there, thus, variables E and  are omitted. Secondly, we considered the solution of Equation 4 in the limit of tf ® ¥. Thirdly, we considered the neutron flux at point r. Equation 1.8 gives it for

.

 

            In order to get an equation for the neutron importance, we apply the same reasoning which leads to the differential form of the transport equation. Our equation will express the following obvious statement: the importance of a neutron is equal to the sum of the importances of its descendants. We can call it the principle of importance conservation.

            Let a neutron be at point r at time t and having energy E, flying in direction . After a time lap of dt, we find it at point (r + vdt), unless it induces some nuclear reaction. The probability of the former eventuality is (1 – vSt(r,E)dt). According to the definition of the kernel  (see Equation 1.2), the probability that, following an (elastic or inelastic) scattering, the energy of the neutron changed to a value in the interval (E¢E¢+dE¢), its flight was deflected to a direction in the solid angle d around , is given by

 

Finally, the total number of fission neutrons of energy in the interval (E¢,E¢+dE¢), of flight direction in the solid angle d around , is

.

 

It follows from the conservation principle of the importance that

 

 

.

 

We get in the limit of dt ® 0:

 

.

 

 

If the total time derivative is expressed by partial derivatives, we get an equation similar to Equation 1.1:

(1.10)

where

 

(1.11a)

and

.

(1.11b)

We have got an equation for n+, in which time t passes “backwards”: if we want to compute n+, we start from the final condition specified by Equation 1.9 for the moment tf, and we solve Equation 1.10 backwards in time.

The adjoint transport equation

            We will show that the operators , , and  defined by Equations 1.11 are adjoint to operators , , and  of Equations 1.2. This means that we have to prove the validity of the following equations:

,

 

 

where n and y+ are arbitrary functions of the neutron and detector fields, respectively, which satisfy the respective boundary conditions. If these equations hold, it follows from this that

.

 

Let us look at the operators  and  terms by term:

·     The kernels of the integral operators are the same except that variables E and E¢ exchanged: we have  in  (see Equation 1.2b) while the corresponding term of  is . Since  and  appear symmetrically, this means that the integral operators are adjoint to each other. As an exercise, we advise to check the same thing for operators  and .

·     The term containing St(r,E) is trivially self-adjoint.

·     The leakage terms are identical except the sign. We prove our statement by taking into account the following identity:

grad(gh) = g×grad h + h×grad g

 

which holds for any differentiable functions g(r) and h(r). We may now write

 

.

 

We apply the Gauss theorem in the first term of the right hand side:

.

 

This surface integral vanishes for every : if  points outwards (see Figure 1.1), we have y+ = 0 by Equation 1.7a, if  points inwards, we have n = 0 by Equation 1.7b. This means that

.

 

This completes the proof of that the adjoint operator of grad is –grad.

            We can now state the following conclusion: the neutron importance satisfies the adjoint transport equation. For this reason, the quantity n+(r,E,,t) is sometimes called the adjoint function.

            The formal relation between the transport equation and its adjoint leads certain authors to call the function n+ adjoint flux. We do not subscribe to that practice. It follows from our derivation and the physical meaning of n+ that there are no particles whose transport would be described by the adjoint transport equation. In addition to this, we have a simple formal objection, too. We have seen that the functions of the detector field are the form vS i.e. their dimension is 1/s. At the same time, the dimension of the flux is 1/(cm2×s). This formal reason alone is sufficient to ban the expression “adjoint flux”.

The adjoint diffusion equation

            Their are three ways to derive the adjoint diffusion equation: firstly, we can give some physical meaning to the adjoint function n+(r,E,t) and apply the same derivation as above; secondly, we can start from the adjoint transport equation and derive the adjoint diffusion equation just as the diffusion equation can be derived from the transport equation; thirdly, we can formally find the operators adjoint to those defined in Equations 1.3. We choose the third way which is the simplest.

            We shall show that operators adjoint to Equations 1.3 are the following:

 

(1.13a)

and

.

(1.13b)

The statement is trivial except for the leakage term which is the first term of . The proof will be based on the identity

 

which holds for any (differentiable) function h(r) and vector g(r). Let n and y+ be arbitrary functions of the neutron and detector fields, respectively, which satisfy the diffusion boundary conditions. Then we may write:

 

.

 

 

We have applied the Gauss theorem for the first term of the right hand side. The surface integral is defined for the outer (extrapolated) surface of the reactor where y+ vanishes. The value of the surface integral is thus 0. In the second term of the right hand side, we apply again the above identity and the Gauss theorem:

 

 

.

 

n satisfies the diffusion boundary condition, thus, the first term on the right hand side vanishes. We conclude now that

.

 

This means that the operator div[vDgrad] is self-adjoint as stated in Equation 1.13a.

            As a conclusion of the present section, we give some hints concerning the derivation of the adjoint diffusion equation starting from the adjoint transport equation. When the diffusion equation is derived from the transport equation, the first step is the series expansion of the neutron density in terms of the spherical functions:

 

Here, the diffusion-theoretical density  is the integral of the transport-theoretical density  with respect to . Starting from this series development, we get Fick’s law:

.

 

In case of the adjoint function, we could formally generalize this as follows:

 

If we put the series development in this form,  would be the integral of the transport-theoretical adjoint function  with respect to . Unfortunately, this is not acceptable for two reasons. Firstly, in contrast to the neutron density , the neutron importance  is not a density with respect to : it is the importance of a neutron flying in some selected direction . Thus, we cannot assign any physical meaning to the integral of . Secondly, the transport-theoretical final condition 1.9 can be transformed into a diffusion-theoretical final condition if  is the average of the transport-theoretical adjoint function  with respect to :

.

(1.15)

In order to show this, assume that the detector function appearing in the final condition is isotropic. (This is the usual case.) Then Equation 1.9 has the form

.

(1.16a)

Its diffusion-theoretical equivalent should trivially read as

.

(1.16b)

Equations 1.16a and 1.16b can hold simultaneously if we use the definition 1.15.

            These considerations lead to the series development

 

If we apply it in Equation 1.10 and 1.11, we derive the adjoint Fick’s law:

(1.17)

and the adjoint diffusion equation corresponding Equations 1.13. When comparing Equation 1.17 with the Fick’s law related to the neutron density, we find a remarkable difference: the negative sign is absent in Equation 1.17. This can be interpreted as follows: vector J+ of the adjoint current points in the direction of the increase of neutron importance.

2. Time dependence of the adjoint function

            Two kinds of eigenvalue equations can be derived from the time-dependent equation 1.1: kinetic (or dynamic) and the static eigenvalue equations. Both are obtained from Equation 1.1 by setting S = 0. This corresponds to a reactor left alone without external neutron source.

            The kinetic eigenvalue equation describes a reactor which develops in a persistent mode i.e. in which the time dependence of the neutron density is separable from the dependence of the other variables:

.

 

By substituting this in Equation 1.1, we can easily show that the time dependence can be only an exponential function:

 

with some suitable w satisfying the kinetic eigenvalue equation:

.

(2.1a)

The static eigenvalue equation is the well known

(2.2a)

equation which is solved by many computer codes.

            It is simple to write the analogous adjoint equations:

(2.1b)

and

.

(2.2b)

The negative sign in the kinetic eigenvalue equation is due to the negative sign of the time derivative in Equation 1.10. If necessary, we apply a subscript “s” in order to identify quantities belonging to the static eigenvalue equation.

            In all what follows, we try to keep general notations by using operators because this allows to make statements which are independent of the model actually used. If it is necessary to specify the other variables in addition to the time variable t, we shall use vector x for that purpose. Its meaning depends on the actual model. For example:

·      transport theory: x = (r,E,);

·      one-group transport theory: x = (r,);

·      diffusion theory: x = (r,E);

·      one-group diffusion theory: x = (r);

·      multigroup diffusion theory: x = (rg) where g is the group subscript,

and so on. When calculating inner products, we have to integrate for all components of vector x. For arbitrary functions y(x) and y+(x), we calculate the integral

.

 

For multigroup models, the integral with respect to energy E is substituted by a sum with respect to the group subscript g. For example, we have in multigroup-diffusion theory

.

 

            The eigenfunctions belonging to different eigenvalues are bi-orthogonal. In order to show this, we write the kinetic eigenvalue equations:

                 and       .

 

We take the inner products of these equations with the eigenfunction of the other equations, and we subtract the resulting equations:

.

 

If , the inner product  vanishes. Under quite general conditions, it is possible to choose the subscripts such that, for j = l, the equality  holds. We can show the bi-orthogonality for the static eigenfunctions too:

      for      .

 

Note that the orthogonality holds only when the production operator is included in the inner product.

Source importance

            Consider a subcritical reactor, which, as a special case, can be also a non-multiplying medium. For such a system, all kinetic eigenvalues are negative: wl < 0, for l = 1, 2, ... As a consequence of this, the neutron density tends to zero for t ® ¥ independently of the initial density. A steady state neutron density can prevail only in the presence of an external source:

(2.3)

where S(x) is a steady state external neutron source.

            The neutron importance has an analogous time behavior. Let us develop the importance in a series of the eigenfunctions:

         for        t £ tf.

 

By virtue of the final condition (see Equation 1.9), we may write

.

 

By the orthogonality proven above, we can obtain the coefficients of this expansion:

 

leading to the result

.

(2.4)

Since all eigenvalues are negative, the limit of this importance is zero for t ® –¥ independently of the final condition. This behavior allows to get an equation for the neutron importance similar to Equation 2.3.

            Assume that the external neutron source injects neutrons into the system since t = –¥. We ask now the following: what will be the total contribution of the descendants of these neutron to the meter reading at time tf. The neutron importance goes to zero for t ® –¥, thus, the contribution of the descendants of the neutrons injected sufficiently early with respect to tf will be negligible. Let the time interval (tt+dt) be close to the moment tf. The contribution of the descendants of the neutrons injected during this interval is by definition

.

 

If this is integrated for all moments, the result is the total contribution:

(2.5)

where

.

(2.6)

It follows from Equation 2.4 that this time integral is convergent. The quantity defined by Equation 2.7 is called source importance.

            Equation 1.10 allows the derivation of an equation for the source importance. If we integrate Equation 1.10 for t in the interval (–¥tf], we get

.

 

By taking into account Equations 2.6 and 1.9, this may be rewritten as

.

(2.8)

This equation is similar to Equation 2.3. We have thus found that there is an adjoint equation with source in which the “source” is the detector function.

            We give an example of applying our last result. Assume that a spatially extended neutron source injects neutrons into the system, and we are interested in the radiation dose in some dose point outside the system. In terms of our former notations, S(x¹ 0 in a large volume while  = 0 if x is other than the dose point. Such problems of biological shielding are often treated by Monte Carlo programs. Unfortunately, the calculation will not converge since the probability of that a neutron hits the dose point is zero if the neutrons are started from the neutron source S(x). The problem is lessened a little by not considering the detector a geometrical point but extending it to a small volume around the dose point. Then the probability of hitting the detector is no more zero but it is still small. The real solution is the adjoint Monte Carlo method which allows to solve Equation 2.8 by Monte Carlo simulation. Although we emphasized in Section 1 that there are no real particles the flux of which could be the adjoint function, we can consider the source importance, at least formally, the flux of some imaginary particles whose transport is described by the operator . Now we start the “adjoint particles” from  as the “adjoint source”, then the pay-off is computed according to Equation 2.5 when the “diffusing adjoint particles” cross the real neutron source. Since the real neutron source is extended in space, the adjoint particles cross it with a finite probability. Of course, even the adjoint Monte Carlo method is helpless if the neutron source S(x) is also concentrated in a point.

Static eigenfunction

            In practical applications, it is the static eigenfunction that is computed most frequently. This means the solution of the static eigenvalue problem defined by Equations 2.2a and 2.2b. The static eigenvalues constitute a monotonously decreasing sequence, the first one being the effective multiplication factor (keff). The static eigenvalue equations corresponding to it are

(2.9a)

and

.

(2.9b)

The subscript “s” means “static” equations. There are several computer programs that solve Equations 2.9 in various models (i.e. finite difference diffusion, finite element diffusion). Consequently, we may assume in the following that ns and  are known in every case considered.

3. Variational principles

            In addition to what was discussed above, the adjoint function plays a further role, too: it allows to formulate the transport equation (and all of its approximate versions) in form of variational principles. We apply variational principles frequently in theoretical physics. Their role is twofold: for one hand, they contain the basic equations in an concise form, and, for the other hand, they allow to compute certain quantities accurately (up to first order) by means of approximate solutions of the basic equations. Since the transport equation is not self-adjoint, the variational principles necessarily contain booth the neutron density and the adjoint function.

The Lagrangian of transport theory

            In transport theory, the Lagrangian is the meter reading defined by Equation 1.8. Let n+ be a (for the time being arbitrary) detector function which satisfies the boundary condition formulated in Equations 1.7. Rewrite R(tf) as follows:

 

where ti is the moment to which initial condition 1.4 applies. R(t) for ti < t < tf is given by the formula

.

 

If Equation 1.1 is also taken into account, the Lagrangian becomes

 

.

(3.1)

According to the variational principle, we ask the question: how to choose n+ in order to make R(tf) stationary against small variations of the neutron density? In order to find the answer, we change n by dn. The neutron density is given for t = ti (see Equation 1.4), consequently, dn(ti) = 0, leading to

.

(3.2)

This vanishes for all dn if

.

 

In other words: the Lagrangian is stationary if n+ satisfies Equation 1.10 i.e. the adjoint equation.

            We have found that the variational principle based on the Lagrangian 3.1 contains the adjoint equation. We can simply see that this principle contains the transport equation 1.1, too. We integrate in Equation 3.1 by parts:

.

(3.3)

Now we take the variation of this with respect to n+:

.

(3.4)

This is stationary for all dn+ if n satisfies Equation 1.1. Furthermore,

 

needs to hold as well, i.e. n+ has to satisfy the final condition 1.9.

            As a conclusion, we can state the following: the Lagrangian 3.1 is stationary with respect to both the neutron density and the adjoint function, and both equations can be obtained from the variational principle based on it.

Stationary formulae

            The most important application of the variational principle formulated above is perturbation theory which allows to compute small reactivity changes with high accuracy. Perturbation theory can be applied to the computation of small changes in the neutron flux, too, but we will not go into the details of this. Before presenting perturbation theory, we present a simple example illustrating the power of the method.

            Equation 1 (i.e. the Helmholtz equation) is self-adjoint. Consequently, we can show that the formula

(3.5)

is stationary against the variations of F . Assume that some function F(r) satisfies Equation 1. Then, its small variations dF(r) do not change B2 in first order:

 

 

.

 

Let us take the example of a sphere with radius R = 1. The Laplace operator has the following form in spherical geometry:

.

 

Equation 3.5 has the following form in this case:

.

(3.6)

            The buckling is known for this case: B2 = p2 = 9,87. In order to illustrate the method, we take some trial functions satisfying the boundary condition F(1) = 0 and certain trivial requirements. One of the simplest functions which are differentiable twice is

.

(3.7a)

Formula 3.6 yields the buckling  for it. Its error is 40%. This trial function has the drawback that it has no extreme for = 0 its first derivative being equal to –2. We modify the trial function:

.

(3.7b)

Its derivative vanishes for = 0. Formula 3.6 leads to  which is a much better approximation, the error being only 6.5%.

            We could imagine further trial functions. Now the question arises: which one is the best. We know from functional analysis that formulae like 3.5 have a minimum property: from among the trial functions invented, that one leads closest to the exact solution for which the buckling is the smallest. This observation allows to refine the trial functions further. Let us consider for example a linear combination of the function 3.7b and of its square:

 

where c is not specified for the time being. Now we proceed as follows: substitute F3 in Equation 3.6, then find the minimum of the buckling in terms of c. We get two values minimizing B2:

c1 = 1.098                    and                   c2 = –1.41.

 

The corresponding bucklings are

for            and                   for .

 

The error of the former value is only 0.1%. The other buckling belongs to the second harmonic which is an approximation for 4p2 = 39.5. Such methods are known as Rayleigh-Ritz methods in mathematics.

4. Perturbation theory

            Perturbation theory is applied when we have to compute the effects of small changes directly. Most frequently, we are looking for the change in reactivity. There are cases when we are interested in the perturbation of the neutron flux but we discuss only reactivity perturbations in the following. We call perturbation every small change in the operators describing the reactor. We give two examples for illustrating what can be considered small:

(1)   When the perturbation consists in moving a control rod a little, the cross sections change only in a small volume of the reactor. The changes in the cross sections can be locally large.

(2)   When the concentration of boric acid dissolved in the moderator changes a little, the absorption cross section is perturbed. Here the perturbation concerns the whole volume of the reactor but the change is small in every point.

Perturbation of the reactivity

            The reactivity is obtained by solving the static eigenvalue equation 2.9a. The first idea is to express the multiplication factor from the equation itself:

 

leading to the reactivity:

 

where y+(x) is an arbitrary weight function (belonging to the detector field). The perturbation of the reactor means that the operators change in Equation 2.9a:

                        changes to                    ,

 

                         changes to                    .

 

As a result of this, both the reactivity and the static neutron density change:

r                          changes to                    r + dr,

 

ns(x)                     changes to                    ns(x) + dns(x).

 

The reactivity perturbation is then obtained as the difference

.

(4.1)

            This formula yields the required results but it is not satisfactory for two reasons: firstly, its application requires that the density perturbation dns(x) is computed, secondly, there is now warranty that the effects of the approximations applied in solving Equation 2.9a were minimal. Both difficulties can be eliminated if the weight function y+ is chosen adequately. The variational principles discussed in the previous section suggest that the best choice is the static adjoint eigenfunction obtained as a solution of Equation 2.9b:

                    or                     .

(4.2)

We show that this is stationary against small variations of ns(x):

 

.

 

In the last step, we took into account Equation 4.2 and the fact that the adjoint function satisfies Equation 2.9b. It follows from this that this choice of the weight function eliminates – at least up to first order – the eventual approximations applied in computing the neutron density. Unfortunately, also the adjoint function can be computed only approximately. Thus, we could think that our method has the drawback that we have to know the adjoint function accurately. Things are not as bad as that: we could show by a similar derivation that Equation 4.2 is stationary against small variations of the adjoint function, too. We may thus draw the conclusion that the error of Equation 4.2 is of the order of

 

i.e. the error is only of second order.

            It follows from our considerations that the effects of small perturbations of the neutron density in Equation 4.1 are negligible, and it is sufficient to take into account the perturbations of the operators:

.

(4.3a)

This leads to the usual perturbation formula for the reactivity:

(4.3b)

where

.

(4.4)

            Equations 4.3 are the fundamental formulae of perturbation theory. They demonstrate how we could realize our goal formulated in the introduction: we can determine the reactivity perturbation by solving Equations 2.9 for the unperturbed state of the reactor. In terms of these solutions, Equations 4.3 yield the reactivity change related to an arbitrary perturbation. It is an important aspect of the method that the eigenvalue equations 2.9 need not be solved for the perturbed reactor state at all.

The characteristic of a control rod

            The characteristic of a control rod is a well-known application of perturbation theory. Assume that the rod moves in an empty tube parallel to the z-axis, and its horizontal coordinates are (x0y0). Let the reactor be critical when the rod is its upper position. Since the rod moves in an empty tube, the cross section are (practically) zero in points inside the tube when the rod is its upper position. We solve the eigenvalue problem 2.9 for this reactor state. We choose the one-group diffusion model according to which the solution is

(4.5)

where H is the reactor height and A is a coefficient depending on the horizontal rod position.

            When the rod is inserted to some depth, we can say (approximately) that only the absorption cross section changed in the volume occupied by the control rod since the scattering cross of the rod’s material is generally small. When the rod is inserted to the height z, the perturbation of the operator is

(4.6)

where Sa is the absorption cross section of the control rod. This means that, in Equation 4.3b, we have to integrate the numerator only for the inserted volume of the control rod while the integral in the denominator is extended to the whole volume of the reactor. Thus, the dependence of the reactivity on z is determined by the numerator. The denominator, Sa, and the other quantities independent of z can be condensed in a single constant C. With this, we get

.

(4.7)

If we plot this as a function of z, we get a curve like that presented in Figure 4.1. It has a typical shape recalling letter “S” for control rods for which the absorption cross section is independent of z. Such S-curves are regularly taken by the reactor operators. The curve has a relatively long section near zH/2 which is linear in a good approximation called the linear section of the control rod. The slope of this section is

.

 

This defines the physical meaning of coefficient C: it is the slope of the rod’s characteristic on its linear section.

Figure 4.1. Characteristic of a control rod according to Equation 4.7

(C = 0,05 $/cm; H = 100 cm)

5. The point kinetic equation and the effective delayed neutron fraction

            The point kinetic equation is a simple tool for discussing and interpreting time dependent phenomena. Its derivation is based on the assumption that the neutron density is separable:

(5.1a)

where y(x,t) is a function to be specified later. The point kinetic equation is related to the factor j(t). Such a separation is obviously arbitrary. It has some use only if the first factor is somehow restricted. Generally, the following is required:

(5.1b)

where n+ is the solution of Equation (2.9b). If we keep the time dependence of y(x,t), the parameters of the point kinetic equation (such as L and beff) will be also time dependent. In order to avoid this, the best choice is

(5.2)

where n is the solution of Equation (2.9a). With this, condition (5.1b) is automatically satisfied.

            We omit the derivations because they are lengthy and not too interesting. We restrict ourselves to the final result concerning beff. As a matter of fact, the main interest of the derivations consists in obtaining this formula, which is fundamental for reactor kinetics. We shall need some notations.

            Introduce subscript m which runs through the fissile and fertile isotopes (m = 235U, 238U, 239Pu etc.). fi(E) is the spectrum of delayed neutron group i, fp(E) is the spectrum of the prompt neutrons. All these spectra may be assumed as independent of m. All others, in turn, do depend on m: bim, nSfm. By means of these notations, we define the production operator describing the prompt neutrons emitted by isotope m:

(5.3)

where bm is the nuclear-physical delayed neutron fraction for isotope m:

.

(5.4)

In an analogous way, we define the production operators for the delayed neutrons (i = 1, 2, ..., 6):

.

(5.5)

The total production operator is thus given as the sum of all these operators:

.

(5.6)

            The static fission spectrum appearing in the operators (1.2c) and (1.3b) is a linear combination of these spectra. The latter should be averaged for all isotopes involved. The averaging factors are the integrals

(5.7a)

where

.

(5.7b)

Strictly speaking, both F and Fm depend on r but this is not shown in the formulae. Now the average static fission spectrum is

.

(5.7c)

            After some lengthy derivations, we get the following formula for the effective delayed neutron fraction:

(5.8)

where

.

(5.9)

            In order to apply these formulae, we have to solve the transport equation (or one of its approximations) and the corresponding adjoint equation. This is usually done for the critical state of the reactor. Formula (5.9) involves integration for the spatial variable r as well as summation (or integration) for the energy variable. Thus beff depends on the spatial and energy dependence of the neutron density and the adjoint function. The main effect is the latter i.e. the energy dependence: beff needs to be computed with a large number of energy groups. We show an example illustrating this statement. If the core composition of the considered reactor is simple like in a critical assembly, the spatial effects are much less important than the spectral ones. Therefore, the best is to treat them by asymptotic reactor theory. However, formula (5.9) should be applied explicitly for energetic reactors in which the fuel composition can change from fuel pin to fuel pin.

            Let us write Equation (5.9) explicitly for the energy-dependent diffusion model:

(5.10a)

where V is the volume of the reactor. Usually, such calculations are performed in multigroup approximation. We rewrite this formula in a multigroup picture for completeness:

.

(5.10b)

G is the total number of energy groups, vg is the characteristic speed of energy group g. The asymptotic approximation means the following:

                and                   .

 

If this is substituted in Equation (5.10a), the integrals for the variable r are the same in the numerator and the denominator leading to

.

 

This can be further simplified by virtue Equations (5.7a) and (5.7b):

.

(5.11)

In these inner products, the integration is done only for variable E.

            The spectra y(E) and y+(E) appearing in Equation (5.11) are computed by means of some asymptotic slowing down program. y(E) is the solution of an equation whose operator form is (cf. Equation (2.9a)):

(5.12a)

with

.

(5.12b)

The form of the destruction operator depends on the slowing down model actually used and on the way of treating the leakage. In order to have something in hand, we give its form for the simplest case i.e. for the case when the leakage is treated in diffusion theory:

(5.12c)

where B2 is the buckling of the fundamental mode. Equations (5.12) constitute an eigenvalue equation for y(E) which has a non-trivial solution only when keff and B2 take some special values. If, for example, we set keff = 1, B2 should be equal to the material buckling of the core composition. For any selected value of keff, the slowing down program computes the corresponding value of B2 and vice versa. The eigenfunction can be normalized in an arbitrary way. The most comfortable normalization is

(5.13)

since it results in the following slowing down equation:

.

(5.14)

            The advantage of this “trick” (which, by the way, is possible only in asymptotic theory) consists in that the slowing down equation has become an inhomogeneous integral equation which can be solved for any value of the buckling B2. If, for some given value of B2, we compute the solution y(E), Equation (5.13) yields the multiplication factor keff as a function of B2. This has a further by-product: the solution allows to compute the fission fractions Fm on the basis of Equation (5.7b) which reads as

(5.15)

in the actual model. According to Equation (5.7c), the spectrum f(E) depends on these fractions. Thus, strictly speaking, we ought to perform an iteration for this dependence. This is rarely done since its effect is rather small in most cases.

            Equation (5.11) can be applied only if the adjoint spectrum y+(E) is known. This can be eventually a severe problem because adjoint slowing down programs are not common place. Fortunately, this difficulty can be easily overcome in asymptotic theory. Let us write the adjoint slowing down equation (cf. Equation (2.9b)):

with

.

 

The adjoint function can be normalized similarly to Equation (5.13):

(5.15)

leading to the adjoint equation

.

(5.17)

We see from this that the adjoint function does not depend on the spectrum f(E). This allows to compute the quantity

 

rather easily. Equation (5.17) remains the same when we put the following factor in its right hand side since it is equal to 1:

.

 

This equation, in turn, is adjoint to the

(5.18)

equation. If its solution is normalized as

(5.19a)

(cf. Equation (5.13)), Equation (5.18) becomes

(5.19b)

(cf. Equation (5.14)).

            As a conclusion of our considerations, we formulate the following algorithm of computing the effective delayed neutron fraction:

1.    We first compute  i.e. that value of B2 for which keff = 1. The spectrum y(E) obtained in this step is used to compute Fm and F according to Equations (5.7a) and (5.7b), see also Equation (5.15).

2.    Keeping the value of B2 constant, Equation (5.19b) is solved for i = 1, 2, ..., 6, and the ki factors are computed by means of Equation (5.19a).

3.    The effective delayed neutron fraction for group i and isotope m is computed by using equation

(5.20)

(cf. Equation (5.11)).

4.    The effective delayed neutron fraction for group i is

 

while the total beff is their sum:

.

 

It is noted finally that Step 1 can be generalized: we first find that value of B2 for which keff is an arbitrary csosen positive number (not necessarily equal to 1).

            We give a numerical example obtained by a 40-group B1 slowing down program using the Grueling-Goertzel approximation. The example will illustrate how important the number of energy groups can be. beff was computed in four approximations:

Approximation 1: asymptotic slowing down in 38 groups out of which one is the thermal group;

Approximation 2: asymptotic slowing down in 6 groups out of which one is the thermal group;

Approximation 3: diffusion calculation of a bare reactor in 6 groups, using finite differences; the geometric buckling of the bare reactor was equal to the B2 obtained in Approximation 1;

Approximation 4: reflected reactor calculated in 6-group finite difference diffusion approximation.

The considered reactor was critical in all cases. For Approximations 1 and 2, the method described was applied while Equation (5.10b) was used for Approximations 3 and 4. Only 235U was considered as fissile isotope. The nuclear-physical value is b = 0.0065 for it. The obtained beff values are summarized in the following table.

            Approximation 1 yielded a large effect since the ratio of the effective and nuclear-physical delayed neutron fractions is

 

i.e. 11%. The difference between the values obtained by Approximations 1 and 2 is also large (2,9%). This means that we loose almost one third (» 2.9%/11%) of the total effect (11%) when the number of energy groups is reduced from 38 to 6. The differences between Approximations 2, 3 and 4 are much less ( 0,6% at maximum). It is remarkable that the effect of the reflector is practically negligible (Approximation 3 vs. Approximation 4).

 

Approximation

beff

1

0.007217

2

0.007016

3

0.006973

4

0.006984

 



[1] It is important to emphasize that the characteristics (r,E,) of the neutron are specified in the same moment t as the measurement is done. In case of the neutron importance to be defined later, these two moments will be different.

[2] Subscript “f” recalls final – just as subscript “i” recalls initial in Equation 1.4.