Examining the Effect of Pileup on the Accuracy of Sharp Indentation Testing

One essential problem related to instrumented indentation is the effect of pileup, which could introduce significant errors on the measured hardness and elastic modulus. In this work, we have assessed some critical issues associated with instrumented indentation by means of numerical simulation. Dimensional analysis is adopted to acquire the principal governing parameters of the process, such as the ratio of yield strength to elastic modulus (Y/E) and the ratio of indentation depth to maximum penetration depth (h/hmax). The dimensional analysis provides a general understanding on the relationships among these dimensionless variables and their influence on the test results. Three kinds of indenters, that is, conical, Berkovich, and Vickers indenter, are comparatively studied. The dependence of pileup on the dimensionless parameters and its effect on the simulation results are also given. These findings are instructive for measuring the mechanical properties by instrumented indentation, especially when the effect of pileup is involved.


Introduction
Instrumented indentation, particularly nanoindentation, has been widely used for characterizing the mechanical properties of materials at small scales [1][2][3][4][5][6].One of the great advantages of this technique is that the mechanical properties, such as hardness and elastic modulus, can be directly determined from the indentation load-displacement curves, while the measurement of the actual hardness impression is not necessary.Recent developments of the nanoindentation technique make it capable of recording experimental data at load and displacement as small as micro-Newton and nanometer levels.
The commonly used definition of hardness is as follows: where  is hardness,  the load applied to the test surface, and  the projected contact area. was calculated by some researchers using the following equation [1,4]: where ℎ  is the contact depth and ℎ is the resulting penetration depth.If pileup effect is negligible, ℎ  could be given by where  is a constant that depends on indenter geometry and  is the elastic stiffness of the contact [4].One should notice that (3) cannot account for the plastic phenomenon of pileup, which may lead to significant change of the true contact area.
The elastic unloading curve of nanoindentation is approximately expressed in the power-law relation [4,7]: where  and  are empirically determined fitting parameters, depending on the geometry of the indenter.ℎ  is the final displacement of the indenter after complete unloading.Then (5) The relationship between the effective elastic modulus and the projected contact area and unloading stiffness is given by [1,4] where  is a dimensionless factor that accounts for deviations in stiffness caused by the lack of axial symmetry for Berkovich or Vickers indenters.The effective elastic modulus is defined by where  and ] are elastic modulus and Poisson's ratio of the specimen and   and ]  are those of the indenter.The elastic modulus of the specimen, , could then be calculated by (7).
In this work, the indenter is taken as rigid body for numerical simulation, and then  =  eff (1 − ] 2 ).
From ( 1) to (7), the hardness as well as the elastic modulus of the material could be determined in the following sequence: (a) Fit the unloading curve of the indentation; find values of  and , and determine  by (5).
(c) Determine the projected contact area  by (2).
The difficulty related to the above method is the determination of the contact depth ℎ  and the true projected contact area , especially when the phenomenon of pileup occurs.The resulting contact area might be greatly underestimated by using these equations if pileup effect exists.Bolshakov and Pharr [7] found by numerical simulation that pileup was significant for conical indenter if the ratio of final displacement to the maximum depth of penetration is larger than 0.7 (i.e., ℎ  /ℎ max > 0.7).Y.-T Cheng and C.-M Cheng [2,[8][9][10][11] did a series of analysis on conical indentation by finite element method (FEM).They analyzed their data by the method of dimensional analysis, which provided us a general understanding of this process.The dependencies of  and  on some critical dimensionless variables, such as /, ], and  ( is the yield strength of the tested material and  is the half angle of the indenter), were investigated.However, their numerical simulations were limited to 2D, which may represent the conical indentation fairly well but are not so representative for some 3D indenters such as Berkovich indenter and Vickers indenter.Dao et al. [12] conducted comprehensive analysis on the data derived from indentation by finite element simulation, where they also used 2D axial symmetric model instead of the real configuration.Giannakopoulos et al. [13] and Larsson et al. [14] investigated Vickers indenter and Berkovich indenter separately, but detailed comparison among these sharp indenters, such as Vickers, Berkovich, and conical, was not provided.Although the geometry functions of Berkovich indenter and Vickers indenter can be designed to be the same as that of the conical indenter, the true projected contact areas at the same indentation depth may still be different considering the effect of pileup.
Apart from geometric difference of indenters, there are several other factors that may influence the measurement of instrumented indentation.The thickness and properties of the test materials and the substrate effect have been studied systematically by many researchers [16][17][18][19][20][21].The stress state within the material also affects the testing results of indentation [22,23].In fact, Larsson et al. successfully related stress state to indentation testing results and found a way of determining residual stress within the material [24,25] by sharp indentation.
In this paper, we will focus on the differences of indenter geometries and their effect on the indentation results by means of numerical simulation.Dimensional analysis will be adopted to acquire the principal governing parameters.The effect of pileup on the simulation results will also be given.
Three kinds of indenters, namely, conical indenter, Berkovich indenter, and Vickers indenter, will be discussed separately.The half angle of the conical indenter is set to be 70.3 ∘ , while the half angles of Berkovich indenter and Vickers indenter are chosen to be 65.3 ∘ and 68 ∘ , respectively.The relationship between the projected contact area  and the contact depth ℎ  (i.e.,  1 (ℎ  )) can be obtained through simple calculation and all of the three indenters have equivalent geometry functions:

Dimensional Analysis
2.1.Dimensional Analysis for Loading.For a simple instrumented indentation test, the procedure is composed of two steps, that is, loading and unloading.During the loading process, the produced load  and projected contact area  are functions of all the independent governing parameters, which include the material parameters (Young's modulus , yield strength , and Poisson's ratio ]), the indenter geometry parameter (those in (2)), and the indentation depth ℎ.For a given indenter,  and  can be written as Among these governing parameters, we may choose two basic dimensions, that is, [] and [ℎ].The dimensions of the rest of the variables can then be described as Since the influence of ] on the measurement of hardness and elastic modulus is relatively small [1] (e.g., ] = 0.25 ± 0.1 produces only about 5% uncertainty in the calculated value of ), we will not incorporate the effect of ] for simplicity in this work.Conducting dimensional analysis by applying the Π theorem [26], the following equations can be obtained: where Then the above equations can be rewritten as Note that   and   are different functions from  Π and  Π , respectively.It is easy to see from ( 12) and ( 13) that  and  are all proportional to ℎ 2 .The hardness can then be calculated by Clearly, the hardness depends only on the material parameters  and , and it is independent from the indentation depth.Moreover, the dimensionless hardness, that is, /, is a function of only one parameter / from (14).

Dimensional Analysis for Unloading.
During unloading, the imposed load  depends on the following independent variables: , , ℎ, and maximum indentation depth ℎ max , while the final indentation depth ℎ  depends only on , , and ℎ max (the effect of ] is not considered here).We have Again, the dimensional analysis gives Equation (17) indicates that the final indentation depth is proportional to the maximum indentation depth and is dependent on /.Since the slope of the unloading curve is essential for calculating  and , we take derivative with respect to ℎ on both sides of ( 16) and obtain the values at ℎ = ℎ max : Or equivalently, The dimensionless form of the elastic contact stiffness /(ℎ max ) depends on /.

Numerical Analysis
Finite element analysis (FEA) of the indentation process has been carried out by using the commercial software Abaqus.The above-mentioned scaling relationships are carefully examined.These dimensionless functions, such as   ,   ,  ℎ , and   , are evaluated.The effect of pileup on the calculated hardness and elastic modulus will also be discussed.
As mentioned before, three kinds of indention, namely, conical, Vickers, and Berkovich, are simulated.All of the indenters have the same projected area function (8).The dimension of the tested specimen was set to be larger than ten times of the contact length in every direction to minimize the influence of the boundary.Similar to experimental test, the bottom surface of the specimen is set to be fixed and all of the remaining surfaces are set to be free.The indenter is free along the indentation direction only and the loads are given through applying displacements on the indenter.The contact between the indenter and specimen is set to be frictionless and the whole process is considered to be time-independent.
The indenters were taken as rigid body, while the tested materials were regarded as elastic-perfect-plastic.The yield strengths of the materials are set to range from 80 MPa to 8 GPa, which covers most of the engineering materials.Two elastic modulus values, that is, 100 GPa and 200 GPa, are considered.Poisson's ratio is set to be 0.38, corresponding to a material of niobium that will be used for comparison with the experimental observations.Actually, the values of Poisson's ratio taken from 0.2 to 0.4 will not lead to any significant difference [8].
A total number of 60,000 brick elements are meshed in the material and about 3,000 quadrilateral elements are meshed in the indenter.An example of the mesh of half of the conical indentation model is given in Figure 1.Intense mesh is given in the region of contact and sparse mesh in the far field to compromise between accuracy and efficiency.Mesh dependence of the simulation results was conducted first and Similar to experiments, the P-h curves are extracted directly from the calculated results of each simulation.Other properties such as hardness and elastic modulus are derived from the P-h curves following the procedure given in Section 1 of this paper.The errors related to FEA include modeling error, discretization error, and numerical error.Since all of the simulation models are identical in the sense of modeling, loading, meshing, contact condition, and solver and the only differences are the geometry of indenters and the properties of testing materials, the simulated results can accurately reflect the role of these different inputs in instrumented indentation.Actually, one advantage of the "numerical experiment" is that it excludes all other factors that may influence the results and focuses on one or two concerned terms, which is very hard to achieve in real experiments.the geometry functions of the three indenters are the same, their pileup patterns are rather different.These differences indicate that axial symmetric conical indenter may not be representative for the 3-side Berkovich indenter or 4-side Vickers indenter, especially when pileup plays a role.The evolution of force  with respect to indentation depth ℎ is plotted in Figure 3.The loading part of each indentation curve is fitted by the power law function.The fitting exponents, , for all the simulations are between 1.93 and 2.02, corresponding well with the prediction by (12).The true contact area,   , between the indenter and the material is easy to obtain from the simulation.The true projected contact area, , can be derived from   through simple geometric relationships.Figure 4 displays the A-h curves for the three kinds of indentations.Although they behave differently in Figure 2, the three kinds of indenters present identical Ah curves given the same  and .The loading part of the A-h curves is also fitted by power law function and it gives the values of  ranging from 2.0 to 2.1, in accordance with the prediction by (13).A reference area   , defined as the projected area calculated from indentation depth ℎ by ( 8) is given in Figure 4. Obviously,   is significantly smaller than the true projected areas when the indentation depth is larger than 1 m, signifying essential error may be introduced if the effect of pileup is ignored.

Results and Discussion
The dimensionless functions, that is,   = /(ℎ 2 ),   = /ℎ 2 ,  ℎ = ℎ  /ℎ max , and   = /(ℎ), were plotted with respect to / in Figures 5, 6, 7, and 8, respectively.It is obvious from Figure 5 that these derived quantities lie in a single curve for each indenter.Thus   is a function of /, as predicted by (12).Nevertheless, these functions are not exactly the same for each indenter.When / is larger than 0.01 (corresponding to materials like sodalime glass, high strength steels, titanium alloys, and so on), the differences among the data derived from the three indenters become evident.Similarly, the dimensionless function for projected contact area   is a function of / and is different for each indenter.But, different from   , the discrepancy of   becomes significant when / is smaller than 0.02 (corresponding to materials like aluminum alloys, low carbon steels, and so on).The dimensionless function  ℎ decreases monotonously with increasing / and appears independent of the geometry of the indenter (see Figure 7).ℎ  /ℎ max is a very useful quantity, because it is easy to obtain from the P-h curve and many researchers [7][8][9]14]  for pileup or sink-in.The dimensionless contact stiffness   becomes smaller as / increases but the differences among different indenters become more significant.Since the contact stiffness is essential for the derivation of the elastic modulus (see ( 6)), its scattering among different indenters may lead to significant variation of the calculated elastic modulus.
The projected contact area  is an important quantity that directly influences the measured hardness.The dependence of the normalized projected contact area /  on ℎ  /ℎ max is given in Figure 9.The true projected contact areas (solid symbols) and those calculated by the Oliver and Pharr (O/P) model [1] (hollow symbols) are all presented.The red dashed line, where /  is 1.0, denotes no pileup or sinkin.All of the hollow symbols are below the dashed line, showing that the O/P model considers sink-in only and could lead to significant error if ℎ  /ℎ max is larger than 0.8.The solid symbols intersect with the dashed line at around ℎ  /ℎ max = 0.85 for all the three indenters, indicating that direct penetration depth and area function could be used to calculate projected contact area  at this point (8).Our simulation results also reveal that less than 10% error will be introduced if pileup or sink-in effect is ignored when 0.8 < ℎ  /ℎ max < 0.9 (corresponding 0.008 < / < 0.02) for Berkovich indentation.However, this error increases to more than 20% for conical or Vickers indentation.This observation demonstrates that Berkovich indenter is less sensitive to pileup or sink-in than Vickers and conical indenters within this region.If pileup and sink-in are ignored, the maximum error could be larger than 50% for the conical indenter.
The derived elastic modulus of the material depends on the contact stiffness , the projected contact area , and the correction factor  (6).The differences of  and  could lead to variation of  for different indenters.Figure 10 displays the dependence of calculated  on ℎ  /ℎ max .The solid symbols represent  derived from the true projected contact area, while hollow symbols represent  derived from the O/P model.Obviously, all of the values of calculated  are larger than the input value and the O/P method could introduce more error if improper correction factor  is used.Figure 10 shows that the error of calculated  increases with increment of ℎ  /ℎ max , demonstrating that the correction factor  should not be a constant but a function of ℎ  /ℎ max ; that is,  =  + (ℎ  /ℎ max ).For the non-work-hardening materials in this work,  and  are 1.021 and 0.358, 0.832 and 0.526, and 0.759 and 0.512 for the Berkovich, Vickers, and conical indenters, respectively.Another noteworthy feature in Figure 10 is that Berkovich indentation will yield higher elastic modulus than Vickers and conical indentations, indicating again that replacing Berkovich and Vickers by conical indenter is not appropriate.An approximate relationship of  among different indenters can be given as  Berkovich = 1.06Vickers = 1.13Conical .This estimation leads to less than 5% error when ℎ  /ℎ max < 0.95 (corresponding to / > 0.002).
Figure 11 illustrates the variation of the normalized hardness  O/P / true with respect to ℎ  /ℎ max . O/P and  true are derived by using the O/P model and the true projected contact area, respectively.The O/P hardness, as expected, is larger than the true hardness for all cases.When ℎ  /ℎ max < 0.8,  O/P / true ≈ 1.1 and this approximation will introduce no more than 5% error for all the three indenters.However, this ratio increases sharply when ℎ  /ℎ max > 0.8.For conical indentation,  O/P / true rises monotonously up to 1.54, but the maximum values for Berkovich and Vickers indentation are less than 1.3.Thus, conical indentation fails again to represent Berkovich and Vickers indentations.
The hardness is normalized by yield stress of the material and plotted against ℎ  /ℎ max in Figure 12.In the case that  is calculated by the O/P method, similar trends of / are Advances in Materials Science and Engineering Figure 9: Influence of ℎ  /ℎ max on the projected contact area.The projected contact area is normalized by   , the area calculated by the indenter geometry function.The solid icons stand for the true projected contact area derived from finite element models, while the hollow icons represent the projected contact area calculated by Oliver and Pharr model [1].observed for all the three indenters in Figure 12, although there are still differences in the magnitudes.When ℎ  /ℎ max approaches unity, the maximum values of / are achieved and are larger than 4.2 for all indenters.When true hardness is applied, different patterns of / are found for different indenters.The relationship between hardness and yield stress,  = /, was firstly proposed by Tabor [27], where  is the constraint factor.Tabor analyzed the rigid/plastic deformation of Vickers indentation and claimed that the hardness of metals is about 2.9∼3.0 times its flow stress at 8∼10% strain; that is,  ≈ 2.9 ∼ 3.0.Since rigid/plastic deformation corresponds to ℎ  /ℎ max = 1, our numerical results could be compared with .Our simulation shows that / ≈ 3.3 for Vickers indentation when ℎ  /ℎ max approaches unity, slightly larger than Tabor's value.Different from Vickers and Berkovich, the value of / for conical indentation plateaus at about 2.7, which is close to that calculated by Bolshakov and Pharr [7].Clearly,  is not a constant but varies with ℎ  /ℎ max .The value of  also depends on the indenter geometry.For soft material, that is, ℎ  /ℎ max close to 1 (or / close to zero from Figure 7),  reaches 2.7, 3.3, and 3.6 for conical, Vickers, and Berkovich indentation, respectively.However,  might be lower than 2.0 for some hard materials.
It should be noted that ℎ  /ℎ max = 0.8 (or / = 0.02) is a critical point, below which the true projected contact area, hardness, and Tabor factor are almost proportional to the values derived from Oliver and Pharr model.In that case, the results by traditional method could still be used if a correction factor is applied (hardness divided by, e.g., 1.1).However, this factor could not be used if ℎ  /ℎ max > 0.8 (or / < 0.02), where the relation between true values and Oliver-Pharr results becomes different for different indenters.As / is less than 0.02 for most of the engineering materials, the corrections should be more careful.
Our experimental results [15] on niobium could be used to compare with some of the simulation results.The commercially pure niobium disks with diameter of 10.0 mm and thickness of 1.0 mm were processed by high pressure torsion (HPT) at room temperature.The Von Mises equivalent strain at the edge of the disk reaches about 90, while it is theoretically zero at the disk center.The large deformation by HPT substantially enhances the yield strength of the material while reducing its ductility.One distinct characteristic of HPT-ed materials is their vanished work-hardening ability, which is in keeping with the assumption of elastic-perfectplastic material in our simulation.In fact, the material undergoes deformation even at the center of the disk during HPT.The Berkovich hardness is 1.2 GPa at the niobium disk center [15], corresponding to about 400 MPa for yield stress  from Figure 12.Considering that the strength of as received niobium is less than 200 MPa, the enhancement in strength and at the same time reduction in ductility at the disk center would be significant.Therefore, assuming elasticperfect-plastic material property at the disk center will not lead to essential error.The HPT process does not change the elastic modulus remarkably.Thus, the value of / is larger at the edge than at the center of the disk.Figure 13 presents the indentation geometry at the center (a) and edge (b) of the HPT processed niobium disk.It can be observed that significant pileup appears at the center (a2), while neither sink-in nor pileup is observed at the edge (b2).The values of ℎ  /ℎ max are 0.93 and 0.81 for Figures 13(a) and 13(b), respectively.The corresponding values of / are estimated to be 0.0057 and 0.019 from Figure 7.This experimental observation of pileup or sink-in is well predicted by Figure 9.
Advances in Materials Science and Engineering

Summary and Conclusion
Three kinds of indentation, namely, conical, Berkovich, and Vickers, were comparatively studied by means of numerical simulations.The effect of pileup on the derived quantities was carefully examined.Dimensionless analysis was adopted to acquire the principal governing parameters of the process.
The main findings are summarized in the following: (1) The dimensionless functions   ,   ,  ℎ , and   , which represent loading, projected contact area, final indentation depth, and contact stiffness, respectively, are all functions of / only for the elastic-perfectplastic material considered in this work.However, these functions are not the same for different kinds of indenters.
(2) When ℎ  /ℎ max > 0.85, the effect of pileup becomes significant.Failure to account for the effect of pileup may lead to more than 50% error in calculating elastic modulus and hardness.
(3) The derived elastic modulus could be different if different indenters are adopted, even if the effect of pileup is considered.An approximate relationship of  of the sample measured using different indenters can be given by  Berkovich = 1.06Vickers = 1.13Conical .This estimation leads to less than 5% error when ℎ  /ℎ max < 0.95 (corresponding to / > 0.002).
The above conclusions are acquired based on a group of specific materials whose strain hardening is ignored and the effect of Poisson's ratio is not considered.Moreover, the timedependent (or strain-rate dependent) characteristics of the indentation process for some materials are not reckoned in for present work.Therefore, application of these results It can be observed that, at the center of the disk, significant pileup effect is present (a2), while, at the edge, neither sink-in nor pileup effect is observed (b2) [15].
should be conditional.The role of strain hardening and strain rate hardening might be significant in determining the properties of the testing materials and should be further investigated.

2
Advances in Materials Science and Engineering could be acquired by calculating the slope of the initial portion of the unloading curve: =  ℎ        ℎ=ℎ max =  max  (ℎ max − ℎ  ) −1 .

1 𝜇mFigure 1 :Figure 2 :
Figure 1: Finite element meshes of the numerical model of conical indentation (only half model is shown for better view).The image in the red box provides the details of the mesh in the region of contact.

Figure 2 Figure 3 :
Figure2displays the displacement field in the loading direction at a penetration depth of 3 m for the three indentations.The current  and  are 800 MPa and 100 GPa, respectively.Only half model is shown for better view.Note from the legend that positive values of displacement (U3) represent penetration into the material, while negative values indicate pileup.Clearly, pileup exists in each of the situations, demonstrating that corrections are required when using traditional process equations (2) and (3).Moreover, although

Figure 4 :Figure 5 :
Figure4: True projected area-displacement curves for three kinds of indenters.  is the reference area that is calculated from the indentation depth by using the area function of the indenters.The solid lines are the power law fitting of the loading parts. = 100 GPa,  = 800 MPa.

Figure 6 :
Figure 6: The variation of dimensionless function   = /ℎ 2 with respect to / for different kinds of indenters.

Figure 7 :
Figure 7: The variation of dimensionless function  ℎ = ℎ  /ℎ max with respect to / for different kinds of indenters.

Figure 8 :
Figure 8: The variation of dimensionless function   = /(ℎ) with respect to / for different kinds of indenters.

Figure 10 :
Figure 10: Influence of ℎ  /ℎ max on the calculated elastic modulus.The modulus is normalized by the input  of the numerical model.The solid icons stand for the elastic modulus derived from true projected contact areas (6), while the hollow icons represent the modulus derived by Oliver and Pharr model[1].The solid and dashed lines are linear fitting of solid icons and hollow icons, respectively.

Figure 11 :
Figure 11: Effect of ℎ  /ℎ max on  O/P / true , the ratio of hardness derived by Oliver and Pharr model to that calculated by using true projected contact area.

Figure 12 :
Figure12: Influence of ℎ  /ℎ max on /.The solid icons stand for the hardness derived from true projected contact areas, while the hollow icons represent those derived by Oliver and Pharr model[1].

Figure 13 :
Figure13: Appearance of the indentation geometry by AFM at the center (a) and edge (b) of the HPTed niobium disk.It can be observed that, at the center of the disk, significant pileup effect is present (a2), while, at the edge, neither sink-in nor pileup effect is observed (b2)[15].