Numerical Simulation of the Depth-Sensing Indentation Test with Knoop Indenter

Depth-sensing indentation (DSI) technique allows easy and reliable determination of two mechanical properties of materials: hardness and Young’s modulus. Most of the studies are focusing on the Vickers, Berkovich, and conical indenter geometries. In case of Knoop indenter, the existing experimental and numerical studies are scarce. The goal of the current study is to contribute for the understanding of the mechanical phenomena that occur in the material under Knoop indention, enhancing and facilitating the analysis of its results obtained in DSI tests. For this purpose, a finite element code, DD3IMP, was used to numerically simulate the Knoop indentation test. A finite element mesh was developed and optimized in order to attain accurate values of the mechanical properties. Also, a careful modeling of the Knoop indenter was performed to take into account the geometry and size of the imperfection (offset) of the indenter tip, as in real cases.


Introduction
Depth-sensing indentation (DSI) tests are typically used to evaluate the hardness and Young's modulus of materials.They can also be used to extract the uniaxial mechanical properties of bulk and composite materials, such as the yield stress and the strain hardening parameter (see, e.g., [1][2][3][4]).The most common hardness testing methods were developed in the early twentieth century.They are typically performed using spherical, conical and pyramidal indenters with Vickers and Berkovich geometries.In addition, the Knoop pyramid geometry is also used [5].
Hardness tests with the Knoop indenter have been valuable in the mechanical characterization of some materials, such as thin coatings [6,7] and biological materials (e.g., dental tissue [8]).Another important application of this hardness tests is in the field of the gradient materials obtained by severe plastic deformation (see, e.g., [9][10][11]), for which it is required to determine the mechanical properties in thin samples and/or thin surface layers of the samples.
In fact, the Knoop indenter geometry leads to wider and shallower indentations for a given applied load than the Vickers and Berkovich geometries.This makes the Knoop hardness test particularly attractive in some cases as for the determination of the near-surface properties and the characterization of brittle materials.Moreover, the results obtained in the Knoop hardness test are sensitive to the indenter orientation, making it a useful tool to analyze the materials anisotropy (e.g., [12]).
Although Knoop indenter is relatively common to use, it has not received suitable attention with respect to the study of some of the peculiarities of the test itself.One of the first studies using numerical simulation of the Knoop indentation test was performed by Rabinovich and Sarin [13], in their study, the Knoop indentation was analyzed in the context of linear elasticity.Giannakopoulos [14] presented analytical results on the response of frictionless and adhesionless contact of flat, linear elastic and visco-elastic isotropic surfaces penetrated by pyramidal indenters including the Knoop geometry.Giannakopoulos and Zisis [15] studied the Knoop indentation of elastic and elastoplastic materials with and without strain hardening.More recently, Giannakopoulos and Zisis [16] presented a finite element study on the adhesionless contact of flat surfaces by Knoop indenter.Only a few experimental studies, such as those by Riester et al. [17,18] were performed in order to clarify the analysis procedure of DSI tests with the Knoop geometry.Recently, Ghorbal et al. [19] performed an experimental work on ceramic materials in order to compare the conventional Knoop and Vickers hardness.
In this context, the present study intends to contribute to clarify some aspects of the Knoop indentation, particularly those related to the analysis of depth sensing indentation (DSI) results.For this purpose, three-dimensional numerical simulations of the indentation tests were performed, using various pyramidal indenter geometries, from Vickers to Knoop.A systematic study is accomplished using materials with different mechanical properties whose ratio between the residual indentation depth after unloading (h f ) and the indentation depth at maximum load (h max ) is in the range 0.022 < h f /h max < 0.984.The correction factor, β, required to determine the Young's modulus, is evaluated as a function of the indenter geometry, from Vickers to Knoop.In addition, numerical simulations using flat indenters with equivalent lozenge geometries were performed, in order to better understand the role of the pyramidal geometry.

Theoretical Aspects
The ability of the ultramicrohardness equipment to register the load versus the depth indentation, during the test, enables us to evaluate not only the hardness, but also other properties, such as the Young's modulus.Based on the Sneddon relationship [20] between the indentation parameters and Young's modulus, Doerner and Nix [21] have proposed an equation that relates the Young's modulus with the compliance, C, of the unloading curve at the point of the maximum load and the contact area, A c , such as: where β is the geometrical correction factor for the indenter geometry.The specimen's Young's modulus, E s , is obtained using the equation: where E and ϑ are the Young's modulus and the Poisson's ratio, respectively, of the specimen (s) and of the indenter (i).When performing numerical simulations, the indenter can in a first approximation considered as rigid for simplicity (e.g., [22]), and so 1 − ϑ 2 i /E i = 0.The accuracy of the Young's modulus results, obtained with the above equations, depends on the correct evaluation of the contact area and compliance.The contact area, A c , can be evaluated by two different procedures.One procedure uses the contour of the indentation in the finite element mesh, in order to make the results independent of the formation of pile-up or sink-in.The other is the usual experimental procedure, which makes use of the compliance, C, evaluated by fitting the unloading part of the load-indentation depth curve, (P − h), using the power law [22]: where T and m are constants obtained by fit and h 0 is the lower value of the indentation depth, h, used in the fitted region, corresponding to a load value P 0 , during unloading.The upper part of the unloading curve taken into account in the fits is 70% [22].Once the value of compliance, C, is known, the contact indentation depth, h c , that allows the calculation of the contact area according to the , where h max is the indentation depth at the maximum load, P max , and ε is a parameter equal to 0.72 for conical or pyramidal (Berkovich [23], Vickers [22] and Knoop [18] indenters).
An approach that can be used for evaluating the correction factor, β, was proposed by Joslin and Oliver [24], combining the hardness and Young's modulus equations: where H = P/A c is the hardness, P the maximum applied load and S the stiffness (S = 1/C).The ratio between the maximum applied load and the square of the stiffness, P/S 2 , is an experimentally measurable parameter that is independent of the contact area and so of the penetration depth [24].

Numerical Simulation and Materials
Three-dimensional numerical simulations of the hardness tests were carried out, using a finite element (FE) in-house code DD3IMP.This FE code, which has already been tested in the case of Vickers, Berkovich and conical indentation of bulk materials and thin films (see, e.g., [22,25,26]), allows simulating the hardness tests with any type of indenter shape, taking into account contact with friction between the indenter and the sample [22,27,28].The mechanical model, which is the basis of the DD3IMP code, considers the hardness test as a quasi-static process that occurs in the large plastic deformations domain.In DD3IMP, the contact with friction problem is modeled using a classical Coulomb's law.To relate the static equilibrium problem to the contact with friction, an augmented Lagrangean method is applied to the mechanical formulation.This leads to a system of non-linear equations, where the kinematic (material displacements) and static variables (contact forces) are the final unknowns [27,28].To solve this problem, the code uses a fully implicit Newton-Raphson-type algorithm.All non-linearities, induced by the elastoplastic behavior of the material and by the contact with friction, are treated in a single iterative loop [27,28].
In the current study, the friction between the indenter and the deformable body was assumed to have a friction coefficient of 0.16.This is a commonly used value and leads to a better description of the indentation process than if frictionless contact is assumed [22,26].

Indenters
The Knoop indenter has a pyramidal geometry, with a lozenge-shaped base having one diagonal (L) 7.11 times longer than the other (m).The angles between the edges (apical angles) are 172.5 o for the long edges and 130 o for the short edges, as shown in Figure 1.unloading curve taken into account in the fits is 70% [22].Once the value of compliance, , is known, the contact indentation depth, ℎ , that allows the calculation of the contact area according to the geometry of the indenter ( = (ℎ )), is determined by the equation ℎ = ℎ −  , where ℎ is the indentation depth at the maximum load,  , and  is a parameter equal to 0.72 for conical or pyramidal (Berkovich [23], Vickers [22] and Knoop [18] indenters).
An approach that can be used for evaluating the correction factor, , was proposed by Joslin and Oliver [24], combining the hardness and Young's modulus equations: where  = / is the hardness, P the maximum applied load and  the stiffness ( = 1/).The ratio between the maximum applied load and the square of the stiffness, /  , is an experimentally measurable parameter that is independent of the contact area and so of the penetration depth [24].

Numerical Simulation and Materials
Three-dimensional numerical simulations of the hardness tests were carried out, using a finite element (FE) in-house code DD3IMP.This FE code, which has already been tested in the case of Vickers, Berkovich and conical indentation of bulk materials and thin films (see, e.g., [22,25,26]), allows simulating the hardness tests with any type of indenter shape, taking into account contact with friction between the indenter and the sample [22,27,28].The mechanical model, which is the basis of the DD3IMP code, considers the hardness test as a quasi-static process that occurs in the large plastic deformations domain.In DD3IMP, the contact with friction problem is modeled using a classical Coulomb's law.To relate the static equilibrium problem to the contact with friction, an augmented Lagrangean method is applied to the mechanical formulation.This leads to a system of non-linear equations, where the kinematic (material displacements) and static variables (contact forces) are the final unknowns [27,28].To solve this problem, the code uses a fully implicit Newton-Raphson-type algorithm.All non-linearities, induced by the elastoplastic behavior of the material and by the contact with friction, are treated in a single iterative loop [27,28].
In the current study, the friction between the indenter and the deformable body was assumed to have a friction coefficient of 0.16.This is a commonly used value and leads to a better description of the indentation process than if frictionless contact is assumed [22,26].

Indenters
The Knoop indenter has a pyramidal geometry, with a lozenge-shaped base having one diagonal (L) 7.11 times longer than the other (m).The angles between the edges (apical angles) are 172.5 for the long edges and 130 for the short edges, as shown in Figure 1.The ideal (in the absence of pile-up or sink-in) indentation contact area,  , of the Knoop indenter as function of the indentation contact depth is given by:  = 2ℎ tan tan = 65.4ℎ, The ideal (in the absence of pile-up or sink-in) indentation contact area, A c , of the Knoop indenter as function of the indentation contact depth is given by: The Knoop indenter geometry was modeled using parametric Bezier surfaces, which allows a satisfactory description of the indenter tip, namely an imperfection such as occurs in the real geometry, similar to the case of offset in the Vickers indenter [25,29].The modeled Knoop indenter, shown in Figure 2, has a tip imperfection consisting of a plane normal to the indenter axis with an area equal to 0.0032 µm 2 (this value is the same as the experimentally found, for the Vickers indenter, by Antunes et al. [29]).where ℎ is the ideal indentation contact depth,  = 65 and  = 86.25 are the semi-apical angles of the Knoop indenter.The Knoop indenter geometry was modeled using parametric Bezier surfaces, which allows a satisfactory description of the indenter tip, namely an imperfection such as occurs in the real geometry, similar to the case of offset in the Vickers indenter [25,29].The modeled Knoop indenter, shown in Figure 2, has a tip imperfection consisting of a plane normal to the indenter axis with an area equal to 0.0032 μm (this value is the same as the experimentally found, for the Vickers indenter, by Antunes et al. [29]).Due to the tip imperfection, the indenter area function does not match the ideal area function above mentioned (Equation ( 5)).Table 1 shows the area function of the Knoop indenter (ratio, R = L/m = 7.11) used in the numerical simulations.This table also shows four others area functions of pyramidal indenters, used in this study, with different values of the ratio, R, between the diagonals of the indenter (R = 1, 2.5, 4 and 5.5), where the ratio R = 1 corresponds to the Vickers indenter.For pyramidal indenters, other than Vickers and Knoop, the angles  and  were chosen such that the tangents of  and  follows a quasi-linear evolution with R, as shown in Figure 3.In this way it is possible to study the extent to which the deviation of Vickers geometry towards the Knoop geometry influences the indentation results obtained with DSI experiments.Due to the tip imperfection, the indenter area function does not match the ideal area function above mentioned (Equation ( 5)).Table 1 shows the area function of the Knoop indenter (ratio, R = L/m = 7.11) used in the numerical simulations.This table also shows four others area functions of pyramidal indenters, used in this study, with different values of the ratio, R, between the diagonals of the indenter (R = 1, 2.5, 4 and 5.5), where the ratio R = 1 corresponds to the Vickers indenter.For pyramidal indenters, other than Vickers and Knoop, the angles θ 1 and θ 2 were chosen such that the tangents of θ 1 and θ 2 follows a quasi-linear evolution with R, as shown in Figure 3.In this way it is possible to study the extent to which the deviation of Vickers geometry towards the Knoop geometry influences the indentation results obtained with DSI experiments.where ℎ is the ideal indentation contact depth,  = 65 and  = 86.25 are the semi-apical angles of the Knoop indenter.The Knoop indenter geometry was modeled using parametric Bezier surfaces, which allows a satisfactory description of the indenter tip, namely an imperfection such as occurs in the real geometry, similar to the case of offset in the Vickers indenter [25,29].The modeled Knoop indenter, shown in Figure 2, has a tip imperfection consisting of a plane normal to the indenter axis with an area equal to 0.0032 μm (this value is the same as the experimentally found, for the Vickers indenter, by Antunes et al. [29]).Due to the tip imperfection, the indenter area function does not match the ideal area function above mentioned (Equation ( 5)).Table 1 shows the area function of the Knoop indenter (ratio, R = L/m = 7.11) used in the numerical simulations.This table also shows four others area functions of pyramidal indenters, used in this study, with different values of the ratio, R, between the diagonals of the indenter (R = 1, 2.5, 4 and 5.5), where the ratio R = 1 corresponds to the Vickers indenter.For pyramidal indenters, other than Vickers and Knoop, the angles  and  were chosen such that the tangents of  and  follows a quasi-linear evolution with R, as shown in Figure 3.In this way it is possible to study the extent to which the deviation of Vickers geometry towards the Knoop geometry influences the indentation results obtained with DSI experiments.In order to better understand some aspects related with the influence of the indenter geometry on the mechanical properties evaluation, numerical simulations using lozenge-shaped flat indenters were also performed.The flat indenters were modulated with five values of the ratio, R = L/m, between the diagonals of the lozenge, as for the pyramidal indenters (R = 1, 2.5, 4, 5.5, and 7.11).

Finite Element Mesh
The test sample used in the numerical simulations has both radius and thickness equal to 40 µm.Its discretization was performed using three-linear eight-node isoparametric hexahedrons.Due to geometrical and material symmetries in the X = 0 and Z = 0 planes, only a quarter of the sample was used in the numerical simulation, as shown in Figure 4.
Metals 2018, 8, x FOR PEER REVIEW 5 of 15 In order to better understand some aspects related with the influence of the indenter geometry on the mechanical properties evaluation, numerical simulations using lozenge-shaped flat indenters were also performed.The flat indenters were modulated with five values of the ratio, R = L/m, between the diagonals of the lozenge, as for the pyramidal indenters (R = 1, 2.5, 4, 5.5, and 7.11).

Finite Element Mesh
The test sample used in the numerical simulations has both radius and thickness equal to 40 µm.Its discretization was performed using three-linear eight-node isoparametric hexahedrons.Due to geometrical and material symmetries in the X = 0 and Z = 0 planes, only a quarter of the sample was used in the numerical simulation, as shown in Figure 4.The FE mesh is composed by 17,850 elements.The size of the elements in the indentation region is 0.055 μm.This refinement has proven to provide accurate values of the indentation contact area, when measured using the contour of the indentation, in case of Vickers and Berkovich geometries, with equal value of offset area (see, e.g., [22,26]).In fact, the Young's modulus values obtained from the indentation contact area, evaluated using the contour of nodes in the FE mesh in contact with the indenter at maximum load presents an error less than 1%, when compared with the values used as input in the numerical simulation (e.g., [26]).In the present study, the size of the central region of the finite element mesh, especially refined, is larger than in these cases, and thus the total number of elements is approximately three times greater, in order to take into account the elongated geometry of the Knoop indenter.

Materials
Three-dimensional numerical simulations of depth-sensing indentation with pyramidal indenters were carried out using 45 fictitious materials, whose mechanical properties are shown in Table 2.In order to cover a wide range of materials used in engineering applications, five values of yield stress (0.2, 2, 6, 10, and 20 GPa), three values of Young's modulus (70, 200 and 400 GPa) and of strain hardening parameter of the Swift law (0.01, 0.15 and 0.3), were taken into account.
The plastic behavior of the materials is described by the von Mises yield criterion and the Swift hardening law:  = ( +  ) where  and  are the equivalent stress and plastic strain, respectively, and ,  and n (strain hardening coefficient) are the material parameters (the yield stress is:  =  ); the parameter  was considered to be equal to 0.005.The elastic behaviour is isotropic and described by the generalised Hooke's law; the Poisson's ratio,  is 0.3, for all simulations.The FE mesh is composed by 17,850 elements.The size of the elements in the indentation region is 0.055 µm.This refinement has proven to provide accurate values of the indentation contact area, when measured using the contour of the indentation, in case of Vickers and Berkovich geometries, with equal value of offset area (see, e.g., [22,26]).In fact, the Young's modulus values obtained from the indentation contact area, evaluated using the contour of nodes in the FE mesh in contact with the indenter at maximum load presents an error less than 1%, when compared with the values used as input in the numerical simulation (e.g., [26]).In the present study, the size of the central region of the finite element mesh, especially refined, is larger than in these cases, and thus the total number of elements is approximately three times greater, in order to take into account the elongated geometry of the Knoop indenter.

Materials
Three-dimensional numerical simulations of depth-sensing indentation with pyramidal indenters were carried out using 45 fictitious materials, whose mechanical properties are shown in Table 2.In order to cover a wide range of materials used in engineering applications, five values of yield stress (0.2, 2, 6, 10, and 20 GPa), three values of Young's modulus (70, 200 and 400 GPa) and of strain hardening parameter of the Swift law (0.01, 0.15 and 0.3), were taken into account.
The plastic behavior of the materials is described by the von Mises yield criterion and the Swift hardening law: σ = k(ε + ε 0 ) n where σ and ε are the equivalent stress and plastic strain, respectively, and k, ε 0 and n (strain hardening coefficient) are the material parameters (the yield stress is: σ y = kε n 0 ); the parameter ε 0 was considered to be equal to 0.005.The elastic behaviour is isotropic and described by the generalised Hooke's law; the Poisson's ratio, ϑ is 0.3, for all simulations.

Indentation Geometry and Equivalent Plastic Strain Distributions
A study of the indentation geometry with the Knoop indenter (R = 7.11) was performed using only three of the materials in Table 2, covering different values in the possible range of h f /h max .Table 3 shows the mechanical properties of these materials, which have two Young's modulus values, 200 and 400 GPa, yield stress of 0.2, 6, and 20 GPa, and two values of the strain hardening parameter, n ≈ 0 and n = 0.3.Table 3 also includes the value of the ratio between the indentation depth after unloading and the indentation depth at maximum load, h f /h max .This ratio, which can be easily obtained from the experimental load-unloading curve, is independent of the maximum indentation depth, for a given material in cases of conical [30] and Vickers [22] indentations.Its range of values is from 0 to 1, which corresponds to materials with purely elastic and rigid-plastic behaviors, respectively.All numerical simulations of the hardness test were performed up to the same maximum indentation depth, h max = 0.2 µm.Figure 5 shows the indentation profiles, obtained for the three materials at maximum load (Figure 5a,c) and after unloading (Figure 5b,d), along the two diagonals, the shorter, m, and the longer, L, respectively.Figure 5a shows that for the shorter diagonal, the indentation profile obtained at the maximum load depends on the material mechanical properties.In case of materials M2 and M3 the indentation surface sink-in.This behavior can be associated with the value of the ratio h f /h max that is equal to 0.40 and 0.25 (for the materials M2 and M3, respectively), and with the high value of the strain hardening parameter, n.In case of M1 material, with a ratio h f /h max equal to 0.97, i.e., close to 1, the indentation surface does not present sink-in or pile-up.On the other hand, the indentation profile obtained at maximum load for the longer diagonal, almost does not depends on the materials mechanical properties (Figure 5c).After unloading, the indentation profiles corresponding to the materials M2 and M3 shown elastic recovery in both diagonals (Figure 5b,d).Moreover, the amount of elastic recover increases with the decrease of the ratio ℎ /ℎ .In case of M1 material, with a ratio ℎ /ℎ = 0.97 and without strain hardening, the indentation profile along the short diagonal shows pile-up (Figure 5b).In fact, for other indenter geometries (conical and Vickers), the pile-up appears for values of the ratio ℎ /ℎ higher than 0.8 in materials without strain hardening (see, e.g., [22,30]).It should be noted that, for a given material, the ℎ /ℎ ratio correlates with the value of the H/ ratio between the hardness and the Young's modulus, and slightly depends on the strain hardening parameter.
Figure 6 shows the comparison between the Knoop indentation diagonals at maximum load.To make possible this comparison, the indentation profiles were normalized by considering the value of the R-ratio between the indenter diagonals (R equal to 7.11).Figure 6a shows that in the material M1 the two diagonal exhibit a similar behavior.In the case of the materials M2 and M3 sink-in is observed along the short diagonal.Moreover, the sink-in slightly increases with the value of the ratio ℎ /ℎ .After unloading, the indentation profiles corresponding to the materials M2 and M3 shown elastic recovery in both diagonals (Figure 5b,d).Moreover, the amount of elastic recover increases with the decrease of the ratio h f /h max .In case of M1 material, with a ratio h f /h max = 0.97 and without strain hardening, the indentation profile along the short diagonal shows pile-up (Figure 5b).In fact, for other indenter geometries (conical and Vickers), the pile-up appears for values of the ratio h f /h max higher than 0.8 in materials without strain hardening (see, e.g., [22,30]).It should be noted that, for a given material, the h f /h max ratio correlates with the value of the H/E ratio between the hardness and the Young's modulus, and slightly depends on the strain hardening parameter.
Figure 6 shows the comparison between the Knoop indentation diagonals at maximum load.To make possible this comparison, the indentation profiles were normalized by considering the value of the R-ratio between the indenter diagonals (R equal to 7.11).Figure 6a shows that in the material M1 the two diagonal exhibit a similar behavior.In the case of the materials M2 and M3 sink-in is observed along the short diagonal.Moreover, the sink-in slightly increases with the value of the ratio h f /h max .After unloading, the indentation profiles corresponding to the materials M2 and M3 shown elastic recovery in both diagonals (Figure 5b,d).Moreover, the amount of elastic recover increases with the decrease of the ratio ℎ /ℎ .In case of M1 material, with a ratio ℎ /ℎ = 0.97 and without strain hardening, the indentation profile along the short diagonal shows pile-up (Figure 5b).In fact, for other indenter geometries (conical and Vickers), the pile-up appears for values of the ratio ℎ /ℎ higher than 0.8 in materials without strain hardening (see, e.g., [22,30]).It should be noted that, for a given material, the ℎ /ℎ ratio correlates with the value of the H/ ratio between the hardness and the Young's modulus, and slightly depends on the strain hardening parameter.
Figure 6 shows the comparison between the Knoop indentation diagonals at maximum load.To make possible this comparison, the indentation profiles were normalized by considering the value of the R-ratio between the indenter diagonals (R equal to 7.11).Figure 6a shows that in the material M1 the two diagonal exhibit a similar behavior.In the case of the materials M2 and M3 sink-in is observed along the short diagonal.Moreover, the sink-in slightly increases with the value of the ratio ℎ /ℎ .The effect of the indenter geometry on the equivalent plastic strain distribution of the indentations was also studied.Figure 7a,c   For each material, the maximum values of the equivalent plastic strain are quite similar for the Knoop and Vickers indenters.However, in the case of Knoop, the equivalent plastic strain distributions are asymmetric with respect to the indenter axes.The maximum values of the equivalent plastic strain are observed along the longest diagonal.For M1 material (ℎ /ℎ = 0.97), the maximum value of the equivalent plastic strain is slightly higher in the Vickers indentation The effect of the indenter geometry on the equivalent plastic strain distribution of the indentations was also studied.Figure 7a,c shows the equivalent plastic strain distributions obtained at the maximum load in the numerical simulations of the materials M1 (h f /h max = 0.97) and M3 (h f /h max = 0.25) using the Knoop indenter.For comparison, the same figure also shows the same distributions obtained with the Vickers indenter (Figure 7b,d).The effect of the indenter geometry on the equivalent plastic strain distribution of the indentations was also studied.Figure 7a,c   For each material, the maximum values of the equivalent plastic strain are quite similar for the Knoop and Vickers indenters.However, in the case of Knoop, the equivalent plastic strain distributions are asymmetric with respect to the indenter axes.The maximum values of the equivalent plastic strain are observed along the longest diagonal.For M1 material (ℎ /ℎ = 0.97), the maximum value of the equivalent plastic strain is slightly higher in the Vickers indentation For each material, the maximum values of the equivalent plastic strain are quite similar for the Knoop and Vickers indenters.However, in the case of Knoop, the equivalent plastic strain distributions are asymmetric with respect to the indenter axes.The maximum values of the equivalent plastic strain are observed along the longest diagonal.For M1 material (h f /h max = 0.97), the maximum value of the equivalent plastic strain is slightly higher in the Vickers indentation (≈0.686) than for Knoop (≈0.674).The maximum plastic strain region is located just at the surface in the edge regions of the indentation (Figure 7a,b).In case of the M3 material (h f /h max = 0.25), as shown in Figure 7c,d, the maximum Metals 2018, 8, 885 9 of 15 value of equivalent plastic strain is also somewhat higher for the Vickers indentation (≈0.159), than for the Knoop (≈0.149).However, for this material, the region with maximum equivalent plastic strain is located beneath the indentation surface, whatever the indentation geometry.

Indentation Contact Area and Young's Modulus
The results of the numerical simulation of the hardness tests with the Knoop indenter, for all materials in Table 2, were used to study the influence of the mechanical properties on the evaluation of the indentation contact area and, consequently, on the Young's modulus.As mentioned above, the indentation contact area was determinate using two different procedures: one of them uses the load-unloading curve as in the experimental DSI procedure for evaluating this area, A h c , and the other considers the contour of nodes of the FE mesh in contact with the indenter at maximum load, A FE (numerical contact area).
Figure 8 shows the evolution of the indentation contact areas, A h c and A FE , as a function of the ratio h f /h max .These contact areas are normalized with respect to the reference area, A REF , corresponding to the obtained with the area function of the indenter (Equation ( 5), with h c equal to the indentation depth, which does not take into account the pile-up or sink-in formation).Figure 8a shows that the contact area, A h c , is independent of the strain hardening parameter and Young's modulus, whatever the value of the ratio h f /h max .Moreover, the ratio A h c /A REF is always less than 1. Figure 8b shows that the numerically calculated contact area, A FE , is nearly independent from the strain hardening parameter and Young's modulus, for the ratio h f /h max < 0.6.However, for h f /h max > 0.6 the normalized contact area depends on the strain hardening parameter.The ratio A FE /A REF is even higher than 1 after a value of h f /h max that depends on the strain hardening parameter (h f /h max equal to 0.85, 0.90 and 0.95 for n equal to 0, 0.15 and 0.3, respectively), indicating pile-up formation.For both cases, A h c /A REF and A FE /A REF , similar behaviors were previously observed for the case of the conical, Vickers and Berkovich indenters (see, e.g., [22,26,30]).
(≈0.686) than for Knoop (≈0.674).The maximum plastic strain region is located just at the surface in the edge regions of the indentation (Figure 7a,b).In case of the M3 material (ℎ /ℎ = 0.25), as shown in Figure 7c,d, the maximum value of equivalent plastic strain is also somewhat higher for the Vickers indentation (≈0.159), than for the Knoop (≈0.149).However, for this material, the region with maximum equivalent plastic strain is located beneath the indentation surface, whatever the indentation geometry.

Indentation Contact Area and Young's Modulus
The results of the numerical simulation of the hardness tests with the Knoop indenter, for all materials in Table 2, were used to study the influence of the mechanical properties on the evaluation of the indentation contact area and, consequently, on the Young's modulus.As mentioned above, the indentation contact area was determinate using two different procedures: one of them uses the loadunloading curve as in the experimental DSI procedure for evaluating this area,  , and the other considers the contour of nodes of the FE mesh in contact with the indenter at maximum load,  (numerical contact area).
Figure 8 shows the evolution of the indentation contact areas,  and  , as a function of the ratio ℎ /ℎ .These contact areas are normalized with respect to the reference area,  , corresponding to the obtained with the area function of the indenter (Equation ( 5), with ℎ equal to the indentation depth, which does not take into account the pile-up or sink-in formation).Figure 8a shows that the contact area,  , is independent of the strain hardening parameter and Young's modulus, whatever the value of the ratio ℎ /ℎ .Moreover, the ratio  / is always less than 1. Figure 8b shows that the numerically calculated contact area,  , is nearly independent from the strain hardening parameter and Young's modulus, for the ratio ℎ /ℎ < 0.6 .However, for ℎ /ℎ 0.6 the normalized contact area depends on the strain hardening parameter.The ratio  / is even higher than 1 after a value of ℎ /ℎ that depends on the strain hardening parameter (ℎ /ℎ equal to 0.85, 0.90 and 0.95 for n equal to 0, 0.15 and 0.3, respectively), indicating pile-up formation.For both cases,  / and  / , similar behaviors were previously observed for the case of the conical, Vickers and Berkovich indenters (see, e.g., [22,26,30]).Figure 9 shows the evolution of the Young's modulus  and  normalized by the input value used in the numerical simulation,  , as a function of the ratio ℎ /ℎ .The average ratio of  / is close to 1.419, except for values of ℎ /ℎ approaching 1, for which the ratio  / increases for low values of the strain hardening parameter of the materials (n = 0 and 0.15).This is a consequence of the pile-up formation during the indentation of these materials, as was also observed Figure 9 shows the evolution of the Young's modulus E h c and E FE normalized by the input value used in the numerical simulation, E REF , as a function of the ratio h f /h max .The average ratio of E h c /E REF is close to 1.419, except for values of h f /h max approaching 1, for which the ratio A h c /A REF increases for low values of the strain hardening parameter of the materials (n = 0 and 0.15).This is a consequence of the pile-up formation during the indentation of these materials, as was also observed for Vickers indenters [25].The ratios of E FE /E REF are slightly higher than those of E h c /E REF , being in average close to 1.469.Slight differences between both ratios, E h c /E REF and E FE /E REF , were already Metals 2018, 8, 885 10 of 15 observed for the Vickers indenter [25].Under these conditions, in the experimental indentation tests, a value for the correction factor β of 1.419 is required to be used for the contact area evaluation (h c ) from the unloading curve, in case of the Knoop indentation.This value is essentially different from other indenter geometries such as conical, Vickers and Berkovich, which require β values around 1.034, 1.055, 1.081, respectively (see, e.g., [22,26]).
Metals 2018, 8, x FOR PEER REVIEW 10 of 15 for Vickers indenters [25].The ratios of  / are slightly higher than those of  / , being in average close to 1.469.Slight differences between both ratios,  / and  / , were already observed for the Vickers indenter [25].Under these conditions, in the experimental indentation tests, a value for the correction factor  of 1.419 is required to be used for the contact area evaluation (ℎ ) from the unloading curve, in case of the Knoop indentation.This value is essentially different from other indenter geometries such as conical, Vickers and Berkovich, which require  values around 1.034, 1.055, 1.081, respectively (see, e.g., [22,26]).Equation ( 4) was used in order to confirm the above correction factor . Figure 10 shows the ratio / versus / obtained by numerical simulation of all materials in Table 2, for the case of the Knoop and Vickers indenters.The reduced Young's modulus,  , was determined considering the input Young's modulus and Poisson ratio, , determined using the contact area evaluated from the contour of the indentations.Equation ( 4) was used in order to confirm the above correction factor β. Figure 10 shows the ratio P/S 2 versus H/E 2 r obtained by numerical simulation of all materials in Table 2, for the case of the Knoop and Vickers indenters.The reduced Young's modulus, E r , was determined considering the input Young's modulus and Poisson ratio, H, determined using the contact area evaluated from the contour of the indentations.for Vickers indenters [25].The ratios of  / are slightly higher than those of  / , being in average close to 1.469.Slight differences between both ratios,  / and  / , were already observed for the Vickers indenter [25].Under these conditions, in the experimental indentation tests, a value for the correction factor  of 1.419 is required to be used for the contact area evaluation (ℎ ) from the unloading curve, in case of the Knoop indentation.This value is essentially different from other indenter geometries such as conical, Vickers and Berkovich, which require  values around 1.034, 1.055, 1.081, respectively (see, e.g., [22,26]).Equation ( 4) was used in order to confirm the above correction factor . Figure 10 shows the ratio / versus / obtained by numerical simulation of all materials in Table 2, for the case of the Knoop and Vickers indenters.The reduced Young's modulus,  , was determined considering the input Young's modulus and Poisson ratio, , determined using the contact area evaluated from the contour of the indentations.The straight-lines in Figure 10 pass through the origin of the axes as indicated by Equation (4) (independently of the indenter geometry, curves match for H/E 2 r = 0, i.e., for materials with rigid-plastic behaviour, which corresponds to the ratio h f /h max = 1).The β factor is evaluated from the slope, ρ, of the straight line, related with β through ρ = (π/4β 2 ).Using this procedure, for the values of the R-ratio of 1, 2.5, 4, 5.5, and 7.11, the β correction factor obtained were 1.054, 1.141, 1.232, 1.351, and 1.473, respectively.The value of 1.473 for the Knoop indenter is very close to that mentioned above (1.469),obtained from the ratios E FE /E REF (Figure 9b).

Flat Indenter
In order to understand if the value of the β factor is affected by the plastic deformation beneath the indenter, five flat-ended punches were considered having the same R-ratio values between the diagonals (L and m) as for pyramidal indenters (Table 1).Also, the flat-ended punches were modeled using parametric Bezier surfaces and the finite element mesh is shown in Figure 4.This study allows an understanding of to what extent the evolution of the β factor values with R is related to the pyramidal geometry of the indenter and/or to the loss of symmetry of the flat punches, when R evolves from 1 to 7.11.
Using flat indenters, numerical simulations were performed up to a 0.025 µm displacement imposed so that only elastic deformation occurs in the material beneath the indenter.Five values of Young's modulus were used: 30, 200, 400, 600 and 800 GPa. Figure 11 shows the evolution of the load as a function of the indentation depth, h e , obtained in the numerical simulations with the flat indenters for values of the R-ratio equal to 1 (as for the Vickers indenter), 4 and 7.11 (as for the Knoop indenter).The figure shows that, for a given indentation depth, the load increases with the increase of the Young's modulus and the R-value.
As above mentioned, Equation (1) comes from the Sneddon's equation [20] that relates the applied load, P, with the elastic deflection of the surface of the material, h e , and can be written as follows: where E and ϑ are de Young's modulus and Poisson ratio of the material, respectively; a is the radius of the rigid Sneddon cylindrical flat indenter, or an equivalent value for other rigid indenter geometry with the same area; and β is a parameter that takes into account the geometry of the indenter (β = 1, for cylindrical flat indenter).The Young's modulus results obtained by the data such as in Figure 11 indicate that, in case of the flat punches, the value of the β parameter in Equation ( 6) is different from 1 and depends on the R-ratio, whatever the Young's modulus, as shown in Table 4.The values of β for R = 1 (as for the Vickers indenter) is in agreement with those obtained in previous studies (see, e.g., [22,26]).In the case of R = 7.11 (as for the Knoop indenter), β is equal to 1.372.The increase of the value of β with R is certainly related with the loss of symmetry of the flat punches with the increase of the ratio R. Indeed, the same was observed in case of the conical, Vickers and Berkovich indenters, whose values of the factor β increase in this order (1.034, 1.055, 1.081, respectively), although on a smaller scale [26].Figure 12 shows that the evolution of  with R for the flat punches is a quasi-linear function, quite close to that observed for the case of the pyramidal indenters (obtained with Equation (4): 1.054, 1.141, 1.232, 1.351 and 1.473).Moreover, the values of  for the pyramidal indenters are slightly higher than those determined for the flat indenters, except the R-ratio equal to 1.  Figure 12 shows that the evolution of β with R for the flat punches is a quasi-linear function, quite close to that observed for the case of the pyramidal indenters (obtained with Equation (4): 1.054, 1.141, 1.232, 1.351 and 1.473).Moreover, the values of β for the pyramidal indenters are slightly higher than those determined for the flat indenters, except the R-ratio equal to 1. Figure 12 shows that the evolution of  with R for the flat punches is a quasi-linear function, quite close to that observed for the case of the pyramidal indenters (obtained with Equation (4): 1.054, 1.141, 1.232, 1.351 and 1.473).Moreover, the values of  for the pyramidal indenters are slightly higher than those determined for the flat indenters, except the R-ratio equal to 1.The dissimilarity between the β values obtained for the flat and the pyramidal indenters is certainly because, in the case of pyramidal indenters, not only elastic deformation occurs, but also plastic deformation appears, which significantly distorts the material surface.Moreover, the increase of the value of R-ratio leads to the lowering of the symmetry of the plastic strain distribution along two axes of the indenter (see Figure 7).

Correlation with Experimental Results
In order to check the performance of the correction factor β proposed in this study, experimental results of five materials were used to calculate the Young's modulus.They are fully dense brittle materials, covering a wide range of mechanical properties, whose experimental data were obtained from bibliography [19].Table 5 presents the Young's modulus results of these materials determined by Grobal et al. [19], using two different methods proposed by: (i) Grobal et al. [19] (E G ); (ii) and Marshall et al. [31] and Riester et al. [18] (E M ).The results obtained by the method proposed in the current study (E), using β equal to 1.419 as defined above for the Knoop indenter, are also shown in this table.For comparison, nominal values of the Young's modulus, E nom , evaluated by ultrasonic method [19] are presented and used to calculate the errors.To evaluate the Young's modulus, Grobal et al. [19] make use of the plot of the contact stiffness, C, as a function of 1/ √ A C , whose values of the slopes of the fitted straight lines ( √ π/2βE r in Equation ( 1)), are also used in current study to assess the value of E, for all materials.The results in Table 5 show that the proposed β coefficient enables relatively good accuracy of the Young's modulus.The average of the absolute value of the error is equal to 4.27%.This is higher than that obtained by the method of Grobal et al. [19] and slightly lower than that achieved using the method by Marshall et al. [31].

Conclusions
A finite element study using the three-dimensional numerical simulation of hardness tests of elastic-plastic materials is performed.Pyramidal indenters with geometry from Vickers to Knoop were used.Also flat-ended punches were considered.This allowed us to obtain important information about the geometry of the indentation surface (sink-in and pile-up formation) and the distribution of plastic deformation beneath the pyramidal indenters.Both types of tests, with pyramidal and flat punch indenters, allow an assessment of the values of the geometrical parameter β to be used in the Sneddon's equation [20], for flat-ended punches, and the Doerner and Nix equation [21], for pyramidal indenters.This permits to determine the reduced Young's modulus of the indented material when using depth-sensing indentation (DSI) equipment.In the case of the Knoop indenter, the correction factor β of at about 1.419 is required, in order to accurately obtain the Young's modulus.This makes the procedure for analyzing the results of the Knoop indenter more expeditious than it currently is.
h c is the ideal indentation contact depth, θ 1 = 65 o and θ 2 = 86.25 o are the semi-apical angles of the Knoop indenter.

Figure 3 .
Figure 3. Evolution of the tangents of  and  as function of the R-ratio values.

Figure 3 .
Figure 3. Evolution of the tangents of  and  as function of the R-ratio values.Figure 3. Evolution of the tangents of θ 1 and θ 2 as function of the R-ratio values.

Figure 3 .
Figure 3. Evolution of the tangents of  and  as function of the R-ratio values.Figure 3. Evolution of the tangents of θ 1 and θ 2 as function of the R-ratio values.

Figure 4 .
Figure 4. Finite element mesh used in the numerical simulations: (a) global view; (b) detail of the central region where indentation occurs

Figure 4 .
Figure 4. Finite element mesh used in the numerical simulations: (a) global view; (b) detail of the central region where indentation occurs.

Figure 5 .
Figure 5. Surface indentation profiles obtained along the two diagonals, the shorter, m (z-axis) and the longer, L (x-axis), respectively: (a,c) obtained at maximum load; (b,d) after unloading.

Figure 5 .
Figure 5. Surface indentation profiles obtained along the two diagonals, the shorter, m (z-axis) and the longer, L (x-axis), respectively: (a,c) obtained at maximum load; (b,d) after unloading.

Figure 5 .
Figure 5. Surface indentation profiles obtained along the two diagonals, the shorter, m (z-axis) and the longer, L (x-axis), respectively: (a,c) obtained at maximum load; (b,d) after unloading.
shows the equivalent plastic strain distributions obtained at the maximum load in the numerical simulations of the materials M1 (ℎ /ℎ = 0.97) and M3 (ℎ /ℎ = 0.25) using the Knoop indenter.For comparison, the same figure also shows the same distributions obtained with the Vickers indenter (Figure 7b,d).

Figure 7 .
Figure 7. Equivalent plastic strain distributions obtained at maximum load in the numerical simulations using the Knoop and Vickers indenters: (a,b) Material M1; (c,d) Material M3.

Figure 6 .
Figure 6.Surface indentation profiles at maximum load, obtained from the results of Figure 5a,c, where z is multiplied by R = 7.11, in order to easily compare the profiles along the two diagonals, the longer, L, and the shorter, m: (a) Material M1; (b) Material 2; (c) Material M3.

Figure 6 .
Figure 6.Surface indentation profiles at maximum load, obtained from the results of Figure 5a,c, where z is multiplied by R = 7.11, in order to easily compare the profiles along the two diagonals, the longer, L, and the shorter, m: (a) Material M1; (b) Material 2; (c) Material M3.
shows the equivalent plastic strain distributions obtained at the maximum load in the numerical simulations of the materials M1 (ℎ /ℎ = 0.97) and M3 (ℎ /ℎ = 0.25) using the Knoop indenter.For comparison, the same figure also shows the same distributions obtained with the Vickers indenter (Figure 7b,d).

Figure 7 .
Figure 7. Equivalent plastic strain distributions obtained at maximum load in the numerical simulations using the Knoop and Vickers indenters: (a,b) Material M1; (c,d) Material M3.

Figure 7 .
Figure 7. Equivalent plastic strain distributions obtained at maximum load in the numerical simulations using the Knoop and Vickers indenters: (a,b) Material M1; (c,d) Material M3.

Figure 8 .
Figure 8. Normalized contact area results obtained in the numerical simulation of the materials with the values of the strain hardening parameter, yield stress and Young's modulus shown in Table 2, using the Knoop indenter.(a) Contact area  / ; (b) Contact area  / .

Figure 8 .
Figure 8. Normalized contact area results obtained in the numerical simulation of the materials with the values of the strain hardening parameter, yield stress and Young's modulus shown in Table 2, using the Knoop indenter.(a) Contact area A h c /A REF ; (b) Contact area A FE /A REF .

Figure 9 .
Figure 9. Normalized Young's modulus results obtained in the numerical simulation of the materials with the values of the strain hardening parameter, yield stress and Young's modulus shown in Table 2, using the Knoop indenter: (a) Young's modulus  / ; (b) Young's modulus  / .

Figure 10 .
Figure 10.Evolution of the ratio / as a function of / , obtained in the numerical simulations of all materials in Table 2, using the Knoop and Vickers indenters.

Figure 9 .
Figure 9. Normalized Young's modulus results obtained in the numerical simulation of the materials with the values of the strain hardening parameter, yield stress and Young's modulus shown in Table 2, using the Knoop indenter: (a) Young's modulus E h c /E REF ; (b) Young's modulus E FE /E REF .

Figure 9 .
Figure 9. Normalized Young's modulus results obtained in the numerical simulation of the materials with the values of the strain hardening parameter, yield stress and Young's modulus shown in Table 2, using the Knoop indenter: (a) Young's modulus  / ; (b) Young's modulus  / .

Figure 10 .
Figure 10.Evolution of the ratio / as a function of / , obtained in the numerical simulations of all materials in Table 2, using the Knoop and Vickers indenters.

Figure 10 .
Figure 10.Evolution of the ratio P/S 2 as a function of H/E 2 r , obtained in the numerical simulations of all materials in Table 2, using the Knoop and Vickers indenters.

Figure 11 .
Figure 11.Evolution of load as a function of the elastic indentation depth obtained in the numerical simulations using flat indenters with different ratio R: (a) R = 1; (b) R = 4; (c) R = 7.11.

Figure 12 .
Figure 12.Evolution of  as a function of the R-ratio obtained in numerical simulations using the five flat and pyramidal indenters with different values of the R-ratio.

Figure 11 .
Figure 11.Evolution of load as a function of the elastic indentation depth obtained in the numerical simulations using flat indenters with different ratio R: (a) R = 1; (b) R = 4; (c) R = 7.11.

Figure 11 .
Figure 11.Evolution of load as a function of the elastic indentation depth obtained in the numerical simulations using flat indenters with different ratio R: (a) R = 1; (b) R = 4; (c) R = 7.11.

Figure 12 .
Figure 12.Evolution of  as a function of the R-ratio obtained in numerical simulations using the five flat and pyramidal indenters with different values of the R-ratio.

Figure 12 .
Figure 12.Evolution of β as a function of the R-ratio obtained in numerical simulations using the five flat and pyramidal indenters with different values of the R-ratio.

Table 1 .
Area function of the indenters used in the numerical simulations.

Table 1 .
Area function of the indenters used in the numerical simulations.

Table 1 .
Area function of the indenters used in the numerical simulations.

Table 2 .
Mechanical properties of the materials used in the numerical simulations.

Table 3 .
Mechanical Properties of The Materials Used in the Study of The Knoop Indentation Geometry.

Table 4 .
Values obtained for the correction factor β in the numerical simulations with flat indenters.