© A. Gandini, published by EDP Sciences, 2021

1 The steps of the HGPT method development

The interest in perturbation/sensitivity methodologies for reactor physics studies started during the stage of Augusto Gandini as an associate researcher by the Reactor Physics Division of the Argonne National Laboratory in years 1961–62 [1]. At that time the Division Director Robert Avery posed the question on the possibility of developing perturbation methods apt for analyzing functionals of the neutron flux, like reaction rate ratios, as we were using perturbation methods for reactivity calculations. It was in later years that a research activity, focused on this subject [2], started together with other researchers, among which Massimo Salvatores, at the Casaccia Center of the Italian National Committee of Nuclear Energy (CNEN). The result of this activity was the methodology that became known as Generalized Perturbation Theory (GPT). Its first relevant application was presented in 1966 at the ANL Conference on Fast Critical Experiments and their Analysis [3]. The method considered adopted heuristic concepts based on conservation principles [4], extending the method originally proposed for reaction rate ratios by L.N. Usachev to any functional of the real and adjoint neutron fluxes, such as reactivity coefficients, effective prompt neutron lifetimes, effective delayed neutron fractions¹. A number of codes implementing this methodology were written at that time by the Fast Reactor Physics Group at CNEN for the calculation of sensitivities and data correlations [3,5]. For a number of years these codes have been a reference for GPT sensitivity analysis applications.

This heuristic methodology, initially limited to the neutron domain, was successively further extended to the nuclide field [6,7] enabling, in particular, sensitivity studies relevant to fuel depletion evolution.

One of a main applications of the GPT methodology has been the differential data adjustments based on integral quantities data measurements in critical facilities (ZPRs, MASURCA, etc.). To this purpose a probabilistic inference based methodology was adopted. This same methodology led to the definition of correlation coefficients between integral quantities relevant to experimental facilities with respect to those relevant to reference reactors [8–10].

In further developments the heuristically based GPT methodology has been applied:

• to nonlinear problems [11], in particular, to the coupled neutron/nuclide field for reactor cycle analysis;

• in the estimation of spatial shifts of power peak points following a perturbation [12];

• to reactor design optimization [13];

• to the development of the EGPT methodology [14], by which, for the analysis of reactivity coefficients, the calculation of the importance functions, implying the solution of inhomogeneous equations, is replaced by the calculation of functions, solution of simpler homogeneous ones governed by a properly modified operator;

• to the analysis of subcritical (ADS) reactors [15]. This led to the definition of ‘generalized reactivity’, properly taking into account the intensive control variable required (for instance, the neutron source strength) for maintaining the established power level. This ‘generalized reactivity’ appears, in particular, in the derivation of the point kinetic equations governing the normalized power of a subcritical system [16,17]:

• to the development of a method for sensitivity analysis in the reliability domain [18];

• to the development of a method by which the information obtained on-line through a system of neutron measuring devices such as self-powered neutron detectors (SPNDs) inserted in the core of a nuclear power reactor allows the on-line detection of a possible hot spot during plant operation [19,20].

2 The HGPT method

Since the beginning of nuclear reactor physics studies, perturbation theory has played an important role. As well known, it was first proposed by Wigner [21] to study fundamental quantities such as the reactivity worths of different materials in the reactor core. It is also well known that this first formulation, today widely used for reactor analysis, makes a consistent use of the classical adjoint flux concept.

The HGPT approach, fundamentally based on importance conservation principles, is generally intended to be used for defining perturbation expressions relevant to a variety of responses in stationary, as well as time dependent, linear and nonlinear fields. Compared with direct calculations, it allows to obtain with relatively limited effort the sensitivity coefficient of the responses of interest, and may then be very well used, as it has been extensively done, in a variety of studies.

The concept of importance was first defined by Kadomtzev [22], in the radiation transport field, and then considered by Usachev [4], as the contribution to a response in a critical system by a neutron through its progeny, and by Lewins [23], as interpretation of the adjoint function adopted in variational techniques.

To be also reminded that the HGPT approach and the variational one [23,24] differ only in the procedure for arriving at the sensitivity/perturbation expressions, these resulting equivalent to each other as demonstrated for different fields by Greenspan [25] and Marques Alvim et al. [26].

The preference, or the merit of each method for arriving at the perturbation/sensitivity expressions of interest is beyond the scope of this presentation. We limit to note that, whereas with the variational approach one makes use of the mathematical concept of adjoint function and of its properties, and then realizes that it may be associated with the importance, with the HGPT approach one starts with defining this latter quantity and heuristically, via conservation principles, arrives at the equation governing it and, finally, at the sensitivity/perturbation expressions. The importance may be subsequently associated with the adjoint function, while the reciprocity relationships implied with the HGPT theory may be associated with its properties.

2.1 The importance function

In the HGPT method the importance function is uniquely defined in relation to a given system response, for example, a neutron dose, the quantity of plutonium in the core at end of cycle, the temperature of the outlet coolant.

The HGPT method was first derived in relation to the linear neutron density field. Then it was extended to other linear ones. For all these fields the equation governing the importance function was obtained directly by imposing that on average the contribution to the chosen response from a particle [a neutron, or a nuclide, or an energy carrier] introduced at a given time in a given phase space point of the system is conserved through time (“importance conservation principle”). Obviously such importance will result generally dependent on the time, position, and, when the case, energy and direction, of the inserted particle.

Consider a linear particle field density represented by vector f (e.g., the multigroup neutron density field) and a response Q of the type² (1)where s ⁺ is an assigned vector function and where <> indicate integration over the phase space. Weighting all the particles inserted into the system, let's assume a source s, with the corresponding importance (f ^*) will obviously give the response itself, that is,(2)which represents an important reciprocity relationship.

From the first derivations mentioned above the rules for determining the equation governing the importance function f ^* were learned. They imply, with respect to the equation governing the real function f:

• change of sign of the odd derivatives,

• transposing matrix elements,

• reversing the order of operators,

• substitution of the real source s with s ⁺.

The first three rules will be generally called “operator reversal” rules.

The HGPT method was then extended to any field governed by linear operators for which the rules for their reversal were known. In particular, it was extended to the derivative fields, obtained from expanding to first order, around a given starting solution, a number of important nonlinear equations, as those governing:

• the coupled neutron/nuclide field, relevant to core evolution and control problems,

• the temperature field, relevant to thermohydraulics.

2.2 General formulation

Consider a generic physical model defined by a number of parameters p_j (j = 1, 2, ..., J) and described by an N-component vector field f obeying equation(3)

Vector f(q,t) generally depends on the phase space coordinates q and time t. Vector p represents the set of independent parameters p_j (j = 1, 2, ...) fully describing the system and entering into equation (3). Their value generally determines physical constants, initial conditions, source terms, etc. Equation (3) can be viewed as an equation comprising linear, as well as nonlinear, operators and is assumed to be derivable with respect to parameters p_j and (in the Frechet sense) component functions f_n (n = 1, 2, ..., N).

Consider now a response of interest, or functional Q given by equation (1). In the following, we shall look for an expression giving perturbatively the change δQ of the response Q in terms of perturbations δp_j of the system parameters. In particular, expressions giving the sensitivity coefficients relevant to each parameter p_j will be obtained.

Expanding equation (3) around a reference solution gives, setting ,(4) where O₂ is a second, or higher order term, and a (Jacobian) operator given by the expressionwhere by we have indicated a Frechet derivative [27].

Since parameters p_j, and then their changes δp_j, have been assumed to be independent from each other, it must follow(5)which represents the (linear) equation governing the derivative functions f_/j. The source term m_/j is here intended to account also, via appropriate delta functions, for the initial and, if the case, boundary conditions.

Consider now functional(6)

Introducing the importance (f ^*) associated with field f_/j, if we use it as weight of the source term m_/j, and integrate space- and time-wise, according to the source reciprocity relationship, equation (2), the resulting quantity will be equivalent to functional Q_j, that is,(7)where the importance f ^* obeys the (index-independent) equation(8)

H ^* being obtained by reversing operator H. As said above, this implies transposing matrix elements, changing sign of the odd derivatives, inverting the order of operators.

We can easily see that the sensitivities s_j (j = 1, 2, ..., J) of system parameters can be written(9)where the first term at the right-hand side represents the so called, easy to calculate, direct term.

The overall change δQ due to perturbations δp_j (j = 1, 2, ..., J) of system parameters can be written, to first order,(10)

2.3 The HGPT applied to critical and subcritical system analysis

Due to its generality, the HGPT approach may be applied without ambiguity to the analysis of critical, as well as subcritical systems. Indication on this latter possibility were made from the late 60's [28–31].

Inherent with the HGPT theory is the concept of the control to be associated with a system generally subject to a power, or flux level constraint. Cacuci says [32]: “Classical perturbation theory certainly cannot be used for optimal operation and control of an ADS”. But, while this is true in relation to the classical perturbation theory, “a la Wigner”, in which the control role is fictitiously played by the coefficient (λ) multiplying the fission source, the same thing cannot be said with respect to the HGPT theory which may account for the real control (the insertion of a control rod, the intensity of the extraneous neutron source strength, etc) [33,34]. Introducing the general frame of optimal control theory does not change this conclusion. As well known, control theory intrinsically uses the sensitivity coefficients of responses of interest, among which the target quantity, during the various steps of an optimization search [35]. The issue here is that of adopting a correct, unbiased sensitivity theory accounting for the real control adopted, in a critical, as well as in a subcritical reactor system. Which is just what HGPT theory is doing.

The above concept of control-related HGPT theory in relation to subcritical systems has been considered in some depth since 1997 [15,16]. Special attention has been given to a response represented by the very control variable, in particular, the source strength itself. This gives rise to a peculiar perturbation expression by which it is possible to evaluate control changes (for instance, at the end of reactor cycle life, in an evolution study) following a perturbation of any system parameter.

2.3.1 The HGPT methodology at quasi-static conditions

The methodology relevant to long term nuclide/neutron core cycle evolution analysis (see Appendix A) may be very well applied to source driven, subcritical systems.

Consider the quasi-static equations governing the neutron density n, the nuclide density c andthe control function ρ during the core evolution during burn-up, in the form(11) (12) (13)where B and E depend on fuel and neutron densities c and n, respectively.

Since we generally consider systems at quasi-static, that is, stationary conditions, the time derivative at second member of equation (11) may be neglected in the course of the integration process.

Any response, functional of variables n, c, and ρ, could be considered for analysis. We think instructive to limit here consideration to the response defined by the expression(14)which corresponds to the relative source strength required at t_F to assure the power level imposed. We may assume that, at unperturbed conditions, ρ (t) = 1 in the interval (t_o, t_F). If some system parameter (for instance, the initial enrichment, or some other material density) is altered, as in an optimization search analysis, it may be of interest to evaluate the corresponding change of ρ at the end of cycle, to make sure that given upper limit specifications of the source strength are not exceeded.

Along with the HGPT methodology, the equations for the corresponding importance functions result(15) (16) (17)

and being coupling operators defined in Appendix A.

Equation (17) corresponds to an orthonormal condition for n ^*.

In order to determine the “final” value n ^*(t_F) required for starting the integration of equation (15), in consideration of the nature of the above governing equations, we shall first write n ^* and ρ^* in the form³ (18) (19)with and being finite functions, vanishing at t_F.

Replacing into equation (15), integrating in the interval (t_F − ε, t_F + ε), and then making ε → 0, we obtain the equation(20)

Let us now define as obeying equation(21)

We note that corresponds to the importance relevant to functional <c(t_F), Sn(t_F)>, that is, to the system power W. From the source reciprocity relationship (Sect. 2), we may write(22)

From constraint, equation (17), we easily obtain(23)

and then(24)

From this “final” value, a recurrent calculation scheme may be defined starting from t_F and proceeding backward.

Along with the HGPT methodology, the sensitivity coefficient relevant to the k'th parameter p_k is found as(25)

with given by equation (23). The first term at right side accounts for effects on ρ(t_F) due to parameter changes at t_F, in particular, if p_k ≡ W, it gives the (trivial) result . The second, integral term accounts for analogous effects on ρ(t_F) produced by parameter changes at times t<t_F.

Rather than on the source term, a control on the neutron absorption in the multiplying region could be of interest. In this case, the (intensive) control variable ρ would represent the average penetration of the control elements, or the average density of the soluble boron in the coolant, and then would enter into the (transport, or diffusion) operator B. The orthonormal condition for the neutron importance n ^* would now be, rather than equation (17),(26)

In this case, the sensitivity coefficient with respect to a given parameter p_k would always be given by equation (25), with obeying equation (21), but with(27)

In general, a control strategy, by which an automatic resetting of the imposed overall power is actuated, might imply a control intervention on both the neutron source strength and the absorbing elements within the multiplying region. In this case, ρ (which remains a unique, intensive control variable) would affect both operator B and the neutron source [in this latter case, via an appropriate ρ- and parameter dependent coefficient α(ρ|p), assumed unity at unperturbed conditions]. The distribution between these two control mechanisms could be described by appropriate parameters(subject to perturbation analysis). The sensitivity coefficient, in this case, with respect to a given parameter p_k would always be given by equation (25), with obeying equation (21), but with(28)

2.3.2 Stationary case

To study a given subcritical system at stationary conditions (which may be interpreted at the beginning of its cycle life), we may consider the same system above in which the neutron source and the nuclide density are assumed time-independent during an arbitrarily small time interval (t_o, t_B). We assume that at t_o the neutron density (n_o), as well as the control (ρ_o) have already reached stationary conditions. So, also these two quantities are time-independent in the same time interval. Their governing equations can then be written, in case the power level is controlled by the source strength,(29) (30)

Also here we shall assume that at unperturbed conditions ρ_o = 1.

The same equations derived previously are applicable to this case, with the advertence of replacing t_F with t_B and setting the coupling operators and appearing in equations (15) and (16) equal to zero. The sensitivity coefficient of the response ρ(t_B) [=ρ(t)=ρ_o, that is, constant in the whole interval (t_o, t_B)] with respect to the j'th parameter p_k can then be obtained. Since in this case c ^*, as well as and vanish, recalling equation (25), we obtain(31)

where(32)and obeys equation(33)

If, rather than via the source strength, the power level reset control is assumed to be regulated via neutron absorption, so that the control ρ_o would enter into operator B, the sensitivity coefficient would be given always by equation (31), but with(34)

We might as well consider a (fictitious) control mechanism affecting the fission source, rather than the neutron absorption, that is, we might choose as control a coefficient multiplying the fission matrix (F) and, therefore, entering into the Boltzmann, or diffusion, operator B (=A + ρ_oF). The sensitivity coefficient would be given again by equation (31), but with(35)

2.3.3 Reactivity of subcritical systems

For resetting the power level, we have considered above different control mechanisms to which the following types of equations governing the neutron density may be associated:(36) (37) (38)

where the control and parameter dependence is indicated. Coefficient α is given and reflects the mixed strategy chosen. Equations (36)–(38⁴) may be generally represented by equation(39)

The sensitivity expression (31) may be generalized so that(40)with obeying equation (33).

A corresponding perturbation expression may now be obtained. Assuming that the power W_o appearing in equation (40) is not subject to perturbation, we may write:(41) where and .

As said previously, δρ_o corresponds to the control change necessary to reestablish the power level existing before the perturbation δm_(n,o). We may as well say that the perturbation δm_(n,o) [and δ(S^T c_o)] would produce a power level change equivalent to that produced by a control change δK_ρ given by the equation(42)

In the case of a (fictitious) control on the neutron fission, setting λ in place of ρ to distinguish this peculiar case, we may explicitly write(43)

The first term at the right side closely resembles the reactivity expression for critical systems⁵. So, we shall call a quantity δK_λ as given by expression (43) a “generalized reactivity”. The second term may be defined the “source reactivity”, whereas the last one a “direct effect”. To account for a generic ρ-mode control mechanism, we shall extend this definition to δK_ρ, similarly defined by equation (42), that is(44)

and call it generalized ρ-mode reactivity.

3 Conclusions

In concluding this presentation I like to remind that also the classical perturbation theory, a la Wigner, may be interpreted along the schemes of the HGPT methodology. The adjoint neutron flux (defined up to an arbitrary multiplicative factor) may in fact be interpreted, recalling Soodak [36], as an importance function corresponding to an arbitrary functional linear with the real flux, defined at an asymptotic time. It may be, in particular, a functional proportional with the asymptotic power.

The HGPT may then be defined as a heuristic methodology on the basis of which perturbative formulations can be derived in any linear, or linearizable field interpretable (really or fictitiously) as a density one.

For its heuristic character, based primarily on the importance conservation principle and on reciprocity relationships, the HGPT methodology offers also the advantage of being easily grasped and used in a variety of situations.

For instance, concerning its application to subcritical systems, it has been shown [37] how, by properly adjusting the external source strength for compensating a slight calibrated control rod insertion so to reset the given nominal power, the subcriticality index (1−K_sub) can be easily estimated online during reactor operation without significantly altering the power level.

Appendix A

The densities of the nuclides that make up the fuel refer to average values relating to the macrozones into which the core has been divided. They are represented with the vector c_z. [z = 1, 2, …, Z (number of sub-areas)]. The density of neutrons, dependent on phase space and time, is given in the multigroup form by vector n. In this discussion, the approximation of diffusion has been adopted (which is widely used for burn-up calculations). An intensive control variable, ρ(t), is associated with the densities n(r,t) and c_z(r), defined in the interval (t_o, t_F), such as to keep the total assigned power W(t) constant. The variable ρ(t) can represent, for example, the overall degree of insertion of the control rods into the core (not their relative movement, which can generally be described by parameters), or the average density of a neutron poison in the coolant. In the calculations a fictitious control parameter is normally adopted, for example a coefficient that multiplies the fission source. The application of the methodology to different control modes is however possible, this implying, as will be seen below, different conditions on the functions that describe the behavior of the system [1]. In a subcritical system (ADS), ρ(t) can represent the intensity of the source (by adjusting the intensity of the accelerator current). The nonlinear equations which the variables n, c and ρ must satisfy can therefore be formally written, in the most general case, as

(A.1)

(A.2) (A.3)where <> sys indicates integration over the entire volume of the multiplying zone, B represents the operator in diffusion or transport (depending on c and, generally, on ρ), E the evolution matrix of the nuclides (depending on n), vector p the system parameters. s_n and s_c are source terms, while S and Γ are the matrices⁶

(A.4)

(A.5)where is the microscopic fission cross section of the isotope m in group g, Γ is a diagonal matrix with elements γ_m representing the amount of energy per fission of the m-th element while V is the diagonal matrix of the velocities v_g of the neutrons. The quantities γ_m, v_g, W and are generically represented as system parameters (p_j). The terms in which the Dirac deltas appear represent initial conditions.

If, for simplicity, we wish to replace the condition of constant power with the condition of constant fission rate, it will be sufficient to replace the matrix Γ with a unitary matrix.

In the following we will consider the general case of both critical and subcritical multiplying systems.

The source term s_c in the second member of equation (A.2) is generally given by a sum of delta functions defined at specific times to take into account fuel loading and shuffling operations.

In equation (A.1), control over the source was chosen for the subcritical case. If the control falls within operator B, the methodology would become similar to that relating to the critical reactor. In quasi-static problems, such as those affecting burn-up studies, the time derivative is negligible. Its notation is still maintained to allow, as we will see later, the determination of the correct operator that holds the importance function.

A very general form of a Q (linear) response can be written in the form(A.6)with given quantities.

From the linearization procedure, and recalling the coordinate complementation rule [1], we can obtain the linear equations that govern the derived functions and the importance functions. The system of equations of derived functions results:

(A.7)

where and are the coupling operators and , respectively, ξ_c is a coefficient equal to 1 within the limits of the space in which the control parameter ρ is defined and otherwise equal to zero, and where ωρ is given by the expression(A.8)

The importance functions: they are instead governed by equations

(A.9)

where , are the adjuncts of the operators , Ω _n and Ω _c , respectively. By placing(A.10)the equation corresponding to the first line, relating to the importance of neutrons, is:(A.11)while the equation corresponding to the second line, relating to the importance of nuclides, is:(A.12)

The equation corresponding to the third line, relating to the importance associated with the power, is, recalling (A.8),(A.13)

In case it is , this equation indicates how the neutron importance function is orthogonal in the phase space to the distribution of the neutron density in critical systems, while in the subcritical ones it is orthogonal to the neutron source distribution. For critical systems we may say that the importance function results void of the fundamental mode.

It should be noted that the role of this function is relatively marginal compared to that of . In fact, while the latter takes into account the effects of a perturbation on the power level during evolution and therefore of the necessary modulation of the flux to maintain the level imposed by the control parameter, the function only takes into account the effects of any alterations of the spatial and energetic distribution of the neutron density on the quantities (average micro_reaction_rates in each macrozone) that appear in the evolution matrix. It is instructive to observe how in the case of the approximation to a single group and zero-dimensional geometry (as in the cell calculation) the value of is equal to zero.

Assuming = 0, multiplying the first row on the left by n and integrating in the phase space, it results(A.14)

The general expression of the variation δQ resulting from a perturbation of the system parameters can be written in the form(A.15)

References