Arnaud Peschard

Alexandre Govin
email: govin@emse.fr

Philippe Grosseau

Bernard Guilhot

R Guyonnet

Effect of polysaccharides on the hydration of cement paste at early ages

Keywords: . Hydration, Polymers, Polysaccharides, Retardation, Cement Paste

Introduction

Organic admixtures have been widely used in concrete and mortar for several decades, providing advanced specific properties needed during the construction. The European market of coatings is changing to the advantage of the remodelling market. Mortars for this application require a water retention value of at least 95% to prevent water absorption. To obtain this value, polysaccharides which can be equally classified as water-reducers and setretarding agents, are frequently used. These admixtures can be either native bio-polymers (e.g. starch) or substituted polymers (e.g. cellulose ethers). Dewacker et al. [START_REF] Dewacker | Cement mortar systems using blends of polysaccharides and cold-water-soluble unmodified starches[END_REF] suggested that a blend of polysaccharides such as cellulose and starch was a good water retention agent. These admixtures are also set-retarding agents which improve workability duration and modify cement hydration. Yamamuro [START_REF] Yamamuro | Property of new polysaccharide derivative as a viscosity agent for self compacting concrete[END_REF] showed that a polysaccharide derivative containing an ionic functional group and a hydrophobic group increases the viscosity of cement suspensions. In the Shotcrete process, Ghio et al. [START_REF] Ghio | The effects of polysaccharide gum additives on the shotcrete process[END_REF][START_REF] Ghio | The rheology of fresh cement paste containing polysaccharide gums[END_REF] showed that concretes formulated with polysaccharides are easier to pump and to spray at high shear rates than unadmixed concrete. Moreover, concrete is more cohesive and viscous at low shear rates.

Tanaka et al. [START_REF] Tanaka | Cement compositions containing a sulfated polysaccharide and method[END_REF] patented an additive that contains at least one sulfated polysaccharide to improve fluidity and workability and that supplies higher final compressive strength.

Hayakawa et al. [START_REF] Hayakawa | An experimental study on the concrete containing cellulose ether[END_REF][START_REF] Hayakawa | Effects of cellulose ether on bond between matrix and aggregate in concrete[END_REF] reported that cellulose ether improves binding between cement matrix and aggregates.

Previous investigations on physico-chemical aspects were restricted to the effect of monosaccharides or sugar acids on cement or single phase hydration [START_REF] Young | The influence of sugars on the hydration of tricalcium aluminate[END_REF][START_REF] Collepardi | Influence of gluconate, lignosulfonate or glucose on the C 3 A hydration in the presence of gypsum with or without lime[END_REF][START_REF] Collepardi | Influence of gluconate, lignosulfonate, and glucose admixtures on the hydration of tetracalcium aluminoferrite in the presence of gypsum with or without calcium hydroxide[END_REF][START_REF] Milestone | Hydration of tricalcium silicate in the presence of lignosulfonates, glucose and sodium gluconate[END_REF]. Influence of such admixtures on Calcium aluminate phase (C 3 A) was studied by Young [START_REF] Young | The influence of sugars on the hydration of tricalcium aluminate[END_REF] who showed that, interaction of sugars on C 3 A prevents rapid formation of the cubic phase C 3 AH 6 and promotes formation of the hexagonal phase C 4 AH 13 . Collepardi et al. [START_REF] Collepardi | Influence of gluconate, lignosulfonate or glucose on the C 3 A hydration in the presence of gypsum with or without lime[END_REF][START_REF] Collepardi | Influence of gluconate, lignosulfonate, and glucose admixtures on the hydration of tetracalcium aluminoferrite in the presence of gypsum with or without calcium hydroxide[END_REF] found that addition of glucose, gluconate or lignosulfonate stabilizes ettringite in the C 3 A-gypsum system. It was also noted that glucose retards gypsum consumption and ettringite formation.

Many studies on the hydration of silicate phase in presence of admixtures were devoted to the determination of calcium hydroxide content. According to Milestone [START_REF] Milestone | Hydration of tricalcium silicate in the presence of lignosulfonates, glucose and sodium gluconate[END_REF] the retarding effect is as follows: sugar acids > sugars > lignosulfonates. Scanning electron microscopic (SEM) examination showed that in presence of these admixtures, calcium silicate phase (C 3 S) particles are covered by flaky calcium silicate hydrate (CSH), instead of acicular hydrates [START_REF] Milestone | Hydration of tricalcium silicate in the presence of lignosulfonates, glucose and sodium gluconate[END_REF].

Bensted et al. [START_REF] Ramachandran | Conduction calorimetric investigation of the effects of retarders on the hydration of Portland cement[END_REF] investigated the effect of lignosulfonate on the hydration of cement and found that retardation is linked either to a preferential complexation with calcium and silicon or to an introduction of organic admixtures into a precipitated gel. This latter hypothesis postulates that the retarding effect would come from the formation of a less permeable gel coating around cement grains. However the mechanism of cement-admixture interaction is yet to be fully understood.

The present study was undertaken to determine the relative efficiencies of polysaccharides and to reveal the mechanism of interaction between cement particles and the admixture at the early stages of hydration.

Experimental

Raw materials

White Porland cement (C 1 ): CPA CEM I 52.5, and grey cement (C 2 ): PMES 42.5 were used in this study. The designation follows French standard NF P 15-301. Chemical composition and mineralogical phases of the two cement samples were calculated by the Bogue approximation [START_REF] Taylor | Cement Chemistry[END_REF], and these are listed in Table 1.

Five powder polysaccharides, a cellulose ether (CE), a starch ether (SE), a native starch (NS), a white dextrin (WD), and a yellow dextrin (YD) both prepared from NS, were tested. Properties of these polysaccharides are given in Table 2.

Methods of investigation

Isothermal microcalorimetry, TGA and FTIR spectrometry were used to monitor kinetics and mechanism of cement hydration.

Each experiment was conducted in triplicate. In each case, cement and polysaccharides were mixed for 2 periods of 2 minutes using a shaker-mixer (Wab, Turbula, Germany).

Unless otherwise stated, polysaccharide to cement weight ratio (P/C) was equal to 0.5 (w/w).

Deionized water to cement ratio (W/C) was approximately 0. Cement pastes were obtained by using a mixer (Technotest, France) for one minute at low speed (60 rpm) and four minutes at high speed (120 rpm). Samples were stored in a water vapour saturated desiccator and under nitrogen atmosphere. Cement hydration was stopped by breaking and grinding hydrated samples and washing them with anhydrous ethanol. Samples were ground in an agate mortar until each particle diameter became lower than 100 µm. This step was reproduced at different lengths of time ranging from 1 to 24 hours.

Thermal analysis was performed on a TG-DSC 111 thermogravimetric apparatus (Setaram, France). All experiments were conducted from 20°C up to 800°C (heating rate of 10°C/min ) under dynamic argon atmosphere on a approximately 25 mg samples.

The amount of portlandite CH(%) was obtained directly from TG curves by Eq.1

H MW CH MW (%) CH WL (%) CH × = (1)
where WL CH (%) corresponds to weight loss in percent occurring during CH dehydration, MW CH and MW H are molecular weights of portlandite and water respectively.

Infra-red spectra were collected with a spectrometer (Biorad FTS 185 Digilab, USA) equipped with a Harrick scattering reflexion accessory. Analysis was conducted on hydrated samples diluted at 5% in KBr. Spectra were converted by the Kubelka-Munk algorithm.

Admixtures can be classified in retarders or accelerators category using isothermal microcalorimetry [START_REF] Wilding | A classification of inorganic and organic admixtures by conduction calorimetry[END_REF]. It enables the hydration process to be monitored continuously. It was conducted on a C80 calorimeter (Setaram, France). Temperature of the calorimetric block was maintained at 25°C and the experiments were conducted on a 250 mg sample at a W/C ratio of 0.4. Cement and water mixing was done outside the calorimeter according to the procedure described by Sauvat et al. [START_REF] Sauvat | A study of ordinary Portland cement hydration with wood by isothermal calorimetry[END_REF] and Ramachandran et al. [START_REF] Ramachandran | Conduction calorimetric investigation of the effects of retarders on the hydration of Portland cement[END_REF]. Consequently, only the exothermic hydration peak was encountered.

Results and discussion

Hydration of C 1

Investigations were first focused on C 1 since this kind of cement is frequently used in coating formulas. The calorimetric curve of neat cement exhibits two separate exothermic peaks. According to Sauvat et al. [START_REF] Sauvat | A study of ordinary Portland cement hydration with wood by isothermal calorimetry[END_REF] they correspond to C 3 A and C 3 S hydration peaks respectively whereas Bensted et al. [START_REF] Bensted | Some applications of conduction calorimetry to cement hydration[END_REF] suggested the reverse order. Retarding ability of polysaccharides on cement hydration is obtained from curves of heat evolution rate. Three different times were defined as proposed by Sauvat et al. [START_REF] Sauvat | A study of ordinary Portland cement hydration with wood by isothermal calorimetry[END_REF], t 1 , t 2 , t 3 which correspond to minimum flux time, second peak time and third peak time respectively. Heat released in 24 hours of hydration was calculated by curve integration from t 1 to t 24h (Table 3). Intensity and shape of the exothermic peaks are not significantly modified by CE and NS addition. Intensity of heat flow for cement formulated with SE is lower. The cement admixed with YD presents a larger wetting peak which induces a shift of the minimum heat flow to 10 hours. Moreover, the heat released during the first 24 hours of hydration is smaller. These results confirm previous studies regarding effects of set-retarding agents on cement hydration [START_REF] Sauvat | A study of ordinary Portland cement hydration with wood by isothermal calorimetry[END_REF][START_REF] Ramachandran | Conduction calorimetric investigation of the effects of retarders on the hydration of Portland cement[END_REF].

Ramachandran et al. [START_REF] Ramachandran | Conduction calorimetric investigation of the effects of retarders on the hydration of Portland cement[END_REF] found that sucrose is an efficient retarder which extends induction period for several hours. Similar result was obtained with YD in the present study.

Previte et al. [START_REF] Previte | Some insights on the mechanism of saccharide set retardation of Portland cement[END_REF] suggested that alkali-stable and high molecular weight sugars are more efficient to retard hydration. It should be noted that oligomers used by these authors had a maximum chain length of 3 units. To extend this previous study, the present results focus on polysaccharides with a much higher degree of polymerization. Polysaccharides with a strong set-retarding ability i.e., SE and YD have a lower average molecular weight than CE and NS.

A weak average molecular weight favors alkali-degradation of polysaccharides [START_REF] Whistler | Alkaline degradation of polysaccharides[END_REF][START_REF] Knill | Degradation of cellulose under alkaline conditions[END_REF].

Several authors have shown that products from decomposition of sugar are more effective than native sugars [START_REF] Ramachandran | Conduction calorimetric investigation of the effects of retarders on the hydration of Portland cement[END_REF][START_REF] Bruere | Set-retarding effects of sugars in Portland cement pastes[END_REF].

Differential Thermal Analysis (DTA) and TGA were used to study the effects of sugars and hydroxycarboxylic acids on hydration of different cement phases. As Bensted et al.

[21] and Odler et al. [START_REF] Odler | Possibilities of quantitative determination of the Aft-(ettringite) and Afm-(monosulphate) phases in hydrated cement pastes[END_REF] suggested, it is possible to assign the first two peaks at 120°C and 145°C to ettringite and gypsum decomposition respectively (Fig. 1a). Beyond 4 hours of hydration, decomposition peaks of CSH and ettringite overlap. It is also possible to monitor gypsum consumption until it completely disappears since no other phase interferes with it.

Decomposition peak of CH occurs between 420°C and 460°C.

The C 1 +SE formula revealed a smaller amount of ettringite and a smaller gypsum consumption than C 1 (Fig. 1b). However the admixture YD accelerated gypsum consumption and ettringite formation (Fig. 1c) as compared with C 1 . It can be at the origin of the wetting peak broadening on the calorimetric curve. Between 4 and 8 hours of hydration, the ettringite phase seems to be stabilized and gypsum is no longer consumed for both formulas. Moreover, decomposition peak of portlandite is less intense for C 1 +SE and C 1 +YD. Further information is obtained by calculation of portlandite content as specified by Eq. 1. Portlandite formation depends strongly on the type of polysaccharide introduced in cement (Fig. 2). YD stops portlandite formation during 8 hours. Beyond this period the formation rate increases and finally portlandite content is very close to those of pure cement. Conversely, CE and NS have had no significant effect on gypsum consumption, on ettringite as well as on portlandite formation (Fig. 2).

A smaller amount of portlandite could be associated with a greater carbonation as proposed by Silva et al. [START_REF] Silva | Evidences of chemical interaction between EVA and hydrating Portland cement[END_REF] for cement admixed with ethylene/vinyl acetate copolymer.

Weight losses occurring during decarbonation were compared. Carbonation does not significantly increase upon polysaccharide incorporation. Thus, at a given time the presence of admixtures YD and SE induces a smaller content of portlandite.

CSH formation was monitored by infra-red spectroscopy since monitoring by thermal analysis at early ages is not possible. CSH is usually observed by a gradual shift of band centered at 925 cm -1 to 970 cm -1 , which is assigned to the asymmetrical SiO 4 stretching frequency [START_REF] Bensted | Some applications of infrared and raman spectroscopy in cement chemistry, Part 3-Hydration of portland cement and its constituents[END_REF]. It indicates polymerization of silicate units SiO 4 4-with formation of CSH phase [START_REF] Mollah | A Fourier transform infrared spectroscopic investigation of the early hydration of Portland cement and the influence of[END_REF]. CSH has already been formed for C 1 and C 1 +CE at 8 hours of hydration (Fig. 3).

On the contrary, C 1 +SE and C 1 +YD revealed no shift attributed to CSH formation.

Results obtained from these studies are consistent. CE has no significant effect on hydrate formation. YD extends induction period, favours ettringite formation at very early ages and retards CSH and CH formation. SE slows down ettringite formation and gypsum consumption and also delays considerably CSH and CH formation.

Influence of SE and YD concentration on cement hydration

Results obtained by isothermal microcalorimetry for an increasing amount of polysaccharides are listed in Table 4. It shows that t 1 and t 2 are not significantly modified by an increase of SE concentration. Time t 3 is not defined for SE concentration over 0.5%. As YD concentration rises, t 1 and t 2 are shifted to later times. In both formulas, evolution of heat released is inversely proportional to concentration. Results of C 1 +YD are consistent with those of Ramachandran et al. [START_REF] Ramachandran | Conduction calorimetric investigation of the effects of retarders on the hydration of Portland cement[END_REF] who proposed a linear relation between concentration of retarders and duration of induction period.

Study on the admixture YD

Influence of cement composition on retarding ability of YD

In order to determine the effect of cement composition and especially C 3 A content on set-retarding ability of polysaccharides, cement C 2 was studied. For this experiment admixture YD was used, as it has a strong retarding effect on cement hydration. As illustrated in (Fig. 4a), gypsum was totally consumed for pure C 2 after 16 hours of hydration. Introduction of YD revealed a decomposition peak at 145°C attributed to gypsum up to 24 hours (Fig. 4b).

Portlandite formation rate is higher for C 1 than C 2 (Fig. 5). It can be associated with a higher C 3 S content for C 1 . C 3 A amount seems to be an important parameter since 0.5% of YD totally inhibits portlandite formation up to 48 hours. Present results are similar to those of Singh [START_REF] Mollah | A Fourier transform infrared spectroscopic investigation of the early hydration of Portland cement and the influence of[END_REF]. Monosi et al. [START_REF] Monosi | Combined effect of lignosulfonate and carbonate on pure Portland clinker compounds hydration. III. Hydration of tricalcium silicate alone and in the presence of tricalcium aluminate[END_REF] found that lignosulfonate totally stops C 3 S hydration up to two weeks without C 3 A. The inhibiting effect lasts only one day with C 3 A. It reveals the importance of C 3 A content in cement which acts as a sink for retarders [START_REF] Ramachandran | Conduction calorimetric investigation of the effects of retarders on the hydration of Portland cement[END_REF].

C 2 hydrated with 0.5% YD suggests an excess of admixture with respect to cement. P/C 3 A ratio for this formula is 0.2 instead of 0.045 for C 1 +0.5%YD. In order to verify the fact that an excess of YD with respect to C 2 was introduced, an experiment on C 1 +3% YD was carried out (ratio P/C 3 A equal to 0.25). Results obtained with this formula are similar to C 2 +0.5%YD, i.e., portlandite formation was completely inhibited.

Influence of dextrinization on set retarding ability of starch compounds

WD and YD were obtained from NS by dextrinization which is a combination of an acid and a thermal treatment. It allows starch solubilization by a decrease of average molecular weight. Soluble fraction value increases with degree of conversion of starch to dextrin (Table 2). Thermal analysis is performed on C 1 with 0.5% of admixture. YD delays portlandite formation whereas NS and WD have no effect (Fig. 6).

Average molecular weights and soluble fraction of starch derivatives are given in Table 2. At a given P/C ratio, YD concentration in liquid phase is higher than WD and NS concentration. Portlandite rate formation is more dependent on soluble fraction than on molecular weight of compounds. Dextrinization of NS is a predominant parameter on the retarding ability of starch derivatives.

Conclusions

Results obtained by thermal analysis show that the introduction of a small amount of polysaccharide (0.5% by weight of cement) induces a strong modification of cement hydration. Retarding ability is dependent on polysaccharide variety and is as follows: NS < WD < CE < YD < SE.

A mechanism which describes cement-polysaccharide interaction is difficult to establish. The first step of C 3 A hydration leading to ettringite is accelerated by YD and retarded by SE. These polysaccharides then stabilize ettringite and gypsum content. This observation will induce a delay in conversion of ettringite to monosulfoaluminate phase. From these results, the action of polysaccharides on C 3 A hydration seems more restricted to growth of hydrates than to nucleation. Consequently the origin of retardation could be linked to an adsorption of admixtures on the first hydrates forming a less permeable coating. As far as silicate hydration is concerned SE and YD (most effective polysaccharides) slow down and delay formation of both CSH and portlandite respectively.

Polysaccharide retarding effect depends on cement composition. C 3 A content might be considered as a key parameter. Extent of retardation is higher for cement with low C 3 A content which is in accordance with authors [START_REF] Ramachandran | Conduction calorimetric investigation of the effects of retarders on the hydration of Portland cement[END_REF][START_REF] Singh | Effect of gluconates on the hydration of cement[END_REF][START_REF] Monosi | Combined effect of lignosulfonate and carbonate on pure Portland clinker compounds hydration. III. Hydration of tricalcium silicate alone and in the presence of tricalcium aluminate[END_REF]. The retarding effect depends also on the chemical structure of admixtures. With respect to starch compounds, YD with a higher soluble fraction than NS and WD is more effective to retard cement hydration. A study with polysaccharides which have only one variable parameter, like degree of polymerization or degree of substitution, would produce valuable results in terms of admixture selection and properties prediction. Experiments on pure cement phase would provide useful information. It could be a first step to model pure phase hydration in presence of polysaccharide.

( 8 . 5 ±

 85 0.5)10 3 a determined by Gel Permeation Chromatography b determined by Pyrolysis-Gas Chromatography-Mass Spectrometry S : soluble fraction value (producers data)

FigFigFig 3 FTIR

 3 Fig 1b

Fig 1 Fig 2

 12 Fig 1 DTG curves of a C 1 , b C 1 +SE, c C 1 +YD.

Fig 4

 4 Fig 4 DTG curves of a C 2 , b C 2 +YD.

Fig 5

 5 Fig 5 Effect of cement composition on the retarding ability of YD.

Fig 6

 6 Fig 6 Effect of starch dextrinization on the evolution of CH formation in C 1 .

Table 1

 1

	Chemical composition	C 1	C 2
	SiO 2	22.2	21.9
	TiO 2	0.2	0.4
	Al 2 O 3	4.4	3.8
	Fe 2 O 3	0.3	4.3
	CaO	67.2	62.8
	MnO	0.01	0.04
	MgO	0.5	2.2
	Na 2 O	0.4	0.4
	K 2 O	0.05	0.39
	P 2 O 5	0.11	0.26
	SO 3	3	2.6
	Bogue approximation	C 1	C 2
	C 3 S	66	51
	C 2 S	14	24
	C 3 A	11	2.6
	C 4 AF	1	13.2

Table 1

 1 Chemical composition of cement and potential phases as determined by Bogue

	approximation.
	A.Peschard, A.Govin * , P.Grosseau, B.Guilhot, R.Guyonnet

Table 3 Formula

 3

		t 1 (hours)	t 2 (hours)	t 3 (hours)	Q 24h (J/g)
	C 1	2.4 ± 0.3	7.2 ± 0.3	11.2 ± 0.6	99 ± 5
	C 1 +CE	3.5 ± 0.5	9.5 ± 0.7	13.5 ± 0.7	97 ± 5
	C 1 +SE	2.2 ± 1.5	8.2 ± 2.1	14.5 ± 4.1	60 ± 8
	C 1 +NS	2.9 ± 0.2	7.7 ± 0.1	11.4 ± 0.1	94 ± 7
	C 1 +YD	7.9 ± 0.4	18.2 ± 1.1	-	9 ± 1

Table 4 Formula

 4 Fig 4b DTG curves of C 2 +YD A.Peschard, A.Govin * , P.Grosseau, B.Guilhot, R.Guyonnet Captions

		t 1 (hours)	t 2 (hours)	t 3 (hours)	Q 24h (J/g)
	C 1	2.4 ± 0.3	7.2 ± 0.3	11.2 ± 0.6	99 ± 5
	C 1 + 0.5% SE	2.2 ± 1.5	8.2 ± 2.1	14.5 ± 4.1	60 ± 8
	C 1 + 1.2% SE	3 ± 0.2	7.8± 0.4	-	45 ± 5
	C 1 + 1.5% SE	3 ± 0.2	7.9 ±0.3	-	28 ± 5
	C 1 + 2% SE	3.8 ± 0.3	7.9 ± 0.4	-	8 ± 1
	C 1 + 0.1% YD	2.2 ± 0.2	8.0 ± 0.3	12.2 ± 0.5	78 ± 5
	C 1 + 0.3% YD	4.6 ± 0.3	15.0 ± 0.5	-	48 ± 5
	C 1 + 0.4% YD	7.2 ± 0.3	16.4 ± 0.8	-	37 ± 4
	C 1 + 0.5% YD	7.9 ± 0.4	18.2 ± 1.1	-	9 ± 1

Table 1

 1 Chemical composition of cement and potential phases by Bogue approximation.

Table 2

 2 Properties of polysaccharides.

Table 3

 3 Times and heat released in 24 hours as a function of polysaccharide added.

Table 4

 4 Times and heat released in 24 hours for increasing concentration of SE and YD.

Acknowledgements

The authors acknowledge the financial support of Weber & Broutin company (Servas, France) and Rhône-Alpes Region. The authors are grateful to J.F. Halet for his assistance.


Pinoli Rahmani Saïd
email: said.rahmani@emse.fr

Debayle Johan

Geometrical Stochastic Modeling and Characterization of 2-D Crystal Population

Keywords: Boolean models, stochastic geometry, mean geometric covariogram, granulometry, morphology

Introduction

The presented work is motivated by a pharmaceutical issue. Indeed, the production of drugs involves a crystallization process in solution depending on several parameters and their quality is closely linked to its geometrical characteristics (size, shape, quantity ...). Therefore, it is necessary to control the distribution of the "crystal geometry" during the process. To meet this goal, an in-situ optical acquisition system provides images of the crystal population crystals' population, giving access to the projected particles during the process (see Figure 1). These images highlight the overlapping of the crystals (caused by the 2-D projection) and it is consequently difficult to individualize the particles for further characterization. Therefore, after a binarization process [START_REF] Presles | Size and shape estimation of 3-D convex objects from their 2-D projections: application to crystallization processes[END_REF] (see Figure 1), the objective of this work is to investigate stochastic geometrical models so as to represent these binary images and to get the geometrical characteristics of the crystals.

Modeling and objectives

Several materials can be modeled by random sets. In fact, the heterogeneity of the materials can be apprehended by a probabilistic approach [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF], [START_REF] Jeulin | Random texture models for material structures[END_REF]. Especially granular or fibrous media [START_REF] Galerne | Modèles d'image aléatoires et synthèse de texture[END_REF], [START_REF] Peyrega | Prediction des propriétés acoustiques de matériaux fibreux hétérogènes à partir de leur microstructure 3D[END_REF] can be represented by unions of overlapping particles (the grains) centered on random positions (the germs), thus giving rise to the germ-grain model.

Ξ = ⋃ 𝑥 𝑖 ∈Φ 𝑥 𝑖 + Ξ 𝑖 ( 1
)
Where Φ is a point process which generates the germs x i , and where the grains 𝚵 𝒊 are convex random sets independent and identically distributed.

Notice that this definition assumes the independence between the particles Ξ 𝑖 and their positions x i there is a more general definition authorizing the correlation between germs and grains [START_REF] Chiu | Stochastic geometry and its applications[END_REF].The proposed approach consists in representing the population of crystals by such a germ-grain model (Figure 2); that is to say that the point process Φ represent the spatial distribution of the particles centers and the convex random sets Ξ 𝑖 the particles themselves. The goal is to adjust the model to the real data by matching measurements computed both on real and simulated images. The real data are here obtained by an imaging acquisition system, meaning that we have only realizations of Ξ ∩ 𝑊 where 𝑊 is a bounded window, and we want to estimate the characteristics of Φ and Ξ 0 . To meet this objective, we focus on two points: firstly, we estimate the characteristics of Ξ from a realization of Ξ ∩ 𝑊 [START_REF] Molchanov | Statistics of the Boolean model for practitioners and mathematicians[END_REF], secondly we establish relationships between the characteristics of Ξ and the local characteristics of the model (Φ and Ξ 0 ). We also introduce an additional assumption: we assume the process Φ comes from a known type (Poisson process, Cox process ...). We will use the homogeneous Boolean model, a germ-grain model in which Φ is a homogeneous Poisson point process. This model is widely used because we have an analytical formula for the Choquet capacity.

𝑇 Ξ (𝑋) = 1 -exp(-𝜆 𝔼[𝐴(Ξ 0 ⊕ 𝐾)]) (2)
Several methods have been presented in the literature so as to connect the global characteristics of the Boolean model to the characteristics of the primary grain. In the plane and the space, the Miles's formulae [START_REF] Miles | Estimating aggregate and overall characteristics from thick sections by transmission microscopy[END_REF] or the minimum contrast method [START_REF] Stoyan | Stochastic Geometry and its application[END_REF] provide theoretical relationships to estimate the average values of the Minkowski functionals of the primary grain from global characteristics. Generally, the primary grain is assumed to have a known and deterministic shape, that is to say, the realizations of the primary grain are homothetic. So as to estimate the variations of the scaling factor from the expectation of the Minkowski functionals. For example, for a disc in the plane, the moments of the first and second orders of the radius of the primary grain are respectively proportional to its average perimeter and its average area. However, if we consider that the shape of the grain can vary, several issues remain unresolved: firstly the Minkowski functionals of a random convex set are not enough to characterize its shape and also their average will not provide sufficient information to characterize its variations. For instance, for a Boolean model whose grains has a shape that depends on several parameters (rectangle, ellipse ...); the estimation of geometric variations of the grain is not straightforward.

The proposed work provides estimators of the geometrical variations of the primary grain of the Boolean model without any assumption concerning its shape.

Method and results

It's well known that the second order moment of the area of the primary grain can be expressed as the integral of the mean geometric covariogram [START_REF] Galerne | Computation of the perimeter of measurable sets via their covariogram. Applications to random sets[END_REF]:

𝔼[𝐴(Ξ 0 ) 2 ] = ∫ ℝ 2 𝛾 Ξ 0 (𝑢)𝑑𝑢 (3)
Where 𝛾 Ξ 0 (𝑢) = 𝔼[𝐴(Ξ 0 ∩ Ξ 0 + 𝑢)] is the mean geometric covariogram of Ξ 0 . Thus, under the hypothesis 𝔼[𝐴(Ξ 0 )] < ∞ it is possible to obtain an estimator for 𝔼[𝐴(Ξ 0 ) 2 ] from an estimator of 𝛾 Ξ 0 (𝑢). In the special case of the Boolean model the mean geometric covariogram can be obtain by the covariance 𝐶 Ξ (𝑢), also called 2-points probability function by the following relation associate to Miles's estimator of intensity.

𝛾 Ξ 0 (𝑢) = 1 𝜆 ln (1 + 𝐶 Ξ (𝑢) -𝑝 Ξ 2 (1 -𝑝 Ξ ) 2 ) (4)
Where 𝑝 Ξ = 𝐶 Ξ (0) is the fraction area of Ξ. Let's notice that the Boolean is stable by convex dilatation; that is to say, for each compact convex set 𝐾 the random set Ξ ⊕ 𝐾 is also a Boolean model of same intensity 𝜆 and of primary grain Ξ 0 ⊕ 𝐾 . Consequently for any compact convex set 𝐾 the quantities 𝜆, 𝔼[𝐴(Ξ 0 ⊕ 𝐾)] and 𝔼[𝐴(Ξ 0 ⊕ 𝐾) 2 ] are estimable. Especially, considering a family of homothetic convex compact sets (𝑟𝐾) 𝑟>0 the second order moment of the area of dilated grain 𝔼[𝐴(Ξ 0 ⊕ 𝑟𝐾) 2 ] is a polynomial function in 𝑟, and can be expressed by Steiner's formula as follows:

𝔼[𝐴(Ξ 0 ⊕ 𝑟𝐾) 2 ] = 𝔼[𝐴(Ξ 0 ) 2 ] + 4𝑟𝔼[𝐴(Ξ 0 )𝑊 0,𝐾 ] + 𝑟 2 (4𝔼 [𝑊 0,𝐾 2 ] + 2𝐴(𝐾)𝔼[𝐴(Ξ 0 )]) + 4𝑟 3 𝐴(𝐾)𝔼[𝑊 0,𝐾 ] + 𝑟 4 𝐴(𝐾) 2 (5)
where 𝑊 0,𝐾 is the mixed area between Ξ 0 and 𝐾 [START_REF] Michielsen | Integral geometry morphological image analysis[END_REF], [START_REF] Schneider | Convex bodies: the Brunn-Minkowski theory[END_REF]. The quantities 𝔼[𝐴(Ξ 0 )] and 𝔼[𝑊 0,𝐾 ] can be estimated by the minimum contrast method [START_REF] Heinrich | Asymptotic properties of minimum contrast estimators for parameters of Boolean models[END_REF] or by the Miles's formulae [START_REF] Chiu | Stochastic geometry and its applications[END_REF]. Consequently the quantities 𝔼[𝐴(Ξ 0 ) 2 ], 𝔼[𝐴(Ξ 0 )𝑊 0,𝐾 ] and 𝔼[𝑊 0,𝐾 2 ] can be estimated by a polynomial approximation of 𝔼[𝐴(Ξ 0 ⊕ 𝑟𝐾) 2 ].

In practice the choice of the values of 𝑟 play an important role, and a precaution for the edge effect must be taken into account; if Ξ(𝜔) ∩ 𝑊 is a realization of Ξ in a bounded windows 𝑊, the dilated model Ξ ⊕ 𝑟𝐾 is only known on the eroded windows 𝑊 𝑟 = 𝑊 ⊖ 𝑟𝐾 and its realization in the windows 𝑊 𝑟 is (Ξ(𝜔) ∩ 𝑊) ⊕ 𝑟𝐾) ∩ 𝑊 𝑟 .

The choice of the convex compact set 𝐾 provides different characteristics of the primary grain; for a ball, the mixed area 𝑊 0,𝐾 is proportional to the perimeter 𝑈(Ξ 0 ), and for a segment oriented by 𝜃 ∈ [0,2𝜋] , it is proportional to the Feret's diameter of Ξ 0 in the direction 𝜃.

This method can be generalized to obtain 𝑛 𝑡ℎ order moments of (𝐴(Ξ 0 ), 𝑊 0,𝐾 ), using 𝑛-points probability functions (article in preparation).We have evaluated the performance of the proposed method for a disk and a segment by simulating several realizations of a test model in a squared window 500 × 500. The test model is isotropic model with a point process intensity 𝜆 = 100 500 × 500 ⁄ and with rectangular grains whose side lengths are independent and follow the Gaussian distributions 𝑁(40,10) and 𝑁(30,10). The results are presented in Figure 3, where the relative errors on the geometrical variations of the grains are shown as a function of the number of realizations.

Conclusions and prospects

A method has been presented so as to geometrically model and characterize a population of crystals. It provides estimators of the variations of the primary grain's morphological characteristics. Especially, using dilatation by a disk or a segment, the proposed method can be used to better characterize the geometry of a primary grain whose shape depends on two parameters (rectangle, ellipse...). We emphasize that our method can be used for any germ-grain model in which we can estimate the mean covariogram. In the future we are looking for more complex germ-grain models than the Boolean model. We are also interested in the influence of the model parameters and the observation window on the accuracy of the estimations. The prospect of describing a convex random set by the characteristics of the random process associated to its Feret's diameters seems to be also promising.

Figure 1 .

 1 Figure 1. Acquisition of an image of a crystals' population and its segmentation.

Figure 2 .

 2 Figure 2. Acquisition of an image of a crystals' population, and realization of a geometric stochastic model representing it.

Figure 3 .

 3 Figure 3. A realization of the test model (a) and the relatives errors of the second order moments of (area, perimeter) on (b) and of (area, Feret's diameter) on (c).


Ngoc Lam

Nicolas Marti

Antoine Blanche

Nicolas Ranc

Véronique Favier
email: veronique.favier@ensam.eu

André Chrysochoos

Nicolas Saintier

Fabienne Gregori

Brigitte Bacroix

Guillaume Thoquenne

Fabienne Grégori

Very High Cycle Fatigue for single phase ductile materials: slip band appearance criterion

Keywords: Very high cycle fatigue -VHCF, slip band, anisotropic elastic crystallinne, finite element

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The VHCF regime is associated with stress amplitudes lower than the conventional fatigue limit and as a result, numbers of cycles higher than 10 9 . The present project is a fundamental research project. However, it is motivated by industrial problems. Some mechanical components, such as pistons, rotating axes, in the transportation and more generally engineering industry, have been designed previously using fatigue resistance data at lower numbers of cycles (<10 7 cycles ; the regime of High Cycle Fatigue, HCF) whereas they must endure cyclic loading for a number of cycles higher than >10 9 cycles [START_REF] Bathias | Gigacycle fatigue in mechanical practice[END_REF]. Many experiments in the literature prove that failure occurs at stress amplitudes lower than the conventional fatigue limit. There is also growing demand for the development of robust life prediction and damage assessment methodologies that allow the safe extension of service lifetimes for transportation and power generation systems beyond their original design lifetimes; the methodologies must also allow the design of new efficient systems intended to operate at even longer lifetimes.

In the High Cycle Fatigue (HCF), much work, including analyzing and modeling proposed by Pr. H. Mughrabi [START_REF] Mughrabi | Cyclic stress-strain response and high-cycle fatigue behaviour of copper polycrystals[END_REF][START_REF] Mughrabi | Specific features and mechanisms of fatigue in the ultrahigh-cycle regime[END_REF] has been performed, and a fair understanding has emerged over the years. Cyclic strain localization in Persistent Slip Bands (PSBs) is generally observed and accepted as the critical onset of fatigue damage in the HCF regime in ductile single-phase materials and also in many precipitation-hardened alloys. Once formed, PSBs, accompanied by the development of extrusions and notch-like deepenings, commonly referred to intrusions, lead to cracks. For ductile single-phase materials, referred as type I materials, only few investigations (mainly on copper [START_REF] Stanzl-Tschegg | PSB threshold and fatigue limit of polycrystalline copper in the VHCF-regime[END_REF][START_REF] Anja | Fatigue damage in copper polycrystals subjected to ultrahigh-cycle fatigue below the PSB threshold[END_REF][START_REF] Stanzl-Tschegg | Mechanisms of strain localization, crack initiation and fracture of polycrystalline copper in the VHCF regime[END_REF][START_REF] Cretegny | AFM characterization of the evolution of surface deformation during fatigue in polycristalline copper[END_REF][START_REF] Huang | The study of fatigue in polycrystalline copper under various strain amplitude at stage I: crack initiation and propagation[END_REF][START_REF] Lukas | High Cycles Fatigue life of metals[END_REF]) are available regarding the shape of the fatigue life curve and the damage evolution. From these experiments, it is suggested that crack initiation takes place at the surface owing to the accumulation of very small but irreversible plastic deformation over very large number of cycles [START_REF] Mughrabi | On multi-stage fatigue life diagrams and the relevant life controlling mechanisms in ultrahigh-cycle fatigue[END_REF]. However, there is no possibility for the formation of a pronounced PSB structure at very low stress amplitudes below the "PSB threshold" (associated with the HCF regime). Therefore, it is necessary to conduct systematic investigations of type I materials after cyclic deformation in the VHCF regime in order to clarify the role of the above mentioned mechanisms and parameters.

In this work, we focused on pure copper which is face-centered cubic (fcc) single-phase metal. This material is a good candidate because its HCF behavior has been largely studied in literature. To perform experiments up to a very high number of cycles in a reasonable time, ultrasonic equipment at a testing frequency of 20 kHz was used. It takes only fourteen hours to perform 10 9 cycles at a testing frequency of 20 kHz, whereas it takes one hundred days at 100 Hz. The main objective is to characterize the earliest signs of plasticity, i.e. the appearance of the earliest slip markings (SMs) (associated with slip bands (SBs)) on the specimen surface. Their evolution with stress amplitude and with number of cycles was also investigated. The evolution of the intrinsic dissipation with the stress amplitudes was estimated. Quasi 3D finite element simulations taking into account the polycrystalline nature of the material were performed to investigate the role of the elastic anisotropy in slip band appearance as suggested by [START_REF] Pommier | Arching' effect in elastic polycrystals: implications for the variability of fatigue lives[END_REF][START_REF] Sauzay | Cubic elasticity and stress distribution at the free surface of polycrystals[END_REF][START_REF] Winter | Dislocation microstructures in fatigued copper polycrystals[END_REF]. A special attention was paid to the role of grain orientation and of the initial surface state following the work of Mareau and Favier et al [START_REF] Mareau | Micromechanical modeling of the interactions between the microstructure and the dissipative deformation mechanisms in steels under cyclic loading[END_REF][START_REF] Mareau | Influence of the free surface and the mean stress on the dissipation in steels under cyclic loading[END_REF] in HCF regime.

Material and experimental procedure

Material and specimens

Commercial polycrystalline copper CuOF 99.95% in initial form of a hot rolled plate of 14 mm of thickness was tested. The specimens were extracted from the centre of the plate and the microstructure of material is presented in Fig. 1. The mean grain size was about 30 μm with a lot of annealed twins. The material is quasi-isotropic textured. The present study was conducted by using an ultrasonic fatigue machine at 20 kHz with cylindrical and flat hourglass shaped specimens. The specimen dimensions were calculated so that the free resonant frequency of specimen in the first longitudinal mode is 20 kHz. All specimens were treated at 250°C for 60 min to relieve the residual stress (ASM specialty Hand book) without change of microstructure. Then, a procedure of mechanical and electrolytical polishing was executed till remove all hardened layers on surface specimen. As a result, before testing, the specimens were mirror polished and without residual stresses. The dynamic Young modulus in calculations was E d = 130 GPa [START_REF] Stanzl-Tschegg | PSB threshold and fatigue limit of polycrystalline copper in the VHCF-regime[END_REF]. The stress/strain amplitude along the specimen was estimated by an one-dimensional elastic calculation. It was in a good agreement with the experimental measurement obtained by the strain gages. The geometry and the stress amplitude distribution in specimens are presented in Fig. 2.

Experimental procedure

To perform experiments up to the VHCF regime, an ultrasonic equipment at a testing frequency of 20 kHz was used. The S-N curve of the material in the VHCF regime was first established. The fatigue tests were carriedout with cylindrical specimens at various stress amplitudes with a constant stress ratio (σ max /σ min ) R = -1.

In a previous work [START_REF] Favier | Microplasticity and energy dissipation in very high cycle fatigue, Fatigue Design[END_REF], the self-f f heating over ultrasonic fatigue tests was measured by an infrared camera. Figure 3 shows the spatial average of the temperature increase (heating) measured at the surface specimen and the deduced intrinsic dissipation. The temperature increase remains lower than 20°C for stress amplitudes lower than 60 MPa. However it reaches 70°C for higher stress amplitudes. This stronger temperature increase may change properties of copper. That is why in order to prevent from a heating effect at high stress amplitudes, the specimens used for establishing the S-N curve were cooled by a cold air gun during the test till rupture. To study the formation and development of slip markings on the specimen surface at very low stress amplitudes and after very high numbers of cycles, the flat specimens were used. This type of specimen is appropriate to study the microstructure around the SMs by Electron Back Scatter Diffraction (EBSD) technique and to measure the roughness of SMs by Atomic Force Microscope (AFM) and temperature and intrinsic dissipation by an infrared camera. The experiments were periodically interrupted after several specified numbers of cycles for surface observations. The observations were exclusively made at the centre of the cross section. Indeed, the stress amplitude at the border is lightly greater (~ 5 %) than one at the centre of the cross section.

Results

S-N Curve

The S-N curve of the studied pure copper obtained by using ultrasonic fatigue machine at 20 kHz is shown in Fig. 4. It displays two regimes. In the first regime, failure occured between 3.6x10 6 and 5x10 7 cycles with stress amplitude ranging from 115 MPa to 97 MPa. The first regime is associated with the high cycle fatigue (HCF) regime. In the second regime, failure took place at numbers of cycles greater than 5x10 7 cycles. The stress amplitude is lower than 97 MPa. The S-N curve slope in VHCF is more than 10 times smaller than in HCF (-0.4 MPA/Log(Number of cycles) versus -5.3 MPa/Log (Number of cycles)). In other words, the number of cycles very quickly increases with decreasing the stress amplitude. At Δσ/2 = 91.2 MPa, we did not observe any failure up to 5.4x10 9 cycles at which the test was stopped. Our results are in very good agreement with Stanzl-Tschegg et al. [START_REF] Stanzl-Tschegg | PSB threshold and fatigue limit of polycrystalline copper in the VHCF-regime[END_REF]'s experiment results obtained on similar pure polycrystalline copper. Their results revealed a "fatigue limit at 1x10 10 cycles" also called "fatigue life threshold" equal to at Δσ/2 = 92.2 MPa.

Three types of slip band in fatigue

The appearance and evolution of slip markings in the VHCF regime for a wide range of stress amplitude from 38 MPa to 115 MPa have been studied. It was found that he lower the stress amplitude, the higher the number of cycles necessary to observe the earliest SMs. Beside, we distinguished 3 types of SMs: type I, type II and type III (Fig. 5) which differ from their location within the microstructure and their morphology. All these SMs were shown to be Persistent Slip Markings (PSM), namely when the surface was electropolished after fatigue, the SMs reappeared at the same sites on the specimen surface when fatigue was resumed.

Type I

Fig. 6 shows SEM images of slip markings type I. They are parallel to each other. They are long, straight and crossing the grain. The width of each SM varies from 0.5 to 3 μm. The distance between the SMs depends on grain size and it is larger for larger grains. It is about 5~7 μm for a grain size of 30 μm f f (30 is the mean grain size of the studied copper). The surface is smooth between the SMs. A quantitative measurement with AFM shows that the height of the extrusions associated with type I SMs is about 0.6 μm and the width is 2 μm (Fig. 7). New parallel slip markings appear inside the grain and the height of SM rises with increasing number of cycles. The majority of SM type II was formed at twin boundaries. The average height is close to SM type I height (Fig. 9) but the maximum height can reach 1-2μm depending on the stress amplitude and the number of cycles. The initial SMs become longer, larger and higher with increasing the number of cycles. Figure 10 shows the evolution of a type II slip band at Δσ/2 = 55 MPa after 10 7 cycles, 10 8 cycles and 2x10 8 cycles. This SM was formed at a grain or twin boundary and it grew only in the left grain. No new SM appeared in the close neighborhood of the initial SM.

Type III

At the very earliest stages, slip markings of type III are invisible under optical microscopy (MO). SEM observations revealed the presence of clusters of very fine slip markings formed inside the grains and near the grain boundaries. Two slip marking directions were usually observed. The first one more marked and occupies the majority of the cluster surface. It is very probably associated with the slip band of the primary slip system. The second is less marked and is called "secondary" slip band (Fig. 11). Fig. 12 shows an AFM image of slip markings of type III at Δσ/2 = 50 MPa after 2x10 8 cycles. Their roughness is less than 50 nm.

Relation between stress amplitude and formation of different types of f f SM

The appearance of the different types of slip markings depends essentially on stress amplitude. At high stress amplitudes, i.e. higher than "fatigue limit at 1x10 10 cycles" at Δσ/2 = 92.2 MPa, the majority of the first slip markings is of type I. The amount of SM type I decreases with decreasing stress amplitudes and SM type II become more and more dominant. At very low stress amplitudes such as 44 MPa, type III is dominant.

The numbers of cycles necessary to form the earliest SMs and the type of SMs as a function of the stress amplitude are showed in Fig. 13. At a given stress amplitude, the test was regularly interrupted to observe the slip markings on the specimen surface. In Fig. 13, the experimental data indicate the last time observation when no slip markings appeared and the next time observation when the earliest slip markings are detected using SEM. Our results are in good agreement with the experimental results obtained by Stanzl-Tschegg and Schönbauer [START_REF] Stanzl-Tschegg | PSB threshold and fatigue limit of polycrystalline copper in the VHCF-regime[END_REF]. From 62 MPa to 34 MPa, the stress amplitude varies linearly with the logarithm of the number of cycles necessary to form the earliest slip markings. The S-N curve and the curve for the appearance of the earliest slip bands are compared in Fig. 14. Up to 10 8 cycles, at a specified number of cycles, the stress amplitude required to break the specimen is twice larger than the stress amplitude to form the earliest SMs on the specimen surface.

Intrinsic dissipation

In a previous paper [START_REF] Blanche | Dissipation measurements in dynamic high cycle fatiguein copper during VHCF loading[END_REF], the average intrinsic dissipation was estimated from thermal field measurements over fatigue tests. Fig. 15 exhibits the intrinsic dissipation (divided by in order to have volume energy rate (W/m 3 ) expressed in (°C/s)) at 10 6 cycles for various stress amplitudes ranging from 20 to 90 MPa. Results stressed that dissipated energy exists whatever the applied load. Higher the stress range, higher the intrinsic dissipation. We can note that even below about 30 MPa, at which no slip marking appear even at number of cycles higher than 10 10 and sometimes considered as an "irreversible threshold" in literature [START_REF] Stanzl-Tschegg | PSB threshold and fatigue limit of polycrystalline copper in the VHCF-regime[END_REF], a non-zero intrinsic dissipation can be detected (Fig. 15 (b)). It was also found that dissipation always increased with the number of cycles revealing continuous changes within the material. This result is consistent with the gradual changes of the specimen surface over cycles.

Calculation of the stress field at the surface of a polycristal l

Slip markings of type II and III are localized at or close to grain and twin boundaries. These observations reveal the key role of grain/twin boundaries acting as stress concentration factors. At such low stress amplitudes (with regard to the macroscopic yield stress and fatigue strength), the material is elastic (except at some locations) and the anisotropy of the grain elastic behaviour is responsible for stress heterogeneity and stress concentrations at grain and twin boundaries. To estimate the stress field within the material and to relate it to the sites of SMs, quasi 3D finite element simulations were performed. The real polycrystalline microstructure observed at the specimen surface using EBSD was geometrically represented and discretized (Fig. 16). Four zones of interest (called A, B, C, D) were chosen because they contained SMs of type II. Considering the four zones, the total number of grains is equal to 192. The total number of SM type II is equal to 11. The 2D representations were then extruded in the third direction. Quasi 3D simulations were performed using the ABAQUS software and considering the orientation and the cubic elasticity of each grain. The plastic behaviour is neglected since we are just interested in the earliest SM appearance for which the material is mainly elastic. A uniform compression/tensile uniaxial stress along E 2 axis (specimen axis) was applied on the studied zone boundaries normal to the loading direction. The others surface boundaries were free. The Von Mises stress field clearly exhibits stress concentrations at some grain boundaries (Fig. 16 (e)). Copper has a face-centered cubic structure and 12 {111} <110> slip systems. To correlate slip markings and slip systems, the resolved shear stress on each slip system s was calculated according to the following equation [START_REF] Bathias | Gigacycle fatigue in mechanical practice[END_REF] where σ is the local stress tensor and R is the Schmid tensor. For each grain, the maximum resolved shear stress was identified and the location at which it was found was noted. For grains displaying one SM type II, the maximum resolved shear stress was systematically found higher close to grain/twin boundaries than within the grain. In some cases, two very close values were found. To determine the active slip system, we used the marking observed on the specimen surface and measured the π angle. The π angle is the angle between the marking of the slip band and the loading axis (Fig. 17). This measured angle was compared with the calculated one assuming the active slip system. The slip system having one of the highest resolved shear stress and the closed angle π to the experimental one was said to be the active slip system according to the Schmid criterion. Fig. 18 displays the number of grains (over the 192 grains) having a maximum resolved shear stress divided by the applied stress max / equal to a given value. The values of max / range from 0 to about 1. They were divided into intervals of 0.02 except for the first interval which includes values between 0 and 0.35. As mentioned above, 11 grains over the 192 grains exhibited one slip marking. These grains are indicated by the following notation: G Number of the grain / Label of the zone (A, B, C or D). For example, G32/A means a slip marking in the numbered 32 grain belonging to the A zone. In addition, the maximum Schmid factor was calculated for each grain. Note that the maximum Schmid factor is also the ratio of the maximum resolved shear stress over the applied stress max / for which the maximum resolved shear stress is calculated assuming a uniform stress state ( = ), in other words isotropic elasticity. To sum up, the number of grains having max / of a given value and calculated assuming isotropic and anisotropic elasticity was plotted as a function of max / (Fig. 18). Note that, while the Schmid factor is bounded between 0 and 0.5, max / calculated assuming anisotropic elasticity reaches 1 because of the stress concentration or/and triaxiality. The Schmid factor ( max / ) for the grains having one slip marking ranges from 0.44 and 0.5. However, many other grains (115 grains) which do not have any slip marking also have a maximum Schmid factor between 0.44 and 0.5. This result demonstrates that considering the maximum Schmid factor is not a good criterion to predict the grains having one SM type II. In addition, it cannot predict the location of the marking since max / is uniform within a grain.

The max / values, calculated assuming anisotropic elasticity, for the grains having one slip marking ranges from 0.78 and 0.9. In that case, only 4 other grains have their max / value in this range and do not have any slip band. This result show that considering the maximum max / is a good criterion to predict the grains having one type II slip band: it predicts more than 70% of grains having one slip band type II. In addition, it predicts the location of the band.

The previous analysis also revealed that the 11 slip bands are located at or very close to a grain boundary whose disorientation is very high from 58.5 ° to 60 °. 9 of these bands appeared at or very close to a twin boundary. The associated slip plane is a plane parallel to the plane of the twin boundary (G1/A, G2A, G32/A, G20/B, G18/C, G28/C, G47/C, G36/G10 and G10/D). For The last two slip bands (G27/C and G8/D), the grain boundary is not a twin boundary, although the disorientation between the two grains is close to 60°. This result proves that twin boundaries are preferred location for appearance of slip markings (slip bands) in fatigue at very low stress amplitudes.

Discussion and conclusion

In LCF and HCF regime, the crack initiation in high purity polycrystalline metal mainly occurs along PSBs crossing the grains (type I). PSBs appear in well-oriented grains for plastic slip and are generally along one dominant direction given by the primary slip system. An analysis based on the maximum Schmid factor (and so which does not take into account the elastic anisotropy) predicts well the grains which have PSBs and the active primary slip system [START_REF] Blochwitz | Analysis of activated slip systems in fatigued nickel polycrystals using the EBSD-technique in the scanning electron microscope[END_REF][START_REF] Bartali | Identification and analysis of slip systems activated during lowcycle fatigue in a duplex stainless steel[END_REF]. In the VHCF regime, it was shown that slip marking of type II perfectly correspond to slip bands associated with a primary slip system. An analysis based on the maximum resolved shear stress (Schmid criterion) is relevant as soon as the stress heterogeneities due to elastic anisotropy are considered. Such a criterion predicts the location and the active primary system of the slip band type II. Similar investigations on slip markings type III are ongoing. Literature results demonstrate the role of twin boundaries for crack initiation [START_REF] Neumann | Crack initiation at grain boundaries in FCC materials -Strength of Metals and Alloys[END_REF][START_REF] Thompson | The origin of fatigue fracture in copper[END_REF]. Very few works were dedicated to the relationships between slip bands and twin boundaries. They rather focused on the role of disorientation [22]. In this paper, the key role of twin boundaries on the slip band type II was proved for pure copper specimens loaded in the VHCF regime. Slip markings were found to continuously change over the cycles. Consistently, the dissipation always increased with the number of cycles even at stress amplitudes at which no slip marking was visible at 10 10 cycles. These results raise the question of the existence of a fatigue limit for pure copper.

Fig. 1 .

 1 Fig.1. Microstructure of material

Figure 3 .

 3 Figure 3. (a) Average temperature and (b) Average intrinsic dissipation during fatigue test at different constant stress amplitude fatigue test for CuOF 99.95% [16].

Fig. 2 .

 2 Fig.2. Ultrasonic fatigue specimens (a) cylindrical specimen (b) flat specimen

Fig. 4 .

 4 Fig.4. S-N curve of pure copper obtained with ultrasonic fatigue tests

Fig. 5 .

 5 Fig.5. Three types of slip markings in VHCF (a) Type I ; (b) Type II ; (c) Type III and its zoom (d)

Fig. 6 .Fig. 7 .

 67 Fig.6. Slip markings type I appeared at Δσ/2 = 89 MPa -10 6 cycles

Fig. 8 .Fig. 9 .

 89 Fig.8. Slip marking type II at Δσ/2 = 57.5 MPa -10 7 cycles (a) SEM image and (b) EBSD analyze at the same zone shows the appearance of SM type II (arrows) at grain boundary

Fig. 10 .Fig. 11 .Fig. 12 .

 101112 Fig.10. Evolution of SM type II with the numbers of cycles at 55 MPa (a) after 10 7 cycles ; (b) after 10 8 cycles and (c) after 2x10 8 cycles

Fig. 15 -

 15 Fig. 15-(a) Dissipation vs. Stress amplitude Fig. 15-(b) Zoom between 20 and 60 MPa

Fig. 16 .

 16 Fig.16. (a) Coordinate system ; (b) SEM image and (c) EBSD image of some type II SB (d) Simulation on Abaqus ; (e) Von-Mises stress distribution

Fig. 17 .

 17 Fig.17. Definition schema of angle π .Fig.18displays the number of grains (over the 192 grains) having a maximum resolved shear stress divided by the applied stress max / equal to a given value. The values of max / range from 0 to about 1. They were divided into intervals of 0.02 except for the first interval which includes values between 0 and 0.35. As mentioned above, 11 grains over the 192 grains exhibited one slip marking. These grains are indicated by the following notation: G Number of the grain / Label of the zone (A, B, C or D). For example, G32/A means a slip marking in the numbered 32 grain belonging to the A zone. In addition, the maximum Schmid factor was calculated for each grain. Note that the maximum Schmid factor is also the ratio of the maximum resolved shear stress over the applied stress max / for which the maximum resolved shear stress is calculated assuming a uniform stress state ( = ), in other words isotropic elasticity. To sum up, the number of grains having max / of a given value and calculated assuming isotropic and anisotropic elasticity was plotted as a function of max / (Fig.18). Note that, while the Schmid factor is bounded between 0 and 0.5, max / calculated assuming anisotropic elasticity reaches 1 because of the stress concentration or/and triaxiality.

.

Fig. 18 .

 18 Fig.18. max/ calculated assuming isotropic or anisotropic elasticity for the 192 grains contained in zones A, B, C and D. 11 grains contained one slip marking type II.

Acknowledgements

The authors are grateful for financial support from Agence Nationale de la Recherche France ANR-09-BLAN-0025-01 and to company Griset for supplying copper. The authors thank also very much Professor Haël Mughrabi for the scientific discussions.


Lucile Pentecoste

Pascal Brault
email: pascal.brault@univ-orleans.fr

Anne-Lise Thomann
email: anne-lise.thomann@univ-orleans.fr

Pierre Desgardin

Thomas Lecas

Tayeb Belhabib

Marie-France Barthe

Thierry Sauvage

A.-L Thomann

M.-F Barthe

Low Energy and low fluence helium implantations in tungsten: Molecular dynamics simulations and experiments

Keywords: Molecular Dynamics, Helium retention, plasma-wall interactions

Introduction

Due to its low sputtering yield and high melting point, tungsten is the material chosen to compose the first wall of the divertor of ITER fusion reactor. This part of the chamber will face the most severe conditions of the plasma, such as high fluxes of light particles with a wide range of energies and high heat loads [1][2][3] . Studies on plasma material interactions are currently carried out to anticipate the evolution of this metal while exposed to fusion plasma. Several reactors are developed to investigate the interactions of such fusion plasma on W such as Nagdis-I and II in Nagoya or Pisces in San Diego 4,5 .

Gas implantation, diffusion and retention into metals are well-explored phenomena 6,7 Experimental studies have evidenced that He ion fluxes on W substrates lead to various surface morphology changes such as blistering, fuzz formation and even flaking [8][9][10][11] . Those features appearance is strongly dependent on the flux and incident energy of the ions; it also depends on the substrate temperature during the implantation process. At high energies, keV to MeV, mean penetration depths are high (100 nm to microns) and defects induced by elastic and inelastic collisions lead to the formation of damages in the substrates. However at low ion energies (below W displacement thresholds: 90 eV transferred from the impinging ion) the formation of the characteristic damage features is still under investigation.

Experiments performed in Nagdis-II by Nishijima et al 8 showed damages as holes appeared for high fluxes (10 19 cm -2 s -1 ) and high doses or number of impinging atoms (10 23 cm -²) for energies ranging from 5 to 30 eV at high temperature (1300 to 2950 K). It was also reported by Lee et al 12 that a limited number of He could remain in the material for energies below W displacement threshold for implantations performed with fluxes ranging from 10 14 to 10 15 cm -2 .s -1 and doses from 10 17 to 10 19 cm -2 . Simulations using molecular dynamics were also performed to study the phenomena involved in the implantation process at low energies. Henriksson et al. studied the behaviour of helium with doses below 10 14 cm -2 and with energies varying from 50 to 200 eV at 0 and 300 K 13 . Li et al. investigated implantations of 10 15 cm -2 impinging helium ions on W at 200 eV at temperatures from 300 to 1500 K 14,15 .

Both observed formation of growing and exploding clusters which liberate their He content to the surface. In the case of slightly higher energies (300 eV to 1 keV), the mean penetration of the He is observed to be deeper and no "cluster rupture" was identified due to simulations stopped too early in the course of He interactions 15 .

In the present article, the objective is to describe the phenomena observed during the simulation of implantation process carried out at doses ranging from 10 14 to 10 17 cm -2 and with a kinetic energy of 300 eV at a temperature of 300 K. the simulations were performed by molecular dynamics and SRIM 16 calculations. MD simulations initial conditions are chosen to closely fit with experiments. The aim is to explain the saturation of the retained He atom number which does not occur because of the "cluster ruptures" but is however evidenced for implantation experiments performed at energies below W displacement threshold.

Then after a review of the details of the simulation and experimental methods, the MD simulation and experimental results are discussed and compared.

Model and Experiments

Model, Calculations: starting parameters and data treatment

To investigate the phenomena occurring during the implantation of Helium into the tungsten, molecular dynamic simulations are carried out on an Alineos quadriprocessor 8 cores high performance computer using LAMMPS [START_REF] Plimpton | Lammps Sandia[END_REF] . Besides, the SRIM code, which is widely referred to in ion/material interaction studies, is used in addition to compare the statistics of the resulting implanted atom concentration and depth profiles.

1 1 1                    k s T T dt   .
Where dt is the integration timestep (here 0.1 fs), Ts is the targeted temperature, Tk is the current temperature and  is the time needed for reaching the targeted temperature Ts. This thermostat allows proper energy relaxation of the substrate for mimicking energy dissipation in the experiment on a reasonable time scale: it is set here to 1 ps, a value consistent with a model of electron-phonon coupling 19,20 that allows estimation of the relaxation time  by the formula:

Z k Lne T m B e D F e 2 2    
ΘD is the surface Debye temperature, Ts the substrate temperature, L the Lorentz number, n the electron density, e the electron charge, kB the Boltzmann constant, Z the valence, me the electron mass, κ the thermal conductivity, εF the Fermi energy and Tk is the kinetic temperature. For W,  = 1 ps and is used as the relaxation time of the Berendsen thermostat.

Substrate thickness is chosen in order to limit to 0.5% the implanted He atoms crossing the substrate through the fixed W atoms layer. This allows observing all the implanted He atoms, their localization and surrounding, considering the channelled atoms as well. Also, the high number of thermostated W atoms combined with the energy dissipation handled by the thermostat helps keeping the desired temperature of the substrate at 300 K such as in the experiments. Moreover the thickness is chosen in order to keep unperturbed thermostated layers at the substrate bottom and so observe a smooth transition toward the fixed bottom layers. Thus none of the perturbations of the W lattice are missing.

Helium ions are treated as atoms, as they are neutralized at the W surface, so classical MD simulations are then applicable. He atoms are injected one by one and in a direction normal to the surface from a randomly chosen position above the box. This allows avoiding any interaction before reaching the substrate's surface. The frequency of the helium atom injection is set between 2.5 and 5 at.ps -1 . Such impingement frequency allows no lateral interaction between induced perturbations in the W crystal during the dissipation time, due to the large enough substrate width. The kinetic energy of the impinging He is set to be 300 eV. The lattice displacement threshold for a W atom is of 90 eV, and the maximum energy transferred from a He atom to a W atom is 25 eV when the He collides with an energy of 300 eV. This value is obtained from the elastic collision energy transfer formula: Etr = Ei*4m1m2/(m1+m2)². The number of helium atoms impinging the W surface ranges from 2000 to 10000 for a surface of 53.8x53.8 Å 2 . It corresponds to the experiment doses ranging from 10 14 to 10 17 cm -2 .

Handling molecular dynamics only requires the knowledge of inter-atomic potentials and a set of initial conditions. We use the Juslin and Wirth potentials 21 : W-W interactions are described by a modified Ackland-Therford potential and W-He interactions are modelled from DFT study. Finally, He-He interactions in the substrate are described by Beck's potential 22 , which would allow He cluster growth.

Analysis of He retention rate, of snapshots during implantation simulations and of implanted He atom depth profiles will give insights into the interaction dynamics of the He with W. SRIM calculations are also carried out for determining retention rates and depth profiles using SRIM 2013 with 2008 stopping power parameters 16,23 . SRIM software uses Binary Collision Approximation and is commonly used in simulating implantation processes. It considers an amorphous substrate with a density matching the W structure (19.3 g cm -3 ) and ignores either the interaction between He atoms or effect of the accumulation in the substrate. Moreover SRIM calculations do not account for recombination of interstitials with the vacancies, nor for defect clustering and irradiation-induced amorphization. Thus comparison with MD simulations, where these effects are considered, would allow to highlight their role on the retention of He in W.

Experiments

Polycrystalline tungsten samples are 7x7x0.3 mm pieces of 99.95% purity. The samples were mechanically polished using diamond pastes with grain size ranging from 6 down to 0.1 µm. They were then annealed under high vacuum (10 -7 to 10 -6 mbar) and at high temperature (1600 °C) for one hour. This annealing removes the polishing induced defects and reveals the polycrystalline structure of the substrates in the near surface region as observed in the figure 1. The most important feature is the large grain size minimizing grain boundaries size obtained after annealing. from 10 11 to 10 14 cm -².s -1 . Implantations are performed at 300 K for various durations to obtain a range of ion impinging number from 10 14 cm -2 to 10 17 cm -2 . The 3 He isotope was chosen instead of 4 He because it can be quantified using the Nuclear Reaction Analysis technique (NRA).

The total amount of 3 He retained in the samples is measured by NRA in a dedicated apparatus DIADDHEM at the CEMHTI 24 laboratory, Orléans.

The nuclear reaction used to quantify the He content is 3 He(d, 1 H)α [25][26][27][28] . A 900 keV deuterium beam is used to create nuclear reactions with 3 He atoms present in the W substrate. The nuclear reaction products are an 4 He ( particle) and a 1 H (proton) with energies of 2-3 and 15 MeV, respectively. The concentration of retained 3 He is determined by counting the emitted protons for a given deuterium beam charge. The total proton number is normalized by the beam charge and compared to the one measured with the same conditions in a standard sample containing 10 16 3 He cm -2 . The 3 He concentration is then determined using a simple proportionality relation.

Even if MD simulations are concerned by 4 He, the mass difference between the two isotopes 3 He/ 4 He is not expected to induce a variation in the defects formed into the tungsten lattice.

The maximum energy transferred by a 3 He atom (19 eV) remains still below the displacement threshold of a lattice W atom. J. Amano et al 29 determined experimentally using FIM technique the diffusion coefficients for both isotopes and concluded there was no strong effect of the mass on the enthalpy of migration and only a difference less than 15% in the diffusion coefficient. Thus, we consider that both isotopes have the same effects on W substrate when implanted and that the comparison between MD and NRA results is possible.

Results and Discussions

Retention rate

MD simulation shows a decrease of the He retention rate when increasing the He atom number. This is different from SRIM calculations which lead to a constant retention rate (Table1). Beyond the nature of interactions considered in each code, an important difference between MD and SRIM calculations lies in the fact that each SRIM impact is on unperturbed amorphous material without any dose effect. Conversely, MD simulations consider helium accumulation under the substrate surface. Thus, it appears to be an important parameter impacting the retention rate when increasing the incident He atoms number due to He-He and He-W interactions leading to a modified W lattice. Henriksson et al. 13 reported such a comparison between MD and SRIM. They observed such a large difference, for low energy implantations, which is decreasing when impinging atom energy reached 1 keV. The SRIM code is known as a powerful tool in the case of high energy implantations as it is based on binary collision approximations. However those interactions are less important in the case of low energy implantations where many body interactions come into play with moving atoms.

The difference found between MD and SRIM highlights the limitation of SRIM code for the low kinetic energies and shows the importance of considering relevant many-body interactions.

Implantation process

Snapshots of the MD simulations are presented in figure 2a At the beginning of the implantation process, helium atoms slow down at various depths and diffuse in the lattice of the substrate (Figure 2b). Early after the beginning, small He cluster start to form and they do not diffuse anymore when they have a size greater than 6 He atoms, in agreement with previous studies. 30,31 They behave as traps for new incoming He atoms (figure 2 a-d Big clusters form between 30 and 60 Å depth below the outermost surface and then coalesce to form broad bubbles (figure 2 f to j). Helium concentration in this 30-60 Å slice is continuously increasing. In the meantime, the W lattice is continuously expanding along the z axis. This is mainly due to "loop punching event": [13][14][15] the W atoms are pushed away as a

He bubble is growing. The W atoms are thus moving along (111) planes as "crowdion cascade" towards the surface where they are stacking in crystal positions above the initial surface. This results in increasing the height of the current surface 15 as illustrated in figure 5. This phenomenon as already been observed experimentally by Gilliam et al. 11 From this step, the simulation is stopped as the implantation begins on a "refreshed" substrate surface already containing implanted He atoms.

For investigating the processes leading to the rupture of the substrate, we study the evolution of the depth concentration profile. The shape of the depth profile obtained with MD simulation is displaying a stratification of the He concentration. Such stratified depth profile was observed by Li et al. 15 for simulations of 200 eV helium implantation in (001) tungsten at 300 K. They observed a near-surface zone below three or four W layers such as the one we defined as A (fig 7 -b) where single helium atoms are located and small clusters formed. They revealed that such small clusters led to changes of the morphology near the substrate surface. However, the biggest clusters formed in their case were located 45 Å below the surface. These clusters can also undergo a "rupture" which consists in emptying from most of their He atoms content, ejecting them away from the surface. In our case, the highest concentration of helium is located deeper (60 -70 Å) than for 200 eV impinging atoms, so clusters are unable to explode and eject all or part of their helium content. The clusters accumulate, coalesce and eventually separate the substrate in two parts.

Depth profile evolution

Evolution of the helium depth distribution is plotted in growth. Regions A and C surrounding the bubble the B zone) do not seem to evolve much after the maximum of concentration is reached (frame g). In the case of the C region, it is due to the impossibility for the helium atoms to cross the bubble where they are trapped.

On the other hand, a slight increase of the concentration in the A zone is observed and corresponds to balance between incoming He atoms which are rapidly slowing down thus contributing to the formation and growth of small clusters, and those which enter temporarily and escape from the substrate. These latter are recorded as backscattered ions.

In order to check if the periodic boundaries used in the simulations are not responsible for such behaviour (cluster formation and saturation) in response to the implantation process, the simulation was performed on a two times larger substrate with the same proportional number of impinging helium atoms. The cluster formation and substrate rupture proceeded similarly.

The highest concentration of helium reachable in the W substrate corresponds to the time when the substrate is starting to break. We define this as the saturation level. After the flaking phenomenon, the surface is considered as partially regenerated leading to a new surface already containing some helium atoms that is exposed to the helium flux. Due to calculation time limitations when the saturation occurs, the simulation is stopped and the implantation process is no more considered on the flaked substrate surface.

The saturation level is reported in Table 2 which lists the number of retained helium atoms and the corresponding implanted dose as a function of the impinging He atom number. The runs 3 and 4, give similar saturation levels since the saturation phenomenon appears between 4000 and 5000 impinging helium atoms, corresponding to doses of 1.38 and 1.73 10 16 He cm -2 . The surface is twice larger for run 5 so the dose is the same as for the run 3, the number of retained atoms is about twice the number for the smaller surface substrate. It means the phenomena are not surface size dependent. In order to check the hypothesis of saturation phenomenon, experiments are carried out performing helium implantations in conditions matching the simulations parameters: 300 eV Helium ion implantation at room temperature and number of impinging He atoms below and above the predicted saturation level.

Helium retention

The helium retention corresponds to the number of atoms which remains in the substrates after the implantation experiment. Experimentally, this number is measured by NRA. In the case of the simulations, Molecular Dynamics results allow performing statistics on each impinging helium atoms upon its spatial coordinates. The number of atoms with a z coordinate below the surface of the tungsten substrate is determined and reported to the of the substrate box. Finally, SRIM software gives directly a number of backscattered ions. The implantation probability is readily deduced from this number.

The measured implanted doses are plotted versus the impinging He atom number in figure 9. Also, are plotted the MD and SRIM simulations results. Two behaviors are observed depending on the impinging He atoms number. Before saturation, i.e. below impinging 10 16 He cm -2 , the implanted dose increases linearly.

Experimentally, the implanted dose corresponds to 83% of the number of impinging He atoms calculated from the ion current measured on the substrate holder. This ratio is higher than the calculated values of 56% for SRIM and 45 to 49% predicted by the MD simulations depending on the dose. experimental high retention rates were also evidenced experimentally in the literature but at higher energies, thus far from our conditions. E Zenobia et al 32 performed implantations with low impinging He atoms number (5x10 16 He cm -2 ) with a kinetic energy of 30 keV, they observed a retention close to 100%. Meanwhile, Gilliam et al 11 considered the retention to be equal to 100% below the "critical dose" at which they observed the formation of blisters during implantations of 1.3 MeV He ions impinging tungsten.

On the other hand, MD simulations are considering the crystallinity of the sample, moving substrate atoms and helium accumulation contrary to SRIM calculations which assume an amorphous substrate with W density but is taking into account electronic effects. So retention values are not expected to be similar. MD implanted dose is even lower than the one calculated by SRIM, due to counting He atoms liberated from bubbles too close from the surface.

A reason which could explain such differences with experiments of the implanted dose for low impinging He atoms number could be the nature of real samples which are polycrystalline. Various crystal orientations are present promoting or not the channeling and grain boundaries exist where He accumulation is different from the grain "bulk". It has been recently demonstrated that (111) orientation of W allows deeper He implantation. 33 At the same time it is reported that "Helium clustering behavior and bubble formation are largely independent of the surface orientation at depths below about 3 nm". 33 So our (100) simulations will provide a relevant picture of bubble growth leading to the substrate fracture (fuzz phenomenon). Various experimental and simulation studies presented the effect of the grain boundaries on He retention in W. On the one hand, Sefta et al. 34 MD study compared the helium implantation processes happening either there were or not grain boundaries.

They showed the propensity to accumulate in those locations which could explain the difference between the experimental values and the simulated ones. On the other hand, the comparison of various grain size made by Zenobia et al 32 , mentioned in the previous section, showed the Helium retention in nano-grain sized tungsten was lower than in polycrystalline W with micrometric grain size, but at considerably higher energies, compared to the present work. They assumed the density of grain boundaries to be responsible for the low retention level. Hu et al carried out large scale MD simulations on a model polycrystalline W at 933 K. 35 They observed that helium segregates towards grain boundaries as well as the outer surface. Thus very different results are reported in the literature regarding the effect of grain boundaries, but also in different situations. In the present experimental work W, substrates exhibits very large grains with small area of grain boundaries, thus these boundaries may have a limited impact on the He retention process.

Moreover, He atoms were found to stop migrating as soon as a small cluster of around 7 atoms is formed, at 300 K. At this temperature, they cannot reach the boundaries anymore and the global implantation process can be considered as implantation in a single crystal grains, which is described in the present MD simulations.

Above an impinging He atoms number of 10 16 He cm -2 , a leveling-off of the implanted dose appears, defining a saturation regime. This trend is observed both experimentally and using MD simulations at an implanted dose slightly above 5.10 15 He cm -2 . The fact that no saturation step is predicted by SRIM shows again that this phenomenon is related to the effect of accumulation of He in the W lattice and to the interaction between implanted He atoms. Such saturation is expected when invoking He mobility and the resulting ease to escape away. This behavior was observed experimentally by Lhuiller et al 36 , who assumed that the low implanted doses measured NRA resulted from the high mobility of helium in tungsten lattice. This was also observed earlier at low temperature by Soltan et al 37 and J. Amano et al. 29 Lee et al. 12 claimed the saturation would happen because of the stress fields increasing in the lattice due to the growth of a localized dense array of HenVm complexes preventing the deep penetration of the He ions and promoting the cluster ruptures and the ejection of the He atoms outside the material.

Our MD results attributing saturation to a flaking mechanism do not contradict these hypotheses, but focus on the role of "exploding" (or expanding) bubbles in liberating large amount of He. This can be related to highly strained tungsten lattice. Strain can originate from bubble growth, due to merging clusters which leads to swelling of the W lattice 38,39 .When an helium atom implanted dose threshold is reached (depending on incident energy and experiment temperature), the flaking and exfoliation occur 10,11 . He diffusion and He atom escaping the substrate observed in the near-surface zone (zone A) contributes also to the saturation phenomenon, the way Lhuillier suggested it 36 . On the other hand, as mentioned before, as soon as helium atoms start forming clusters, their mobility dramatically decreases even close to the surface 40,41 . This also leads to the bubble growth, expansion and flaking phenomenon when these clusters merge with the deeper ones.

One should be careful in directly comparing experiments and MD simulations: experiments are carried out on polycrystalline W while present MD simulations are considering a perfect single crystal W. Nevertheless, the MD study of Sefta et al. 34 mimicking He interactions on model polycrystalline W at 60 eV show that polycrystalline state only affect bubble size distribution and location for the same implanted dose. Flaking then should not be affected by features when the bubbles coalesce close to saturation level. Moreover experiments of Gilliam et al. 11 showed that the crystalline state only plays a role for low dose while at higher dose, i.e. close to observed saturation, there is no more difference between single crystal and polycrystalline W substrates.

Conclusion

Molecular dynamics and SRIM simulations are carried out to give an insight into helium behaviour during implantation process into tungsten leading to bubble formation.

Experiments are compared to MD simulations, paying attention to the close matching of initial conditions. MD simulations describe three steps in the helium behaviour:

-He penetrates the substrate and diffuses interstitially, until it encounters another He atom.

-Two He atoms form a cluster which stops diffusing and become a trap site for other incoming He atoms. The stratification observed on the calculated depth profiles is due to the formation of clusters at different depths of the impinging helium atoms and on the way by which the energy is dissipated in the substrate.

-Formed He clusters grow and coalesce beneath W substrate surface until He concentration becomes high enough to lead to tungsten flaking. On the other hand, we do not observe near-surface "cluster rupture" leading to only He atom ejection towards vacuum 300 eV 3 He implantation experiments in polycrystalline W samples corroborate the impinging He atom number and the retention doses at which the saturation phenomenon occurs, which is found around 5.10 15 He cm -2 . This is also consistent with previous experimental studies at low energy. However, below a certain impinging He atom number, there is a difference between the retained doses with the two methods, simulation and experiments; this is attributed to the polycrystalline nature of the samples. Further work will allow studying influence of the impinging helium atom kinetic energy and substrate temperature. An interesting point would be to introduce a kinetic energy distribution function more consistent with the experiments and include grain boundaries in the simulation box.

A

  (100) body-centred cubic tungsten crystal box consisting of 17 x 17 x 204 unit cells is built, for which dimensions are then 53.8 x 53.8 x 645.7 Å. Periodic boundaries conditions are set in the x [100] and y [010] directions, and free motion is allowed along z [001]. The temperature of the substrate is set to 300 K, thus matching the room temperature experiments. The bottom of the box is composed of a fixed layer of 6 unit cells (20.5 Å) which mimics the unperturbed bulk. The rest of the W atoms are maintained at the temperature 300 K by the mean of a Berendsen thermostat 18 . The thermostat is used because the fluxes of impinging atoms used in MD simulations are unphysically high as compared to the experiments in order to reduce the large calculation times and properly accounting He interactions with W and already He present in W. Then over a set time, velocities of the W atoms are re-evaluated and brought back to values matching the initially set global temperature of the substrate. The substrate atom W current velocities are modified using the multiplicative coefficient :

2

 2

Figure 1 :

 1 Figure 1: SEM images of the polished W samples surfaces magnified 5000 times (a)prior and (b) after 1600°C annealing.

  to 2l. It displays the evolution of the substrate, the formation of He clusters, growing and forming bubbles.The vertical dotted line gives the initial z=0 position of the W surface. When increasing injection of He atoms, swelling of the surface is observed and thus the outermost W surface is evolving towards higher z values. The bottom of the substrate is truncated at a depth of 150 Å for clarity, but extends to 645 Å for MD simulations.

Figure 2 :

 2 Figure 2: Snapshots of the implantation process 9971total impinging He atoms at 300 eV and at 300 K. The W atoms are figured with small dots and He atoms with large balls for clarity. (a) The substrate before any impinging atoms reaches it. (b-e) Helium atoms diffuse and form clusters, the substrate surface moves toward higher z. (f-j) clusters number and size grow with a preferential zone (B) at the mean penetration depth of the Helium atoms. (k-l) Clusters in the intermediate zone (B) coalesce which leads to the rupture of the substrate in two parts: the flaking.

  ). Two typical trajectories are presented in figures 3a and b. In both cases the mobility of the growing clusters is reduced as soon they reach a size of 3 He atoms.

Figure 3 :

 3 Figure 3: Two different trajectories of He atoms leading to cluster growth in Big

Figure 4 :

 4 Figure 4: Coalescence of two He clusters at times (a) 600 ps (b) 720 ps and (c) 800 ps

Figure 5 :

 5 Figure 5: The overall surface is expanded displaced W atoms, due to bubble growth,

Figure 6 :

 6 Figure 6: details of the substrate flaking. (a) He bubbles are ready to coalesce on the width of

1 .

 1 Final depth profileIn Figure7(a) a of the helium atom positions is presented for 2000 impinging and 901 retained atoms. While in figure7(b), the helium depth profiles determined both by SRIM and MD simulations are potted.

Figure 7 :

 7 Figure 7: a) Snapshot of the MD simulated He atoms distribution. The W atoms have been removed for clarity; b) Corresponding plot of the depth distribution comparing MD simulations (full line) and the SRIM calculations (dotted lines). The area labelled A, B and C define three different implantation zones.

  Fig 8 for various number of impinging He atoms.

Figure 8 :

 8 Figure 8: Depth profile plots as a function of number of impinging He atoms . The vertical dotted line displays the initial surface at t=0 and the full one corresponds to the actual swelled surface. The separating lines show the evolution of the B and C zones. Each depth profile corresponds to a frame presented in fig 2. Labels on left part of the plot correspond to snapshots of figure 2.

Figure 9 :

 9 Figure 9: Implanted dose (measured or calculated) as function of the number of impinging He atoms per unit area. The empty diamonds are the experimental measurements of the dose obtained by NRA. The full triangles are results from the molecular dynamic simulations. The full line corresponds to the 100% implantation probability and the dotted line corresponds to the implantation probability obtained by SRIM calculations (Pi=0.56).

Table 1 :

 1 MD and SRIM results for the retention rates of He ions incident on W normal to the surface.

	Number of He	Retention rate
		MD Simulations	SRIM calculations
	2000	0.45	0.560
	4000	0.430	0.56
	5000	0.38	0.56
	9971	0.2	-
	10000 (2x larger substrate)	0.48	0.56

Table 2 :

 2 Number of injected and retained He atoms for various MD simulation runs

	Run	Impinging	Corresponding	Retained	Implanted dose	Retention
	#	atom	dose(He cm -2 )	atom number	(He cm -2 )	rate
		number
	1	2000	6.91x10 15	906	3.13x10 15	0.45
	2	4000	1.38x10 16	1721	5.94x10 15	0.430
	3	5000	1.73x10 16	1887*	6.52x10 15	0.38
	4	9971	3.37x10 16	1939*	6.70x10 15	0.2
	5	10000 (large)	1.73x10 16	4754*	8.24x10 15	0.48

(*) Substrate rupture occurred at this number of retained He atoms.

Impinging He number (at.cm -2 ) NRA Measure MD calculations 100% SRIM-Pi=0.56 10 15 10 15 10 17 10 16 10 16 10 18

	(at.cm -2 )
	Dose
	Implanted

Acknowledgements

Conseil Régional du Centre is gratefully acknowledged for granting "PIMAT" project. We would like to thank C. S. Becquart and J. Boisse for providing us with modified He-W potentials.