Finite Element Simulation of Spherical Indentation Experiments

The response of linear elastic materials during indentation can be predicted by exact elasticity solutions as long as the indenter is described as a solid of revolution [1]. However, modelling the indentation process in the elasto-plastic region is a more complex problem. Since the constitutive equations are non-linear and a number of plastic properties must be included to describe material behaviour, an analytical solution is not easily obtained. As a result, understanding of elasto-plastic indentation behaviour can be better achieved through finite element (FE) simulation.

Computer simulations of indentation tests are often associated with efforts to develop material characterisation techniques based on indentation data [2–9]. In such context, FE modelling of conical, pyramidal and spherical indentations has been extensively used to simulate experiments and provide output suitable for the development of characterisation schemes applicable to materials of various constitutions and compositions. However, critical issues relating to the accuracy and reliability of the modelling have not always been discussed in great depth.

The subject matter of this paper is spherical indentations on ductile materials. The amount of published numerical work in this area is considerable. For this reason, brief reference is here made to a rather limited number of relevant publications. The aim is to highlight adopted tools and assumptions regarding materials and their properties.

General purpose software has been routinely used and the axisymmetry of the problem taken into account to reduce it into a two-dimensional one [2–3, 5–10]. The ball indenter has been modelled as rigid [2–3, 7–8, 10, 12] or elastic [2, 5–6, 9, 11]; the indented specimen elasto-plastic input was either obtained from tensile tests [8, 13] or conventionally built from known material parameters [3, 5–7, 9]; friction between indenter and specimen was sometimes ignored [3, 7, 9, 12] but more often accounted for [2, 5, 6, 8, 10, 13]. Various mesh schemes were applied to test and optimise the performance of the modelling [2, 9, 10].

The development of an FE analysis of the problem is here described in some detail. Particular attention is given to contact modelling as well as the shape, size and distribution of elements. The model was validated by comparing its predictions to exact elastic solutions as well as experimental results from stainless steel and carbon steel samples, whose stressstrain curves had been obtained by tensile testing. These curves were used as FE input to simulate the process of indenting these materials by a sphere of 150 mm nominal radius. Friction was found to have significant influence on pile-up. It was also shown that the agreement between numerical and experimental results could be enhanced if the creep occurring during indentation is taken into account.

Three types of steel were tested and analysed. They are two grades of mild steel, 43 A and 50 A, and a stainless steel, grade 316 L, supplied in plates of 8 mm thickness. Different grades of steel were considered in order to cover a range of hardness. Steel chemical composition is detailed in Table 1.

Sections taken from the supplied plates were metallographically prepared and examined by optical microscopy. The rolling direction was determined by observing the orientation of inclusions. The grain size was found equal to 10 to 15 μm. This was suitable as several grains were thus included within an indent, which might typically be 70–120 μm in diameter. Thus, the indentation data reflected the average behaviour of the bulk metal as determined by a tensile test, rather than the behaviour of a single grain.

Tensile tests provide reliable and realistic information on elasto-plastic characteristics, which is essential input to the planned FE analysis. Tensile test samples with a circular cross-section were machined from the as-received plates according to the ASTM specifications. They were cut in their rolling direction with a gauge length of 25 mm. Tensile testing was carried out using an Instron machine with the crosshead speed set at 1 mm/minute. The elongation was measured using an extensometer that can cover up to 40% strain. Elongation beyond this value was measured from the crosshead motion.

For a more accurate measurement of load and elongation, a digital data logger was connected to the Instron tensile tester for calibration and the test itself. The data produced by the Instron, corresponding to force and elongation values, were converted into engineering stress-strain points. Then, the true stress-logarithmic strain $(\sigma-\varepsilon)$ relation, which is the final information sought from tensile testing, was obtained using the equations

where $\sigma_0$ and $e$ are the engineering stress and strain values, respectively.

The material was assumed to have the same elasto-plastic behaviour in compression as in tension. Thus the $\sigma-\varepsilon$ curve obtained from the tensile testing could be used as the input into the FE simulation of indentation, which causes a predominantly compressive stress field. Since the strains developing in an indentation test are expected to be higher than the fracture strain in tensile test, the $\sigma-\varepsilon$ curve was extrapolated over a range of higher strain values.

The maximum calculated logarithmic strain is at the engineering strain corresponding to maximum load where necking commences. The strain up to this point is considered uniform, but beyond that the strain becomes localised around the neck region hence its measured value depends on the considered gauge length, in which case eqns (1) and (2) are not valid. To determine the $\sigma-\varepsilon$ curve at higher strain values, the power law function

was fitted to the $\sigma-\varepsilon$ curve; in eqn (3), $K$ is the strength coefficient and $n$ the strain hardening exponent. It was not always found possible to fit a single power law function to the whole post-yield region of the curve. In such a case, this part of the curve was split into two and a different function was fitted to each of them. Typical engineering and the true stress-logarithmic strain curves with their extrapolated parts are plotted in Figure 1 and Figure 2.

The moduli of elasticity were taken equal to 210 and 185.1 GPa for mild steel and stainless steel, respectively. The values of material parameters for all specimens are summarised in Table 2. Although a single power law function is well fitted to the stress-strain curve for plates S6 and S7, S5 showed different plastic behaviour despite being of the same grade. Hence, obtaining the stress-strain curve for each plate is necessary for a precise description of material behaviour in FE simulations.

The indentation process was simulated by developing an elastic-plastic FE model using a version of ANSYS [14], a general purpose type of analysis software. Isotropic hardening plasticity and the von Mises yield criterion were assumed. Large strain and other non-linear features of ANSYS were adopted. The material data input for each specimen was the corresponding true stress-logarithmic strain curve. The diamond spherical indenter tip was modelled as both rigid and linear elastic; in the latter case, its modulus of elasticity and Poisson’s ratio were taken equal to 1141 GPa and 0.07, respectively. These values were provided by the manufacturers of the instrumented indentation tester.

The simulation of indentation process is a 3D axisymmetric problem. It assumes a spherical indenter tip pressed into a cylinder of material that should be representative of a homogenous infinite half space. Loading the indenter is prescribed as applied displacement on its top surface. Running a 3D problem is expensive and requires large memory and disk space. Advantage can be taken of the axisymmetry of the problem, so that it can be simplified to a 2D problem instead, and only half of the indenter and the material cylinder section by a plane through the axis of symmetry are modelled. This reduces the number of elements and the number of degrees of freedom per node; hence, the running time and disk space required.

A four-node quadrilateral element was used to mesh the deformable indenter and the sample. Each node of this element has two degrees of freedom.

The free mesh option in ANSYS was initially used to mesh both the indenter and the sample. Over the region where the two bodies are expected to be in contact, a regular mesh was mapped on both bodies for generating elements of the same size. This mesh regularity was expected to improve numerical performance and thus enhance accuracy over this region of interest. However, the indentation process involves high geometrical non-linearity. With increasing penetration depth, the induced strains increase and element distortion increases excessively and reaches a stage when the solution fails to converge even at shallow indentation due to violating element geometry limits.

The automatic free and mapped mesh options in ANSYS create a large number of elements with little control on their shape and interior angles that deviate from 90°. The steeper the transition from a fine mesh, in the contact region, to a coarse one, at the far edges, the greater the irregularity in element shape.

Decreasing the element geometric irregularity by introducing a smooth variation of mesh density produces a mesh with larger number of elements, hence longer running time and convergence problems. This necessitates the imposition of more control on all elements of the mesh, which can only be achieved by generating the mesh manually. Whilst this task is very complicated and tedious to do using ANSYS, another type of modelling software, PATRAN [15], offers a powerful facility of manual mesh generation with full control on all elements size and shape. In addition, it has the capability of doubling the element size during the transition from the fine mesh zone to the coarse mesh one without violating any element geometry limits. This is done by using trapezoidal elements as interlinks between the twosize zones, as illustrated in Figure 3. As a result, fewer elements are created for the whole mesh compared with that of ANSYS.

Although this mesh decreased the running time for every load step, convergence difficulties persisted at deep penetration. This is due to the high strains imposed on the elements within the indentation zone of the material, especially those closer to the contact centre, which undergo excessive deformation and almost get squashed as a result. To counteract this excessive change in the aspect ratio, the initial element geometry of the specimen was modelled as rectangle instead of square. The aspect ratio was taken 2.0 with the longer dimension in the direction of indentation. This helped the model reach the sought penetration depth without any convergence problems, thus this kind of mesh was adopted for all subsequent FE simulations.

Contact between the specimen and the indenter was modelled using contact elements generated over the region of possible contact between the two surfaces. The contact element, initially used to represent node-to-surface contact, was a triangle with the base being the target line between two nodes on one of the surfaces and the opposing vertex being the contact node on the other surface. The indenter surface, which should be convex, is considered ‘contact’ and the flat sample surface is the ‘target’.

Contact between two bodies is detected when virtual penetration occurs within a certain limit. A compatibility control method is used to ensure that one surface does not penetrate into another surface by more than an acceptable tolerance. The node-to-surface model uses as control parameter a normal contact stiffness, whose optimum value can be ascertained by trial and error, that is, by testing the performance of the model for a range of its values.

It was noted that, with this contact model, the solution was extremely slow to converge with the number of iterations per load step increasing considerably after only a few microns of simulated indenter penetration. The convergence problem became even worse when sliding friction was included between the two surfaces. This prompted the adoption of an alternative contact modelling, according to which contact is represented by a pair of contact elements, one of them used to represent contact and sliding between the target surface, that of the indenter, and a deformable surface, that of the specimen. The geometric description of the two elements is illustrated in the diagram of Figure 4.

The specimen top surface was overlaid with the contact element and the indenter surface was meshed with the target element. The penalty method with Lagrange multiplier was used to ensure that penetration compatibility is satisfied. Convergence issues were resolved with this modelling, which was subsequently adopted in all simulations producing the results presented in this paper.

Initially, the mesh of the indented specimen was generated with its refined part in the indentation zone made up of element size 1.0 × 2.0 µm². The dimensions of the mesh that covers the modelled half of the cylinder were 672 × 448 mm². This is referred to as mesh A and is illustrated in Figure 5.

To investigate the sensitivity of the results to mesh changes, two other meshes, labelled B and C, were generated by enlarging the specimen model dimensions and the refined zone covering the indentation region. In all three meshes, the indenter is meshed using the free mesh option in ANSYS. Elements on the indenter perimeter from the tip up to the height that corresponds to the expected indentation depth are taken to be 1.0 µm wide, which is the same width as the corresponding elements of the specimen over the contact region. The characteristics of all meshes are summarised in Table 3.

The developed FE model was first applied to simulate the indentation process of a perfect sphere with a radius of 150 µm on a linear elastic material. It was thus possible to validate the model against Sneddon’s elastic solution of a sphere pressed on an elastic flat surface [1]. The contact was modelled as frictionless. Loading was prescribed as displacement controlled increased by 0.1 µm at every load step. Peak load was attained when maximum depth reached 20 µm. At every load step the indent profile geometry, contact radius and depth were directly measured and compared with the corresponding predictions of Sneddon’s solution. The predicted load-indentation depth curve by FE for each mesh A, B and C, described in Table 3, as well as Sneddon’s solution are shown in Figure 6.

As observed in the plots of Figure 6, best agreement between numerical predictions and Sneddon’s elastic solution is achieved with meshes B and C. The same conclusion is reached by plotting the ratio of contact depth, $h_c$, to indentation depth, $h_t$, which should be equal to 0.5 according to Sneddon’s solution. Since there is no significant difference between the results produced by these meshes and, as can be seen in Table 3, mesh B has less than half the number of elements than mesh C, mesh B was adopted for the FE simulations of contact discussed in this paper.

The initial elasto-plastic indentation simulations were carried out in order to assess the effect of variations in critical input parameters on the results. For this purpose, idealised material properties and a perfect spherical indenter of 150 $\mu\mathrm{m}$ radius were adopted. The stress-strain curve produced by tensile testing of specimen S5 was modified to an idealised curve with a single value for the strain hardening exponent and the strength coefficient. The modified stress-strain curve was based on the actual curve for steel plate S5. It corresponded to a material with Young's modulus $E$ =210 GPa, elastic limit $\sigma_e=340 \mathrm{MPa}, \sigma_{0.2}=375 \mathrm{MPa}, n=0.132$, and $K=782.91 \mathrm{MPa}$, satisfying the continuity condition at yield point where

The part on the stress–strain curve joining the elastic limit to the yield stress was created using a quadratic function. The modified stress-strain curve was very close to the original one, characterised by three different regions.

Load-indentation curves were obtained for all three meshes detailed in Table 3 and it was noted that the elasto-plastic solutions almost overlap each other. Sensitivity to mesh variation was therefore found negligible. This suggested that mesh B would be satisfactory for further elasto-plastic FE simulations.

Contact tolerance and contact stiffness were two arbitrarily specified numerical parameters associated with contact modelling. Significant changes to their recommended values had no noticeable effect on the indentation output.

Next, the effect of friction between indenter and specimen was examined. Five different values of the friction coefficient of contact were tried, namely 0, 0.05, 0.1, 0.2 and 0.5. The resulting P–h curves can be compared by reference to Figure 7.

As noticed in Figure 7, the curves start to deviate at around two-thirds of the maximum depth and the maximum difference at the peak load is less than 4%. However, when the pile-up was assessed, the effect of friction was found to be more pronounced. The profiles of the indent for the five friction coefficients can be compared by referring to Figure 8.

Finally, the effect of changing the elastic modulus of the indenter was assessed. Due to anisotropy of diamond, its Young’s modulus, E_i , varies with crystal orientation. E_i can be determined by indirect methods. One such method is by measuring the speed of sound or compression waves in a diamond crystal. An approximate value of 1050 GPa is given for E_i by Wilks and Wilks [16], which was said to vary with the crystal orientation by approximately 10%. This value is very close to that initially adopted.

To assess the effect of any variation of E_i on the P-h curves resulting from FE simulation, the latter was run for four values of Ei ranging from 800 to 2000 GPa in addition to the case of an infinitely rigid indenter. Comparison of the resulting load-depth curves showed that there is no significant difference between them for the extreme cases of Young’s modulus values, and their difference is even more negligible for the realistic range lying between 800 and 1500 GPa. As a result, the adopted value of 1141 GPa was proved to be acceptable.

The performance of the developed FE model in the region of elasto-plastic deformation was validated through comparisons of its predictions with experimental data obtained by instrumented indentation testing [17]. FE simulations were performed on plates S5, S6, S7 and stainless steel. The actual properties of these specimens, described in Section 3, were entered as material input to the elasto-plastic FE analysis. The indenter model was changed from perfect sphere to the real profile determined by digitising and averaging SEM photographs of the actual indenter tip [18]. The program was run with friction coefficient m = 0.2. Plots of both experimental and numerical results for plates S5 and stainless steel are shown in Figure 9 and Figure 10, respectively.

The good agreement observed in these figures was also noted in similar plots for the other two steel plates. The observed variation in the experimental curves is considered to lie within acceptable limits and it can be understood to be due to variations in material microstructure, surface finishes, and equipment error, among others. The experimental results appear to indicate a stiffer material than that anticipated from the predictions but it should be noted that the FE simulation does not consider any creep effect. Taking that effect into consideration, as discussed in the next section, brings the experimental points closer to the corresponding FE predictions.

The reliability of any output from indentation experiments depends on ensuring ideal measurement conditions, which help eliminate any undesirable effects. Thus, the recorded indentation data would correspond to the material response predicted by the FE simulation. One effect that cannot be eliminated completely is that of loading rate, which manifests itself as creep; the latter is observed as a plateau in the experimental load-indentation curve at peak load. The extent of creep depends on the material and it is also influenced by the hold/dwell time, which is the duration of maximum load application.

If creep is not considered, this can have a significant effect on the calculation of hardness and the reduced modulus from the unloading curve. It was found by Chudoba and Richter [19] that creep not only affects the measured depth but also the slope of the unloading curve at peak load, which is a key quantity in calculating the reduced modulus. It was reported that when creep was not allowed to take place, the resulting unloading curve showed bowing towards higher depth.

Noticeable creep deformation was reported during indentations on tungsten and aluminium using the Berkovich indenter [19–21]. Hence, it was considered essential to study the effect of creep using a spherical indenter on mild steel in order to eliminate its effect from any subsequently developed material characterisation scheme. For this purpose, spherical indentations were performed at three locations on steel plate S5. The peak load was chosen to be the maximum applicable by the instrument, that is, 20 N.

Since the extent of creep depends on the loading rate, its effect was examined by varying the loading time, that is, the time required to reach the peak load. Three values for the loading time, namely 8, 40, and 200 s, were adopted. When the peak load was reached, it was held constant for a dwell time of 180 s then partial unloading by 80% followed at the same loading rate. The obtained load-indentation curves are shown with the corresponding FE prediction in Figure 11.

The capability of continuously recording the creep deformation with time made it possible to monitor the development of such deformation over the dwell time. This is shown for the various loading times in Figure 12. As can be clearly seen in that figure, the creep effect can be removed by holding the peak load constant for a sufficiently long period. Creep deformation takes place at a high rate during the first few seconds then it diminishes within less than a minute. Its value reaches 150 nm/s at the first few seconds of the hold time and reduces to less than 3 nm/s after 15–20 s.

It can also be observed that the loading time has a noticeable effect on the three main segments of the load-displacement curve, namely, the loading curve, hold time at peak load (amount of creep), and unloading curve. In the case of loading time of 8 s, the loading curve appeared stiffer and the amount of creep deformation was large, that is, 1.387 µm. On the other hand, a more compliant loading curve and less creep were noticed when loading time was longer. In the case of loading times of 40 and 200 s, the corresponding measured creep was 0.6894 and 0.4983 µm, respectively.

Even though the long loading time results in less creep, allowing enough dwell time makes the depth at peak load reach the same final value and the unloading curves almost overlap regardless of the loading time as long as it is not very short, that is, at least 40 s.

Applying the required dwell time in the case of the load time of 8 s, will bring the maximum depth close to the corresponding value of the other two cases. Nevertheless, the unloading curve shows bowing at peak load, which will incur an error in calculating the contact stiffness, and thus the modulus. This indicates that short-time loading causes erroneous indentation data, which they seem to include undesired deformations. This explanation of bowing can be supported by the findings of Chudoba and Richter [19] who reported this bowing phenomenon when dwell time is zero and loading time is 92 s. The material behaviour in the case of short-time loading can be understood as a dynamic response to high strain rate loading, whereby the material seems stronger, compared with its response to static loading.

Referring to Figure 9 and Figure 10, it can be observed that the creep value is not only dependent on the loading rate but also on the depth reached, which is in agreement with the findings of Chudoba and Richter [19]. Figure 12 shows that a dwell time of 120 s is enough to allow creep deformation to take place in full and it can compensate for the value of the loading time, as long as the latter is 40 s or more. It also shows, together with Figure 11, that a loading time of 200 s is adequate to produce reliable indentation data provided enough dwell time is allowed, and that increasing the loading time will not result in noticeable change in the amount of creep. Therefore, a value of 120 s for the dwell time and loading time of 200 s were adopted in the subsequent experiments for material characterisation.

This creep phenomenon clearly implies that since the loading rate effect was not considered in the FE model, only those points corresponding to the end of the creep plateau of P–h curves should be considered when comparison with corresponding FE predictions is made, as the other points on the loading curve exclude creep, thus are expected to represent a stiffer response.

An axisymmetric FE model was developed for simulating the elasto-plastic spherical indentation of test materials with known, experimentally obtained stress-strain relations. The specimen mesh was generated manually in order to address the high geometric nonlinearity induced by deep indentation. Mesh sensitivity was carried out to ensure the capability of the adopted mesh with regards to accuracy and efficiency.

The FE model was first validated against the Sneddon’s elastic solution for a sphere onto a flat surface. Then it was applied to the indentation of idealised elasto-plastic materials, in order to assess the effect of various FE control and physical parameters. These included: the contact tolerance and contact stiffness, friction coefficient between indenter and specimen and the elastic modulus of the diamond indenter. It was found that none of these parameters, within their range of variability and load range of application, had a noticeable effect on the simulated P–h curves. However, friction was shown to have a pronounced influence on the predicted pile-up, thus the contact area, which is an important parameter in characterisation procedures based on spherical indentation data.

Time-dependent deformation, so-called creep, was also investigated. The experimental conditions influencing creep were examined. It was found that loading rate and the hold/ dwell time at peak load govern not only the extent of creep but also the initial contact stiffness. It was also noticed that the higher the loading rate, the larger the creep effect, thus a longer dwell time will be required to allow for this deformation to take place as long as the full load is not applied abruptly, i.e. in less than 10 s, which can affect the initial unloading contact stiffness. Based on the experimental observations on the test materials, it was found adequate to adopt, for practical purposes, a loading time of not less than 40 s followed by a dwell time of not less than 100 s to account fully for any time-dependent deformation.

Load-indentation curves were simulated using the actual imperfect geometry of the spherical indenter and compared with those obtained from experiments. Indentation experiments were performed in a gradually amplified cyclic loading pattern and the peak load was held for a specified time to remove the effect of creep deformation. The points corresponding to the end of the creep plateau of the experimental P–h curves were found in good agreement with the corresponding FE predictions.