In the case of internal boundaries we have a condition of continuity of the angular flux ф(г, E, ft), which in the case of the multi-group diffusion equation is equivalent to continuity of flux and current for each group. Because of Fick’s law the current continuity means

D № D,^

dn dn

where n is in a direction normal to the surface and the suffixes 1 and 2 refer to the regions on the two sides of the boundary.

In the case of continuous refuelling, after a running-in phase an equilibrium is reached, and the reactor composition does not change with time. In this equilibrium condition all burn-up stages of the fuel are present at the same time in the reactor (uniformly graded exposure). The mean reactor composition is then an average of the compositions of all burn-up stages of a fuel element. If Nk(t) is the concentration as a function of time of isotope к in a fuel element, the average concentration Nk of this isotope in the reactor is

where T is the residence time of the fuel elements in the reactor. This expression supposes a very high number of fuel elements in the reactor and continuous refuelling. This condition is exactly satisfied for pebble-bed reactors and still sufficiently well satisfied in the case of prismatic fuel with on-load refuelling.

The existence of an equilibrium condition supposes a constant power level of the reactor so that the time t of eqn. (9.2) actually corresponds to an irradiation tф. Periodic power variations due to load following do not disturb the equilibrium condition if the time period is short compared with the decay constant of the isotopes considered. In eqn. (9.2) the concentrations Nk(t) must then be calculated using an average power level. This approximation is valid for most isotopes. The most important exception is given by l35Xe. The concentration of this isotope is usually calculated separately considering the operational requirement of the reactor (see Chapter 12).

These equilibrium calculations are usually the first calculations performed when designing and optimizing a reactor. In this case one has to iterate in order to obtain a critical reactor. Usually a calculation is started with a guess on the initial composition and the average composition is obtained. This can be done for each reactor region separately (e. g. inner and outer core) or for the whole reactor (point model calcula­tions). Once the equilibrium composition is known a criticality calculation is performed. In the case of the point model

У. 2/fc

fee* =—- !——— — (9.3)

к

where the summation is extended over all the isotopes k. In the case of a multi-region reactor calculation, a diffusion calculation is performed.

If кея is different from the desired value (usually slightly above 1 for control reasons) one iterates on the concentration of the fertile material of the fresh fuel (or in some cases on the fuef residence time) until this value is reached.

The %ak and %fk appearing in (9.3) are obtained as

lk=jjk(t)dt

in analogy to (9.2). Here the 2k(f) are macroscopic cross-sections including self — shieldings. These 2k(f) must be calculated from the concentrations obtained by solving the depletion equation (9.1), taking into account the cross-section variations due to spectrum changes during the fuel lifetime.

Here one can distinguish between different cases. If the fuel elements are small compared to the neutron migration length, the neutron spectrum is determined by the average reactor composition and remains constant in the equilibrium phase, indepen­

dently of the burn-up of the single fuel element which is being studied. In this case the constants of eqn. (9.1) can only be calculated when the average composition given by

(9.2) is known.

A double iterative procedure is then necessary. First with a guess of the average composition the cross-sections are calculated and eqn. (9.1) is solved. Now it is possible to obtain the average composition from eqn. (9.2) from which new cross-sections for eqn. (9.1) are calculated. Once this iteration (inner iteration) has converged, an iteration (outer iteration) on the fertile concentration (or on the residence time) will start to obtain criticality, and the whole cycle of calculations is repeated.

In order to deal with reprocessing, at the end of each inner iteration the fresh composition for the next iteration consists of all the U (or Pu) isotopes contained in the discharged fuel (whose concentration is multiplied by a reprocessing efficiency) plus enough feed fuel to have either the same moderator to fuel ratio of the first given loading, or the same life average moderator-to-fuel ratio. This method can also be applied to treat two intimately mixed sorts of fuel (feed and breed) out of which, in some cases, only one type is reprocessed (fuel segregation). Those two types of fuel could possibly have two different residence times. If the fuel elements are so big that their neutron spectrum is determined only by their composition, the above described inner iterations are no longer necessary. This does not normally happen, but in big fuel blocks the neutron spectrum is usually influenced both by the composition of the block under study and of the surrounding blocks. An exact treatment of this would require a two-dimensional burn-up calculation. Usually first survey studies are performed with the simpler methods here described and later more detailed calculations are performed in which the effect of neighbouring fuel blocks, control rods, reflector blocks, etc., can be properly treated. An interesting way of taking into account the effect of neighbour­ing blocks is given by the use of collision probability methods in the multi-cell option of the WIMS code. These energy-dependent probabilities, which are used to couple the individual cell calculations, are usually obtained from diffusion theory calculations on the regions of interest (see also §4.11).

If the parameters of eqn. (9.1) can be considered as constant, analytical solutions are possible with a considerable saving in computing time in the calculation of equilibrium cycles. This supposition requires beside the Constance of the spectrum, also a Constance of the self-shieldings. As outside the resonances self-shieldings are very near to unity in HTR fuel, this assumption is often satisfied. Besides the concentration of fertile material, and hence its self-shielding, changes very little during burn-up.

A notable exception can be given by the Pu resonances in the low enriched-uranium fuel cycle where, as the Pu concentration changes with time, the approximation of constant self-shieldings can be rather crude.

Particularly interesting in the case of analytical solutions is a method developed by Blomstrand which, using Fourier transform of eqn. (9.1), relates the average composition to the feed composition, avoiding the inner iterations described above. In the BASS and BABS codes"21 written according to this method it is possible to specify directly the characteristics of the average equilibrium core composition (e. g. its moderator to fissile ratio S) and obtain the required feed composition.

For more detailed equilibrium calculations space dependence has to be considered. In space-dependent equilibrium burn-up codes it is possible to treat various core regions iterating on the feed fuel composition, or on the burn-up, of these regions in order to obtain a flat power distribution (e. g. same maximum power in all regions). The one-dimensional FLATTER code,’131 which uses the Blomstrand method combining the BASS code with a one-dimensional diffusion calculation, can divide the core in a number of burn-up sub-regions (independent BASS calculations). These sub-regions can be combined in two or three macro-regions whose initial composition or burn-up can be changed in order to achieve power flattening.

The energy E appearing in the Breit-Wigner formula (3.1) is the relative energy between target nucleus and neutron. This means that the thermal motion of the nuclei of reactor material should be considered and for the calculation of all reactions taking place in a reactor the relative energy should be used. In practice one uses effective cross-sections in which only the neutron energy appears. These effective cross-sections are, of course, temperature dependent, and must be defined in such a way as to give the proper reaction rate. Therefore we must have:

vacW(v, T) = j VrdO-(Vrd)P(vrei, T) dvrc,

where v = neutron velocity,

vr<!, = relative neutron-nucleus velocity, va = velocity of absorber nucleus,

P(prd, T) = probability of having the relative velocity vrel,

T = temperature.

Considering that and that

P(va, T) = M(va, T) Maxwellian distribution of nucleus velocities

we have

where E = energy corresponding to velocity v,

A = mass of the absorber in units of neutron mass.

One can demonstrate that if cr(vIci) follows the 1/u law, also cr^s will follow the 1/u law, independently of the temperature T. Of course <7ея-мт if T-»0°K. In the case of resonances, temperature has the effect of lowering the maximum of the resonance and broadening it at the same time (Doppler broadening).

Substituting the Breit-Wigner expressions in (3.5) one can obtain the Doppler broadened formulation of the resonance cross-section (the suffix eff used in (3.5) will be dropped here).

(Ta =-^сгоф(С, к) absorption cross-section,

Г„ . .

CTS = уг (Тоф + <T s, pot + <T s, int scattering cross-section,

with

crs, int = interference scattering, crs, pot = potential scattering.

і/»(£, K) and xU> x) are the so-called shape functions

. 2(E — Eo)

where к = ————-

Г

С ————— Г75 = тг — ratio of natural to Doppler width,

1 D

E0 = resonance energy.

We have seen that for narrow resonance approximation the flux is given by (7.4)

2Po. + 2,i 1

2,o + 2,i + 2o E

2„ 1

2,i + 2o E

If the resonance is neither narrow nor wide one can use the intermediate form, first proposed by Goldstein and Cohen,<2>

(7.8)

where

When A = 1, eqn. (7.8) gives the NR expression and when A = 0 eqn. (7.8) gives the NRIM approximation.

With NR approximation for the moderator eqn. (7.3) takes the form

where the integral operator Кф is defined by

From eqn. (7.9) it is possible to define the iterative sequence

(7.10)

where фт, фа . . ., ф(п) are successively closer approximations to the correct solution. The first guess can be taken from the NR or NRIM approximation and the iteration can be repeated until ф(а> converges to the solution of eqn. (7.9). However, these iterations become quite laborious. The intermediate resonance (IR) method consists of inserting a

given by фт is equal to I<2) given by ф(2

Usually a reactor optimization is first performed on the equilibrium composition and only in a second stage the running-in calculations are performed.

All ordinary burn-up codes can be used for running-in calculations, except those which have been particularly developed for equilibrium calculations. Point model calculations are sometimes used for a first assessment (e. g. the codes GARGOYLE’9’ and HELIOS’10’) but they cannot give any information on the power distribution so that one — and eventually two-dimensional codes have to be used.

The one-dimensional RINA code of Ispra’111 uses a nodal approach for the solution of the radial diffusion equation in one energy group. For a given refuelling interval the programme searches for a single value of the initial fissile concentration in one or more zones, necessary to obtain an imposed reactivity at the end of the interval, or searches for the refuelling time when the initial composition is fixed in each region. Burnable poisons can be considered.

Two-dimensional calculations in which an R-Z flux map is synthetized from

one-dimensional calculations are used in the VSOP code system/12’l3) Usually a first survey of the various possible running-in procedures is performed with fast zero or one-dimensional codes. Each of these calculations is accompanied by a fuel-cycle-cost calculation and by an evaluation of the constraints (temperatures, burn-up, control-rod requirement, temperature coefficient, etc.).

In this way an optimization of the running-in strategy can be obtained. In a later stage more detailed calculations can be performed with two — or three-dimensional codes.

Let us consider the time-independent Boltzmann equation (4.4) with isotropic sources (again the fission source has been for simplicity included in S).<l0)

,ft)= I ls(E’^E, Cl’

+ 4l~S(r, E)

С1)ф(г, E’, ft’) dE’ dil’

expanding the scattering cross-section of scatterer і in Legendre polynomials

<tAE’^E, CI’^CI) = 2 s„(E’-»E)Pi(#io) (4-42)

і

where p0 = ft ■ Cl’ — cos во and defining

XAE’ -*• E) = 4тг 2 NME’ -*■ E)

і

where N, is the atomic density of scatterer i, we obtain

a • viKr, e, a) + х, ф(г, в, п) = — J — 2 f s„(E’-*• е)р,(р0)ф(г, e, ft) dE’ dcr

4тг i J

+ 7— S(r, E). (4.43)

47Г

The one-dimensional form of this equation is sufficient for all interesting applications of the Bn method.

If z is the space variable, instead of the variable ft it is sufficient to introduce the variable p = cos в where в is the polar angle between the direction ft and the z-coordinate.

Because of symmetry reasons the flux is then independent of the azimuthal angle <p and we can use the variable

1A(z, E, p) = j^ 1jt(z, E, ft) d<p

so that we have

Making use of the addition theorem of Lengendre polynomials it is possible to express Pi(po) in terms of p and p’

P,(po) = Pl(p)Pl(p’) + 2’Z,(j РГ (P)РГ (p’) cos [m(<p — <p’)]. (4.45)

m=l 1» ~T tTl )l

Carrying out the integration with respect to <p’ from 0 to 27t the second term of (4.45) disappears so that (4.44) becomes

+ ^S(z,£).

We multiply this equation by e‘Bz and integrate over all z obtaining an equation for the Fourier transform of the flux densityt

Equation (4.46) becomes

The Fourier transform of the flux density can be expanded in Legendre polynomials

Ф(В, Е, р) = їф,(В, E)P,{p) (4.49)

l

dividing both sides of the equation by

1 — iyp, where у =

multiplication by Pj(p) and integration over p, because of the orthogonality of the Legendre polynomials, gives

In the Bn approximation the expansion in Legendre polynomials of the scattering cross-section (4.42) is truncated at the nth term. The very important peculiarity of the Bn method is that after having made this truncation no truncation in the expansion of ф

tThe reader is reminded that not all texts use the same definition of a Fourier transform. The definition (4.47) implies an inverse transform:

ip(z, E, ju.) = J e, Bz ф(В, E, ці) dB. Sometimes the following definitions are used:

ф(В, E, ці) = J e IB‘ ip(z, E, ці) dz, ф(г, E, ft) = 2~ J е>Вг Ф(В, E, ці) dB.

A symmetric definition is also possible

ф(В, Е,р)=уІe~,B‘ ф(г, E, p) dz, іl/(z, E, fi)= j е‘Вг ф(В, E, ці) dB

is necessary. The precision with which the components фі of the neutron flux density are calculated is only influenced by this truncation. In general the B„ method gives more exact values for the components ф; than the P, method because neutrons emerging from collisions are more isotropic than the flux (or neutrons entering a collision). If the scattering, for example, is linear in the laboratory system the В, method is exact, and if scattering is isotropic in laboratory system, B0 is exact. The B„ approximation starts from the angular distribution that neutrons would have with isotropic scattering (in the case n = 0) and calculates deviations from this distribution. Besides the anisotropy of the scattering cross-section is a known quantity so that it is easier to assess the value of the approximation which is being made with the B„ method.

The difficulty of the B„ method consists in the fact that in order to obtain the flux it is in principle necessary to perform a Fourier inversion. The only way of avoiding this inversion is to assume separability of the space and energy dependence of the neutron flux

ф(г, Е, B) = WB, E, (l)eiBr

in this case the energy dependence of the Fourier transform is the same as the one of the flux.

While this approximation would be too crude for complete reactor calculations, it may still be quite useful in the case of spectrum calculations where it is often used to obtain within-group fluxes (see §8.1).

equation (see §4.11). Integral transport theory is also used in the General Atomic MICROX code<10> which solves the neutron slowing down and thermalization equations in 200 energy groups for a two-region lattice cell. The fluxes in the two regions are coupled by collision probabilities computed with flat flux approximation. Grain structure of the fuel is taken into account with the method developed by Walti (see §7.13).

The leakage in this code is treated with a diffusion formulation derived from an approximation of the В, equations. The term

which appears in the Вi equations [see eqn. (4.50)] is substituted by

where фітсі(Е) is a fixed reference current spectrum (since in the Pi approximation the ф, component can be interpreted as a neutron current). With this assumption it is possible to eliminate ф, and have an equation for ф0 where the leakage term can be expressed in the form О(Е)В2ф0(Е) (see ref. 10, pp. 35-39). In this way the Bt equations are transformed into the multi-group diffusion equation and the transport cross-section is calculated weighting the scattering matrices on the reference current spectrum. It must be noted that this approximation is rather similar to the one made for obtaining the multi-group diffusion equation from the Pi equations. The use of energy-dependent bucklings becomes then possible. As expected this approximation is rather crude for hydrogen but sufficient for heavy nuclides. As in the case of the MUPO code this allows a much faster calculation. The term D(E)B2(E) refers to the cell averaged leakage and is added to the removal cross-sections of each region treated in this code.

Another peculiarity of the MICROX code is the replacement of the elastic and inelastic slowing-down kernels of each individual nuclide by synthetic scattering kernels of mathematically more tractable form. This permits a very rapid evaluation of the slowing-down integrals. The scattering kernel 2s(u->u’) is expressed as

2s(m -» и’) = ^ Nkak, s(u)g(u -> u’)

where и = lethargy, Nk = atomic density of nuclide k, ak, s = scattering cross-section of isotope k, and the distributions gk are defined by

gk(u -» u’) = #(u’ — u)R'[Zkpke-pt, u’-u>

where #(x) is the Heaviside function (= 1 for x > 0; =0 for x < 0). The parameters Zk and pk that characterize the function gk are complex numbers chosen in such a way as to conserve the number of neutrons in a scattering collision, to give the correct mean lethargy gain per collision, and to have the minimum mean square deviation between the synthetic and the actual kernel (for details see ref. 10).

Another method of performing spectrum calculations in heterogeneous cells is given by the WIMS code. The last versions of WIMS have actually evolved to a modular system including almost all types of neutronics design calculations, but we will only deal here with the spectrum calculations."11

The programme includes a sixty-nine group library. Between 4 eV and 9.118 keV are thirteen resonance groups whose cross-sections are calculated by means of the equivalence theorem (see § 7.11) using a library of resonance integrals. This library of resonance integrals is obtained by solving the slowing-down equations for a homoge­neous mixture of moderator and resonance absorber with the SDR code in some 12,000 energy mesh points. This sixty-nine group library is then condensed to fewer groups using the SPECTROX technique developed by Leslie for the calculations of heavy water systems."21 In this technique each of the principal regions of the lattice are coupled by means of collision probabilities.

The collision probability equations can be modified to account for flux shape in large moderator regions supposing that the flux rise in the moderator at any energy is proportional to the current leaving the moderator at that energy, the constant of proportionality being taken from one-group diffusion theory. After having condensed the data to a lower number of groups the WIMS code can perform a cell transport calculation with either collision probability or S„ methods.

A unique feature of WIMS is the possibility of expressing the effect of the environment on the neutron spectrum of the region considered through collision probabilities (MULTICELL approximation)."31 The multi-cell model describes the core of a reactor in terms of a number of different cell types, each being weighted in proportion to its frequency of occurrence in the system. The cells are effectively combined for a single collision probability transport solution by specifying the matrix of probabilities that a neutron leaving each cell type will enter a cell of the same or any other type."41

The simplest assumption is that these probabilities are directly proportional to the contiguous surfaces between the cells. The multi-cell WIMS method assumes that the emission density is uniform in space and angle across the region associated with each cell type. The spatial distribution is, however, often significantly non-uniform when regions are more than one or two mean free paths in thickness. For these situations a coarse mesh correction has been developed"51 in terms of a multiplicative factor on the surface-to-surface probabilities. The simplest correction factor is a constant factor /3 (=0.7 say) to reduce the overall leakage."61 Improved values of /3 can be obtained by an iterative process with a space-dependent few-group diffusion calculation for the whole reactor. This can form a possible alternative to the buckling iteration. The advantage of the multi-cell approach is that very large meshes can be used in regions of constant spectrum, allowing a finer mesh in the regions where spectrum variations are important.

A multi-cell layout for an HTR core is represented in Fig. 8.1.031

Heterogeneous cell calculations are performed in France with the APOLLO code"71 which calculates the space and energy-dependent flux for any one-dimensional medium with collision probability methods. Various options are possible. Collision probabilities can be calculated with the flat flux assumption in each region, or in large media the flux can be considered to vary linearly within each region.

Anisotropic collision probabilities can be used to treat linear anisotropic scattering. The multi-cell formulation of WIMS has been included in APOLLO to provide a simplified treatment of two-dimensional geometries. Two libraries are available with 186 or 99 groups.

Assuming Ak = const, and independent of feedback effects, simple analytical solutions can be found. These solutions are usually uninteresting for power operation, where Ak is a function of temperature, but can adequately represent the reactor operation at very low power. In this case, neglecting the source S, we have a homogeneous system of m + 1 linear differential equations the solution of which is of the form

<Hf) = X A, e

і =0

where m is the number of delayed neutron groups being considered. The coefficients A, are determined by the initial conditions. The exponents to, are obtained as the roots of the characteristic equation

ІШ + 1 у ti)3,

1ш + 1 /to + 1 frl to + A, where p = Д/c/fceff is the reactivity.

Introducing the period T = 1 /to this equation takes the form

which is the well-known inhour equation (so called because it can be used to define a reactivity unit, the inhour, which is the reactivity which produces a reactor period of 1 h). Equation (12.14) has m + 1 poles, of which m occur for values of to at which to + A, vanishes (to = — A,) and the (m + l)st pole occurs at to = — 1//.

The roots are all real; m of them are negative and lie between the poles. The remaining root to0 is algebraically the largest one and has the same sign as p.

Some time after the beginning of an excursion, the exponentials due to the first negative m roots will have vanished, and only the last exponent to0 will remain to describe the reactor behaviour. The period T = l/to0 is called the stable reactor period. In the case in which S is not negligible (e. g. reactor start-up) an analytical solution can still be found adding to the solution (12.13) of the homogeneous system a solution of the inhomogeneous system.

6.1. Slowing-down in graphite

Fission neutrons are emitted with an energy spectrum ranging from ~ 10 to 0.5 MeV.

Collisions with moderator and heavy metal atoms slow down the neutrons until their energy reaches an equilibrium with the thermal motion of the medium. For light nuclei, like carbon, only elastic scattering takes place at energies of interest for nuclear reactors and the collision can be studied as a classical mechanics problem.

The collision can be considered both in the centre of mass frame of reference and in the laboratory frame of reference, the first one being only a means to calculate the data in the laboratory system. The problem of the slowing down by elastic collision is treated in most reactor physics textbooks, so that we only quote here the most significant results. Defining іjj the scattering angle in the centre of mass frame of reference, we can obtain (see, for example, ref. 1, § 7.1.1)

(6.1)

where A = atomic weight of scatterer in units of neutron mass,

v’ = neutron velocity in the laboratory frame of reference before collision, v = neutron velocity in the laboratory frame of reference after collision,

E’, E = neutron energies corresponding to the velocities v’ and v.

The neutron loses most energy when its direction of motion is reversed (cos ф = -1) and thus:

(6.2)

It is obvious that even in the case of isotropic scattering in the centre of mass system, the scattering is anisotropic in the laboratory system, except for heavy nuclei where the two systems practically coincide.

In the case of isotropic scattering in the center of mass system (which is usually the case below 1 MeV) all values of E are equally likely in the interval

aE’^E^E’

tThe present definition of a is the most commonly used in the literature. Unfortunately in the past some authors (e. g. Wigner, Weinberg and Nordheim) defined a as 4AI(A + l)2 so that a is replaced by (1 — a) in every formula.

so that the probability of having a final energy in the small interval dE about E is

when aE’ E E’ and zero when E > E’ or E < aE’.

The scattering kernel is then given by

Z,(E’^E) = Xs(E’)g(E’^E).

The average cosine of the scattering angle in the Laboratory system is given by

cos в = (6.3)

The average logarithmic energy loss per collision is given by

f-lnf-f ln^g(E’^E)dE = +T^—na. (6.4)

.C J aE’ & A OL

The average number of collisions from fission energy Ef to thermal energy Elh is

We can now compare the different moderating properties of materials used in nuclear reactors (Table 6.1).

Table 6.1 H D Be C 0 U A 1 2 9 12 16 238 a 0 0.111 0.640 0.716 0.778 0.983 і 1.000 0.725 0.209 0.158 0.120 0.00838 n (2 MeV -> 0.0253 eV) 18 25 86 114 150 2172

We see that graphite, because of its high atomic weight, is a rather poor moderator, but it has the advantage of a very low absorption cross-section, so that if we take as representative the parameter £(2s/2„), graphite is only inferior to heavy water and beryllium.

As we have seen, the average logarithmic energy loss per collision is independent of energy: this suggests the convenience of using a logarithmic energy scale, called lethargy (и), defined as follows:

и = In Щ (6.6)

where E0 is the highest energy considered, usually 10 MeV.

Considering the slowing down in an infinite homogeneous medium we can write the following balance equation:

Х,(Е)ф(Е) = S(E) + fE Xs(E’^E)4>(E’)dE’ (6.7)

which is the slowing-down equation.

The slowing-down equation is actually a simplification of the Boltzmann equation

(4.2) from which it can be obtained supposing time- and space-independence and integrating over all angles. This equation can always be numerically solved using the multi-group method. For spectrum calculations the Boltzmann equation in a more complete form is normally used (P,, B„ or diffusion approximation), without the need of supposing space independence (see Chapter 8).

Considering only slowing down in the moderator we have for elastic scattering

X,(E)<fi(E) = S(E) + 1-!— [ ‘ ХЛЕ’)ф(Е’)Щ^. (6.8)

1-а JE £J

If %S(E) = const and Xa = 0 at the energies below the fission spectrum where S(E) = 0 we obtain as a solution a 1 IE dependence of the neutron flux. This is a very important result because the scattering cross-section of most moderators is constant over a very broad energy range.

In case of off-load refuelling the equilibrium composition and reactivity are not constant, but periodical, the period being the interval between two subsequent reloadings. In reality the period can be longer than this interval because the refuelling of one region is not equivalent to the reloading of another one, but usually this is only considered in detailed calculations and not in survey studies. In a zero-dimensional case the average reactor concentration Nk{t) of isotope к as a function of time during the refuelling period of length Tp is given by

і M-1

M(D = tt 2 Nk(t + nTP)

1*1 n =0

Fig. 9.6. Concentration of isotope к as a function of time.

The calculation is performed with methods similar to those used for the continuous refuelling. The reactivity at the end of each period must have the values specified as the minimum acceptable for control and operation reasons. In order to reach this an iteration on the feed-fuel composition can be performed.

A typical zero-dimensional code of this type is the General Atomic GAFFEE code which uses analytical solutions/14’ A modified version (MOGA) developed in Ispra can deal with concentration-dependent self-shieldings.

For analytical solutions a further approximation is needed because in the case of discontinuous reloading the flux does not remain constant during the period between two subsequent reloading. An analytical solution of eqn. (9.1) requires constant flux. It is a reasonable approximation to use the flux corresponding to the time average of the concentrations during the period.

As for the case of continuous refuelling, also for off-load refuelling the methods used for zero-dimensional equilibrium calculations can be extended to deal with space dependence. Using one group diffusion calculations based on nodal methods the GAFFEE type of calculation has been extended to one and two-dimensional calcula­tions (codes RACE and TRACE’15’).

9.2. Burn-up calculations in pebble-bed reactors

A special problem is posed by the burn-up of pebble-bed reactors. In this case the spherical fuel elements are loaded on the top of the reactor, flow across the core and are discharged from the bottom.

In the pebble-bed reactors now in operation or under construction (AVR and THTR) the fuel elements flow many times through the reactor before being finally discharged. Reactors in which the fuel elements flow only once through the core are being also considered for future applications (the so-called OTTO fuel cycle). In the normal case with many passages of the fuel elements, pebble-bed reactors are the perfect example of continuous reloading.

In any volume element of the core all burn-up stages of the fuel are represented and, because of the small size of the fuel elements, the neutron spectrum is determined by the average composition. The burn-up of the elements flowing out of the bottom of the core is measured (usually by means of a small reactor,(16>) the elements whose burn-up is higher than a given limit are discarded, the others are reloaded. A two-region core for radial power flattening can be obtained loading preferably fresh elements in the outer core and older elements in the inner parts.

The problem is complicated by the radial dependence of the axial fuel velocity, which is usually higher at the core centre and lower near the reflector.

First assessments are usually made with zero-dimensional burn-up codes for continu­ous refuelling, in which the effect of the axial fuel element movement does not need to be considered. Later more detailed space-dependent calculations are needed in which the peculiarities of this reactor type are taken into consideration.

Space-dependent burn-up calculations for pebble-bed reactors can be performed with standard burn-up codes treating a very high number of regions and considering each region as composed of a high number of fuel element classes representing different burn-up stages (and possibly different fuel-element types). The burn-up of each class must be followed independently of the others. The axial fuel element movement can then be simulated with a very high number of reloading-reshuffling operations in which the composition of each axial zone is shifted into the next one. If necessary different spectral regions can be considered. For each of these regions the cross-sections for the burn-up equations are calculated with spectrum calculations.1171 This method can give very accurate results but is not suited for fast calculations.

The simplified methods which have been developed for this reactor type have proved in practice quite fast and accurate. Supposing that the fuel moves only in the axial direction it is possible to write the general burn-up equation’181

v _ ф 2 ]sjj(Tfiy. k + Ф 2 Nscrasysk + 2 NAjCXjk — ЛkNk — 4>Nk(Tak (9.4)

Ol OZ i=l s=r

where Nk = atomic concentration of isotope k, v = axial ball velocity v = f(r, t), ф = flux, ф = f(r, z, f),

<Jfi = isotope і fission cross-section,

(Tai = isotope і absorption cross-section,

Лі = isotope і decay constant,

yik = yield of isotope к due to a fission in isotope i, ysk = probability that a neutron absorption in isotope s produces isotope k, a, k = probability that the decay of isotope і produces isotope k, r = radius,

z = axial coordinate—the balls are loaded in the point z = 0 and move in the direction of increasing z,

H = depth of the pebble bed.

This equation is identical to (9.1) with the addition on the left-hand side of the term

dNk

~z— v dz

taking into account the axial movement of the fuel. In equilibrium condition

and v = f(r), ф = f(r, z) are supposed to be known. In this case it is possible to solve numerically eqn. (9.4), if the boundary condition Nk(r)z.0 is given for all k. The difficulty consists in relating Nk(r)z,0 with Nk(r)z. H because the composition of the discharged balls is not known.

A simple way of solving this problem is to assume that the composition of a fuel element is fully defined by its irradiation. This method is used by the KUGEL code.’19’ This programme calculates the statistical space distribution of the fuel elements in a pebble-bed reactor taking into account the radial dependence of the passage time through the core. A fuel element is followed during all its passages through the core until it is discharged and the probability of finding it in all points of the core is calculated as a function of its irradiation.

The composition of the fuel element as a function of irradiation (here always measured in full power days) is read in by KUGEL in form of a table. In this way the burn-up equations do not appear in the code. This table is obtained from a zero­dimensional burn-up code which calculates numerically the burn-up of a fuel element in a constant flux taking into account the most important fission product chains. The flux level is the one obtained from a zero-dimensional burn-up calculation for continuous refuelling (e. g. the BASS code"2’) which must have been performed in a previous assessment stage. The KUGEL code divides radially the reactor in a certain number of axial channels.

In each channel the pebble velocity is supposed to be constant and it is assumed that the movement is restricted to the vertical direction. The total irradiation time of a fuel element is subdivided in a discrete number of age groups of width т (this age is an irradiation (ф x t) and not a time, but it is convenient to measure it in full power days). The fresh elements belong to group 1. Usually between 100 and 200 age groups are considered. An age spectrum A(f) is defined as the probability of finding a fuel element as function of its age t. If the age is discretized in age groups, the age spectrum is a vector A with as many components ak as age groups k. t At each passage through the core the irradiation of a fuel element increases by a fixed amount h which is

tin the following treatment we define with a small letter the components of a vector defined by the corresponding capital letter, so that ak are the components of vector А, Ь’л the components of vector Bl, etc.

characteristic of each channel і (proportional to neutron flux and inversely proportional to the axial velocity). This means that an element belonging to age group к if loaded in channel і will, at its exit from the core, belong to age group к — t — Sr, where 8, is a “delay” characteristic of this particular channel

(9-5)

t ф

where і = channel index,

tt = passage time in channel i, фі = average flux in channel i, ф = average core flux

(the discretization of the age groups means also that Si will have to be rounded off to the nearest integer).

At the beginning of each passage j it is possible to define a “loading spectrum” L’ whose components /*’ give the probability that a fuel element, at the beginning of its jth passage through the core, belongs to age group k. At the first passage

In each passage the fuel elements are distributed among the various channels according to their age, following a given distribution law. This distribution law should insure a constant level of the pebbles in all channels and a flat radial power distribution.

Let us define W‘ as the probability that a fuel element in passage j is loaded in channel i, and

(9.5)

where j = passage index, J = total number of passages.

In order to have the same level of pebbles in each channel the probability p, should be proportional to channel cross-section and inversely proportional to channel passage time

St It,

2 s, it,

where Si = cross-section of channel /, N = total number of channels.

In order to obtain a radially flat power distribution it is usually necessary to define two core zones. The inner zone includes channels 1 to R and the outer zone channels R + 1 to N. In this case it is possible to define

wi,

і = 1 І = 1

Pt and P2 are the probabilities that a fuel element (in the total of its passages) comes in the inner or the outer zone respectively. We then have:

One could assume that all elements whose age is smaller than T* are loaded into the outer region, those whose age is greater than T* are loaded into the inner region and those with age greater than T are discarded.

Defining

J К*

a? _

<=£ — j к

2 2 ь’

= 1 k = 1

where X* and К are the age groups corresponding to T* and T respectively. It would then result

P, = 1 — 5£, P2 = $.

К’ or T* would be defined by condition (9.8) and no parameter would remain free to allow for power flattening.

In order to avoid that, one can introduce two corrections.

ai Probability for an element of age < T* to be loaded into the inner zone.

a2 Probability for an element of age > T* to be loaded into the outer zone.

It follows:

P, = a,^ + (l-a2)(l-^),

(9.9)

P2 = (l-a,)^ + a2(l-^).

Pi and P2 must still satisfy the conditions (9.8), but now two parameters (ai and a2) are free for adjustments, is only known if the problem of the fuel element distribution as function of age has been solved, i. e. it can be only obtained through an iterative process. As first guess for the iteration it is assumed

(This formula would be valid for infinite axial velocity.) The loading spectrum L’ is subdivided into two parts Lx and L2‘ belonging to the inner and outer core zones respectively. For an age group k, corresponding to an age smaller than T we have:

and for an age group k corresponding to an age greater than T we have:

Within each core zone the distribution is then proportional to the channel cross-section and inversely proportional to the passage time. The age spectrum B’ at the entrance of each channel і for each passage j can be obtained from the loading spectrum L1 according to the distribution law.

The elements of the age spectrum В/ are obtained in the following way:

ik t 1 к R

2 Silt,
and

■ S It-

Ь[к= 12k—n——1— for the outer zone. (9.12)

2 S, lt,

i=R + 1

Let us define Ql as the age spectrum at the end of channel і for the passage j. Its elements q’ik can be obtained from the elements of the age spectrum В/,

qik+s,= b’,k. (9.13)

Out of the age spectrum Ql the part belonging to an age < T is used to produce the loading spectrum for the passage j + 1,

llk+’ = ‘Zq[k for k = K — (9.14)

1 = 1

The part corresponding to an age greater than T forms the discharge spectrum D,

dk=’Z’Zq[k for all k>K (9.15)

j-i i=i

defined as summation of the discharged elements over all channels N and all passages J.

The integrated loading spectrum

(9.16)

к

gives for each passage j the probability of still finding a fuel element in the reactor at the jth passage. The calculation is continued until yJ+’ =0.

J is then the greatest number of passages which a fuel element can have through the reactor.

Having started with a loading spectrum normalized to 1, also the discharge spectrum must give:

2 dk = 1.0.

к

The average age of a discarded element is

Td = t 2 kdk. (9.17)

к

Until now we have assumed that it is possible to measure exactly the age of a fuel element.

In reality one has a Gaussian distribution of the measured values so that eqns. (9.8) to 9.11) have to be modified in order to take into account this effect.

The probability that a fuel element of age T — e is classified as having an age greater than T is

-5= [ e~h2x2dx, (9.18)

V 7Г J.

where h is the precision index. This same expression gives the probability that a fuel element of age T + e is classified as having an age lesser than T. (It is assumed that h is the same for all age groups.)

Let us define

gn=~y=( e 1,1×1 dx, (9.19)

V7T J„T-T,2

g„ is the probability that a fuel element which is classified as belonging to an age group equal or greater than k, belongs in reality to group к — n (or the probability that a fuel element which is classified as belonging to an age group equal or smaller than k, belongs in reality to group к + n). We have then a vector G whose components are

go, g 1, gl ■ ■ ■ gs

assuming that all g„ with n > s are negligible.

Equations (9.14) and (9.15) are then substituted by

N

/i+, = 2qU for k^K-s, (9.20)

і = I

a = 2 (9.21)

І=I i=l

N

and (i+‘ = 2 q[kG~gK-k) for K-s<k^K,

і — 1

A =2 2 41*0 “ft-*). (9.22)

І~• i=1

N

and Tk+’ = 2 qlkgk-к for K<kSK + s,

і = 1

dk= 2 2^U for к > К + s. (9.23)

j-ii-i

In a similar way eqns. (9.10) and (9.11) are changed in order to take into account the errors in the measurements of T*.

After having calculated the age spectra B! for each channel і and each passage j, it is possible to calculate the total age spectrum A, at the entry of each channel,

J

A = 2 B,‘. (9.24)

i-i

These age spectra are then normalized so that

2 L

к

From the total age spectrum at the entry of each channel it is then possible to obtain the total age spectrum in every other point z,

alk+,= a,.k. (9.25)

In order to calculate s one needs an assumption on the flux distribution along the channel

<p(x) dx

———— 5,.

<p (x) dx

The discretization of the age in groups means that 5 will have to be rounded off to the nearest integer. Following this procedure the KUGEL code calculates a two­dimensional matrix of total age spectra. It is then easy to convert these total age spectra in isotopic compositions.

Let A* be the total age spectrum in a given point x. A previously performed burn-up calculation relates each age group к with a concentration vector Сь (This vector has as many components as isotopes and gives for each isotope the concentration correspond­ing to age group k.)

The composition in every point x is then given by

r=^CA’. (9.27)

к

As already pointed out, this treatment assumes that the composition of a fuel element is fully defined by its irradiation, independently of the actual form of the time dependence of the flux during burn-up. In reality, as the flux is not constant in the axial direction, the fuel elements are not burnt in a constant flux.

The error introduced by this simplification is very small except for those isotopes (like,35Xe and 233Pa) with short half-lives whose concentration is the result of a competition between flux dependent processes (production and absorption) and flux independent processes (decay). The concentration of 135Xe can be calculated at a later stage independently of the burn-up calculations (see Chapter 12).

Assuming that the axial fuel element velocity is radially constant, it is possible to perform an equilibrium burn-up calculation without the need of the above assumption.

In this case one can check the error introduced by this assumption, which is very small, even in the case of a greatly distorted axial flux distribution. Figure 9.7 gives the discharge spectrum for the 300 MWe THTR reactor calculated once without error in the burn-up measurement, and once with an error of ± 125 days.

A low axial fuel element velocity has the effect of shifting the fissile concentration, and hence the power, towards the top of the core. This can be seen in Fig. 9.8 where the THTR axial power distribution is compared with a hypothetical case with infinite velocity.

This fact can be used in order to approximate an exponential axial power. This same principle of separation of the burn-up calculation from the calculation of the fuel element distribution can be used for running-in calculations.<20> The fuel is again classified in age groups. The core is subdivided radially in channels and axially in zones. After each time step a zone is shifted into the next one. At defined time intervals a two-dimensional diffusion calculation is performed. The burn-up measurement and the

Full power days [t] Fig. 9.7. Discharge spectrum for the 300 MWe THTR.

Depth in the core [m] Fig. 9.8. Influence of the passage time of the fuel elements on the axial power distribution of pebble-bed reactor.

radial distribution of the recirculated elements is taken into account in the calculation with the same methods used for the equilibrium case. In the first core the power flatten­ing is obtained with a proper radial distribution of dummy elements and burnable poison.

These simplified burn-up methods cannot be applied in all cases. A very strong axial dependence of the power distribution influences the axial dependence of the neutron spectrum and of the cross-sections. An extreme case is given by the fuel cycles with only one passage through the core (OTTO) where the very strong axial dependence of burn-up is used to approximate an exponential axial power distribution with the aim of temperature flattening. In this case the simplification of performing a burn-up calcula­tion with constant flux and cross-sections, separate from the calculation of the fuel element distributions is no longer acceptable. On the other hand, the fact of having only one passage through the core simplifies the space-dependent burn-up calculations performed with standard methods. A high number of core regions can be used and the axial movement can be discretized in many axial shifts of one region into the next. An iteration is necessary since in order to obtain the axial flux distribution, which is needed for the burn-up calculation, one needs the axial dependence of the composition which is obtained from the burn-up calculation.