C Villani

M Maurice

G Levesque

Kermouche

S Breumier

A Villani

C Maurice

M Lévesque

G Kermouche

Effect of crystal orientation on indentation-induced residual stress field: Simulation and experimental validation

Keywords: High Angular Resolution Electron Back-Scattered Diffraction, Crystal Plasticity Finite Element, Single crystal, Shot-peening

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

Residual surface compressive stresses are known to significantly increase metallic alloys fatigue life, corrosion resistance and wear properties. Such stresses often result from plastic strains induced by impact mechanical surface treatments. For instance, shot peening consists in stretching a thin surface layer of a ductile metallic part by impacting it with hard particles called shot. Deep rolling consists in pressuring a hard roller onto a ductile metallic surface.

Like shot peening, the process stretches a thin surface layer in the treated part. Turning, milling and drilling can also induce a near-surface residual stress-field through tool-part interactions. Non contact surface treatments, like laser peening -for which a laser irradiationinduced plasma transmits a pressure shockwave to the treated part -also induce compressive 1 residual stresses. Predicting the residual stress field resulting from a given manufacturing process is challenging since numerous parameters must be accounted for: process parameters, material behavior, geometry, environment etc.

Recent shot peening simulations can account for shot-shot and shot-target interactions, which allows for simulating the process from the shot exiting the nozzle up to their effect on the treated part's residual stresses. [START_REF] Gallitelli | Simulation of shot peening: from process parameters to residual stress fields in a structure[END_REF][START_REF] Ghasemi | Shot peening coverage effect on residual stress profile by FE random impact analysis[END_REF][START_REF] Hu | A unifying approach in simulating the shot peening process using a 3d random representative volume finite element model[END_REF][START_REF] Jiabin | An evaluation on SP surface property by means of combined FEM-DEM shot dynamics simulation[END_REF][START_REF] Tu | A sequential dem-fem coupling method for shot peening simulation[END_REF]. Most simulation works considered macroscopic and isotropic constitutive theories, which implicitly assume that the grain size is much smaller than the area between a single shot and the target surface. Simulations can usually predict residual stress fields within reasonable accuracy because methods like X-Ray Diffraction and Hole drilling provide spatially averaged residual stresses over many grains. [START_REF] Ghasemi | Shot peening coverage effect on residual stress profile by FE random impact analysis[END_REF][START_REF] Ikushima | Numerical analysis of residual stress distribution on peening process[END_REF][START_REF] Mahmoudi | A comprehensive experimental and numerical study on redistribution of residual stresses by shot peening[END_REF].

Very few authors have attempted to simulate shot impacts on models where grains are explicitely represented and where crystal plasticity constitutive laws are used [START_REF] Moreno | Shot peening fem simulation: A novel approach based on crystal plasticity[END_REF][START_REF] Musinski | On the eigenstrain application of shot-peened residual stresses within a crystal plasticity framework: Application to ni-base superalloy specimens[END_REF][START_REF] Rousseau | Finite element multi-impact simulations using a crystal plasticity law based on dislocation dynamics[END_REF]. The prohibitive computational time required to run such simulations and the challenges associated with measuring intragranular residual stresses might explain why such fundamental work has seldom been attempted. Kobayashi et al. [START_REF] Kobayashi | Mechanism of creation of compressive residual stress by shot peening[END_REF] used X-Ray diffraction measurements to show that the in-depth residual stress profiles produced on steel polycrystals by impact and quasi-static indentation were different. However, the spherical impacter and indenter they used had diameters of 50 mm and 75 mm, which is much larger than typical shot whose diameters can range from 0.1 mm to 2 mm [START_REF] Nordin | Measuring shot peening media velocity by indent size comparison[END_REF]. Juran et al. [START_REF] Juran | Investigation of indentation-, impact-and scratch-induced mechanically affected zones in a copper single crystal[END_REF] used Electron BackScattered Diffraction (EBSD) measurements to demonstrate that the orientation gradients produced on a single crystal by a 2.5 mm diameter ball impact at 100 m s -1 and spherical indentation are similar, for similar imposed loads. The strain rate involved in their study was about 100 s -1 , which is lower than usual shot peening strain rates, ranging from 10 4 s -1 to 10 6 s -1 [START_REF] Meguid | 3D finite element analysis of peening of strain-rate sensitive materials using multiple impingement model[END_REF]. However, these results suggest that the local residual stress field induced by spherical indentation could be used to estimate that induced by a shot impact.

Single crystals plastic flow resulting from indentation has been widely investigated in the past twenty years. Most researches investigated the so-called size-effect in metals [START_REF] Demir | Investigation of the indentation size effect through the measurement of the geometrically necessary dislocations beneath small indents of different depths using ebsd tomography[END_REF][START_REF] Rester | Indentation across size scales -a survey of indentation-induced plastic zones in copper {111} single crystals[END_REF], or developed inverse methods to identify crystal plasticity parameters [START_REF] Alcalá | Extracting uniaxial responses of single crystals from sharp and spherical hardness measurements[END_REF][START_REF] Petryk | On direct estimation of hardening exponent in crystal plasticity from the spherical indentation test[END_REF][START_REF] Renner | Sensitivity of the residual topography to single crystal plasticity parameters in berkovich nanoindentation on fcc nickel[END_REF]. Other works also focused on the single crystal indentation-induced strain mechanisms by analysing either the pile-up/sinking-in patterns [START_REF] Alcalá | The influence of plastic hardening on surface deformation modes around vickers and spherical indents[END_REF][START_REF] Casals | Finite element crystal plasticity analysis of spherical indentation in bulk single crystals and coatings[END_REF][START_REF] Liu | Indentation analysis of mechanical behaviour of torsion-processed single-crystal copper by crystal plasticity finite-element method modelling[END_REF] or EBSD-measured crystal misorientation fields [START_REF] Zaafarani | Three-dimensional investigation of the texture and microstructure below a nanoindent in a cu single crystal using 3d ebsd and crystal plasticity finite element simulations[END_REF].

In particular, the orientation-dependent behaviour of such strain mechanisms has been widely studied [START_REF] Han | Indentation Pileup Behavior of Ti-6Al-4V Alloy: Experiments and Nonlocal Crystal Plasticity Finite Element Simulations[END_REF][START_REF] Liu | Indentation analysis of mechanical behaviour of torsion-processed single-crystal copper by crystal plasticity finite-element method modelling[END_REF]. Similar studies on bicrystal have also been conducted to study the particular mechanisms induced near grain boundaries [START_REF] Dahlberg | Geometrically necessary dislocation density measurements at a grain boundary due to wedge indentation into an aluminum bicrystal[END_REF][START_REF] Liu | Crystal plasticity fem study of nanoindentation behaviors of cu bicrystals and cu-al bicrystals[END_REF][START_REF] Su | Quantifying deformation processes near grain boundaries in alpha titanium using nanoindentation and crystal plasticity modeling[END_REF].

Fewer works focused on the residual stress field induced by indentation. Zheng et al [START_REF] Zeng | Residual stress fields at the surface of sharp pyramid indentations[END_REF] proposed an analytical model to predict the residual stress field produced by pyramidal indentation on silica. Boyce et al [START_REF] Boyce | The residual stress state due to a spherical hard-body impact[END_REF] evaluated the residual stress field produced by a spherical impact on a Ti-6Al-4V polycrystal using X-ray diffraction and compared their experimental results to finite element simulations reproducing the same experiment. However, to the best This paper investigates the residual stress field generated by spherical indentations on single crystals. The study builds on the experimental work of Juran et al. [START_REF] Juran | Investigation of indentation-, impact-and scratch-induced mechanically affected zones in a copper single crystal[END_REF], more specifically on a <001> indented copper single crystal. Residual stress fields quantified by HAR-EBSD are succesfully compared to those predicted by CPFE simulations. CPFE simulations are further used to investigate the crystal orientation and repeated indentations effects on the subsurface residual stresses.

The paper is organized as follows: section 2 describes the main theories and experimental methods we relied on. Section 3 presents the Finite Element model used. Section 4 details the numerical model experimental validation, predicted indentation results for single indentations on different orientations and repeated indentations on the same crystal. Finally section 5 discusses the results and concludes this paper.

Background

Spherical indentation on a copper single crystal

The experimental data used in this work was extracted from the spherical indentation test presented in the work of Juran et al. [START_REF] Juran | Investigation of indentation-, impact-and scratch-induced mechanically affected zones in a copper single crystal[END_REF] and is briefly recalled here for the sake of completeness. The copper single crystal was produced by directional solidification using a standard zone melting method based on a horizontal Bridgman-type apparatus. The samples were cut by Electric Discharge Machining (EDM) into parallelepiped-shaped specimens and mechanically polished. A final electro-polishing step removed remnants of mechanical hardening due to prior polishing. The indented surface was {100} oriented along the surface normal. Quasi-static indents were performed using a load-controlled standard Brinell

High Angular Resolution EBSD measurements

Residual stresses are produced by local lattice elastic distortions that induce subtle changes in the EBSD diffraction pattern. The HAR-EBSD analysis retrieves local elastic strain variations by measuring and suitably interpreting the warping field between a well-chosen reference pattern and any other pattern in the map. Only relative deviatoric strains can be measured, but assuming a plane stress condition, together with Hooke's law, leads to the full strain and stress tensor fields recovery. Many improvements to the original method [START_REF] Wilkinson | High-resolution elastic strain measurement from electron backscatter diffraction patterns: New levels of sensitivity[END_REF] have been made in recent years [START_REF] Maurice | On solving the orientation gradient dependency of high angular resolution ebsd[END_REF][START_REF] Villert | Accuracy assessment of elastic strain measurement by ebsd[END_REF]. The analysis presented in this work was performed using an EBSD pattern iterative remapping algorithm and finite strain theory, as explained in Maurice et al. [START_REF] Maurice | On solving the orientation gradient dependency of high angular resolution ebsd[END_REF].

In-depth EBSD measurements were performed using a Zeiss supra 55-VP Field Emission Gun (FEG) Scanning Electron Microscope (SEM) operated at 20 kV with a probe current of 2 nA. The EBSD analysis was carried-out using an HKL system (Oxford instrument) composed of NordlysII camera and the Channel 5 software suite.

Crystal plasticity framework

The single crystal constitutive behavior was modeled using the Meric-Cailletaud largestrain formulation crystal plasticity framework [START_REF] Méric | Finite element calculations of copper bicrystal specimens submitted to tension-compression tests[END_REF]. The deformation gradient F was decomposed between its elastic and plastic contributions, respectively F e and F p , as:

F = F e • F p , ( 1
)
where • denotes the singly contracted product. The plastic part was related to the slips occurring in the different slip systems through

Ḟp • (F p ) -1 = N s=1 γs m s ⊗ n s , ( 2
)
where γs is the s-th slip system shear strain rate, N is the number of activable slip systems ( 12for face-centered cubic structures) and m s and n s are unit vectors representing respectively the slip direction and the normal to the slip plane.

The stress τ s projected on each slip system was related to the system's shear strain rate according to a Norton law where

γs = f s (τ s ) K n , (3)
with

< x >= 0 if x < 0 x if x ≥ 0 ,
K and n being material constants. An isotropic hardening term r s was added to the flow rule f s as

f s = |τ s | -r s . ( 4
)
r s results from the interactions between the different dislocation slip systems modeled by

r s = R 0 + N q=1 h sq (R q -R 0 ), (5)
with

R q = R 0 + Q(1 -e -v q b ), (6)
where R 0 is the critical resolved shear stress, Q and b are two phenomenological constants, v q is the cumulated plastic slip for the q-th system and h sq is the interaction coefficient between systems (s) and (q). These interaction coefficients represent the averaged dislocation interactions contribution between systems (s) and (q). For Face-Centered Cubic (FCC) structures, the interaction matrix contains six different coefficients that account for the interaction phenomena listed in Table 1.

The coefficients used in the finite element analysis for copper are those identified in [START_REF] Méric | Finite element calculations of copper bicrystal specimens submitted to tension-compression tests[END_REF] and recalled in Tables 1, 2 and 3.

Finite element analyses

K (MPa s 1/n ) n R 0 (MPa) Q (MPa) b 5 10 10 6

15

Table 3: Norton law and isotropic hardening parameters [START_REF] Casals | Finite element crystal plasticity analysis of spherical indentation in bulk single crystals and coatings[END_REF] coupled to ABAQUS to integrate the crystal plasticity constitutive law described in section 2.3.

The copper sample was modeled as a parallelepiped divided in two regions, as shown in Figure 1: a refined zone (Zone A) in the contact area and a border zone (Zone B) with progressively unrefined meshes to avoid border effects. Dimension R (Zone A) was chosen so as to contain all the elements where plastic strains occurred. The outside layer thickness b was chosen to avoid border effects on the residual stress field (see Appendix A). The Dimension b = R was chosen as the converged value.

The substrate was meshed with 3D hexahedric linear reduced integration elements (C3D8R) using an updated lagrangian formulation. Reduced integration was used to avoid locking effects resulting from plastic incompressibility [START_REF] Giannakopoulos | Analysis of vickers indentation[END_REF]. Very little mesh instabilities (e.g. hourglass effect) were observed and were corrected using a low hourglass stiffness (1 MPa). The mesh size was determined through a convergence study based on the residual stress field (see Appendix A). The converged element size was taken as 28 µm in the contact zone. The chosen mesh density also ensured that the maximum contact radius was composed of at least 10 elements.

The indenter was modeled as a rigid half-sphere meshed with tetrahedral elements. A penalty algorithm was used for the contact, which was considered frictionless. Simulations with different friction coefficients were also performed and led to similar residual stress fields (see details in Appendix B).

The simulations were performed in three steps. (i) First, the indenter was moved (displacement- control) toward the substrate along the y-axis down to a depth of 40 µm. This corresponds to the experimental maximum indentation depth. (ii) The indenter was gradually removed from the contact surface to simulate unloading, following the y-axis. (iii) Finally, half the substrate's elements (with coordinates z>0 mm) were deactivated to model the sample's cutting. All displacements were blocked at the bottom of the parallelepiped through all the steps (U x = U y = U z = 0). This reproduces the experimental boundary conditions, as the sample was stuck on a sample holder before indentation.

F (N) Experimental R 0 =10 MPa R 0 =15 MPa + + + + + + + R 0 =20 MPa × × × × × × × R 0 =35 MPa * * * * * * *
To recover the experimentally measured macroscopic behavior, the choice was made to vary only the critical resolved shear stress parameter R 0 . This choice allows to retrieve the experimentally imposed load for a given indentation depth, without altering the simulated material's hardening behavior. The indent size as a function of the indentation force was plotted for several R 0 and compared to the experimental values obtained for the numerous indentations reported by Juran et al. The indent diameter is defined as the largest dent dimension along the x-axis. Figure 2 shows that R 0 = 10 MPa best fits the experimental data.

Simulation of a single indent aligned with the [001] crystal direction was performed for comparison with experimental data. The crystal [100] direction was aligned with the simula- The stress redistribution during successive indents was investigated. Five successive indents were simulated following the indentation strategy presented in Figure 3, in the same simulation run. The fives indents were performed along the y-axis at x = 0 and for different z values. The first three indents' centers were spaced apart by a distance of d imp . This distance corresponds to a single imprint diameter. The order chosen for these three indents is specified by the numbers in Figure 3. The two last indents were performed on the resulting dimples.

Their centers were spaced by a distance of

d imp 2
from the first indent center. The geometry was expanded to contain the plastic strains in the refined mesh zone.

Results

Numerical predictions comparison with experimental data

Figure 4 shows the simulated indentation imprint highlighted by the contact pressure field at maximum indentation depth. This imprint is compared to that experimentally measured.

Both experimental and numerical imprints are square-shaped, as in the work of Juran et al, and have similar dimensions. It can however be noted that the experimental imprints present a rounder shape than the predicted one.

Figure 5 shows the experimentally measured and predicted σ xx residual stress field over the specimen's cross-section after cutting. Please note that a zero-centered color scale was chosen for all the residual stress fields presented in this study to better emphasize the relative importance of tensile and compressive stresses. Figure 5 shows that the HAR-EBSD method can capture local stress variations within the crystal. below the surface (see Section 5).

× ××××× × × × × × ×× × × × × × × × × × × × × × × × × × × × × × × × ×
A large tensile zone is present below the indent, for both the simulation and the experiment. These observations suggest that our simulation methodology is sufficiently robust to further investigate the crystal's orientation effect on the residual stress field during an indentation test, at least qualitatively.

Influence of crystal orientation on the residual stress field resulting from a single indent

Figure 7 shows the predicted residual stress field in an indented sample before cutting.

This figure further confirms that spherical indentation in the crystal's [001] direction produces a large tensile residual stress zone beneath the indent.

The volume fraction f t of elements for which the first stress invariant was positive in the refined mesh zone within the first 200 µm under the surface was computed as:

f t = Ne i=0 A i V i V tot , ( 7
)
where

A i = 1 if tr(σ) 3 ≥ t 0 if tr(σ) 3 < t . ( 8
)
N e is the number of elements in the first 200 µm below the surface, V i the i-th element when considering the volume fraction of elements with a negative stress invariant, f c , using a threshold t =-20 MPa. Appendix C details the methodology used to generate these figures.

Please note that the figures do not complement each other as a stress threshold t was used to calculate both volume fractions.

These figures show that many crystal orientations can lead to tensile zones after one indent. The largest tensile zones were found for samples indented close to the [001] orientation while the smallest tensile zones were found for those indented close to the [101] orientation.

However, Figure 8(b) shows that these orientations do not yield the largest compressive zones, which are found close to the [213] orientation.

Further simulations were performed on fourteen orientations ranging from [001] to [111] in the (1 10) plane to visualize the residual stress field evolution leading to the f t and f c values presented in Figure 8.

Figure 9 shows the residual first stress invariant field evolution as a function of crystal orientation. The tensile residual stress zone observed for the sample indented along the [001] direction is confined, almost uniformly, close to the surface, when the indentation direction shifts toward the [111] direction. This is particularly emphazised in Figure 10 that reports the in-depth stress profile taken at the indentation center for three different orientations. For the sake of clarity, the curves presented in Figure 10 were smoothened to remove discretizationinduced oscillations using a gaussian filter as:

Γ f iltered (x) = ∞ -∞ Γ(x -t) 1 s √ 2 e -(t-µ) 2 2s 2 dt, (9)
where Γ f iltered and Γ are respectively the filtered and the raw curves, and s and µ are the gaussian parameters (taken as µ = 0 and s = 4).

Stress redistribution after several indentations

The residual stress fields produced after each of the five indents are presented in Figures 11 and12 for samples indented in the [001] and [111] direction, respectively.

A tensile residual stress zone is still present underneath the indented surface for the

[001] oriented sample after five indents. By opposition, the residual stress field predicted after indenting in the [111] direction is purely compressive, and more compressive than that predicted for the [001] orientation. The depth where compressive residual stresses can be found is also deeper for the [111] oriented sample.

Discussion

Figure 6 shows that more compressive residual stresses are measured close to the indentation axis (x = 0 mm) and beneath the tensile zone deeper in the sample. Such bias in the compressive stress intensities could be induced by the HAR-EBSD parameter choice. Errors in the HAR-EBSD pattern center coordinates, or in the reference pattern stress, could for instance lead to such bias [START_REF] Maurice | On solving the orientation gradient dependency of high angular resolution ebsd[END_REF][START_REF] Villert | Accuracy assessment of elastic strain measurement by ebsd[END_REF]. Also, small misalignments of the cutting plane with the indentation center could also contribute to such bias. This could also result from the simplicity of the procedure used to simulate sample's cutting: the numerical procedure assumes that no residual stresses have been generated through the cutting process.

Our model predicted a significantly spread tensile residual stress zone after spherical indentation, whose size depends on the crystal orientation. A tensile zone after dynamic impact has been reported before in the literature for large diameter balls by Kobayashi et al. [START_REF] Kobayashi | Mechanism of creation of compressive residual stress by shot peening[END_REF]. As recalled in Section 1, the authors performed quasi-static spherical indents and impacts on a steel polycrystal. Only impact tests led to subsurface tensile residual stresses in their study. According to the authors, the tensile residual stress zone resulted from plasticity mechanisms brought by stress waves propagation and high strain rates. In our case, these residual stresses result from the anisotropic nature of crystal plasticity, as it has not yet been revealed using homogeneous macroscopic isotropic constitutive models.

Figures 11 and12 showed that the volume of tensile residual stresses decreases with an increasing number of indents. However, these figures also suggest that there is a link between the relative volume of tensile and compressive residual stresses produced after one and several indents. The larger tensile stresses volumes after one indent result in shallower compressive stresses and a more heterogeneous stress field after several indents. For grains whose dimensions are close to the contact radius, or for textured microstructures, the residual stress state of a treated material could therefore highly depend on its initial microstructure orientations.

Figure 8 showed that grains with orientations close to <001> have a higher probability to present subsurface tensile residual stresses. By opposition, indentation on grains with orientations close to the center of the IPF should exhibit larger compressive stresses volumes.

Knowing the initial microstructure orientations, the IPF map of f t and f c can therefore provide insights in the shot-peening induced residual stress field heterogeneities and residual compressive stress depth. This orientation-dependent behavior could therefore have consequences for surface treatments involving spherical impacts on polycrystals with a coarse grain structure or with a particular texture. This could be particularly the case for the shot peening of welded zones [START_REF] Fonda | Texture variations in an aluminum friction stir weld[END_REF][START_REF] Meo | Finite element analysis of residual stress induced by shot peening process[END_REF].

Works performed by Chen et al. [START_REF] Chen | Influence of shot peening on surfacelayer characteristics of a monocrystalline nickel-based superalloy[END_REF] showed that this effect is also present after multiple impacts on a different material. The authors shot-peened a nickel-based superalloy single crystal along the [001] and [111] orientations at 80% and 400% coverage. Their work showed that compressive residual stress intensities decreased in the [001] single crystal, when compared to that in [111] oriented single crystal, even after a 400% coverage. Their results corroborate the results we presented in this paper. Future work will focus on the study of this orientation-dependent behaviour for materials which are known to present particular textures and for dynamic impacts.

Conclusion

The objective of this work was to investigate the residual stress field generated by spherical indentations on a single crystal and to assess the relevance of CPFE for mechanical surface treatment simulations. The main contributions are as follows:

• HAR-EBSD observations revealed a large tensile residual stress zone after spherical indentation on a [001] oriented single crystal copper.

• A spherical indentation CPFE procedure using Meric-Cailletaud's crystal platicity framework was developed. A good agreement was found with HAR-EBSD observations in terms of residual stresses.

• This analysis should be conducted with one and several high velocity impacts to validate its conclusions on real engineering mechanical surface treatments. CPFE simulations for strain rates close to those encountered in shot peening will require the development of new constitutive law identification methodologies at small scale. Further work will therefore focus on high strain rate micro-mechanical testing, such as micro-pillar compression tests [START_REF] Guillonneau | Nanomechanical testing at high strain rates: New instrumentation for nanoindentation and microcompression[END_REF]. This work will be performed on copper, as it is a well-known material, on which HAR-EBSD observations can be performed with a reasonable noise level and which present no mechanical twinning nor phase transformations during shot-peening. Such material presents very few industrial applications for shot-peening but will allow to develop this methodology on a simple case, in the perspective of using it on more complicated materials in the future. [001] direction and with an angle β = 0 • (see Figure 1). Convergence found for these conditions was assumed to hold for any other crystal orientations. Convergence of f t to a value of approximately 13 % is observed for a total of 66 517 elements.

in the crystal's [001] direction and with an angle β = 0 • (see Figure 1). The value of f t stabilizes around 13 % for 66 517 elements. Evolution with the mesh density of the tensile stress volume fraction computed using the three axial stress components is shown in Figure 14. It further confirms the residual stress field convergence for each stress components. Also the stress profile taken along z=0 at a depth of 40 µm is represented in Figure 15 for different mesh densities. This results reveals that the stress variations seem only to converge around 114 518 elements. Simulations in this work were therefore performed with this mesh density.

This corresponds to 28 µm sized elements in Zone A (Figure 1). Convergence found for these conditions was assumed to hold for any other crystal orientations.

Indenter mesh

A convergence study on the indenter's mesh size was performed to remove its influence on the predicted residual stress and misorientation fields. The indenter being asymmetrically meshed (with respect to the specimen's symmetry axes), a low indenter mesh density in the contact zone leads to asymmetric fields. Convergence was therefore assessed by characterizing the residual stress field symmetry in the (x,y) plane at z=0 at several depths y using an error criterion defined as:

y = 2 N y Ny/2 i=0   σ(x i ) -σ(x N -i ) max i (σ i )   2 , ( 10
)
where N y is the number of data point at depth y.

The evolution of y with mesh density at different depths is presented in Figure 16. The figure reveals that the indenter mesh has indeed an influence on the field symmetry for low mesh densities. The indenter was meshed with 15 040 elements to reduce this influence.

Border dimensions

Border effects can induce local variations of the residual stress field, depending on Zone B (Figure 1) relative size to the contact radius. The value of f t was therefore computed Figure 17 shows the evolution of f t with Zone B to Zone A dimension ratio 2b R . Border influence on the value of f t is observed for low Zone B dimension values. Convergence can however be observed for borders twice larger than the refined zone. allows to describe all possible crystal orientations. Simulated orientations were chosen by meshing the IPF with tetrahedral elements, using gmsh software [START_REF] Geuzaine | Gmsh: A three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF], so as to cover the whole domain as homogeneously as possible. The chosen mesh is represented on Figure 19(b). Each node of the mesh represents a simulated indentation direction.

Data availability

The raw data required to reproduce these findings cannot be shared at this time as all experimental results were extracted from the work of Juran et al. [START_REF] Juran | Investigation of indentation-, impact-and scratch-induced mechanically affected zones in a copper single crystal[END_REF]. However, finite element result files, simulation input files as well as meshing tools used for this study can be found at https://doi.org/10.5281/zenodo.2550820 and can be reused.

Figure 1 :

 1 Figure 1: (a)Schematic view of the indentation model geometry and mesh. C3D8R hexahedral andC3D4 tetrahedral elements were used respectively for the copper sample and the indenter mesh. The indenter was modeled as rigid. The sample followed a CPFE law for copper. Contact was modeled as frictionless using a penalty method. (b)Representation of the substrate crystal orientation with respect to the global coordinate system. Simulations with an in-plane rotation φ = 20 • and 45 • were also performed to rule out any effect of the mesh on the residual stresses, using the same criterion as for the convergence study presented in Appendix A.

Figure 2 :

 2 Figure 2: Variation of the indent diameter as a function of the applied force for different values of the yielding criterion R 0 . The indent diameter is defined as the largest dent dimension along the x-axis. Results obtained for R 0 = 10 MPa seem to best fit experimental data, especially for the largest force values.

Figure 3 :

 3 Figure 3: Indentation strategy chosen to study the influence of several indents. d imp represents the first residual imprint diameter. Five successive indents were performed following the x-axis. First three indents were spaced apart from d imp . Two last indents were performed on the resulting dimples.

Figure 4 :

 4 Figure 4: (a) Topography of the resulting indent shape: comparison of the experimental (left) and numerical (right) residual indent enhanced by the contact pressure field. (b) Quantitative comparison of the predicted and experimentally measured imprint shape using the imprint diameters and the curvature radius at the imprint corners. Both imprint present similar dimensions. The experimental imprint is rounder than the predicted one.

Figure 5 :

 5 Figure 5: In-depth residual stress field after indentation and sample cutting at the indent center: HAR-EBSD results (left) and Finite Element results (right). Both fields present similar shapes and amplitudes. Tensile residual stresses are observed in the first 200 µm under the surface.

Figure 6

 6 presents residual σ xx profiles taken along two path located 70 µm and 493 µm below the surface for both simulation and experiment. This comparison suggests that there is a reasonable agreement between the experimentally measured and predicted residual stress fields. The measured residual stresses are more compressive under the surface at the indent center (x = 0 mm) than those predicted. Such a shift in the compressive stress values is also found deeper into the sample, at 493 µm

Figure 6 :

 6 Figure 6: Experimentally measured and predicted σ xx profiles after an indentation, along two paths 70 µm and 493 µm below the surface. The Finite Element model captures well the global shape and amplitudes of the field. The most compressive stresses are found experimentally under the surface at x = 0 mm. A deeper profile located 493 µm below the surface reveals a gap between predicted and experimentally measured stresses.

  figure for f t for 43 crystal orientations. Figure 8(b) shows the inverse pole figure obtained

Figure 7 :

 7 Figure 7: Finite Element in-depth σ xx field after indentation in the [001] direction, before sample cutting. Tensile stresses are present down to 230 µm and therefore do not result from sample's cutting.

Figure 8 :

 8 Figure 8: Inverse Pole Figure representing the volume fraction of elements with a positive (a) and negative (b) residual first stress invariant after one simulated indent. Each black mark represents a simulation result. Largest tensile zones and smallest compressive zones are found when indenting close to the [001] orientation. Please note that the figures do not complement each other as a stress threshold t was used to calculate both volume fractions.

Figure 9 :

 9 Figure 9: Evolution of the residual first stress invariant field with indentation orientation varying from [001] to [111]: (a) [2 15 2], (b) [4 15 4], (c) [8 15 8], (d) [13 15 13]. The tensile stress zone is present for several orientations and is confined, almost uniformly, close to the surface for orientations close to [111].

Figure 10 :

 10 Figure 10: Evolution of the in-depth residual stress profile for orientations [001], [4 15 4] and [111]. The profiles are extracted on a path at the geometry center, as sketched at the left. A gaussian filter was applied to remove discretization-induced oscillations. These profiles emphasize the respective reduction and augmentation of the tensile and compressive zones close to the [111] orientation.

  258

Figure 10

 10 Figure10shows that, as the orientation approaches the [111] direction, both the com-

Figure 11 :

 11 Figure 11: Redistribution of the residual first stress invariant field during the indents presented in Figure 3, when indenting in the [001] direction.Figure's numerotation corresponds to the indent number following Figure 3's convention. Tensile stresses are still present after 5 indents. The final residual stress field is highly heterogeneous

Figure 12 :

 12 Figure 11: Redistribution of the residual first stress invariant field during the indents presented in Figure 3, when indenting in the [001] direction.Figure's numerotation corresponds to the indent number following Figure 3's convention. Tensile stresses are still present after 5 indents. The final residual stress field is highly heterogeneous

•

  Single indentation simulations in 43 different crystal orientations revealed that the tensile and compressive stress fields highly depend on the initial crystal orientation. Indentation close to the [001] orientation lead to more tensile stresses and less compressive stresses than for any other crystal orientations. Conversely, more compressive stresses and less tensile stresses were found when indenting close to the [213] orientation. Simulations of five successive indentations in the [001] and [111] directions showed that there is a link between the volume of tensile and compressive stresses induced by one and several indents. More tensile stresses and less compressive stresses produced by one indent leads to shallower compressive stresses and a more heterogeneous stress field after five indents. Accounting for crystal plasticity anisotropy when simulating surface treatments could therefore reveal local tensile residual stresses and local stress heterogeneities in favorably oriented crystals. Provided that the grain sizes are sufficiently large, the initial material microstructure could therefore have an influence on shot peening efficiency. This conclusion should naturally depend on the shot-peening parameters. The underlying deformation mechanisms leading to this orientation-dependent behaviour have not been investigated in the present work. Simulations inputs and results files generated in the present work are therefore provided to the community for further interpretations (see section 7). Stress heterogeneities revealed by this study could affect the fatigue life prediction of shotpeened parts. Local tensile stresses could favor crack propagation during fatigue. Shallower compressive stresses than those expected using macroscopic isotropic laws could result in overestimations of shot-peened part's fatigue life.

Figure 13 :

 13 Figure13: Evolution of f t with the number of elements. Indentation was performed in the crystal's [001] direction and with an angle β = 0 • (see Figure1). Convergence found for these conditions was assumed to hold for any other crystal orientations. Convergence of f t to a value of approximately 13 % is observed for a total of 66 517 elements.

Figure 14 :

 14 Figure 14: Evolution of f t , computed using (a) σ xx , (b) σ yy and (c) σ zz with the number of elements.Indentation was performed in the crystal's [001] direction and with an angle β = 0 • (see Figure1). Convergence found for these conditions was assumed to hold for any other crystal orientations. These figures confirms the previous converged element size for each stress components

Figure 15 :

 15 Figure 15: Evolution of the stress profile along z=0 at a depth of 40µm with teh mesh density. Convergence of the profile can qualitatively be observed for a similar mesh density as previously. This confirms the convergence of the residual stress field

Figure 16 :

 16 Figure16: Evolution of the residual stress field dissimetry y with the indenter mesh density at several depth y. Indenter mesh density on the field symetry can be observed for low mesh densities.

Figure 17 :

 17 Figure 17: Evolution of f t with Zone B to Zone A dimension ratio 2bR . Stabilization of the value of f t can be observed for borders twice larger than the refined zone.

Figure 18 :

 18 Figure18: Evolution of f t with the static friction coefficient value. A slight decrease (9 % of the frictionless value) is observed for high friction coefficients. Variation of f t with friction is therefore considered negligible.

Figure 19 :

 19 Figure 19: (a) Representation of a vector's (green arrow) spherical coordinate angles α and φ. A point on an Inverse Pole Figure (IPF) corresponds to the projection of the green vector on the hemispherical plane. Its coordinates on the IPF can be expressed using α and φ. (b) Mesh used to generate the IPF. Each of the 43 nodes represents the crystal orientation axis aligned with the indentation axis.

Table 1 :

 1 Hardening interaction matrix parameters[START_REF] Méric | Finite element calculations of copper bicrystal specimens submitted to tension-compression tests[END_REF]

	Finite Element Analyses were performed using ABAQUS Finite Element software to solve
	the global mechanical equilibrium. the Zmat module of Z-set Finite Element software was
	h 1	h 2	h 3	h 4	h 5	h 6
	Self-hardening Coplanar Hirth lock Collinear Glissile Lomer junction
	1	4.4	4.75	4.75	4.75	5
					)
		159,300	122,000		81,000

C 11 (MPa) C 22 (MPa) C 44 (MPa

Table 2 :

 2 Elastic coefficients[START_REF] Méric | Finite element calculations of copper bicrystal specimens submitted to tension-compression tests[END_REF]

Acknowledgement

This work was supported by the Natural Sciences and Engineering Research Council (NSERC) discovery grants Canada Research Chairs programs. The authors thank Romain Quey for the productive discussions which allowed to perform the Inverse Pole Figure mapping.

Appendix A. Convergence study

Global mesh convergence

This work mainly focused on spherical indentation-induced residual stress fields. Particular attention was given to the tensile residual stress quantity. Mesh size convergence was therefore based on the volume fraction of element with a positive first stress invariant, f t .

Simulations were performed with eight different substrate mesh sizes. The total number of elements ranged from 20 277 to 280 277.

Figure 13 shows the evolution of f t with the total number of elements for an indentation

Appendix B. Results sensitivity to friction

Most CPFE studies on single crystal indentation modeled the contact as frictionless [START_REF] Alcalá | Reassessing spherical indentation: Contact regimes and mechanical property extractions[END_REF][START_REF] Casals | The duality in mechanical property extractions from vickers and berkovich instrumented indentation experiments[END_REF][START_REF] Petryk | On direct estimation of hardening exponent in crystal plasticity from the spherical indentation test[END_REF][START_REF] Zaafarani | Three-dimensional investigation of the texture and microstructure below a nanoindent in a cu single crystal using 3d ebsd and crystal plasticity finite element simulations[END_REF]. However, according to Marteau et al. [START_REF] Marteau | Review on numerical modeling of instrumented indentation tests for elastoplastic material behavior identification[END_REF], friction coefficient influence has been observed in the litterature when modelling spherical indentation. According to their review, the friction coefficient can have an influence on the pile-up quantity and very locally, on strains and stresses distribution, for low indentation depths.

Friction coefficient influence on the value of f t was therefore studied. Three simulations were performed using different friction coefficients ranging from 0 to 0.4.

Figure 18 shows the evolution of f t with the friction coefficient. A slight decrease is observed for higher friction (9% of the frictionless value). Contact was therefore modelled as frictionless in this study.

Appendix C. Inverse Pole Figure construction

Face-centered cubic crystal orientation can be described by three vectors, corresponding to the three rotated crystal axes. The orientation of one axis can be described by its spherical coordinate angles α and φ (as shown in Figure 19(a)). One point on the IPF represents the projection of the crystal vertical axis y (aligned with the indentation direction) on the hemispherical plane of the unit sphere. The second crystal axis x is chosen using spherical coordinates α x = α y + π/2 and φ x = φ y , so that the third vector always remains in the hemispherical plane.

Considering the face centered cubic crystal symmetry, covering the whole IPF therefore