The paper presents the recent developments of the acceleration techniques for the Method Of Characteristics (MOC) in the code APOLLO2. The main contribution concerns the introduction of a multidomain DP N technique where all regions belonging to the same macro-domain are coupled by an integral operator that is strictly equivalent to the MOC. Different macro-domains are coupled via currents that are defined with the DP N formalism. This new Integral DP N (IDP N ) operator is build by using transmission and escape probabilities factors that respect symmetries/anti-symmetries and complementary properties that are enforced to preserve the physics of the problem and to save computational effort. These factors are computed using the numerical tracking of the MOC operator. The paper presents results on realistic assembly calculations that demonstrate the effectiveness of the IDP N operator as a synthetic acceleration tool.

INTRODUCTION

In this paper we introduce an evolution of the DP N operator currently used in the two dimensional MOC module of the APOLLO2 code [START_REF] Sanchez | APOLLO II: A user-oriented, portable, modular code for multigroup transport assembly calculations[END_REF] (called TDT) to accelerate the free iterations for the MOC. [START_REF] Halsall | CACTUS, a Characteristics Solution to the Neutron Transport Equation in Complicated Geometries[END_REF][START_REF] Smith | CASMO-4 characteristic methods for two-dimensional PWR and BWR core calculations[END_REF][START_REF] Dahmani | An efficient preconditioning technique using Krylov subspace methods for 3D characteristics solvers[END_REF][START_REF] Hong | CRX: a code for rectangular and hexagonal lattices based one the method of characteristics[END_REF][START_REF] Zim | Fundamentals and recent developments of reactor physics methods[END_REF] Much work has been dedicated in recent years to the theme of acceleration techniques for the MOC. The reason for all this interest is the capability of the MOC to solve unstructured geometry problems without using any homogenization procedure. As with every method that uses the source iteration algorithm, the MOC can be very slow for collision dominated problems. Therefore acceleration techniques are essential in order to permit a wider practical application of the MOC. The DP N method, originally introduced by Sanchez and Chetaine [START_REF] Sanchez | A Synthetic acceleration for a two dimensional characteristic method in unstructured meshes[END_REF], is the operator to be used in the framework of the synthetic acceleration procedure for TDT. Generalizations of this technique have been recently introduced in [START_REF] Santandrea | Acceleration Techniques for the characteristics method in unstructured meshes[END_REF][START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF] to improve the physical approximations of the initial formulation and to speed up the solution procedure. Further improvements are presented here and results are shown to illustrate the speed up factors of the TDT code (which is currently used for industrial calculations as a flux solver of the APOLLO2 [START_REF] Sanchez | APOLLO II: A user-oriented, portable, modular code for multigroup transport assembly calculations[END_REF] code system).

The principal contribution of this paper is a new way of using the DP N acceleration. Much alike the classical Interface Current (IC) [START_REF] Sanchez | New development in the Interface-Current Module TDT[END_REF] treatment of the Collision Probability Methods (CPM), we show that the MOC admits a similar Integral DP N (IDP N ) treatment. In this approach one subdivides the system into macro-domains in each of which one applies an integral formalism that does not introduce any supplementary approximation with respect to the MOC. The coupling between regions belonging to different macro-domains is done via interface DP N currents that result from some simplifying assumptions on the interface flux anisotropy. These assumptions are the only approximations of the method with respect to the MOC. One can expect this operator to behave better than the classic DP N operator, where the DP N approximations were made at the interfaces between every region and not only at the interfaces between macro-domains. The main differences with respect to the standard IC method are that our method is not based on the CPM but on the MOC, that the interface approximations are given by the DP N method and that the resulting equations have a response-matrix form that is different from that of the IC one. Since the numerical derivation of the IDP N method is consistent with the MOC, one can expect that this method will be a good acceleration tool. We have not investigated, as it has been done for the IC method, the behavior of the IDP N method as a stand-alone operator. Only its capability as a means to accelerate the MOC free iterations is considered here. Moreover, even if the IDP N method theory is valid also for three dimensional unstructured geometries, our practical implementations are limited to two dimensional systems, due to the actual capabilities of TDT.

The paper is organized as follows. In section 2 we recall the basis of the MOC and introduce some notation for the following sections. In section 3 the IDP N synthetic operator is introduced and its detailed properties are discussed together with the computational algorithm. In section 4 the external synthetic acceleration is revised and discussed. Section 5 is devoted to the analysis and discussion of results of calculations for a realistic assembly. Conclusions follow. Appendix A and B are dedicated to the description of improvements for the treatment of the equations for the boundary conditions, and to show many symmetries/anti-symmetries and reciprocity relations of the IDP N . We also demonstrate that this operator is robust in the limit of nearly voided media.

THE METHOD OF CHARACTERISTICS

The DP N acceleration equations are used as a low-order operator for the acceleration of the MOC. [START_REF] Sanchez | A Synthetic acceleration for a two dimensional characteristic method in unstructured meshes[END_REF][START_REF] Santandrea | Acceleration Techniques for the characteristics method in unstructured meshes[END_REF][START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF] The equations are obtained by the same projective technique that is used to derive the transport MOC equations but over different representation functions. We present here the transport MOC equations and in the next section the low-order DP N equations. Details can be found in [START_REF] Santandrea | Qualification Physique et amélioration du schéma numérique d'une méthode des caractéristiques pour l'équation de transport des neutrons[END_REF][START_REF] Sanchez | Treatment of boundary condition in trajectory based deterministic transport methods[END_REF].

The method of characteristics provides an iterative solution for the discrete ordinate formulation of the transport equation in a domain composed of unstructured homogeneous regions. The MOC equations are based on two approximations. The first consists of introducing a region-based flat spatial representation for the source term :

( ) ( ) ,
,

s i i i q A S φ ⎡ ⎤ ⋅ Σ + ∈ ⎣ ⎦ r ∼ Ω Ω , c e l l , i r (1) where ( ) ( ) { } A A ρ = Ω Ω is a vector of real spherical harmonics, { } , , s i s i diag ρ Σ = Σ
contains for each region the scattering coefficients associated to each spherical harmonic and i S denotes the external source. The components of vector i φ are the angular flux moments ( ) ( )

(4 ) 1 , 4 i d A π φ π = ∫ Ω Ω Ω i ψ (2)
where is the region averaged angular flux. All angular integrals are calculated with the discreteordinates quadrature formula.

( )

i ψ Ω
The second approximation is to use an expansion for the angular flux on region boundaries. For unstructured meshes one uses a set of parallel trajectories for each discrete ordinate direction Ω. The angular flux in direction Ω is then assumed to be constant over the cross sectional area associated with each trajectory,

( ) ( ) ( ) , , ,
t t t z z ψ θ ψ ⊥ Ω = + ∑ r r r r ∼ Ω Ω ⊥ Ω , (3)
and it is computed analytically along the trajectory. In this formula z and ⊥ r are, respectively, the coordinate in direction Ω and the coordinates over a plane orthogonal to Ω and ( )

t θ ⊥
r is the characteristic function associated to the cross sectional area of trajectory t.

The MOC equations comprise a balance equation and a transmission equation. Only the latter concerns us here. The transmission equation, which results from the analytical integration along each trajectory, gives the angular flux exiting the trajectory in terms of the incoming value and of the internal source. The homogeneity of the region as well as the flat source approximation are used to obtain:

( ) ( ) ( ) ( ) ( 1) , , ( ) 1 (
) ( ) ( ), i i i i i i R t R t n n n i i i i E R e t e t q ψ ψ Ω -Σ -Σ - + - Σ - = + Σ ( 4
)
where n is the iteration index, R i (t) is the length of trajectory t within region i and Ω its direction.

A transport iteration consists of a sweep over all the trajectories crossing the domain. This provides the fluxes leaving the domain and, via use of the region balance equations, the updated values of the region averaged angular fluxes. The latter are then used in equation ( 2) to update the region averaged flux moments that are compared to the previous ones to check convergence. If the iteration has not converged, the sources are updated via Eq. ( 1) and the transport iteration is repeated. Details of the full MOC equations can be found in [START_REF] Halsall | CACTUS, a Characteristics Solution to the Neutron Transport Equation in Complicated Geometries[END_REF][START_REF] Sanchez | Treatment of boundary condition in trajectory based deterministic transport methods[END_REF].

THE INTEGRAL DP N SYNTHETIC OPERATOR

A technique similar to the MOC is used to derive the DP N acceleration equations but with approximations (1) and ( 3) replaced by low-order ones. The source region expansion is as in (1) but with a smaller number of angular moments. The second approximation is more severe. The region boundary is decomposed into surfaces, a surface α being typically the common interface with a neighboring region, and a spatially constant angular approximation is assumed for the boundary fluxes:

( )

( , ) , , s A s S α α ψ ψ Ω Ω ± ± ∈ r r ∼ (5)
where S α stands for a generic surface and α ψ ± indicates exiting (+) or entering (-) surface averaged angular moments. The vector contains all the spherical harmonic of polar order k ( )

A Ω N ≤ . The index "N" indicates that a "Double P N " spherical harmonics expansion is used in (5) (one per side), this gives its name to the acceleration. Together with equation ( 5) the following currents are defined:

(2 ) ( ) ( , ) , s S J dS d n A A α ρ ρ α υ π ψ Ω Ω Ω Ω ± ± ± ± ± = ρ υ υ α α ψ = ∑ ∫ ∫ r i (6)
where coefficient A are defined below in [START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF]. We note that for the lowest moment this quantity is exactly the neutron current.

The first equation type of the DP N synthetic method expresses the current exiting a surface α in function of the interface fluxes of neighboring regions. This is done by using the same set of tracking trajectories used for the MOC. In the standard DP N , method approximation (5) is applied to the surfaces of all regions for all entering fluxes. This allows writing a transmission equation that couples all ,α ψ ± quantities. In this work we consider a partition of the system into macro-domains composed of homogeneous regions and we apply [START_REF] Hong | CRX: a code for rectangular and hexagonal lattices based one the method of characteristics[END_REF] to only the surfaces between macro-domains. We use the integral transport equation over a given trajectory "t". The angular flux exiting a region contained in macro-domain I is:

( ) ( , ) , 1 , j ij j i e x i j I j i j q e e e τ τ τ ψ ψ - - + - ∂ - ∈ < = + - Σ ∑ (7)
where τ(i,∂) is the optical length between the boundary of the macro-domain I and the entering boundary of region i along trajectory t, τ j the optical length over a region j, Σ j the total cross section in region j and τ ij the optical distance between the entering point in region i and the exiting point of region j. Here we consider a point as "exiting" or "entering" with respect to the given trajectory. Figure 1 gives a sketch of the situation.

Using the projection defined in (6) over an arbitrary surface, and applying the DP N approximation (1) to the macro-domain interface entering fluxes, gives: , , ,

,

i i I i I J T E q T ρ ρ υ υ ρ υ υ α α β β α υ β υ ψ ψ + + - - + - ∈∂ ∈ = + = ∑ ∑ Eq + (8)
where the source q υ i has been defined implicitly in [START_REF] Sanchez | APOLLO II: A user-oriented, portable, modular code for multigroup transport assembly calculations[END_REF] and E and T are escape and transmission coefficients respectively. The sum over β is done on region surfaces that belong to the surface ∂I of the macro-domain I, and the sign "+/-" indicates that one has to consider the entering/exiting current. Equation [START_REF] Santandrea | Acceleration Techniques for the characteristics method in unstructured meshes[END_REF] gives us an expression for all currents appearing in the system, regardless of their position with respect to macrointerfaces. The only approximation in ( 8) is to assume a reduced DP N flux anisotropy in the entering flux. In order to allow coupling between macro-domains, another DP N expansion is assumed for the currents exiting a macro-domain interface, that is to say:

, J A ρ ρ υ α α α υ ψ + + + = ∑ υ . ( 8a
)
The definitions and properties of the coefficients A,E and T are quite similar to those of the standard DP N method but, in order to stress some important differences, we recall them here. The matrix coefficients A α appearing in ( 6) and (8a), the transmission and escape coefficients in equation ( 8), are:

( ) ( ) (2 ) 1 , 4 A d S d nA A ρυ ρ υ α α π π Ω Ω Ω Ω + + = ⋅ ∫ ∫ ( ) ( ) ( ) ( ) ,( ) ( ) , ( ) 1 , 4 1 ( ), 4
where (β→α) (respectively (i→α)) indicates that the integration domain has to include only the trajectories going from surface β to surface α (respectively from region i to surface α), τ αβ is the optical distance between α and β, τ α ,(i) is the optical distance along direction Ω between the exiting side from region "i" and the surface α, and, finally, E(R i ) is defined in (4) has been used. It is important to note that for consistency all these coefficients are computed using the MOC trajectories for the integral over surface α. Moreover it is easy to see that the coefficients T and E obey symmetries/anti-symmetries and complementary relations. By "complementary" relation we mean a relation that links the two set of coefficients (E and T). The symmetry/anti-symmetry relation for the T coefficients are:

, T T s s T ρυ υρ ρυ υ ρ α β α β β α + - + - + - = = (10)
where s ρ stands for the parity of the harmonic ρ. It is important to stress how (10) can be used to greatly reduce the computational effort necessary to obtain the transmission coefficients. Actually, instead of computing the coefficients T by using [START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF] and sweeping over all angles, we strongly reduce the integration domain thanks to [START_REF] Sanchez | New development in the Interface-Current Module TDT[END_REF], and compute, sweeping the track, the following integrals:

( ) ( ) [ ] [ ] ( ) 0, 0, / 2 1 ˆ, 4 T d S d nA A e ρ αβ τ ρυ υ α β α β α ρ υ π π π + - - + - → ≤ ∈ × = ⋅ ∫ ∫ Ω Ω Ω Ω Ω (11)
where the conditions listed below the integral sign mean that we are considering only the lowest diagonal angular harmonics and the angular directions belonging to the first two octants. Thanks to [START_REF] Sanchez | New development in the Interface-Current Module TDT[END_REF] it is easy to see that:

ˆˆ, T T s s T ρυ ρυ ρυ υ ρ α β α β β α + - + - + - = + (11a)
for ρ<υ, while for ρ>υ the transmission terms are obtained directly from [START_REF] Sanchez | New development in the Interface-Current Module TDT[END_REF]. Calculations made with [START_REF] Santandrea | Qualification Physique et amélioration du schéma numérique d'une méthode des caractéristiques pour l'équation de transport des neutrons[END_REF] imply the use of less than half of the trajectories that would be needed were we to use [START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF]. A set of complementary relations are also valid for the escape coefficients:

( ) , , [ i I I i i E T T T A ρυ ρυ ρυ ρυ ρυ β β αβ α α β α β α β β δ δ δ Σ + + + + - + - ∉∂ ∈∂ ∈ = - - + ∑ ] . α + (12)
Equation ( 12) expresses the escape coefficients from a region "i" through a surface α as a sum over all the surfaces of the region. Three terms inside the sum must be distinguished. The first one, in which appears the symbol I β δ ∉∂ , comes from the surfaces that are completely internal to the macro-region I. The second one has to be included when the boundary surface β belongs to the border of the macro-region I. The last one has to be considered when β coincides with α. Expression [START_REF] Sanchez | Treatment of boundary condition in trajectory based deterministic transport methods[END_REF] gives the escape coefficients only in terms of the transmission ones. Therefore, the transmission coefficients are the only quantities directly computed from tracking.

Equation ( 8) must be completed to include the boundary conditions, that is to say, by expressing the entering boundary currents in terms of the exiting ones. A detailed description of these equations can be found in [START_REF] Santandrea | Qualification Physique et amélioration du schéma numérique d'une méthode des caractéristiques pour l'équation de transport des neutrons[END_REF]. In this paper we introduce an important improvement for the boundary currents equations that is detailed in appendix A. This treatment is common to both the standard DP N and IDP N methods.

Since the source q contains flux moments, supplementary equations are needed to close the problem.

These are obtained from the DP N region balance relations [START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF]:

( ) , . i i i i i i J J V V B υ υ α α q ρ ρ υ υ α υ υ φ + - - + = Σ Σ ∑ ∑ (13)
Here B is an operator that takes into account the lack of numerical orthogonality for the spherical harmonic functions and is defined as:

( ) ( ) ( ) 4 1 Ω . 4 B d A A ρυ ρ υ π π = Ω ∫ Ω (13a)
Equation ( 13) becomes ill-conditioned for voided or weakly absorbent regions and is therefore cast into a more stable form. The algebraic manipulations are parallel to the ones used in the standard DP N method [START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF], but many substantial differences appear due to the full coupling of the fluxes for regions belonging to the same macro-region. Another benefit of the procedure is to reduce the size of the problem to be solved iteratively. In fact, the iteration problem then deals only with the surface currents.

Substituting the expression (8) of the currents into formula ( 13), one obtains the following system per macro-region I :

,

I I I I I e x t M I Gq φ ψ - = + (14)
where the expressions for matrices I, G, M are given in appendix B. It is important to note that the coefficients of the "incoming" matrix I can be directly obtained from the escape coefficients E given by ( 12), as detailed in equation (B2). This relation allows further computational savings. The matrix G of ( 14) is defined in (B3a) as

, BV G= -E Σ
where the operator E can be defined as a function of the transmission factors as showed in (B5), and B is the matrix (13a). An analytical definition of the operator G is given in appendix B (formulas (B7) and (B7a)), and there we show also that this operator is well defined also for nearly voided media. Moreover we also prove that:

ij ij ji G G s G s ρυ υ ρ ρ ρ υ = = υ . ( 15
)
Equation ( 15) enables us to reduce the computational effort by adopting a strategy similar to the one sketched in [START_REF] Santandrea | Qualification Physique et amélioration du schéma numérique d'une méthode des caractéristiques pour l'équation de transport des neutrons[END_REF] and (11a) for the transmission factors.

The numerical counterpart of (B7) and (B7a) in the framework of the MOC is straightforward, and it is the way we have chosen to compute the operator E. The other possibility to compute E I , the escape matrix for a given macro-region I, is described in appendix B, and makes use of reciprocity relationships between E I and the transmission factors between all surfaces of I. We have opted for the first way since the storage required by E I is by far lower than the one required by the transmission factors. System ( 14) is then used to express the flux moments as:

ˆ, ext I Gq

φ ψ - = + (14a)
where

1 Î M I - =
and . These expressions are inserted back into the transmission equation [START_REF] Santandrea | Acceleration Techniques for the characteristics method in unstructured meshes[END_REF] through the source term (1). This gives:

1 Ĝ M G - = ( ) ( ) 1 1 ŝ s d ext A T E I A E G I q ψ ψ + - - - = + Σ + Σ + , ( 16
)
where

{ } , , s i s i diag ρ Σ =
Σ and I d is the identity matrix. It is worth noting that the matrix T+EΣ s Î, that appears in ( 16), is sign-symmetric (as shown in appendix B), a property that is exploited to reduce the computational effort. Equation ( 16) is the final iterative form used to solve the problem and, since it is written for currents between macro-regions, it has a reduced dimension with respect to the standard DP N method. This means that problem (16) should be faster to solve than its equivalent in the standard DP N method. Once the currents have been obtained by solving (16), volume fluxes are obtained using (14a).

All transmission and escape coefficients can be shown to be bounded. For the transmission coefficients, this fact is clear from their definition [START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF]. On the contrary, to show boundedness for the escape coefficients, it is necessary to develop and implement numerically an asymptotic expansion of equation [START_REF] Sanchez | Treatment of boundary condition in trajectory based deterministic transport methods[END_REF].

In fact, it is clear that expression [START_REF] Sanchez | Treatment of boundary condition in trajectory based deterministic transport methods[END_REF] is greatly ill-conditioned for nearly voided systems. This development is given in appendix B.

Problem (16) has the same response-matrix form as the standard DP N operator. In [START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF] it has been shown that a numerous set of preconditioning and optimization techniques can be used for this problem. This is an important point for our acceleration technique since all previous work, done in the framework of the standard DP N operator, can be re-used for the IDP N acceleration.

All the geometry partitions that we have used in this article have been obtained by the Spectral Bisection Algorithm (SBA) [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF], that is a standard way to compute domain decompositions.

EXTERNAL SYNTHETIC ACCELERATIONS

In this work we have also applied the new IDP N acceleration as a multigroup acceleration technique to improve the rate of convergence in multigroup calculations. Until now, most of the work that has been done to accelerate the external iterations used some non-linear strategy, [START_REF] Zim | Fundamentals and recent developments of reactor physics methods[END_REF] some classical linear extrapolations (like the Chebyshev one) or some Krylov subspace techniques. [START_REF] Warsa | Krylov Subspace Iterations for Deterministic k-Eigenvalue Calculations[END_REF] Recently, Suslov [START_REF] Suslov | An algebraic collapsing acceleration method for the acceleration of the inner (scattering) iterations in long characteristics transport theory[END_REF] proposed a linear algorithm using a mono-group synthetic operator in the external framework. We implemented this idea in APOLLO2 with slight modifications, and used the standard DP N and the new IDP N operator to test its efficiency. In this section we describe the algorithm and discuss implementation details.

The multigroup external equations are driven by fission and, written with the iteration index, read :

( ) ( ) ( ) ( , ,1 ( 1)
1

p j p j p d u p L H H F ψ ψ ψ λ - - - = + ) 1 -, ( 17
)
where L is the one-group transport operator, H d is the group-to-group and slowing down scattering operator, H u is the up-scattering operator and F is the fission operator. The upper index p denotes the external iteration index, while j is the thermal iteration index. Denoting by J tot the number of thermal iterations per external sweep, it is assumed in (17) that ( ) ( )

, tot p pJ ψ ψ =
for every external iteration p. One can then write the converged solution as

( ) ( ) ( ) , p j p j f ψ ψ ∞ = + , , (18)
where

( ) ( ) ( ) F L H ψ ψ λ ∞ - = ∞ . ( 19
)
We substitute formula ( 18) in ( 19), taking into account [START_REF] Zmijarevic | [END_REF], to obtain:

( ) ( ) ( ) ( ) ( ) ( ) ( ) , 1 , , , 1
, given by (20). This last problem is solved using once again the power iteration solution strategy, that is to say by fixing the terms coming from fission and compute a multigroup solution with the help of "j" thermal iterations. The new estimated eigenvalue for the next iteration is given by the usual condition:

p j p p j p j p j s u u p F L H f H f f H F ψ ψ ψ λ λ λ - - ⎛ ⎞ - = + + ∆ + - ⎜ ⎜ ⎝ ⎠ ⎟ ⎟ (20)
( ) ( ) ( ) ( ) ( ) ( ) ( ) , , , , , ,
p j p j k k p j p j k F f F f ψ λ ψ - + == + 1 ,
where the index k refers to the external iterations necessary to solve the problem (20).

RESULTS AND COMPARISONS

In this section we will prove the effectiveness of the IDP N acceleration, presenting some results related to the Atrium Benchmark.

[16] Here we will compute a version of this benchmark that makes use of a sixteen group P 0 transport corrected cross section library that has been obtained by a previous self-shielding phase in APOLLO2. [START_REF] Zmijarevic | [END_REF] The Atrium benchmark is a severe realistic assembly calculation, that has strong transport effects and slow convergence rate. The spatial computational mesh used for our calculation is shown in figure 2, and it consists of 3560 regions. A full description of the case can be found in [16], and we here only insist on some basic features. The Atrium benchmark represents a BWR assembly. The fuel and gadolinium rods are surrounded by an external as well as an internal box, where homogeneous water flows.

As it is typical for these reactors, liquid water concentration varies inside the assembly. In fact, for the considered version of the benchmark, we have a 0% vapor fraction inside both the internal and the external box of the assembly, while we have a 40% vapor fraction for the water in direct contact with fuel rods. Fuel rods are characterized by three different enrichments, and strong flux heterogeneities appear next to gadolinium rods, that are distinguishable in figure 2 by the very fine mesh refinement necessary to represent flux gradients. We used the transport code APOLLO2 where all the developments we present here are now implemented. The tracking data used to compute the benchmark are shown in table 1. We used a product angular formula with sixteen uniform azimutal angles, and two optimal Bickley polar angles. [START_REF] Leonard | Optimal Polar Angles and Weights for the Characteristics Method[END_REF] The transversal spatial step is 0.025 cm. This results in a total of 643,617 horizontal intersections with computational regions.

A first comparison deals with two calculations that have used a DP 1 or IDP 1 multigroup initialization and a thermal synthetic acceleration, a procedure that is also discussed in [START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF]. Table 2 gives the performances of the old DP 1 operator and of the new IDP 1 one. Free iterations are also shown in order to illustrate the efficiency of the acceleration methods. All the results in the table, with the exception of the free iterations ones, were obtained using a mesh grid that was constructed using the recursive SBA. [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] For the standard DP 1 method, the chosen grid is used to make a two-level multigrid acceleration method [START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF], while for the IDP 1 operator we use this same grid to compute macro-domain integral probabilities. We have always used the adapted ILU0 preconditioning technique [START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF] to speed up the DP N and IDP N equations. Due to their great impact over the total computational time the building times are also shown. All computing times are expressed in seconds on a 2.8Ghz 512Kb cache PC.

Building times strongly depend on the preliminary manipulations required: algebraic elimination process for the standard DP 1 method, and computation of the transmission probabilities for the IDP 1 method.

As the number of regions per macro-region increases, the building time also increases, but the burden of the accelerated iterative resolution strongly decreases. Therefore the appropriate number of regions to be included in each macro-region should be optimized. This phenomenon was already described in [START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF] and we

will not insist on it here. Our numerical tests show that a good number of regions per macro one for the Atrium benchmark is 70, so we always retained the grid mesh resulting from the recursive SBA with this fixed maximum number of regions. It is important to stress that, while for the standard DP N method the grid mesh has no physical impact, for the IDP N it should be beneficial to take advantage of a "physical" grid, that is a grid with reduced flux anisotropy on grid interfaces. This is the case for example for the classical IC method where the grid is given by the region decomposition of the assembly. In fact one can suppose that flux anisotropy is reduced at region interfaces. In this work we will always use the "numerical" grid given by the SBA, even if it certainly does not guarantee the best result.

As one can see from table 2, the best time with consistent tracking, 26.4 sec., is obtained for the IDP 1 operator. The standard DP 1 operator spends a lot of time in the initial building phase, and requires 151 internal transport iterations. In spite of the non-"physical" SBA grid mesh, the spectral behavior of the IDP 1 operator is quite better since it converges with only 102 iterations. The thermal synthetic acceleration, used with the IDP 1 method, reduces the number of total transport iterations by almost a factor of 40 with respect to free iterations. Unfortunately, the IDP 1 acceleration (as well as the DP 1 one) requires an expensive building phase that accounts for half of the total time; this explains why the accelerated calculation is only 15 times faster than the free iterations.

In order to reduce the expensive building phase, we introduced a low-order tracking described at the third row of table 1. The low-order tracking is now used to compute a synthetic IDP 1 operator: this choice entails a lower building cost but can cause a poorer convergence rate. If one compares the data for the reference and low-order tracking, one can see that the number of tracks that are used to build the IDP 1

operator is now decreased by a factor of four with respect to the reference case. We will call this approach "inconsistent" since it is generally required that the acceleration operator be obtained in a "consistent" way with respect to transport, in order to avoid problems of poor convergence or even of divergence of the accelerated algorithm. Using an inconsistent acceleration can therefore lead to lack of robustness. At this time we do not suggest a method to prevent these drawbacks. In the results we show, and more generally on the basis of our experience, we can state that the behavior of the acceleration is not degraded too much if one maintains the low-order tracking formula at an acceptable level.

The new IDP 1 results are given in fourth row of table 2. The use of a low-order tracking should entail poorer IDP 1 performances. The results in table 2 show that this effect is very limited, since the total number of iterations goes from 102 to 107. On the contrary, thanks to the faster building time, the total computational time decreases from 26.4 to 23.6 seconds. Note that while most of the IDP 1 building time is directly proportional to tracking, the standard DP 1 building time strongly depends on the geometrical structure of the problem, and weakly on tracking. The use of a low-order tracking has practically no effect on the standard DP 1 operator, and therefore we do not present inconsistent DP 1 results in table 2.

A last batch of results (table 3) deals with the impact of the external synthetic acceleration. We used this algorithm in conjunction with the standard DP N operator or the new IDP N operator. Both have been computed inconsistently with respect to transport: we used the low-order tracking described by the third row of table 1. In table 3, one can see that the external acceleration is very effective for the synthetic IDP 1 operator, since it reduces the total number of internal transport iterations from 106 to 48. Unfortunately this yields only a slight time reduction, the total computational time goes from 23 to 22.3 seconds. This is due to the fact that the synthetic external acceleration "transfers" a great amount of work from the transport to the synthetic operator. And in fact the total number of synthetic operator problems to be solved is multiplied by more than a factor two, with respect to the cases presented in table 2.

The only possibility to further reduce the computational time is to solve the IDP N problem faster, without increasing building times. We suppose that this could be possible by adopting a full multilevel algorithm (and not a purely two-level one as in this paper). We will address this issue in future works.

CONCLUSIONS

In this paper a new DP N synthetic operator has been introduced. Its main approximation consists in supposing reduced flux anisotropy over a limited number of surfaces. In fact, after having subdivided the system into some macro-regions, the IDP N operator supposes limited flux anisotropy only over those surfaces that divide macro-regions. Thanks to this preliminary hypothesis, and by making use of the basic MOC equations, one can find an integral operator whose unknowns are the flux components over macrointerfaces. Many symmetries/anti-symmetries and complementary relations help to minimize the computational effort, and taking into account important physical properties in the numerical operator.

Calculations made with APOLLO2 on the Atrium Benchmark show that the IDP 1 method is more efficient than its standard counterpart, and that it can be retained as an efficient acceleration tool. The preliminary time necessary to build the operator is reduced and the spectral performance is improved. The behavior of the new IDP N operator in conjunction with the external synthetic acceleration approach has also been evaluated. Even if the overall computational time is reduced, further improvements of the IDP N operator are necessary to make this reduction significant. In particular, we suggest adopting a true multilevel algorithm that will help to reduce building and solution times. We will address this issue in future works.

APPENDIX A: BOUNDARY CONDITIONS TREATMENT

In this section we give a detailed description of the boundary current treatment. The complete form of the boundary current equations has been given elsewhere. [START_REF] Santandrea | Acceleration Techniques for the characteristics method in unstructured meshes[END_REF] What we describe here is an improvement in the treatment of these equations for the DP N problem. We shall show that for all local boundary conditions (i.e. all boundary conditions where the neutrons exiting the system re-enter through the same surface), exiting currents can be eliminated from the iterative part of the problem with almost no computational and memory cost. Let us write the transmission problem (8) in the following manner:

i i T T E S ψ ψ ψ ψ βψ + - - ∂ ∂ - + ∂ ∂ ⎧ = + + ⎪ ⎨ = ⎪ ⎩ , ( A1
)
where we have distinguished the transmission contributions coming from internal currents "ψ -" and from external currents "ψ - ∂ " . Each entering external current can be written as

( ) t α α ψ βψ -= + S . (A2)
For local boundary conditions t(α)=α, but in general surfaces do not coincide. For both cases equation (A1) can be re-written as

i i T T E ψ ψ ψ + - + ∂ ∂ = + + (A3)
where T T β ∂ ∂ =

.

The system of equations (A3) is a reduced form of the initial system (A1): it contains less unknowns and equations but shares the same profile for the remaining equations. This means that contrarily to standard elimination processes, where fill-in always appears, here no elements are added to the initial operator. In the case of no local boundary conditions, the further elimination of the boundary exiting currents appearing in (A3) would cause fill-in and, for this reason, we decided not to eliminate them. For local boundary conditions the situation is different. In fact one can extract the boundary external currents from (A3), and write them as

( ) ( ) ( , ) ( ) ( ) , i i ij j i j i T T E , q ψ ψ ψ ∂ + - ∂ ∂ + ≠ = + + ∑ , (A4) where ( ) ,
i ψ ∂ + is the outgoing external modes related to the i-th external surface, T (i) is the transmission matrix giving the contribution of internal currents to the i-th external boundary current, is the transmission matrix describing the link between the j-th external surface and the i-th one and describes the escape.

It is clear that the only case when is different from zero is when the i-th and the j-th surface share the same region. For convenience we will suppress from now on the tilde notation over the T ( , )

i j T ∂ ( ) i E q ( , ) i j T ∂ ∂ symbol since
each external transmission term will always be post-multiplied by the albedo operator. The objective is to show that building the formal inversion

,

i j i j I T - ∂ ⎛ ⎞ - ⎜ ⎝ ⎠
∑ ⎟ is not computationally expensive and can be done without any fill-in in the LU decomposition necessary to compute this inverse. To see this, let us eliminate the first external current from the system (A4), to

i i i i i j ij j i i j i j D I T T T T T T T T E T E q ψ ψ ψ ∂ ∂ ∂+ ∂ - ∂ ∂ ∂ ∂+ ∂ ≠ > ⎡ ⎤ ⎡ ⎤ ⎡ - = + + + + + ⎣ ⎦ ⎣ ⎦ ⎣ ⎤ ⎦ ∑ , obtain ( ) 1 ( ,1) (1, ) ( ) ( ) ( ,1) (1) ( ,1) (1, ) ( , ) ( ) ( ) ( ,1) (1) , , 1 i i
which is valid for all i>1. Rearranging these equations to a unit diagonal gives:

( ) ( ) ( , ) ( ) ( ) , 1 , 1 , 1 1
,

i i i j j i j i j T T E q 1 i ψ ψ ψ ∂ + - ∂ ∂ + ≠ > = + + ∀ ∑ > , (A6) where , ( ) ( ) ( ) ( ) ( ) ( ) 1 ,1 1 1 1
i i i i T D T T T - ∂ = + ( ) ( ) ( ) ( ) ( ) ( ) 1 , , 1 1 , 1 ,1 i
j i j i j i T D T T T - ∂ ∂ ∂ = + , ∂ and
( ) ( ) ( ) ( ) ( ) ( ) 1 ,1 1 1 1 i i i i E D E T E - ∂ = + .
System (A6) is formally similar to system (A4), with the exception that the number of unknowns is reduced by one in all equations but the first, and that the profile of the operator has been correspondingly reduced since each coupling term due to the first external current has been eliminated from all equations. A similar reasoning can be applied to all the other currents, so that at the final step one obtains an upper diagonal operator. Before the l-th step one has to eliminate the l-th current from the equations greater than l.

After the elimination of this current from the i-th equation of the system (where clearly i>l) this equation reads ( )

( ) ( , ) ( ) ( ) , , ,
i i i j j l l j i j l T T E i l q ψ ψ ψ ∂ + - ∂ ∂ + ≠ > = + + ∑ (A7)
where, after defining

( ) ( ) , , i l l i l i D I T T ∂ ∂ = - , one has ( ) ( ) ( ) ( ) ( ) ( 1 , 1 , 1 1
i i i l l i l l l T D T T T - - ∂ - - = + ) l l , 1 l ∂ - , (A8a) ( ) ( ) 1 , ( , )
,

i j i j l l i T D T - ∂ = . (A8b)
Equations (A8a-b) allow the calculation of all matrices

( ) i l
T from all matrices at the preceding step ( ( )

1 i l T -).
As before, the elimination process gives rise only to the cancellation of boundary terms, without any fill-in. This back-substitution process can be easily repeated, until each external current is expressed in function only of internal currents and collision sources. In fact when treating the j-th equation in the backward procedure, one can use the following situation:

( ) ( ) ( , ) ( ) ( ) , , , ( ) ( ) ( ) ,
, ,

j j jk k j j i k j i i i i back i back T T E j q for i j T E q ψ ψ ψ ψ ψ ∂ + - ∂ ∂ + > ∂ + - = + + > = + ∑ (A9)
where one can note that all equations beyond the j-th one have already been diagonalised. The backsubstitution of these equations into the j-th one gives:

( ) ( ) , , ( ) ( ) ( , ) ( ) ( ) ( , ) ( ) , , , , , j
j j back j back j j j i i j j i i j i i b a c k j i i b a c k i j i j T E T T T E T E Vq ψ ψ ∂ + ∂ - ∂ > > ⎡ ⎤ ⎡ = + + + ⎢ ⎥ ⎢ ⎣ ⎦ ⎣ ∑ ∑ ⎤ ⎥ ⎦
, which gives the recursive definition of the back matrices. Finally the system (A4) is reduced to a system of the form:

( ) ( ) ( ) , i i i back back T E q ψ ψ ∂ + - = + , ( A10
)
where all external elements have disappeared. Eqs (A10) are substituted into (A4) to completely eliminate the boundary currents from the problem and finally obtain:

, ( ( ) ( i back back T T T T E E ) ) q ψ ψ + +∂ - + ∂ ∂ - ∂ ∂ ∂ = + + + ∑ ∑ ,
where the sums over is done on boundary currents. The balance equations have also to be consistently modified. In fact the original balance equations (15a), written without the hat sign, can be split in the following manner:

∂ H GVq H H GV φ ψ βψ ψ - + - ∂ ∂ = + = + + q ) , ( A11
)
where H ∂ stands for the boundary part of the balance operator H. By making use of (A10) one obtains:

( ) ( ˆĤ T H G H E Vq φ β ψ β - ∂ ∂ = + + + , if , if . j j j I j I j j I j I T T E E q J J T A E q I ρυ ρυ υ ρυ ρυ υ α β α β β α α β υ υ ρ ρ α α ρυ υ ρυ υ ρυ υ α β β α α α β υ υ υ ψ α ψ ψ α + - - - - + - + - + - - - - + ∈∂ ∈ ∈∂ ∈ ⎧ - + - ⎪ - = ⎨ - + ∈∂ ⎪ ⎩ ∑ ∑ ∑ ∑ ∑ I ∉∂ (B1)
From equation (B1) we derive the following expression for the leakage term of the balance equation:

( ) ( ) ( ) ( ) ( ) , , , , , , , , , , , 1 1
I i I ij I I i I i i i I I I j j j i i I i E J J T T T A E E E q ρυ β ρυ ρ ρ α α ρυ ρ υ ρ υ ρ υ α α α β υ α β α β α β α β α β υ α ρυ ρυ ρυ α α α α α α δ δ δ δ δ ψ + - + - - + - + + - ∉∂ ∈∂ ∈ ∈ ∂ ∈ ∉∂ ∈∂ ∈ ∈ - ⎧ ⎫ ⎡ ⎤ = - + - ⎨ ⎬ ⎣ ⎦ Σ Σ ⎩ ⎭ ⎧ ⎫ ⎡ ⎤ + - + ⎨ ⎬ ⎣ ⎦ Σ ⎩ ⎭ ∑ ∑ ∑ ∑ , . j j I υ υ ∈ ∑ - (B2)
Equation (B2) implicitly defines the matrix I of equation [START_REF] Suslov | An algebraic collapsing acceleration method for the acceleration of the inner (scattering) iterations in long characteristics transport theory[END_REF], and another auxiliary matrix E I used to construct the matrices G and M of equation [START_REF] Warsa | Krylov Subspace Iterations for Deterministic k-Eigenvalue Calculations[END_REF]. The incoming matrix I can be better re-written in function of the escape coefficients [START_REF] Sanchez | Treatment of boundary condition in trajectory based deterministic transport methods[END_REF]:

( ) ( ) ( ) ( ) , , , , , , , 1 1 ,
I I i i i I I i i i I T T T A s s T T T A s s E ρυ ρυ ρυ ρυ ρυ α α α β β α β α β α β α α , ρ υ ρυ ρυ ρυ ρυ ρ υ ρυ α α α β β α β α β α α β α δ δ δ δ δ δ - + - - - + - - + - + + + - + + ∉∂ ∈∂ ∈ ∉∂ ∈∂ ∈ ⎡ ⎤ = - - + - = ⎣ ⎦ Σ ⎡ ⎤ - - + - = ⎣ ⎦ Σ ∑ ∑ - (B3)
where [START_REF] Sanchez | New development in the Interface-Current Module TDT[END_REF] has been used together with the fact that A s s A 13), one obtains the following expression for the matrix M I of equation ( 15) :

( ) ( ) , [ j ] I i ij I ij s j ij ij i B V M V E ρυ ρυ ρυ υ ρυ δ δ δ = + -Σ Σ , (B4)
where i,jOEI. Finally, the matrix G that appears in ( 15) is given by: ,

I I i BV G = - Σ E (B4a)
where E I is defined in (B1) and B in [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. Equations (B3) and ( 10) allow us to show that the matrix T+EΣ s Î appearing in ( 16) is sign-symmetric that is to say that it satisfies the relation ( ) ( )

ˆt s s T E I s T E I s + Σ = + Σ .
Indeed, first note that (B3) and ( 10) can be written in matrix form as t and . Moreover let us suppose that (and then also

T sTs = t E sI = -s 1 I I M - Σ 1 I I M - Σ ) is sign-symmetric. This implies that ( ) ( 1 t s I s I T E M I s T E M I s - + Σ = + Σ ) 1 - , (B5)
which is the expected result. We still need to demonstrate Interchanging i with j in the left hand side amounts to interchanging α with β in the right hand side.

Using [START_REF] Sanchez | New development in the Interface-Current Module TDT[END_REF], one then sees that each term in (B6) is transformed either into itself multiplied by s ρ s υ , as it is the case for the terms into the third line of (B6), either in other terms of the same row multiplied by s ρ s υ , as it is the case for the terms in the first two lines. We will find this result in another way in formula (B8). Equations (B4) and (B8) can be used to reduce the computational effort as well as the memory storage, since they show that only halves of the matrices E I and T+EΣ s Î have to be computed and stored.

Thanks to (B6) an analytical definition of operator G can be given:

( ) ( ) ( , ) ( ) ( , ) 1
( ) ( ) 4 lm s i i s j ij ij s l m l T r V j i V m T r V G E d dr n A A E e E i τ ρυ ρυ ρ υ τ τ π - ∈ Ω ∩ → ∂ ∈ Ω ∩ = - = Ω Ω ⋅ Ω Ω ≠ ∑ ∫ ∫ j.(B7)
In equation (B7) the symbol (j→i) stands for the angular directions going from region j to region i, the symbol T(r,Ω) stands for the set of all the chords intersected by the trajectories that pass through the point r with angle Ω, so that T(r,Ω)∩V i (T(r,Ω)∩V j ) stands for the subset of chords that belong to region i (or region j) and symbol E is defined in (4). Therefore the sum in (B7) takes into accounts the case of the intersections with non-convex regions. The optical chord lengths τ l ,τ m and τ lm are related respectively to the chord in region i, to the chord in region j and to the distance between the exiting point of segment l and the entering point of segment m. An analogous expression can be written for the diagonal terms:

( ) ( ) (2 ) ( , ) , ( , ) 1 4 {
[ ( )]/ ( ) ( ) }, i lm s i s i i ii ii s i V l i l m l T r V l m T r V l m B V G E d r d nA A l E E e E ρυ ρυ ρυ ρ υ π τ π τ τ + ∂ - ∈ Ω ∩ ∈ Ω ∩ ≠ = - = ΩΩ⋅ Ω Ω Σ - Σ + ∫ ∫ ∑ ∑ i τ (B7a)
where the second sum inside parentheses comes again from non-convex regions. Inspecting (B7) and (B7a), one can readily realize the symmetric role of indexes i and j as well as of the angular indexes. This assures us

that ij ij ji G G s G s ρυ υ ρ ρ ρ υ = = υ . (B8)
Next, we show that the escape coefficient for the source in Eq. ( 9) remains bounded at the vacuum limit. From the expression for i E ρυ α we obtain where we used ( )

2 2 0 1 lim 2 i i i R i i i i R e R O -Σ Σ → Σ - → - + Σ Σ
i . Formula (B9) shows that the escape coefficients remain bounded, but unfortunately does not give any complementary relation. This means that, in case of nearly voided media a new integral, given by (B9), must be computed.

Another term that must be treated specifically at the vacuum limit is the diagonal term of operator G in equation (B3a). One can show that this diagonal term has the following limit that is to be computed explicitly for the nearly voided cases.

Figure 1: Sketch for the transport equation [START_REF] Sanchez | A Synthetic acceleration for a two dimensional characteristic method in unstructured meshes[END_REF]. The macro-region "I" contains many micro regions "r". For a trajectory along direction Ω, many intersections have to be considered from the entering to the exiting point. In order to express each exiting angular flux along this trajectory, the exponential of the optical length between these regions must be computed. For example, to express thanks to equation [START_REF] Sanchez | A Synthetic acceleration for a two dimensional characteristic method in unstructured meshes[END_REF], a sum over the preceding collision sources (weighted with and ) and over the entering flux must be done. Table 3: Results for the Atrium assembly calculation using the standard DP 1 operator and the new integral IDP 1 operator with double tracking. All cases are DP N -initialized eigenvalue calculations and external iterations are accelerated with the external synthetic algorithm. The computational times are expressed in seconds. The reference k-effective is 1.13811.

Figure 1: Sketch for the transport equation [START_REF] Sanchez | A Synthetic acceleration for a two dimensional characteristic method in unstructured meshes[END_REF]. Macro-region "I" contains many micro regions "r". For a trajectory along direction Ω many intersections have to be considered from the entering to the exiting point.

In order to express each exiting angular flux along this trajectory the exponential of the optical length between these regions must be computed. For example to express thank to equation ( 7) a sum over the preceding collision sources, weighted with and , and over the entering flux must be done.

  between two consecutive transport iterates. Moreover, λ is the new eigenvalue obtained as the result of the multigroup correction problem (20). Problem (20) is an eigenvalue problem but with the presence of an external inhomogeneous term. An accelerated external iteration consists of a free MOC external iteration followed by one synthetic multigroup solution correction

  The remaining step is to complete the inversion of the original by substituting back, from the last until the first, the new computed equations. Let us call N bd the total number of boundary equations in (A4). After the N bd steps described above, one has an upper diagonal operator where the last equation gives the N bd -th external current in function only of internal unknowns and collision sources. The backsubstitution of this equation into the (N bd -1)-th one allows again to express the (N bd -1)-th external current in function only of internal currents and collision sources.

  coefficients of the incoming matrix I and escape coefficients E. Re-inserting (B1) into the balance equation (

FIGURE 2

 2 FIGURE 2

Figure 2 :

 2 Figure 2: Computational mesh for the BWR atrium assembly.

Figure 2 :

 2 Figure 2: Computational mesh for the BWR Atrium assembly.

TABLE 1

 1

	trajectories	tracks	azimuthal angles	polar angles	Transversal
					Integration step
	162	643,617	16 (uniform)	2 (Bickley)	0.025cm
	143	162,610	10 (uniform)	2 (Bickley)	0.05 cm

Table 1 :

 1 Tracking parameters for the calculation of the Atrium assembly with cyclic trajectories. The number of tracks is the total number of intersections of trajectories with regions. The second and the third rows gives, respectively, the data of the reference tracking and of the DP N low-order one.

			TABLE 2
		building time Total time Nb. External it. Nb. Internal it.
	DP 1 Consistent tracking	22.7	41.4	7	151
	IDP 1 Consistent tracking	12.8	26.4	5	102
	IDP 1 Inconsistent tracking	9.44	23.6	5	107
	Free Iterations		412	13	3848

Table 2 :

 2 Results for the Atrium assembly calculation using the standard DP 1 operator and the new integral IDP 1 operator. All cases are DP N -initialized eigenvalue calculations. The computational times are expressed in seconds. The reference k-effective is 1.13811.

TABLE 3

 3

		building time Total time Nb. External it. Nb. Internal it.
	DP 1	19.5 sec	40.3	6	96
	IDP 1	9.4 sec	22.3	3	48
	Free Iterations		412	13	3848

APPENDIX B: MATRIX COEFFICIENTS OF THE INTEGRAL DP N METHOD AND ASYMPTOTIC LIMITS OF ESCAPE PROBABILITIES

In this appendix we give some numerical details about the DP N equations. First of all we detail the expression of matrices M, G and H appearing in [START_REF] Warsa | Krylov Subspace Iterations for Deterministic k-Eigenvalue Calculations[END_REF]. Next we consider the vacuum limit for the escape coefficients in equation [START_REF] Santandrea | Analysis and improvements of the DP N acceleration technique for the method of characteristics in unstructured meshes[END_REF]. Substituting ( 8) into (13) yields:


Adenot

Dr Sercombe Jérôme

Jérôme Sercombe
email: sercombe@cea.fr

Romain Vidal

Christophe Gallé

Frédéric Ade

F Adenot

Experimental study of gas diffusion in cement paste

Keywords:

Experimental study of gas diffusion in cement paste

Introduction

Durability of concrete structures often depends on the rate of ingress of gaseous or aqueous elements into the concrete porous network. Migration of gaseous elements in concrete can be driven by pressure or concentration gradients, i.e., by viscous flow or diffusion. Carbonation is a typical example where gas diffusion, i.e., carbon dioxyde, plays an important role with respect to the overall kinetics of the deterioration process [START_REF] Saetta | The carbonation of concrete and the mechanism of moisture, heat and carbon dioxyde flow through porous materials[END_REF][START_REF] Chaussadent | State of the art and considerations about the carbonation of reinforced concrete[END_REF][START_REF] Bary | Coupled moisture-carbon dioxyde-calcium transfer model for carbonation of concrete[END_REF]]. Drying is another example where the kinetics of gas diffusion, i.e., vapour water diffusion, can be of some importance [4,[START_REF] Xi | Moisture diffusion in cementitious materials[END_REF][START_REF] Mainguy | Modelling of isothermal hydric transfers in porous media. Application to the drying of cement-based materials[END_REF]. In consequence, many mathematical descriptions used to estimate the durability of concrete structures depend on the diffusion coefficient of gas species in the porous network of concrete.

Gas diffusion through porous media is often divided in three independent modes or mechanisms with a distinct diffusion coefficient for each of them [START_REF] Mason | Gas transport in porous media : the dusty-gas model[END_REF]: free-molecule or Knudsen diffusion, molecular or ordinary diffusion and surface diffusion. Usually surface diffusion is neglected, mostly because it is basically an assumption which does not rest on experimental evidence. Theoretically, ordinary diffusion, in which the different species of a mixture move relative to each other under the influence of concentration gradients, occurs predominantly when molecule-molecule collisions dominate over molecule-pore wall collisions. On the contrary, Knudsen diffusion, in which molecules of different species move entirely independently of each other, occurs predominantly when molecule-molecule collisions can be ignored compared to molecule-pore wall collisions.

The prevalence of Knudsen or ordinary diffusion depends, first, on the mean free path of the gas molecule (approximately 0.1 µm for a gas molecule at atmospheric pressure and 20 • C) which depends itself on the total pressure, the temperature of the gas mixture and on the molecular weight of the gas species, second, on the size (diameter if the pores are assumed to be cylindrical) and degree of connectivity of the unsaturated pores (accessible to gas). Since the pore diameters in cement pastes, mortars or concretes are widely distributed from the nanometer to the millimeter scale, it is difficult to state, a priori, which is the dominant mechanism in cementitious materials.

Furthermore, the pore size distribution and degree of connectivity of the unsaturated pores in cementitious materials strongly rely on the mix-design properties (w/c ratio, type of cement, proportion of aggregates,...) and on the curing/conservation conditions (Relative Humidity, temperature) which modify the moisture content and saturation level (i.e., the percentage of the total pore volume filled with water) of the materials. In this respect, some studies have pointed out the importance of moisture content, w/c ratio, aggregate content, type of cement, curing/conservation conditions on the diffusion of gas through mortars or concretes [START_REF] Kobayashi | Oxygen diffusivity of various cementitious materials[END_REF][START_REF] Ohama | Pore size distribution and oxygen diffusion resistance of polymer-modified mortars[END_REF][START_REF] Sharif | Determination of the effective diffusion coefficient in concrete via a gas diffusion technique[END_REF][START_REF] Klink | Characterization of the gas transport properties of porous materials by determining the radon diffusion coefficient[END_REF].

To the authors knowledge, the prevalence of Knudsen or ordinary diffusion in cementitious materials has never been assessed clearly in spite of its importance with respect to the mathematical description of diffusion-based phenomena (carbonation, drying, ...). In this paper, the impact of the moisture content, the total porosity and pore size distribution, of the total gas pressure and of the gas molecular weight on gas diffusion in pure cement pastes is therefore investigated experimentally. Con-clusions are then given concerning the prevalence of Knudsen or ordinary diffusion at low to moderate pressures.

Modes of gas transport in porous media

Usually, four modes of gas transport can be considered in porous media [START_REF] Mason | Gas transport in porous media : the dusty-gas model[END_REF], as illustrated schematically in Figure [START_REF] Saetta | The carbonation of concrete and the mechanism of moisture, heat and carbon dioxyde flow through porous materials[END_REF]. Three of them are related to concentration or partial pressure gradients (ordinary diffusion, Knudsen diffusion, surface diffusion), and one to the total gas pressure gradient (viscous or bulk flow). In the discussion which follows, no total pressure gradient (no bulk flow) is considered since this is the condition which prevails in the experiments presented in this paper. Surface diffusion is neglected since its contribution to the overall transport cannot be assessed precisely. Binary gas mixtures only are considered.

Molecular or ordinary diffusion

Molecular diffusion defines the mechanism by which the different species of a mixture move relative to each other under the influence of concentration gradients, and where molecule-molecule collisions dominate over molecule-pore wall collisions, see Figure [START_REF] Saetta | The carbonation of concrete and the mechanism of moisture, heat and carbon dioxyde flow through porous materials[END_REF]. This case is encountered when the mean free path λ of gas molecules (the mean free path of a gas molecule at atmospheric pressure is approximately 0.1 µm) is much smaller than the characteristic length scale of the pores (the pore diameter if the pores are assumed to be cylindrical). For a binary mixture, in which there is no pressure gradient, the purely diffusive fluxes J 1D and J 2D in a porous material can be written as :

J 1D = -D 12 ∇c 1 (1)
J 2D = -D 21 ∇c 2 (2)
where D 12 and D 21 are the gas diffusion coefficients, ∇c 1 and ∇c 2 the concentration gradients of gas species 1 and 2 in the pore volume, such that D 12 = D 21 and ∇c 1 + ∇c 2 = 0. The measurable binary gas diffusion coefficients D 12 and D 21 through the porous medium are usually related to the free gas diffusion coefficient

D f ree 12 = D f ree

21

of the gas species in open space as follows [START_REF] Mason | Gas transport in porous media : the dusty-gas model[END_REF] :

D ij = g τ D f ree ij (3)
where g and τ are respectively the percentage of open pores in which ordinary diffusion can take place (unsaturated and interconnected pores) and the tortuosity, i.e., the average length of the diffusion path in the material. Complex relations based on pore geometrical arguments have been derived for the tortuosity factor, which usually depends on the total porosity and the saturation level of the porous medium [START_REF] Carman | Flow of gases through porous media[END_REF][START_REF] Millington | Gas diffusion in porous media[END_REF].

The free gas diffusion coefficient D f ree 12 in a binary gas mixture at low to moderate pressures is well known and results from the molecular theory of gases. The following approximate expression can be obtained in many textbook on the subject [START_REF] Reid | The properties of gases and liquids[END_REF][START_REF] Gosse | Transport properties of gases at moderate pressures[END_REF]:

D f ree 12 = 1.858 × 10 -7 1 M 1 + 1 M 2 1/2 T 3/2 P σ 2 12 Ω D (4)
with T the absolute temperature (K), P the pressure of the gas mixture (bars), M i the molecular weight of gas i (g/mol), σ 12 the characteristic length ( Å) and Ω D , the diffusion collision integral (-). Ω D can be estimated according to the formulae given in the appendix. Expression (4) shows that the free gas diffusion coefficient in a binary mixture is inversely proportional to the pressure of the gas mixture, i.e., D f ree 12 = f (1/P ). According to eq. ( 3), a similar dependency can be expected for the diffusion coefficient D ij through the porous medium in case ordinary diffusion is predominant.

Free-or Knudsen diffusion

In the free-molecule or Knudsen diffusion mechanism, the gas molecules collide more frequently with the pore walls than with other gas molecules, see Figure [START_REF] Saetta | The carbonation of concrete and the mechanism of moisture, heat and carbon dioxyde flow through porous materials[END_REF]. This case is encountered when the mean free path λ of the gas molecules is of the order of the characteristic length scale of the pores (diameter D if cylindrical). Since molecule-molecule interactions are negligible, the Knudsen diffusion flux J iD is independent of the gas composition and is therefore defined for each gas species i by the following relation [START_REF] Mason | Gas transport in porous media : the dusty-gas model[END_REF] :

J iK = -D iK ∇c i (5)
with D iK as the Knudsen diffusion coefficient of gas species i in the porous medium, ∇c i the concentration gradient of gas species i. Due to the complexity of the geometry of most porous network, attempts have been made to related the macroscopic Knudsen diffusion coefficient D iK to that of a single cylindrical pore of average radius r (mean radius of the pores in which diffusion can take place, assumed cylindrical and interconnected) through a relation similar to eq. ( 3) :

D iK = g τ D pore iK (6)
with D pore iK the Knudsen diffusion coefficient through a single pore of radius r and g /τ the porosity/tortuosity factor as previously defined. For a long cylindrical pore, D pore iK can be expressed as follows in function of the mean pore radius r and the gas species molecular weight M i [START_REF] Mason | Gas transport in porous media : the dusty-gas model[END_REF][START_REF] Carman | Flow of gases through porous media[END_REF]:

D pore iK = 2r 3 8RT πM i (7)
Note that D pore iK is independent of the pressure of the gas mixture and inversely proportional to √ M i . According to eq. ( 6), this should also be the case for the Knudsen diffusion coefficient D iK of gas species i in the porous medium.

Combination of diffusion mechanisms for binary mixtures

The combination of ordinary and Knudsen diffusion for a binary mixture of gas species i and j is based on momentum-transfer arguments [START_REF] Mason | Gas transport in porous media : the dusty-gas model[END_REF]. According to Newton's second law of motion, if a gas species i is not accelerated on the average, then the momentum transferred to it by collisions with the walls and with species j must be balanced by the external force acting on it, i.e., the gradient of the partial pressure of the species i, ∇p i . Equations ( 1) and ( 5) might be rewritten as follows :

(∇p i ) molecule = RT D ij J iD (8)
(∇p i ) wall = RT D iK J iK (9)
Combining ( 8) and ( 9) with the equilibrium condition ∇p i = (∇p i ) molecule +(∇p i ) wall , the following expression for the combined Knudsen-ordinary diffusion of one component i of a binary mixture at uniform total pressure is obtained :

(∇p i ) = RT D ij J iD + RT D iK J iK (10)
From eq. ( 10), the total flux of species i, J i = J iD + J iK can now be written as :

J i = -D i ∇c i (11)
with the combined Knudsen-ordinary diffusion coefficient D i of species i in the porous material given by :

D i = 1 D iK + 1 D ij -1 = g τ   1 D pore iK + 1 D f ree ij   -1 (12)
Relation [START_REF] Reinecke | Knudsen diffusion, gas permeability, and water content in an unconsolidated porous medium[END_REF] refers to a parallel model, as indicated in reference [START_REF] Carman | Flow of gases through porous media[END_REF]. It differs from the approach adopted by Klinkenberg [START_REF] Klinkenberg | The permeability of porous media to liquids and gases[END_REF] to combine Knudsen diffusion and viscous flow [START_REF] Reinecke | Knudsen diffusion, gas permeability, and water content in an unconsolidated porous medium[END_REF] (serie model). It is of some interest to note that, according to the abovementioned pressure dependency of D pore iK and D f ree ij , the combined Knudsen-ordinary diffusion coefficient D i tends towards the Knudsen diffusion coefficient D iK when the total pressure of the binary gas mixture tends towards 0. This is consistent with the fact that the mean free path tends to infinity when pressure drops to zero. The mean free path becomes therefore much greater than the characteristic size of the pores. Molecule-pore wall collisions are thus predominant. After this period, they were kept for almost 10 years in sealed dessicators where the Relative Humidity (RH) was controled by saturated salt solutions, as indicated in Table 1. To study the diffusion properties of the "dry" materials, one probe of each hardened cement paste was also kept in sealed dessicators containing silica gel (RH ≈ 3%). In 2003, the central part of the probes were cut in 20 mm thick slices, three per probe. The end parts of the probes, 10 mm thick, were also kept in order to estimate the total porosity and the density of the cement pastes. The porosity has been measured by complete drying at 60 • C and re-saturation of the samples at atmospheric pressure. The main physical properties of the studied cement pastes are given in Table [START_REF] Chaussadent | State of the art and considerations about the carbonation of reinforced concrete[END_REF]. After the sawing operation and before testing, the 10 and 20 mm thick samples have been re-conditioned in the same atmospheric conditions as before sawing.

Experimental program

Experimental setup

The main element of the experimental set-up is the diffusion cell, as shown in Fig.

( A diffusion test comprises the following steps. The vacuum is first made in the two compartments of the diffusion cell and in the gas flow lines. One compartment is then filled with pure hydrogen or xenon, the other with pure nitrogen. Periodic analyses of the gas mixture in the nitrogen compartment are then realized. Since an analysis results in a small gas expansion due to the non-negligible volume of the gas flow lines situated between the nitrogen compartment and the chromatograph, the nitrogen compartment is immediately re-filled with pure nitrogen after each measure. This ensures that the gas pressure in the two compartments remains equal during the test. Gas analyses are performed till the steady diffusion state is reached.

The duration of a diffusion test depends on the initial water content of the cement paste sample. For saturated samples, it can last several months. For dry samples, the steady state is usually reached within a few minutes. The mass of the samples are measured before and after each diffusion test in order to check the constancy of the cement paste water content.

Two VARIAN 3400 chromatographs are used for the gas analyses : one with argon as the carrier gas (for the nitrogen-hydrogen gas mixture) and one with helium (for the nitrogen-xenon gas mixture). The change in carrier gas is made to increase the sensitivity of the measures. A molecular sieve is used to remove moisture from the gas mixtures before passing through the capillary columns of the chromatographs for separation. The operating temperature of the chromatographic ovens is of 180 • C, which allows to perform analyses in less than 10 minutes. The chromatographs use thermal conductivity detectors to identify and quantify the elements of the gas mixture. Since the identification is not limited to nitrogen, hydrogen and xenon, any lack of air tightness in the diffusion cell is immediately seen through the detection of small quantities of oxygen. For the gas species of interest (hydrogen, xenon), a precise calibration of the thermal conductivity detectors is made prior to the tests which allows for the quantification of very small percentage in the gas mixture (as low as 0.01 %).

Measurement of the diffusion coefficient

Gas analyses in the nitrogen compartment are performed regularly till the steady diffusive flow is obtained. The analyses give access to the evolution in time of the quantity of hydrogen (or xenon) that has passed through the cement paste sample,

Q i (t) with i = H 2 or Xe.
Considering one-dimensional diffusion through a cement paste sample of thickness L, and using the steady-state assumption, equation [START_REF] Klink | Characterization of the gas transport properties of porous materials by determining the radon diffusion coefficient[END_REF] for the total diffusive flux of a gas species i can be re-cast as follows :

J i = -D i c u i (t) -c d i (t) L ( 13
)
with c u i (t) and c d i (t) respectively as the concentration of gas species i in the upstream and downstream compartments of the diffusion cell. At the beginning of the diffusion test, the upstream compartment is filled with the gas species i (H 2 or Xe), and hence c u i (t = 0) = c • i . The downstream compartment is filled with nitrogen and hence c d i (t = 0) = 0. The test is then performed till steady state is reached. This is usually achieved for small concentrations of gas species i (less than 5%) in the downstream compartment of the diffusion cell. Hence the change in concentration gradient during the test can be neglected :

c u i (t) -c d i (t) ≈ c u i (t = 0) -c d i (t = 0) = c • i ( 14
)
Integration in time of equation [START_REF] Carman | Flow of gases through porous media[END_REF] with the condition (14) leads to the following expression for the quantity of gas species i that has passed through the cement paste sample :

Q i (t) = Sc • i L D i t (15)
with S the surface of cement paste available to gas diffusion. By plotting Q i ×L/S/c • i in function of time as obtained from the periodic gas analyses, it is then possible to estimate the combined Knudsen-ordinary diffusion coefficient D i , as shown in Figure (4). Linear regression is used to best-fit the experimental data. Only the last 5 -6

points are used since the first measures are characteristic of the transient regime.

This method is valid only for gases which are inert with respect to the cement paste. This is the case of hydrogen, xenon, and nitrogen. and hygral condition (RH). The shape of the curves obtained for the two cement pastes are similar. First, they present a plateau for RHs between 0 and 55% (mean diffusion coefficient of about 10 -6 m 2 /s). The plateau is followed by a considerable decrease of the diffusion coefficient till 82%-93% RH (mean value of about 10 -10 m 2 /s). The measures obtained on fully saturated samples (RH 100%) indicate that the diffusion coefficient decreases again when the RH exceeds 93%, leading values between 10 -11 and 10 -13 m 2 /s. A higher discrepancy in the experimental results has been observed for samples kept at high RHs (82-100%). This higher discrepancy results probably from the increasing importance of any heterogenity on diffusion when the pore network is close to full saturation.

Experimental results and discussion

The porosity of OPC pastes (CEM I) can usually be split up in two major components : 1) the capillary porosity (pore diameters in the 3 -100 nm range), which depends strongly on the w/c ratio and 2) the microporosity, associated with the C-S-H gel. The C-S-H gel porosity, being closely related to the C-S-H internal structure, is known to vary little with the w/c ratio [START_REF] Baroghel-Bouny | Characterization of cement pastes and concretes : methods, analysis, interpretations[END_REF][START_REF] Powers | Studies of the physical properties of hardened Portland cement paste[END_REF]. According to water vapour desorption test results [START_REF] Baroghel-Bouny | Characterization of cement pastes and concretes : methods, analysis, interpretations[END_REF], the capillary porosity is first desaturated when the RH decreases from 100% till 40-50%. For lower RHs, it is essentially the micropores which are desaturated [START_REF] Baroghel-Bouny | Characterization of cement pastes and concretes : methods, analysis, interpretations[END_REF][START_REF] Powers | Studies of the physical properties of hardened Portland cement paste[END_REF]. As shown by the curves of Figure ( 5), the impact of the w/c ratio on diffusion appears to be particularly important for RHs greater than 55% : in this domain, the diffusion coefficient of the CEM I, w/c=0.45 cement paste, is 2 to 8 times greater than that of the CEM I, w/c=0.35 cement paste. This difference can be attributed to the increase of the capillary porosity with the w/c ratio. On the contrary, the w/c ratio has a very limited impact on diffusion when the RH is less than 55%, see Figure [START_REF] Xi | Moisture diffusion in cementitious materials[END_REF]. This result is consistent with the fact that it is mainly the C-S-H pores which are water-desaturated at low RHs.

Surprisingly, the desaturation of C-S-H pores, which represent a considerable part of the total porosity (usually around 50%), does not lead to an increase of hydrogen diffivity in the cement paste, as shown by the plateau obtained for RHs in the range 3-55%, see Figure [START_REF] Xi | Moisture diffusion in cementitious materials[END_REF]. This result can possibly be explained by the small size of the C-S-H pores (1 -3 nanometers in diameter) which may make them improper to gas diffusion, or by the important porous network that is already available to gas diffusion as a result of the desaturation of capillary pores.

The impact of the pore size distribution on hydrogen diffusion can further be stretched by plotting the results obtained for the CEM I (OPC) and CEM V (BFS-PFA) cement paste samples. It should be recalled that, in spite of their relatively close total porosities (see Table 2), the pore spaces of the two cement pastes are known, according to microstructural studies, permeability tests [START_REF] Gallé | Gas permeability of unsaturated cement-based materials: application of a multi-scale network model[END_REF][START_REF] Gallé | Effect of drying on cement-based materials pore structure as identified by mercury intrusion porosimetry. A comparative study between oven-, vacuum-, and freeze-drying[END_REF], and radionuclides diffusion tests [START_REF] Richet | Migration of radioelements within cementitious materials. Impact of leaching on transport mechanism[END_REF], to be very different. In particular, MIP (Mercury Intrusion Porosimetry) test results show that the percolation threshold for CEM V cement pastes is found for much smaller pore access diameters than for CEM I cement pastes. This indicates that the pore connectivity of CEM V materials is not as good as in CEM I materials (increased tortuosity). In particular, the capillary porosity of mature CEM V cement pastes is usually thought to be markedly discontinuous [START_REF] Taylor | Cement chemistry[END_REF].

The hydrogen diffusion coefficients of the w/c 0.45 CEM I and CEM V cement pastes are plotted in Figure [START_REF] Mainguy | Modelling of isothermal hydric transfers in porous media. Application to the drying of cement-based materials[END_REF] The impact of pore connectivity on gas diffusion in cement paste can probably best be seen by comparing the diffusion coefficients obtained on the dry CEM I and CEM V samples, see the 3% RH column in Table [START_REF] Bary | Coupled moisture-carbon dioxyde-calcium transfer model for carbonation of concrete[END_REF]. The diffusion of hydrogen in the dry CEM V cement paste is almost 20 times slower than in the dry CEM I cement paste.

Since the evolution of the diffusion coefficients between the dry state and 55% RH is very small for both cement pastes (RH range where the micropores are desaturated), these results can be attributed to main differences in the capillary pore system of the two cement pastes. They are consistent with the expected discontinuity of capillary pores in the CEM V cement paste. They confirm that pore connectivity has a great impact on gas diffusion in cement pastes.

Impact of the gas molecular weight on gas diffusion in cement paste

To study the impact of the gas molecular weight on gas diffusion in cement pastes, diffusion tests have been performed on the CEM I w/c 0.45 cement paste samples with hydrogen-nitrogen and xenon-nitrogen gas mixtures. Xenon was selected first because of its high molecular weight (131 g/mol) as compared to hydrogen (2 g/mol)

and second because of its inert nature. According to Eq. ( 4), the free diffusion coefficient of xenon in nitrogen at 1 bar is equal to 1.24 10 -5 m 2 /s. That of hydrogen in nitrogen at the same pressure is equal to 7.53 10 -5 m 2 /s (see the appendix for the application of Eq. ( 4) to hydrogen-nitrogen and xenon-nitrogen gas mixtures).

In theory, if ordinary diffusion predominates and if the porosity/tortuosity factor is the same for xenon and hydrogen, then the hydrogen to xenon diffusion coefficient ratio in cement paste equals that of the same gas species in open space :

D Xe D H 2 = D f ree Xe/N 2 D f ree H 2 /N 2 = 1.2410 -5
7.5310 -5 = 0.165 [START_REF] Reid | The properties of gases and liquids[END_REF] Applying in the same manner Eq. ( 7) and ( 8) to xenon and hydrogen, the hydrogen to xenon diffusion coefficient ratio in cement paste, in case Knudsen diffusion predominates, is given by :

D Xe D H 2 = M H 2 M Xe = 2 131
= 0.124 [START_REF] Gosse | Transport properties of gases at moderate pressures[END_REF] Note that the diffusion coefficients ratios of Eqs. ( 16) and ( 17) are very close.

The diffusion tests have been performed first with hydrogen-nitrogen and then with xenon-hydrogen. In a few cases, due to small variations in the water content of the samples, a second hydrogen-nitrogen diffusion test has also been undertaken. The mean experimental xenon and hydrogen diffusion coefficients obtained from the tests are given in Table (4). Between 3 and 6 measures were made to obtain each value, except for the 3% RH which is based on a single measure. Theoretical estimates of the xenon diffusion coefficient in cement paste in case of predominant ordinary or Knudsen diffusion are also given in Table (4), based respectively on Eqs. ( 16) or [START_REF] Gosse | Transport properties of gases at moderate pressures[END_REF] and on the mean hydrogen diffusion coefficients obtained during the tests.

The evolution with RH of the xenon diffusion coefficient in the CEM I cement paste is presented graphically in Figure [START_REF] Mason | Gas transport in porous media : the dusty-gas model[END_REF], together with that of the hydrogen diffusion coefficient and the theoretical estimates in case of predominant Knudsen or ordinary diffusion mechanisms. Figure [START_REF] Mason | Gas transport in porous media : the dusty-gas model[END_REF] shows that the shape of the curves obtained for xenon and hydrogen are very similar. First, they decrease slowly till 70% RH and then sharply. The measured diffusion coefficients for xenon are overall 5 to 20 times smaller than that of hydrogen. In the tests performed on saturated samples (100% HR), no xenon flow has been observed, even after a year. It appears therefore that the diffusion of xenon in cement paste becomes very difficult when the percentage of unsaturated pores tends towards zero.

The 3% RH measure appart, the xenon to hydrogen diffusion coefficients ratio is almost constant in the RH range 55-82%, with a mean value of 0.1, close to the theoretical estimate of Eq. ( 17) in the case of Knudsen diffusion (0.124). This means that the unsaturated porous space available to gas diffusion at RHs less than or equal to 82% is fully accessible to hydrogen and xenon. This is not the case when the RH is greater or equal to 93% since the xenon to hydrogen diffusion coefficients ratio drops respectively to 0.05 and 0 for RHs of 93 and 100%.

These results can possibly be attributed to the formation of discontinuities in the porous network accessible to gas diffusion. When pore humidity is high, the narrow connections (necks) that exist between large pores which are still unsaturated can be filled with liquid water [START_REF] Xi | Moisture diffusion in cementitious materials[END_REF]. Hence, the diffusion of a gas throught nearly saturated materials depends on the solubility and on the diffusion coefficient of the gas in the liquid phase that saturates the pore necks. The solubility of xenon in water is five times greater than that of hydrogen and cannot therefore explain the lack of xenon diffusion through the nearly saturated CEM I cement paste samples. Its important size (the atomic radius of xenon is almost six times greater than that of hydrogen), which would limit its diffusion in saturated narrow pore connections, could however be at the origin of these results.

Impact of the total gas pressure on gas diffusion

The total gas pressure dependency is one of the major difference that exists between ordinary and Knudsen diffusion in porous media. In case of ordinary diffusion, the diffusion coefficient of a gas species in a porous medium is inversely proportional to the total pressure, see Eqs. ( 4) and [START_REF] Xi | Moisture diffusion in cementitious materials[END_REF]. In case of Knudsen diffusion, it does not depend on the total pressure of the gas phase, see Eqs. ( 7) and [START_REF] Kobayashi | Oxygen diffusivity of various cementitious materials[END_REF]. To determine which is the dominant diffusion mechanism in cementitious materials, hydrogen diffusion tests have been performed at four different pressures (0.1, 0.5, 1 and 1.9 bars) on samples of the two CEM I cement pastes in equilibrium with air at 3, 55, 70 and 82% RH. One sample per RH and cement paste type has been tested at the four different pressures in a random order.

The hydrogen diffusion coefficients measured at the four different pressures and on the different cement pastes are given in Table [START_REF] Xi | Moisture diffusion in cementitious materials[END_REF]. As shown by these results, whatever the RH and the cement paste, the hydrogen diffusion coefficient decreases with increasing pressure. This trend is consistent qualitatively with an ordinary diffusion mechanism. However, quantitatively, the evolution of the hydrogen diffusion coefficient is far from being proportional to 1/P . This can be best seen by plotting the D i (P )/D i (1 bar) diffusion coefficient ratio in function of the pressure ratio 1/P (bar/bar), as illustrated in Figure [START_REF] Kobayashi | Oxygen diffusivity of various cementitious materials[END_REF]. In case of a purely ordinary diffusion mechanism, D i (P )/D i (P • ) = P • /P , as shown by the continuous line in Figure [START_REF] Kobayashi | Oxygen diffusivity of various cementitious materials[END_REF]. In case of a purely Knudsen diffusion mechanism, D i (P )/D i (P • ) = 1, as indicated by the dashed line in Figure [START_REF] Kobayashi | Oxygen diffusivity of various cementitious materials[END_REF].

Figure [START_REF] Kobayashi | Oxygen diffusivity of various cementitious materials[END_REF] shows that all the test results (circles and triangles) are close to the pure Knudsen diffusion mechanism line, indicating that the total pressure of the gas phase has little impact on the hydrogen diffusion coefficient in cement paste. Obviously, the pressure dependency of the diffusion mechanism at hand in cement paste does not vary with the hygral state of the cement paste, i.e., its internal RH. This means that it is the same diffusion mechanism which dominates in cement paste when the capillary pores and/or the micropores are accessible to the gas phase. The complex pore connections and tortuosity in cement paste are very probably at the origin of this result. In fact, the diffusion path through cementitious materials can be seen as a succession of meso-and micropores connected by necks which makes therefore the overall gas flow dependent on the size of the smallest pores. In the studied pressure range 0.1-1.9 bars, the dominant diffusion mechanism is therefore of Knudsen type and reflects the importance of the interactions between the gas molecules and the pore structure of the material.

Conclusions

In this paper, experimental results on gas diffusion in three cement pastes have been presented. The main conclusions of this study are :

1. Gas diffusion depends essentially on the saturation level of the cement paste (which determines the percentage of pores available to gas), and hence on the RH of air when the material is in equilibrium with its external environment. The hydrogen diffusion coefficient was found to vary between 10 -6 m 2 /s for dry samples and 10 -13 m 2 /s for saturated samples.

2. The total porosity and the pore size distribution/connectivity of cement paste have a strong impact on hydrogen diffusion. Capillary pores are of great importance with regards to gas diffusion in cementitious materials since the diffusion coefficient decreases strongly in the RH range 55-100%. On the contrary, as shown by the lack of variation of diffusion properties when the RH is smaller than 55%, micropores play a negligible role with regards to gas diffusion in cementitious materials.

3. Gas diffusion in CEM I cement paste depends approximately on the square root of the gas molecular weight, except when the saturation level of the material is high, i.e., when the RH is greater or equal to 93%. This limit is probably due to the discontinuity of the unsaturated pore network at high RH which makes gas diffusion through nearly saturated cementitious materials dependent on the gas solubility and gas diffusion properties in liquid water.

4. In the pressure range 0.1 -1.9 bars, gas diffusion in CEM I cement paste depends slightly on the total pressure of the gas phase. This result indicates that Knudsen diffusion mechanisms are dominant in these materials. The pore size (capillary pores or micropores) and the saturation level have no impact on this tendency. These results need to be confirmed for higher total gas pressures. 1.3 10 -6 1.9 10 -6 6.6 10 -7 4.5 10 -10 1.3 10 -10 Xenon 3.0 10 -7 1.9 10 -7 5.5 10 -8 6.3 10 -11 7.1 10 -12 Xenon Eq. ( 16) 2.2 10 -7 3.2 10 -7 1.1 10 -7 7.4 10 -11 2.2 10 -11 Xenon Eq. ( 17) 1.7 10 -7 2.4 10 -7 8.1 10 -8 5.6 10 -11 1. 7 10 -11 Table 5: Hydrogen diffusion coefficients (m 2 /s) obtained on the CEM I cement paste samples tested at 0.1, 0.5, 1 and 1.9 bars Pressure (bars) 0.1 0.5 1 1.9 CEM I w/c 0.35 HR 3% 2.01 10 -6 1.27 10 -6 1.05 10 -6 8.04 10 -7 CEM I w/c 0.35 HR 55% 1.72 10 -6 1.41 10 -6 1.18 10 -6 9.07 10 -7 CEM I w/c 0.35 HR 70% 7.4 10 -8

List of Tables

List of Figures

7.24 10 -8 5.45 10 -8 4.47 10 -8 CEM I w/c 0.35 HR 82% 2.28 10 -10 1.79 10 -10 1.15 10 -10 7.98 10 -11 CEM I w/c 0.45 HR 3% 2.31 10 -6 1.39 10 -6 1.08 10 -6 1.09 10 -6 CEM I w/c 0.45 HR 55% 3.27 10 -6 2.95 10 -6 2.45 10 -6 1.81 10 -6 CEM I w/c 0.45 HR 70% 1.14 10 -6 9.93 10 -7 8.02 10 -7 6.21 10 -7 CEM I w/c 0.45 HR 82% 3.64 10 -9 3.52 10 -9 2.69 10 -9 2.70 10 -9

2. 1 .

 1 Materials: characteristics, conditioning and sampling The experimental program has been carried out on hardened cement pastes made of French industrial CEM I (OPC) and CEM V (BFS-PFA) cements. Cement pastes were prefered to concretes or mortars since diffusion tests could be performed on relatively thin specimen (20 mm) which can be considered homogeneous. Six cylindrical probes (40 mm diameter and 80 mm high) per cement type (CEM I and V) and w/c ratio (0.35 and 0.45) have been prepared in 1994. After demoulding, they have been cured at 20 • C for 9 months in a Ca(OH) 2 water saturated solution.

  ). It consists of two stainless steel compartments of equal dimensions (internal diameter 60 mm, internal length 90 mm, total volume of about 90 ml), separated by the cement paste sample, glued with an epoxy to two half stainless steel rings.The ring system is designed to ensure the gas tightness on the outer surface of the cement paste sample. To this end, it covers partly the flat faces of the cement paste samples, leaving a diffusion surface of 20 mm in diameter, see Fig. (2). The airtightness of the cell is ensured by rubber O-rings which are in contact with the two steel rings. The experimental set-up used for the diffusion tests, as sketched in Fig. (3), consists of a diffusion cell, vacuum pumps, pressure sensors, gas flow lines and a chromatograph. The vacuum pumps are used at the beginning of the test to eliminate the gas phase initially in the diffusion cell and in the gas flow lines. The pressure sensors are used to control and monitor during the test the pressure of the gas mixture in the two compartments of the diffusion cell. The chromatograph is used to analyze periodically the gas mixture in one of the two compartments.

4. 1 .

 1 Impact of w/c ratio and cement type on hydrogen diffusionAll the diffusion tests presented hereafter have been performed with a hydrogennitrogen gas mixture at a total pressure of 1 bar. The mean hydrogen diffusion coefficients measured on the three cement pastes (CEM I with w/c ratios equal to 0.35 and 0.45, CEM V with a w/c ratio equal to 0.45) pre-conditioned in specific atmospheric conditions (3, 55, 70, 82, 93 and 100% RH) are given in Table[START_REF] Bary | Coupled moisture-carbon dioxyde-calcium transfer model for carbonation of concrete[END_REF].In general, three measures per cement type and hygral condition (RH) have been made. The results show the pronounced evolution of the diffusion coefficient with the moisture content of the samples.To get a better picture of this evolution, the hydrogen diffusion coefficients of the two CEM I cement pastes are plotted in Figure[START_REF] Xi | Moisture diffusion in cementitious materials[END_REF] in function of the RH. The continuous lines of Figure (5) join the mean diffusion coefficients obtained for each cement paste

  in function of the RH. The continuous lines of Figure[START_REF] Mainguy | Modelling of isothermal hydric transfers in porous media. Application to the drying of cement-based materials[END_REF] join the mean diffusion coefficients obtained for each cement paste and hygral condition (RH). The shape of the two curves of Figure[START_REF] Mainguy | Modelling of isothermal hydric transfers in porous media. Application to the drying of cement-based materials[END_REF] are close but the diffusion coefficient drop in the case of the CEM V curve seems translated towards higher RHs. Both curves present a plateau for RHs between 3 and 55%. The mean hydrogen diffusion coefficients for the two materials are however very different at this stage, being close to 10 -6 m 2 /s and 5.10 8 m 2 /s for the CEM I and CEM V cement pastes, respectively.The plateau is followed by a slow decrease of the diffusion coefficient till 70% RH for the CEM I material (mean value of 6.10 -7 m 2 /s) and 82% RH for the CEM V (mean value of 4.10 -9 m 2 /s). A rapid decrease of the diffusion coefficients over 3 orders of magnitude is then observed for both materials when the RH increases. The average diffusion coefficients obtained for the two cement pastes at 93% RH are finally very different, being equal to 4.10 -10 m 2 /s for the CEM I and 3.10 -12 m 2 /s for the CEM V. Measures made on saturated samples (100% RH) indicate a further decrease of the hydrogen diffusion coefficient in both materials. In particular, no hydrogen flow has been observed during two of the three 7-months long diffusion tests performed on CEM V saturated samples.

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: Transport processes in porous media

Figure 4 :

 4 Figure 4: Estimation of the Knudsen-ordinary diffusion coefficient from experimental

Figure 5 :

 5 Figure 5: Effective diffusion coefficients versus RH for the CEM I cement paste

Figure 6 :

 6 Figure 6: Effective diffusion coefficients versus RH for the CEM I and CEM V

Figure 7 :

 7 Figure 7: Hydrogen and Xenon effective diffusion coefficients versus RH for the

Figure 8 :

 8 Figure 8: Hydrogen diffusion coefficients ratio versus total pressure ratio for the

  4 10 -11 CEM I w/c 0.45 1.4 10 -6 1.6 10 -6 5.8 10 -7 3.5 10 -10 1.3 10 -10 CEM V w/c 0.45 7.5 10 -8 3.0 10 -8 5.0 10 -9 4.5 10 -9 2.1 10 -12

Figure 2 :Figure 3 :Figure 5 :Figure 6 :

 2356 Figure 2: Sketch of the diffusion cell

17

 17

Figure 7 :

 7 Figure 7: Hydrogen and Xenon diffusion coefficients versus RH for the CEM I cement paste samples with a w/c ratio of 0.45.

Table 1 :

 1 Salt solutions and related Relative Humidities

Table 2 :

 2 Properties of the tested cement paste materials

Table 3 :

 3 Mean hydrogen diffusion coefficients (m 2 /s) obtained on the CEM I and CEM V cement paste samples (total pressure 1 bar).

Table 4 :

 4 Mean xenon and hydrogen diffusion coefficients (m 2 /s) obtained on the CEM I w/c 0.45 cement paste samples and theoretical estimates of the xenon diffusion coefficient in case of predominant ordinary (Eq. 16) or Knudsen (Eq. 17)

	diffusion mechanisms.

Table 5 :

 5 Hydrogen diffusion coefficients (m 2 /s) obtained on the CEM I and CEM V cement paste samples tested at 0.1, 0.5, 1 and 1.9 bars

Table 6 :

 6 Parameters used in Eq. (4) to calculate the free diffusion coefficients of hydrogen and xenon in nitrogen.

Table 1 :

 1 Salt solutions and related Relative Humidities

	RH (T=20 • C)	54.5%	69.9% 81.8% 93.2% 100%
	Salt solution	Mg(NO 3 ) 2 ,H 2 O	KI	KBr KNO 3 H 2 O

Table 2 :

 2 Properties of the tested cement paste materials

	Cement type	CEM I CEM I CEM V
	w/c ratio	0.35	0.45	0.45
	wet density (g/cm 3 )	2.08	2.0	1.92
	porosity (water)	27%	34%	35%

Table 3 :

 3 Mean hydrogen diffusion coefficients (m 2 /s) obtained on the CEM I and CEM V cement paste samples (total pressure 1 bar). /c 0.35 8.6 10 -7 7.8 10 -7 3.8 10 -8 6.1 10 -11 8.

	HR (%)	3	55	70	82	93
	CEM I w

Table 4 :

 4 Mean xenon and hydrogen diffusion coefficients (m 2 /s) obtained on the CEM I w/c 0.45 cement paste samples and theoretical estimates of the xenon diffusion coefficient in case of predominant ordinary (Eq.[START_REF] Reid | The properties of gases and liquids[END_REF] or Knudsen (Eq. 17) diffusion mechanisms.

	HR (%)	3	55	70	82	93
	Hydrogen

Table 6 :

 6 Parameters used in Eq. (4) to calculate the free diffusion coefficients of hydrogen and xenon in nitrogen.

		M i (g/mol) ( Å) σ i	i /k (K)
	H 2	2.016	2.827 59.7
	N 2	28.0	3.568 113.0
	Xe	131.3	4.061 225.3

Appendix

The data used in this paper to calculate the xenon-nitrogen and hydrogen-nitrogen free diffusion coefficients are summarized in Table [START_REF] Mainguy | Modelling of isothermal hydric transfers in porous media. Application to the drying of cement-based materials[END_REF].

The characteristic length σ 12 of Eq. ( 4) has been estimated as follows :

[17] gives the following approximate expression to determine the diffusion collision integral Ω D :

with

and

Molecular diffusion

Viscous flow

Surface diffusion