Gabriele Grassi
email: gabriele.grassi@cea.fr

A Non-Linear Space-Angle Multi-Grid Acceleration for the Method of Characteristics in Unstructured Meshes

A new space-angle multi-grid technique has been developed to accelerate the free inner transport iterations based upon the Method of Characteristics (MOC). We present a two-level scheme: it consists of a ne level on which the MOC transport calculation is performed and a more coarsely discretised phase space in which a low-order problem is solved as an acceleration step. A ux-volume homogenisation technique is employed to dene the coarse-level cross sections. This entails the non-linearity of the scheme. Restriction and prolongation operators are dened between the two levels. After each ne transport iteration, a low-order transport problem is iteratively solved on the homogenised grid. A coarser angular representation is used within a MOC-like framework. We employ discontinuity factors to reconstruct the scalar incoming and outgoing currents on each region of the coarse discretisation. The solution of the above-mentioned low-order problem is used to correct the angular moments of the ux resulting from the previous free transport sweep. A complete description of the low-order operator and of the grid-to-grid transfer operators is given. A further application of the method to the acceleration of outer transport iterations is also presented. In order to test the eectiveness of our method, numerical tests for given benchmarks geometries have been performed. Results are discussed.

Introduction

Transport problems for nuclear reactors are multi-group k-eigenvalue problems, usually with upscattering. The multi-group formalism yields a system of one-group xed-source equations coupled via ssion and group-to-group transfer sources. The Method of Characteristics (MOC), implemented in the TDT solver for the APOLLO2 code [START_REF] Sanchez | APOLLO II : A User-Oriented, Portable, Modular Code for Multigroup Transport Assembly Calculations[END_REF], provides an iterative solution for the discrete ordinate formulation of the one-group linear transport equation [START_REF] Lewis | Computational Methods of Neutron Transport[END_REF]. In the MOC approach, the phase space is discretised as follows: the geometrical domain D is decomposed into a set of homogeneous regions {D i } on which a at-source approximation is done; a set of discrete angular directions and associated weights {Ω n , ω n } n=1,N is chosen. Then, a set of parallel trajectories is tracked for each direction. The method essentially consists of two main equations: a balance equation for the angular ux on each region i, and a propagation equation giving the angular ux leaving the region i in terms of the incoming angular ux and the internal source. The MOC has proved to be an advantageous tool for the solution of transport equation in unstructured meshes and, therefore, for realistic application to reactor analysis [START_REF] Sanchez | A Synthetical Acceleration for Two-Dimensional Characteristics Method in Unstructured Meshes[END_REF][START_REF] Halsall | CACTUS: A Characteristics solution to the neutron transport equations in complicated geometries[END_REF][START_REF] Goldberg | The Characteristics Method in General Geometry[END_REF].

However, since reactor transport problems are often characterised by collision-dominated regimes, any iterative solution of the transport equation -including MOC -generally requires a great number of iterations to converge [START_REF] Adams | Fast Iterative Methods for Discrete-Ordinates Particle Transport Calculations[END_REF]. This is why the ecient use of MOC needs acceleration techniques. The main idea of a multi-grid approach to speed up the convergence of transport iterations is to build a hierarchy of low-order problems, each solved in a more coarsely discretised phase space and consequently involving fewer unknowns. The solution of a given low-order problem provides a correction to accelerate the ux moments resulting from the previous iteration on the immediately higher level.

In this paper, a new non-linear space-angle multi-grid acceleration method dealing with full anisotropy of scattering is proposed. We extend and generalise some previous results for a spatial multi-grid method [START_REF] Grassi | A Non-linear Spatial Multigrid Acceleration for the Method ofr Acceleration of the Characteristics in Unstructured Meshes[END_REF]. To illustrate our acceleration, we focus upon a two-level scheme. This choice is not restrictive, as the algorithm can be extended to any number of levels with no signicant complications. The rst step is to dene a more coarsely discretised phase space. Coarsening aects both spatial and angular variable. For the spatial variable we act as follows. From the MOC spatial discretisation, that we will henceforth call ne grid, an associated coarse grid is obtained from the previous via a simple agglomeration procedure (see Fig. 1), in which each coarse mesh D I is obtained by fusing together adjacent ne meshes. Consequently, the measure of each macro-region is given by the sum of the measures of all the ne regions in it:

D I = ∪ i∈I D i ⇒ V I = i∈I V i . (1)
This requires a table of correspondences between the two spatial grids, and the implementation of a homogenisation technique in order to derive the material properties for the coarse meshes (typically the cross sections). In our case, the coarse-mesh total cross sections are obtained from a ux-volume homogenisation procedure. An immediate consequence of the latter is the non-linearity of the scheme.

For the angular variable, a low-order set of discrete angular directions and associated weights is employed: {Ω m , ω m } m=1,M with M < N .

(2)

On the coarse level, isotropic sources are used. The latter are derived from the high-order ones by an appropriate projection procedure which eliminates the angular dependence of the sources on the coarse level. Consequently, no table of correspondences between the two angular sets is needed. The scheme here presented involves a MOC-based transport sweep followed by an acceleration step consisting in iteratively solving a low-order transport problem for the zero-th ux moment. This problem is dened in order to provide, in a reference situation, the same result as the original high-order heterogeneous problem as far as averaged scalar properties are concerned. The MOC formalism is kept: outgoing angular uxes are computed using the same propagation equation as for the ne case. However, discontinuity factors (DFs) are introduced in the balance equation in order to reconstruct outgoing and incoming currents on each macro-region. Finally, a prolongation operator, using shape factors depending on the ne ux, is employed to reconstruct the ux moments on the ne grid. This method was initially meant to speed up convergence of inner iterations. Afterwards, it has been extended to the acceleration of outer iterations.

Both applications are presented here.

The paper is structured as follows: in section 2 the basic MOC equations are reviewed; in section 3 we give the operators of restriction and prolongation, and discuss the introduction of DFs; the low-order operator is introduced in section 4, where a complete discussion of the equations is given; in section 5 we present the application of our acceleration scheme to outer iterations; tests are given in section 6 for eigenvalue problems for a realistic BWR assembly and for the C5G7MOX Benchmark; closing conclusions are presented in the last section.

Method of Characteristics

The basis for reactor core, assembly-level transport calculations is the one-group linear Boltzmann equation. The most straightforward method to iteratively solve it is the source iteration (SI) which computes updated angular uxes assuming that the emission density is known:

           (Ω • ∇ + σ) ψ (n+1) (x) = q (n) (x)
x ∈ X,

ψ (n+1) (x) = ψ (n) _ (x) x ∈ ∂ _ X . (3)
In this equation, n denotes the iteration index, X = {x : r ∈ D, Ω ∈ 4π} is the phase space with its boundary ∂ _ X, and σ is the total cross section. The emission density reads as follows:

q = Hψ + S (4)
where S stands for the external source (including ssion contributions and transfers from the other groups), and

(Hψ) (x) = K k σ s k (r) l≤|k| A kl (Ω) φ kl (r) (5)
is the within-group-scattering term, written using the classical expansion of the collision term on spherical harmonics A kl (Ω).

In this section we briey review the basics of the method of the characteristics for unstructured meshes. More details about the angular and spatial approximations as well as the derivation of the equations can be found in [START_REF] Sanchez | A Synthetical Acceleration for Two-Dimensional Characteristics Method in Unstructured Meshes[END_REF][START_REF] Sanchez | Treatment of boundary conditions in trajectory-based deterministic transport methods[END_REF].

Approximations

The method is based upon two main spatial approximations. First a at-source approximation is made on each region:

q (r, Ω) = i χ i (r) q i (Ω) , (6)
where χ i (r) is the characteristic function of region i. Further, the angular ux is assumed to be constant across sectional area associated to each trajectory

ψ (r, Ω) = t Ω χ t (r ⊥ ) ψ t (z, Ω) , (7)
where χ t (r ⊥ ) is the characteristic function associated to the cross sectional area of trajectory t, z is the coordinate along the trajectory, and ψ t (z, Ω) is the angular ux along the trajectory.

The phase space in which the discretised transport problem is to be solved is dened as

X = {D i } × {Ω n } . (8)

Basic Equations

MOC balance and propagation equations are obtained by a projective technique involving Eq. ( 3)

and its integral form. The propagation equation is given by the integral transport equation across a region i and along a trajectory t parallel to the discrete direction Ω n :

ψ +,i (t, Ω n ) = ψ -,i (t, Ω n ) + β i (t, Ω n ) [q i (Ω n ) -σ i ψ -,i (t, Ω n )] , (9)
where the escape coecient is

β i (t, Ω n ) = 1 -e -σiRi(t,Ωn) σ i , (10)
and R i (t, Ω n ) is the length of the trajectory within the region. The average emission density q i (Ω n ) is written as:

q i (Ω n ) = C i (Ω n ) + S i (Ω n ) , (11)
where S i (Ω n ) is the average external source and C i (Ω n ) is the within-group transfer. The latter is given by:

C i (Ω n ) = Ki k=0 σ s k ,i l≤|k| A kl (Ω n ) φ kl i , (12)
where K i is the degree of anisotropy in region i. The moments φ kl i of the angular ux are obtained as follows:

φ kl i = n w n A kl (Ω n ) ψ i (Ω n ) . ( 13
)
The balance equation is obtained by integration of Eq. ( 3) over the volume of a region i:

t Ωn,t∩i w ⊥ (t, Ω n ) [ψ +,i (t, Ω n ) -ψ -,i (t, Ω n )] + σ i ψ i (Ω n ) V i = q i (Ω n ) V i , (14)
the sum in t being done for all trajectories with direction Ω n that intersect region i, and the weight w ⊥ (t) represents the orthogonal area associated with the trajectory.

Boundary Conditions

We consider two kinds of boundary conditions that are used for numerical solutions of the transport equation in the TDT solver:

• geometrical motions (e.g. translation, rotation or planar symmetry) are used to reduce the size of a domain that has these symmetries. They may be thought as "exact" boundary conditions. For this case, a special set of periodic trajectories, often called cyclic trajectories is used.

• albedo conditions introduce a physical approximation to represent the spatial and angular distribution of the particles re-entering the geometrical domain D.

For the albedo cases, the general approach in TDT is the following. The boundary ∂D is decomposed into a set of surfaces {α} and the entering and exiting angular domains {2π} ± into a set of angular subdomains {ρ} ± . The current of particles leaving (-) each surface α through angular domain ρ is evaluated. In the iterative scheme, the latter is used to compute the entering (+) current with the help of an albedo condition:

J ρ α,+ ⇒ albedo J ρ α,-. (15)
For a trajectory t entering through surface α with angular direction Ω ∈ ρ, the incoming ux is given by:

ψ _ (t, Ω) = 1 c ρ α J ρ α,-, (16)
where

c ρ α = α dS ρ dΩ |Ω • n| . (17)
In the following part of the paper, we will consider only isotropic albedo cases. In this case we have:

ψ _ (t, Ω) = 1 c α J α,-, (18)
with

c α = α dS 2π+ dΩ |Ω • n| = n ω n t Ωn,t∩α w ⊥ (t, Ω n ) . ( 19
)
3 General Framework

In this section we present the prolongation and restriction procedures and discuss the introduction of the discontinuity factors.

Prolongation and Restriction

In classical linear multi-grid methods [START_REF] Wesseling | An Introduction to Multigrid Methods[END_REF][START_REF] Hackbusch | Multi-Grid Methods and Applications[END_REF], the low-order problem is commonly solved for an additive error to the high-order calculation. Here, instead, the acceleration step consists in solving a simplied transport problem providing a scalar ux {φ h I } on the coarse level. This ux is directly used to reconstruct the ne angular moments resulting from the previous transport iteration (n + 1). The ne moments are consequently accelerated as follows:

φ kl i (n+1) acc = γ kl i (n+1) φ h I (∞) , ∀ i ∈ I . ( 20
)
To do this, we employ shape factors {γ kl i } depending on the unaccelerated ux. They are dened as:

γ kl i = φ kl i φI , ( 21
)
where φI is the volume-averaged scalar ux:

φI = 1 V I i∈I φ i V i . (22)
The prolongation procedure [START_REF] Santandrea | Acceleration technique for the method of characteristics in unstructured meshes[END_REF] is equivalent to applying a multiplicative correction to the high-order calculation [START_REF] Smith | Full-Core, 2-D, LWR Core Calculations with CASMO-4E, PHYSOR 2002[END_REF]. We note that, in a reference situation where the converged ne ux is available, the solution of the low-order problem must preserve the integral of the scalar ux on each mesh of the coarse spatial grid, i.e.

φ h I = φI . ( 23
)
This requirement is fundamental for the scheme to converge. Indeed, our choice for homogenisation and the introduction of the DFs are strictly tied to the fullment of this condition.

On the coarse level, isotropic sources are used. They are obtained from the high-order ones by a space-angle restriction procedure. The latter is applied whether to the emission density q i when computing the DFs, or to the external source S i when solving the low-order problem. Let be u i one of the above-mentioned sources, we have:

ūI = 1 4π n ω n i∈I u i (Ω n ) V i . (24)
For the spatial variable, this projection consists in volume-averaging on each mesh of the coarse grid.

Further, the integration over all directions eliminates the angular dependence on the coarse level. This permits to avoid building a table of correspondences between the two angular grids. We note that, by construction, the scalar source are preserved on both levels.

Homogenisation

This part has been inspired by Koebke's earlier works [START_REF] Koebke | A new approach to homogenization and group condensation[END_REF][START_REF] Koebke | Advances in homogenization and dehomogenization, Internal Topical Meeting on Advances in Mathematical Methods for the Solution of Nuclear Engineering Problems[END_REF] on homogenisation. The main idea is that an equivalent homogenised problem can be built for a given heterogeneous problem by introducing additional degrees of freedom such as discontinuity factors. Koebke's method dealt with a nodal approximation of the transport equation, and discontinuity factors were applied to the interface scalar uxes in order to full the so-called equivalence theorem. Following an approach analogous to Koebke's, we explain why the DFs have been introduced into our acceleration, and how this permits to respect condition (23).

The availability of an accurate heterogeneous solution of Eq. ( 3) is assumed in the following discussion.

The heterogeneous solution satises the integral neutron balance equation on whatever coarse mesh D I :

∂D I J (r) • dS + D I σ (r) φ (r) dr = D I Q (r) dr . (25)
We write an analogous equation for the homogenised low-order problem:

∂D I J h (r) • dS + D I σ h (r) φ h (r) dr = D I Q h (r) dr . ( 26
)
Let the scalar source Q h be dened as the volume-averaged ne one. The homogenised total cross section is taken piece-wise constant:

σ h (r) = I χ I (r) σ h I . ( 27
)
If the leakage term is preserved, we are able to ensure the preservation of the total reaction rate as well, for any arbitrary value of the homogenised total cross sections σ h . Thus, we are allowed to choose a ux-volume weighting procedure:

σ h I = D I σφdr D I φdr , (28)
which yields the preservation of the integral of the scalar ux and, therefore the fullment of condition (23). However, once homogenised total cross sections are chosen as in Eq. ( 28), additional degrees of freedom are necessary to preserve the leakage term. Some approaches can be found in litterature.

As an example, one can consider the coarse mesh nite dierence (CMFD) method [START_REF] Smith | Full-Core, 2-D, LWR Core Calculations with CASMO-4E, PHYSOR 2002[END_REF][START_REF] Joo | Dynamic Implementation of the Equivalence Theory in the Heterogeneous Whole Core Transport Calculation[END_REF] that acts as follows: the leakage is expressed in terms of the net currents on each surface of a rectangular homogenised region, diusion is used as a low-order operator and coupling correction coecients on each surface are introduced in order to match the average net currents. The novelty in this work is that a transport-like operator is employed for the homogenised problem. Furthermore, we decided to write the leakage in terms of the outgoing (+) and the incoming (-) currents as follows:

∂D I J (r) • dS = J +,I -J -,I , (29)
and to ensure the preservation of each one of them. We indicate with { Ĵ±,I } the currents given by the low-order transport operator and with {J h ±,I } the corrected ones. A discontinuity factor is dened for each one of them, in order to ensure the following:

J ±,I = J h ±,I = f ±,I • Ĵ±,I . (30)
The discontinuity factors {f ±,I } are calculated as the ratio between a reference value and the value produced by the low-order operator in a reference situation.

Summarising, the introduction of the DFs as additional degrees of freedom permits to dene an equivalent low-order problem which is able to reproduce the same result as for the reference heterogeneous problem when considering averaged properties. However, the problem in practical applications is that these factors, as the diusion-like coupling coecients for the CMFD formulation, should be computed a priori, since we do not have access to the converged ne ux before the iterative procedure is achieved. This is why a dynamic calculation of the DFs is performed: DFs are computed after each ne iteration, taking the latter as a reference situation.

V-cycle scheme

Our two-level method can be well described as a V-cycle. A general sketch is given in Fig. 2. After each ne transport iteration, an associated low-order problem is iteratively solved on the coarse grid, using the homogenised parameters derived from the previous transport sweep on the ne level. A maximum number M n of iterations may be xed a priori. The coarse solution is used to reconstruct the ne ux moments with the help of the shape factors. Then, the source is updated for the next transport sweep.

Acceleration of Inner Iterations

At each transport sweep, the region-averaged angular uxes are computed in terms of sources and boundary conditions:

ψ (n) _ (t, Ω) , q (n) i (Ω) → ψ (n+1) i (Ω) . (31)
Then, an associated low-order problem is solved in order to speed up the convergence of the inner iterations. In this section we point out the main aspects of our acceleration. The low-order problem and its iterative solution are discussed.

Low-Order Problem

On the homogenised coarse grid, we consider a transport problem analogous to the ne case. The MOC-formalism is used to iteratively solve it. Let be p the iteration index, we have:

           ∇ • Ω + σ h ψ(p+1) = q h,(p) in X ψ(p+1) = ψ(p) _ on ∂ _ X . (32)
In this equation q h is the coarse emission density and σ h is the homogenised total cross section obtained as:

σ h I = 1 V I i∈I σ i γ (n+1) i V i , (33)
where γ i is the shape factor for the scalar ux as in Eq. ( 21). The more coarsely discretised phase space is written as X = {D I } × {Ω m }.

As earlier mentioned, because of homogenisation and source projection, the currents computed by the low-order transport operator are not consistent with those resulting from the high-order calculation.

In particular, in a reference situation where the converged ux is available, the leakage term is not preserved on the two levels. Thus, we cannot directly employ ψ to compute the scalar ux needed for prolongation in Eq. [START_REF] Aen-Nea | Expert Group on 3-D Radiation Transport Benchmarks, Benchmark Specication for 2-D/3-D MOX Fuel Assembly Transport Calculations without Spatial Homogenization (C5G7MOX)[END_REF]. Nevertheless, this drawback can be overcome by modifying the iterative approach for the low-order problem. Unlike for the ne level, each iteration is performed to compute the global incoming and outgoing currents per macro, without evaluating region-averaged angular uxes.

The same sweep philosophy is used. Continuity of the interface angular uxes along each trajectory is still valid. Afterwards, the scalar ux {φ h I } is computed in a post-treatment phase involving the discontinuity factors previously computed.

The same spatial approximations as for the ne level are made. Thus, the source is assumed to be at on each homogenised mesh:

q h (r, Ω) = I χ I (r) q h I (Ω) , (34)
where χ I (r) is the characteristic function of region I. Further, we make the assumption that the angular ux is constant across sectional area associated to each trajectory as in Eq. ( 14):

ψ (r, Ω) = t Ω χ t (r ⊥ ) ψt (z, Ω) . ( 35
)

Equations

A MOC-like formalism is kept on the coarse grid. Indeed, like for the ne case, two main equations are employed: a propagation equation giving the outgoing (+) angular ux along each trajectory in terms of the incoming (-) angular ux and the internal source, and a balance equation for the scalar ux.

The propagation equation is obtained by the integral transport equation across a macro-region I and along a trajectory t parallel to the discrete direction Ω m :

ψ+,I (t, Ω m ) = ψ-,I (t, Ω m ) + β h I (t, Ω m ) q h I -σ h I ψ-,I (t, Ω m ) . ( 36
)
The escape coecient is dened in an analogous way to the ne case:

β h I (t, Ω m ) = 1 -e -σ h I R I (t,Ωm) σ h I , (37)
where σ h I is the homogenised total cross section dened in Eq. (33), and R I (t, Ω m ) is the length of the chord intersected in the macro-region I by the trajectory t. By analogy to the ne level, the coarse emission density q h I is written in terms of an external source SI , which is obtained by restriction (24), and a collision term C h I :

q h I = C h I + SI . ( 38
)
We will come back later to the collision term. Total outgoing and incoming currents are built up at each coarse iteration as follows:

Ĵ±,I = m ω m t Ωm ω ⊥ (t) ψ±,I (t, Ω m ) . (39)
These currents are corrected with the help of the discontinuity factors, and then used in the balance equation giving the coarse scalar ux:

f +,I • Ĵ+,I -f -,I • Ĵ-,I + σ h I φ h I V I = Q h I V I , (40)
where the scalar source Q h I is obtained by integration of the emission density q h I over all directions as follows:

Q h I = m ω m q h I . (41)
It remains to explain how the collision term is dened and how the emission density is updated.

Collision Term

As already mentioned, the low-order problem is dened in a reference situation where the ne converged ux is used and the source on the coarse level is obtained by projection:

q h I = qI . (42)
Now, in our iterative scheme, the coarse emission density is written in terms of the within-group scattering contribution and the external source as in Eq. ( 38). This means that, in the above-mentioned reference situation, further iterations must not modify the value of the coarse scalar ux in order not to compromise the convergence of the scheme, which entails:

C h I = CI = 1 V I i∈I C i V i . (43)
This can be achieved by choosing a ux-volume homogenisation procedure for the scattering cross sections:

σ h s I = 1 V I i∈I σ s0,i γ i V i , (44)
and simply dening the within-group collision term as:

C h I = 1 4π σ h s I φ h I V I . ( 45
)
As a consequence, the coarse emission density is iteratively updated as follows:

q h,(p) I = 1 4π σ h s I φ h,(p) I + SI , ( 46
)
where p is the iteration index for the coarse level calculation.

Boundary Conditions

Like for the ne case, two kinds of boundary conditions have been considered: geometrical motions using cyclic trajectories, and isotropic albedo conditions.

In the case of geometrical motions, the same approach as for the traditional MOC [START_REF] Sanchez | Treatment of boundary conditions in trajectory-based deterministic transport methods[END_REF] is kept. Therefore, the angular ux ψ (x) along a periodic compound trajectory requires the simultaneous computation of the two angular ux dened by a unit incoming angular ux in the absence of volumetric sources, ψbd (x), and by the volumetric sources with zero incoming angular ux, ψq (x). Once these two uxes have been computed, the nal angular ux is given by:

ψ (x) = ψq (x) 1 -ψbd (x) ψbd (x) + ψq (x) . (47)
In the albedo condition case, the boundary is decomposed into a set of surfaces {β}, and table of correspondences is established between the ne and the coarse surfaces. In particular we have:

S β = α∈β S α . ( 48
)
By analogy with the ne level, the incoming ux associated to a trajectory t entering through surface β with angular direction Ω m is given by:

ψ_ (t, Ω m ) = 1 c β J h β,- (49)
where

c β = β dS 2π+ dΩ |Ω • n| = m ω m t Ωm,t∩β w ⊥ (t, Ω m ) . ( 50
)
For the rst iteration we take:

J h,(0) β,-= 1 S β α∈β S α J (n+1) α,- , (51)
then, the incoming currents through each macro surface β are updated with an albedo condition derived from the ne level:

J h β,+ ⇒ albedo J h β,-. (52)

Computing DFs

Here we focus upon the leakage term. We discuss how discontinuity factors (DFs) are computed in the implementation of our scheme. Their evaluation is situated between a ne iteration and the solution of the associated low-order problem. After a transport sweep, the incoming (-) and the outgoing (+) currents on each coarse-mesh D I are available:

J (n+1) ±,I = n w n t Ωn,t∩I w ⊥ (t) ψ (n+1) ±,i∈∂I (t, Ω n ) , (53)
where ψ

(n+1)
±,i∈∂I is the angular ux entering (-) or leaving (+) a portion of a ne region's boundary ∂D i having a non-null intersection with the boundary ∂D I of the coarse region D I . It is worth noting that, in practice, a modication of the implementation of the ne transport iterations in TDT code is needed to compute these currents. More details about this point can be found in appendix A. We remind the reader that, when dening DFs, we want to preserve the leakage term on each coarse-mesh D I for a reference situation. In this work, we chose to achieve it by preserving the total incoming and outgoing currents.

The reference is given by the currents at iteration (n + 1) obtained from the emission density and the boundary conditions at iteration (n). Consequently, we consider the following xed-source problem:

           Ω • ∇ + σ h ψ = q on X ψ = ψ_ on ∂ _ X (54)
where the source is given by projecting {q

(n) i } as in Eq. ( 24) and the incoming ux is obtained as shown in subsection 4.4 using {J

α,-}. The total cross section is given by Eq. ( 33), and the more coarsely discretised phase space X is the same as in Eq. (32). Within the MOC framework, only one sweep is needed to solve Eq. (54). Starting from the boundary of the domain the angular ux leaving a macroregion D I is computed for a characteristic trajectory t via Eq. (36), in which qI is used. The outgoing ux is used as incoming ux for the following region along the above-mentioned trajectory. The total incoming (-) and outgoing currents (+) are given by Eq. (39). Then, the DFs are given as the ratio between the heterogeneous and the homogeneous currents as follows:

f ±,I = J ±,I Ĵ±,I ref .

(55)

These coecients are taken as constant in the low-order calculation and, after each coarse iteration, the incoming and outgoing currents are corrected as follows:

J h,(p)

±,I = f ±,I • Ĵ(p) ±,I , (56)
where p is the iteration index.

Note on Albedo Conditions -Boundary DFs

As above mentioned, on the coarse level the boundary is decomposed into a coarser set of surfaces {β}. The outgoing current through each surface α is reconstructed and updated as follows:

(J α,+ ) (n+1) acc = ζ (n+1) α J h β,+ (∞) , ∀ α ∈ β . (57)
The shape factors {ζ α } are dened using the unaccelerated leaving currents as:

ζ α = J α,+ Jβ,+ , (58)
where Jβ,+ is the average current on the macro surface β:

Jβ,+ = 1 S β α∈β S α J α,+ . (59)
We note that if the converged ne ux is available, the partial currents through surfaces β must be preserved on both levels, i.e. the following condition must be fullled:

J h β,+ = Jβ,+ . (60)
Now, the introduction of DFs ensures the preservation of the global incoming and outgoing currents on each region of the coarse discretisation. However, this does not entail that the current of particles leaving the boundary of the domain is preserved as well. Boundary discontinuity factors are therefore introduced for each macro surface β and the corrected currents are given as follows:

J h β,+ = f β Ĵβ,+ . (61)
This permits to respect condition (60). Like for the DFs, the computation of f β is performed between a transport sweep and the solution of the associated low-order problem. Reference values are represented by averaged currents at the ne iteration (n + 1) and sources at iteration (n):

ψ (n) _ (t, Ω) , q (n) i (Ω) → J(n+1) β,+ . ( 62
)
After each transport ne iteration, the currents for each macro surface β are available. The low-order operator provides Ĵβ,+ . Then, the boundary discontinuity factors are computed as follows:

f β = Jβ,+ Ĵβ,+ ref . ( 63
)
These coecients are taken as constant in the low-order calculation and, after each coarse iteration, the currents leaving the domain through surface β are corrected as follows:

J h,(p) β,+ = f β Ĵ(p) β,+ , ( 64
)
where p is the iteration index for the coarse-level calculation.

Acceleration of Outer Iterations

In the previous part of our work, we focused upon the one-group problem. A spatial-angle multi-grid acceleration for the inner iterations has been presented. In this section, an extension of the low-order problem equations to the k-eigenvalue multigroup problem is proposed.

Criticality Calculations

The multi-group formulation of Boltzmann equation for multiplying systems is typically solved by the method of the power iteration [START_REF] Lewis | Computational Methods of Neutron Transport[END_REF][START_REF] Nakamura | Computational Methods in Engineering and Science[END_REF][START_REF] Wachspress | Iterative Solutions of Elliptic Systems and Applications to the Neutron Diusion Equations of Reactor Physics[END_REF] as follows:

F (p+1) (r) = 1 k (p) ef f AF (p) (r) . ( 65
)
In this equation F (r) is the spatial distribution of ssion neutrons produced in the reactor:

F (r) = g (νσ f ) g (r) φ g (r) , ( 66
)
and A is the scalar multi-group transport operator. Power iteration is often referred to as outer iteration.

In TDT code, the eigenvalue is updated as follows:

k (p+1) ef f = k (p) ef f D F (p+1) (r) dr D F (p) (r) dr . ( 67
)
Each power iteration requires the iterative solution of a xed-source multi-group upscattering problem for the angular ux. Let be l the iteration index, we have:

L g ψ (l+1) g = H g→g ψ (l+1) g + g <g H g →g ψ (l+1) g + g >g H g →g ψ (l) g + 1 4π Q g f , (68)
In this equation, L g . = Ω • ∇ + σ g is the one-group transport operator, while H g →g is the slowing-down operator from group g to group g, and Q g f is the ssion source for group g. This problem is solved using a Gauÿ-Seidel iterative procedure. If no up-scattering is present, which is the case for fast groups, only one iteration is needed. Neglecting the iteration index, on each region Eq. (68) reads:

J g +,i (Ω n ) -J g -,i (Ω n ) + σ g i ψ g i (Ω n ) V i = C g i (Ω n ) V i + g =g C g →g i (Ω n ) V i + 1 4π Q g f,i V i , (69)
where C g →g i is the group-to-group collision contribution which is written as follows:

C g →g i (Ω n ) = k σ g →g s k ,i l<|k| A kl (Ω n ) φ kl,g i , ( 70
)
and the ssion source per group Q g f,i is written as:

Q g f,i = 1 k ef f x∈i χ g x F i,x , (71)
the sum being done over all the isotopes in i. The ssion source F i,x for an isotope x in a mesh i is dened as follows:

F i,x = g (νσ f ) g x φ g i .
(72)

Extension of the Low-Order Equations

The basic idea is to exploit the coarse scalar uxes computed while accelerating inner iterations, before performing a new power iteration on the ne level. To do this, after each ne outer iteration, we dene an associated low-order eigenvalue multi-group problem on the coarse level. The solution of this problem {φ h,g I } is then used, together with the previously computed shape factors per group {γ kl,g i }, to reconstruct the ux moments on the ne level:

φ kl,g i acc = γ kl,g i φ h,g I , ∀ i ∈ I . (73)
The so-obtained scalar ux is used to correct the ssion source:

(F x,i ) acc = g (νσ f ) g x φ g i acc , (74)
and consequently to accelerate the multiplication eigenvalue. The low-order problem must be dened so that, in a reference situation where the converged ne solution is available, the scalar balance equation per macro-region D I and per group g is respected on both levels. On the coarse level, the balance equation reads:

f g +,I • Ĵg +,I -f g -,I • Ĵg -,I + σ h,g I φ h,g I V I = Q h,g I V I . (75)
The DFs {f g ±,I } come from the previous ne outer iteration, as well as the homogenised total cross sections. The source Q h,g I is obtained by integrating the coarse isotropic emission density q h,g I over all

directions {Ω m }. The latter term is written as:

q h,g I = C h,g I + g <g C h,g →g I + g >g C h,g →g I + 1 4π Q h,g f,I , (76)
respectively in terms of the within-group, the down-scattering, the up-scattering terms and the ssion source per group g. By analogy with the within-group case, the group-to-group collision term is written as:

C g →g,h I = 1 4π σ g →g,h s I φ g ,h I V I , (77)
with

σ g →g,h s I = 1 V I i∈I σ g →g s0,i γ g i V i . (78)
On the other hand, the ssion term is expressed as follows:

Q h,g I = 1 k h ef f x∈I χ g x F h x,I , (79)
where the sum is done over all isotopes x in I. The ssion source is dened as:

F h x,I = g (νσ f ) h,g x φ h,g I , (80)
with

(νσ f ) h,g x = 1 V I i∈I (νσ f ) g x γ g i V i . (81)
Finally, k h ef f is the eigenvalue of the low-order multigroup problem. This value is computed as in Eq. 65.

A particular case of this method is the acceleration of thermal iterations. This consists in solving the Eq. ( 68) only for thermal groups with ssion source derived from the last ne outer iteration as follows:

Q h f,I = Qf,I = 1 V I i∈I Q f,i V i . (82)
In this case only the thermal uxes are reconstructed on the ne level.

Tests and Results

The aim of an acceleration method is to diminish, for a given class of problems, the overall computational cost, i.e. to reduce the total number of transport sweeps and the total computing time. Now, there is a trade-o between the accuracy and the computational cost of the low-order problem. An accurate approximation of the transport operator on the coarse level may drastically cut o the number of iterations but require an excessive cost in terms of computational eort and memory storage. On the other hand, a too poor approximation, which is easily solved, may not suciently reduce the global computational cost.

In order to test the eectiveness of our method, some tests have been performed. Results are presented in this section. All the calculations have been run on a DEC-Alpha 1000 MHz machine.

is signicantly reduced, the ratio between the number of tracks on the ne level and those on the coarse one is only ∼ 5.9. This is why the acceleration is expensive in terms of computing time.

In Table IV we give some results for the case using a coarser tracking for the acceleration. An initialisation-like procedure has been also implemented. The initialisation consists in performing a preliminary multi-grid-accelerated transport calculation using the low-order parameters (set of angular directions and spacing between parallel trajectories) for tracking on both levels. The same above-mentioned iterative strategy is kept for the initialisation problem although this calculation is not necessarily converged to the same precision since the maximum number of outer iterations is a priori xed to 5. The so-obtained approximated ux is then used as an initial guess. For results in Table IV, the ratio between the number of tracks on the ne level and those on the coarse one is ∼ 84.46. Consequently, the computing time for solving the acceleration equations is signicantly diminished, e.g. for the non-initialised case, with thermal acceleration, it passes from 89.52 to 4.59 sec. Moreover, the use of a signicantly reduced number of tracks on the coarse level permits to diminish the memory storage needs. Finally, the fastest calculation is obtained for the initialised case in which both inner and thermal iterations are accelerated.

C5G7MOX Benchmark

Here some results for the C5G7MOX Benchmark [START_REF] Aen-Nea | Expert Group on 3-D Radiation Transport Benchmarks, Benchmark Specication for 2-D/3-D MOX Fuel Assembly Transport Calculations without Spatial Homogenization (C5G7MOX)[END_REF] are presented. We compare our results to an APOLLO2 calculation using a DP 1 acceleration [START_REF] Santandrea | Acceleration technique for the method of characteristics in unstructured meshes[END_REF] for the inner iterations and a synthetic acceleration for the outer iterations [START_REF] Santandrea | Private Communication[END_REF]. The MOC calculation uses a 19188-region discretisation for the geometrical domain, and the trajectories have been tracked for all the directions in a product angular quadrature formula comprising 8 uniformly spaced azimuthal angles in (0, π) and 2 polar angles derived from a Bickley-Naylor formula in (0, π/2). On the coarse level we employ a 1414-mesh discretisation and a low-order angular representation. All tracking parameters for both levels can be found in Table V, while results are given in Table VI. The following iterative strategy has been adopted for both methods:

one transport sweep per group followed by an acceleration step for each outer iteration. Then, the acceleration for the outer iterations requires to solve a multigroup k-eigenvalue problem on the coarse level. The initialisation step consists in performing an unconverged multigrid-accelerated transport calculation by using the coarser angular representation on both levels, then we employ the solution of this problem as an initial guess. With reference to Table VI, we note that our method is slightly less eective than DP 1 for this benchmark: one more outer iteration is needed, which entails more transport sweeps (35 vs 28). On the other hand, the total computing times are comparable. It is worth noting that the two methods have a dierent allocation of the total computing time. In fact, the implementation of our method does not require extra-time in the building phase, dierently from the DP N technique.

However, the solution of the multigroup problem by our two-level method is more expensive (103.22 vs 37.50 sec). In this phase the most time consuming part is found to be the sweeping on the coarse level (64.34 over 103.22 sec). This is mostly due to the acceleration of the outer iterations which requires, as already mentioned, the solution of a low-order multigroup k-eigenvalue problem on the coarse level after performing an outer iteration on the ne level.

7 Conclusions and Future Work

In this work, we have presented a space-angle two-level acceleration for the MOC numerical solution of the neutron transport equation. Unlike the classical linear multi-grid methods for which a low-order problem is solved for the error at a given transport iteration, in our approach the acceleration step supplies a coarse homogenised scalar ux, which is used to reconstruct the moments of the ne angular ux. In a reference situation where the converged ne ux is available, the coarse solution must be able to reproduce the same averaged properties as the ne one. The implementation of our method has required a homogenisation/dehomogenisation technique, which has led us to the main approximations of our scheme: the discontinuity factors and the shape factors. Some parallels and distinctions with the CMFD approach proposed in literature have been drawn. Furthermore, we note that, since our method uses a MOC-like formalism to solve the low-order problem, no constraints stand against the fact of considering an unstructured grid on the coarse level. This is a clear advantage with respect to CMFD or other non-linear diusion-like acceleration techniques which generally require to superimpose a rectangular acceleration mesh on the heterogeneous geometry. This aspect makes our method attractive, for example, to whole core transport calculations. For the latter, an accurate geometrical description of elements close to the boundary (e.g. bae/reector representation) requires indeed an unstructured grid. Finally, a transport-like operator is employed on the coarse level, the correction provided by our acceleration step can be more accurate than that provided by CMFD for a class of reactor transport problems in which diusion modes are not dominant and transport eects due to higher modes become more important (e.g. accident situations).

Tests performed for the Atrium assembly have given encouraging results. They showed that our approach may be very eective in cutting o the number of transport sweeps. However, the solution of the acceleration problem may be too expensive if using the same set of trajectories on both ne and coarse levels. This is because the acceleration step is still based upon a MOC-like approach and, consequently, its cost depends on the number of tracks (i.e. the total number of intersections of trajectories with regions where the number of trajectories depends upon the angular quadrature formula and the spacing between the trajectories). Then a space-angle method, making use of a low-order tracking, permits to signicantly reduce the computing time and the memory storage needs for acceleration. Nevertheless, as a general consideration, we note that this further approximation of the transport operator on the coarse level may lightly worsen the total number of transport sweeps.

Furthermore, tests performed for the C5G7MOX benchmark showed that our method can perform well for reactor transport problems. It does not need extra work in the building phase, which is a clear advantage. However, the computational eort for solving the acceleration equations still remains high for such a case. The aim of this paper was to introduce a new method and show that it can be eective for a class of transport problems. Further work is needed for improving and optimisation.

Work is under way to investigate the extension of the method to more levels for reactor-core-size transport problems. We note that the use of more grids will not further diminish the number of transport sweeps, but only the computational eort needed for the acceleration phase. In this case, one can imagine to employ, on the rst coarse level, an accurate approximation of the transport operator and the phase space in order to drastically cut o the number of iterations, and use more coarsely dened problems on the lower levels in order to reduce the computation time.

A Modifying MOC on the Fine Grid Within the general framework of MOC, the tracking is built by computing and storing a set of geometrical data. In practical applications, a record is created for each trajectory:

{r i , l i | i = 1, M } , ( 83
)
where M is the number of intersections with regions; and r i , l i are the order number of region i and the track length of the trajectory across this region, respectively. Therefore, sweeping a trajectory consists in covering, in the order given by record (83), the M chords given by intersection of the trajectory with the ne regions. This approach has been kept to solve the low order problem, with the creation of an analogous record for each trajectory:

{r I , l I | I = 1, M } , ( 84
)
where M is the number of intersections with macro regions; and r I , l I are the order number of region I and the track length of the trajectory across this region, respectively.

In the implementation of MOC, one free iteration consists in computing the average angular ux ψ i on each region D i and, if needed, the currents leaving the domain through the discretised boundary α.

However, in our scheme, global incoming and outgoing currents on each macro region J ±,I have to be evaluated in order to compute DFs with Eqs. (53) and (55). In practice, the aim is to know when the boundary of a coarse mesh is crossed. This means that the ne discretisation has to know the coarse one, and that the tracking on the ne grid has to be reformulated to take into account both grids and the table of correspondence between them. To do this, we have introduced two more records for the trajectory on the ne grid:

• a coarse record giving, for each macro region I crossed by the trajectory, its order number:

{r I | I = 1, M } , (85)
• a record giving, for each of these macro-regions, the number of ne regions in it:

{n I | I = 1, M } . (86)
Consequently, the sweeping of a trajectory on the ne grid is modied. In fact, instead of a single cycle over the M ne chords, a double cycle is performed: a cycle over the M coarse chords, and for each one of them an inner cycle over all the ne chords in it. This allows us to know when the boundary of a coarse mesh is crossed, and to compute the currents we are interested in. Note that, if the same set of trajectories is employed on both grids, record (85) is derived directly from record (84).

(

Figure 2 :

 2 Figure 2: Visualisation of V-Cycle scheme for the acceleration of the inner iterations.

Table I :

 I Atrium Benchmark. Tracking parameters for the ne level. The number of tracks is the total number of intersections of trajectories with regions.

	Regions	Azimuthal angles	Polar angles	Spacing (cm)	Trajectories	Tracks
	3052	20 (uniform)	3 (Bickley)	0.02	1084	954,466

Table IV :

 IV Results for the Atrium assembly benchmark using the space-angle two-level acceleration method. All times are expressed in seconds. Two batches of results are presented: for the not-initialised and the initialised case.

			Free	Space-Angle MG
				Only Inner	Inner & Thermal
		building time	5.85	4.22	4.25
		solving time	555.24	62.39	45.93
	Not initialised	(ne+coarse+other) total time	561.09	(46.05+2.36+13.98) 66.61	(28.61+4.59+12.74) 50.18
		# transport sweeps	1403	103	64
		building time		4.11	4.15
		initialising time		3.33	4.80
		solving time		45.12	27.65
	Initialised	(ne+coarse+other)		(30.24+2.68+12.20)	(18.64+4.93+4.08)
		total time		52.56	36.60
		# transport sweeps		64	40

k ef f = 1.12854

Table VI :

 VI Results for C5G7MOX Benchmark. Our space-angle acceleration method is applied to both inner and outer iterarions. We compare our results with an APOLLO2 reference calculation. All times are expressed in seconds.

		APOLLO2 Reference	Space-Angle MG
	k ef f	1.18647	1.18637
	building time	85.89	17.70
	initialisation time	15.86	21.74
	solving time	37.50	103.22
	(ne+coarse+other)		(23.94+64.34+14.94)
	total time	139.25	142.66
	# transport sweeps	28	35
	# outer iterations	4	5

Acknwoledgement

The author would like to thank R. Sanchez and S. Santandrea for helpful discussions and suggestions during the preparation of this paper.

Atrium Benchmark

The case of a BWR assembly benchmark [17] is here presented. For our ne level MOC calculation, we use a 3052-region discretisation and a 6-group P1 cross sections library [START_REF] Aggery | Etude du Benchmark d'un assemblage REB de type MOX Atrium 10[END_REF]. Because of the enforced specular boundary conditions on all sides, cyclic trajectories have been used. The trajectories have been tracked for all the directions in a product angular quadrature formula comprising 20 uniformly spaced azimuthal angles in (0, π) and 3 polar angles derived from a Bickley-Naylor formula [START_REF] Leonard | Optimal Polar Angles and Weights for the Characteristics Method[END_REF] in (0, π/2). The tracking parameters for the ne level can be found in Table I. On the coarse level we have 136 meshes, with a ratio of ne to coarse meshes of ∼ 22.4. Two angular quadrature formulae have been considered on this level. All tracking parameters can be found in Table II. In the rst row we employ the same angular set as on the ne spatial grid, which is equivalent to consider a spatial multigrid acceleration. In the second row, a low-order angular approximation is used. All the calculations have been converged to the same solution (eigenvalue k ef f = 1.12854) with a relative precision of 10 -5 on inner iterations and eigenvalue, and a relative precision of 10 -4 on thermal iterations and ssion rates. For all the accelerated calculations, the following iterative strategy has been adopted: only one transport iteration per group is performed for each outer iteration, and the maximum number of inner and thermal iterations for the low-order problem is a priori xed to 15 and 4 respectively. The results are shown in Tables III andIV.

Those tables give global computation parameters such as the total computing time and the total number of transport sweeps. Total calculation time is divided into the dierent phases of the calculation: building (which includes reading geometry, tracking, and computing MOC parameters), initialisation (if present), and the solution of the k-eigenvalue multi-group problem. For the latter, the contributions from nelevel and coarse-level calculation are given between parenthesis. All the calculations are compared to the free-iteration results.

Table III refers to the case with the same set of trajectories on both levels: the number of transport sweeps is drastically reduced (ratio∼ 13.62), while the gain for the total computing time is not great (ratio ∼ 3.83). In fact, the computation cost required to solve the acceleration equations is signicant. Therefore, our gain in time does not show much in the overall calculation. This aspect is more evident when accelerating the thermal iterations. In this case, the further reduction of the total number of inner iterations is not followed by a consistent computing time gain. The total calculation time is greater than for the case without thermal acceleration (156.35 vs 146.30 sec). This is due to the nature of our low-order problem. In fact, since it employs a MOC-like iterative solution, its computational cost depends upon the total number of tracks. Therefore, although the number of regions on the coarse level List of Figures


L Di Pace

F Dacquait

P Schindler

V Blet

F Nguyen

Y Philibert

B Larat

Development of the PACTITER code and its application to safety analyses of ITER Primary Cooling Water System

Keywords: ITER, safety analysis, ACPs, PACTITER, PACTOLE

The PACTITER code derives from the PACTOLE code, developed by the CEA for predicting Activated Corrosion Products (ACPs) in PWR primary circuits. The operating conditions, material compositions and water chemistry of the various Primary Heat Transfer Systems (PHTS) of the International Thermonuclear Experimental Reactor (ITER) made mandatory the adaptation of the PACTOLE code.

PACTITER was developed on the basis of dedicated experiments, namely devoted to determine copper solubility and stainless steel release in the ITER primary cooling systems conditions, which are rather different from those in PWR (i.e. water chemistry and temperatures). The PACTITER code has been extensively used in support of the ITER Generic Site Safety Report (GSSR) in the field of accident analysis and worker collective dose assessment.

Introduction

For almost 30 years, the PACTOLE code has been developed by the Commissariat à l'Energie Atomique (CEA) to predict the formation, activation, migration and deposition of corrosion products in the PWR primary system [START_REF] Beslu | A computer code PACTOLE to predict activation and transport of corrosion products in a PWR[END_REF][START_REF] Nguyen | First numerical simulations of contamination of the PWR primary circuit by activated corrosion products with the PACTOLE V3.0 code[END_REF]. The other main purpose of the code application is estimating the effects of changes in primary system design and operating parameters. The PACTOLE code has been validated by comparison between simulations and more than 30 years of contamination measurements on PWRs [START_REF] Robin | Validation et qualification du code PACTOLE V2 (Validation and qualification of the PACTOLE V2 code)[END_REF][START_REF] Robin | PACTOLE: A Computer Code to Predict the Activation and Transport of Corrosion Product in PWRs[END_REF].

After an initial use of PACTOLE code, made by simply adapting the input to the ITER Primary Heat Transfer System (PHTS) characteristics, it was necessary to modify it and develop PACTITER to take into the presence of Cu-alloy as a new material in contact with water in the divertor cooling circuit. Afterwards, the modified PACTOLE code, called PACTITER, has been used in support of the Non-Site Specific Safety Report -2 (NSSR-2) (1996)(1997)(1998) and to the Generic Site Safety Report (1998)(1999)(2000)(2001)(2002)(2003) for accident analysis and worker collective dose assessment. Indeed, the Activated Corrosion Products (ACPs) in the ITER PHTS can be of some concern as contributor to the source term of potential released activity to the environment in case of accident (LOCA for instance) and to the Occupational Radiological Exposure (ORE) during the normal operation of ITER. Some dedicated experiments have been performed to acquire physical-chemical input data and validation data in PHTS conditions, particularly, in order to determine the corrosion product release term into the fluid due to the stainless steel corrosion in PHTS conditions, and the copper solubility in the ITER water chemistry conditions.

PACTOLE/PACTITER phenomenological model description

In a PWR primary system or an ITER PHTS, the formation process of ACPs is very complex. It involves many different mechanisms that react among each other. The first mechanism is the uniform and generalized corrosion of metallic alloys. For stainless steel materials, this leads to the generation of a dual oxide layer: an inner compact layer (chromite) and an outer porous layer (ferrite). The inner layer is a passive oxide layer, which limits ion exchanges between metallic alloys and primary coolant but does not eliminate them: ions are released in the primary coolant. The quantities of released materials are small (of about several mg/dm²/month) and do not alter component soundness.

The primary coolant transports ions generated by the corrosion-release phenomenon or by oxide dissolution. When the coolant becomes supersaturated in corrosion products, ions can precipitate on the walls or in the bulk of the fluid to form particles. Particles are also generated by erosion processes. Transported by the primary coolant, particles are deposited inside the circuits or they can agglomerate. Dissolution and precipitation depend on the corrosion product equilibrium concentrations, which depend on coolant chemical treatment (pH, H 2 concentration (or Redox potential), temperature).

Two types of radioactive corrosion product formation coexist. On the one hand, the activation of corrosion products occurs when they are deposited on surfaces under neutron flux. On the other hand, the corrosion of structural materials under neutron flux is accompanied by a release of radioactive corrosion products.

The PACTOLE and the PACTITER codes are based on a control volume approach, the primary circuit is represented by an arrangement of several volumes in which transient mass balance equations are solved :

( ) ∑ ∑ - = - + ∂ ∂ Source Sink m m e s i J J m m t m & & (1)
Where m i is the mass of the i th isotope in a considered medium defined as a form for the corrosion products, t is the time,

e s m m & & -
is the convective term (balance between input and output) and except for nuclear reactions, which are mass production/loss within the considered medium, J m is the exchange mass rate between two different media.

Seven different media are taken into account in a control volume and depicted in Figure 1. The transfer mechanisms, which are modelled in the mass balance equations, between all the considered media, are summarised in Figure 2.

The corrosion and release processes are formulated in the mass balance equations as functions of their time production rates. This formulation yields results on the consequences of material surface treatments, on the corrosion and release processes and, thus on the circuit contamination.

Application of PACTOLE to ITER: the development of the PACTITER code

The PACTITER code is currently being used for predicting the Activated Corrosion Product (ACP) inventory in the various ITER (International Thermonuclear Experimental Reactor) Primary Heat Transfer System (PHTS) circuits of the Tokamak Cooling Water System (TCWS) for accident analysis and personnel safety (ORE assessment).

At the very beginning (1994)(1995)(1996), the PACTOLE code was used for predicting activated corrosion products inventory and distribution in the first wall/shielding blanket cooling loops of the ITER PHTS [START_REF] Di Pace | Activated corrosion products in ITER first wall and shielding blanket heat transfer system[END_REF][START_REF] Nisan | Evaluation of the Activated Corrosion Product Source Term for ITER Heat Transfer System[END_REF], by simply adjusting the input deck to the case.

The major modifications in the PACTOLE input deck were directed at simulating the particular ITER PHTS operating conditions, which are quite different from those of PWRs. In PWRs, the coordinated B-Li chemistry is used to control a constant pH 300 °C to a target value of 7.2 characterized by the use of boric acid (H 3 BO 3 ) and Li hydroxide (LiOH) differently from ITER PHTS whose chemistry is closer to a Boiling Water Reactor water chemistry mainly controlled by its conductivity (less than 0.3 µS/cm).

The other main difference is the temperature range, which is 50-240 °C for ITER operation.

The other peculiarity of the ITER PHTS is the presence of copper alloys in the divertor cooling loops, that forced to consider a new element, copper, which was not included in PACTOLE. That required the realization of the new code PACTITER to include copper as the sixth element in addition to Fe, Ni, Cr, Co and Mn.

The three following nuclear reactions involving Cu, found to be the most important ones for safety analysis, were then included in PACTITER:

-63 Cu(n, γ) 64 Cu;

-65 Cu(n, 2n) 64 Cu;

-63 Cu(n, α) 60 Co.
First, a Cu solubility curve as a function of temperature was defined on the basis of data from a survey of the existing literature [START_REF] Eriksson | Aqueous Corrosion of Copper Alloys -Literature Survey[END_REF].

In parallel, an experimental campaign, performed by CEA Fontenay-aux-Roses, was launched in early 1996 to obtain Cu solubility data for ITER divertor PHTS conditions [START_REF] Di Pace | Development of the PACTITER code and its application to the assessment of the ITER Divertor cooling loop corrosion products[END_REF].

The ACP assessment for all ITER cooling loops were completed in early 1998 in support to NSSR-2, Volume III "Radiological and Energy Source Terms". The ACP assessment was then repeated to take into account design modification in the PHTS design (e.g. the heat exchangers with primary coolant on the shell side) [START_REF] Di Pace | Final report on Activated Corrosion Products Evaluation for the ITER-FDR PHTS[END_REF]. The comparison with experimental corrosion data related to ITER operating conditions and Cu-alloys provided a satisfactory result [START_REF] Di Pace | Development of the PACTITER code and its application to the assessment of the ITER Divertor cooling loop corrosion products[END_REF].

Further to the decision taken to update the ITER Final Design Report defined in 1998 (ITER FDR) to reduce size and cost of the machine still satisfying the overall programmatic objective of the ITER EDA (Engineering Design Activity), the machine design was reviewed including that of PHTS during the extension (1998-2001) of the EDA phase. The number of cooling loops of the PHTS was reduced from a total of 20 to 6. Again, the code PACTITER was applied to estimate the ACP inventory of the ITER PHTS according to the new design in support of accident analysis and worker collective dose assessment, with results documented in the Generic Site Safety Report (GSSR) [10].

ITER safety requirements for ACPs

ACPs will be present in ITER in the various in-vessel and vacuum vessel coolant loops as well as in any coolant loops related to test blanket modules, auxiliary heating or diagnostics equipment. These products impact occupational exposure, routine effluents to the environment, and potential releases during accidents. That has made the ACP inventory evaluation an important task for ITER public and occupational safety. ACPs present in fission reactors cooling loops do not pose a significant public hazard in these installations, however they are a significant hazard to site personnel. In fact, in case of severe accidents in LWRs, ACPs are negligible in terms of radioactive inventory, mobilisable inventory and environmental source term, compared to the analogous inventories from the reactor core (fission products and actinides). On the contrary for fusion devices such as ITER, ACPs could be important in terms of the whole mobilisable inventory in case of accidents compared to other major source terms, tritium and dust.

The enveloping ACPs inventory adopted for ITER accident analyses [11] was assumed to be 10 kg/loop as deposit and 0.6% of deposit mass in the suspended forms (soluble ions or crud particulates), irrespective of the cooling loop type, first wall//shielding blanket, divertor or vacuum vessel.

The other important safety requirements are the project release guidelines for the ACPs, tritium and dusts as shown in Table 1. For the occupational radiation exposure (ORE) doses, the limits and project guidelines have been adopted [12] as shown in Table 2.

ITER ACP assessment by PACTITER

Many calculations were carried out by the PACTITER code to assess ACP in the PHTS of ITER TCWS, following the several evolutions of the design.

The PACTITER geometric model of one first wall/shielding blanket cooling loop of the ITER PHTS (according to the Final Design Report 2001) is depicted in Figure 3.

The ACP assessment carried out for a single First Wall/Shield cooling loop provided a value of 1.8 g of suspended cruds and 10.9 g of dissolved ions. The total ACP deposit mass was 1.4 kg/loop. The hazard of the ACP is typically dominated by Mn-56 that arises from Mn and Fe in the base alloy steel, by Co-57 and Co-58, that arise from Ni in the base steel and by Co-60 that arises principally from Co in the base steel. The ACP deposit inventory is more than seven times lower than the limit of 10 kg/loop while the ACP inventory in the suspended forms is about 5 times lower than the limit. The ACP assessment for the divertor cooling loop provided values closer to the limits, 8.9 kg for the deposit mass and 42 g in the suspended forms. The ACP inventory for VV loops was estimated to be 6 kg/loop for mass deposit and 25.6 g/loop in the suspended forms.

The estimated base metal release rates for the different cooling loops of the ITER PHTS are shown in Table 3.

Mass release rates of the VV loop stainless steels are lower due to the lower coolant operating temperature (100 -103 °C) and to the lower coolant velocity (< 0.5 m/s) in most regions of the loop.

ACP mass and activity inventory were used in accident analyses involving LOCA with consequent releases to the environment. Among several accident sequences analysed [13], a large DV ex-vessel coolant pipe break provided the largest ACP environmental release, equal to 1.2 % of the project release guidelines (see Table 1).

The ACP radioactive inventory distribution, assessed by the PACTITER code for ITER PHTS, was used to determine the dose rates around its system components (e.g. piping, main pumps, valves, heat exchangers) during scheduled and unscheduled maintenance.

Starting from the radiation field and taking into account the anticipated work effort, the related annual collective ORE was estimated to be 65 pers•mSv [14] for the first wall/shield loops (30 pers•mSv) and divertor loop (35 pers•mSv). The contribution from VV loops was not estimated being the related ACP radiation hazard negligible compared to the other loops, as written before.

Experimental activity in support of PACTITER development and validation

Two main experimental activities were set up and carried out aiming at the development and validation of PACTITER.

The first one, launched in 1996 was performed at CEA Fontenay-aux-Roses, to get Cu solubility data in the ITER divertor PHTS conditions.

The second one was carried out at CEA Cadarache in two phases (2001 and 2004) to study the ACP release mechanism in close ITER coolant chemistry, temperature and velocity conditions. Both are described in more detail in the following sections.

Copper solubility data in ITER divertor PHTS conditions

The experimental campaign, performed by CEA Fontenay-aux-Roses, was launched in early 1996 to get Cu solubility data in the ITER divertor PHTS conditions [START_REF] You | Thermodynamic stability and solubility of copper and copper oxides[END_REF]. This experimental study was devoted to determine the valence states of the copper (Cu 0 , Cu + or Cu ++ ) and to obtain solubility in steady state conditions. In order to check the stable phases of Cu in ITER divertor coolant conditions, autoclave tests were performed on copper and copper oxides. These included four tests at 100 and 250 °C on CuO and Cu 2 O, pure copper and Cu-Cr-Zr (the reference Cu alloy for ITER divertor high heat flux components).

The major findings of this part of the study were: The solubility measurements, illustrated in Figure 4 [8], are quite different from those obtained from the literature [START_REF] Eriksson | Aqueous Corrosion of Copper Alloys -Literature Survey[END_REF]. The copper concentrations determined by tests are four or five orders of magnitude larger than the copper solubility obtained by thermodynamic data. Such difference could be explained by:

• the formation of copper colloids during tests;

• the presence of another Cu species, such as copper oxides previously formed during the exposure to vapours at room temperature or before establishing reducing conditions;

• the oxidation of Cu due to the existence of oxygen traces in water (copper is known as a very good getter of oxygen).

ITER ACP release mechanism

In order to get realistic values of release kinetics and thus to improve the understanding of ACP generation, experiments have been carried out in 2001 and 2004 using the CORELE loop (see Figure 5) in thermalhydraulics operating conditions envisaged for ITER PHTS and with the chemistry given in Table 4.

Having circulated in the irradiated tubes, the coolant (deionized water) is purified by passing through mixed beds of ion-exchange resins (Figure 5). No release can exist in the loop [cold part in polypropylene, hot part in zircaloy], except in the test sections.

Then, all the released ions are radioactive and trapped in the measurement ion-exchange resins. The nature and quantity of these radio nuclides are then determined by gamma spectrometry of the resins.

All the tubes considered in this study are made of 316 L Stainless Steel whose composition (given in Table 5) is close to an ITER Grade 316L SS.

The detected 58 Co, 59 Fe, 54 Mn and 51 Cr radio nuclides are then assumed to be representative of the main elements of the steel from whom they are generated (see Table 6). In the case of 59 Fe and 54 Mn, both can be considered as tracers of iron justifying the fact that the release rate of Fe is taken as the mean value of the 59 Fe release rate and the 54 Mn release rate.

For each chemical species, the release rate R (mg/dm²/month) is calculated using the following expression:

T S A w A R w ⋅ ⋅ ⋅ = (2)
where :

• A : Activity (Bq) of radionuclide in measurement resin at time t corresponding to the beginning of the experiment i.e. when the temperature in CORELE is reached and the coolant flows through the test sections. Actually, this activity is measured at the end of the experiment by disconnecting the measurement resins from the circuits of the loop and placing them in a dedicated device. This device enables the gamma spectrometry of the resins in a low environmental background. The measured activity is then corrected from the radioactive decay in order to correspond to time t.

• w : weight % of the element in the tube (mg/mg of tube) (manufacturer's data).

• A w : specific activity of tracing radionuclide in the tube at time t (Bq/mg of tube). This activity is deduced from the measurement at the date of irradiation corrected from the radioactive decay in order to correspond to time t.

• S : Exchange surface of tube with primary coolant (dm²). T : Duration of the test (month).

The global release rate is the sum of the different release rates relative to each radionuclide measured in the resins. They are summarized for each test in Table 7.

The two first tests were performed with the same irradiated tubes although the last one was performed with new freshly irradiated pair of tubes. The Li content, introduced as LiOH, was required due to the fact that the used water was not perfectly deionized and slightly acidic.

Release of chromium

Chromium was detected during the first test and was only trapped in measurement resin,

Influence of temperature

Globally, it can be seen that the lower the temperature, the lower the release rate in the tested range. These results could be explained by the existence of a diffusion barrier (possibly as a chromite) at the interface fluid/stainless steel. At low temperature (say less than 100°C), the solid diffusion of atomic oxygen in this layer could be considerably reduced, thus limiting the further formation of oxides (corrosion), and then the associated releases. As an example, the thickness of a Cr 2 O 3 layer is commonly of the order of a nanometer (L=10 -9 m). Considering a solid diffusion coefficient, D, of 10 - 21 m²/s (common value) leads to a characteristic diffusion time (L²/D) of 20 minutes compared to more than 300 hours of duration of the tests. Dividing D by 1000 when T is decreased from 150°C (423 K) to 100°C (373 K) could be possible if the activation energy would be higher than 180 kJ/mol.

Influence of velocity

From the comparison of the different tests conducted at about 1 m/s and 4 m/s (at the same temperature), one can conclude that the lower the velocity, the slightly higher the release rate. Note that the difference in velocity is obtained by a reduction in the hydraulic diameter thanks to an insert inside the test section.

The velocity is a key parameter in the release global process as it is responsible for erosion of the deposits and for mass transfer limitation in the hydraulic diffusion layer.

The erosion depends directly on the velocity whatever the temperature and the wall roughness in the tests conditions. No particles having been detected: one must conclude that, in these last 2004 tests conditions, erosion is not a relevant mechanism. On another hand, the mass transfer coefficient in the diffusion layer is given through the Nusselt number calculated, for instance, by the Dittus-Boelter correlation (see Table 8). It appears that an increase in coolant velocity due to the reduction in the hydraulic diameter globally decreases the mass transfer -and then the release rate.

Finally, two types of mass transfer limitations can be put forward in these ITER-like PHTS conditions :

the thermally activated solid diffusion of oxygen in a surface barrier (possibly as a chromite) which reduces the potential source of release by corrosion, the hydraulic diffusion which depends solely on Reynolds number instead of the velocity.

Conclusions

The application of the PACTITER code to the prediction of ACP for fusion reactors has been ongoing for more than 10 years. The extensive work of development, qualification and validation through dedicated experiments, and application of PACTITER to ACP assessment, mainly for ITER safety analyses (release in accident conditions and occupational dose to the staff), has produced valuable results. The ultimate aim is to further reduce the uncertainties in the source terms considering the specific material and operating scenarios of ITER machine. This effort will also fulfill some of the requirements for ITER licensing in Europe.

Acknowledgements

This work, supported by the European Communities under the Contract of Association between EURATOM on one side and CEA / ENEA on the other side, was carried out within the framework of the European Fusion Development Agreement. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Captions

Figures: Table 1 -ITER project release guidelines [10] Table 2 -ITER limits and project guidelines for ORE doses

•

  In reducing conditions, stable phase is Cu 0 , dissolution can occur according to the chemical reaction:Cu 0 (s) + H + (aq) ⇔ Cu + (aq) + 1/2 H 2(g) ;• Kinetics of reduction from Cu I and/or Cu II to Cu 0 is very slow.The other part of the experimental campaign was devoted to determine the total concentration of soluble copper in equilibrium with Cu-Cr-Zr alloy in chemical conditions close to those of the ITER divertor primary coolant. The solubility measurements were carried out in a flow device made of titanium. The measured pH during all the tests was in the range of 6.4 ±0.2.

  thus indicating the presence of a soluble form of chromium. This detection could thus be linked to the temperature of the test (150°C) and to the slight presence of oxygen which could favor the formation of a chromite layer (Cr 2 O 3 or MCr 2 O 4 , where M= Fe, Ni). The existence of such MCr 2 O 4 layer is coherent with experimental results indicating the release of Fe and Ni (upon the assumption that 58 Co is representative of Ni). It is worth noting that when the temperature is lower (100°C), no soluble chromium is detected when the O 2 content in the fluid is less than 24 ppb.

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 12345 Figure 1 : Media in a Control Volume for a PACTOLE representation

Figure 1

 1 Figure 1

< 1 g

 1 -T as HT and 0.1g-T as HTO and 1 g-metal as AP and 5 g-metal as ACP per year.Incidentsdeviations from normal operation, event sequences or plant conditions not planned but likely to occur due to failures one or more times during the life of the plant but not including Normal Operation < 1 g-T as HT or 0.1 g-T as HTO or 1 g-metal as AP or 1 g-metal as ACP or equivalent combination of these per event Accidents postulated event sequences or conditions not likely to occur during the life of the plant. < 50 g-T as HT or 5 g-T as HTO or 50 g-metal as AP or 50 g-metal as ACP or equivalent combination of these per event.Note: (*) HT: elemental tritium (including DT); HTO: tritium oxide (including DTO); AP: divertor or first-wall activation products; ACP: activated corrosion products.

Table 3 -

 3 ITER PHTS cooling loop material average release rate

Table 4 -

 4 Operating conditions for simulated ITER tests in the CORELE loop

Table 5 -

 5 CORELE test -Composition of the used SS316L tubes

Table 6 -

 6 CORELE test -Activation products of the main elements of the tubes

Table 7 -

 7 CORELE test -Measured release rates

Table 8 -

 8 CORELE test -Thermal-hydraulic parameters

Table 1 -

 1 ITER project release guidelines[10]

	Events or Conditions	Project Release Guideline(*)
	Normal Operation
	event sequences and plant conditions planned and
	required for ITER normal operation, including
	some faults, events or conditions which can occur
	as a result of the ITER experimental nature

Table 2 -

 2 ITER limits and project guidelines for ORE doses

	Dose Limits
	ICRP recommended limit for annual	20 mSv averaged over 5 years not to
	individual worker doses	exceed 50 mSv in any year
	Project Guidelines
	Annual individual worker doses	5 mSv/year
	Individual worker dose for any given	0.5 mSv/shift
	shift

Table 4 -

 4 Operating conditions for simulated ITER tests in the CORELE loop

	Operating conditions 150 °C / 150 bar	100 °C / 150 bar	100 °C / 120 bar
	Velocity V (m/s)	4.12	1.02	3.96	0.95	3.82	1.01
	Duration (hours)	353		332		335
	O 2 concentration
		4.15 <O 2 < 22.5	9<O 2 <24		3 <O 2 < 4
	(ppb)
	H 2 concentration
	(cm 3 /kg)	23.7<H 2 <26		23.5<H 2 <25.5	23.5<H 2 <25
	Li (ppm)	0.21		0.195		0.23
	pH 25°C measured	6.6		6.85		6.6

Table 5 -

 5 CORELE test -composition of the used SS316L tubes .2 0.093 16.77 1.640 2.06 0.29 0.007 0.012 0.070 0.440 0.031

	weight %	Fe	Ni	Co	Cr	Mn	Mo	Cu	S	C	N	Si	P
	From	67.4 11
	manufacturer

Table 7 -

 7 CORELE test -Measured release rates

	Temperature	Velocity	Species	Release rate	Global release rate
	°C	m /s		( mg / dm² / month )	( mg / dm² / month )
			Ni	4.7 ± 0.9
		4.15	Fe	16. ± 3	22 ± 3
	150		Cr	1 ± 0.4
			Ni	8.1 ± 0.7
		1.03	Fe	27 ± 5	35 ± 5
			Cr	0.27 ± 0.12
		3.98	Ni Fe	< DL < DL	< DL
	100		Ni	0.25 ± 0.1
		0.99	Fe	2.30 ± 1.3	2.6 ± 1.3
		3.96	Ni Fe	0.5 ± 0.3 0.27 ± 0.03	0.8 ± 0.3
			Ni	< DL
	100	1.01	Fe	2.7 ± 1.5	2.7 ± 1.5
	DL : Detection Limit