Multiscale Modelling of Nanoindentation

Nanoindentation is an interesting technique used to probe the local mechanical properties of a material. Although this test has been widely used and developed over the world during the past few years, it remains a lot of uncertainties regarding the interpretation of nanoindentation data. In this study, we propose to simulate the nanoindentation test of FCC single crystals like Cu or Ni using three numerical models. At the lowest scale, molecular dynamics simulations give details of the nucleation of the first dislocations induced by the indentation. At an intermediate scale, discrete dislocation dynamics simulations are performed to study the evolution of the dislocation microstructure during the loading. Finally, at the upper scale, 3D finite element modelling using crystal plasticity constitutive equations give a continuum description of the indentation induced plasticity. It is shown how the different models are interconnected together.


Introduction
Last few years have seen a large development of the nanoindentation techniques. Those advances are mainly motivated by the recent interest for the analysis of thin films behaviour in computer technology or more generally for testing the near-surface properties of materials. The high precision of the latest instruments in the indentation positioning makes it now possible to use the indentation technique to map the mechanical properties of a surface with a submicron resolution. Such an instrument can be seen as a mechanical properties microprobe. Despite these numerous studies, a nanoindentation test is still difficult to interpret as demonstrated by the numerous publications on this topic [1][2][3]. As an example, the indentation size effect which denotes a reverse dependence of the hardness with the penetration depth remains unexplained although the attempts proposed using some strain gradient plasticity concepts [4][5][6]. Numerical simulations could help to fully understand the physics of this test. In the case of crystalline materials, many numerical tools are available. In this paper, it is proposed to simulate nanoindentation of face centred cubic single crystals using three models involving three different scales. At the lowest scale, Molecular Dynamics (MD) simulations give access to the positions of all the atoms contained within a simulation box of few hundred of Angströms wide. Such simulations can explain the early stage of the nanoindentation process, i.e. the nucleation of the first dislocations in a perfect defect free crystal. At an intermediate scale, Discrete Dislocation Dynamics (DDD) simulations can describe the interaction of the dislocation lines and the effect induced on the macroscopic response. For DD simulations, local rules for the nucleation of dislocation loops are introduced thanks to the results obtained by MD. Finally at the upper scale, Finite Element (FE) simulations are performed within the frame of crystal plasticity. For these simulations, the constitutive equations are based on a physical description of the plasticity for which the internal variables are the dislocation densities on all the slip systems. In this modelling, the parameters of the behaviour law are identified by DD simulations.

Molecular Dynamics Simulations
MD simulations are performed using the embedded atom method (EAM) with the inter-atomic potential dedicated to Nickel developed by Voter. The simulated volume is typically a parallelepiped box of 224x284x285 Å 3 (1.675.080 atoms). The indenter is represented by a spherical repulsive potential of radius R=120 Å. The position of the center of the sphere is imposed step by step (displacement control). The normal to the indented surface is chosen as (111). Periodic boundary conditions are applied on the ( 01 1 ) and ( 1 2 1 ) sides of the box. The atoms lying on the bottom surface are forbidden to move along the (111) direction but the other in plane directions are allowed (see Fig.1).  A simulation campaign has been performed using both molecular static and MD simulations which confirmed that the indentation induced dislocations do not depend upon the loading rate. As shown in the loading curves given in Fig.2(b), after an initial elastic stage, discontinuity is observed in the forcedisplacement curves which is the signature of dislocation pop-in. The process for the formation of the dislocations beneath the indenter is rather complex as illustrated by the snapshots printed in Fig. 2. In this figure, only the atoms for which the number of neighbours does not correspond to that of a perfect fcc crystal are plotted. This way, dislocation leading and tail partials appears in red (dark) and atoms in the stacking fault in yellow (clear).
After an elastic deformation of the volume (Fig 2(a)), defects appear beneath the indenter ( Fig. 2(b)). These defects consist of three extrinsic defects lying in the three {111} planes activated by the indentation. Each defect is made of an assembly of three adjacent stacking fault planes. This structure is unstable and rapidly transforms into intrinsic defects (Fig. 2(c)) then into three half-loops of dislocations touching the free surface ( Fig. 2(d)). Each loop is lying in the two slip system sharing a common Burgers vector. These loops expend in the volume and cross-slip events finally emit a closed prismatic loop which can easily glide on the cylinder based on the Burgers vector ( Fig. 2(e-f)). In Fig. 2(f), two loops are completely formed and a third one on the left is prohibited by the periodic boundary conditions which reintroduce the prismatic loop escaping from the right side of the box. When changing the size of the simulation box, we observe that three identical prismatic loops are simultaneously nucleated when indenting along the (111) direction. Moreover, when the indenter is further pushed in the crystal, more and more prismatic loops are introduced with a size directly related to the contact area between the indenter and the surface. This nucleation mechanism will now be used in the DD simulations presented below. For DD, the details of the formation of the prismatic loops are not accounted for and only the final shape and position of the prismatic loops are conserved. This means that at this upper scale, pure perfect prismatic dislocations will be introduced in the simulation box each time the load is sufficient to make them progress in the volume.

926
The Mechanical Behavior of Materials X Fig. 2: Snapshots of the defects induced by nanoindentation obtained by MD simulations.

Discrete Dislocation Dynamics Simulations
DDD simulations are performed using the edge-screw model as described in [8][9]. The boundary conditions corresponding to nanoindentation are enforced by a coupling with the finite element method (see Fig. 3). In this coupling, the stress field applied to each dislocation segment is the superposition of the field generated by all the segments already included in the simulation box as if they were in an infinite medium with an elastic stress field computed by FE. The latter is the solution of a boundary problem where the boundary conditions corresponding to the indentation problem are compensated by the image forces and the displacements induced at the boundary by the dislocations. The algorithm used is the following: the position of the indenter is imposed step by step in quasi-state manner. A Hertz pressure corresponding to a given penetration depth is applied on the top surface. This induces a force on the indenter which is purely elastic and close to the Hertz theory (See Fig. 3(b)). Then prismatic loops are introduced beneath the indenter in order to decrease the force toward the mater curves obtained experimentally. When the experimental curve is reached, dislocations are allowed to move until 50nm (a (b (f) (e (d (c they all reach a stable position. Next iteration is then applied: the load is increased, dislocations are introduced and the dislocation structure is equilibrated. Fig. 4 shows snapshots of the dislocation structure developed during the nanoindentation along (001) of single crystal of Cu up to 50 nanometers with a spherical indenter of radius R=420nm [9]. Fig. 4: Snapshots of the dislocation structure induced by (001) indentation in Cu [9].
It was observed that the obtained dislocation microstructure is contained within half a sphere of diameter 1.5 micrometers which is in good agreement with Johnson's equation [10]. An estimate of the hardness can be computed from the simulation by computing the yield stress induced by the dislocation density included in the half-sphere. The density was measured to ρ=10 16 m -2 from which a hardness of H=1.4GPa was computed. Experimentally, the hardness found for a (001) indentation of a Cu single crystal was H=1.4 ±0.2 GPa which shows that the DDD modelling is consistent with the reality.

Finite Element simulations
FE simulations are performed using ABAQUS v6.4 with a set of constitutive equations based on dislocation densities. In this model, the plastic flow linking the stress on a given system (s) to the strain rate on the same system is given by the following three equations [11]: The flow law relates the shear strain rate ) The evolution of the dislocation densities is given by a relation involving both a storage and an annihilation term: where K is a dislocation mean free path, y c an annihilation distance and [b] an interaction matrix. Note that this equation gives a saturation value for the dislocation density when the storage term is equal to the annihilation term. Although these three equations are closely related to dislocation motion, the value to give to the parameters is hard to quantify since an average has been done over the entire population of dislocations of a system (s). Thus, DDD simulations have been performed in order to compute these values in the case of Cu single crystal [12].
The obtained constitutive equations have been introduced in ABAQUS v6.4 using the open VUMAT and UMAT procedures; the first one being for implicit integration and the later one for explicit algorithm. Simulation of nanoindentation test of a Cu single crystal along (001) has been performed using the explicit scheme with the same geometry as the one chosen for DDD in previous section, i.e. an indenter radius of R=420nm and a maximum penetration depth of 50 nanometres. The mesh used here consists of 10.308 eight nodes elements refined beneath the indenter as shown in Fig. 5. The indenter is supposed to be infinitely rigid and the contact is frictionless. A mapping of the Von Mises stresses developed beneath the indenter is shown in Fig5(b) for the maximum penetration depth. One can observe the four-fold symmetry enforced by the four Burgers vectors (eight slip systems) activated when the indenter is pressed along the (001) direction. This observation is confirmed by the displacement field obtained for the same depth as shown in Fig.  6(a) where the same four-fold symmetry is observed for UZ, the displacement projected along the indentation axis. Fig. 6(b) shows that a pile-up phenomenon of 10 nanometres height is observed along the [110] directions.

Concluding Remarks
The three numerical models presented here are inter-connected together since the FE simulations are based on constitutive equations which are identified by DDD simulations, and DDD simulations use a nucleation criterion derived from the MD simulations. Work is now in progress in order to use these three tools to understand the origin of the indent size effect.