Richard Sanchez
email: sanchezr@snu.ac.kr

Daniele Tomatis
email: daniele.tomatis@cea.fr

Igor Zmijarevic
email: igor.zmijarevic@cea.fr

Kevin Dugan
email: kevin.dugan@cea.fr

Analysis of Alpha Modes in Multigroup Diffusion

Keywords: neutron kinetics, alpha modes, clustering

The alpha eigenvalue problem in multigroup neutron diffusion is studied with particular attention to the theoretical analysis of the model. Contrary to previous literature results, the existence of eigenvalue and eigenflux clustering is here investigated without the simplification of a unique fissile isotope or a single emission spectrum. A discussion about the negative decay constants of the neutron precursors concentrations as potential eigenvalues is provided. An in-hour equation is derived by a perturbation approach recurring to the steady state adjoint and direct eigenvalue problems of the effective multiplication factor and is used to suggest proper detection criteria of flux clustering. In spite of to prior work, the in-hour equation results for a necessary and sufficient condition for the existence of the eigenvalue-eigenvector pair. A simplified asymptotic analysis is used to predict bands of accumulation of eigenvalues close to the negative decay constants of the precursors concentrations. The resolution of the problem in one-dimensional heterogeneous problems shows numerical evidence of the predicted clustering occurrences and also confirms previous theoretical analysis and numerical results.

I. INTRODUCTION

The kinetic neutron equation with precursors is an efficient tool for the analysis of the neutron time evolution and, therefore, finds multiple applications in reactor control and the study of accident scenarios. The time-dependent flux depends on the initial conditions for neutrons and precursors, the boundary condition for neutrons, the presence of external sources and the time-dependent changes of the cross sections, such as those produced by rod motion and changes of Boron concentration or to small stochastic perturbation induced by the coolant flow and, in the long time, by radioactive decay and nuclide depletion and creation from fission.

However, if the cross sections and boundary conditions remain constant in time and the sources vanish, then the state of the system tends asymptotically in time to an exponential behavior, which is independent of the earlier changes of the system and, specially, of the initial conditions. [START_REF] Wing | An introduction to transport theory[END_REF] The exponential behaviors that a given system can adopt are the solutions of an eigenvalue equation and can be used to describe the fast evolution of the system as well as to characterize the reactivity of the system. The solutions of this equation are known as time-dependent modes or, more simply, alpha modes, where "alpha" refers to the most frequently adopted symbol for these eigenvalues.

Alpha modes have been applied to formally derive different forms of the well-known in-hour equation, to obtain solutions of the kinetic equations by expansion techniques, [START_REF] Foulke | Application of the natural mode approximation to space-time reactor problems[END_REF] to develop numerical solution methods [START_REF] Borsari | Time-varying modes in space-dependent nuclear reactor kinetics by a perturbation approach[END_REF] and also as weighting fluxes to homogenize the kinetic equation. [START_REF] Dugan | Cross-Section Homogenization for Reactivity-Induced Transient Calculations[END_REF] These applications as well as the theoretical interest of alpha modes have * Corresponding author been the object of intense study. However, with the exception of a few scattered mathematical results, [START_REF] Porsching | On the spectrum of a matrix arising from a problem in reactor kinetics[END_REF][START_REF] Asahi | Theory of Omega-d Modes[END_REF][START_REF] Devooght | Spectrum of the multigroup-multipoint diffusion operator with delayed neutrons[END_REF] a detailed description of these modes have not yet been given for the diffusion equation. The purpose of this paper is to give a detailed analysis of the alpha modes for a slab geometry using multigroup diffusion theory. We give ample numerical evidence for all the modes predicted from mathematical analysis as well as from physical arguments.

General equations for the alpha eigenvalue problem are discussed in Sec. II, including a perturbation expression in terms of reactivity. In Sec. III we summarize our results and observations for the multigroup one-dimensional slab. Some numerical results are illustrated in Sec. IV, while in Sec. V we compare with Asahi's analytical results for the one-group problem. [START_REF] Asahi | Theory of Omega-d Modes[END_REF] Conclusions follow in Sec. VI.

For simplicity we shall use a notation based on a continuous formulation and let to the reader the change to a fully discretized operator which is used in the numerical application.

II. THE ALPHA EIGENVALUE EQUATIONS

Our starting point are the time-dependent kinetic diffusion equations coupled to the precursors equations in a heterogeneous domain D:

             1 v ∂ t + B pr ψ = p χ p λ p C p , (∂ t + λ p )C p = F p ψ. (1a) (1b)
Here ψ(x, t) is the scalar flux, x = (r, E) stands for the phase space variables, C p (r, t) is the concentration for precursor p and B pr = L -P pr is the diffusion kinetic operator with prompt neutron production

L = -∇ • D∇ + Σ -H, (2)
where D(x, t) is the diffusion coefficient, Σ(x, t) is the total macroscopic cross section, H stands for the scattering operator and P pr = i χ i (1 -β i )F i is the prompt fission operator with the sum in i over the fissile isotopes. Also in Eqs. (1b) and [START_REF] Foulke | Application of the natural mode approximation to space-time reactor problems[END_REF], F p = i β i→p F i and F i are the production operators for precursor p and neutrons from isotope i, respectively. The latter is defined as F i ψ = νΣ f,i ψdE. We note N p and N i the number of precursors and fissile isotopes, respectively, and take the prompt and delayed fission spectra, χ i and χ p , to be normalized so that χdE = 1.

By introducing the reduced scalar product over the energy E,

f, g = f gdE, (3)
where f denotes the complex conjugate of f , the two latter formulas can be written as F i ψ = νΣ f,i , ψ and 1, χ = 1. Finally, λ p and β i→p are the decay constant for precursor p and its fission yield from isotope i, respectively, and β i = p β i→p is the total precursor yield for isotope i. We shall follow the usual convention of ordering the λ p in increasing values. The kinetic equations in (1) are supplemented with initial conditions for both ψ and C p and with boundary conditions for ψ. The cross sections, external source and boundary conditions may change in time. But at any time t = t * one might consider a new kinetic problem with the initial and boundary conditions at time t * and such that for t > t * the cross sections and boundary conditions remain constant and the external source vanishes. This defines a kinetic problem with constant cross sections and boundary conditions and with no external source which will evolve in time according solely to its initial conditions at t * .

The alpha eigenvalue equations for this t * -kinetic problem are obtained by introducing an exponential time behavior in (1), ψ(x, t) ∼ e αt ψ α (x), C p (r, t) ∼ e αt C p,α (r), to obtain:

             α v + B pr ψ α = p χ p λ p C p,α , (α + λ p )C p,α = F p ψ α . (4a) (4b)
Here ψ α (x) and C p,α (r) are respectively the scalar flux and the concentration for precursor p of an α-mode, and all other quantities and operators are like those defined for the timedependent kinetic equation but with the cross sections evaluated at time t * . More generally, by introducing the vector function Ψ α = {ψ α , C α }, where C α has components {C p,α , α = 1, N p }, the set of alpha equations can be written in the form of a classical linear eigenvalue problem:

M Ψ α = α Ψ α , (5)
where the structure of matrix operator M, not explicited here, can be inferred from Eqs. [START_REF] Dugan | Cross-Section Homogenization for Reactivity-Induced Transient Calculations[END_REF].

When appropriate we shall use Hilbert spaces for a precise formulation of results. We note H X the Hilbert space of complex-valued functions defined over phase space X and satisfying the boundary conditions of problem (4). This space is endowed with the scalar product

( f, g) = X f (x)g(x)dx, (6)
where dx = drdE is the volume element in X. We also note H D the Hilbert space of complex-valued functions defined over the geometrical domain D with scalar product

( f, g) D = D f (r)g(r)dr.
Thus, we consider problem (5) in the Hilbert space of complexvalued functions H X × (H D ) N p with scalar product

(( Ψ, Ψ )) = (ψ, ψ ) + p (C p , C p ) D .
An useful form of the alpha-equations can be obtained by solving Eq. (4b) for the precursor concentrations:

C p,α = F p ψ α α + λ p , α ∈ C\Λ, (7)
where Λ = {-λ p , p = 1, N p }. Next, replacing this result in Eq. (4a) yields an expression for the alpha-eigenvalue problem in terms of only the flux:

α v + B ψ α = -K α ψ α , α ∈ C\Λ, (8)
where B = L -P,

P = i χ i F i (9)
is the production operator with the steady-state fission spectrum

χ i = (1 -β i )χ i + p β i→p χ p (10)
and

K α = p α α + λ p χ p F p (11)
comes from the delayed contribution to production. We consider this problem in H X . We note that Eq. ( 8) is a non linear eigenvalue problem for (α, ψ α ). Moreover, as we discuss next, except for very pathological cases, Λ is in the resolvent set of the linear operator in Eq. ( 5) and therefore Eqs. ( 4) and [START_REF] Gozani | The concept of reactivity and its application to kinetic measurements[END_REF] are equivalent.

Let us analyze in which conditions α = -λ p can be in the spectrum of Eqs. [START_REF] Dugan | Cross-Section Homogenization for Reactivity-Induced Transient Calculations[END_REF]. If this is true, then one must have from Eq. (4b)

F p ψ α = 0. ( 12
)
We note that one can write F p = hπ p where h(r) = (F p 1)(r) and π p = h -1 F p is a projector π p : H X → H D . Equivalently, constraint [START_REF] Sanchez | Duality, Green's functions and all that[END_REF] can be written as

ψ α ∈ ker(π p ),
where ker(π p ) is the null space of π p . For the other precursors one obtains

C p ,α = 1 λ p -λ p F p ψ α , p p, (13)
which yields a source equation for the flux:

Bψ α = χ p λ p C p,α , (14)
where

B = B - λ p v - p p λ p λ p -λ p χ p F p .
Hence, a solution of Eq. ( 14) with the constraint in ( 12) is a pair (ψ α , C p,α = Bψ α /(χ p λ p )) 2 iff 3 (i) 0 ψ α ∈ ker(π p ), (ii) Bψ α ∈ E p , where we have defined

E p = { f (x) = χ p (E)h(r); h(r) ∈ H D }. In other words -λ p is an alpha- eigenvalue iff B(ker(π p )/{0}) ∩ E p ∅.
Although this iff condition does not apparently require operator B to have an inverse, an equivalent formulation can be written as B-1 E p ∩ (ker(π p )/{0}) ∅. Thus, a solution might exists if operator B is degenerated or if it is invertible over E p ; the former condition seems unrealistic but the second looks plausible because v is large and Σ contains scattering, capture and fission and therefore dominates H. However, in both cases the solution of equation ( 14) must satisfy the orthogonality condition [START_REF] Sanchez | Duality, Green's functions and all that[END_REF], which most likely will not be feasible. Therefore, with the exception of unrealistic cross sections, the existence of an alpha eigenvalue equal to minus a precursor decay constant can be excluded.

Clearly the precedent conclusion applies to the case of the one-group diffusion equation. But we give here a counter example for which -λ p is an alpha eigenvalue. Consider the case when the diffusion coefficient does not depend on the energy, D = D(r), and take ψ α (x) = f (E)h(r) with π p f = 0 and h(r) an eigenfunction of the Laplacian equation -∇ • D∇h = γh which satisfies the boundary conditions of the alpha problem. Then Bψ α = h Bγ f , where Bγ equals B with the replacement -∇ • D∇ → γ. Therefore, by setting χ p = Bγ f /c, where c(r) = Bγ f dE, the pair (ψ α , C p,α = (c/λ p )h) is a solution of the alpha eigenvalue problem for α = -λ p . A problem with the definition of χ p is that it depends on space, but this can be avoided by considering a homogeneous problem. However, even then, one would have to "adjust" the cross sections, perhaps with nonphysical values, so as to ensure that χ p remains positive.

Perturbation formulation

Reactivity is a measure of the disequilibrium of a system and there are several ways to introduce a notion of reactivity in neutron kinetics [START_REF] Gozani | The concept of reactivity and its application to kinetic measurements[END_REF][START_REF] Cacuci | Eigenvalue-dependent neutron energy spectra: Definitions, analyses, and applications[END_REF], in particular that derived with the help of the adjoint α equation. [START_REF] Henry | The application of inhour modes to the description of non-separable reactor transients[END_REF] In this work we define a generalized reactivity as ρ λ = 1 -1/λ, where λ is any eigenvalue of the quasi-static eigenvalue equation obtained by setting α = 0 in Eqs. [START_REF] Gozani | The concept of reactivity and its application to kinetic measurements[END_REF]:

B λ ϕ λ = 0, (15)
2 For simplicity, we assume that χ p is strictly positive. 3

if and only if

where

B λ = L -(1/λ)P. ( 16
)
Here λ ∈ σ(L -1 P), where σ denotes the spectrum of the operator, is an eigenvalue of the quasi-static diffusion operator and P is the steady-state production operator in Eq. ( 9). The boundary condition for Eq. ( 15) is that of the parent alpha eigenvalue equations. We note E λ and N λ , respectively, the subspace of eigenfunctions of eigenvalue λ and its dimension. We recall that the eigenvalue λ = k with maximum absolute value in σ(L -1 P) is real, non degenerate and has a positive eigenfunction ϕ k (a physical flux).

In order to analyze the solutions of the alpha-eigenvalue equation it is advantageous to write Eq. ( 8) in terms of the quasi-static operator B λ in [START_REF] Planchard | Méthodes mathématiques en neutronique[END_REF] and cast the alpha-eigenvalue equation as

B λ ψ α = K λ,α ψ α , (17)
where λ ∈ σ p (L -1 P) and

K λ,α = ρ λ P -( α v + K α ). ( 18
)
The appeal of formulation [START_REF] Gozani | Subcritical reactor kinetics and reactivity measurements[END_REF] is that (i) it offers a perturbation expression for the alpha eigenvalue problem in terms of the better known steady-state eigenvalue problem, and (ii) the operator B λ on the left of the equal sign is singular and therefore the "source" term K λ,α ψ α must satisfy a solvability condition.

According to the Riesz-Schauder theory, [START_REF] Yosida | Functional analysis[END_REF] the solvability condition requires the source term to be in the orthogonal complement E ⊥ of the eigenspace E † λ associated to the eigenvalue λ ∈ σ((PL -1 ) † ) of the adjoint operator B † λ . 4 More explicitly, (ϕ

† , K λ,α ψ α ) = 0 ∀ϕ † ∈ E † λ . (19)
Thus, the solvability condition results in a system of N λ onepoint-like "kinetic" alpha equations. For each function ϕ † ∈ E † λ we have an inhour-type equation:

ρ λ = α         Λ + p Γ p α + λ p         , (20)
where

Λ = (ϕ † , 1 v ψ α ) (ϕ † , Pψ α ) , Γ p = (ϕ † , χ p F p ψ α ) (ϕ † , Pψ α ) . ( 21
)
As an aside, we note that in Eqs. (21) one can write Γ p = i β i→p γ i,p with the effectiveness factors [START_REF] Keepin | Physics of nuclear kinetics[END_REF]

γ i,p = ϕ † k , χ p νΣ f,i , ψ α dr/(ϕ † k , Pψ α ) and (ϕ † k , Pψ α ) = i ϕ † k , χ i νΣ f,i
, ψ α dr, recovering thus the formula derived for an infinite homogeneous medium in Ref. [START_REF] Dall | Neutron spectrum kinetics in the infinite homogeneous reactor[END_REF]. Also, because of the projection by ϕ † the original eigenvalue α is dispersed among the N p + 1 roots of Eq. (20).

Under the solvability constraint (19), the general solution of ( 17) is an arbitrary multiple of

ψ α = ϕ λ + ψ ⊥ , (22)
where ϕ λ = n a n ϕ n is written using an arbitrary basis {ϕ n , n = 1, N λ } in E λ and ψ ⊥ ∈ E ⊥ . The latter condition entails

a n = (ϕ † n , ψ α ), (23)
where the

{ϕ † n , n = 1, N λ } is the dual basis in E † λ , i.e., such that (ϕ † n , ϕ n ) = δ nn .
We end this section with a comment regarding the effect of spatial symmetries on the structure of the in-hour equation. Consider a problem with a spatial involution, that is a symmetry s such that s 2 = 1. Then ψ α obeys the relation s 2 ψ α = ψ α and, since operator s 2 has eigenvalues 1 and -1, one can write ψ α as the sum of an even (+) and an odd (-) components ψ α,± = (s ± 1)ψ α with the property that drψ α,+ (r)ψ α,-(r) vanishes. Because both the steady-state transport equation and the alpha eigenvalue equation necessarily share the same geometry one concludes that the in-hour equation does not exist when ϕ † and ψ α have opposite symmetries.

III. SUMMARY OF RESULTS

The spectrum of the alpha mode equations is conditioned by the nature of these equations and by the degree of complexity accepted in the delayed neutron term. In the continuous case the eigenmodes, Ψ α (x), belong to a functional space of infinite dimensions, while in the discretized case they become vectors in a finite-dimensional space. On the other hand, the degree of complexity of the fission delayed contribution depends on the number of fissile isotopes N i and the number of precursors N p . The operator in Eq. ( 11) can be written as a sum of projectors with r-dependent coefficients:

K α = p α α + λ p a p π p , (24)
where a p (r) = F p χ p and π p = (1/a p )χ p F p is a projector over the subspace generated by the {χ p , p = 1, N p }. Therefore, for α {0∪Λ}, K α is a finite linear combination of projectors with finite range equal to the span of the functions {χ p , p = 1, N p }.

Moreover the null space of K α contains the set of all the functions that are orthogonal to the finite set

{νΣ f,i , i = 1, N i }.
Simplifications of the fission delayed contributions are typically introduced to facilitate the theoretical analysis of the alpha spectrum. A simplifying assumption can be applied to the delayed fission spectra, to the fission production or to both. The first consists of replacing all the delayed neutron spectra with an averaged one, χ p ∼ χ, which leads to the simplified form K α = χF α , where F α = p [α/(α+λ p )]F p . In the second one replaces the fission production for all fissile isotopes with a single one, νΣ f,i ∼ w i (r)νΣ f , where by convenience we assume i w i (r) = 1. This approximation yield the simplified form K α = χ α F , where now χ α = p [α/(α + λ p )]β p χ p with β p (r) = i w i (r)β i→p . Note that these approximations are automatically satisfied if one assumes one single precursor (χ = χ 1 ) or one single fissile isotope (w 1 = 1 and β p = β 1→p ). Finally, when both approximations are used one has the simplified operator

K α = f α (r)χF , (25)
where

f α (r) = p α α + λ p β p (r). ( 26
)
Note that for a single fissile isotope β p (r) = β p , while for a homogeneous medium β p = i w i β i→p . For both cases f α (r) is independent of r and is an analytical function of α in C\Λ.

There is a sizable number of theoretical studies of the alpha spectrum in transport theory but, to our knowledge, there are only three communications pertaining to diffusion theory. Porsching [START_REF] Porsching | On the spectrum of a matrix arising from a problem in reactor kinetics[END_REF] used algebraic techniques to analyze the spectrum of the alpha mode discretized diffusion equations in a homogeneous medium with a single fissile isotope. However his equations were based in a synthetic expression for the time-dependent flux and, in regard to the present standard alpha eigenvalue problem, Porsching's results apply only to the one-group diffusion equation. A generalization to the multigroup diffusion equation was carried out by Devooght [START_REF] Devooght | Spectrum of the multigroup-multipoint diffusion operator with delayed neutrons[END_REF] for the fully discretized diffusion equation in a heterogeneous medium with no upscattering and with a single fissile isotope. To simplify the delayed fission term, Devooght assumes that all the fissile neutron spectra are identical, i.e., not only that χ p ∼ χ but also χ 1 ∼ χ, which lead him to lump prompt and delayed fission in a single term and work with the function f (α) = 1f α , which is also analytical in C\Λ. This function was used in the general study of the spectrum. The analysis of the fast spectrum, based on a perturbation result for the case of degenerated operators discussed in Kato, [START_REF] Kato | Perturbation theory for linear operators[END_REF] was carried out for the case with no precursors, while a technique early introduced by Wings [START_REF] Wing | An introduction to transport theory[END_REF] was applied to the investigation of the degeneracy of the delayed spectrum. The result of the analysis is that the spectrum consists of N g + N p clusters of N r eigenvalues each, where N g and N r are the number of groups and the number of spatial components, respectively. The delayed spectrum has N p clusters of real eigenvalues with the p-th cluster confined to the open interval5

I p = (-λ p , -λ p-1 ) (27)
with the convention that -λ 0 = ∞. The only theoretical result for the one-group continuous diffusion equation in one-dimensional heterogeneous slab geometry was given by Asahi [START_REF] Asahi | Theory of Omega-d Modes[END_REF] a few years earlier than Devooght's paper but seemingly unknown to the latter. Asahi demonstrated that for an arbitrary geometry the eigenvalues were real and bounded and that the -λ p values were not in the spectrum. He gave also a more detailed analysis for the case of the one-dimensional slab, classifying the eigenfunctions by the number of nodes and showing that there are a countable infinite number of eigenvalues in each interval I p which accumulate at -λ p as well as in the open interval

I L N p = (-∞, -λ N p ). ( 28
)
The rightmost eigenvalue in each interval has an eigenfunction with positive flux (zero nodes) and the number of nodes of the eigenfunctions fluxes increases monotonously by one for the other eigenvalues in the interval. Finally, as shown in Appendix A, an analytical solution can be obtained for the one-group homogeneous case which confirms the structure of the prompt and delayed spectrum as well as the emergence of flux clustering. This result leads to a new derivation of Asahi's results for the one-group homogeneous finite slab (see Appendix B).

In the second part of this section we give general observations in the nature of the spectrum of the alpha eigenvalue problem derived from robust arguments and from numerical exploration of multigroup diffusion in a slab geometry. Our results match the known theoretical predictions in their restricted domain of validity. [START_REF] Porsching | On the spectrum of a matrix arising from a problem in reactor kinetics[END_REF][START_REF] Asahi | Theory of Omega-d Modes[END_REF][START_REF] Devooght | Spectrum of the multigroup-multipoint diffusion operator with delayed neutrons[END_REF] Our discretization was done for the one-dimensional piecewise heterogeneous slab geometry by writing the linear alpha eigenvalue equations ( 5) in its multigroup form and by introducing a numerical discretization for the spatial dependence of both ψ and C. We observe that the number of degrees of freedom (DOF) of the discretized system of equations is

N DOF = N r (N g + N p ). ( 29
)
A natural ordering is build in the kinetic equations. This ordering stems from the observation that most of the neutrons resulting from fission are prompt and that, therefore, the contribution of delayed neutrons to the overall neutron production in the kinetic equations can be considered as small. In addition, the frequencies of emission of delayed neutrons are also much smaller than the neutron collision frequency. These two conditions can be expressed as

β p
1 and λ p vΣ.

In the following we shall assume that all cross sections are O(1) and we shall define

β = O(β) and v = O(1/v).
Physical data tell us that β ∼ 10 -3 -10 -4 and v ∼ 10 -5 -10 -9 . This ordering has direct implications in the distribution of the alpha spectrum, which can be sorted out into two subsets: delayed and prompt eigenvalues.

Delayed spectrum

These modes appear in the range |α| ∼ O(λ p ), in N p clusters of N r eigenvalues each, for a total of N r N p eigenvalues (see Fig. 1). The eigenvalues in the p-th cluster gather on the right of -λ p :

-λ p < α p,N r < ... < α p,2 < α p < -λ p-1 ,
where α p,i is the i-th eigenvalue of the cluster in the interval I p . Also, we have written α p,1 simply as α p to indicate the special nature of the dominant mode in the cluster. While the other eigenvalues in the cluster are very close to -λ p and their values do not change much with the reactivity, the value of the dominant eigenvalue α p increases with the reactivity towards -λ p-1 . The delayed eigenvalues are simple and real negative with the exception of α 1 which becomes positive for positive reactivity. This latter eigenvalue admits a positive eigenvector (both fluxes and concentrations are positive). The remaining dominant eigenvalues (α p , p > 1) have positive fluxes but the concentrations change of sign.

(α) -λ N p . . . -λ p -λ p-1 . . . -λ 1 0 α 1 • •• •• •• •• • N r N r α p,r , r = N r , . . . , 2
α p
Fig. 1. Clusterings of eigenvalues approaching the negative values of the precursor's decay constants.

The p-th cluster can be analyzed using the expression

α + λ p = p (30)
and considering the eigenvalue problem for p → 0. The analysis shows that there is a leading concentration for each of the N p clusters which because of the special variation degenerates into N r independent distributions. Hence, this part of the spectrum can also be named "precursor" spectrum. We have also made the conjecture that the eigenmodes concentrations in each cluster are linearly independent and this has been supported by all our calculations. In order to analyze the clustering results we resort to an approach different from the ones adopted in the literature. We recognize that the behavior of the delayed spectrum is directly related to the changes undergone by the precursors concentrations and, therefore, replace Eq. ( 8) with an "equivalent" equation for the concentrations. We use Eq. ( 7) to write the delayed contribution in Eq. ( 8) in terms of the concentrations and by expliciting the ψ α in the resulting equation we write

ψ α = -(α/v + B) -1 α p χ p C p,α , α σ(-vB). (31)
Next operating on this equation with (α + λ p ) -1 F p yields a finite system of equations for the concentrations, which we choose to write:

(α + λ p )C p,α (r) = -α p (T pp ,α C p ,α )(r), α σ(-vB), (32)
where T pp ,α is an integral operator over E and r with kernel

t pp ,α (r → r) = dEνΣ f,p (x) dE χ p (E )g α (x → x), (33)
where νΣ f,p = i β i→p is the fission production term for precursor C p , x = (E, r) and g α (x → x) is the Green's function of operator (α/v + B) -1 . We note that Eq. ( 32) is not equivalent to Eq. ( 8) for the simple reason that the energy variable does not appear in the former equation. This can be stated more precisely for the discretized case where the number of DOF for Eq. ( 32) is N r N p whereas Eq. ( 8) has N r N g N p DOF. This shows that, except for the one-group case, the two equations are not equivalent. However, if α is an eigenvalue of Eq. ( 32) with eigenvector C α (r) and α σ(-vB), then α is also an eigenvalue of the full alpha problem in Eq. ( 8) with the eigenvector given by Eq. ( 31). Thus, we can safely use Eq. ( 32) to analyze the delayed alpha spectrum. We introduce the assumption in Eq. ( 30) and use asymptotic analysis to study the behavior of the spectrum near -λ p . A full analysis would consist of expanding the concentrations and the eigenvalue as a sum of increasing powers of p , but here we use a rough approximation and keep only the leading terms in the equation. However, the discussion is awkward when working with Eq. (32) and we shall, instead, take a step back and work with Eqs. ( 7) and [START_REF] Gozani | The concept of reactivity and its application to kinetic measurements[END_REF]. We notice first that, consistently with Eq. ( 7), we have C p,α = F p ψ α / p ∼ -1 p and, for p p, C p ,α = F p ψ α /(λ p -λ p + p ) ∼ O(1). Hence, we can neglect the contribution from the p p precursors in the sum in the right-hand-side of Eq. ( 8). Following the derivation of Eq. ( 32), this implies that the system of equations in (32) uncouples and, in particular, yields a single equation for C p,α . To simplify our notation we write this equation as follows

p C p,α (r) = λ p (T α C p,α )(r), α σ(-vB) (34)
with p playing now the role of eigenvalue. Still, this last equation is difficult to analyze and we should simplify the operator by keeping only the leading contributions. Going back to Eq. ( 8), we observe that because F p ∼ β the order of the right-hand term in this equation is -1 = β / p instead of -1 p . So, in order to have 1, we have to assume p β . Having made sure that the right-hand-side of Eq. ( 8) is ∼ -1 , we observe that on the left-hand-side (α/v + Σ -H -B)ψ ∼ O(1) so we have α/v + B ∼ -∇ • D(E, r)∇. The spectrum of this operator consists of a countable infinite number of real positive eigenvalues growing towards +∞; these eigenavlues have finite multiplicity and the corresponding eigenfunctions form a complete basis in L 2 . [START_REF] Planchard | Méthodes mathématiques en neutronique[END_REF] It follows that he inverse of this operator is compact and definite positive and its Green's function is of the form g α (E, r → r), where E appears only as a parameter. This implies that operator T is also compact, and therefore it has real eigenvalues, but it is not necessarily definite positive. However, it is easily verified that for a homogeneous medium as well as in the one-group case T is definite positive and therefore its eigenvalues are positive and accumulate at zero. This proves that the eigenvalues that the eigenvalues α which converge towards -λ p do this from the right (α > -λ p ).

Prompt spectrum

The remaining part of the spectrum lies on the complex plane (α) < -λ N p and consists of simple eigenvalues that are real or complex conjugated in the range |α| ∼ O(vΣ); the exception is the first dominant eigenvalue α N p +1 which is real and, with increasing reactivity tends to -λ N p . The eigenvectors have negligible concentration values and they are dominated by energy modes.

To analyze these eigenvalues we use again an asymptotic approach by assuming that α ∼ O(v). Hence, one can neglect the contribution of the delayed neutrons and write the alpha eigenvalue equation by keeping only the leading terms as

( α v + B pr )ψ α ∼ 0. ( 35
)
This shows that the prompt spectrum is closely related to the spectrum of the diffusion equation with prompt fission, α ∈ σ(-vB pr ).

Flux clustering

We use the expression flux clustering to refer to a set of distinct eigenvalues which have very close fluxes. Flux clustering for the alpha eigenvalue transport problem was discussed earlier by Gozani [START_REF] Gozani | Subcritical reactor kinetics and reactivity measurements[END_REF] and by Henry. [START_REF] Henry | The application of inhour modes to the description of non-separable reactor transients[END_REF] Gozani applied his analysis also to diffusion theory, but the most interesting result was Henri's analysis of the two-group P1 transport equation, which is closely related to diffusion; by using the analytical solutions of this problem, Henri was able to give a satisfactory argument, albeit non rigorous, in support of flux clustering.

The flux components of the N p+1 dominant alpha modes can be shown to cluster about the dominant mode of the quasistatic equation ϕ k . This includes the dominant modes in the delayed spectrum (α p , p = 1, N p ) plus the dominant mode of the prompt spectrum (α N p +1 ). It can also expected that non-dominant delayed modes can exhibit some form of flux clustering, like sharing the the same number of spatial fluctuations in the slab case, and that this type of flux clustering might even apply to prompt modes. However, in our numerical calculations we have not observed yet any such clustering, except for the fact that for each energy group one might be able to identify N r modes that are associated, in the sense that they have a similar spatial behavior.

Our analysis of flux clustering is based on the perturbation approach discussed in Sec. 1, whose main results are the inhour equation and the flux perturbation formula, Eqs. (20) and Eq. ( 22), respectively. For simplicity, hereafter we shall only consider the case when λ is the dominant eigenvalue of the quasi-static diffusion operator, which we denote by k. This eigenvalue is simple, real positive and its eigenvector ϕ k is the so-called fundamental flux. Hence, the system of in-hour equations in (20) simplifies to a single equation which we write as

ρ k = ω         Λ α + p Γ p,α ω + λ p         , (36)
where, to stress the fact that the equation has N p + 1 roots of which only one is the original alpha eigenvalue, we have replaced α with ω and we have used the lower index α to indicate that the coefficients Λ and Γ p depend on the eigenvector ψ α .

The in-hour equation can have complex roots and we refrain here from considering this case. But, as shown in Appendix C, if the coefficients Λ α and the Γ p,α in the in-hour equation have the same sign, which we take to be positive, then all the N p + 1 roots are real. Therefore, since the associated steady-state adjoint flux ϕ † k is positive, we shall consider only the case when the alpha eigenvector ψ is also positive. Note that this implies that the eigenvalue α is real. This is in fact the more clear-cut and best understood example of flux clustering.

A substantial argument in support of flux clustering can be made by looking for alpha modes with fluxes close to the fundamental static flux ϕ k . We proceed by writing the flux perturbation formula, Eq. ( 22), in the form

ψ α = (1 + )ϕ k + ψ ⊥ , ( 37
)
where | | 1 and where we have adopted the normalizations ϕ k = ψ α . Next, we use this expression to compute the coefficients in the in-hour equation Eq. ( 36) and observe that to leading term in the coefficients depend only on the static flux ϕ k . The reason is that ψ ⊥ ∼ √ . As an aside, we point out that another reason for neglecting ψ ⊥ is that the coefficients depend on integrals of a product of a positive function times ψ α and that the impact of ψ ⊥ in the value of the integrals is further reduced because, as opposed to ϕ k , ψ ⊥ must necessarily change of sign.

Hence, to leading order in , we might conclude that Eq. ( 36) is independent of ψ α and also of small variations in . Since this equation has N p + 1 real roots we can safely assume that each one of these roots corresponds, to order , to an alpha eigenvalue. However, there are two conundrums regarding the precedent conclusion. The first one is that if the in-hour equation were to depend only on ϕ k , then this flux will be a solution of the equation, but the only case in which this might happen is for k = 1 which entails, of course, α = 0. The second problem arises from the fact that the in-hour equation is only a necessary condition for (α, ψ α ) to be an eigenpair of the alpha eigenvalue problem. Hence, the problem has to resolved with a simultaneous consideration of both the in-hour equation and the original alpha eigenvalue equation. For the present approach this involves showing that solutions of the form in Eq. ( 37) are truly solutions of the alpha eigenvalue problem. This is what, in a limited way, Henry did for the two-group P1 equation. For a more general approach may consist perhaps of using a direct perturbation analysis between the alpha eigenvalue equation and the equation satisfied by ϕ k . Finally, we observe that the leading eigenvalue α 1 dominates the alpha spectrum and consequently it has a positive flux. It remains to show, however, that the fluxes of the other N p dominant eigenvalues α p , p = 2, N p are also positive, a fact that is confirmed by previous speculations and by all our calculations. In view of formula (37), this could perhaps be demonstrated by a dominance argument.

IV. NUMERICAL RESULTS

We have used finite differences to discretize the multigroup diffusion alpha eigenvalue equation in 1D slab geometry as well as the eigenvalue equations for the quasistatic equation and its adjoint, and implemented three different criteria, all based in the L 2 norm ψ re f -ψ α 2 = g,r (ψ re f -ψ α ) 2 g,r < clus , to explore flux clustering about alpha modes (ψ re f = ψ α , ∀α ), about the modes of the quasistatic equation (ψ re f = ϕ λ , ∀λ) and by using perturbation analysis ψ ⊥ 2 , |1 -a| < clus , where the latter formula has been obtained from Eqs. ( 22) and ( 23) for the case of a non degenerated eigenvalue (N λ = 1). Calculations are performed in double precision by a Python module (v.2.7), using the numerical numpy library (v.1.11.0). The module solves the multi-group, (S N ), multi-region α and λ eigenvalue problems with reflection, periodic and void boundary conditions; the full set of eigenpairs is computed.

We have analyzed a large number of problems for the infinite homogeneous case as well as for some finite cases with 1, 8 and 281 energy groups and for increasing values of N r . Flux weighting was used to homogenize and collapse the cross sections for a PWR fuel assembly at different conditions (see Table I). Thus, the corresponding α-sets are close. The calculations where done with 8 groups of precursors and 6 fissile isotopes (N p = 8 and N i = 6). We present here a sample of the results.

Case Tfuel Tmod Boron conc. By increasing the number of regions the delayed eigenvalues cluster very close to the negative decay constants and are only clearly visible for the one-group calculation, as shown in Table IV. The table also shows that the dominant eigenvalue of the clustering, α p , does not change with the number of regions.

ρ nb. ( • C) ( • C) (PPM) (pcm
As we show next, this follows from the fact that an infinite homogeneous medium sustains an uniform flux solution, regardless of the spatial discretization. The equation describing the energy spectrum of the uniform solution results by setting to zero spatial derivatives in Eq. ( 8):

α v + Σ -H -P ψ α = -K α ψ α , α ∈ C\Λ,
This equation is independent of the spatial discretization, be it continuous or not, and therefore its solutions are present for any value of N r . Moreover, according with the analysis in Sec. III.3, these uniform solutions appear in a cluster which contains all the N p + 1 leading modes. Contrary to the few group cases, complex eigenvalues appear with 281 groups in the infinite homogeneous problem (i.e. one spatial cell and zero current at the boundary), and for different reactivities, see Fig. 2. The distribution of the prompt eigenvalues is in general weakly influenced by the reactivity.

V. COMPARISON WITH ASAHI'S RESULTS

Asahi applied a continuous space description to the onegroup, one precursor problem and computed results for a symmetric reflected critical slab reactor with a central core of half width a = 10 cm surrounded by two reflectors of width b = 10 cm. [START_REF] Asahi | Theory of Omega-d Modes[END_REF] He computed the first 15 eigenvalues in the delayed spectrum and the first 16 in the prompt spectrum by progressively increasing the number of zeros of the associated fluxes. Asahi presented his results with 6 and 3 digits, respectively, for the delayed and prompt eigenvalues.

We have used our numerical approach to recalculate Asahi's reactor by subdividing each cm into M equal regions, for a total of N r = 40 × M spatial DOF, and by gradually increasing M until convergence. For M = 32 the corresponding delayed and prompt eigenvalues (15+16), out of a total of N DOF = 2560 eigenvalues, converged to 6 digits and their associated fluxes exhibit the same number of zeros as those from Asahi's. Table V shows our 31 eigenvalues and the relative difference in % to the Asahi's values. Note that there are large differences in the prompt eigenvalues. We think that some of Asahi's values lacked precision. Our results prove that, as the number of DOF increases, our numerical technique approaches the more and more closely the full spectrum of the operator.

N r = 1 N r = 2 N r = 4 -1.3821081-3 -1.3821081-3 -1.3821081-3 -1.2466691-2 -1.2466692-2 -1.2466698-2 -1.2466699-2 -1.4012729-2 -1.4012729-2 -1.4012729-2 -2.

VI. CONCLUSIONS

We have presented theoretical arguments about the alpha eigenvalue problem in multigroup neutron diffusion, in support of the general behavior observed in our calculations including the distribution of the prompt and delayed spectrum as well as flux clustering. Our analysis and numerical results confirm the few theoretical and numerical results found in the literature regarding the diffusion problem.

We have used a perturbation formulation, based on the steady state eigenvalue problem of the effective multiplication factor, to derive a generalized form of the in-hour equation. As opposed to previous derivations of the in-hour equation, which have been derived by projection of the original alpha equation, in our formulation this equation arises as a necessary and sufficient condition for the existence of the alpha eigenvalueeigenvector pair.

We have also conducted many numerical simulations for one-dimensional homogeneous and heterogeneous configurations in multigroup diffusion theory. In all cases, our numerical results supports the theoretical predictions. The characterization of the alpha modes in multigroup neutron transport is deferred to a future study.

APPENDIX A: ANALYSIS OF THE ONE-GROUP, HO-MOGENEOUS CASE

The analysis for the one-group case is made especially simple because χ i = χ = 1 and one can add all fissile isotopes into a single one with νΣ f = i νΣ f,i . Moreover, for a homogeneous medium we can also define β p = i β i→p νΣ f,i /νΣ f and cast Eq. ( 8) into the simplified form

f (α)ψ = -Bψ, (38)
where

f (α) = α v + p β p α α + λ p νΣ f , ( 39
)
and B is the diffusion operator with steady-state neutron fission spectrum (unity in the present case). Finally, using the boundary conditions one can construct the eigenfunctions of problem Bψ ζ = ζψ ζ and then any (α, ψ ζ ), where α is a root of

f (α) = ζ (40)
is an alpha eigenpair. Consider a homogeneous slab of thickness a. For the case with zero flux boundary conditions one has the eigenfunctions ψ(z) ∼ sin(B n z) and the eigenvalues

ζ n = νΣ f -Σ a -DB 2 n ,
where B n = (n + 1)π/a. Note also that

ζ n = -ρ n νΣ f ,
where ρ n = 1 -1/λ n is the generalized reactivity associated to the classical reactor eigenvalues λ n . There are N p + 1 roots of (40) and they split into N p + 1 different groups; new roots are added progressively to each group as n increases. We show next that the higher N p groups cluster on the right of each -λ p , but for the group at the left of -λ N p , the eigenvalues decrease without bounds towards -∞.

We examine first the behavior of the α's that are close to -λ p so that α + λ p = with | | 1 (this assumption will be proved by the following analysis). For these α's the critical condition (40) can be written to leading order in as -β p λ p νΣ f -1 = ζ n , where now plays the role of the eigenvalue α. The solution is n = -

β p λ p νΣ f ζ n → α n = -λ p 1 - β p νΣ f ζ n .
For increasing n, ζ n becomes negative and increases as (n + 1) 2 so that the factor between parentheses remains positive but diminishes very fast towards 1. Remember that we have assumed n to be small, which implies |ζ n | β p λ p νΣ f . Clearly, this can only happen if ζ n < 0.

Consider, finally, the case with α -λ N p Here we can neglect precursor contributions and write α = vζ, which entails that ζ is negative. Note that in this case the first ζ of this group is the smallest root of (40) for n = 0, which we know that is -λ N p . Here we have

α n -α n-1 = vD π a 2 (2n + 1).
This implies that the distance between two successive eigenvalues in this group increases with n and that lim n→∞ α n = -∞.

It follows that the eigenvalues in this group do not cluster.

APPENDIX B: ASAHI PROBLEM REVISITED

The simple technique in the previous Appendix can be extended to a general analysis of a piecewise homogeneous slab geometry. The basic idea is to use the Laplacian eigenfunctions in each homogeneous piece and impose interface and boundary conditions in order to obtain a general solution for the problem which is then replaced in Eq. (38). Unfortunately, as we show next, this leads to a complicated system of trigonometric equations for the determination of the alpha spectrum.

As an example of this analytical approach we consider the Asahi problem discussed in Sec. V. This problem contains a central core of width 2a and a bilateral reflector of width b. The domain is symmetric with respect to the middle plane. We shall treat only the right half of the domain for z ∈ [0, a + b], with vacuum boundary condition on the right at z = a + b and either the current or the flux vanishing on the left at z = 0 . The latter conditions are necessary because the initial problem is symmetric and therefore accepts even and odd solutions. Given the boundary conditions, the solutions are

ψ c (z) ∼ f c (B c z) for z ∈ [0, a] and ψ r (z) ∼ f r (B r (a + b -z)) for z ∈ [a, a + b].
The core function f c can be either a cosine or a sine, depending on whether the alpha eigenfunction is symmetric or antisymmetric, while f r is a sine. We impose continuity of the flux and current at the core-reflector interface z = a:

A c f c (B c a) = A r f r (B r b), A c D c B c f c (B c a) = -A r D r B r f r (B r b),
where the A's are two constants and f is the derivative of f . Dividing the second equation by the first we get the dispersion relation:

D c B c f c (B c a) = B r D r f r (B r b). (41)
Here f r is the cotangent function while f c is the tangent function for the symmetric solution or minus the cotangent for the antisymmetric one.

We have now a solution that satisfies boundary and interface conditions with B c and B r constrained by Eq. (41). Next, we our solution in the alpha eigenvalue Eq. (38

) to obtain          f (α) = νΣ f -Σ a,c -D c B 2 c , α v = -Σ a,r -D r B 2 r , (42a) (42b)
where f (α) has been defined in (39).

Equation (41) and the two latter equations form a system of nonlinear equations for the triplet (α, B c , B r ). We now proceed to eliminate B c and B r and obtain a final equation for α. From the two equations in (42) we obtain:

B 2 r = - α + vΣ a,r vD r , B 2 c = - f (α) + v(Σ a,c -νΣ f ) vD c , (43)
and use of these expressions in the dispersion relation results in a nonlinear equation F(α) = 0 for α. In practice one might proceed with a zero search for this nonlinear equation, where given a value of α ∈ R, one evaluates B r and B c via relations (43) and, then calculates F(α) as the difference between the left minus the right hand sides of Eq. (41).

Clearly, B r = 0 is a singular point at which the solution vanishes in the reflector and in the core. Note that B 2 r is negative for α > -vΣ a,r , in which case the cotangent f r becomes a hyperbolic cotangent and the reflector flux, which can be taken as positive, decreases from the interface to vanish at the boundary. This case, for which the dispersion relation is easy to solve at least graphically, corresponds to the delayed spectrum calculated by Asahi, whose eigenfunctions have no zeros in the reflector. In the opposite case, with B > 0 or, equivalently, with α < -vΣ a,r , the flux in the reflector can change of sign. This is the case of Asahi's prompt spectrum. The numerical analysis fo Asahi's problem is further simplified because of the simple form of f (α) for the case of a single precursor, but we do not pursue this development further here.

APPENDIX C: ROOTS OF THE GENERALIZED NORDHEIM EQUATION

Equation (20) is the condition for the existence of the eigenpair (α, ψ α ). We write this equation in the form ρ = ρ(ω):

ρ = ω         Λ + p Γ p ω + λ p         ,
where the Λ(ψ α ) and Γ p (ψ α ) are defined in Eq. (23). The above equation can be written as a polynomial in ω of order N p + 1 and therefore has as many roots but not all necessarily real. We assume that Λ and the Γ p 's are bounded and that Γ p Γ p for p p . To simplify the analysis we shall assume that these coefficients are real and positive. Under these assumption we show that for any ρ ∈ R the equation has exactly N p + 1 real distinct roots. As an aside, we recall that for ρ = ρ λ one of the roots must be ω = α.

Let R 0 = R\Λ be the real line excluding the N p points Λ = {-λ p , p = 1, N p }. We note that ρ(ω) is C ∞ over R 0 . It is clear that every one of these points is a vertical asymptote. Let ε > 0, then for ε → 0 we find that ρ(-λ p ± ε) → ∓λ p Γ p /ε. Therefore, ρ(ω) goes to -∞ at the right of the asymptote and to ∞ at its left. Moreover, the derivative

∂ ω ρ = Λ + p λ p Γ p (ω + λ p ) 2
shows that, ρ(ω) is a monotonously increasing function in R 0 . The roots of ρ(ω) = 0 give the locations at which the function changes of sign. Note that there is a single root in each interval (-λ p+1 , -λ p ) for p = 1, N p -1. Clearly ω = 0 is another root and, since the function monotonously increases, there are no roots for strictly positive ω. Next, we check the behavior as ω → ±∞ to find that, for large ω, ρ(ω) → ωΛ. Therefore there is a last root at some value smaller than -λ N p . The general aspect of the graph is depicted in Fig. 3. Let us discuss what happens when not all the coefficients are positive. Assume that only one, say Λ, is negative. This changes the behavior for ω → ±∞ and now we find that ρ(ω) → ∓∞ for ω → ±∞. This agrees with the fact that presently the derivative has two zeros. Hence, ρ(ω) has a positive minimum ρ min > 0 for ω < -λ N p and a positive maximum ρ max < ρ min for ω > -λ 1 . The number of roots is N p + 1 except for ρ ∈ (ρ max , ρ min ) where it is N p -1. Note that for ρ = ρ max or ρ = ρ min there is a double real root.

Finally, we consider the case when the coefficient Γ p is negative for one intermediary p value, while all other coefficients remain positive. Here a positive minimum ρ min > 0 appears in the interval (-λ p+1 , -λ p ) while a positive maximum ρ max < ρ min appears in the interval (-λ p , -λ p-1 ). The discussion is similar to the previous case. We left the analysis of more complicated cases, where more than one coefficient is negative, to the interested reader.

Fig. 2 .

 2 Fig. 2. Distribution on the complex plane of the prompt eigenvalues with 281 groups for the cases in TableI.

Fig. 3 .

 3 Fig. 3. Sketch of generalized Nordheim in-hour equation for N p = 6.

TABLE I .

 I Physical conditions of the cross section preparation in the PWR fuel assembly.For all cases in TableI, the neutron flux components of the ψ α corresponding to the dominant eigenvalues α p clustered within clus = 10 -2 . The dominant eigenvalues for cases 3 and 4 are reported in Tables. II and III, respectively.

	α s	rel. diff. (%)
	281G	8G	1G
	-1.4397459518e-03	+4.003 +2.054
	-1.4020587820e-02	+0.056 +0.171
	-3.7787664954e-02	-0.011	-0.133
	-8.2605497512e-02	+0.023	-0.063
	-2.0561160002e-01	-0.016 +0.056
	-5.9302199351e-01	-0.019	-0.016
	-1.4945472270e+00	-0.034	-0.011
	-3.4222886970e+00	-0.017 +0.014
	-2.6351110679e+02 -35.019	-0.972
	TABLE III. Comparisons of the fundamental and delayed
	reactor periods with 1, 8 and 281 energy groups for case 4
	(k ∞ = +0.999910).
						)
			1	286	286	0	+6880.6
			2	286	286	400	+3405.6
			3	286	286	600	+1690.3
			4	286	286	800	-9.3
			5	650	306	0	+4915.5
			6	650	306	400	+1624.2
			7	650	306	600	+0.1
			8	650	306	800	-1611.3
			9	1100	306	600	-1302.5
			10	1800	306	600	-3037.1
				α s		rel. diff. (%)
				281G		8G	1G
			+5.8029093620e+02 -25.514 +0.214
				-1.2552000351e-02	-0.033 +0.002
				-2.9640277111e-02	-0.235	-0.074
				-4.3562900973e-02	-0.139	-0.005
				-1.4048345489e-01	-0.282	-0.030
				-3.2792305224e-01	-0.680	-0.118
				-7.0546045788e-01	-0.379	-0.054
				-1.7081304062e+00	-0.301	-0.041
				-3.6271274472e+00	-0.139	-0.027
			TABLE II. Comparisons of the fundamental and delayed re-
			actor periods with 1, 8 and 281 energy groups for case 3
			(k ∞ = +1.017200).

Table I .

 I

	α 1,(•)	r.d. (%)	α 2,(•)	r.d. (%)
	-0.000009		--6.593 × 10 1	+0.046
	-0.095135	-3.307e-04 -1.340 × 10 3	+0.027
	-0.098363	-2.959e-04 -3.957 × 10 3	-0.084
	-0.099205	-4.805e-04 -8.107 × 10 3	-0.040
	-0.099536	-2.496e-04 -1.382 × 10 4	+0.153
	-0.099697 +3.656e-04 -2.098 × 10 4	-0.076
	-0.099788	-3.637e-04 -2.790 × 10 4 +14.352
	-0.099843	-2.955e-05 -2.961 × 10 4	+0.033
	-0.099879 +2.620e-04 -3.205 × 10 4	-0.160
	-0.099904 +3.194e-04 -4.114 × 10 4 +15.417
	-0.099922 +3.322e-04 -5.168 × 10 4	-0.037
	-0.099936	-2.921e-04 -5.736 × 10 4	-0.071
	-0.099946	-9.027e-05 -6.038 × 10 4	+0.962
	-0.099954	-1.342e-04 -6.950 × 10 4	-5.892
	-0.099960 +1.894e-04 -8.366 × 10 4	-0.044
			-9.781 × 10 4	+0.006
	TABLE V. Eigenvalues with the relative difference (%) to
	Asahi's results: delayed spectrum (left), prompt spectrum
	(right).

It is understood that the adjoint operator is provided with appropriate "adjoint" boundary conditions[START_REF] Sanchez | Duality, Green's functions and all that[END_REF].

In reality, the result obtained by Devooght is that the -λ p also belongs to the p-th clusters and that, therefore, the precursor clusters contain N r + 1 eigenvalues each. However, he recognized that this pathological result occurred because he had artificially replaced the precursor concentrations C p (r) with the energy dependent functions C ♣ p (r, E) = χ p (E)C p (r) resulting in the associated eigenvalue problem having N r N p (N g -1) more eigenvalues than the original one


K Herbreteau
email: kevin.herbreteau@cea.fr

N Marie

F Bertrand

J-M Seiler

P Rubiolo
email: pablo.rubiolo@lpsc.in2p3.fr

FAST-RUNNING TOOL DEDICATED TO POWER EXCURSION SIMULATIONS IN SODIUM-COOLED FAST REACTOR

Keywords: Sodium-cooled Fast Reactor, Transient OverPower, Severe Accident, Fast-running tool

Within the framework of the Generation IV Sodium-cooled Fast Reactor (SFR) within the CEA (French Commissariat à l'Energie Atomique et aux Energies Alternatives) is involved, the innovative reactor designs under severe accidents conditions have to be assessed. Such accidents have been mainly simulated with mechanistic calculation tools (such as SAS and SIMMER). As a complement to these codes, which provide reference accidental transients calculations, a new physico-statistical approach is currently developed at CEA; its final objective being to derive the variability of the main results of interest to quantify the safety margin. This approach involves fast-running tools to simulate extended accident sequences, by coupling models of the main physical phenomena with advanced statistical analysis techniques. They enable to perform a large number of simulations in a reasonable computational time and to describe all the possible scenario progressions of the hypothetical accidents. In this context, this paper presents a physical tool (numerical models and result's assessment) dedicated to the simulation of the beginning of the primary phase of the Unprotected Transient OverPower accidents (i.e. before large pin degradation). At the beginning of this primary phase, the increase of nuclear power induces a strong temperature gradient in the fuel pellets leading to specific mechanical behaviours, such as swelling and thermal expansion, before their meltdown. The fuel pin thermal evolution during slow power increases, such as control rod withdrawal accidents, and fast power increases have been performed. Validation on some CABRI experiments was carried out and focused on the axial distribution of melting points in the fuel column, the molten fraction at the end of the transient and the coolant temperature evolution during the transient. The results are consistent with experimental data and an application case on the CFV core of ASTRID was performed in order to give tracks for the safety of the reactor.

INTRODUCTION

The current main objective of the French Generation IV project is to design a new reactor based on improved technologies in terms of safety and reliability at an industrial scale. Among other concepts, the Sodium-cooled Fast Reactor (SFR) has been selected for its ability to secure the nuclear fuel resources and to manage radioactive waste by minor actinides transmutation. A major innovation of the new SFR French concept concerns the core, which is featured by a very low (even negative) reactivity variation caused by a potential sodium voiding. More specifically, this feature has a strong impact on the severe accident sequences especially regarding the material heating-up and on the global core power evolution [START_REF] Bertrand | Comparison of the behavior of two core designs for ASTRID in case of severe accidents[END_REF]. In this framework of safety studies devoted to the Unprotected Transient OverPower (UTOP), a physical tool, called OCARINa (Outil de Calcul Analytique Rapide pour les Insertions de réactivité dans un réacteur à neutrons rapides refroidi au sodium), is being developed. As a complement to reference transients which are usually simulated with mechanistic codes such as SIMMER [START_REF] Kondo | SIMMER-III : computer program for LMFR core disruptive accident analysis[END_REF], physical tools, coupled to advanced statistical technics, enable to get the variability of the main results of interest for the safety. This general approach combining a mechanistic code and evaluation tools has already been conducted for another accidental initiator family (USAF -Unprotected SubAssembly Fault [START_REF] Marie | Physico-statistical approach to assess the core damage variability due to a Total Instantaneous Blockage of SFR fuel sub-assembly[END_REF] and ULOF -Unprotected Loss Of Flow [START_REF] Droin | Modélisation du transitoire ULOF[END_REF]).

Up to now, there is no tool dedicated to the simulation of the UTOP transient for heterogeneous cores. The SAS code is used to describe the primary phase (i.e. before large pin degradation) of a severe accident but it cannot handle two molten cavities in the fuel pin for heterogeneous fuel pin and is not able to describe the degradation phase (i.e. after large pin degradation). Moreover, the SIMMER code is currently used to describe the degradation phase but the mechanical models are not able to describe the primary phase. As an example, the fuel axial expansion is not modeled whereas it has a strong negative reactivity effect. These limitations prove that the development of this physical tool fills a gap in severe accident studies.

A focus on the phenomena occurring during UTOP transients are described in section 2. They are modelled in accordance to the level of details required to catch all the decisive phenomena. The physical models used so far (thermal models) are presented in section 3. Their validation against some CABRI experiments and an application on the CFV core are then presented in the same section.

UTOP TRANSIENT

Core design is obviously driven by performance and safety objectives. An important design and safety purpose is the prevention of a power excursion in case of sodium voiding. The core concept (called CFV) is an axial heterogeneous core of 1500MWth on the contrary to more classical homogeneous cores used in former SFR. The low sodium void effect of the CFV core mainly results from the presence of a sodium plenum above the fissile zones, combined to the presence of a fertile plate in the inner zone of the core encompassed by two fissile zones (displayed in Figure 1). The larger height of the outer fissile zone also enables the void reactivity effect to be lowered as well, due to neutron leakage enhancement.

Different initators

The UTOP is characterized, in the first place, by an external reactivity insertion leading to an increase of the nuclear power. The increase could be local or global with a specific kinetic (slow or fast power variation). Four types of UTOP, according to different initiators, are usually considered:

• Inadvertent Control-rod Withdrawal (IRW)
It is the most common transients occurring during reactor operation as they manage the regulation of the core power. In a fast neutron reactor, in such accident, a linear power increase of about 1-3 %PN/s [W/s] (PN stands for the nuclear power on normal conditions) in pins located near the withdrawn control rod could be generated.

• Compaction of the core

The reactivity inserted is assumed to have a sinusoidal shape due to a succession of compaction (positive reactivity insertion) and stacking (negative reactivity insertion).

• Fall of the core

The fall of the core is assumed to be radially homogenous; all the control-rods are stuck at their mechanisms. The reactivity inserted is therefore equal to the neutronic weight of the control-rods.

• Gas ingress into the core

This event may occur due to the entrainment of a part of the cover-gas plenum by the primary pumps that could push the gas into the core. The associated reactivity insertion depends on the size of the bubble and on the location of the bubble. The CFV core has a global negative void effect but it is positive close to the fissile zones.

The aim of the tool is to describe all of those UTOP transients on a reactor scale. Currently, powerimposed transients are modeled considering a single pin.

Overall phenomenology

Considering any of the previous postulated event, the multi-physic (thermic, mechanic, neutronic) phenomenology remains the same (Figure 2). Due to the fast kinetic of UTOP transient (except powerramp transients), some physical phenomena have not the time to settle and there is no need to model them such as axial thermal conduction (see section 3).

The main phenomena up to clad failure can be summarized as the following:

1. Reactivity insertion: the core becomes super-critical; 2. The nuclear power increases; 3. The fuel temperature increases, the Doppler effect tends to slow down the power increase; 4. If the Doppler effect isn't efficient enough to lower the reactivity down to zero, the fuel could melt; 5. The fuel melts in the center of the pellet: a molten cavity is created; 6. The fuel phase change combined with thermomechanical phenomena induces significant mechanical strain in the fuel pin clad;

7. Fuel pin clad failure; 8. Molten fuel is ejected out of the pin into the coolant increased by the eventual fission gas present.

Up to now at the current development stage of the tool, only a simplified sequence is modeled. No neutronic coupling is taking into account: the power variation is imposed during the transient.

Few words about CABRI

The CABRI core is a 25MW (Mega-Watt) thermal nuclear reactor operating with LWR-like fuel. The test channel is in the center of the core. The power transients are triggered with transient rods filled with gaseous Helium-3. Fast depressurization of these rods leads to a reactivity insertion in the core and a power transient up to 20GW (Giga-Watt).

PHYSICAL MODELS

The UTOP transient is characterized by fuel heating-up to its possible melting. Thus, a two-dimensional meshing of the pin is considered the most reasonable for an accurate simplified simulation in order to take into account the evolution if the fuel property versus the temperature which have been evidenced having an influence on the transient, especially through the thermal conductivity of the fuel [START_REF] Philipponneau | Thermal conductivity of (U,Pu)O2-x mixed oxide fuel[END_REF]. The sodium flow is considered as unidirectional from the bottom (inlet) to the top (outlet). Variables related to sodium and to the wrapper are space-averaged over each axial mesh as illustrated in Figure 3.

One or several subassemblies of the reactor core are modeled by a representative subassembly (noted SA) including the pins, its sodium coolant flow and its surrounding hexagonal wrapper (as done in SAS-SFR [START_REF] Lemasson | Simulation with SAS-SFR of a ULOF transient on ASTRID-like core and analysis of molten clad relocation dynamics in heterogeneous subassemblies with SAS-SFR[END_REF]). For example, for a single CABRI pin geometry, the pin modeled is presented in Figure 4.

Thermal modeling

The knowledge of the temperature in the fuel, the clad and the coolant is a mandatory step for the transient modeling. In fact, the temperature has an impact on the mechanical evolution of the fuel (thermal expansion, phase change, strains, creeps etc.) and on the nuclear power evolution via the Doppler effect. A finite difference method is used with an implicit coupling between the temperature of the sodium and the temperature in the fuel/clad [START_REF] Latrobe | A comparison of some implicit finite difference schemes used in flow boiling analysis[END_REF]. The coupling is made with the flux exchanged between the sodium and the clad. Many UTOP transients are very fast and require an unconditional stable numerical scheme. That's why a completely implicit numerical scheme has been developed. An explicit numerical scheme has been tested (the coupling is made with the flux exchanged between the sodium and the clad at the previous time step) and was not stable for fast transients.

Equations

The energy balances used are:

• One-dimensional radial heat transfer in the fuel:

�ρC p � f ∂T f ∂t = div�k f grad ����������⃗ T f � + q ′′′ (1)
• One-dimensional radial heat transfer in the clad:

�ρC p � cd ∂T cd ∂t = div�k cd grad ����������⃗ T cd � (2)
• One-dimensional axial energy balance in the coolant (liquid phase only):

∂ρ c H c ∂t + ∂G c H c ∂z = 2π S c (r cd φ cd + r w φ w ) (3)
• Energy balance of the wrapper:

�ρC p V� w ∂T w ∂t = -S w φ w (4)
Where ρ, C p , k, T stand respectively for density [kg/m 3 ], specific heat [J/kg/K], thermal conductivity [W/m/K] and temperature [K] of the fuel (f), the clad (cd), the coolant (c) and the wrapper (w).

In equation ( 1 Due to the kinetic of the transient and regarding the characteristic time of axial conduction in the fuel and the clad defined by:

τ axial conduction = ρC p k L 2
Where L stands for a characteristic length of the medium (two pellets for the fuel, the same value for the clad), it leads to a characteristic time of axial conduction (~10 5 s) much higher than the transient time (~10 -3 s for fast transient to 100 s for power ramps). That is why the axial conduction in the fuel and the clad are neglected.

Boundary conditions

An axisymmetric condition and an adiabatic condition beyond the wrapper are taken as boundary conditions.

The heat exchange between the fuel and the clad through the gap is described by a global heat exchange coefficient h gap [W/m 2 /K]. Concerning the CABRI tests, its value depends on the kinetic of the transient. For slow transient, the value h gap = 10 4 W/m 2 /K is taken and it corresponds to an average value determined by a GERMINAL calculation [START_REF] Melis | GERMINAL -A computer code for predicting fuel pin behavior[END_REF]. For fast transients, the value h gap = ∞ (in practice 10 9 W/m 2 /K) is taken, corresponding to a perfect solid contact between the pellet and the clad. It is justified by the resulting fast increase of the pellet temperatures which make the solid contact quasi-instantaneous by thermal expansion and phase change.

Finally, the heat exchange between the clad and the sodium is given by the Lyon-Martinelli [START_REF] Agosti | Heat Transfer Correlations for Liquid Metal Cooled Fast Reactors -Short Handbook[END_REF] transfer correlation: Nu = 7 + 0.025Pe 0.8

Where Nu and Pe stand for Nusselt and Peclet numbers.

Numerical resolution and analysis

After discretization, the energy balance equations and the boundary conditions provide a tridiagonal linear system per axial mesh:

MT = v
With T stands for the unknown vector containing the temperature of every node in an axial mesh (fuel, clad, coolant and wrapper), v stands for a known vector and M is a tridiagonal matrix. The inversion of the matrix is assured with the Thomas algorithm [START_REF] Thomas | Elliptic Problems in Linear Differential Equations over a Network[END_REF]. Moreover, in order to provide an accurate result at each matrix inversion, a numerical analysis was performed focusing on:

• the condition number к(M) defined, under the Frobenius norm [START_REF] Golub | Matrix computation, 3 rd edition[END_REF], by the ratio between the highest and the lowest eigen values of the matrix. For a well-conditioned matrix, к(M) = 1;

• the residual error, defined by r = |v -MT|, must be negligible.

The block diagram of the thermal resolution is shown in Figure 5. It begins with a steady state calculation followed by the transient calculation. The thermal properties are calculated with the temperature at the previous time step because at the current time step, the temperature in unknown. That's why a convergence loop is used to keep the implicit numerical scheme even on solid properties.

Figure 5. Block diagram of the thermal resolution. 3.2 Validation of the thermal modeling

This numerical scheme was testes on six CABRI tests [START_REF] Cranga | Transient material behaviour in CABRI-1 experiment failure under fully and semi-restrained fuel pin conditions[END_REF], [START_REF] Haessler | The CABRI-2 Programme-Overview on Results[END_REF], [START_REF] Fukano | Fuel pin behavior under slow ramp-type transientoverpower conditions in the CABRI-FAST experiments[END_REF]. Four fast transients: E5, E7, LT2, AGS0 and two slow transients: PF1, E9 (Table I). All the results are summarized in the same Table I. Some results obtained for fast transients (E5) and one slow transient (PF1) are shown in Figure 6 and7. The unvalidated results (Table I) are also shown in Figure 8. Globally, the thermal modeling gives good agreements with experimental data inside the uncertainty ranges so the model is considered validated regarding the objectives of the physical tool. It shows that the heat-up of all the medium during a power excursion can be accurately described with a simplified model. However, the differences with experimental data for the E9-test can be explained by the lack of mechanical modeling during the transient. Indeed, the evolution of the pellet-clad exchange coefficient is not yet modeled in the physical tool (see section 3.1.2). The average value taken for the whole transient must be better for the PF1 test than the E9 test, which could explain why good results are obtained for the PF1 test. Anyhow, the pellet-clad exchange coefficient evolution during the transient constitutes one of the following steps of modeling.

Pin burn-up (at%)

5 4, 6 12,4 2,9 6,4 4,6 Pin design OPHELIE-6 (annular) OPHELIE-6 (annular) QUASAR (solid) RIG-2 (solid) SCARABIX (annular) OPHELIE-6 (annular)
Condition number к ~ 1 ~ 1 ~ 1 ~ 1 < 5 < 4

Residual error r

< 2,5.10 -4 < 2,5.10 -4 < 2,5.10 -4 < 2,5.10 -4 < 6.10 -8 < 6.10 -8

Coolant

Application on the CFV design

UTOP power-imposed simulations have been performed focusing on an inner CFV core subassembly, described in Figure 9.

The idea is to suppose that a CFV inner core pin was inserted into the CABRI core. The same power evolution as E5 and E7 transients (Figure 10) are imposed. The linear power distribution at steady state is exposed on Figure 11. Due to the lack of neutronic and mechanic modeling, it is assumed in a first approximation that there is no modification of the axial form factor of the neutron flux during the transient like in CABRI. The linear power distribution magnitude is updated during the transient depending on the power evolution. At each time step, the linear power q ′ (z, t) is equal to:

q ′ (z, t) = q ′ (z, t 0 ) P relative (t)

Where q ′ (z, t 0 ) stands for the linear power at the beginning of the transient and P relative (t) stands for the relative power evolution (Figure 11).

It has to be noted that the power of the CABRI core was set regarding the energy per unit of mass injected into the fuel. The power integrated of a CFV pin is lower than pins used for E5 and E7 transients due to the fertile zones. So the energy injected into the CFV pin is lower than the energy injected into the CABRI pins and the results can't be directly compared. Regarding the heat exchange coefficient between the fuel and the clad through the gap, a constant value was taken for the fertile zones and another one for the fissile zones. Those values were provided by another GERMINAL calculation [START_REF] Melis | GERMINAL -A computer code for predicting fuel pin behavior[END_REF].

The results obtained for the sodium temperature at the Top of the Fuel Column (TFC) show (Figure 12) that for E5-like transient, there is a good margin regarding sodium boiling (around 250K). Also for E7like transient, there is a margin regarding sodium boiling but a smaller one (around 50K). Due to the lower energy injected into the pin, it is not guarantee that sodium boiling do not appear for the same energy injected as the CABRI tests. Now focusing on the molten radius at the power peak (Figure13), the E5-like leads to small amount of molten fuel in the upper fissile zone (less than 1% of the total volume of fuel). It means that it may not be axial motion of molten fuel through the pin according to previous experimental data [START_REF] Papin | Behavior of Fast Reactor Fuel During Transient and Accident Conditions[END_REF] and the central hole could be efficient enough to decrease strains on the clad and avoid its failure [START_REF]Le code PHYSURA, Objectifs, Description et Validation[END_REF].

On the contrary, E7-like transient leads to 30% of molten fuel in the pin with more than 80% of molten fuel in the upper fissile zone. This amount allows axial motions considering the experimental criteria of 40% of molten fuel used in the PHYSURA code [START_REF] Sato | Transient fuel pin behaviour and failure conditions in CABRI-2 in-pile tests[END_REF]. Due to the solid pellets constituting the fertile zones, the molten fuel cannot flow through the central hole by gravity flow or squirting effect [START_REF]Le code PHYSURA, Objectifs, Description et Validation[END_REF], preventing the release of strains on the clad. In that configuration, the clad will endure strong stresses and may break. Moreover, the melting pools surrounding the inner fertile zone could lead to the collapse of the fertile pellets. Those conclusions have to be confirmed simulating as well thermomechanical modeling of the fuel and the clad.

CONCLUSION

This paper deals with the preliminary development of a physical tool which will be used for safetyinformed design and to address core melting prevention and mitigation situations of sodium-cooled fast reactors. This tool will allow emphasizing main dominant phenomena and trends of significance for safety assessment. It is more specifically dedicated to the primary phase of an Unprotected Transient OverPower (i.e. before large pin degradation) in a heterogeneous core for which none adequate tool exist. This study focuses on the modeling, validation on CABRI experiments and an application on the CFV core of the heat-up of the fuel, clad, coolant and the wrapper up to large fuel fusion. The modeling of this tool is described. At that stage, it handles challenging implicit coupled resolution of:

• a one-dimensional transient heat equation in the solid materials of the SA (or single pin);

• a one-dimensional energy balance on the sodium flow along the SA height (one-phase flow modeling) and on the wrapper.

Validation against experiments has been performed using CABRI experimental database. It provided good results, especially regarding the coolant temperature evolution and the molten radius along the pin. It shows that the thermal evolution of the fuel, the clad and the coolant during a power excursion are accurately described with those analytical models. Two realistic power transients were then imposed on a CFV inner core and it gives tracks for the following steps of modeling such as sodium boiling and the collapse of the inner fertile zone.

The forthcoming work will focus on the thermomechanical modeling of the fuel and the clad. Indeed, most of mechanical phenomena occurring during a UTOP transient are thermal activated such as the thermal expansion of the fuel its expansion when it melts. Such phenomena will induce strong stresses on the clad and could lead to its failure. Its prediction (height of the failure and the failure time during the transient) therefore constitutes the main next step of modeling.

Figure 1 .

 1 Figure 1. Radial cut of CFV core geometry.

Figure 2 .

 2 Figure 2. Phenomenology tree of the UTOP transient up to clad failure.

Figure 3 .

 3 Figure 3. 2D meshing used for the modeling.

Figure 4 .

 4 Figure 4. Axial cut of a subassembly (or a single pin in the CABRI channel).

  ), 𝑞𝑞 ′′′ stands for the volumetric power source of the fuel [W/m 3 ]. At each time step, it is actualized depending on the imposed transient power. Then in equation (3), H c , G c , S c , r cd , φ cd , r w , φ w stand respectively for the coolant enthalpy [J/kg], the coolant flow rate [kg/m 2 /s], the flow cross section [m 2 ], the outer radius of the clad [m], the heat flux exchange between the clad and the coolant [W/m 2 ], the inner radius of the wrapper [m] and the heat flux exchange between the coolant and the wrapper [W/m 2 ]. Finally in equation (4), V w and S w stand for the control volume of wrapper [m 3 ] and the exchange surface between the coolant and the wrapper [m 2 ].

Figure 6 .

 6 Figure 6. Molten radius for E5 transients (Solidus front -solid, Liquidus front -dotted, Experimental double-phase front -squared). The fuel is in a double-phase state between the solidus and the liquidus front.

Figure 7 .

 7 Figure 7. Evolution of the Top Fuel Column (TFC) sodium temperature during PF1 transient (Experimental data -squared).

Figure 8 .

 8 Figure 8. Molten radius for E9 transients (Solidus front -solid, Liquidus front -Experimental double-phase front -squared). The fuel is in a double-phase state between the solidus and the liquidus front.

Figure 9 .

 9 Figure 9. CFV subassembly features [3].

Figure 10 .

 10 Figure 10. Relative power evolution imposed during the transient (E5 -solid, E7 -dotted).

Figure 11 .

 11 Figure 11. Linear power axial distribution at steady state (inner core, BFC -Bottom of the Fuel Column).

Figure 12 .

 12 Figure 12. Evolution of the Top Fuel Column (TFC) sodium temperature during the transient (E5 -solid, E7 -dotted, Boiling temperature = 1153K at 1 bar).

Figure 13 .

 13 Figure 13. Melting radius at the power peak (Solidus front E7 -solid, Liquidus front E7 -dotted, Solidus front E5 -circled, no liquidus front E5). The fuel is in a double-phase state between the solidus and the liquidus front.

Table I . Validation database used for the thermal modeling. (𝐏𝐏 𝐦𝐦𝐦𝐦𝐦𝐦 -Maximum power, 𝐏𝐏 𝟎𝟎 -Power at steady state, 𝐭𝐭 𝐏𝐏𝐦𝐦𝐦𝐦𝐦𝐦 -Time of the power peak) (𝐕𝐕 -Validated (see the section above), 𝐍𝐍𝐍𝐍 -Data not available, 𝐍𝐍𝐕𝐕 -Not yet validated)

 I

	Transient	E5	E7	LT2	AGS0	PF1	E9
	Type	TOP	TOP	TOP	TOP	Power ramp	Power ramp
	𝐏𝐏 𝐦𝐦𝐦𝐦𝐦𝐦 𝐏𝐏 𝟎𝟎	(𝐭𝐭 𝐏𝐏𝐦𝐦𝐦𝐦𝐦𝐦 )