It is known that the diffusion (as well as the transport) equation has a well defined solution in V provided the entering current is given along the boundary dV. From the Green’s function and from the operators in (4.11) we set up the following iteration scheme. To formalize this, we write the solution as

Yk(r) = E Gkk'(r’ ^ r)4′ (r’)dr’.

V=1J dV

that can be put into the concise form

G

Jk = E Rkk’ Ik’, (4.25)

k’=1

where we have suppressed that the partial currents depend on position along the boundary and the response matrix R includes an integration over variable r’.

When volume V is large, we subdivide it into subvolumes (nodes) and determine the response matrices for each subvolume. At internal boundaries, the exiting current is the
incident current of the adjacent subvolume. Thus in a composite volume the partial currents are connected by response matrices and adjacency. We collect the response matrices and adjacency into two big response matrices:

J_ = RI; I = HJ, (4.26)

and because the adjacency is an invertible relationship, we multiply the first expression by H and get

I = HRI. (4.27)

Since there is a free parameter keff in matrix R, it makes the equation solvable. At external boundaries there is no adjacency, but the boundary condition there provides a rule to determine the entering current from the exiting current. With these supplements, the solution of equation (4.27) proceeds

I(m+1) = HR(m) l(m). (4.28)

The iteration starts with m = 0 with an initial guess for the keff and the entering currents I. Let us assume that the needed matrices are available, their determinations are discussed in the subsequent Subsection. The iteration proceeds as follows. We sweep through the subvolumes in a given sequence and carry out the following actions (in node m):

• collect the actual incoming currents of subvolume m.

• determine the actual response matrix to calculate the new exiting currents and contributions to volume integrals[17].

• determine the new exiting currents (J) from the entering currents and the response matrices using equation (4.26) and the contributions to the volume integrals.

After this, pass on to the next node. When the iteration reaches the last node, the sweep ends and the maximal difference is determined between the entering currents of the last two iterations. At the end of an iteration step, the parameter keff is re-evaluated from the condition that the largest eigenvalue of HR should equal one. If the difference of the last two estimates is greater than the given tolerance limit, a new iteration cycle is started, otherwise the iteration terminates. If we have a large number of nodes, the improvement after the calculations of a given node is small. This shows that the iteration process is rather slow, acceleration methods are required.

It has been proven (Mika, 1972) that the outlined iteration is convergent. The goal of the iteration is to determine the partial current vector. The length of vector I is Nnode x nf x G. From the point of view of mathematics, the iteration is a transformation of the following type:

A(keff )xm = aXm+1, (4.29)

where m is the number of the iteration, matrix A(keff) makes the new entering current vector xm+1 from the old entering current vector xm. In the case of neutron diffusion or transport, operator A(keff) maps positive vectors into positive vectors. In accordance with the Krein-Ruthman theorem, A(keff) has a dominant eigenvalue and the associated eigenfunction[18]. When keff is a given value, the power method is a simple iteration technique to find a good estimate of x = lim,^^, x,. Solution methods have been worked out for practical problems in nuclear reactor theory: for the solution of the diffusion and transport equations in the core of a power reactor. The original numerical method is described elsewhere, see Refs. (Weiss, 1977), (Hegedus, 1991).

Note that the iteration (4.29) is just an example of the maps transforming an element of the solution space into another element. Thus in principle one can observe chaotic behavior, divergence, strange attractors[19] etc. Therefore it is especially important to design carefully the iteration scheme. The iteration includes derived quantities of two types: volume integrated and surface integrated. When you work with an analytical solution, the two are derived from the same analytical solution. But when you are using approximations (such as polynomial approximation), it has to be checked if the polynomials used inside the node and at the surface of the node are consistent. In an eigenvalue problem, parameter keff in equation (4.29)should be determined from the condition that the dominant eigenvalue a in (4.29) should equal one. First we deal with the general features of the iteration.

As has been mentioned, one iteration step (4.29) sweeps through all the subvolumes. The number of subvolumes (Gado et al., 1994) varies between 590 and 7980, the number of unknowns is 9440 and 111680. At the boundary of two adjacent subvolumes, continuity of Ф and DdnФ (the normal current) is prescribed.

In node m in iteration i. In the derivation of the analytical solution we have assumed the node to be invariant under the group Gy. Actually, not the material properties are stored in a program because the material properties depend on:

• actual temperature of the node;

• the initial composition of the fuel (e. g. enrichment);

• the actual composition of the fuel as it may change with burn-up;

• the void content of the moderator;

• the power level;

In the calculations, app. 50-60% of the time is spent on finding the actual response matrix elements, because those depend on a number of local material parameters (e. g. density, temperature, void content). We mention this datum to underline how important it is to reduce the parametrization work in a production code.