A. Kinematics of crystalline deformation

The quantitative description of kinematical theory for the mechanics of elastic–plastic deformation of crystals is derived from the earlier work done by Taylor^{Reference Taylor29} and Hill.^{Reference Hill30} It is found that material flows through the crystal lattice via dislocations gliding along the corresponding slip systems. The constitutive models used here follows the earlier works by Asaro,^{Reference Asaro31} Asaro and Rice^{Reference Asaro and Rice32} and Huynh.^{Reference Si, Lu, Huynh, Tieu and Liu33–Reference Liu, Lu, Tieu and Yu35}

The crystalline material under load undergoes crystallographic slip due to dislocation motion on the active slip systems and elastic deformation including stretching and rotating of the crystal lattice,^{Reference Asaro31–Reference Si, Lu, Huynh, Tieu and Liu33,Reference Peirce, Asaro and Needleman36,Reference Asaro and Needleman37} which leads to a hypothesis that the crystallographic slipping may be supposed to firstly occur from the reference configuration to the intermediate configuration, followed by the elastic stretching and rotation from the intermediate configuration to the current configuration. Therefore, three configurations have been proposed to describe the deformation process. The initial configuration corresponds to the original undeformed state of the element and the current configuration corresponds to the deformed state of the element. The intermediate configuration is obtained from the current configuration when relaxing the lattice elastic stretching and rotation. The total deformation gradient F can be decomposed into two different components, namely F* and F^P as

(1)$${\bi{F}} = {{\bfi\partial \bi{x}} \over {\bi\partial X}} = {{\bi{F}}^*} \cdot {{\bi{F}}^{\rm{P}}}\quad ,$$

where X represents the position of material points in the reference configuration and x denotes the position of material points in the current configuration. Here, F is the total deformation gradient, F^∗ describes the elastic deformation gradient including the stretching and rotation of the crystal lattice, and F^P indicates the crystallographic slip on the slip system, which can be written as

(2)$${{\bi{F}}^{\rm{P}}} = \mathop\sum\limits_{\alpha = 1}^N {{{\bi{F}}^{\left( \alpha \right){\rm{P}}}}} = \mathop\sum\limits_{\alpha = 1}^N {{\bi{I}} + {\gamma ^{\left( \alpha \right)}}\left( {{\bf{s}}_0^{\left( \alpha \right)} \otimes {\bf{m}}_0^{\left( \alpha \right)}} \right)} \quad ,$$

where ${{\bi{F}}^{\left( \alpha \right){\rm{P}}}}$ is the contribution of the α-th slip system to F^P, γ^(α) is the shear strain of α-th slip system, ⊗ indicates the dyadic product of two tensors or vectors, I is a second-order unit tensor, and N is the number of activated slip systems. The Green strain tensor can be written as

(3)$${\bi{E}} = {1 \over 2}\left( {{{\bi{F}}^{\rm{T}}}{\bi{F}} - {\bi{I}}} \right)\quad .$$

In the initial configuration, the crystal slip system α consists of a slip direction vector ${\bi{s}}_0^{\left( \alpha \right)}$ lying in a slip plane with normal vector of ${\bi{m}}_0^{\left( \alpha \right)}$. Both ${\bi{s}}_0^{\left( \alpha \right)}$ and ${\bi{m}}_0^{\left( \alpha \right)}$ are unit vectors and should satisfy the orthogonal relationship, namely

(4)$${\bi{s}}_0^{\left( \alpha \right)} \;\cdot\; {\bi{m}}_0^{\left( \alpha \right)} = 0\quad .$$

Both ${\bi{s}}_0^{\left( \alpha \right)}$ and ${\bi{m}}_0^{\left( \alpha \right)}$ will not be changed during crystallographic slip from the initial configuration to the intermediate configuration but converted with the lattice when the lattice is stretched and rotated from the intermediate to the current configuration. The slip direction vector ${\bi{s}}_0^{\left( \alpha \right)}$ can be given in the deformed configuration by

(5)$${{\bi{s}}^{\left( \alpha \right)}} = {{\bi{F}}^{\left( \alpha \right)}} \cdot {\bi{s}}_0^{\left( \alpha \right)} = {{\bi{F}}^*}\left[ {\mathop \sum \limits_{\alpha = 1}^N {\bi{I}} + {\gamma ^{\left( \alpha \right)}}\left( {{\rm\bf{s}}_0^{\left( \alpha \right)} \otimes {\rm\bf{m}}_0^{\left( \alpha \right)}} \right)} \right]{\bi{s}}_0^{\left( \alpha \right)} = {{\bi{F}}^*} \cdot {\bi{s}}_0^{\left( \alpha \right)}\quad ,$$

and the normal to the slip plane $m_0^{\left( \alpha \right)}$ after deformation can be written as

(6)$${\bi{m}^{\left( \alpha \right)}} = {{\bi{s}}^{\left( \alpha \right) - 1}} = {\left( {{\bi{F}^*} \cdot {\bi{s}}_0^{\left( \alpha \right)}} \right)^{ - 1}} = {\bi{m}}_0^{\left( \alpha \right)} \cdot {{\bi{F}}^{* - 1}}\quad .$$

It should be noted that ${s^{\left( \alpha \right)}}$ and ${m^{\left( \alpha \right)}}$ are not in general unit vectors but remain orthogonal as

(7)$${\bi{s}^{\left( \alpha \right)}}\;\cdot\;{\bi{m}^{\left( \alpha \right)}} = 0\quad .$$

In the crystalline coordinate system, for an FCC crystal, ${s^{({\rm{\alpha }})}}$ and ${m^{\left( \alpha \right)}}$ are given by

(8)$${{\bi{s}}^{\left( {1 - 3} \right)}} = {1 \over {\sqrt 3 }}\left( {1,1,1} \right),{{\bi{s}}^{\left( {4 - 6} \right)}} = {1 \over {\sqrt 3 }}\left( { - 1,1,1} \right),{{\bi{s}}^{\left( {7 - 9} \right)}} = {1 \over {\sqrt 3 }}\left( {1, - 1,1} \right),{{\bi{s}}^{\left( {10 - 12} \right)}} = {1 \over {\sqrt 3 }}\left( {1,1, - 1} \right)\quad ,$$

(9)$$\eqalignb{ & {{\bi{m}}^{\left( 1 \right)}} = {1 \over {\sqrt 2 }}\left( {0, - 1,1} \right),{{\bi{m}}^{\left( 2 \right)}} = {1 \over {\sqrt 2 }}\left( {1,0, - 1} \right),{{\bi{m}}^{\left( 3 \right)}} = {1 \over {\sqrt 2 }}\left( { - 1,1,0} \right),{{\bi{m}}^{\left( 4 \right)}} = {1 \over {\sqrt 2 }}\left( {1,0,1} \right)\quad , \cr & {{\bi{m}}^{\left( 5 \right)}} = {1 \over {\sqrt 2 }}\left( {1,1,0} \right),{{\bi{m}}^{\left( 6 \right)}} = {1 \over {\sqrt 2 }}\left( {0, - 1,1} \right),{{\bi{m}}^{\left( 7 \right)}} = {1 \over {\sqrt 2 }}\left( {0,1,1} \right),{{\bi{m}}^{\left( 8 \right)}} = {1 \over {\sqrt 2 }}\left( {1,1,0} \right)\quad , \cr & {{\bi{m}}^{\left( 9 \right)}} = {1 \over {\sqrt 2 }}\left( {1,0, - 1} \right),{{\bi{m}}^{\left( {10} \right)}} = {1 \over {\sqrt 2 }}\left( {0,1,1} \right),{{\bi{m}}^{\left( {11} \right)}} = {1 \over {\sqrt 2 }}\left( {1,0,1} \right),{{\bi{m}}^{\left( {12} \right)}} = {1 \over {\sqrt 2 }}\left( { - 1,1,0} \right)\quad . \cr}$$

The velocity gradient L can be uniquely decomposed into a symmetric part D (stretching rate tensor) and a skew-symmetric Ω (spin tensor) part as

(10)$${\bi{L}} = {\bi{D}} + \iOmega \quad ,$$

(11)$${\bi{D}} = {1 \over 2}\left( {\bi{L} + {\bi{L}^{\rm{T}}}} \right)\quad ,$$

(12)$$\iOmega = {1 \over 2}\left( {\bi{L} - {\bi{L}^{\rm{T}}}} \right)\quad ,$$

where D is, in general, called the rate of deformation. Both tensors D and Ω can also be decomposed into the elastic stretching and lattice rotation component (D* and Ω*) and the plastic component (D^P and Ω ^P) as follows,

(13)$${\bi{D}} = {\bi{D}^*} + {\bi{D}^{\rm{P}}}\quad ,$$

(14)$$\iOmega = {\iOmega ^*} + {\iOmega ^{\it{P}}}\quad ,$$

and the tensor D^P in Eq. (13) can be expressed as

(15)$${\bi{D}^{\rm{P}}} = {1 \over 2}\left( {{\bi{L}^{\rm{P}}} + {\bi{L}^{{{\rm{P}}^{\rm{T}}}}}} \right) = \mathop\sum\limits_{\alpha = 1}^N {{\bi{p}^{\left( \alpha \right)}}{{\dot \igamma }^{\left( \ialpha \right)}}} \quad ,$$

where p(α) can be defined as

(16)$${\bi{p}^{\left( \alpha \right)}} = {1 \over 2}\left( {{\bi{s}^{\left( \alpha \right)}} \otimes {\bi{m}^{\left( \alpha \right)}} + {\bi{m}^{\left( \alpha \right)}} \otimes {\bi{s}^{\left( \alpha \right)}}} \right)\quad .$$

Similar to D^P, the tensor Ω ^P in Eq. (14) can be expressed as

(17)$${\iOmega ^{\it{P}}} = \mathop\sum\limits_{\alpha = 1}^N {{\bi{w}^{\left( \alpha \right)}}{{\dot \igamma }^{\left( \ialpha \right)}}} \quad ,$$

where q^(α) can be defined as

(18)$${\bi{w}^{\left( \alpha \right)}} = {1 \over 2}\left( {{\bi{s}^{\left( \alpha \right)}} \otimes {\bi{m}^{\left( \alpha \right)}} - {\bi{m}^{\left( \alpha \right)}} \otimes {\bi{s}^{\left( \alpha \right)}}} \right)\quad ,$$

where ${{\bi{s}}^{\left( \alpha \right)}} \otimes {{\bi{m}}^{\left( \alpha \right)}}$ is the Schmid factor. Here, p^(α) and w^(α) are symmetric and asymmetric part of the Schmid factor, respectively.

B. Rate-dependent hardening model

The crystal plasticity constitutive model used in this study follows the line described by Asaro^{Reference Asaro31} and it is incorporated into the implicit finite element code ABAQUS/Standard through the user material subroutine (UMAT). The prominent functions of UMAT are to provide the material with the Jacobian matrix for the constitutive model, and to update the stresses and the solution dependent state variables. In this study, we adopt the UMAT framework initially developed by Huang^{Reference Huang38} and use the formulations established by Bassani and Wu^{Reference Bassani and Wu39} as the hardening model, which has been specifically described by Liu et al. and Lu et al.^{Reference Liu, Lu, Tieu and Yu35,Reference Lu, Deng, Tieu, Su, Zhu and Liu40} The formulations of the rate-dependent hardening model used in the present study can be expressed as

(19)$${\dot \igamma ^{\left( \ialpha \right)}} = \dot \gamma _0^{\left( \alpha \right)}{\mathop{\rm sgn}} \left( {{\itau ^{\left( \ialpha \right)}}} \right){\left| {{{{\itau ^{\left( \ialpha \right)}}} \over {\itau _{\it{c}}^{\left( \ialpha \right)}}}} \right|^{\rm{n}}}\quad ,$$

with

(20)$${\mathop{\rm sgn}} \left( x \right) = \left\{ {\matrix{ {\! - 1,\quad x\; \lt\; 1} \cr \;\;{1,\quad x \;\ge\; 1} \cr} } \right.\quad ,$$

where $\dot \gamma _0^{\left( \alpha \right)}$ is the reference (initial) value of the shear strain rate, which is a constant for all the slip systems. Here, n is the strain rate sensitive exponent. Both $\dot \gamma _0^{\left( \alpha \right)}$ and n are the material parameters. Here, $\tau _{\rm{c}}^{\left( \alpha \right)}$ is the critical resolved shear stress of the slip system α, which represents the strength of the material.

The rate of change of the critical resolved shear stress is expressed by^{Reference Hill30}

(21)$${\dot \itau _{\it{c}}}^{\left( \ialpha \right)} = \mathop\sum\limits_{\beta = 1}^N {{f_{\alpha \beta }}{h_{\alpha \beta }}{{\dot \igamma }^{\left( \ibeta \right)}}} \quad ,$$

(22a)$${h_{\alpha \alpha }} = \left[ {\left( {{{\rm{h}}_0} - {{\rm{h}}_{\rm{s}}}} \right)\sec {h^2}\left( {{{\left( {{h_0} - {h_{\rm{s}}}} \right){\igamma ^{\left( \ialpha \right)}}} \over {{\tau _1} - {\tau _0}}}} \right) + {h_{\rm{s}}}} \right]\left[ {1 + \mathop \sum \limits_{\matrix{ {\beta = 1} \cr {\beta \ne \alpha } \cr} }^N {f_{\alpha \beta }}\tanh \left( {{{{\igamma ^{\left( \ibeta \right)}}} \over {{\gamma _0}}}} \right)} \right],\quad {\rm{for}}\,\alpha = \beta \quad ,$$

(22b)$${h_{{\rm{\alpha \beta }}}} = q{h_{\alpha \alpha }},\quad {\rm{for}}\;\alpha \ne \beta \quad ,$$

where h _αβ is the hardening modulus including the self-hardening of each system (α = β) and latent hardening (α ≠ β), the constant q is the latent hardening parameter, the constant γ₀ is the reference value of slip, γ is the shear strain, τ₀ is the initial critical resolved shear stress, τ₁ is the breakthrough stress where large plastic flow initiates, h ₀ is the hardening modulus just after initial yield, h _s is the hardening modulus during easy glide, and f _αβ is the magnitude of the strength of a particular slip interaction between two slip systems α and β. The factor f _αβ depends on the geometric relation between two slip systems. There are five constants for f _αβ, namely α₁ (no junction), α₂ (Hirth lock), α₃ (coplanar junction), α₄ (glissile junction), and α₅ (sessile junction).

C. Finite element implementation

The aforementioned crystal plasticity constitutive model is implemented into the implicit finite element code ABAQUS/Standard by using the UMAT which are able to provide the material Jacobian matrix ${{\partial \Delta \sigma } \mathord{\left{\vphantom {{\partial \Delta \sigma } {\partial \Delta \varepsilon }}} \right\kern-\nulldelimiterspace}/ {\partial \Delta \varepsilon }}$, for the constitutive model and to update the stresses and the solution dependent state variables.

The commercial software Abaqus6.9, Pro/ENGINEER Wildfire 4.0, and CUBIT 13.1 have been used to simulate the deformation procedure of nanoindentation. A 3D model has been built to describe the behavior of nanoindentation, which is shown as Fig. 1. The Voronoi diagram is the most popular method to generate polycrystalline material structure^{Reference Zhu, Thorpe and Windle41} in the present. In this study, the 3D Voronoi diagram has been used to generate two 3D cells. Firstly, the input file of 3D nanoindentation setup was generated by Abaqus6.9. Secondly, the input file was implemented into a MATLAB program to construct two 3D cells (grains). Each Voronoi cell corresponds to one seed and the number of seeds will be controlled to determine the average size of the cells. The generated two Voronoi cells are then assigned with different crystallographic orientations and implemented into the CPFEM model. In the CPFEM simulation, each Voronoi cell represents a virtual grain. The detailed description for constructing poly-crystal CPFEM model has been given in our earlier work.^{Reference Huynh34}

FIG. 1.

3D nanoindentation model setup: (a) the whole model; (b) the upper surface of Cu single crystal; (c) the upper surface of Cu bicrystals (or Cu–Al bicrystals) and (d) through thickness view of Cu bicrystals (or Cu–Al bicrystals).

The indenter is a Berkovich indenter with a 200 nm radius round tip. 13,024 eight-node brick elements and 14,463 nodes with reduced integration (element id: C3D8R) are used in the CPFEM model. A refined mesh is generated in the contact area [Figs. 1(b) and 1(c)] directly underneath the indenter to obtain the accurate contact solution while a coarser mesh is created in the rest region to decrease the elements number in the model and to keep the computer time within limits. The importance of having an appropriate mesh density in the contact area has been proposed by Liu et al.^{Reference Liu, Varghese, Ma, Yoshino, Lu and Komanduri20,Reference Liu, Wang, Yoshino, Roy, Lu and Komanduri23} The size of the smallest element is about 0.1 μm in all three directions.

The X, Y, and Z coordinates represent the rolling direction, transverse direction, and normal direction, respectively. In this model, the specific Goss orientation is set as the initial orientations, namely the (011) slip plane parallels to the Z axis, and the [100] slip direction parallels to the X axis, to validate the model by comparing the numerical result with the experimental observation. Figure 1(c) shows the upper surface of the sample for Cu bicrystal or Cu–Al bicrystal. For Cu bicrystals, the orientations of Grain 2 are always set as Goss orientations while the orientations of Grain 1 are set as Cube, Brass and R, respectively. For a Cu–Al bicrystal, Grain 2 is set as Cu and its orientations are always set as Goss orientations. Grain 1 is set as Al and its orientations are set as Cube, Brass and R, respectively. The misorientation between Grain 1 and Grain 2 for both Cu bicrystals and Cu–Al bicrystals is listed in Table I. Clearly, all the GBs between Grain 1 and Grain 2 are high angle GBs.

TABLE I.

Misorientations for Cu bicrystal and Cu–Al bicrystal between Grain 1 and Grain 2.

The dimensions of the sample in the FEM model are 60 µm × 60 µm × 30 µm. The height of the sample is much larger than the maximum indentation depth (1 µm) to avoid the influence from the sample.^{Reference Fischer-Cripps42} In addition, the length and the width are bigger enough as the stress contour never reaches four surrounding surfaces of the sample.

All nodes on the bottom surface and four surrounding surfaces of sample are constrained along three axes considering the real indentation condition (the metal sample was mounted inside a piece of resin and just the top surface is free). Liu et al.^{Reference Liu, Wang, Yoshino, Roy, Lu and Komanduri23} found that both the indentation depth and the load–displacement curve are not affected by the coefficient of friction. Therefore, a frictionless contact pair is defined by two contact surfaces with associated nodes between the indenter and sample. The time step increment is set for the convergence of modeling. In this study, a fixed time step increment 0.01 s is chosen and a total time step increments of 24,062 in the simulation are used, including contact, loading, and unloading.

According to the report from Bassani et al.^{Reference Wu, Bassani and Laird43} and Huang,^{Reference Huang38} the factor f _αβ for single-crystal copper is taken as: α₁ = α₂ = α₃ = 8, α₄ = 15 and α₅ = 20 in this study. Three elastic moduli can be used as C₁₁ = 168,400 MPa, C₁₂ = 121,400 MPa, and C₄₄ = 75,400 MPa, respectively.^{Reference Wu, Bassani and Laird43} Other material parameters in the hardening models (Eqs. 19–22) are listed in Table II.^{Reference Wu, Bassani and Laird43,Reference Liu, Lu and Tieu44} Franciosi et al.^{Reference Franciosi, Berveiller and Zaoui45} and Lu et al.^{Reference Lu, Deng, Tieu, Su, Zhu and Liu40} have reported that the factor f _αβ for aluminum can be chosen as: α₁ = α₂ = α₃ = 1.75, α₄ = 2 and α₅ = 2.25. Other material parameters in the hardening models are listed in the Table III.^{Reference Liu, Lu, Tieu and Yu35,Reference Lu, Deng, Tieu, Su, Zhu and Liu40} All of these parameters were evaluated by fitting the simulated stress–strain curve with the experimental results of single-crystal aluminum under plane strain compression.^{Reference Liu, Maurice, Driver and Hansen46,Reference Akef, Fortunier, Driver and Watanabe47} Aluminum has an FCC structure with elastic moduli C₁₁ = 112,000 MPa, C₁₂ = 66,000 MPa, and C₄₄ = 28,000 MPa. There is only one set of slip systems for FCC metals, {111}$\left\langle {110} \right\rangle$, which includes 12 different slip systems shown in the Table IV. In the deformed single-crystal copper and aluminum, slips occur on the {111} planes and in the $\left\langle {110} \right\rangle$ directions. All of aforementioned parameters for copper were evaluated by both experiments and simulations^{Reference Huang38,Reference Casals, Ocenasek and Alcala48,Reference Alcala, Casals and Ocenasek49} and all parameters for aluminum were validated in the simulations of nanoindentation, rolling, tensile, and ECAP deformation.^{Reference Si, Lu, Huynh, Tieu and Liu33,Reference Liu, Lu, Tieu and Yu35,Reference Lu, Deng, Tieu, Su, Zhu and Liu40,Reference Huynh, Lu, Si and Tieu50,Reference Liu, Tieu, Lu, Zhu and Deng51}

TABLE II.

Parameters in the constitutive model for copper.^{Reference Wu, Bassani and Laird43,Reference Liu, Lu and Tieu44}

TABLE III.

Parameters in the constitutive model for aluminum.^{Reference Liu, Lu, Tieu and Yu35,Reference Lu, Deng, Tieu, Su, Zhu and Liu40}

TABLE IV.

Notation of the slip systems for the FCC materials considered in this study.