Manoj Ducher
email: manoj.ducher@impmc.upmc.fr

Fabio Pietrucci

Etienne Balan

Guillaume Ferlat

Lorenzo Paulatto

Marc Blanchard

Van der Waals contribution to the relative stability of aqueous Zn(2+) coordination states

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The reactivity of zinc in aqueous solutions has important consequences in many scientific fields, including biology, geochemistry, or material sciences. Zinc most often occurs in its ionic form Zn 2+ with [Ar] 3d 10 4s 0 electron configuration. This filled d shell configuration does not exert a strong control on the coordination state and Zn has been shown to occur in 4-, 5-or 6-fold coordination states in various organic (e.g. proteins) or inorganic phases.

In water, experimental studies point to an average coordination number of 6 at room temperature, with a first solvation shell made of H 2 O molecules only. [1][2][3][4][5][6][7][8][9][10][11] In comparison, early theoretical studies of the aqueous environment of Zn based on molecular clusters, 3,[12][13][14][15][16] lead to contrasted results, depending on their theoretical level (e.g. MP2 vs. hybrid-DFT) and numerical parameters (e.g. basis quality). 13,14,16 With increasing computational resources, larger systems were simulated using classical, 11,17,18 hybrid (QM/MM) [19][20][21] and DFT based 22,23 molecular dynamics (MD) and pointed to a Zn coordination number of 6, consistent with experimental observations. However, a recent theoretical investigation of aqueous Zn coordination, 24 using the GGA-BLYP functional and based on a thermodynamic integration of forces, suggests that at room temperature, tetraaquo and pentaaquo-zinc complexes have lower free-energies than hexaaquo-zinc complexes. These results thus significantly differ from the other investigations performed at a similar theoretical level.

In the present study, we address this question by combining long-time FPMD simulations with free-energy calculations and by using either a standard GGA or a van der Waals (vdW) included functional. The comparison of the GGA and vdW structural results with X-ray absorption experiments allows to highlight serious structural flaws in the GGA case as a result of an overestimation of 4-fold coordination states. In addition, we show that the energetics is crucially dependent upon the choice of the functional: the energetics' order of the various complexes is actually reversed from GGA to vdW and it is consistent with experiments only in the latter case. Thus, we attribute the inconsistencies in the previous controversial results using standard GGA simulations of Zn 2+ coordination in aqueous solutions, 24 to the neglect of van der Waals interactions.

Methods

The calculations presented here are carried out within the DFT framework, 25,26 using planewaves, periodic boundary conditions and pseudopotentials. Even though DFT is intended to solve exactly the Schrödinger equation in its theoretical conception, its practical implementation requires some approximations. In particular, an exact analytical expression for the electronic exchange-correlation energy is unknown and several approximate formulations exist. The widely used one in condensed matter physics is the Perdew-Burke-Ernzerhof (PBE) parameterisation of generalised gradient approximation (GGA) functional. 27 By construction, GGA functionals essentially miss van der Waals effects and recent studies emphasised on their importance. [28][29][30] DFT calculations with non-local functionals are intended to provide a better description of van der Waals interactions through the addition of a non-local correlation term E nl C to the exchange-correlation energy E XC :

E XC = E GGA X + E LDA C + E nl C (1)
where

E GGA X is the GGA exchange energy, E LDA C
is the LDA correlation energy. Energetic and structural improvements are observed using non-local functionals especially in water. [31][32][33][34] Here we use the vdW-DF2 non-local functional 31 to evidence improvements over standard GGA calculations in this particular system.

The simulated system contains 1 Zn atom and 64 water molecules in a cubic box of length 12.65 Å so as to reproduce the experimental density of water (1.0 × 10 3 kg/m 3 ). In the initial configuration, the Zn coordination number is 6. The oxidation state of Zn being +2, a negative compensating uniform background charge is applied. All calculations are performed using the PWscf code included in the DFT package Quantum ESPRESSO, 35 in which Born-Oppenheimer type FPMD is implemented. The time step used, 0.97 fs, is chosen to be smaller than the smallest vibrational period in the system (here the O-H stretching period of 9-10 fs). The simulations are performed in the NVT ensemble, the temperature being maintained at 300 K using a Berendsen thermostat 36 with a coupling time constant of 0.3 ps. Ultrasoft pseudopotentials for Zn, O and H from the GBRV library are used for the description of the ionic cores. 37 The cut-off energies for wave-functions and charge-density expansions in finite basis-sets are fixed at 50 Ry and 450 Ry respectively. Given the large system size, the Brillouin zone sampling for electronic integration is restricted to the Γ-point.

The instantaneous evolution of Zn 2+ coordination is investigated with FPMD trajectories of 30 ps length, amongst which the equilibration and production lengths are 6 and 24 ps respectively.

The thermodynamic properties of the system are then investigated by computing the 4

where the individual contributions s Zn-O i are calculated using a switching function of the interatomic distances (r Zn-O i ):

s Zn-O i = 1 - r Zn-O i ro n 1 - r Zn-O i ro m . (3)
The parameters r 0 , n and m and were set to 2.8 Å, 10 and 20 respectively in order to have a smooth function decaying from 1 to zero between the first and second peak of the Zn-O radial distribution function (see Fig. S1 in the Supporting Information). This choice leads to non-integer numbers for the 3 characteristic coordination numbers (Fig. 1). Nevertheless, in the following, we will describe the coordination states of Zn using CN=4, CN=5 and CN=6 notations, which correspond to the presence of 4, 5 and 6 water molecules within a radius of 2.8 Å from Zn 2+ respectively. Metadynamics and umbrella-sampling simulation lengths are 180 ps and 100 ps respectively and are performed using the plugin PLUMED, 40 which is an open source library for free-energy calculations.

Metadynamics explore the free-energy landscape by depositing a bias potential V B (s, t) in the CV space that facilitates the crossing of free-energy barriers. By doing so, configurations having a low probability of being accessed in FPMD, can be explored in metadynamics in a reasonable time (few ps). The bias potential is obtained by depositing every 60 MD steps, Gaussian hills of height 2.6 kJ mol -1 and width of 0.07 CV unit. The profile of the freeenergy F (s) is estimated from the opposite of the final bias potential. The latter is obtained 5

Page 5 of 24

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation by adding the Gaussians deposited during the trajectory at each value of the coordination number, as the time-average of V B (s, t) between the time after which the system has explored all the relevant CV space (t 1 = 86 ps) and the time before which important distortions occur due to the exploration of irrelevant high-energy configurations (t 2 = 138 ps) 41 :

F (s) ≃ - 1 t 2 -t 1 t 2 t 1 V B (s, t) dt . (4)
By comparing the estimates in the two halves of the interval [t 1 , t 2 ], we expect a statistical error bar of ±5 kJ mol -1 (Fig. S2 in the Supporting Information).

In umbrella-sampling, instead, the CV space is cut into 20 windows and a harmonic restraining potential of type k(s -s c ) 2 with a spring constant value (k) of 0.594 Ry is added at the centre (s c ) of each window. k was chosen to give an ideal overlap of 0.9 between consecutive probability distributions. For each window, a trajectory of 5 ps (with the first ps discarded as equilibration) is performed using a configuration extracted from metadynamics as initial structure. The free-energy profile is reconstructed by the weighted histogram analysis method (WHAM). 42 Free-energy profiles calculated with our in-house WHAM code are similar to those calculated with Grossfield [START_REF] Grossfield | the weighted histogram analysis method[END_REF] and GROMACS 44 codes. By comparing the profiles estimated using only either half of the production trajectory in each window, we expect a statistical error bar of ±2.5 kJ mol -1 (Fig. S3 and S4 in the Supporting Information).

Theoretical EXAFS spectra have been computed from the MD configurations using the Feff8 code [START_REF] Ankudinov | [END_REF] and the methodology detailed in Ref. 46 Snapshots were extracted every 10 MD steps (after equilibration) from each FPMD trajectory (GGA and vdW-DF2), resulting in two series of about 2,400 configurations. For each configuration, an individual EXAFS spectrum was calculated with SCF potentials within a radius of 3.0 Å around the absorbing atom (Zn), a half-path distance of 7.0 Å and considering all multiple-scattering paths up to 7 legs. The amplitude reduction factor S 2 0 was set to 1, the value used for the experimental extraction. 47 An average calculated spectrum (χ calc ), directly comparable to the experimental one (χ exp ), is then obtained from the average over the individual spectra. The agreement between the calculated and experimental spectra is estimated by δ, defined as a sum of squared differences between the k 3 -weighted spectra:

δ = k k 3 χ calc (k) -k 3 χ exp (k) 2 (5)
In the following comparisons, two different experimental spectra taken from two recent studies, 22,47 are used in order to illustrate the error bar associated to the EXAFS signal extraction. These were acquired in diluted solutions of zinc-containing salts (0.1 M Zn(NO 3 ) 2 in Ref. 47 and 0.4 M Zn(ClO 4 ) 2 in Ref. 22 ) so that the contribution from counter-ions in the Zn first neighbours shell can be safely neglected.

Results and discussion

FPMD simulations

The instantaneous Zn coordination number, calculated using Eq 2 and 3, is shown along both the GGA-PBE and vdW-DF2 trajectories in Fig. 1 (after the 6 ps equilibration time).

Representative snapshots of the different Zn complexes observed along the FPMD trajectories are also shown (Fig. 1).

In the GGA-PBE case, after a relatively long duration (12 ps) over which the Zn coordination is close to 6 (hexaaquo complex), it suddenly falls to 5 and then to ∼ 4 a few picoseconds later. The tetraaquo complex then formed remains stable until the end of the GGA-PBE trajectory, without any return to the CN=6 state (Fig. 1). This suggests that, within the GGA scheme, tetraaquo complexes are more stable than hexaaquo ones and/or the presence of an energy-barrier between CN=4 and CN=6. This is well in line with previous GGA-BLYP simulations (using both long FPMD simulations of 60 ps and free-energy calculations), 24 which suggested that all 3 complexes co-exist at 300 K, with tetraaquo com- ). Also added are the results obtained from the portion of the GGA-PBE trajectory where the coordination number is 6 (CN=6) and the FPMD (GGA) results from Liu et al. 24 and Cauët et al. 22 The contrasted average Zn speciation between the two trajectories is reflected in the zinc-oxygen radial distribution functions (RDF) g Zn-O (r) and in the corresponding running coordination numbers n Zn-O (r) (Fig. 2a). The RDF first peak is located at r Zn-O = 2.02 and 2.14 Å for GGA-PBE and vdW-DF2, respectively. The latter value is within the experimental range of observed Zn-O bond lengths, from 2.06 9 to 2.15 Å, 7 while the GGA-PBE value is too short. In addition, for the GGA-PBE case, RDF specific to each coordination type g CN Zn-O (r) (CN= 4, 5, 6) have been calculated by using the corresponding portions of trajectory (e.g. from 16 to 30 ps for CN=4, see Fig. 1): the first peak location is 1.97, 2.05 and 2.11 Å for CN=4, 5 and 6, respectively. This shows that the origin of the too short RDF first peak location obtained with GGA-PBE is the overestimation of contributions from the CN=4 and CN=5 complexes. Further, the g CN =6 Zn-O (r) obtained with GGA-PBE is quite comparable to the vdW-DF2 g Zn-O (Fig. 2b): this indicates that although modifying the energetics of the system, the inclusion of van der Waals interactions has only a minor effect on the geometry of the hexaaquo complex. We also add in the comparison (Fig. 2b) the results from previous GGA studies. The data from Cauët et al. 22 (who observed only hexaaquo complexes) are quite similar to our GGA-PBE results restricted to the CN=6 trajectory and to the vdW-DF2 full trajectory (Table 1) while the data from Liu et al. 24 (who pinpoint the existence of pentaaquo and tetraaquo complexes in however unreported proportions) show a shorter RDF first peak at 2.05 Å (and a running coordination number of 5.2 at 2.8 Å). This shows that the position of the RDF first peak is a reasonably good marker of the average coordination number, mostly irrespectively of the functional used.

The GGA-PBE and vdw-DF2 scenario are now tested by comparing EXAFS calculated spectra to experimental ones 22,47 (Fig. 3). The spectrum obtained from the GGA-PBE trajectory shows significant discrepancies with the experimental ones (Fig. 3a): in particular the too large period (in k space) essentially reflects a too short mean Zn-O distance (typically 2.02 Å versus 2.078 Å experimentally 47 ). A significant improvement is provided with the vdW-DF2 scheme, the residual between the spectra (δ) being reduced by a factor of about 3-4 (Fig. 3b). A slight frequency misfit is still visible, reflecting in this case a slightly too long mean Zn-O distance. As described previously for the RDF, the GGA-PBE trajectories can be decomposed according to the Zn coordination number (6, 5 and 4), allowing to produce representative spectra for each type of complex (panels c, d and e of Fig. 3). Clearly, the spectrum obtained for hexaaquo complexes (Fig. 3c) shows a pretty good match of the experimental spectra; typically comparable (or slightly better) than that obtained with the vdW-DF2 full trajectory (which corresponds mostly to hexaaquo complexes as well). By contrast, the δ residual significantly increases for lower coordination types, with a severe frequency misfit most pronounced in the tetraaquo case (Fig. 3e). This decomposition thus allows confirming the main origin of the GGA-PBE deficiency, namely a significant overestimation of contributions from 4-fold and 5-fold complexes. While a quantitative assessment of the relative proportion of each complex from the EXAFS spectra is beyond the scope of the present paper, it seems reasonable to discard the presence of tetraaquo complexes in These results are consistent with those of Liu et al. 24 (also reported in Fig. 4 and Table 2), obtained by thermodynamic integration of forces using the GGA-BLYP functional. Beyond the overall good agreement between the techniques, our results show a more pronounced barrier height between CN=4 and CN=5 (17-18 kJ mol -1 ) as compared to those of Liu et al. 24 (9 kJ mol -1 ), which is likely a result of the different GGA functionals used (PBE versus BLYP). In any case, three distinct techniques used to calculate the free-energy profiles, provide the same overall picture (which is in contradiction with experiments). It shows that within the GGA scheme, unbiased MD simulations starting from equilibrated tetraaquo or pentaaquo complexes have a small probability to sample the hexaaquo state since it requires to overcome a barrier of at least 12 kJ mol -1 . The reverse situation (sampling CN=4 or 5 starting from CN=6) is more likely to happen since the corresponding barriers are lowered to 2-6 kJ mol -1 . As a matter of fact, GGA FPMD studies with relatively short trajectory lengths (∼ 15 ps), 22,23 found only hexaaquo complexes while longer trajectories (≥ 30 ps) display the 3 complexes (this work and Liu et al. 24 ).

Table 2 Summary of the free-energy differences (∆F ) between hexa, penta and tetraaquo-zinc complexes and associated barrier heights calculated from metadynamics using the GGA-PBE functional, and from umbrella-sampling using either the GGA-PBE or the vdW-DF2 functional. The results are compared with those of Liu et al. 24 obtained from free-energy calculations using the GGA-BLYP functional and thermodynamic integration of forces. Values are given in kJ mol -1 .

Metadynamics

Umbrella sampling Umbrella sampling Thermodynamic Integration using GGA-PBE using GGA-PBE using vdW-DF2 using GGA-BLYP, from Ref Back to our results in closer details, the free-energy profiles obtained from metadynamics and umbrella sampling agree within 10 kJ mol -1 (Fig. 4), as expected from other similar simulations of chemical reactions in solution. 48 However, two differences regarding the location of the extrema are noticed: first, the coordination number minimum is shifted from According to the present free-energy calculations, van der Waals interactions are key ingredients to stabilise hexaaquo zinc state with respect to tetraaquo and pentaaquo states at the GGA level. However, in a previous study of aqueous Zn molecular clusters, Diaz et al. 16 have pointed out that a stabilisation of hexaaquo zinc complex with respect to the tetraaquo one could also be obtained by improving the description of the electronic structure (using calculations at the MP2 level). To further document the respective contribution of electronic exchange and non-local correlation, we have performed additional calculations on the same Zn clusters using hybrid functionals (PBE0 and B3LYP) as well as GGA and vdW functionals (Table S1 in the Supporting Information). On these complexes GGA functionals predict the tetraaquo complex to be the most stable. A better description of the electronic exchange as treated by hybrid functionals leads to the same picture with however a relative stabilising effect on the hexaaquo complex. Only the inclusion of non-local correlation leads to the expected energetic order between the two complexes. Although using a simplified Using FPMD simulations, free-energy calculations and simulated EXAFS spectra, we have shown that the description of aqueous Zn solutions as obtained from standard GGA functionals suffers from a severe overestimation of the stability of tetraaquo complexes (at the expense of the hexaaquo complexes). Further, the origin of the contrasted results obtained from different FPMD GGA studies [22][23][24] has been unveiled. Indeed, because of the presence of large energy barriers, the simulated system can remain kinetically trapped for relatively long times in the (metastable) hexaaquo state: such a situation will fortuitously give a good account of the experimental situation although being numerically unconverged, leaving the problem unnoticed. 22,23 The incorrect GGA behaviour is revealed by long trajectories and/or enhanced-sampling techniques (here metadynamics and umbrella-sampling) and by a direct confrontation with experiments. The MD-EXAFS combination used here provides a clear enough test of the GGA insufficiencies.

Strikingly, the use of a recent functional including van der Waals contributions, vdw-DF2, changes drastically the relative energetic stability of the aqua complexes, providing a picture fully consistent with experiments in which Zn appears mostly as a hexaaquo complex. This shows that even weak interactions such as the van der Waals forces must be accounted for a reliable description of the chemical and thermodynamical properties of Zn-bearing aqueous solutions. Although the importance of these contributions has been progressively recognised in the recent years in other systems, we are not aware of any other cases in which these forces induce such a marked and qualitative change. From the structural point of view however, it should be noticed that the local geometry of the 6-fold coordination complex is only weakly affected by the use of vdW-DF2 as compared to standard GGA.

The present results may have implications for other hydrated cations, in particular regarding the relative stability of various coordination states. It can be anticipated that the use of non-local functionals incorporating van der Waals contributions could improve the description of main group ions, 50,51 such as Ca, Na and K, which play important biological

6 Page 6 of 24 ACS

 624 Paragon Plus EnvironmentJournal of Chemical Theory and Computation

7 PageFigure 2

 72 Figure2(a) Zinc-oxygen radial distribution function g Zn-O (r) and running coordination number n Zn-O (r) obtained from the whole GGA-PBE and vdW-DF2 trajectories. (b) Zoom on the first peak region (1.7-3.0 Å). Also added are the results obtained from the portion of the GGA-PBE trajectory where the coordination number is 6 (CN=6) and the FPMD (GGA) results from Liu et al.24 and Cauët et al.22

11 Page 11 of 24 ACSFigure 3

 11243 Figure 3 Calculated average EXAFS spectra of molecular configurations extracted from (a) the whole GGA-PBE, (b) the whole vdW-DF2 FPMD trajectories, (c) the portion of GGA-PBE trajectory where CN=6, (d) the portion of GGA-PBE trajectory where CN=5 and (e) the portion of GGA-PBE trajectory where CN=4. These calculated spectra are compared with the experimental spectra from d'Angelo et al.47 and Cauët et al.22

14 Page 14 of 24 ACS

 1424 Paragon Plus EnvironmentJournal of Chemical Theory and Computation tween hexaaquo and pentaaquo complexes) is reversed from 6 to -4 kJ mol -1 . The value of the barrier height between CN=6 and CN=5 is 11 kJ mol -1 , explaining why in vdW-DF2 FPMD trajectory, hexaaquo complexes are mostly observed. Thus the inclusion of vdW forces induces an inversion of the energetic stability of the Zn-aqua complexes which is at the origin of the stabilisation of hexaaquo complexes in FPMD trajectory as evidenced in the previous section (3.1). However, the absolute value of |∆F 6-5 | is relatively small (as compared to |∆F 6-4 |), thus the presence of a fraction of pentaaquo complexes in aqueous Zn solutions cannot be ruled out neither from the energetic calculations nor from the EXAFS spectroscopic measurements.

  cluster model of Zn environment in water, these additional results point once again towards a dominant role of van der Waals contribution in the relative stability of aqueous Zn 2+

Page 4 of 24 ACS Paragon Plus Environment Journal of Chemical Theory and Computation free

  -energy profiles associated to the changes in Zn coordination with oxygen atoms. These calculations use metadynamics38 and umbrella-sampling39 strategies which are based on a exploration of the configuration space restricted to the vicinity of a path defined by a

	collective variable (CV). Here, the CV is the coordination number (CN) s,
	64
	s =	s Zn-O i
	i=1

Table 1

 1

	summarises structural data extracted from GGA-PBE

Page 9 of 24 ACS Paragon Plus Environment Journal of Chemical Theory and ComputationTable 1

 1 Structural information extracted from GGA-PBE and vdW-DF2 FPMD trajectories. The GGA-PBE trajectory is decomposed in three portions as shown in Figure1depending on the observed Zn coordination number (CN) and refered below as CN=6, CN=5 and CN=4. Information on Zn-O radial distribution functions g Zn-O (r) as well as the number of water molecules n H 2 O , average Zn-O distances r Zn-O in the first and second solvation shells are reported. Average OZnO angle < OZnO , Zn-H distance r Zn-H , O-H bond length r O-H and HOH angle < HOH within the first solvation shell are also given. Root mean square σ of these quantities are also reported. Some of these quantities reported by Cauët et al.22 are given for comparison.

			GGA-PBE vdW-DF2	GGA-PBE CN=6	GGA-PBE CN=5	GGA-PBE CN=4	Ref 22
		g Zn-O (r) 1 st max 2.02	2.14	2.11	2.05	1.97	2.12
		g Zn-O (r) 1 st min 2.8	2.8	3.0	2.6	2.35	2.8
		n H2O	4.70	5.98	5.97	5	4	6
		σ n H 2 O r Zn-O	0.84 2.13	0.14 2.19	0.16 2.17	2.10	1.99	2.18
	1 st solvation sphere	σ r Zn-O < OZnO σ < OZnO r Zn-H	0.20 98 12 2.67	0.13 90.1 7.7 2.78	0.15 89.9 7.9 2.77	0.12 96 11 2.67	0.07 101 13 2.58	0.11 2.76
		σ r Zn-H	0.16	0.15	0.17	0.14	0.12	0.17
		r O-H	1.010	0.996	1.000	1.008	1.019	0.97
		σ r O-H	0.033	0.027	0.030	0.033	0.033	0.03
		< HOH	107.2	105.9	106.0	106.9	108.3	106.6
		σ < HOH	5.4	5.0	5.0	5.3	5.5	5.5
		g Zn-O (r) 2 nd max 3.9	4.3	4.2	4.0	3.9	4.2
	2 nd solvation sphere	g Zn-O (r) 2 nd min 4.9 n H2O 12.2 1.2 σ n H 2 O r Zn-O 4.12	5.1 14.7 1.5 4.37	4.9 12.0 1.1 4.33	4.7 11.84 0.93 4.10	4.5 9.82 0.92 3.95	4.7 14 4.22
		σ r Zn-O	0.33	0.33	0.30	0.25	0.27	0.27
				10

Page 10 of 24 ACS Paragon Plus Environment Journal of Chemical Theory and Computation

ACS Paragon Plus EnvironmentJournal of Chemical Theory and Computation

ACS Paragon Plus Environment Journal of Chemical Theory and Computation ACS Paragon Plus Environment Journal of Chemical Theory and Computation Acknowledgement The authors thank P. d'Angelo for providing us with the EXAFS spectra of Ref. 47 This work was performed using HPC resources from GENCI-IDRIS (Grants 2016-i2016041519 and x2017081875). This work was granted access to the HPC resources of The Institute for scientific Computing and Simulation financed by Region Île de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" programme. This work was supported by the French national programme EC2CO-Biohefect. ACS Paragon Plus Environment Journal of Chemical Theory and Computation

Supporting Information Available

This material is available free of charge via the Internet at http://pubs.acs.org/.

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation


Bacroix

S Dufrenoy

T Chauveau

I Lemaire-Caron

R Brenner

B Bacroix

Keywords: residual stresses, hcp structure, titanium alloys, self-consistent model, X-ray diffraction (XRD)

Comparison of 2 methodologies developed for the determination of residual stresses through X-ray diffraction: application to a textured hcp titanium alloy

Introduction

Residual stresses are stresses which are retained within a body when no external forces are acting.

They arise because of misfits (or incompatibilities) between different regions of the investigated material sample or component, which are a direct consequence of what these have been subjected to. Residual stresses can be classified according to their origin (e.g. [START_REF] Withers | Residual stress and its role in failure[END_REF][START_REF] Bretheau | Les contraintes résiduelles : d'où viennent-elles ? Comment les caractériser ?[END_REF]), the scale over which they equilibrate (e.g. [START_REF] Noyan | An X-ray diffraction study of the residual stress-strain distributions in shot-peened two-phase brass[END_REF]), or their effect on the macroscopic behavior of a structural piece (e.g. [START_REF] Turner | A study of residual stresses in Zircaloy-2 with rod texture[END_REF]).

As they can have beneficial or detrimental effects on the life time of mechanical components, indeed, the three aspects (origin, distribution and influence on mechanical behavior) are important.

In the case of polycrystalline materials, the most common origins of the misfits are non-uniform plastic flow, steep thermal gradients or phase transformations. And due to the polycrystalline nature of the materials, these stresses indeed vary from point to point within the microstructure, and thus, several orders of stresses are generally distinguished [START_REF] Bretheau | Les contraintes résiduelles : d'où viennent-elles ? Comment les caractériser ?[END_REF][START_REF] Baczmański | Plastic incompatibility stresses and stored elastic energy in plastically deformed copper[END_REF] (see below).

In many industrial applications, in which materials and components are subjected to severe and complex thermomechanical treatments (including forming, heat treatment, welding, machining, …), the origin of the existing residual stresses is usually a combination of several misfit sources (e.g. plasticity and thermal contraction due to rapid cooling [START_REF] Turner | Modeling the mechanical response of rolled Zircaloy-2[END_REF]), which are not easy to know a priori. The so-called macroscopic 1 st order stresses (which are averaged on a quite large volume and thus neglect the underlying microstructure) are of primary importance, since they are indeed the ones which have the greatest influence on the life time of mechanical parts of a structure. They are more and more often determined through X-Ray Diffraction (XRD), and this is why some standards (e.g. [START_REF] Afnor | Non-destructive Testing -Test Method for Residual Stress analysis by X-ray Diffraction[END_REF]) have been developed to define the best protocol to be able to evaluate these stresses from measurements with sufficient quantitative accuracy. However, all these standards rely on crude assumptions concerning the material, which are rarely verified for metals which present a hexagonal close packed (HCP) crystallographic structure (such as Zr, Ti, Mg, ….), and which are widely present e.g. in the aeronautic or energy industries. The two main important assumptions which are mentioned in these standards are that (i) the investigated material can be considered isotropic, and (ii) the 2 nd order stresses -of inelastic origin -can be neglected. Unfortunately, these two assumptions are rarely met for the above mentioned metals, as well as for dual-phase Ti alloys with complex microstructures, which are now widely used in the aeronautic industry.

The aim of the present work is thus to propose a modification of the classical methodology, which minimizes the error made on the estimation of the 1 st order stress, and thus improves the recommended standard, although it still based on a priori assumptions. In order to validate this new method and to compare it to the classical methodology, both are applied to a theoretical case for which the residual stresses are a priori completely known. It consists into the simulation of a tensile test of an anisotropic Ti sample with a polycrystalline mean field model. The paper is thus organized according to the following outline. The classical definitions of the stresses defined at different scales are first recalled in Section 1; then, the polycrystalline model, tested material and performed simulations are presented in Section 2. The classical standard and proposed methodology are detailed in Section 3. Results obtained with both methods are then presented in Section 4 and compared to those of the polycrystalline model. The paper ends up with some general recommendations about the use of these methods and perspectives.

Definition of the various order stresses for a polycrystalline material

Let us consider a polycrystalline sample, schematically represented by a 2 dimensional section in Fig. 1. This sample has been e.g. machined on the upper surface (perpendicular to the figure and identified by the direction normal to it, labelled ND), which means that it has been subjected to heterogeneous stress and temperature fields throughout its thickness. As a result, the internal structure and mechanical state within the superficial layer has been modified compared to the lower part of the sample. This induces misfits between various layers, and thus residual stresses. During such a machining operation, plastic strain and heating/cooling are the main processes responsible for the appearance of residual stresses. Due to the polycrystalline structure of the material, plastic deformation, as well as thermal dilatation/contraction, generally varies from grain to grain, producing in turn misfits between grains, i.e. at the intergranular level [START_REF] Pang | In situ generation of intergranular strains in an Al7050 alloy[END_REF][START_REF] Clausen | Lattice strain evolution during uniaxial tensile loading of stainless steel[END_REF]. And we also know that in such materials, deformation is rarely homogeneous within the grains, resulting in turn into heterogeneities at the intragranular level. The various stresses usually defined in the field of residual stresses are recalled below. As a result of inhomogeneous boundary conditions imposed to the whole sample, the stress state depends on the location. We can thus define stresses at different levels:

the local stress at point of the sample, the average stress for grain , the average stress on a given investigated volume V containing a large number of grains and finally the stress averaged on the whole sample, which can be the imposed stress Σ during a mechanical test. The investigated volume V can be for example a part of a superficial layer, whose thickness will both depend on the measurement method, as well as the investigated process. It has to be distinguished from the diffracting volume Ω which is defined latter. In the context of residual stresses, it is common to decompose these various stress as follows

= + = 1 . = + = 1 . = + = 1 . = Σ (1)
We see that , , are fluctuations around average values, which are classically called 3 rd , 2 nd and 1 st order stresses respectively. When the whole sample is unloaded, then Σ = = 0 and = . In this case, correspond to the 1 st order residual stress tensor (averaged on the investigated volume) that we want to evaluate. The 2 nd order stresses, , which correspond to fluctuations within the grains around , are then the 2 nd order residual stresses. If the whole sample has been deformed homogeneously, then, in the unloaded state, we necessarily have = = 0 and thus = . The investigated volume can be considered in this case to be a

Representative Volume Element (RVE) for the material. We will only have 2 nd order residual stresses in this case. In what follows, we will neglect the intragranular stress fluctuations ( ) and consider only 1 st and 2 nd order stresses.

Simulation of a tensile test

We start simulating a tensile test in a textured polycrystalline sample (supposed to be initially stress free), in order to estimate the stress distribution within the polycrystalline sample for various macroscopic stress levels. In other words, we create an artificial "experimental" data set, for which the 2 nd order residual stresses will be due to plasticity only. The considered material is commercially pure titanium (grade 2), whose initial texture in presented in Fig. 2. The texture is typical from rolled and annealed titanium, and is composed of two main texture components, which are the socalled tilted 0001 〈112 0〉 and tilted 0001 〈11 00〉 orientations [START_REF] Zhu | Texture Evolution and Associated Nucleation and Growth Mechanisms during Annealing of a Zr Alloy[END_REF]. This texture statistically respects the orthotropic symmetry arising from the rolling process and is far from isotropy (the maximum is larger than 6 for the {0002} pole figure). The Orientation Distribution Function is calculated from pole figures using the LaboTex ® Software and a dataset of 50000 orientations is extracted, to represent the texture for further calculations.

The original elasto-plastic SC model, initially developed by Kröner [START_REF] Kroner | Berechnung der Elastischen Konstanten des Vielkristalls aus den Konstanten des Einkristalls[END_REF] for the case of elasticity, and further extended to the treatment of thermo-elasto-plasticity by Hutchinson [START_REF] Hutchinson | Elastic-plastic behaviour of polycrystalline metals and composites[END_REF] and Turner et al. [START_REF] Turner | Modeling the mechanical response of rolled Zircaloy-2[END_REF] has been used here, since it has been shown to be well adapted to the calculation of 1 st order stresses arising from plasticity or thermal treatment for HCP materials [START_REF] Tome | Role of internal stresses in the transient of irradiation growth of Zircaloy-2[END_REF][START_REF] Gloaguen | Measurement and prediction of residual stresses and crystallographic texture development in rolled Zircaloy-4 plates: X-ray diffraction and the self-consistent model[END_REF]. Compared to other polycrystalline models used un the context of residual stresses (see e.g. [START_REF] Baczmański | Examination of the residual stress field in plastically deformed polycrystalline material[END_REF][START_REF] Jeong | Residual stress evolution due to cool down in viscoplastic metal matrix composites[END_REF]), it has also the advantage of being solely based on mechanical interactions between one grain and the surrounding material, without the addition of any fitting parameter. It is used here to simulate a tensile test along RD up to 5% total strain, i.e. in elasto-plasticity (thermal effects are not considered here). It is then possible to extract the stresses and within each grain of the material at various stages of the loading process. The use of this model is based on the following procedure:

(i)
The material is considered to be heterogeneous, i.e. composed of a discrete set of grains of various orientations composing the texture; here, the sole consideration of the distribution of orientations is taken into account and not their position;

(ii) The simulation is performed on a so-called Representative Volume Element (RVE), on which the boundary conditions of the simulated test are imposed;

(iii) The localization step allows establishing relationships between the macroscopic quantities (at the level of the RVE) and the microscopic ones (at the level of each grain) and the knowledge of the local constitutive law allows getting the unknown macroscopic quantities through the homogenization step.

Only the principal equations relating stress and strain tensors (or stress and strain rate tensors) will be recalled below, and the reader is referred to the work of Turner, Tome et al. [START_REF] Turner | A study of residual stresses in Zircaloy-2 with rod texture[END_REF][START_REF] Turner | Modeling the mechanical response of rolled Zircaloy-2[END_REF][START_REF] Turner | Self-consistent modeling of visco-elastic polycrystals: Application to irradiation creep and growth[END_REF][START_REF] Zecevic | An elasto-plastic self-consistent model with hardening based on dislocation density, twinning and de-twinning: Application to strain path changes in HCP metals[END_REF] for further details concerning this model. It is important to note that during plastic deformation, the model not only allows the quantifications of these tensors but also to calculate the grain shapes and texture evolutions due to slip on several systems within each grain. For the performed simulations, because of the limited amount of imposed strain, the grain shape will be considered to be constant, which reduces the computing time, but texture evolution, although expected to be limited, will be considered.

If the model is used in the elastic range, then the local law relating elastic strain and stress ! at the level of the grains reads

= " : = $ : (2)
With " and $ the 4 th order elastic stiffness and compliance tensors respectively.

The localization step allows to write

= % : & = ' : ( (3)
where & and ( are the macroscopic quantities (at the level of the RVE) and % and ' are the 4 th order tensors of stress concentration or strain localization respectively. The average of these 4 th order tensors over the grains contained within the RVE (represented below by the symbol 〈. 〉 ) * ) is equal to the identity tensor, since by definition we have

〈 〉 ) * = & and 〈 〉 ) * = + (4)
We can then write macroscopic constitute laws -a priori unknown -of the same form as the local ones, & = " , : (

( = $ -:& (5)
where " , and $are called the 4 th order macroscopic stiffness and compliance tensors.

Then, if we assume that each grain is a homogeneous ellipsoidal inclusion embedded into an infinite medium, the average localization strain tensor reads ' = ./ + $ * : " , 01 : 2" -" , 45 01 [START_REF] Turner | Modeling the mechanical response of rolled Zircaloy-2[END_REF] with $ * , the Eshelby tensor, which depends on the shape of the grain g and on the infinite medium.

This equation leads to the determination of all local and global quantities from macroscopic imposed boundary equations. Now, the extension of this model to the case of elasto-plasticity (EP) requires considering incremental laws. The crystal plasticity law assumes that plastic strain takes places within each grain g by slip on specific slip systems, which are potentially active when : 6 7 = 8 9 7 with :; 7 ≥ 0 (7)

6 7 is the Schmid tensor, 8 9 7 the critical resolved shear stress and :; 7 the shear rate for system s. The components of the Schmid tensor are classically defined by

= > 7 = 1 2 2? 7 @ > 7 + ? > 7 @ 7 4 (8)
with ? AAB 7 and @ AB 7 being unit vectors representing the slip plane normal and slip direction of system s and ? 7 and @ > 7 being the individual components of these vectors.

By decomposing the total strain into the sum of the elastic and plastic parts, we get ; = $ : ; + C = :; = D : ;

The Hooke's law which relates the stress rate to the elastic strain rate now writes ; = " : E ; -C = :; F = G : ;

in terms of strain and stress rates, in which Gand D , are called the 4 th order macroscopic effective moduli.

After Yoshida et al. [START_REF] Yoshida | Micromechanical modeling of the workhardening behavior of single-and dual-phase steels under two-stage loading paths[END_REF], we consider the following saturating hardening law for the slip systems 8; 9 = C ℎ 7 :;

with ℎ 7 = IJ + 1 -J K 7 Lℎ M N1 + O P Q R S TUV W X W01 (12)
In this equation,

: Y = Z ∑ |:; 7 | 7 ] M
^ is the total strain, J, ℎ M and @ are hardening parameters and 8 `ab is the reference critical shear stress of the so-called easy glide systems. The selected slip systems families and associated parameters, which have been identified on experimental tensile curves obtained by Benmhenni et al. for commercially pure titanium [START_REF] Benmhenni | Micromechanical modelling of monotonic loading of CP alpha-Ti: Correlation between macroscopic and microscopic behaviour[END_REF], are listed in Table 1.

The single crystal elastic constants are taken equal to: c 11 = 162.4 GPa, c 1f = 92 GPa, c 1h = 69 GPa, c hh = 180.7 GPa, c kk = 46.7 GPa [START_REF] Fisher | Single-crystal elastic moduli and the hcp -> bcc transformation in Ti; Zr and Hf[END_REF]. It is worth mentioning that the HCP single crystal presents transverse isotropy, which means that the Young modulus in tension is constant in the basal plane. It is also seen in Fig. 3 (left) that it varies from 10.4 GPa in the basal plane to 14.4 GPa along the <c> axis. In other words, the single crystal elastic anisotropy is limited but not completely negligible. When the SC model is used to calculate the macroscopic Young's modulus of the textured polycrystal as a function of the tensile direction, the range is still reduced to 10.8 to 12.4 GPa but small variations are now observed within the (X,Y) plane (Fig. 3 right). Thus, we can say that the material is macroscopically moderately elastically anisotropic. Now, the boundary conditions imposed to the RVE to simulate a tensile test are the following Ε ; = m Ε ; 11 ? ? ? ? ? ? ? ? o and Σ ; = m ? 0 0 0 0 0 0 0 0 o

In other words, we solely impose a macroscopic uniaxial stress tensor together with the magnitude of the tensile strain rate Ε ; 11 . The simulated curve is shown in Fig. 4 (left). The principal strains (width and thickness strains) are also presented in Fig. 4 (right). It is worth noting that this plot highlights the plastic anisotropy of the material, since these two strains would be equal for an isotropic sample. It is seen on this curve, that these two strains are indeed almost equal at the very beginning of the deformation, i.e. in the elastic range, as expected from the calculation of the Young's modulus in this range.

Apart from the initial state characterized by the point labelled 0 and which corresponds to zero stress everywhere in the RVE, 4 points have been identified along the curve, at which 1 st and 2 nd order stresses will be calculated. For that purpose, the RVE will represent the investigated volume described above (i.e. a part of a larger sample), subjected to the calculated macroscopic stress state, which will thus be equivalent to a residual stress state Σ = . It is worth noting that the classical situation in which residual stresses are usually measured corresponds actually to point 4, i.e. a state at which residual stresses due to an heterogeneous plastic strain field have been previously

introduced into the material at the level of each grain ( *s as well as at the level of the RVE ( , which thus remains elastically loaded, although being a part of an unloaded sample. For each of these points, the macroscopic stress as well as the distribution of stresses inside the grains will be extracted from the calculations. For each of these points, we can write At point 0:

= = 0 (14a)
At point 1(elastic range): ≠ 0 and = % :

with % the elastic stress concentration tensor. At point 4 (elastic range):

! = % : + *s (14f)
Thus, because of the rewriting of Eq.14c, at each step of the loading path, the stress state within grain g, can be expressed by the general expression

! = % : + v* (15)
i.e., the sum of two terms, the 1 st one associated with the macroscopic 1 st order stress state , and the second one associated with the fluctuations -of Non Elastic (NE) origin -arising from the overall thermomechanical history of the sample and related in the present case solely to elastoplastic incompatibilities (EP). However, in the context of stress measurements associated with severe processes such as machining or additive manufacturing, plasticity generally occurs concurrently with thermal expansion or phase transformations, and the term v* is then a consequence of all these processes occurring heterogeneously within the material. However, unlike in the case of elasto-plasticity, the quantitative prediction of this term is still a challenge, which necessitates complex numerical tools (see e.g. [START_REF] Jeong | Residual stress evolution due to cool down in viscoplastic metal matrix composites[END_REF][START_REF] Guillemot | 3D Heat Transfer Analysis for a Hybrid Approach to Predict Residual Stresses After Ball-End Milling[END_REF]). In the present case, reduced to elastoplasticity, , as well as *s (and thus v* can be extracted from our simulations, for all the points highlighted on the curve (in Fig. 4). Due to the definition of the 2 nd order stresses, we also have 〈 *s 〉 ) * = 0 ⟹ 〈 v* 〉 ) * = 0.

Classical method for the determination of residual stresses

As already mentioned, residual stresses are necessarily associated with residual elastic strains ensuring the compatibility of the total strain. But, as seen previously from the above detailed equations, they may simultaneously depend on elastic, thermal and plastic properties of the material. Diffraction (of neutrons or X-rays) is now a widely-used technique to measure these elastic strains. It is also worth noting that this technique is ideally suited for comparison with mean -field models (such as the one described above), since both methodologies allow the determination of the global response of a large number of grains, without explicit account of their position within the sample. Several types of experimental set-up exist to measure these strains, from classical laboratory devices to neutron sources which allow the characterization of large volumes [START_REF] Webster | Non-destructive measurement of residual stress by neutron diffraction[END_REF] or third generation synchrotron sources such as the European Synchrotron Radiation Facility (ESRF)

which provide intense beams of high energy X-rays that are much more penetrating than characteristic laboratory X-ray devices [START_REF] Daymond | A new stroboscopic neutron diffraction method for monitoring materials subjected to cyclic loads: Thermal cycling of metal matrix composites[END_REF]. The experimental devices and set-up depend also to a certain extent whether we need to examine thin films or thick samples [START_REF] Ligot | Stress analysis of polycrystalline thin films and surface regions by X-ray diffraction[END_REF]. Recently, a new device has been developed, which allows a rapid determination of a 2D X-ray image, composed of several peaks [START_REF] Dufrenoy | Influence of surface integrity on fatigue limit: Ti-10V-2Fe-3Al titanium application[END_REF]. Whatever the precision, penetration depth and experimental set up, all these devices are based on the same type of analysis, which relies on the measurement of elastic strains, averaged, not on the whole investigated volume, but on a so-called diffracting volume Ω (see Fig. 1), only composed, for a given diffracting direction, of all grains in diffracting orientation.

For a given X-ray source of wavelength λ, and an imposed diffraction vector K, characterized by the plane indices {hkl} and the two angles x and y with respect to the sample reference frame (see Fig. 5), it is possible to measure the actual lattice parameter of the screened plane family, averaged on all grains in diffraction position, i.e. those composing the diffracting volume Ω. The Bragg's law reads in this case

2〈 Oz{ 〉 | sin • Oz{ = @€ (16)
with • Oz{ the Bragg angle associated with the plane family {hkl} and the X-ray source, and 〈 Oz{ 〉 | the lattice parameter in the direction of the diffraction vector averaged on all diffracting grains. From this, we can get the elastic strain, again averaged on Ω, provided that the initial lattice parameter (under zero stress) M Oz{ is known. This average elastic strain reads

〈• 'ƒ 〉 |"… †. = 〈 Oz{ 〉 | -M Oz{ M Oz{ (17)
This strain (which indeed constitutes one single component of the strain tensor) is measured along the direction characterized by the direction (x, y), i.e. in the reference frame linked to diffraction.

We denote by A the rotation matrix between the two references frames (sample and diffraction), such that ‡′ > = ‰ z ‰ >{ ‡ z{ [START_REF] Zecevic | An elasto-plastic self-consistent model with hardening based on dislocation density, twinning and de-twinning: Application to strain path changes in HCP metals[END_REF] where ‡ z{ and ‡′ > are the components of a given tensor, expressed in the reference frames linked to the sample and to diffraction respectively. This rotation matrix reads

‰ = m Ž••yŽ••x •'@yŽ••x -•'@x -•'@y Ž••y 0 Ž••y•'@x •'@y•'@x Ž••x o (19)
And we can thus write

〈• 'ƒ 〉 |"… †. = 〈 '' ′ ! 〉 |"… †. (20)
Now, in order to see how this quantity is related to the residual stresses, we first suppose that the measured elastic strains are the sole consequence of a prior heterogeneous elasto-plastic strain field within the sample, similar to the one simulated above with the SC model. We can then calculate the same quantity with the model as the one measured by first expressing the relationship between calculated elastic strain and stress tensors for one grain g, in the diffraction reference frame (by using Eq. 15)

′ ""•. = $′ : ′ ""•. = $′ : 2%′ : ""•.

-

+ v*""•. - 4 † (21)
then by averaging this elastic strain on all diffracting grains

〈 ′ ! 〉 |""•. = 〈$′ : %′ 〉 | : ""•. - + 〈$′ : v*""•. - 〉 | (22)
and finally by extracting the component of interest

〈 '' ′ ! 〉 |""•. = 〈$′ ''˜™ ! : %′ ˜™š› ! 〉 | : z{ - ""•. + 〈$′ ''˜™ ! : > - v*""•. 〉 | (23)
Now, if we want the 1 st order residual stress to be expressed in the reference frame linked to the sample, we can write † Of course, in this equation, as in Eqs. 22 to 24, the stress tensors are not the only calculated quantities: the tensor %′ is also calculated with the model.

〈• 'ƒ 〉 |""•. = 〈 '' ′ ! 〉 |""•. = 〈$′ ''˜™ ! : %′ ˜™š› ! 〉 | : ‰ z ‰ {W W ""•. + 〈$′ ''˜™ ! : > - v*""•. 〉 |
or more globally

〈• 'ƒ 〉 |""•. = 〈• 'ƒ 〉 |""•. + 〈• 'ƒ 〉 |""•. (24a) (24b)
Assuming now that the same relationship holds between the measured strain and the 1 st and 2 nd order stress tensors underwent by the sample and constituting grains, a priori unknown, we can also

write ‡ 〈• 'ƒ 〉 |"… †. ≅ 〈• 'ƒ 〉 |"•ž. = 〈$′ ''˜™ ! : %′ ˜™š› ! 〉 | : ‰ z ‰ {W W "•ž. + 〈$′ ''˜™ ! : > - v*"•ž. 〉 |
Or more simply

〈• 'ƒ 〉 |"… †. ≅ 〈• 'ƒ 〉 |"•ž. = 〈• 'ƒ 〉 |"•ž. + 〈• 'ƒ 〉 |"•ž. (25a) (25b)
Equation (25a) thus relates a quantity which can be measured by X-ray diffraction 〈 '' ′ ! 〉 |"… †. with (i) the 1 st and 2 nd order stress tensors that we want evaluate, and v* -

(of various origins) and

(ii) some additional unknown tensor %′ (which describes the sole elastic behavior of the material).

The matrix A is defined experimentally and the stress components $′ ''˜™ ! derived from the elastic constants of the single crystal are usually taken from the literature.

Unfortunately, in the general case of an anisotropic and strongly textured material, it is difficult to get a quantitative estimation of all stress tensors from the sole measurements (there are indeed too many unknowns in both terms of Eq. (25a)); this difficulty can be overcome by making some additional assumptions based on the following considerations:

(i) The second term 〈• 'ƒ 〉 |"•ž. = 〈$′ ''˜™ ! : > - v*"•ž.
〉 | depends on the whole thermomechanical history of the material and thus the only way to get a precise estimate of all included 2 nd order stress tensors is to simulate this history with sufficient precision, which is, as already mentioned, an extremely difficult task, still out of reach for many complex processes, such as e.g. welding; however, as these 2 nd order stresses arise from fluctuations ‡ We distinguish here the calculated (calc.) quantities from the modelled (mod.) ones: the first ones refer to the predicted quantities using the SC model; the second ones refer the approximated expression we will fit on the measurement points. This expression can be more or less simple and may also need the use of the SC model, as seen below.

around average values, we assume that they can be neglected (1 st assumption) and concentrate on the determination of the 1 st order stress tensor. In fact, in some studies, the second term is even completely omitted, without any justification [START_REF] Ligot | Stress analysis of polycrystalline thin films and surface regions by X-ray diffraction[END_REF]. Eq. ( 25) can then be re-written as:

〈• 'ƒ 〉 |"•ž. = 〈$′ ''˜™ ! : %′ ˜™š› ! 〉 | : ‰ z ‰ {W W "•ž. = $′ hh > aŸ ‰ z ‰ {W W "•ž. (26)
In this expression, the matrix A depends solely on the measurement direction but the tensor $′ aŸ depends in quite a complex way on the single crystal elastic constants as well as on the diffraction volume (and thus the texture of the material). If we want to avoid requiring to a model to assess it, we can make the following assumptions:

(ii) In the general case, the stress tensor has 6 independent components, but it often reduces to 3, since it is usually considered that the stress state measured on the surface of the material is a plane stress state, because of the free surface (2 nd assumption);

(iii) If the material can be considered to be elastically isotropic (3 rd assumption), the $′ aŸ tensor does not depend anymore on the diffracting plane and volume, but solely on the macroscopic elastic constants which are known a priori (called 1¡ and f¡ in Eq. ( 27) below, defined from ¢ and £, the macroscopic Young Modulus and Poisson ratio) [START_REF] Afnor | Non-destructive Testing -Test Method for Residual Stress analysis by X-ray Diffraction[END_REF]. In this case, Eq. ( 26) reduces to the well-known sin f x linear fit (because 〈• 'ƒ 〉 | becomes linear as a function of sin f x, see Eqs ( 27) and ( 28)) and only one peak and three y values (to extract the 3 stress components) are needed. For hexagonal materials, the 123 3 peak is generally recommended, because of the high multiplicity of this crystallographic plane family.

〈• 'ƒ 〉 |"•ž. = 1 + £ ¢ 2 11 "•ž. cos f y + 1f "•ž. sin 2y + ff "•ž. sin f y4 sin f x - £ ¢ 2 11 "•ž. + ff "•ž. 4 = 1 2 f¡ . § '"•ž. . sin f x + 1¡ Tr § "•ž.
with 1¡ = -£ ¢ © and

1 f f¡ = 1ª« * § '"•ž. = 11 "•ž. cos f y + 1f "•ž. sin 2y + ff "•ž. sin f y (27)
(iv) In any other case, the $′ aŸ tensor has to be calculated with a polycrystalline model describing the elastic response of the material, like e.g. the one used above to simulate the tensile curve [START_REF] Pang | In situ generation of intergranular strains in an Al7050 alloy[END_REF][START_REF] Kroner | Berechnung der Elastischen Konstanten des Vielkristalls aus den Konstanten des Einkristalls[END_REF], and the modelled average strain then reads

〈• 'ƒ 〉 |"•ž. = 1 2 $′ |1 aŸ . § '"•ž. . sin f x + $′ |f aŸ Tr § "•ž. ( 28
)
The terms $′ |1 aŸ and $′ |1 aŸ -called the Diffraction Elastic Constants (DEC) -which now depend on the diffracting plane as well as on diffracting volume Ω (relative to the texture of the material) replace the macroscopic quantities, and can also be calculated with the SC model in it elastic version. It is worth mentioning though that, if the SC model has proven its excellent predictive capacity in elasticity, some other models, such as Voigt, Reuss (or the average of the two previous ones, also called Hill estimation) can also be selected since they provide explicit relations for the macroscopic Hooke's law and are thus much simpler to use while still providing an acceptable accuracy [START_REF] Belkhabbaz | Prediction of the overall behavior of a 3D microstructure of austenitic steel by using FFT numerical scheme[END_REF].

The 3 macrostress components are then estimated by the least square method, i.e. by minimizing the following expression, in which the sum is generally performed -for one selected diffraction peak -on all measured points (i), associated with varying angles (x, y).

∑ ¬〈• 'ƒ 〉 |"… †. -〈• 'ƒ 〉 |"•ž. - f = Min (29)
In this equation, 〈• 'ƒ 〉 |"•ž. is represented either by Eq. 27 if we assume elastic isotropy or by Eq. 28 if we do not assume elastic isotropy. In this case, the material can further be untextured (i.e. macroscopically isotropic) or textured. In the first case, it is easy to show that the average elastic strain is again linear with respect to sin 2 x, but with a slope which varies with the diffraction plane.

In any other case, the variation generally becomes non linear. In most cases, the minimization is realized on a series of measurements performed for one single peak, but some recent procedures rely on multiple peak minimization [START_REF] Dufrenoy | Modeling methodology for stress determination by XRD in polycrystalline materials[END_REF][START_REF] Daymond | Measured and predicted intergranular strains in textured austenitic steel[END_REF]. It has also been suggested by the French standard that the DEC can also be experimentally determined if the material is single-phased using 4-points bending test with XRD, but it is rarely done in practice [START_REF] Afnor | Non-destructive Testing -Test Method for Residual Stress analysis by X-ray Diffraction[END_REF][START_REF] Jeong | Evaluation of biaxial flow stress based on elasto-viscoplastic self-consistent analysis of X-ray diffraction measurements[END_REF].

For industrial applications, the material is most often considered elastically isotropic, which suppresses the necessity to go through a complex calculation. We now consider our predicted data as a virtual experimental data set and evaluate the impact of the proposed simplifying assumptions on the precision obtained on the quantitative estimate of the sole 1 st order macroscopic stress.

Results

For each of the 4 selected points along the predicted tensile curve, macroscopic elastic strains 〈• 'ƒ 〉 |""•. have been calculated with the elasto-plastic model for 2 diffracting peaks 123 3 and {0006} (the 1 st one associated with a high multiplicity and the 2 nd with a very low multiplicity) and several values of x and 2 values of y. For these calculations, both elastic and textural anisotropies have been taken into account and the averaging have been made on all grains in diffracting position, with an additional spread of 5° around them. These calculated elastic strains are represented in Fig. 6 to 9 and the associated diffracting volumes (evaluated at 0 strain) are shown in Fig. 10. For points 2 and 4, for which the elastic strains contain two components (see Eq. ( 24b)), these two components have been also represented separately.

These figures call for the following comments:

(i) It is interesting to note, that for both peaks, the predicted evaluations of the 1 st order strain component 〈• 'ƒ 〉 |""•. appear quasi-linear for the 4 points, in spite of the fact that the material is both (weakly) elastically and (strongly) texturally anisotropic. This is presumably because the elastic anisotropy is small.

(ii) For the {0006} measurements, some points are "missing" in the figures on the curves corresponding to y = 0, since the influence of the texture is quite strong in this case, and for some y values, there is no grain in diffracting position (see also Fig. 10). It is worth mentioning though that this situation is very often encountered in stress measurements, since the most studied materials in this context (steels and titanium alloys) are generally strongly textured (see e.g. the experimental diffracting volume fractions in [START_REF] Jeong | Evaluation of biaxial flow stress based on elasto-viscoplastic self-consistent analysis of X-ray diffraction measurements[END_REF], which are on the same order of magnitude as the ones calculated here). Due to the discrete description of the texture in the present case, the evaluation of the diffracting volume depends strongly on the selected spread (taken here equal to 5°, which is quite small), but a diffracting volume equal to 1% still corresponds here to 500 orientations.

(iii) As soon as 2 nd order stresses are created within the material, some non linearity appears clearly on both total strain 〈• 'ƒ 〉 |""•. and 2 nd order strain 〈• 'ƒ 〉 |""•. profiles (see Points 2 to 4).

(iv) For point 4 -which is similar to a case where 1 st order stresses are classically determinedthe total strain profile deviates significantly from linearity for both peaks, and especially for y = 0.

(v) We can also estimate from these curves, the quantitative influence of the 2 nd order term, whenever it is present; it is seen on all concerned curves that 〈• 'ƒ 〉 |""•. and 〈• 'ƒ 〉 |""•. are of the same order of magnitude, although not always of the same sign; generally, the 2 nd order term is smaller than the 1 st order one, but it can sometimes be larger (for example, for point 4, at y = x = 0 for the 123 3 peak, we have 〈• 'ƒ 〉 |""•. = 0.001 and 〈• 'ƒ 〉 |""•. = -0.0005).

The total strains calculated for the 123 3 peak at y = 0 are now replotted in Fig. 11 for the 4 points. The influence of the imposed macroscopic stress on the maximal strain value is quite clear: this maximal value is close to 0.002 for an imposed stress equal to 200 MPa in the elastic range (after the establishment of plastic residual stresses (Point 4) or not (Point 1)), whereas it reaches 0.004 for an imposed stress close to 500 MPa in the elasto-plastic range, and is reduced to 0.001 in the unloaded condition (which would correspond to the case of a macroscopic residual stress equal to 0). Now, considering that these calculated elastic strains are experimental data, we can estimate the 1 st order stresses for the 4 studied points through the minimization procedure explained above by neglecting the 2 nd order term (1 st assumption) and by assuming that the material is either (i) elastically isotropic (EI) (3rd assumption, Eq. 27), (ii) elastically anisotropic but macroscopically isotropic (IT) (Eq. 28) or (iii) elastically and macroscopically anisotropic (FA) (again Eq. [START_REF] Dufrenoy | Modeling methodology for stress determination by XRD in polycrystalline materials[END_REF]. Two examples of this adjustment are shown in Fig. 12 and all resulting stress determinations are presented in Table 2.The stress state is also assumed to be plane stress, which is completely consistent with the imposed boundary conditions (2 nd assumption). The minimization has been performed on points associated with one single peak (as in recommended standards) -here the 123 3 peak or the {0006} peak, or with two peaks simultaneously (proposed method) (see Eq. 29).

The data presented in Fig. 12 correspond in a sense to two extreme cases. The 1 st one (Fig. 11 left) concerns the estimation of a sole 1 st order stress tensor in the absence of 2 nd order stresses. In this case, the adjustment is necessarily better than in all other cases, since the neglected term is 0.

Also, the influence of elastic or textural anisotropy is quite limited, since all curves are linear and the calculated error parameter is indeed quite small for all performed adjustments. The 2 nd case presented in Fig. 12 (right) corresponds to the worst case for stress estimation, since both terms are present in the predicted value. It is seen that the 4 calculated profiles are non linear in this case (see also Fig. 9). The adjusted profiles look however linear, in spite of the fact that the texture is taken into account. This is again presumably due to the weak elastic anisotropy. Now, from the close inspection of the estimated stress values for each point and associated error parameter (Table 2), we can say that:

(i) For Point 1, for which the only origin of stress is elasticity, the adjustment is quite good for all cases, but always slightly better by taking into account both elastic and plastic anisotropy; the minimization on two peaks is also better than the minimization on one single peak;

(ii) For Point 2, for which the level of stress is quite high, again taking into account both elastic and textural anisotropy improves the quality of the fits, except for the {0006} one for which it has been seen that the diffracting volume was reduced for a lot of measurement points; minimizing on two peaks however still provides a better estimation than the minimization on one single peak;

(iii) Point 3 corresponds to the situation in which the pseudo-experimental elastic strain is composed on the sole term which is neglected in the minimization procedure. Whatever the selected peak, the error appears smaller when the elastic (and thus also textural) anisotropy is neglected.

(iv) For Point 4, taking into account the anisotropy of the material does not improve the quality of the fit, when the minimization is performed on one single peak. However, when the fit is performed simultaneously on two peaks, the quality of the fit is significantly improved, and especially when all sources of anisotropy are taken into account. In such a case, the error parameter is indeed divided by 1.12 when going from the recommended standard (1 st adjustment line, last column ¯ = 124.9) to a double peak minimization (7 th adjustment line, last column ¯ = 56.4).

It is thus clear from Table 2 that the adjustment on one single peak is acceptable only when there is no 2 nd order stress of inelastic origin (point 1 only) but also that in the 3 other cases, the adjustment is considerably improved by performing the minimization on both peaks simultaneously. This is nicely illustrated graphically in Fig. 13. The proposed multipeak minimization performed by taking into account both elastic anisotropy and texture (results highlighted in pink in Table 2) provides a much better estimation of the macroscopic stresses than the procedure recommended by actual standards (results highlighted in grey in Table 2) for the situations where 1 st and 2 nd order stresses are present. This is a very interesting result, since most of the cases for which stresses are estimated in practice correspond to cases for which both 1 st and 2 nd order terms (of various origins and not only of plastic origin) are present. It is also worth noting that although the calculated curves may appear quite close for Points 1 and 4 (see Fig. 11), the corresponding adjusted stress values resulting from a single peak minimization can be quite different (see Table 2).

Conclusions and perspectives

We have detailed in the present paper the classical procedure for the 1 st order stress estimation from XRD and calculated with a model the 1 st and 2 nd order stresses and associated elastic strains arising from an elasto-plastic deformation in a Ti tensile sample. The purpose of this calculation was not to estimate the predictive capacity of the model or to compare experimental and calculated data as in previous studies [START_REF] Jeong | Evaluation of biaxial flow stress based on elasto-viscoplastic self-consistent analysis of X-ray diffraction measurements[END_REF][START_REF] Baczmański | Elastoplastic properties of duplex steel determined using neutron diffraction and self-consistent model[END_REF], but simply to produce a virtual experimental data set. Thus, by considering the calculated data as pseudo-experimental data, we were able to apply on them the method recommended by all standards for stress determination, and to compare the results to a method based on multi-peak minimization.

For the tested case, we have shown that -Neglecting the 2 nd order stresses in the determination of the 1 st order ones induces a significant error when the standard method is employed (i.e. minimization on one single peak); this point has already been underlined by several authors, which compared in the past experimental and calculated profiles [START_REF] Jeong | Evaluation of biaxial flow stress based on elasto-viscoplastic self-consistent analysis of X-ray diffraction measurements[END_REF][START_REF] Baczmański | Elastoplastic properties of duplex steel determined using neutron diffraction and self-consistent model[END_REF]. However, these authors did not quantify the error made by completely neglecting the 2 nd order term as in the present case (aimed at evaluating the capacity of the proposed minimization procedure to estimate the 1 st order stress without being obliged to calculate the 2 nd order term), and introduced rather an additional adjusting parameter to minimize the error made between experimental and calculated profiles [START_REF] Baczmański | Examination of the residual stress field in plastically deformed polycrystalline material[END_REF], and thus to evaluate the predicting capacity of one model in simple cases. The use of an additional adjusting parameter has also been employed by other authors more recently [START_REF] Gloaguen | Measurement and prediction of residual stresses and crystallographic texture development in rolled Zircaloy-4 plates: X-ray diffraction and the self-consistent model[END_REF][START_REF] Gloaguen | Residual stress fields analysis in rolled Zircaloy-4 plates: Grazing incidence diffraction and elastoplastic self-consistent model[END_REF].

-This error is considerably reduced when minimization is performed on two peaks, although the reason for this improvement is not completely clear yet;

-When minimization is performed on two peaks, taking into account elastic anisotropy (IT or FA hypotheses) and texture only very slightly improves the result compared to the IE hypothesis.

As the presented data obviously depend on the selected material (texture) and sampling procedure, this work is now continued to test the validity of the multi-peak minimization in other cases, that is for different texture and different stress origins (and especially with residual stresses due to thermal origin or phase transformation), in single and dual phase Ti alloys. Also, some experimental data obtained with the new portable device [START_REF] Dufrenoy | Influence of surface integrity on fatigue limit: Ti-10V-2Fe-3Al titanium application[END_REF][START_REF] Dufrenoy | Modeling methodology for stress determination by XRD in polycrystalline materials[END_REF] are presently analyzed to confirm the validity of a multi-peak analysis. Indeed, with this device, it is possible to measure elastic strains during a tensile test (in order to know a priori the value of the 1 st order macroscopic stress, just like in the present work), Also, the analysis of the data, which is performed with the software MAUD § is based on

  The previous elastic stiffness " and compliance $ tensors are then replaced by tangent G and D modulus tensors. Again we can write macroscopic laws & ; = G -: ( ; ( ; = D , : & ;

2 (

 2 elasto-plastic range): ≠ 0 but no simple relationship exists between global and local stresses, since the final stress state depends on the whole strain path; from the definitions given in Eq. 1, we can write= + *s ,where *s are intergranular fluctuations around the average value due to the elasto-plastic strain. This can be re-

  #. -!!%&'. ( ) + )) "#. -))%&'. ( ) + !) "#. -!)%&'. ( ) .

Fig. 1 :

 1 Fig. 1: Definition of the investigated V and diffracting Ω volumes. Illustration of the 1 st and 2 nd order stresses in an unloaded sample.

Fig. 2 :

 2 Fig. 2: Pole figures measured by X-ray Diffraction on a rolled and annealed titanium sample. RD and TD stand for rolling and transverse directions respectively. The intensities are multiples of the random intensity.

Fig. 3 .

 3 Fig. 3. Calculated Young's modulus for a Ti HCP crystal expressed in the crystal frame (left) and for the T40 sample expressed in the macroscopic frame (right). The values are expressed in GPa.

Fig. 4 .

 4 Fig. 4. Simulated tensile curve for a textured titanium sample (left) and Evolution of the transverse and thickness strains (the dashed line corresponds to y = x) (right).

Fig. 5 :

 5 Fig. 5: Definition of the two reference frames, linked to the sample (X1, X2, X3) and to the diffracting vector, characterized by the two angles and and the axes (L1, L2, L3).

Fig. 6 :Fig. 7 :Fig. 8 :Fig. 9 :Fig. 10 :

 678910 Fig. 6: Calculated elastic strains 〈 〉 . as a function of sin for Point 1 (elastic loading) on the tensile curve; the macroscopic elastic strain depends only on the macroscopic 1 st order stress tensor:

Fig. 11 :Fig. 12 :Fig. 13 :

 111213 Fig. 11: Calculated elastic strains 〈

Table 2 :

 2 Comparison of calculated and adjusted principal stresses for the 4 investigated points and 3 different assumptions concerning the material. The least square procedure has been applied on 123 3 and 0006 peaks, first separately and then together. The calculated error parameter χ is equal to:

Table 1: Selected slip systems for titanium and associated hardening parameters.

Acknowledgments: Fruitful discussions about residual stresses with D. Aliaga and N. Guillemot from Airbus Helicopters are acknowledged.

Funding:

The authors would like to thank the Conseil Général de Seine St Denis (CG93) for the attribution of a PhD fellowship to S. Dufrenoy.