A Comparative Study of Indentation Size Effect Models for Different Materials

1. Introduction

In traditional mechanical test, standard tensile specimen needs to be cut from materials or structure for obtaining the mechanical properties of materials. However, for service equipment, small size materials and expensive materials in aerospace, nuclear power, microelectronics, medicine and other fields, it is difficult to obtain the mechanical properties parameters by the traditional mechanical test. Due to the “no sample” and “micro damage”, the indentation test has been proposed by Oliver and Pharr, and widely employed to calculate the elasticity modulus and flow stress of different materials [1,2,3,4,5]. The phenomenon of indentation size effect (ISE), which the measured hardness changes with indentation depth, have been observed in various materials such as metallic, ceramics, fibers and thin films [6,7,8,9,10].

Over the past several decades, several theories had been developed to explain the phenomenon of ISE. The major explanations of ISE can be classified into three branches: the energy balance theory, the strain gradient plasticity theory and the dislocation theory [11]. In addition, other factors can be regarded to cause the ISE, which include friction between the surface and the indenter [12], variation of contact surface [13], presence of residual surface stress [14], energy dissipation associated with the contact surfaces [15], variation in solid load bearing contact area inside the indentation region [16], the role of elastic deformation [17], and pre-existing defects [18].

Based on above explanations and experimental data, many ISE models have been proposed to describe the ISE. An empirical power function of applied load and indentation diagonal is introduced by Meyer. According to Hays and Kendall [19], there is a threshold minimum load that must be surpassed for plastic deformation to take place. Li and Bradt [20] formulated the proportional specimen resistance (PSR) model by taking into account energy balance. Gong et al. [14] have proposed a modified proportional specimen resistance (MPSR) model, taking into account the impact of residual compressive or tensile stress on hardness measurements. An elastic plastic deformation (EPD) model, which was introduced by Bull et al. [21], suggested that the elastic stresses at the indentation edge are a direct result of the discrete plastic deformation. Ruiz-Moreno et al. [22] proposed an empirical model based on exponential function. A model, which taking into account the contribution of grain size and dislocation pinning defects, was developed by Hou and Jennet [23].

One of the most popular models of ISE is the Nix-Gao model which considers the variation of geometrically necessary dislocations (GND) density with indentation depth [24]. However, further nanoindentation analysis revealed that the Nix-Gao model tends to overestimate hardness in cases of extremely shallow indentations, due to the intense repulsive force between GNDs [25]. Therefore, Swadener et al. [25]., Feng and Nix [26], and Durst et al. [27] thought that the Nix-Gao model should be modified by considering that the ratio of the size of plastic zone to indentation depth is not fixed. A modified version of the Nix-Gao model was developed by Huang et al. [28] for assessing nano-indentation hardness, taking into account a maximum permissible GND density. Feng-Nix [26] and Haušild [29] had incorporated the exponentially decay function and negative power function into the Nix-Gao model to depict the connection between plastic zone size and indentation depth, respectively. In the meantime, Abu Al-Rub [30] enhanced the Nix-Gao model by incorporating the Taylor hardening law, which takes into account nonlinearity of flow stresses caused by both SSDs and GNDs. Qiu et al. [31] improved the Nix-Gao model by introducing the friction stress that was not related to dislocation activities. Different with above modification, Yuan and Chen [32] improved the Nix-Gao model by a gradient plasticity model with an intrinsic micro-material length. Liu et al. [17] had integrated elastic deformation into the Nix-Gao model to enhance its applicability for predicting ISE at depths below 100 nm.

Above models have contributed to comprehensively understand the ISE. These models often involve different form and number of parameters that can make our wonder which is the best in existed ISE models. So far, Nix-Gao model and four modification of it are evaluated by Haušild [29]. Renjo et al. [1], Maiti et al. [33], Alao et al. [34], and Kathavate et al. [35] used part of models to predict the ISE of specific material. And others only compared Nix-Gao model with one improved by themselves [26,30,31]. However, to the best of our knowledge, a comprehensive and quantitative evaluation of existed ISE modes is still lacked for various materials. In this paper, existed ISE modes are reviewed. The ISE data of different materials are assembled from publications. These reviewed models are then applied to predict the ISE of assembled data. After presenting the predicted results and evaluating the accuracy of each model for different materials, the suitable models are recommended to describe the specific ISE of different materials.

2. Indentation Size Effect (ISE) Models

In this section, sixteen ISE models of isotropic materials will be reviewed in the follow. Readers should note that the ISE models of multilayers and bicrystals materials are not covered in this article [36,37,38,39,40].

2.1. Meyer Model

An empirical power function of applied load and indentation diagonal is introduced by Meyer as follow:

$$P={\alpha }_1{d}^m$$

(1)

where P is the maximum load applied during the indention test, d is the indentation length, m and α₁ are Meyer’s index and material constant, respectively.

For a given indenter, a linear relation is assumed between indentation length d and indentation depth h as follow:

$$h={\alpha }_{2}d$$

(2)

If taking no account of tip bluntness, the contact area A is proportional to the contact area of the indenter d², which can be expressed as:

$$A={\alpha }_3{d}^2=\frac{{\alpha }_3}{{\alpha }_2^2}{h}^2$$

(3)

where α₂ and α₃ are the constants that depend on indenter geometry.

Combining Equations (1)–(3) and hardness H = P/A, the Meyer model is obtained as:

$$H=m_1{h}^{{-m}_0}$$

(4)

where m₁ = α₁/[α₃(α₂)^n-2] ≥ 0 and m₀ = 2-m ≥ 0 are constants depends on material and indenter geometry.

2.2. Hays and Kendall (H-K) Model

According to Hays and Kendall [2], plastic deformation can only occur after a minimum load has been surpassed. The relationship between the applied load P and the indentation diagonal d is described as follows:

$$P=k_0+k_1{d}^2$$

(5)

where k₀ is required load of producing the initial plastic deformation, and k₁ is the load independent constant.

Combining Equations (1)–(3) and hardness H = P/A, the H-K model is obtained as:

$$H=k_2+\frac{k_3}{h^2}$$

(6)

where k₂ = k₁/α₃ ≥ 0 and k₃ = k₀(α₂)²/α₃ ≥ 0.

2.3. Proportional Specimen Resistance (PSR) Model

Taking into consideration energy balance, Li and Bradt [20] given a correlation between the applied load P and indentation diagonal d as

$$P=l_{1}d+l_2{d}^2$$

(7)

where l₁ and l₂ are the elastic and friction resistance at indenter/specimen interface and plastic deformation constant, respectively.

Combining Equations (1)–(3) and hardness H = P/A, the PSR model is determined as:

$$H=b_0+\frac{b_1}{h}$$

(8)

where b₀ = l₂/α₃ ≥ 0 and b₁ = l₁α₂/α₃ ≥ 0.

2.4. Modified Proportional Specimen Resistance (MPSR) Model

Taking into account the impact of residual compressive or tensile stress on hardness measurements, Gong et al. [14] had proposed a modified proportional specimen resistance (MPSR) model as follow:

$$P=a_0+a_{1}d+a_2{d}^2$$

(9)

where l₀ is the induced-residual stress effect during surface preparations, l₁ and l₂ are described previously.

Combining Equations (1)–(3) and hardness H = P/A, the MPSR model is determined as:

$$H=g_0+\frac{g_1}{h}+\frac{g_2}{h^2}$$

(10)

where g₂ = a₀(α₂)²/α₃, g₀ = a₂/α₃ ≥ 0 and g₁ = a₁α₂/α₃ are constants.

2.5. Elastic Plastic Deformation (EPD) Model

Assuming that the elastic stresses at the indentation edge are a direct result of the discrete plastic deformation, the EPD model is introduced by Bull et al. [21] as:

$$H={\left(e_0+\frac{e_1}{h}\right)}^2$$

(11)

where e₀ ≥ 0 and e₁ ≥ 0 are constants depends on material and indenter geometry.

2.6. Exponential (EXP) Model

An empirical model, which based on exponential function, is proposed by Ruiz-Moreno et al. [22] as follow:

$$H=r_1\mathit{exp}\left(-r_{2}h\right)+r_3$$

(12)

where r₁, r₂ and r₃ are constants.

2.7. Hou-Jennet(H-J) Model

Considering the contribution of grain size and dislocation pinning defects on the flow stress (hardness) of a material, a model was developed by Hou and Jennet [23] as follow:

$$H=j_0+\sqrt{\frac{j_1}{h}+j_2}$$

(13)

where j₂ = j₃/D+j₄r ≥ 0, D and r are grain size and density of dislocations respectively, j₀ ≥ 0, j₁ ≥ 0, j₃ ≥ 0 and j₄ ≥ 0 are material constants.

2.8. Nix-Gao (N-G) Model

Taking the effect densities of geometrically necessary and statistically stored dislocations on hardness into consideration, Nix–Gao model [24] is proposed by using the strain gradient plasticity theory to describe the ISE of crystalline materials, then the relation of indentation hardness and depth is given as follow:

$$H=H_0\sqrt{1+\frac{h_0}{h}}$$

(14)

where H₀ ≥ 0 is the classical indentation hardness which depends on density of statistically stored dislocations, indenter’s geometry, magnitude of the Burgers vector and material’s shear modulus in the deformation zone, h₀ ≥ 0 is a length scale characterizing the depth dependence of hardness.

2.9. Nix-Gao-Feng (N-G-Fs) Model

Exponentially decay function that described the relationship between the size of plastic zone and indentation depth, had been used to modify the Nix-Gao model by Feng-Nix [26], then the relation of indentation hardness and depth is given as follow:

$$H=H_0\sqrt{1+{\left(1+r_0\exp \left(-\frac{h}{h_1}\right)\right)}^{-3}\frac{h_0}{h}}$$

(15)

where r₀ ≥ 0 and h₁ ≥ 0 are fitting parameters.

2.10. Nix-Gao-Huangs (N-G-Hs) Model

A modified version of the Nix-Gao model was developed by Huang et al. [28] for assessing nano-indentation hardness, taking into account a maximum permissible GND density. The model given a piece-wise relation between indentation hardness and depth as follow:

$$H=\left\{\begin{array}{cc}{H}_0\sqrt{1+\frac{2h_0}{3h_{nano}}}& h<h_{nano}\\ H_0\sqrt{1+\frac{h_0}{h}-\frac{h_0{h}_{nano}^2}{3h^3}}& h\ge h_{nano}\end{array}\right.$$

(16)

where h_nano ≥ 0 is a new characteristic length for nano-indentation which depends on indenter’s geometry, magnitude of the Burgers vector and maximum allowable GND density in the deformation zone.

2.11. Nix-Gao-Haušild (N-G-H) Model

Haušild [29] assumes that the relation between number of dislocations contained in the effective plastic zone under the indenter and the penetration depth is a power function, then Nix-Gao model is modified as follow:

$$H=H_0\sqrt{1+\left(1-\exp \left(-\frac{h^n}{h_1}\right)\right)\frac{h_0}{h}}$$

(17)

where h₁ ≥ 0 is a parameter (of the similar meaning as in N-G-Hs model) and n is an exponent (shape parameter).

2.12. Nix–Gao-Abu Al-Rub (N-G-A) Model

Abu Al-Rub [30] enhanced the Nix-Gao model by incorporating the Taylor hardening law, which takes into account nonlinearity of flow stresses caused by both SSDs and GNDs. The relation of indentation hardness and depth of this model is given as follow:

$$H=H_y+\left(H_0-H_y\right){\left(1+{\left(\frac{h_0}{h}\right)}^{\frac{\beta }{2}}\right)}^{\frac{1}{\beta }}$$

(18)

where H_y ≥ 0 is the hardness related to the initial yield stress, and β ≥ 0 is the interaction coefficient.

2.13. Nix–Gao-Qius (N-G-Qs) Model

Introducing the intrinsic lattice resistance, the Nix-Gao model is improved by Qiu et al. [31] as follow:

$$H=3{\sigma }_0+H_0\sqrt{{\left(1-\frac{3{\sigma }_0}{H_0}\right)}^2+\frac{h_0}{h}}$$

(19)

where σ₀ is the contribution of friction stress to the tension flow stress.

2.14. Nix–Gao-Yuan-Chen (N-G-YC) Model

Yuan and Chen [32] improvs the Nix-Gao model by a gradient plasticity model with an intrinsic micro-material length, and the relation of indentation hardness and depth is expressed as follow:

$$H=H_0\sqrt{1+\frac{h_0}{h}-{\left(\frac{h_3}{h}\right)}^2}$$

(20)

where h₃ ≥ 0 is the intrinsic micro-material length, and depends on the material property and the indenter geometry.

2.15. Nix–Gao-Lius (N-G-Ls) Model

Liu et al. [17] had integrated elastic deformation into the Nix-Gao model to enhance its applicability for predicting ISE at depths below 100 nm, and this model is expressed as follow:

$$H=H_0\sqrt{1+\frac{h_0}{f^3{h}^2}\left[\frac{h}{\pi }-\left(\frac{\pi -2}{\pi }\right)\frac{{\left(h+\xi \right)}^m}{\left(c+dh\right)}\right]-\frac{h_0}{f^3{h}^3}\frac{{\left(h+\xi \right)}^m}{\left(c+dh\right)}\left[\frac{h}{\pi }-\left(\frac{\pi -2}{\pi }\right)\frac{{\left(h+\xi \right)}^m}{\left(c+dh\right)}\right]}$$

(21)

where f is the ratio of GND storage region radius to the indentation radius, c and d are associated with load-indentation depth curve, ξ and m are geometric parameters accounting for the indenter shape, m = 2 for a conical tip and m = 1.5 for a spherical tip.

2.16. Li-Zhang -Hans (L-Z-Hs) Model

A mechanistic model incorporating surface undulation and indenter tip irregularity, which can simultaneously describe Normal and Reverse ISE, is proposed by Li et al. as follow [41]:

$$H=H_0{\left(\frac{h}{h+\delta }\right)}^m$$

(22)

where δ is the deviation of indentation depth, and caused by surface undulation and indenter tip irregularity. m = 2 for a conical tip and m = 1.5 for a spherical tip.

3. Collection of Experimental and Simulated Results

In recent decades, many researches have consistently shown that isotropic materials exhibit ISE when subjected to nanoindentation [17,26,29,41,42,43,44,45,46,47]. Three main types of ISEs are observed from experimental and simulated results as follow:

Normal ISE: In Figure 1, the hardness H decreases with the increase of indentation depth h, and then tend to a constant.

Reverse ISE: In Figure 1, the hardness H increases with the increasing indentation depth h and then tends to a constant.

Transition of normal to reverse ISE: In Figure 1, the hardness H first increases, then decreases and finally tend to a constant with the increasing indentation depth h.

Readers should note that the ISE of multilayers and bicrystals materials, which the hardness H first decreases and then increases, then decreases, and finally tends to a constant with the increasing indentation depth h, is not covered in this article [36,37,38,39,40].

To evaluate the accuracy and rationality of the sixteen ISE models in section 2, nine test data of MgO single crystals, 16MND5 steel, High entropy alloy (HEA), Ni Carbide Silicon (NiCSi), TC4 titanium alloy, Pulsed electro-deposited Ni (PED Ni), ZrO₂ ceramic, Cu single crystals, Y₂O₃-ZrO₂ ceramic, are collected from the literature. The collected data sets are provided in tabular form and given in Table 1. It is crucial to ensure a sufficient number of experimental or simulated data points to avoid overfitting and accurately determine the model parameters. Otherwise, it becomes difficult to distinguish the accuracy of different ISE models. Considering for the collected data is very time-consuming; these data are uploaded in ResearchGate (DOI: 10.13140/RG.2.2.26190.68167) for use by other researchers.

4. Comparative Study of Indentation Size Effect Models

In the sixteen ISE models, the model parameters can be determined by utilizing an equivalent amount of experimental or simulated data. However, more experimental or simulation data are preferred to obtain high-accuracy ISE curves. Thus, an optimization method is utilized to calculate the model parameters of sixteen ISE models. For three types of ISE data in Table 1, the parameters of sixteen ISE models are determined by optimization method and given in Table 2.

4.1. Qualitative Comparison

By inserting these optimized parameters in Table 2 into the sixteen ISE models, predicted ISE curves can be obtained. The experimental or simulated data and the predicted ISE curves of sixteen ISE models are shown in are shown in Figure 2, Figure 3 and Figure 4.

According to Figure 2, the Meyer, H-K, PSR, MPSR and EPD give a poor prediction to the experimental or simulated data of the normal ISE. This is because Meyer could not predict the finite asymptotic indentation hardness when the indentation is deep, and H-K, PSR, MPSR and EPD cannot reflect the nonlinear variation of the normal ISE. The empirical EXP model is nonlinear, and can predict a non-vanishing asymptotic indentation hardness. However, EXP model underestimates the indention hardness for the relatively lower indentation depth and overestimates the indentation hardness for the relatively higher indentation depth. This means that the EXP model cannot reflect the intrinsic nonlinear characteristics of the normal ISE.

It can be found that N-G, N-G-Hs, N-G-YC and N-G-Ls models can accurately describe the normal ISE of HEA [see in Figure 2e,f], but it gives a poor fit to the normal ISE of MgO and 16MND5 [see in Figure 2a–d]. This is because the indentation depth of MgO and 16MND5 is at the nanometer and micron scale, while the indentation depth of HEA is mainly at the micron scale. The assumption of N-G, N-G-Hs, N-G-A, N-G-YC and N-G-Ls models, which the radius of GNDs zone is equal to the contact radius, breaks down when the indentation depth transits from micron to nanometer scale [29,48].

Figure 2 shows that predicted ISE curves of H-J, N-G-A and N-G-Qs models agrees well with the experimental or simulated data of the normal ISE. The reason is that friction stress that is not related to dislocation activities takes into account has been introduced into these models. It can be found from Equations (13), (18) and (19) that the H-J and N-G-Qs models have the same expression, and the H-J and N-G-Qs models are obtained when β = 2 in N-G-A [Equation (18)]. Moreover, N-G-Fs and N-G-H models also can accurately predict the experimental or simulated data of the normal ISE [see in Figure 2]. This is because N-G-Fs model takes the progressive expansion of the plastic zone with decreasing indentation depth into account, and N-G-H model suppose that the number of dislocations in effective plastic zone increases with the increasing hⁿ in Equation (18) [49].

Figure 3 shows the compared experimental or simulated results of different materials and predicted curves of sixteen ISE models for reverse ISE. For NiCSi [see in Figure 3a,b], the EXP model performs well, but the MPSR and N-G-A models can only reflect the trend of reverse ISE. For the TC4 [see in Figure 3c,d], the MPSR, EXP, N-G-Fs, N-G-H, N-G-A and L-Z-Hs models predict accurately, while N-G-Hs and N-G-YC models can only describe its growth trend. For the PED Ni [see in Figure 3e,f], the MPSR, N-G-H, N-G-A, N-G-Ls and L-Z-Hs models well predicts the reverse ISE, whereas the EXP and N-G-YC models can only embody the trend of reverse ISE. Other models fail to predict the reverse ISE of NiCSi, TC4 and PED Ni.

For the transition of normal to reverse ISE, the prediction of the sixteen ISE models is illustrated in Figure 4. The MPSR, N-G-Fs N-G-Hs, N-G-H and N-G-Ls models give good prediction of the test data of ZrO₂ and Cu [see in Figure 3a–d]. The experimental data of Y₂O₃-ZrO₂ can be well described by the N-G-Fs and N-G-H models [see in Figure 4e,f].

In fact, the sixteen ISE models are proposed to describe the normal ISE of different materials. However, to our surprise, some of these models are able to capture the reverse ISE and the transition of normal to reverse ISE of different materials. Moreover, it can be concluded from the prediction results of different models that the friction stress that is not related to dislocation activities and the Indentation size effect of plastic zone are critical to accurately predict the three types of ISEs. Therefore, these two factors should be considered in future models.

4.2. Quantitative Comparison

In order to quantificationally evaluate the accuracy of the sixteen ISE models, the determination coefficients (DC) are defined as follow:

$$\mathrm{DC}=1-\frac{{\sum }_{i=1}^k{\left(H_{\mathrm{ri}}-H_{\mathrm{p}\mathrm{i}}\right)}^2}{{\sum }_{i=1}^k{\left(H_{\mathrm{ri}}-H_{\mathrm{rm}}\right)}^2}$$

(23)

where k is the number of experimental or simulated data for a specific material, i = 1; 2. . . k. H_ri is the i-th experimental or simulated data, and H_pi is the i-th predicted value according to ISE models. H_rm is the mean value of the experimental or simulated data, which is defined as follow:

$$H_{\mathrm{rm}}=\frac{1}{k}{\sum }_{i=1}^k{H}_{\mathrm{ri}}$$

(24)

The DC is employed to determine the reliability of prediction by the sixteen ISE models. Having a higher DC is desirable, the best outcome is when the experimental or simulated data and predicted data match perfectly, resulting in zero misfit (DC = 1).

Substituting optimized parameters in Table 2 into the sixteen ISE models, the predicted value H_p can be calculated for each experimental or simulated data H_r. Then by using the Equations (23) and (24), the DCs of the sixteen ISE models are determined for different materials and given in Table 3.

Among the sixteen ISE models, N-G-Fs model has the best fits to the normal ISE of MgO and 16MND5 with the maximum DC values. The N-G-H and N-G-Ls models give the maximum DC for HEA. Moreover, it can be found that other ISE models except the H-K, PSR, EPD, N-G and L-Z-Hs provide higher DC values for the normal ISE.

For reverse ISE in Table 3, the highest DC value for NiCSi is given by the EXP model, while the N-G-Fs model provides the maximum DC value for TC4 and the N-G-A model forecasts the maximum DC value for PED Ni. Moreover, when the indentation depth is very small, it can be noted in Table 3 that N-G-YC mode fits the data of NiCSi and PED Ni with high DC values (0.9314 and 0.9835), but it predicts exactly opposite trend compared to the experimental data [see in Figure 3a,b,e,f]. This means that the DC alone cannot fully assess the accuracy of models’ predictions, and the comparison of test data with models’ predictions should be analyzed.

For the transition of normal to reverse ISE in Table 3, the N-G-YC model has the highest value of DC for ZrO₂. One can see that the maximum DC of Cu is obtained by using the N-G-Fs and N-G-H models. The maximum DC of Y₂O₃-ZrO₂ determined by the N-G-H model is 0.9206.

In addition, the mean DCs of the sixteen ISE models are calculated and given at the right of Table 3. It can be seen that the mean DC=0.92060 of the N-G-Fs model is the maximum among the sixteen ISE models, followed by the N-G-H model. It is obvious that a better correlation is obtained by using the N-G-Fs and N-G-H models, but they cannot accurately predict the reverse ISE. Therefore, it can be concluded that the N-G-Fs and N-G-H models should be modified to accurately predict the reverse ISE.

5. Conclusions

In this paper, through applying the sixteen ISE models to predict the three types of ISEs, namely, normal ISE, reverse ISE and transition of normal to reverse ISE. The following conclusions can be drawn.

For the normal ISE, it can be found that the predicted ISE curves of H-J, N-G-A and N-G-Qs models agrees well with the experimental or simulated data. The reason is that friction stress that is not related to dislocation activities takes into account has been introduced into these models. N-G-Fs and N-G-H models also can accurately predict the normal ISE. The reason for this is that the N-G-Fs model considers the gradual expansion of the plastic zone as indentation depth decreases, while the N-G-H model assumes an increase in the number of dislocations in the effective plastic zone as the hⁿ increases. On this basis, it can be concluded that the friction stress that is not related to dislocation activities and the indentation size effect of plastic zone are critical to accurately predict the three types of ISEs. Therefore, these two factors should be considered in future models.

The sixteen ISE models are originally proposed to describe the normal ISE of different materials. However, to our surprise, some of these models are able to capture the reverse ISE and the transition of normal to reverse ISE of different materials. The DCs of the sixteen ISE models are also determined for different materials. For reverse ISE, the highest DC value for NiCSi, TC4 and PED Ni are given by the EXP, N-G-Fs and N-G-A models, respectively. For the transition of normal to reverse ISE, the N-G-YC, N-G-Fs, and N-G-H models produce the maximum DC for ZrO₂, Cu and Y₂O₃-ZrO₂, respectively. Moreover, the mean DC of the N-G-Fs model is the maximum among the sixteen ISE models, followed by the N-G-H model, but they cannot accurately predict the reverse ISE. Therefore, it can be concluded that the N-G-Fs and N-G-H models should be further modified to accurately predict the reverse ISE.