Reload Pattern

Annals of Nuclear Energy 28 (2001) 1377±1411 www.elsevier.com/locate/anucene

Development of parallellized higher-order generalized depletion perturbation theory for application in equilibrium cycle optimization R. van Geemert *, J.E. Hoogenboom Interfaculty Reactor Institute, Mekelweg 15, 2629 JB, Delft, The Netherlands Received 25 October 2000; accepted 20 November 2000

Abstract As nuclear fuel economy is basically a multi-cycle issue, a fair way of evaluating reload patterns is to consider their performance in the case of an equilibrium cycle. The equilibrium cycle associated with a reload pattern is de®ned as the limit fuel cycle that eventually emerges after multiple successive periodic refueling, each time implementing the same reload scheme. Since the equilibrium cycle is the solution of a reload operation invariance equation, it can in principle be found with sucient accuracy only by applying an iterative procedure, simulating the emergence of the limit cycle. For a design purpose such as the optimization of reload patterns, in which many di€erent equilibrium cycle perturbations (resulting from many di€erent limited changes in the reload operator) must be evaluated, this requires far too much computational e€ort. However, for very fast calculation of these many di€erent equilibrium cycle perturbations it is also possible to set up a generalized variational approach. This approach results in an iterative scheme that yields the exact perturbation in the equilibrium cycle solution as well, in an accelerated way. Furthermore, both the solution of the adjoint equations occurring in the perturbation theory formalism and the implementation of the optimization algorithm have been parallellized and executed on a massively parallel machine. The combination of parallellism and generalized perturbation theory o€ers the opportunity to perform very exhaustive, fast and accurate sampling of the solution space for the equilibrium cycle reload pattern optimization problem. # 2001 Elsevier Science Ltd. All rights reserved.

* Corresponding author at current working address: Paul Scherrer Institute, CH-5232 Villigen PSI, Switzerland. E-mail addresses: [email protected] (R. van Geemert), [email protected] (J.E. Hoogenboom). 0306-4549/01/$ - see front matter # 2001 Elsevier Science Ltd. All rights reserved. PII: S0306-4549(00)00135-3

1378

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1. Introduction In order to evaluate the `quality' of a reloading scheme, it is important not to study the associated fuel economy for only the forthcoming cycle. Instead, some method should be available to gain some indication about the scheme's multi-cycle performance. The most obvious and interesting multi-cycle evaluational method is to study the equilibrium cycle behaviour. The equilibrium cycle associated with a reload pattern is de®ned as the limit fuel cycle that eventually emerges after multiple successive periodic refueling with the same reload pattern. In reload pattern optimization, search algorithms are mostly based on the assessment of the e€ects of permutations in the fuel bundle shu‚ing scheme, followed by acceptance decisions in which the permutation yielding the largest improvement in the objective function value (while satisfying the operational and safety constraints) is chosen to determine the next reference pattern in the search procedure. Generally, the exact value of the objective function for an equilibrium cycle can only be determined by implementing a forward iterative procedure to obtain the equilibrium cycle solution from BOC to EOC. This implies that reload pattern optimization methods in which large numbers of di€erent refueling schemes must be evaluated in this way are quite expensive from a computational point of view. Hence it would be helpful to have accelerated iterative methods available for fast assessment of di€erent reloading schemes in such a way that the exact perturbed equilibrium cycle characteristics can be obtained with considerably reduced computational e€ort. A lot of work has already been done in the ®eld of applying perturbation theory to in-core fuel management optimization, but most of it is dedicated to non-equilibrium cycle fuel management (White and Avila, 1990; Maldonado et al., 1995; Moore and Turinsky, 1998; Van Geemert et al., 1998b). Yang and Downar (1989) did develop perturbation theoretical methods for the equilibrium cycle, in order to assess the e€ects of changes in general reactor physics parameters (such as material properties or the enrichment level of fresh fuel assemblies) on the required feed enrichment or on the cycle length for a constrained equilibrium cycle (with the reload operator ®xed). Constrained here means that either the feed enrichment or the cycle length are dictated by the condition that the e€ective multiplication factor ke€(uc)(EOC) of the uncontrolled core (without external reactivity control, such as control rods or soluble boron) at EOC is ®xed at a certain value (usually unity). In their paper, either the feed enrichment is ®xed, whereas the length of the equilibrium cycle is subject to this EOC condition, or vice versa. Unfortunately, Yang and Downar only reported results for zero-dimensional cases (so with a uniform ¯ux in the entire reactor core) and for a ®xed reload pattern. However, for application in loading pattern optimization, it is of course necessary to account for spatial variations in the ¯ux distribution and to consider changes in a far more abstract, discrete equilibrium cycle parameter, which is the reload pattern itself. In studying equilibrium cycle behaviour, one should realize that a change in the shuf¯ing scheme a€ects the entire BOC fuel density distribution. The generalized perturbation theoretical approach for the equilibrium cycle proposed by us will be shown to result in an accelerated iterative scheme that yields the exact perturbation in the equilibrium cycle solution (caused by a limited change in the reload operator)

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1379

in considerably fewer iterative steps and thus with a signi®cantly reduced computational e€ort. 2. General aspects of equilibrium cycle optimization As is generally known (Stevens, 1995, Van Geemert, 1999), no genuine ordering principle exists for the reload patterns or, in other words, it is not possible to map the solution space into regions (of comparable performances) which can be easily de®ned in terms of the variables of the reload operator. However, it seems that it is possible for the equilibrium cycle to establish some rough link between the shape of the power distribution and the equilibrium cycle fuel economy. Utilizing a simpli®ed but rather elegant model (De Jong, 1995), it can be shown that equilibrium cycles are generally more economical when they are characterized by power distribution shapes that are ¯at with regard to both space and time. Via Haling's theorem (stating that an operation cycle's total maximal power peaking factor during the cycle is minimal in the case of a constant power distribution shape throughout the cycle) this general correlation can be translated into the following rough property:loading schemes that are economical in terms of the equilibrium cycle generally give rise to rather constant and ¯at power distributions which thus yield low maximal power peaking factors. Obviously, knowledge of the existence of this general, rough correlation does not reduce the complexity of the loading pattern optimization problem in any way not only because it is only a rough correlation but also since it is impossible to establish an easy mathematical relationship between the way in which the core is con®gured and the resulting time-dependent power distributions. Two extreme examples of loading pattern types are the centre-to-outside loading (COL) and the outside-tocenter loading (OCL). In the case of COL loading, the fresh fuel assemblies are placed in the centre of the core, whereas the burnt fuel assemblies will be placed closer to the periphery as their burnup increases. In the case of OCL loading the reverse happens. Generally, in the equilibrium situation they will both give rise to non-optimal operation cycle behaviour with an additional unacceptably high power peaking factor in the COL case. In practise, power distributions that are ¯at and that change minimally throughout the operation cycle are to be achieved by, for example, checkerboard-like patterns with ring structures in the placement of the di€erent fuel ages. And hence in general, loading pattern design can be done quite satisfactorily by the application of rough engineering knowledge translated into some general rules of thumb for where and where not to place fuel assemblies with di€erent histories. But very obviously, a huge freedom still exists in the design of such patterns. Further optimization starting from a pattern obtained in this way can only be achieved by means of implementing automated optimization procedures featuring inherent mathematical, physical or even biological mechanisms as the driving force of their e€ectiveness. Well-known examples of these are cyclic interchange optimization methods (De Jong, 1995; Van Geemert, 1999), simulated annealing methods (Kropaczek, 1991; Parks, 1987; Verhagen, 1993; Smuc, 1994; Stevens, 1995; Van Geemert, 1998a) and evolutionary optimization algorithms (DeChaine 1996;

1380

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

Parks, 1996; Axmann, 1997; Van Geemert, 1999). Since in principle, the equilibrium cycle can only be determined exactly as the converged outcome of a computationally expensive calculational procedure, these optimization algorithms can be quite timeconsuming if no attempts are made to accelerate the iterative calculations. For equilibrium cycle optimization, it is worthwile to implement perturbation theory approaches at two di€erent levels: (i) Acceleration of the iterative solutions of the quasi-static neutronics equations by application of generalized perturbation theory (GPT). (ii) Acceleration of the perturbed limit cycle iterations by application of an extended time-dependent generalized variational approach, which we call generalized equilibrium cycle depletion perturbation theory (GECDPT). In this approach, the classical depletion perturbation theory (Williams, 1979) is embedded. The GPT approach associated with acceleration technique (i), which is not really speci®c for equilibrium cycle optimization purposes, has already been applied succesfully for non-equilibrium cycle optimization purposes (Maldonado et al., 1995; Moore and Turinsky, 1997, Van Geemert et al., 1998b). The most important line of research adopted in this work was to investigate whether a similar higher-order iterative method, to be referred to as GECDPT, could be set up for the equilibrium cycle. In the single time-step GPT, a generalized response functional is de®ned that consists of a response funtion and a summation of constraint functions formed by inner products of Langrange multipliers (later to be referred to as the adjoint ®elds) and the system equations that govern the time evolution of the reactor. This is done in such a way that the constraint functions vanish when the system equations are exactly satis®ed. GECDPT involves the set-up of an extended generalized response functional that includes the cyclic reload operation invariance equation as a constraint function in addition to the other, more traditional constraint functions (i.e. the ¯ux shape l-eigenvalue equation, the power normalisation equation, the depletion equations, etc.). As it turns out, an adjoint cyclic invariance equation for the adjoint time-dependent nuclide density ®eld surfaces that, if satis®ed, will give rise to the disappearance of the principal term in the development of the generalized functional change induced by a discrete change in the reload operator. This term will be shown to be dominant with respect to the higher-order terms for suciently modest changes in the reload operator, as a result of which the higher-order iterative scheme will quickly converge to the exact perturbed equilibrium cycle. The GECDPT approach follows from the analysis of an extended form of the GPT functional. Because of this, and because the single time-step GPT method has been applied in this study, we will now provide an introduction into the single time-step GPT method, so as to facilitate the introduction of the GECDPT formalism. 3. The single time-step GPT method The time evolution of a PWR core is basically described by the energy-, space-, and time-dependent neutron ¯ux …t† and by the space- and time-dependent nuclide

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1381

density matrix N…t†. Generally, since the neutronic and depletion equations are coupled, a quasi-static approach, featuring a number of time steps ti, must be adopted for numerical simulation of how the core evolves from BOC to EOC. In this approach, it is assumed that the neutronics eigenvalue equation, determining the spatial ¯ux shape i , can be written compactly as …Li

li Fi †

ˆ0

i

…1†

where the eigenvalue l is the reciprocal of the e€ective multiplication factor of the uncontrolled (=without external reactivity control) core ke€ and L and F are the loss and production operator, respectively. The normalisation requirement for the spatial ¯ux shape vector can be written as X …2† I ˆ1 all nodes I

with the index I denoting the di€erent core nodes in the system geometry. In the 1 12group approximation adopted in this study, there is e€ectively only one energy group, but the general GPT and GECDPT formulations presented here can be easily extended for incorporation of multiple energy groups. When determining the response functional, one should explicitly account for the fact that a change in the nuclide density ®eld also in¯uences the neutron ¯ux ®eld. In conformity with the variational approach, this can be realized by treating the neutronics and ¯ux normalisation equations as constraints on the response quantity itself, and appending them to the response function using Lagrange multipliers. Many inner product de®nitions occur, which can be most conveniently written using the Dirac bracket notations:

 X ajb ˆ aI bI ; …3† I

For the eigenvector and the eigenvalue l, de®ned at a certain time step ti, two coupled functionals can be de®ned, in which the quasi-static equations are appended to the response functions using the Lagrange multipliers  , a and  . For the sake  of simplicity, we now abandon the time index (i) 8 > > > R1 ˆ > < > > R2 ˆ l ‡ > > :

D D

 

j …L

lF†

j …L lF† D E  jF

E E

D E a 1j

1

 …4†

We note that, due to the commutativity property of the adjoint operators, we can write: D E D E  j …L lF† ˆ …L lF †  j …5†

1382

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

with L and F the adjoint operators of L and F (occurring in the de®nition of R2,

  respectively. The same applies to the inner product j … L lF † , which can be

 written as …L lF †  j . If l and are exact solutions to the quasi-static equations, then R1 ˆ and R2 ˆ l irrespective of the choices for  , a and  . If the nuclide density distribution N is perturbed ( ˆ N ! N0 ˆ N ‡ N), R1 and R2 will be perturbed as well, since the change in N e€ects a change in the operators L and F:
R1;2 ! R01;2 L0 ; F0 ; 0 ; l0 …6† where the prime variables refer to their perturbed values. Again, if l0 and 0 are exact solutions to the perturbed equations, then R01 ˆ 0 and R02 ˆ l0 . In single-time step GPT, the Lagrange multipliers  , a and  are numerically conditioned in such a way that the following iterative system, 8 D E D E 
…m‡1†  > > ˆ …m † j L0 l…m† F0 …m† a 1 j …m † 1 > > D E <
 j L0 l…m† F0 …m† …7† ; > D E > l…m‡1† ˆ l…m† ‡ > >  j F0 …m† : if convergent (depending on how much L0 and F0 di€er from L and F), converges to the exact perturbed solution (l0 ; 0 ) in a number of steps that becomes very small for small changes L and F. The conditioning of  , a and  follows from a number of Euler stationarity conditions which for the speci®c ®rst-order components in the development for R1 and R2 that contain perturbations in forward quantities that are still to be determined (l and  ), to vanish. Substracting the expressions for the perturbed and the unperturbed response functional and ordering the terms gives the developments D E D E   …m‡1† ˆ  …m† j …L lF† ‡ l…m†  j F D E …L lF †  ‡ 1a j  …m† …1st order† D E D E …8†  j …L lF† …m† ‡ l…m†  j F …2nd order† D E ‡ l…m†  j F …m† …3rd order† and l…m‡1† ˆ

D …L 

lF †



j

…m †

D



E

‡

j F0

D

 …m†

j …L E

lF†

Noting that …m† ˆ  …m† can be written as 1 j  Eq. (8) can be written for this case as:

…m†

E …9†

…m†



with 1 the unity matrix,

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411



…m‡1†

D

ˆ



j …L

lF†

E

E jF E 1 j  …m† …1st order†

‡ l…m†

D

1383



D L lF  ‡ 1a D E  j …L lF† …2nd order† D E ‡ l…m†  j F …m† …3rd order†

…10†

The third and the second term on the right will vanish if the adjoint ¯ux shape equation …L

lF †



ˆ1

1 a ˆ Q

…11†

and the orthogonality condition D E  jF ˆ0

…12†

are satis®ed. For argueing how nthis combined objective o can be realized, it is worthwile to consider the eigenset ln ; n ; n ˆ 0; 1; 2; . . . of the adjoint ¯ux shape l-eigenvalue equation:
…13† L ln F n ˆ 0 n which is, of course, fully analogous to the standard eigenset ln ; the normal ¯ux shape l-eigenvalue equation: …L

ln F †

n

 n; n

o ˆ 0; 1; 2; . . . , of

ˆ0

…14†

Generally, ln ˆ ln . The fundamental adjoint ¯ux shape 0 (conventionally written simply as  , like 0 is written as ) is a very important quantity. The basic reason for this is the validity of a mathematical orthogonality property that can be proved in the following way: after premultiplication of Eq. (14) (and abandoning the   time index

 i)with the lth eigenvector l , recalling the commutativity property l j L n ˆ L l j n and employing Eq. (13), we obtain D E D E …15† ll l j F n ˆ ln l j F n Hence, since ln 6ˆ ll in the case n 6ˆ l, the inner product for the case n 6ˆ l. Thus, the validity of the property D E  j F ˆ cn nl n l



 l

jF

n



should vanish …16†

has been established. The implication of this property is that, if some adjoint ¯ux shape  is to be orthogonal to the ®ssion rate distribution F , it should not contain any component of  i.e. it should be orthogonal to  . So, Eq. (12) basically

1384

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

requires that the Lagrange multiplier  is not allowed to contain any component of the fundamental solution of the homogeneous adjoint equation …L lF †  ˆ 0. This fundamental mode is the only eigenstate which does not yield zero when formed the inner product with F , since the eigensets and n ; n ˆ 0;  1; 2; . . . and n ; n ˆ 0; 1; 2; . . . obey the general orthogonality property l j F n ˆ cn nl . Thus, n j F unequals zero only in the case n=0. More if p is the particular solution to

 speci®cally,  Eq. (11) satisfying Eq. (12) such that p j F vanishes, where is the fundamental solution to the homogeneous equation, then p ‡ b  is also a solution of Eq. (11) for all b. However, due to Eq. (12) the value of b is ®xed to be zero, so that  will not be `contaminated' with the fundamental solution  . In each step in the iterative process of solving Eq. (11), this can be e€ected by application of the ®ltering operation D E  jF  E  :ˆ  D …17†  jF The choice for the multiplier a is determined by the fact that Eq. (14) and Eq. (11) specify j Q ˆ 0. This becomes obvious when analysing the inner product

 that  jQ : D E D E j Q ˆ j L lF j  ˆ …18† D

…L

lF† j



E

D ˆ 0j

Hence, a is given by D E 1j Eˆ aˆD 1j



E

ˆ0

…19†

…20†

because of the normalisation condition (2) on . If Eqs. (11), (12) and (20) are satis®ed by  and a, the expression for the ®rst-order prediction of the change in the response ¯ux  , denoted by …1† , is D E  …1† ˆ j …L lF† …21† Needless to say, this expression yields an accurate result only in case of very small changes L and F. Regarding Eq. (9) since …L lF †  ˆ 0, the ®rst term in the nominator of (9) vanishes, and thus we obtain: D E  j …L lF† …m† D E …m‡1† l ˆ …22†  j F0 …m† With the nodal ¯ux distribution as the response vector for functional R1 and the l-eigenvalue as the scalar response quantity for functional R2, these expressions

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1385

eventually yield a higher-order iterative scheme from which the exact perturbed nodal ¯ux distribution can be obtained at low computational cost: 8 D E D E
 > >  …m‡1† ˆ …1† j L …m† …lF†  …m† ‡ l  j F …m† > > D E <  j …L lF† …m† …23† > D E > …m‡1† l ˆ > >  j F0 …m† : This scheme has the convenient property that it will converge very rapidly when the perturbation operators L and F are relatively small with respect to L and F. This can be made clear by pointing out that from Eq. (11) it follows that





 …L



lF † 1 ˆ



…L

1

lF †

so that the condition number for the scheme (23) becomes



…L lF †



…L lF †

…24†

…25†

Hence, the smaller L and F, the faster the convergence. This will generally be the case for an appropriately conditioned numerical design strategy such as reload pattern design optimization, in which the results of many di€erent small changes in the system need to be evaluated during the optimization procedure. We stress that only when many of such changes are to be evaluated, the evaluational economy will outweigh the extra computational e€ort involved in determining the adjoint ®elds. It should be pointed out here that scheme (23) results from the exact response functional expansion by explicit implementation of the supposedly exact validity of the adjoint Eqs. (11) and (20). Hence if this validity is limited because of a limited number of iterative steps available for solving Eqs. (11) and (20), small errors can be expected. This problem vanishes when one applies the iterative scheme (7) that is based on the original exact response functional expansion, without ommitting any terms. Since the same adjoint ®elds  and a are used as in scheme (23), the terms that are neglected in scheme (7) will become very small all right, such that the same convergence behaviour will emerge. Thus, scheme (7) is the scheme to implement, and the signi®cance of scheme Eq. (23) reduces to a mere theoretical one, since it re¯ects the inner dynamics of scheme (7) explicitly, and clari®es why fast convergence occurs for small changes in L and F caused by a limited change X in the reload operator. 4. Generalized equilibrium cycle depletion perturbation theory Since the time dependent neutronics and nuclide density ®elds implicitly depend on one another, a spatial perturbation in the nuclide density ®eld at BOC (so at

1386

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

t ˆ t‡ 0 ) will perturb the entire neutron/nuclide ®eld for the forthcoming cycle. In classical depletion perturbation theory (DPT) (Gandini, 1975; Williams, 1979), it is conventionally assumed that the perturbed BOC nuclide density distribution is known exactly. Due to the fact that the neutronics ®eld and the nuclide density ®eld are governed by a set of coupled (di€erential) equations and cannot be varied independently, it is impossible to calculate directly the exact way in which a perturbation in the combined neutron/nuclide ®eld will propagate from BOC to EOC. Any data change that changes one ®eld will also change the other ®eld, because the two ®elds are constrained to ``move'' only in a fashion consistent with their coupled ®eld equations. This will be even more complicated when one is dealing with an equilibrium cycle, in which case one wants to evaluate the e€ect of a change in the reload operator on the equilibrium cycle solution. The reload operator can be represented mathematically by a square binary matrix X, the elements of which are to be de®ned by

XIJ ˆ

8 <1 :

0

if the fuel element that has resided in node J is located in node I after reloading otherwise

…26†

P If I XIJ ˆ 0 Pthen obviously the fuel element that has resided in node J is to be discharged. If J XIJ ˆ 0 then apparently none of the older fuel elements which were already present in the core will be place in node I during reloading. In that case, node I is to receive a fresh fuel element. If the fuel elements are characterized by their nuclide density vectors, the invariance of the equilibrium cycle with respect to the reload operation can be mathematically de®ned as:
-
N t‡ 0 ˆ XN tT ‡ …1

X†NF

…27†

NF is a convenient notation for …NF ; NF ; . . . ; NF †T with NF the vector denoting the densities of the di€erent nuclides in the fresh fuel. The symbol 1 indicates the unity operator. The equilibrium cycle can be determined iteratively using this invariance condition. Starting from an initial guess for the equilibrium cycle, successive cycles are simulated, with at each reloading the application of the same reload scheme and the same composition of the fresh fuel assemblies, until a converged limit cycle emerges. For compact representation and notation of equibrium cycle reload operators, one can apply the fact that reload schemes for N-node systems, featuring n age groups for the fuel elements (batches), lead to N/n fuel bundle life history trajectories. These trajectories represent the sequences of fuel positions that act as hosts for a particular fuel element during its time in the reactor core (normally 4 cycles, so n=4), as it travels from one position to the next until it is discharged after n cycles. The collection of all such history trajectories is a convenient and compact representation of the associated reload pattern, and can be written as

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

8 …1 † i1 > > > > ... > > < ... Hˆ >... > > > > ... > : …Nn† i1

!

!

i…21† ... ... ... . .N. … † i2 n

!

!

... ... ... ... ... ...

!

!

9 i…n1† > > > ... > > > = ... ... > > > > . .N. > > … †; in n

1387

…28†

In this denotation, in which for example the ®rst row of H denotes the ®rst fuel element life history trajectory, each arrow represents an individual fuel element repositioning operation. The nonzero elements in the matrix X follow from Xi… j† ;i… j† ˆ 1 for all m ˆ 1; . . . n 1; j ˆ 1; . . . ; m

m‡1

N : n

…29†

In order to account for the e€ect on the equilibrium cycle of a change in the reload operator, an extra term representing the reload operation invariance condition needs to be added to the response functional. Needless to say, this will make the determination of the properly chosen adjoint nuclide density ®eld somewhat more complicated. Naturally, in the more general time-dependent picture a number of additional constraint functions must be included in the functional, corresponding to the depletion processes and the refueling operation which occurs periodically in time. Again, all the system equations can be treated as constraints on the response, and appended to the response function using appropriate Lagrange multipliers. Since the reload operation takes place in the nuclide density space, the reload equation will have the same Lagrange multiplier as the nuclide transmutation equation. First of all, we note that the transmutation equation dictating the time evolution of the nuclide density ®eld can be written as   @ ^ N…t† N…t† ˆ ^…t† ‡ D @t

…30†

with ^ denoting the absorption cross sections matrix and D^ the decay operator. We emphasize here that the operators active in the nuclide type space are marked with a hat, and that an implicit product in the nuclide type space occurs when the operators ^ and D^ appear in combination with the nuclide density ®eld N…t†. The fast ¯ux …t† can be written as a product of the time dependent shape function …r; t† (which is „ normalized such that …r; t†dV ˆ 1 for all t) and a time dependent power normalisation factor …t†. This allows for a convenient implementation of the condition that the reactor power is restricted to remain at a constant level. Using the inner product notation, this condition can be written as: D E wf i F;i j i ˆ Ptotal for all t …31†

1388

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

The response quantity we are interested in is the full collection of all nuclide densities at the end of the operation cycle, so de®ned at t ˆ tT- . In order
to stress that the response is de®ned at t ˆ tT- , we denote the response quantity as N t-T . Applying the Lagrange multipliers Ci , Pi and N …t†, the general response functional R X; N…t†; NF ; i ; i ; li can be written as: RN tT




T 1 … ti‡1  X N …t† j ^i ‡ iˆ0 ti

T 1D X iˆ0

‡

 i

j …Li

… t‡  T N …t† j tT

li Fi †

 @ j N…t† dt @t

^ i‡D E i

T 1  D X Pi wf i F;i j iˆ0

@ N…t† ‡ …X @t

 1† N…t†

E i

Ptotal

 NF …t

 tT † dt



…32†

in which the integrand in the last term is the di€erential equation describing the reload operation at EOC. Forward

integration of this equation from tT yields the reload equation N t‡ t :ˆ XN X†NF . Now, for the sake of convenience, T ‡ …1 0 we think of R as being composed of two parts R1 and R2: R=R1+R2. R1 is associated with the `in-cycle' system equations (describing the `in-cycle' neutronics, power control and burnup) and R2 is associated with the discrete time event of reloading (with respect to which the equilibrium cycle should be invariant): R1 ˆ

T 1 … ti‡1  X iˆ0 ti T 1 X iˆ0

N …t† j ^i

^ i‡D

 D Pi wf i F;i j

E



i

 @ j N…t† dt @t Ptotal



T 1D X C iˆ0

 i

j …Li

li Fi †

E i

…33†

and
R 2 ˆ N tT ‡

… t‡  T N …t† j tT

@ N…t† ‡ …X @t


1† N tT

 NF …t

 tT † dt

…34†

If N…t†, i , li and i are exact
equilibrium cycle solutions to the quasi-static burnup equations, then R ˆ N tT . Now, if the reload operator X and the fresh fuel composition NF (and hence the entire equilibrium cycle solution) are perturbed, this will in¯uence R:   R ! R0 X0 ; N0 …t†; N0F ; 0i ; 0i ; l0i ; …35† where the prime variables refer to their perturbed values. Again, if N0 …t†, 0i , 0i and l0i exactly satisfy the perturbed equilibrium cycle system equations dictated by X0
0 0 0 and NF , then R ˆ N tT , so the perturbed functional value will be equal to the exact perturbed EOC equilibrium cycle nuclide density distribution. Expanding R

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1389

about the unperturbed state and neglecting the second-order terms in the continuous variables, we obtain the ®rst order prediction for R, denoted as …1† R: * + * +! T 1 … ti‡1 X @R1 @R @R1 @R1 …1† j N…t† dt ‡ li ‡ i ‡ j i  Rˆ @li @N…t† @i @ i ti iˆ0 h i ‡ R2 X; NF ; N…t† …36† with h i
R2 X; NF ; N…t† ˆ N tT … t‡  h i   T @ N …t† j N…t† ‡ XN…t† ‡ X N…t† NF ‡ XN…t† …t ‡ @t tT E ‡ XN…t†…t tT † ‡ …1 X0 †NF …t tT † dt

tT †

…37†

h i We note that R2 X; NF ; N…t† as de®ned in Eq. (37) is not a ®rst order approximative expression, but an exact expansion of the change of R2 due to the simultaneous changes X, NF and N…t†. This is why no superscript is added. Many of the contributions to …1† R can be forced to vanish by de®ning the Euler±Lagrange stationarity equations for the adjoint ®elds C i , Pi and N …t†. For example, the functional derivative with respect to i is: @R1 ˆ @i

… ti‡1 D ti

N …t† j ^

i

E j N…t† dt

D Pi wf F;i j

E

…38†

i

For this expression to vanish, the multiplier Pi should be chosen such that: Pi ˆ

1 D wf F;i j

E

… ti‡1 D N …t†j^

i

ti

E

i jN…t†

dt

…39†

Employing the commutativity property of the adjoint operators L* and F*, the functional derivative with respect to i can be written as @R1 ˆ @ i

… ti‡1 D ti

E N …t†j^i jN…t† dt

Li


li Fi C

 i

wf Pi i F;i

In analogy with Eq. (11), this partial derivative will vanish if C inhomogeneous adjoint ¯ux shape equation
Li li Fi C i ˆ Qi with the adjoint source

…40†  i

satis®es the …41†

1390

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

Qi ˆ

… ti‡1 D ti

E N …t†j^i jN…t† dt

wf Pi i F;i

…42†

1 The requirement that @R @li vanishes results in the same orthogonality condition as  Eq. (12), forcing C i to be orthogonal to the neutron production term Fi i at t ˆ ti : D E …43† C i j Fi i ˆ 0

The stationarity condition corresponding to a variation in N…t† is more complex than for the other variables. This is basically because N…t† is the only `forward' quantity that is allowed to vary continuously in time. Due to this, a signi®cant part of its contribution in the variation functional is embedded in an integral, while the other part merely consists of terms de®ned at discrete equidistant moments ti containing the `snapshots' N…ti † of the nuclide density ®elds. To derive the Euler conditions for N …t†, we ®rst consider the variation of R1 with respect to N…t†: T 1 … ti‡1  X  R1 N…t† ˆ N …t†j^i …1 †

h

 @ jN…t† dt @t + * @Fi @F;i  j i ‡ Pi wf i li j @Ni @Ni

i

iˆ0 ti

T 1 X

"* C

iˆ0

 ij

@Li @Ni

^ i‡D

+# i

…44† :Ni

Using the adjoint operators ^  and D^  and partial integration with respect to t yields  … ti‡1  D

E @  ‡ ^ jN…t† dt ˆ N t‡ N …t†j^i i ‡ D t j N i i @t ti …45†   … ti‡1  D E

^  ‡ @ N …t† j N…t† dt N ti‡1 j N ti‡1 ‡ ^  i i ‡ D @t ti If we de®ne * zi ˆ C

 ij

@Li @Ni

@Fi j li @Ni

+ i

‡

Pi wf i

* + @F;i j i ; @Ni

…46†

  @ N …t†jN…t† dt @t D

Ei N ti‡1 jN ti‡1

…47†

h i we can write …1† R1 N…t† as T 1… ti‡1  h i X ^  i R1 N…t† ˆ iˆ0

D

‡ N



ti

t‡ i




zi jN

^ i‡D ‡ t‡ i


E

If we have the N …t† satisfy the in-cycle time boundary discontinuity conditions

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

N ti




ˆ N t‡ i




zi ; i ˆ 0; 1; . . . ; T

2;

the summation will reduce to h i D

E D 

E …1† R1 N…t† ˆ N t‡ z0 jN t‡ N tT jN tT 0 0   T‡1 … ti‡1  X ^  ‡ @ N …t†jN…t† dt ^  i i ‡ D ‡ @t iˆ0 ti

1391

…48†

…49†

Further, if the adjoint transmutation equation, @  N …t† ˆ @t

 ^  i

i

 ^  N …t†; ti < t < ti‡1 ‡D

is satis®ed, the summation eventually reduces to h i D

E D 

E ‡ jN …1† R1 N…t† ˆ N t‡ t t t z N j N 0 T T 0 0

…50†

…51†

 Before we consider R2 N t†Š, impose the stationarity condition for this case and derive the adjoint cyclic boundary condition, we point out some continuity and discontinuity properties and conventions for the forward and adjoint nuclide density ®elds. For the `in-cycle' intervals characterized by ti < t < ti‡1 , the forward timedependent nuclide density distribution is absolutely continuous in the sense that

N ti ˆ N t‡ . The only discontinuity for the forward ®eld occurs at t ˆ tT , where i the reload operation is applied, resulting in the reloaded state at t ˆ t0 . The adjoint time-dependent nuclide density distribution however has a discontinuity between each in-cycle interval and the next as well, in accordance with (48). We stress that, strictly mathematically, t ˆ tT and t ˆ t0 coincide due to which, in the framework of the formalism presented here, t ˆ tT and t ˆ t0 are to be thought of as t0 ˆ t ˆ t‡ T or t ˆ tT ˆ t0 . Now, for the adjoint nuclide density distribution we adopt the convention
that N …ti †  N t‡ i
. Likewise,
at the cyclic boundary (de®ned equivalently by t ˆ t0  ‡ and t ˆ tT ), N t‡ 0 
N t0 . For the forward nuclide density distribution we implement N …tT †  N tT
at the cyclic boundary, which is of course adjoint to the de®nition N …t0 †  N t‡ for the adjoint ®eld and in analogy with the opposed time 0 directions of forward and backward depletion. These conventions are inspired by the fact that forward/backward depletion during a time interval yields certain values for the ®nal/initial nuclide densities for that interval. An important consequence of these conventions occurs in time integrals over the cyclic bounday, involving inner products of adjoint and forward ®elds. For example: … t‡ D T

tT

N …t†j…t

E D

E tT †N…t† dt ˆ N t‡ 0 jN tT

…52†

1392

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

Utilizing these properties, we will now consider R2 from Eq. (35) and derive the adjoint cyclic boundary condition corresponding to ®rst-order stationarity of R2 with respect to a perturbation in the nuclide density distribution. Realizing that
N tT can be written as N tT




… t‡ D T 1…t ˆ N…tT † ˆ tT

E tT †jN…t† dt;

…53†

applying partial integration with respect to t again and using the adjoint reload matrix X , R2 can be written as h i D

E D 

E ‡ jN jN t t t R2 X; NF ; N…t† ˆ N t‡ z ‡ N 0 0 0 T T   … t‡   T @ ‡ …X 1†…t tT † N …t† ‡ 1…t tT †jN…t† dt ‡ …54† @t tT … t‡ D E n   o T N …t†j X N0 …t† N0 F ‡ …1 X†NF …t tT † dt ‡ tT

We note that, in correspondence with Eq. (52), … t‡ D T

tT

n  N …t†j X N0 …t†

 N0 F ‡ …1

D
 0
ˆ N t‡ 0 jX N tT

E tT † dt

o X†NF …t

 E N0 F ‡ …1 ‡ X†NF ;

…55†



and recall that, for the equilibrium case, N t‡ 0 is related to X, NF and N tT as    


0 N t N N t …56† ˆ X N ‡ X ‡ NF N t‡ F F 0 T T Implementing this, and realizing that the ®rst term in Eq. (54) can be placed inside the integral constituting the third term, we obtain D   E

0 0 i j X N t N ‡ … 1 X †N …57† …1† R ˆ N t‡ z F 0 F 0 T … t‡   T @N …t† ‡ ……X ‡ @t tT



1†N …t†

z0 ‡ 1†…t

 tT †jN tT dt


…58†

Now, the integral in Eq. (58) will vanish if its integrand vanishes: @  N …t† ˆ …X …N …t† @t

z0 †

N …t† ‡ 1†…t

tT †

…59†

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

Integration of this di€erential equation from tT to t‡ T yields




  ‡ N t‡ z0 N t‡ N tT T ˆ X N tT T ‡1

1393

…60†

‡ For the equilibrium h cycle, i t‡ T is equivalent to t0 . Thus, the condition that the integral part of R2 N…t† disappears results in the following adjoint cyclic boundary condition to be satis®ed by the adjoint time-dependent nuclide density distribution N …t†:


z0 ‡ 1 …61† N tT ˆ X N t‡ 0

If this condition is satis®ed, …1† R can be written very compactly as D   E D

E …1† R ˆ S0 j X N tT N0 F ‡ …1 X†NF ‡ S0 j XN tT with the sensitivity operator S0 de®ned as
z0 S0 ˆ N t‡ 0

…62†

…63†

We stress that the principal part of this ®rst-order prediction …1† R of the perturbed equilibrium cycle EOC ®eld, D   E
S0 j X^ N tT N0 F ‡ …1 X†NF …64† can be obtained without any new iterative calculations. Further, it is interesting to point out the physical interpretation of Eq. (61). This physical interpretation is set in a reversed time picture in which importances are propagated from one fuel element position preceding the other in correspondence with an adjoint reload operation. This interpretation is illustrated in Fig. 1. Whereas the normal reload scheme indicates in a forward way how disturbances in fresh fuel elements travel through the core via the periodic shu‚ing process, the adjoint reload scheme indicates (in a backward way) the history of an observed disturbance in a discharged fuel element, thus eventually pinpointing from where and when the disturbance originates. Mathematically this implies that the adjoint trajectory representation H* can be obtained from H simply by reversing the orientation of the arrows in expression (28), after which the non-zero elements in X can be obtained with formula (29). 5. Accelerating the limit cycle iterations The general numerical goal de®ned for this study is to develop and implement an accelerated way to perform the equilibrium cycle iteration. In the normal picture this iteration is done by plain simulation of how a limit cycle is reached in real life, i.e. one simply performs the neutronics and depletion calculations during each mth iterated operation cycle, followed by a reloading, each time with the same reload pattern

1394

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

Fig. 1. Forward and adjoint reload operations.

of which one wants to obtain the associated equilibrium cycle. This process ends when the limit cycle has been obtained with sucient accuracy. In the numerical picture proposed here, we assume that one forward equilibrium cycle solution has been obtained for a speci®c reference reload pattern. The reference reload operator and its associated equilibrium cycle solution are assumed to have been adopted as the case for which the adjoint Euler stationarity Eqs. (41) and (50) have been solved. Further, it is assumed that the adjoint ®elds C  and N …t† have been solved such that the additional conditions (43), (48) and (61) are satis®ed. Now, all the new reload operators of which the associated equilibrium cycle solutions need to be evaluated are treated as perturbations on the reference reload operator. In the way to perform the perturbed equilibrium cycle iteration proposed here, a signi®cant calculational acceleration is realized in three ways: . The iteration is started with the estimate (64) as the 0th iterand, which is the ®rst-order prediction of how the equilibrium cycle EOC nuclide density distribution is perturbed. Via Eq. (56), this can be translated into a prediction of how the equilibrium cycle BOC nuclide density distribution is perturbed. . The individual `in-cycle' neutronics equations, for determination of the eigenvalues li and the ¯ux distributions i , i ˆ 0; . . . ; T 1 are performed by solving the iterative scheme (7) instead of applying the power method to Eq. (1). This yields a signi®cant reduction in the number of iterative steps required to obtain a certain prespeci®ed accuracy level. Since the depletion calculations are computationally cheap, these are executed in a normal way ; in this manner, the perturbed time-dependent nuclide densities N0 m …t1 †; N0 m …t2 †; . . . ; N0 m …tT † corresponding with the mth equilibrium cycle iterative step, are obtained. . After the (T-1)th neutronics calculation has been completed, and N0 m …tT † has  been calculated, a correction formula, in which
N …t† is used, is applied in …m‡1† order to obtain the corrected iterand corr N tT , before the normal numerical reload operation is performed.

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1395

We will now derive the accelerative expression to be used. First of all, we recall the general shape of the response functional,  T 1 … ti‡1 
X @ ^ jN…t† dt R ˆ N tT ‡ N …t†j^i i ‡ D @t iˆ0 ti T 1D X iˆ0

C i jLi

li F i j

E

i

T 1  D E X Pi wf i F;i j i iˆ0

… t‡  T ‡ N …t†j

@ N…t† ‡ …X @t

tT

 1† N…t†

Ptotal

 NF …t



 tT † dt

…65†

such that the accelerative expression to be implemented is: m† N …corr




tT ˆ 

…m†




N tT ‡

T 1D X iˆ0

T 1 … ti‡1  X iˆ0 ti

h C i j…m† …Li

 … t‡  T ‡ N …t†j tT



…m†

N …t†j

iˆ0

1†…t

tT † ‡ X…m† N…t††…t

…t

^ i‡D

  @ N…t† dt @t

T 1 iEX h D Ei  …m† j P w    i i F;i i i f

li Fi †

@ ‡ …X @t

 ^i

 
tT † …m† N…t† ‡ X N tT

tT † ‡ …1

X0 †NF …t

E tT † dt

NF

 …66†

Since in each mth iterative step the exact li, i are obtained from the `in-cycle' neutronics eigenvalue equations [in an accelerated way through solution of system (7)], all the individual terms occurring in the summations from i=0 to i=T 1 will vanish:    8 @ > …m † ^ >  ^i i ‡ D N…t† ˆ 0; > > < h i @t …67† …m† …Li li Fi † ;i ; ˆ 0 > > h D E i > > : …m†   j ˆ0 i

F;i

i

Further, recalling Eqs. (54) and (55), the expression reduces to

D

E D  ‡
…m†
E m† …corr N tT ˆ …m† N tT ‡ N tT j…m† N tT N t0 j N t‡ 0 … t‡ 

  @   …m † ‡ …X 1†…t tT † N …t†j N…t† dt ‡ @t tT D E D



E …m† ‡ N t‡ N tT ‡ N t‡ 0 jX N tT ‡ NF 0 jX D E
‡ N t‡ X0 †NF 0 j…1 T

…68†

1396

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

We now take into account that, if Eq. (59) is satis®ed, we have   @ ‡ …X 1†…t tT †† N …t† ‡ 1…t tT † ˆ X z0 …t tT †; @t due to which we obtain
D

E D  ‡
…m†
E D 
E m† …corr N tT ˆ N tT j…m† N tT N t0 j N t‡ ‡ X z0 j…m† N tT 0 D E D



E …m† ‡ N t‡ N tr ‡ N t‡ 0 jX N tT ‡ NF 0 jX D E
0 ‡ N t‡ … †N j 1 X F 0

…69†

…70†

Now, we
wish to end up with an accelerative equation in terms of only the iterand …m† N tT . For achieving this, we recall Eq. (56) and obtain
D


…m†
E m† …corr N tT ˆ N  tT N t X N t‡ z j  0 T 0 D 


E 0 …m† ‡ N t‡ N tT  …m 1 † N t T …71† 0 jX 


If the adjoint cyclic boundary condition N tT X N t‡ z0 ˆ 1 is satis®ed, 0 the scheme reduces to 

D


E m† 0 …m 1 † …m 1 † N tT ˆ …m† N tT ‡ N t‡  N t N t j X  …72† …corr T T 0 As argued before, we should realize that scheme (72) is exactly valid only if the adjoint cyclic boundary condition expressed by Eqs. (59) and (61) is satis®ed exactly. If for any reason one does not want to invest that much time in solving Eq. (61) with very high accuracy, one should simply implement 

D


E m† 0 …m† …m 1 † …corr N tT ˆ …m† N tT ‡ N t‡  N t N t j X  T T 0 D E
‡ " j …m† N tT …73† with " the residual operator
that
is zero only when the adjoint cyclic boundary
condition N tT X N t‡ z0 ˆ 1 is satis®ed exactly: 0


X N t‡ z0 1 …74† " ˆ N  tT 0 As long as " is small enough, this accelerative scheme should, when applied before each standard reload operation, result in an accelerated convergence towards the

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1397

exact
perturbation in the equilibrium cycle EOC nuclide density distribution, N tT . We would like to point out that, in principle, schemes (7) and (72) independently contribute to the acceleration of the total numerical strategy for obtaining the perturbed constrained equilibrium cycle solutions as functions of the reload operator and feed enrichment perturbations X ! X ‡ X and NF ! NF ‡ NF . So if scheme (7) an acceleration a1 and scheme (72) an acceleration a2, the total acceleration e€ect is quanti®ed by the product a ˆ a1 a2 . Further, we point out that in our numerical study we did not implement any additional classical numerical acceleration strategy, such as for example Gauss-Seidel (GS) acceleration. If one would do so, such an additional numerical acceleration strategy should simply add an extra factor to the formula for a, for example a ˆ aGS a1 a2 . For the GPT-based schemes, the more favourable convergence properties arise from their inherently more favourable mathematical structure in the case of modest perturbations. It is interesting to point out that the dependences of a1 and a2 on dX are quite di€erent. From scheme (23) it is quite obvious that the smaller X (and thus the Li and Fi ), the faster the convergence of scheme (7). Now, in order to gain more knowledge about the in¯uence of X on the equilibrium cycle convergence behaviour of the GECDPT (level 2) acceleration, we recall Eq. (56), which expresses the perturbed reload operation in terms of the nuclide density perturbations for the (m+1)th and the mth iterand,    


0 …m † …m‡1† N t‡ N t N  N t …75† ˆ X N ‡ X ‡ NF ; F F 0 T T 

N tT ˆ X …m† N t‡ 0

  NF ‡ X0 …m

1†

N tT




 NF ‡ NF

…76†

and subtract the latter from the former to obtain

0 …2†…m† …2†…m‡1† N t‡ N tT 0 ˆX

…77†


with @…2†…m‡1† N t‡ 0 de®ned as

…m‡1† …2†…m‡1† N t‡ N t‡ 0 ˆ 0

…78†

…m† N t‡ 0





During the depletion that follows the reload operation, the …2†…m‡1† N t‡ 0 is even
tually translated into …2†…m‡1† N tT via numerical solution of the burnup equations. Obviously, since a depletion operator will always have a norm which is smaller than

unity, …2†…m‡1† N tT will have a smaller norm than …2†…m‡1† N t‡ 0 :




…2†…m‡1† N tT



…79†


ˆ adepletion < 1

…2†…m† N t‡

 0 We note that, using de®nition (78) and the perturbed adjoint reload operator X0  , the level 2 accelerative scheme can be rewritten as

1398

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411


D 
…2†…m†
E 2†…m† …corr N tT ˆ X0 N t‡ N t j  0 T ; which implies that




…2†…m‡1†  N t

0   ‡
corr T 1

ˆ
ˆ X N t0 ‡ a2

…2†…m† N t0

…80†

…81†



Realizing that X0  N t‡ into its ®xed component X N t‡ 0 can be decomposed 0
and its variable component X N t‡ 0 , it becomes clear that the dependence of a2 on X is less pronounced than the dependence of a1 on X and of a somewhat more complex nature. From the numerical results that are presented at the end of this paper it can be concluded that, for most (magnitudes of) choices for X, a12 remains signi®cantly below unity (thus realizing and that the maximum speed-up

a speed-up)


occurs for X ˆ 0. (which is inversely proportional to X0  N t‡ 0 6. Parallel numerical solution of the adjoint ®elds N …t† The practical utility of the adjoint equations, like that of the forward quasi-static equations, depends on how easily they can be solved numerically. Whereas the acquisition of the adjoint ®elds i and C i is a rather straightforward task from a numerical point of view, the numerical determination of N …t† dictates a somewhat more involved and computationally more demanding procedure. The only diculty in solving the i and the C i from the inhomogeneous equations (11) and (41) arises
from the property of the operator Li lFi that it is singular in the sense that a fundamental solution (the fundamental adjoint ¯ux i ) exists for the homogeneous equation
…82† Li lFi i ˆ 0 due to which this fundamental solution i must be ®ltered out after each iterative step for determining the i
and the C i , as described in Eq. (17). Due to the singularity property of Li li Fi , the iterative solution of Eqs. (11) and (41) will feature a relatively slow convergence rate. As far as the acquisition of the N …t† is concerned, the calculational ¯ow is very similar to that for the forward equilibrium cycle iteration, except that the adjoint burnup calculations proceed backward in time. As is shown in Eq. (42), the ¯ux adjoint source Qi at ti depends on an integral of N …t† over the future time interval ‰ti ; ti‡1 Š. Further, the ®nal value of N …t† at the end of each time interval is ®xed by the discontinuity condition (48). Its magnitude depends not only on the future behaviour of N …t† but also on the C i and Pi de®ned at the ®nal time of the interval. Thus, the adjoint equations have to be solved backward in time. For solving the adjoint limit cycle iterations for determining N …t†, the ¯ow chart depicted in Fig. 2 can be implemented.

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1399


Fig. 2. Backward ¯ow chart for iterative calculation of N t‡ such that the adjoint cyclic boundary 0 condition (61) is satis®ed.

The inner loop depicted in Fig. 2 also depicts ¯ow chart for
the computational
each adjoint burnup calculation from EOC tT to BOC t‡ 0 . We note here that the N …t† can be thought of as the collection of fully decoupled individual
adjoint nuclide density ®elds NJq …t† corresponding to the response quantities NJq tT that is, the equilibrium cycle EOC nuclide density of nuclide type q in node J): n o N …t† ˆ NJq …t†; J ˆ 1; . . . ; N; q ˆ 1; . . . ; nn …83† The parallellisation of the calculation of N*(t) simply consists of distributing the calculation of the di€erent adjoint ®eld collections NJq …t† over di€erent parallel processors. An example of this is illustrated in Fig. 3.

1400

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

Fig. 3. Parallellisation of the solution of N …t† over N processors; the ®rst index indicates the response node, the second the nuclide type.

We stress that, in this study of methodology development, we allowed ourselves the simpli®cation of working with only 3 nuclide types: 238U, 235U and 239Pu and no ®ssion products taken into consideration. Nevertheless in principle, the DPT formalism presented here allows any number of nuclide types (®ssile nuclides, ®ssion products, absorber nuclides, etc.) to be taken into consideration. We indicate that
the size of the largest adjoint object to be stored, N t‡ is proportional to the 0 square of the number of nodal nodes in the system and to the square of the number of nuclide types taken into consideration. With present-day storage facilities, this should be no practical problem for most cases concerning realistic sizes of cores featuring
octant symmetry. Further, the computational e€ort involved in obtaining N  t‡ 0 is only linearly proportional to the number of nodal nodes in the system and to the number of nuclide types taken into consideration. For two-dimensional nodal core models to be used in reload pattern design, this should be manageable. Since the correction scheme (72) is only `pushing' the perturbed equilibrium limit cycle iteration in the
right direction, one could even simply choose to calculate only the part of N t‡ associated with the nuclides that play the most dominant role in 0 determining the ¯ux distributions and in the burnup process. The remaining part of
N  t‡ 0 then consists of zero elements and hence does not play a role in the correction scheme (72). Then, although the correction e€ected by (72) directly a€ects only the most important nuclides and does not directly a€ect the less dominant nuclides, its e€ect on the less dominant nuclides will be rather quickly conveyed through the transmutation scheme in which the dominant nuclides are generally positioned at the roots.

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1401

7. Example of an application: acceleration of equilibrium cycle feed enrichment minimization for a PWR core The formalism described here has been programmed in MPI-Fortran and implemented on the massively parallel computer available at Delft University of Technology. It has been applied for acceleration of a cyclic interchange optimization procedure, in which a (local) optimum for a prespeci®ed objective function is found by stepwise improvement through the evaluation of all possible octant-symmetric ternary (involving three elements) and binary (pairwise, so involving two elements) interchanges of elements in the reference fuel element trajectory representation H. The optimization case for the test 4-batch refueled PWR core containing 384 fuel elements considered was . Minimization of the required feed enrichment (FE) for the constrained equilibrium cycle, subject to the conditions that the radial power peaking factor f should remain below 1.8 and that the reloading should occur in an octantsymmetric way

With constrained equilibrium cycle we mean that the feed enrichment (symbolized here simply by the fresh fuel composition NF )is dictated by the condition that the leigenvalue at EOC should exactly equal 1. Or, in other words, the EOC-e€ective multiplication factor lT 1 the uncontrolled core (without reactivity control, such as control rods or soluble boron) is constrained to be 1. The overall assumption is that the cycle length and the power level are constant. The feed enrichment satisfying this condition is referred to as the minimal required feed enrichment (MRFE). Should, for the speci®c reload pattern in question and for the chosen power level and cycle length, the feed enrichment be chosen to be smaller than MRFE, the core would become subcritical before EOC and thus the envisaged operation cycle would not be feasible. Obviously, this MRFE is a function of the choice for the reload operator X and thus of the reload pattern H. Since in the search procedure the MRFE needs to be calculated for each candidate reload pattern H, several equilibrium limit cycle iterations (with the FE ®xed) were performed within the loop in which MRFE is determined. In this picture, the limit cycle iteration can be regarded as another type of inner iteration. The geometry of the test PWR core is illustrated in Fig. 4, in which also the reloading structure of the starting reload pattern in the search procedure (chosen intuitively to realize a low power peaking factor for satisfying the constraint f41.8) is indicated. Naturally, since the cyclic interchange optimization algorithm is inherently parallel, the search procedure has been parallellized too in the sense that the workload of evaluating all the interchanges (pattern perturbations) is uniformly distributed among the di€erent available processors. In the two-dimensional 1 12-group nodal core model (van Geemert, 1999) adopted in this study, the radial power peaking factor is de®ned as   PI ; I ˆ 1; . . . ; N …84† f ˆ max N Ptotal

1402

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

Fig. 4. Octant-symmetric test PWR core. The di€erent shades of grey indicate the di€erent fuel ages: the darker, the higher the burnup.

with pI the nodal power in node I and N the total number of nodes. In the MRFE calculation, iterations have taken place at three di€erent levels. The two innermost iteration levels are depicted in the forward ¯ow chart shown in Fig. 5 representing the equilibrium cycle iteration for a speci®c ®xed NF. We now de®ne and specify the following convergence parameters and criteriae adopted in our study. For monitoring the convergence of scheme (1) or (7) in Fig. 5, the convergence parameter "…i1 † for iterative step i1 was de®ned as   I;i1 I;i1 1 …i 1 † j; I ˆ 1; . . . ; N …85† " ˆ max j I;i1 1

For monitoring the convergence of the equilibrium limit cycle iterations, the coni2 † for iterative step i2 was de®ned as vergence parameter "…EC ( …i † ‡
) i2 1† ‡
Nq;I2 t0 t0 N…q;I …i 2 † j; I ˆ 1; . . . ; N; q ˆ 1; . . . ; nn ; …86† "EC ˆ max j i2 1† ‡
t0 N…q;I and for monitoring the convergence of the FE iterations, the convergence parameter i3 † "…FE for iterative step i3 was de®ned as

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1403

Fig. 5. Forward ¯ow chart for the iterative calculation of N…t† such that the cyclic boundary condition is satis®ed and thus the equilibrium cycled is reached, for a speci®c ®xed NF.

(

i3 † "…FE

l 1 …i3 † lT 1 …i3 ˆ max j T lT 1 … i 3 1†

1†

) j; j1

1

lT …i3 †j

…87†

The following convergence criteria were applied for evaluating both the reference and the perturbed reload pattern: " 410 6 ; "EC 410 5 ; "FE 410

5

…88†

For solving the adjoint ®eld equations, similar convergence parameters and criteria were adopted as for the forward equations. In order to establish "FE 410 5 , generally about 11 or 12 outer iterations were required, using the outer iterative scheme

1404

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

FEi ˆ FEi

1‡

1 lT 1 … i 1† …FEi lT 1 …i 1† lT 1 …i 2†

1

FEi 2 †

…89†

Some examples have been isolated of reload operator perturbations of di€erent magnitudes, involving binary, ternary, 6-fold and 10-fold interchanges of elements in the reference fuel element trajectory representation H. For these examples and in the entire optimization procedure, we have calculated the perturbed constrained equilibrium cycle solutions using the formalism described in this paper, so by applying the system (7) 8 D E D E 
…m‡1†  > > ˆ …m † j L0 l…m† F0 …m† a 1 j …m † 1 > > D E <
 j L0 l…m† F0 …m† …90† > D E > l…m‡1† ˆ l…m† ‡ > >  j F0 …m† :

in order to obtain l0 ; 0 ˆ l ‡ l; ‡  , instead of using the power method for solving …L lF† ˆ 0 (Eq. (1)), and by applying the EOC correction 

D

m† 0 …m† N tT ˆ …m† N tT ‡ N t‡  N tT j X …corr 0

…m

1†

N tT


E

prior to each implementation of the reload operation that can be written in terms of nuclide density perturbations as    


0 …m † …m‡1† N t‡ N t N N t …91† ˆ X N ‡ X ‡ NF F F T corr T 0 In Tables 1 and 2, the numbers of iterative steps required for convergence as de®ned in (88) are listed for the di€erent considered cases, consisting of reload operator perturbations of di€erent magnitudes, involving binary, ternary, 6-fold and 10-fold interchanges of elements in the reference octant-symmetric fuel element trajectory representation H. For the large ®ctitious 4-batch refueled PWR core illustrated in Fig. 4 and starting from the arbitrary reload pattern (chosen intuitively to realize a low power peaking factor for satisfying the constraint f41.8) that is illustrated in this picture, we have searched for the loading scheme associated with the lowest MRFE subject to f41.8 $ and the EOC constraint lT ˆ 1. For this, the ternary interchange optimization procedure was applied. The search results are illustrated graphically in Fig. 6 and Fig. 7. The ®nal search result (that is, the reload pattern corresponding to the lowest encountered MRFE subject to f41.8) is plotted in Fig. 8. In Table 3, the improved MRFEs and their associated power peaking factors are listed as a function of the number of adjoint ®eld adjustments (= number of forward MRFE iterations). Also indicated are the required cumulative processing times on the parallel `Vermeer' CRAY-T3E, using 52 processors and applying the GPT-

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1405

Table 1 Average numbers of required `level 1' iterative steps for solving the ¯ux eigenvalue equations for " 410 6 , with and without application of the single time step GPT method Magnitude of permutation

GPT-accelerated

Non-accelerated

Binary Ternary 6-Fold 10-Fold

4 6 8 10

50 50 50 50

Table 2 Average numbers of required `level 2' iterative steps for achieving equilibrium cycle convergence
"EC 410 5 , with and without application of the EOC-correction (72) Magnitude of permutation

GECDPT-accelerated

Non-accelerated

Zero Binary Ternary 6-Fold 10-Fold

3 4 4 5 6

12 12 12 12 12

Fig. 6. Convergence towards the minimized MRFE as a function of the cumulative number of adjoint ®eld adjustments.

1406

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

Fig. 7. Convergence towards the minimized MRFE.

Fig. 8. The octant-symmetric reload pattern corresponding to the minimized MRFE (subject to f41:8). The darker the colour of the fuel assemblies, the higher their burnup.

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1407

Table 3 Listing of the search information Cumulative number of adjoint ®eld adjustments

MRFE (percent)

Power peaking factor

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

1.5800 1.5690 1.5613 1.5570 1.5537 1.5515 1.5499 1.5483 1.5469 1.5458 1.5445 1.5426 1.5421 1.5408 1.5404 1.5401 1.5400 1.5398 1.5395 1.5394 1.5393 1.5392 1.5391 1.5391

1.4579 1.6494 1.7563 1.7476 1.7561 1.7710 1.6214 1.6327 1.6488 1.7665 1.7157 1.7027 1.6828 1.7271 1.7313 1.7669 1.7618 1.7617 1.7607 1.7655 1.7657 1.7660 1.7619 1.7621

Cumulative CPU-time (s) using GPT-based acceleration and 52 parallel processors

± ± ± ± ± ±

740 1480 2220

7400 ± ± ± ± ± ± ± ± ± ± ± ± 17,020 17,760

Estimate of cumulative CPU-time (s) using no GPT and no parallellization

± ± ± ± ± ±

648,100 1,296,200 1,944,300

6,481,000 ± ± ± ± ± ± ± ± ± ± ± ± 14,906,300 15,554,400

based accelerative techniques. The last two columns contain estimates of the cumulative CPU-times that would have been required in case no parallellization or/and no GPT-based accelerative techniques would have been applied. We note that a certain fraction of the cumulative CPU time in the 4th column arises from the computational e€orts invested in solving the adjoint ®eld equations. 8. Analysis of the achieved speed-up In our sequential implementation for the core depicted in Fig. 4 (featuring 52 response nodes in the octant), the solution of the i took about as much CPU-time
as one reference forward MRFE iteration, whereas the acquisition of N t‡ 0 took about 5 times as much time, per response nuclide type. So obviously (and as expected), the computation of the adjoint ®elds certainly requires a non-trivial amount of computational e€ort! However, for a design purpose such as reload pattern optimization, this investment pays o€ extremely well in the case where there are very large

1408

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411


Fig. 9. Speedup behaviour in case of three nuclides included in N t‡ 0 .

numbers of pattern changes to be evaluated. If these evaluations can be accelerated considerably using the obtained adjoint ®elds (as indicated by the results in Tables 1± 3), the gain in computational eciency will be very signi®cant indeed. Suppose for
 ‡ example that we have (nn) nuclide types to be included in N t . The CPU-time 0
(T symbolizes required for calculating N t‡ is then about [5(nn)+1]T MRFE MRFE 0 the average total time required for determining the minimally required feed enrichment associated with a certain reload pattern, without the use of perturbation theory in any of the three iteration levels). However, if the average acceleration due to the use of the adjoint ®elds is a1a2 and a parallellization over Z processors is applied, the gain G in computational eciency for the case of evaluating M perturbed reload patterns can be expressed by 8 za 1 > ; …only level …1† acceleration† > a > 1 > < 1 ‡ ‰1 ‡ ZŠ M …92† Gˆ Za1 a2 > > … level … 1 ‡ 2 † acceleration † > > : 1 ‡ a1 a2 ‰5…nn† ‡ 1 ‡ ZŠ M
For very large M the limiting case limM ! 1 G Z; M; a1;2; …nn† ˆ Za1 a2 will be approached. Obviously, a fully sequential approach (no parallelization) corresponds to Z=1. In Figs. 9 and 10, the logarithm of G is plotted as a function of M for the case a1=10, a2=2 and for di€erent combinations of speed-up mechanisms applied.

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1409


Fig. 10. Speedup behaviour in case of 40 nuclides included in N t‡ 0 .


For Figs. 9 and 10, the numbers of nuclides types accounted for in N t‡ 0 are 3 and 40, respectively. In these pictures one can
see that it is not always by de®nition advantageous to  calculate and apply N t‡ 0 in addition to calculating and applying the i . In the case of only 3 nuclide types accounted for, this already o€ers computational advantage for low values of M (starting from M=400), whereas in the case of the large number of 40 nuclide types it would not become advantageous until M exceeds 4000. However we should recall our earlier remark that for speed-up purposes one
does not necessarily have to include all of the present isotopes in N t‡ 0 , but only
the most dominant ones. In this case, less time has to be spent on calculating N t‡ 0 at the cost of a limited loss of accelerative power. 9. Conclusions We conclude that the GPT-based accelerated iterative methods described in this paper are well applicable to a heuristic equilibrium cycle optimization procedure that is based on assessing the e€ects of many di€erent modest permutations in the loading scheme. By this application, a considerable reduction of the computational costs for the limit cycle iterations is realized. The price to be paid, the necessity of solving an extensive set of adjoint equations a few times during the search procedure, will generally pay o€ easily if the evaluations of several thousands of candidate schemes can

1410

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

be signi®cantly accelerated using the adjoint ®elds. Of interest here is that, for these GPT-based schemes, the more favourable convergence properties arise from the inherently more favourable mathematical structure of the iterative equations in the case of modest perturbations. When setting up the accelerated iterative expression for the equilibrium cycle, an adjoint reload operation invariance equation emerges for the adjoint ®elds. Due to this, the acquisition of the adjoint time-dependent nuclide density distribution consists of an iterative procedure oriented backwards in time, featuring solution of adjoint burnup equations and application of adjoint reload operations until sucient convergence has been realized. Both the solution of the adjoint ®eld equations and the implementation of the applied heuristic optimization algorithm could be parallellized and executed on the parallel `Vermeer' CRAY-T3E available at Delft University of Technology. From the speed-up that could be realized, it can be concluded that the combination of parallellism and generalized perturbation theory o€ers the opportunity to perform exhaustive, fast and accurate sampling of the solution space for the equilibrium cycle optimization problem. Acknowledgements The authors would like to express special gratitude to the Delft Center for HighPerformance Computing (HPAC) for the use of its `Vermeer' CRAY-T3E massively parallel computer. References Axmann, J.K., 1997. Parallel Adaptive Evolutionary Algorithms doe Pressurized Water Reactor Reload Pattern Optimizations. Nuclear Technology 119, pp. 276±292. DeChaine, M.D., Feltus, M.A., 1996. Fuel Management Optimization Using Genetic Algorithms and Expert Knowledge. Nuclear Science and Engineering 124, 188±196. Gandini, A., 1975. Time-Dependent Generalized Perturbation Theory for Burnup Analysis. CNEN RT/ FI(75) 4, Comit. Nazionale per l'Energia Nucl., Rome. de Jong, A.J., 1995. Reload pattern Optimization for Batch Refuelled Nuclear Reactors, report IRI-13195-010, Interfaculty Reactor Institute, Delft University of Technology, Delft, The Netherlands. Kropaczek, D.J., Turinsky, P.J., 1991. In-Core Nuclear Fuel Management Optimization for Pressurized Water Reactors Utilizing Simulated Annealing. Nuclear Technology 95, 9±32. Maldonado, G.I., Turinsky, P.J., Kropaczek, D.J., 1995. Employing nodal generalized perturbation theory for the minimization of feed enrichment during pressurized water reactor in-core nuclear fuel management optimization. Nuclear Science and Engineering 121, 312±325. Moore, B.R., Turinsky, P.J., 1998. Higher order generalized perturbation theory for boiling water reactor in-core fuel management optimization. Nuclear Science and Engineering 30, 98±112. Parks, G.T., 1987. An Intelligent Stochastic Optimization Routine for In-Core Fuel Cycle Design. Transactions of the American Nuclear Society 57, pp. 259±260. Parks, G.T., 1996. Multi-Objective Pressurized Water Reactor Reload Core Design by Nondominated Genetic Algorithm Search. Nuclear Science and Engineering 124, 178±187. Smuc, T., Pevec, D., Petrovic, B., 1994. Annealing Strategies for Loading Pattern Optimization. Annals of Nuclear Energy 21, pp. 325±336. Stevens, J.G., Smith, K.S., Rempe, K.R., 1995. Optimization of Pressurized Water Reactor Shu‚ing by Simulated Annealing with Heuristics. Nuclear Science and Engineering 121, 67±88.

R. van Geemert, J.E. Hoogenboom / Annals of Nuclear Energy 28 (2001) 1377±1411

1411

van Geemert, R., De Leege, P.F.A., Quist, A.J., Hoogenboom, J.E., Gibcus, H.P.M., 1998a. Research reactor in-core fuel managament optimization by application of multiple cyclic interchange algorithms. Nuclear Engineering and Design 186, 369±377. van Geemert, R., Hoogenboom, J.E., Quist, A.J, 1998b. In-Core Fuel Management Optimization by incorporation of a Generalized Perturbation Theoretical Approach in a Heuristic Search Procedure. Proceedings PHYSOR98 conferentie, V1.91±101, Long Island, USA. van Geemert, R., 1999. Nuclear Reactor Reload Pattern Optimization by Application of Heuristic Search and Perturbation Theoretical Methods. PhD thesis, ISBN 90-407-1867-9/CIP, Delft University Press, Delft, The Netherlands. Verhagen, F.C.M, van der Schaar, M., 1993. Simulated Annealing in LWR Fuel Management. TOPNUX'93 Proceedings 2, pp. 37±39. Williams, M.L., 1979. Development of depletion perturbation theory for coupled neutron/nuclide ®elds. Nuclear Science and Engineering 70, 20±36. White, J.R., Avila, K.M., 1990. Developing Feasible Loading Patterns Using Perturbation Theory Methods. Proc. PHYSOR90 conference, Marseille, France. Yang, W.S., Downar, T.J., 1989. Depletion perturbation theory for the constrained equilibrium cycle. Nuclear Science and Engineering 102, 365±380.