- Email: [email protected]
Multiscale modelling of nanoindentation test in copper crystal
Available online at www.sciencedirect.com
Engineering Fracture Mechanics 75 (2008) 3755–3762 www.elsevier.com/locate/engfracmech
Multiscale modelling of nanoindentation test in copper crystal ˇ erny´, Jaroslav Pokluda Jana Hornı´kova´ *, Pavel Sˇandera, Miroslav C Brno University of Technology, Technicka´ 2, CZ-61669 Brno, Czech Republic Received 21 August 2007; received in revised form 24 October 2007; accepted 28 October 2007 Available online 7 November 2007
Abstract The nanoindentation test in the dislocation free crystal of copper is simulated by a nonlinear elastic finite element analysis coupled with both ab initio calculations of the ideal shear strength and crystallographic considerations. The onset of microplasticity, associated with the pop-in effect identified in experimental nanoindentation tests (creation of first dislocations), is assumed to be related to the moment of achieving the value of the ideal shear strength for the copper crystal. Calculated values of the critical indentation depth lie within the range of experimentally observed pop-ins in the copper crystal. The related indentation load is somewhat lower than that observed in the experiment. 2007 Elsevier Ltd. All rights reserved. Keywords: Nanoindentation; Ab initio calculation; Ideal shear strength; Copper crystal; Finite element analysis
1. Introduction The nanoindentation is considered to be a very promising experimental approach to measuring the ideal shear strength since the stressed volume beneath the sharp indenter may be defect-free. The local shear component of the stress reaches its maximum value at certain (not very large) penetration depth of the indenter into the bulk. The value of this stress can approach the ideal shear strength sid (ISS) and, consequently, it can be high enough to nucleate dislocations [1–3]. This process might be detected as a pop-in on the nanoindentation load–displacement curve. Indeed, some authors already reported such behaviour and tried to compare the related stress to the theoretically determined ideal shear stress [4]. A relevant model of the nanoindentation test must be based on joined approaches on three different levels: (i) Atomistic (calculation of ISS); (ii) mesoscopic (crystallography); (iii) macroscopic (finite element analysis of stress–strain response). As a rule, however, the experimental values of the pop-in deformation lie above the values obtained by majority of recent models. The reasons for that discrepancy may have several origins. First, the simple theory must be modified in order to reflect the real stress state at the point of maximum shear under the indenter. In particular, the influence of the compressive normal load on the slip plane on the ideal shear strength value is to be *
Corresponding author. Tel.: +420 54114 2866; fax: +420 54114 2842. ˇ erny´), E-mail addresses: [email protected] (J. Hornı´kova´), [email protected] (P. Sˇandera), [email protected] (M. C [email protected] (J. Pokluda). 0013-7944/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2007.10.016
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
3756
considered in the ab initio analysis. Second, a solution reflecting the three-dimensionality of the nanoindentation test is to be obtained since the real stress state under the indenter must be strongly different from that two-dimensional mostly utilized in recent analysis. Third, the nonlinearity in the stress–strain relation is to be taken into account. Fourth, a correction is needed in order to match the orientation of the appropriate crystallographic plane, i.e. to meet the actual resolved shear stress. Fifth, friction forces between the surfaces of the substrate and the indenter should be also included into the model. Sixth, the anisotropy of the elastic response of the crystal is to be considered. Seventh, the process of dislocation nucleation is to be appropriately reflected in the atomistic model of the shear deformation of the perfect crystal. The authors of the model [4] claim that all the above mentioned corrections, except for the last one, were considered in their analysis. However, no detailed descriptions of related procedures are reported in their paper. Therefore, the relevancy of these results remains rather questionable. In our previous work [5], the first two corrections were included in the model. As a result, the theoretical values obtained lie at a lower limit of the range of experimental pop-ins. In this paper, the third, fourth and fifth corrections (the nonlinearity, the crystallography and the friction) are fully taken into account. Moreover, the influence of anisotropy is also partially included by respecting the computed stress–strain curve of the copper crystal in the [0 0 1] loading direction. 2. Description of procedures used in the multiscale model A three-dimensional isotropic FEM analysis was performed by using the finite element ANSYS code. A frictionless sphere–conical indenter with the radius of 0.5 lm in agreement with the experiment [3] is pressed into a 5 lm thick substrate disc with the radius of 10 lm (see Fig. 1). In the vicinity of the interface, the sphere is meshed with elements approximately 0.5 nm and the mesh of the substrate is refined from 1 lm at the outer edge of the disc to 0.5 nm directly beneath the indenter. The size of the elements near the contact is an order magnitude smaller in comparison to the depth of penetration (see Fig. 2). Owing to the axial symmetry of the 3D indentation model, only 2D section could be analyzed, as depicted in Fig. 2. The Young modulus E of the diamond indenter was taken to be of 1141 GPa and the Poisson ratio m = 0.07. To model a contact area between the indenter tip and the specimen one has to identify where the contact might occur during the deformation. Once identifying the potential contact surfaces, they are to be defined via a target and the contact elements, which will then track the kinematics of the deformation process. In the contact area between boundaries of the indenter and the solid, the indenter was established as the ‘‘target’’ surface, and the other one as the ‘‘contact’’ surface. The contact surface in the vicinity of the indenter tip is created by a set of contact elements that comprise the surface of the deformable body. It is formed by contact pairs – i.e. the ‘‘target surface’’ and the ‘‘contact surface’’. These contact elements have the same geometric characteristics as the underlying elements of the deformable body. The contact detection points are located at the Gauss integration points, which are interior to the contact surface elements. The contact
a
b
c
r
indenter
r
a
h
5μm
α
20μm
Fig. 1. The geometrical configuration of the model. (a) The cross-section of the nanoindenter and the substrate disc. (b) Scheme of the tip of the diamond sphere–conical indenter with the semi-angle a = 68 and the tip radius r = 500 nm. (c) Scheme of the indenter penetration into the copper crystal and the penetration depth h.
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
3757
Fig. 2. An example of the finite element mesh near the contact boundary.
surface element is constrained against penetration into the target surface at its integration points. However, the target surface can, in principle, penetrate through into the contact surface, see Fig. 3. The code ANSYS updates tangential contact stiffness based on current contact normal pressure and maximum allowable elastic slip. Beside the frictionless analysis, we also made an analysis involving the friction. In the Coulomb friction model, two contacting surfaces can carry shear stresses up to a certain magnitude across their interface before they start sliding relative to each other. This state is known as sticking. The Coulomb friction model defines an equivalent shear stress s, at which sliding on the surface begins as a fraction of the contact pressure. Once the shear stress is exceeded, the two surfaces will slide relative to each other. The sticking/sliding calculations determine when a point changes from sticking to sliding or vice versa. As shown in Fig. 4a, for the condition of full slip (no friction), under loading, points of both contact surface (substrate disc) and target surface (indenter) move inward toward the axis of symmetry under the influence of the applied forces Fa. Movement of points within the substrate disc generates ‘‘internal’’ forces Fs (i.e. from the stresses set up in the material) which are proportional to the relative displacement. Movement ceases when the internal forces Fs balance the applied forces. As shown in Fig. 4b, the case of no slip (i.e. full adhesive contact), points on the contact surface want to move inward under the influence of the applied forces Fa but are prevented from doing so by frictional forces Ff (the friction force Ff balances the applied force Fa). In case of partial slip, Fig. 4c, points on the
indenter target surface segment Gaussian integration point
contact surface element
contact surface segment deformed body (substrate disc)
Fig. 3. Gaussian integration points located interior to the contact surface elements.
3758
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
Ff
Fa Fs
Fa
Fa Fs
Ff
Fig. 4. Contact points on the surfaces of both the indenter and the substrate disc. (a) Full slip. (b) No slip. (c) Partial slip.
contact surface want also to move inward under the influence of the applied force Fa. Some points are prevented from doing so by frictional forces, which, due to the local magnitude of the normal forces, are large enough to balance the applied forces. For other points, however, the applied forces are greater than the frictional force and those points do move inwards – slip occurs between the surfaces. For those points that have slipped, the frictional force has already reached its maximum value. Internal forces can still increase with increasing load. Relative movement occurs until the internal force Fs plus the maximum frictional force Ff opposes the applied force Fa. The friction force is now applied by a new point on the target surface, which has come into contact with the point on the contact surface. Now, at full load, the applied force at a point that has slipped is balanced by the sum of the maximum frictional force and the internal force. The nonlinearity and, to a certain extent, also the anisotropy of the stress–strain elastic response was taken into account by a multilinear approximation of the stress–strain curve calculated ab initio for the [0 0 1] uniaxial tension–compression of the copper crystal [6] (see Fig. 5). This curve was utilized in the ANSYS code procedure [7] as an equivalent, von Mises, stress–strain dependence. The dependence of sid on the superimposed normal stress rn acting perpendicularly to {1 1 1} slip planes was also computed by means of the ab initio approach. The h2 1 1i {1 1 1} shear system was selected as the preferred (primary) one, because it exhibits the lowest shear strength. The electronic structure calculations in the ab initio analyses were performed using a plane wave code VASP (Vienna Ab initio Simulation Package) [8]. This code uses ultra-soft pseudo-potentials of the Vanderbilt type [9]. The exchange-correlation energy was evaluated by using the generalized-gradient approximation (GGA) of Perdew and Wang [10]. After proper convergence tests, 18 · 18 · 18 Monkhorst–Pack k-point mesh was found to be satisfactory for integration in the Brillouin zone. The energy cut-off of the plane-wave expansion was increased to 300 eV in order to obtain a reliable stress tensor. The solution was considered to be self-consistent when the energy difference of two consequent iterations was smaller than 10 leV. The crystal energy per atom E = E(u, c) was computed as a function of two independent parameters: the normalized interplanar distance c (the interplanar distance divided by the equilibrium lattice parameter a0) and the normalized plane shift u (see Fig. 6). The latter quan-
Fig. 5. The stress–strain response of the copper crystal to the uniaxial loading in the [0 0 1] direction.
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
displaced atoms
a0 u
111
equilibrium positions
3759
a0 c
211 3 a 8 0
Fig. 6. An illustration of the basic shear variables: The normalized shift u of two adjacent {1 1 1} planes and the normalized interplanar distance c.
pffiffiffiffiffiffiffiffi tity was scaled itsffiffiffiffiffiffiffi values of 0 and 3=2 correspond to the identical fcc state. On the other hand, u pffiffiffiffiffiffiffiffiso that p ffi values of 1=6 and 2=3 lead also to the fcc state but of an opposite stacking order. The corresponding pffiffiffiffiffiffiffiffi value of the normalized interplanar distance for the cubic state is of 1=3. In our model, the {1 1 1} planes remain undistorted (atomic positions within the planes are fixed) along the entire shear deformation path. The stress acting perpendicularly to the slip planes can be evaluated as a first derivative of the crystal energy per atom with respect to the interplanar distance c at any plane shift as 1 dE rn ¼ pffiffiffi ; 3V dc 0
u¼const
where the equilibrium atomic volume V 0 ¼ a30 =4. Similarly, the shear stress applied to slip planes can be computed as follows: 1 dE dE rn s ¼ pffiffiffi : du 3V du rn ¼const
0
The function c = c(u) is to be considered here since the distance c for rn = const varies with the shift u. If no other instability precedes, the maximal value of shear stress (smax) can be considered to be ISS under the given normal stress. Table 1 contains the ideal shear strength values (under the zero normal stress) computed by various authors including that obtained in our present calculations. The computed dependence of the shear strength on the normal stress is depicted in Fig. 7. This dependence is nearly linear and can be simply parameterized as smax ¼ 3:166 0:122rn : The global procedure simulates, step by step, the penetration of the indenter into the copper crystal. In order to identify the appropriate crystallographic plane in which the latter condition for the dislocation emission is firstly reached, the activity in the shear system ð1 1 1Þ½1 1 2 was considered in the calculations. With regard to the cylindrical symmetry of the 3D model, the stress tensor transformation was performed to obtain all possible positions of the crystallographic systems rotating relatively to the main coordinate system. This rotation enabled us to compute the values of s and rn as functions of the angle /. Consequently, the maximum of the ratio #i(/) = si(/)/sid,i(/) was searched on the whole circle going through each node of the mesh. The highest value of that ratio over all nodes and all angles is denoted as #max. These values are related to each individual deformation step characterised by the depth h of the penetration.
Table 1 The ideal shear strength of Cu along with available literature data Calculation
Present
Ogata et al. [11]
Krenn et al. [12]
Sˇandera and Pokluda [13]
smax (GPa)
3.2
2.2
2.7
2.9
3760
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
Fig. 7. The ideal shear strength smax as a function of the normal stress rn.
3. Results and discussion The results obtained by means of the frictionless analysis and that involving the friction revealed to be practically identical. This is not too much surprising with regard to results already reported by other authors (e.g. [4]). The values of #max as a function of the penetration depth h are depicted in the Fig. 8. When the value of #max exceeds 1, the condition for the dislocation emission is safely fulfilled. This corresponds to the depth h 12 nm. Fig. 9 shows the calculated contours of # (as a function of coordinates) at the moment when #max reaches 1.03. This value lies within the black area, which is slightly off the loading axis in agreement with other published models of the nanoindentation tests. The situation in Fig. 9 corresponds to the calculated
Fig. 8. The dependence of the parameter #max on the penetration depth h. The predicted moment of the first emission of dislocation loops corresponds to exceeding the critical value #max = 1 for h 12 nm.
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
3761
Fig. 9. Contours of # in the indentation step corresponding to the network maximum of #max G 1.03 (h 12 nm) beneath the indenter (colour variant of the same picture can be obtained from on-line version).
indentation load of 100 lN and to the indenter displacement of 12 nm. These values have to be comparable with experimentally observed pop-ins (150–270 lN and 10–20 nm [3]). While the computed displacement value lies within the experimental range, the calculated indentation load is somewhat lower. This might be a consequence of the fact that the model does not appropriately reflect the process of dislocation nucleation in the perfect crystal on an atomistic level. Moreover, the dependence of the ideal shear stress on the normal stress component should be refined by including the in-plane relaxation. An improvement of the model in this way is under development. In spite of that slight inaccuracy, the results reveal that the nanoindentation may provide a very good tool for measuring the ideal shear strength. 4. Conclusion The nanoindentation test in the dislocation free crystal of copper was simulated by utilizing a multiscale analysis. The onset of microplasticity, associated with the pop-in effect identified in experimental nanoindentation tests (creation of first dislocation loops), is assumed to be related to the moment when the value of the ideal shear strength for the copper crystal was reached. In particular, the influence of the compressive normal load on the shear plane on the ideal shear strength value, the three-dimensionality of the nanoindentation test, the nonlinearity in the stress–strain relation, the orientation of relevant crystallographic planes, the anisotropy of the elastic response of the crystal and the friction stress between the indenter and the substrate were considered in the model. The mechanical characteristics of the perfect copper crystal were calculated by using the ab initio approach. The three-dimensional isotropic FEM analysis, based on the finite element ANSYS code, was used to simulate the development of the stress–strain field in the substrate. The computed displacement value lies within the experimentally observed range of the pop-in effect while the related indentation load is somewhat lower. This might be particularly a consequence of the fact that the model does not appropriately reflect the process of dislocation nucleation in the perfect crystal on the atomistic level. Nevertheless, the results reveal that the nanoindentation test can serve as a sufficiently precise tool for experimental determination of the ideal shear strength values. Acknowledgements This research was supported by the Ministry of Education and Youth of the Czech Republic under Grants No. OC148/P19, MSM0021630518 and the Czech Science Foundation under the Project No. GA106/05/0274.
3762
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
References [1] Chiu YL, Ngan AHW. A TEM investigation on indentation plastic zones in Ni3Al(Cr,B) single crystals. Acta Mater 2002;50(10):2677–91. [2] Go¨ken M, Kempf M, Nix WD. Hardness and modulus of the lamellar microstructure in PST-TiAl studied by nanoindentations and AFM. Acta Mater 2001;49(5):903–11. [3] Go¨ken M, Kempf M. Pop-ins in nanoindentations – the initial yield point. Z Metallkd 2001;92(9):1061–7. [4] Krenn CR, Roundy D, Cohen ML, Chrzan DC, Morris Jr JW. Connecting atomistic and experimental estimates of ideal strength. Phys Rev B 2002;65(13):134111. ˇ erny´ M, Sˇandera P, Pokluda J. Onset of microplasticity in copper crystal during nanoindentation. Key Engng Mater [5] Hornı´kova´ J, C 2007;348–349:8001–4. ˇ erny´ M, Sˇob M, Pokluda J, Sˇandera P. Ab initio calculations of ideal tensile strength and mechanical stability in copper. J Phys: [6] C Cond Matter 2004;16:1045–52. [7] ANSYS Users Manual. Revision 5.0, I–IV. Swanson Analysis System, Houston, USA, 1992. [8] Kresse G, Furthmu¨ller J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys Rev B 1996;54(16):11169–86. [9] Vanderbilt D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys Rev B 1990;41(11):7892–5. [10] Perdew JP, Wang Y. Accurate and simple analytic representation of the electron-gas correlation-energy. Phys Rev B 1992;45(23):13244–9. [11] Ogata S, Li J, Hirosaki N, Shibutani Y, Yip S. Ideal shear strain of metals and ceramics. Phys Rev B 2004;70(10):104104. [12] Krenn CR, Roundy D, Morris Jr JW, Cohen ML. The non-linear elastic behavior and ideal shear strength of Al and Cu. Mater Sci Engng A 2001;317(1–2):44–8. [13] Sˇandera P, Pokluda J. Ideal strength of cubic crystals under triaxial tension. Metal Mater 1994;32(4):180–93.
Engineering Fracture Mechanics 75 (2008) 3755–3762 www.elsevier.com/locate/engfracmech
Multiscale modelling of nanoindentation test in copper crystal ˇ erny´, Jaroslav Pokluda Jana Hornı´kova´ *, Pavel Sˇandera, Miroslav C Brno University of Technology, Technicka´ 2, CZ-61669 Brno, Czech Republic Received 21 August 2007; received in revised form 24 October 2007; accepted 28 October 2007 Available online 7 November 2007
Abstract The nanoindentation test in the dislocation free crystal of copper is simulated by a nonlinear elastic finite element analysis coupled with both ab initio calculations of the ideal shear strength and crystallographic considerations. The onset of microplasticity, associated with the pop-in effect identified in experimental nanoindentation tests (creation of first dislocations), is assumed to be related to the moment of achieving the value of the ideal shear strength for the copper crystal. Calculated values of the critical indentation depth lie within the range of experimentally observed pop-ins in the copper crystal. The related indentation load is somewhat lower than that observed in the experiment. 2007 Elsevier Ltd. All rights reserved. Keywords: Nanoindentation; Ab initio calculation; Ideal shear strength; Copper crystal; Finite element analysis
1. Introduction The nanoindentation is considered to be a very promising experimental approach to measuring the ideal shear strength since the stressed volume beneath the sharp indenter may be defect-free. The local shear component of the stress reaches its maximum value at certain (not very large) penetration depth of the indenter into the bulk. The value of this stress can approach the ideal shear strength sid (ISS) and, consequently, it can be high enough to nucleate dislocations [1–3]. This process might be detected as a pop-in on the nanoindentation load–displacement curve. Indeed, some authors already reported such behaviour and tried to compare the related stress to the theoretically determined ideal shear stress [4]. A relevant model of the nanoindentation test must be based on joined approaches on three different levels: (i) Atomistic (calculation of ISS); (ii) mesoscopic (crystallography); (iii) macroscopic (finite element analysis of stress–strain response). As a rule, however, the experimental values of the pop-in deformation lie above the values obtained by majority of recent models. The reasons for that discrepancy may have several origins. First, the simple theory must be modified in order to reflect the real stress state at the point of maximum shear under the indenter. In particular, the influence of the compressive normal load on the slip plane on the ideal shear strength value is to be *
Corresponding author. Tel.: +420 54114 2866; fax: +420 54114 2842. ˇ erny´), E-mail addresses: [email protected] (J. Hornı´kova´), [email protected] (P. Sˇandera), [email protected] (M. C [email protected] (J. Pokluda). 0013-7944/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2007.10.016
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
3756
considered in the ab initio analysis. Second, a solution reflecting the three-dimensionality of the nanoindentation test is to be obtained since the real stress state under the indenter must be strongly different from that two-dimensional mostly utilized in recent analysis. Third, the nonlinearity in the stress–strain relation is to be taken into account. Fourth, a correction is needed in order to match the orientation of the appropriate crystallographic plane, i.e. to meet the actual resolved shear stress. Fifth, friction forces between the surfaces of the substrate and the indenter should be also included into the model. Sixth, the anisotropy of the elastic response of the crystal is to be considered. Seventh, the process of dislocation nucleation is to be appropriately reflected in the atomistic model of the shear deformation of the perfect crystal. The authors of the model [4] claim that all the above mentioned corrections, except for the last one, were considered in their analysis. However, no detailed descriptions of related procedures are reported in their paper. Therefore, the relevancy of these results remains rather questionable. In our previous work [5], the first two corrections were included in the model. As a result, the theoretical values obtained lie at a lower limit of the range of experimental pop-ins. In this paper, the third, fourth and fifth corrections (the nonlinearity, the crystallography and the friction) are fully taken into account. Moreover, the influence of anisotropy is also partially included by respecting the computed stress–strain curve of the copper crystal in the [0 0 1] loading direction. 2. Description of procedures used in the multiscale model A three-dimensional isotropic FEM analysis was performed by using the finite element ANSYS code. A frictionless sphere–conical indenter with the radius of 0.5 lm in agreement with the experiment [3] is pressed into a 5 lm thick substrate disc with the radius of 10 lm (see Fig. 1). In the vicinity of the interface, the sphere is meshed with elements approximately 0.5 nm and the mesh of the substrate is refined from 1 lm at the outer edge of the disc to 0.5 nm directly beneath the indenter. The size of the elements near the contact is an order magnitude smaller in comparison to the depth of penetration (see Fig. 2). Owing to the axial symmetry of the 3D indentation model, only 2D section could be analyzed, as depicted in Fig. 2. The Young modulus E of the diamond indenter was taken to be of 1141 GPa and the Poisson ratio m = 0.07. To model a contact area between the indenter tip and the specimen one has to identify where the contact might occur during the deformation. Once identifying the potential contact surfaces, they are to be defined via a target and the contact elements, which will then track the kinematics of the deformation process. In the contact area between boundaries of the indenter and the solid, the indenter was established as the ‘‘target’’ surface, and the other one as the ‘‘contact’’ surface. The contact surface in the vicinity of the indenter tip is created by a set of contact elements that comprise the surface of the deformable body. It is formed by contact pairs – i.e. the ‘‘target surface’’ and the ‘‘contact surface’’. These contact elements have the same geometric characteristics as the underlying elements of the deformable body. The contact detection points are located at the Gauss integration points, which are interior to the contact surface elements. The contact
a
b
c
r
indenter
r
a
h
5μm
α
20μm
Fig. 1. The geometrical configuration of the model. (a) The cross-section of the nanoindenter and the substrate disc. (b) Scheme of the tip of the diamond sphere–conical indenter with the semi-angle a = 68 and the tip radius r = 500 nm. (c) Scheme of the indenter penetration into the copper crystal and the penetration depth h.
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
3757
Fig. 2. An example of the finite element mesh near the contact boundary.
surface element is constrained against penetration into the target surface at its integration points. However, the target surface can, in principle, penetrate through into the contact surface, see Fig. 3. The code ANSYS updates tangential contact stiffness based on current contact normal pressure and maximum allowable elastic slip. Beside the frictionless analysis, we also made an analysis involving the friction. In the Coulomb friction model, two contacting surfaces can carry shear stresses up to a certain magnitude across their interface before they start sliding relative to each other. This state is known as sticking. The Coulomb friction model defines an equivalent shear stress s, at which sliding on the surface begins as a fraction of the contact pressure. Once the shear stress is exceeded, the two surfaces will slide relative to each other. The sticking/sliding calculations determine when a point changes from sticking to sliding or vice versa. As shown in Fig. 4a, for the condition of full slip (no friction), under loading, points of both contact surface (substrate disc) and target surface (indenter) move inward toward the axis of symmetry under the influence of the applied forces Fa. Movement of points within the substrate disc generates ‘‘internal’’ forces Fs (i.e. from the stresses set up in the material) which are proportional to the relative displacement. Movement ceases when the internal forces Fs balance the applied forces. As shown in Fig. 4b, the case of no slip (i.e. full adhesive contact), points on the contact surface want to move inward under the influence of the applied forces Fa but are prevented from doing so by frictional forces Ff (the friction force Ff balances the applied force Fa). In case of partial slip, Fig. 4c, points on the
indenter target surface segment Gaussian integration point
contact surface element
contact surface segment deformed body (substrate disc)
Fig. 3. Gaussian integration points located interior to the contact surface elements.
3758
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
Ff
Fa Fs
Fa
Fa Fs
Ff
Fig. 4. Contact points on the surfaces of both the indenter and the substrate disc. (a) Full slip. (b) No slip. (c) Partial slip.
contact surface want also to move inward under the influence of the applied force Fa. Some points are prevented from doing so by frictional forces, which, due to the local magnitude of the normal forces, are large enough to balance the applied forces. For other points, however, the applied forces are greater than the frictional force and those points do move inwards – slip occurs between the surfaces. For those points that have slipped, the frictional force has already reached its maximum value. Internal forces can still increase with increasing load. Relative movement occurs until the internal force Fs plus the maximum frictional force Ff opposes the applied force Fa. The friction force is now applied by a new point on the target surface, which has come into contact with the point on the contact surface. Now, at full load, the applied force at a point that has slipped is balanced by the sum of the maximum frictional force and the internal force. The nonlinearity and, to a certain extent, also the anisotropy of the stress–strain elastic response was taken into account by a multilinear approximation of the stress–strain curve calculated ab initio for the [0 0 1] uniaxial tension–compression of the copper crystal [6] (see Fig. 5). This curve was utilized in the ANSYS code procedure [7] as an equivalent, von Mises, stress–strain dependence. The dependence of sid on the superimposed normal stress rn acting perpendicularly to {1 1 1} slip planes was also computed by means of the ab initio approach. The h2 1 1i {1 1 1} shear system was selected as the preferred (primary) one, because it exhibits the lowest shear strength. The electronic structure calculations in the ab initio analyses were performed using a plane wave code VASP (Vienna Ab initio Simulation Package) [8]. This code uses ultra-soft pseudo-potentials of the Vanderbilt type [9]. The exchange-correlation energy was evaluated by using the generalized-gradient approximation (GGA) of Perdew and Wang [10]. After proper convergence tests, 18 · 18 · 18 Monkhorst–Pack k-point mesh was found to be satisfactory for integration in the Brillouin zone. The energy cut-off of the plane-wave expansion was increased to 300 eV in order to obtain a reliable stress tensor. The solution was considered to be self-consistent when the energy difference of two consequent iterations was smaller than 10 leV. The crystal energy per atom E = E(u, c) was computed as a function of two independent parameters: the normalized interplanar distance c (the interplanar distance divided by the equilibrium lattice parameter a0) and the normalized plane shift u (see Fig. 6). The latter quan-
Fig. 5. The stress–strain response of the copper crystal to the uniaxial loading in the [0 0 1] direction.
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
displaced atoms
a0 u
111
equilibrium positions
3759
a0 c
211 3 a 8 0
Fig. 6. An illustration of the basic shear variables: The normalized shift u of two adjacent {1 1 1} planes and the normalized interplanar distance c.
pffiffiffiffiffiffiffiffi tity was scaled itsffiffiffiffiffiffiffi values of 0 and 3=2 correspond to the identical fcc state. On the other hand, u pffiffiffiffiffiffiffiffiso that p ffi values of 1=6 and 2=3 lead also to the fcc state but of an opposite stacking order. The corresponding pffiffiffiffiffiffiffiffi value of the normalized interplanar distance for the cubic state is of 1=3. In our model, the {1 1 1} planes remain undistorted (atomic positions within the planes are fixed) along the entire shear deformation path. The stress acting perpendicularly to the slip planes can be evaluated as a first derivative of the crystal energy per atom with respect to the interplanar distance c at any plane shift as 1 dE rn ¼ pffiffiffi ; 3V dc 0
u¼const
where the equilibrium atomic volume V 0 ¼ a30 =4. Similarly, the shear stress applied to slip planes can be computed as follows: 1 dE dE rn s ¼ pffiffiffi : du 3V du rn ¼const
0
The function c = c(u) is to be considered here since the distance c for rn = const varies with the shift u. If no other instability precedes, the maximal value of shear stress (smax) can be considered to be ISS under the given normal stress. Table 1 contains the ideal shear strength values (under the zero normal stress) computed by various authors including that obtained in our present calculations. The computed dependence of the shear strength on the normal stress is depicted in Fig. 7. This dependence is nearly linear and can be simply parameterized as smax ¼ 3:166 0:122rn : The global procedure simulates, step by step, the penetration of the indenter into the copper crystal. In order to identify the appropriate crystallographic plane in which the latter condition for the dislocation emission is firstly reached, the activity in the shear system ð1 1 1Þ½1 1 2 was considered in the calculations. With regard to the cylindrical symmetry of the 3D model, the stress tensor transformation was performed to obtain all possible positions of the crystallographic systems rotating relatively to the main coordinate system. This rotation enabled us to compute the values of s and rn as functions of the angle /. Consequently, the maximum of the ratio #i(/) = si(/)/sid,i(/) was searched on the whole circle going through each node of the mesh. The highest value of that ratio over all nodes and all angles is denoted as #max. These values are related to each individual deformation step characterised by the depth h of the penetration.
Table 1 The ideal shear strength of Cu along with available literature data Calculation
Present
Ogata et al. [11]
Krenn et al. [12]
Sˇandera and Pokluda [13]
smax (GPa)
3.2
2.2
2.7
2.9
3760
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
Fig. 7. The ideal shear strength smax as a function of the normal stress rn.
3. Results and discussion The results obtained by means of the frictionless analysis and that involving the friction revealed to be practically identical. This is not too much surprising with regard to results already reported by other authors (e.g. [4]). The values of #max as a function of the penetration depth h are depicted in the Fig. 8. When the value of #max exceeds 1, the condition for the dislocation emission is safely fulfilled. This corresponds to the depth h 12 nm. Fig. 9 shows the calculated contours of # (as a function of coordinates) at the moment when #max reaches 1.03. This value lies within the black area, which is slightly off the loading axis in agreement with other published models of the nanoindentation tests. The situation in Fig. 9 corresponds to the calculated
Fig. 8. The dependence of the parameter #max on the penetration depth h. The predicted moment of the first emission of dislocation loops corresponds to exceeding the critical value #max = 1 for h 12 nm.
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
3761
Fig. 9. Contours of # in the indentation step corresponding to the network maximum of #max G 1.03 (h 12 nm) beneath the indenter (colour variant of the same picture can be obtained from on-line version).
indentation load of 100 lN and to the indenter displacement of 12 nm. These values have to be comparable with experimentally observed pop-ins (150–270 lN and 10–20 nm [3]). While the computed displacement value lies within the experimental range, the calculated indentation load is somewhat lower. This might be a consequence of the fact that the model does not appropriately reflect the process of dislocation nucleation in the perfect crystal on an atomistic level. Moreover, the dependence of the ideal shear stress on the normal stress component should be refined by including the in-plane relaxation. An improvement of the model in this way is under development. In spite of that slight inaccuracy, the results reveal that the nanoindentation may provide a very good tool for measuring the ideal shear strength. 4. Conclusion The nanoindentation test in the dislocation free crystal of copper was simulated by utilizing a multiscale analysis. The onset of microplasticity, associated with the pop-in effect identified in experimental nanoindentation tests (creation of first dislocation loops), is assumed to be related to the moment when the value of the ideal shear strength for the copper crystal was reached. In particular, the influence of the compressive normal load on the shear plane on the ideal shear strength value, the three-dimensionality of the nanoindentation test, the nonlinearity in the stress–strain relation, the orientation of relevant crystallographic planes, the anisotropy of the elastic response of the crystal and the friction stress between the indenter and the substrate were considered in the model. The mechanical characteristics of the perfect copper crystal were calculated by using the ab initio approach. The three-dimensional isotropic FEM analysis, based on the finite element ANSYS code, was used to simulate the development of the stress–strain field in the substrate. The computed displacement value lies within the experimentally observed range of the pop-in effect while the related indentation load is somewhat lower. This might be particularly a consequence of the fact that the model does not appropriately reflect the process of dislocation nucleation in the perfect crystal on the atomistic level. Nevertheless, the results reveal that the nanoindentation test can serve as a sufficiently precise tool for experimental determination of the ideal shear strength values. Acknowledgements This research was supported by the Ministry of Education and Youth of the Czech Republic under Grants No. OC148/P19, MSM0021630518 and the Czech Science Foundation under the Project No. GA106/05/0274.
3762
J. Hornı´kova´ et al. / Engineering Fracture Mechanics 75 (2008) 3755–3762
References [1] Chiu YL, Ngan AHW. A TEM investigation on indentation plastic zones in Ni3Al(Cr,B) single crystals. Acta Mater 2002;50(10):2677–91. [2] Go¨ken M, Kempf M, Nix WD. Hardness and modulus of the lamellar microstructure in PST-TiAl studied by nanoindentations and AFM. Acta Mater 2001;49(5):903–11. [3] Go¨ken M, Kempf M. Pop-ins in nanoindentations – the initial yield point. Z Metallkd 2001;92(9):1061–7. [4] Krenn CR, Roundy D, Cohen ML, Chrzan DC, Morris Jr JW. Connecting atomistic and experimental estimates of ideal strength. Phys Rev B 2002;65(13):134111. ˇ erny´ M, Sˇandera P, Pokluda J. Onset of microplasticity in copper crystal during nanoindentation. Key Engng Mater [5] Hornı´kova´ J, C 2007;348–349:8001–4. ˇ erny´ M, Sˇob M, Pokluda J, Sˇandera P. Ab initio calculations of ideal tensile strength and mechanical stability in copper. J Phys: [6] C Cond Matter 2004;16:1045–52. [7] ANSYS Users Manual. Revision 5.0, I–IV. Swanson Analysis System, Houston, USA, 1992. [8] Kresse G, Furthmu¨ller J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys Rev B 1996;54(16):11169–86. [9] Vanderbilt D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys Rev B 1990;41(11):7892–5. [10] Perdew JP, Wang Y. Accurate and simple analytic representation of the electron-gas correlation-energy. Phys Rev B 1992;45(23):13244–9. [11] Ogata S, Li J, Hirosaki N, Shibutani Y, Yip S. Ideal shear strain of metals and ceramics. Phys Rev B 2004;70(10):104104. [12] Krenn CR, Roundy D, Morris Jr JW, Cohen ML. The non-linear elastic behavior and ideal shear strength of Al and Cu. Mater Sci Engng A 2001;317(1–2):44–8. [13] Sˇandera P, Pokluda J. Ideal strength of cubic crystals under triaxial tension. Metal Mater 1994;32(4):180–93.