Determination of elastic and viscoplastic material properties obtained from indentation tests using a combined finite element analysis and optimization approach

Conventional indentation tests do not provide an accurate estimation of viscoplastic material properties. In this work, a combined finite element analysis and optimization approach is developed for the determination of elastic–plastic and creep material properties using only a single indentation loading–unloading curve based on a two-layer viscoplasticity model. Utilizing the indentation loading–unloading curve obtained from a finite element-simulated experiment with a spherical and a conical indenter, a set of six key material properties (Young’s modulus, yield stress, work hardening exponent and three creep parameters) can be determined. Non-linear optimization algorithms are used with different sets of initial material properties, leading to good agreements with the numerically simulated target loading–unloading curves.


Introduction
Instrumented indentations at micro or nano-scales have become established techniques for measuring the mechanical properties of a large variety of materials. The Oliver and Pharr approach [1] is the most well-known method for determining the hardness and elastic modulus of the test material. This method involves the determination of the mechanical properties of the surface of a given material from loading-unloading curves obtained from micro-or nano-indentation tests, based on the assumption that the material behaves in an elastic-plastic response during the loading phase and the unloading behaviour of the indenter is fully elastic with no plastic deformation [1]. Various approaches have been proposed using dimensionless analysis and numerical optimisation to extract material properties of power law materials with good estimations [2][3][4][5][6][7][8][9].
Indentation creep tests [10][11][12][13][14][15] have been used to obtain the creep properties of materials, despite the fact that uniaxial tensile creep tests are the standard techniques to obtain creep parameters. There are several advantages of indentation creep tests. For example, only a small amount of the material is needed and the test can be used for the characterization of the local deformation behaviour. In general, creep is usually used to describe a time-dependent material behaviour of metals at high temperature that is a result of visco-elastic deformation when stress or strain is applied. Elasticplastic behaviour usually refers to elastic and time-independent plastic deformation, whereas viscoelastic and viscoplastic behaviours refer to time-dependent elastic and time-dependent plastic deformations respectively.
Recent studies [16][17][18][19][20] have found that the determination of material properties from time-dependent material behaviour based on conventional indentation tests, based on the Oliver-Pharr method, does not provide an accurate estimation of material properties. The contact between an indenter and a material specimen is visco-elastic and not purely elastic, in which creep occurs during the instrumented indentation unloading which leads to an overestimation of Young's modulus.
Research efforts [16][17][18][19][20] have focused on extracting the time-dependent mechanical properties, which are limited to viscoelastic materials and do not consider the plastic deformations of the materials, but some researchers [21][22] have assumed that plastic-viscoelastic procedures occur separately in an indentation test. Tweedie and Van Vliet [20] stated that the plastic deformation can be negligible as the indentation is relatively shallow, hence a purely viscoelastic behaviour can be analysed for the time-dependent response. However, plastic deformation may not be negligible, since time-dependent indentation tests may involve very high localised contact stresses resulting in plastic strains.
The main objective of this study is to determine the elastic-plastic and creep material properties from indentation loading-unloading curves using a Finite Element (FE) approach combined with optimization algorithms for a combined two-layer viscoplasticity material model available in the ABAQUS FE code [23]. The authors' previous optimization approach [9], which was focused on the evaluation of elastic-plastic material properties from instrumented indentation loading-unloading curves using sharp indenters, is extended in this study to elastic-plastic and creep material properties using a spherical indenter. The current investigation builds on previous studies evaluating elasticplastic material properties from indentation loading-unloading curves, using a novel two-layer viscoplasticity model and combined FE and optimization methods to arrive at the mechanical properties of elastic-plastic and creep material to within an error of less than 10%.

Typical indentation loading-unloading curve
The Oliver-Pharr method [1], which is the most well known method for the interpretation of indentation tests, is usually used to determine the hardness and elastic modulus of the test material. A typical loading-unloading curve is shown in Fig. 1, where h m is the maximum depth, h f is the final depth after the indenter is fully unloaded, dP dh is the initial slope of the unloading curve and P m is the maximum indentation load. In general, the main parameters, shown in Fig. 1 have been used to determine the hardness and elastic modulus of specimen from loading-unloading curves.

Two-layer viscoplasticity model
The FE analysis of the bulk material indentation is based on an axisymmetric indenter that is modelled using the ABAQUS Standard FE code. The two-layer viscoplasticity model within ABAQUS [23,24] is chosen as an example to demonstrate both creep and elastic-plastic material behaviour occurring in the indentation test. A two-layer viscoplasticity model is developed for modelling materials in which both significant time-dependent and time-independent behaviour are observed, which for metals typically occurs at elevated temperatures. A one-dimensional idealization of the two-layer viscoplasticity model is presented in Fig. 2, which describes the combined effect of a rate-independent (elastic-viscous network) and a rate-dependent (elastic-plastic network) material behaviours. It is noted that the rate-independent behaviour exhibits permanent deformations after the load application, whereas the rate-dependent behaviour exhibits permanent deformation of the material under load over time. The model consists of an elastic-plastic network that is in parallel with an elastic-visco (Maxwell model) network, where is the elastic modulus of elastic-plastic network, is the elastic modulus of elastic-viscous network, is the initial yield stress, ′ is the power law hardening with work hardening exponent, 1 , and A and 2 are the Norton creep parameters (based on the Norton Law: creep strain rate = A 2 ).
The elastic-plastic network predicts the time-independent behaviour of the material, whereas the elastic-viscous network predicts the time-dependent behaviour of the material. The two-layer viscoplasticity model is based on the von-Mises yield condition in the elastic-plastic network and the Norton power law (secondary creep) for the viscoplastic behaviour in the elastic-viscous network.
The two mechanisms are assumed to be independent, and the total stress σ is the sum of the stress σ v in the elastic-viscous network and the stress σ p in the elastic-plastic network. In this study, the twolayer viscoplasticiy model is combined with a power-law strain hardening for the time-independent , ′ behaviour and the viscoplastic behaviour of the material is assumed to be governed by the Norton Law, also known as the Norton-Hoff law (secondary creep).
The material behaviour in the two-layer viscoplasticity model in ABAQUS covers elastic, plastic, and viscous deformations. The elastic part of both networks in For many engineering materials, the plastic behaviour can be closely approximated by a power law description, which is presented schematically in

Fig.3. Power law elasto-plastic stress-strain behaviour [4]
A simple elastic-plastic, true stress-true strain behaviour is assumed to be Where is the initial yield stress, 1 is the work hardening exponent and the coefficient R can be expressed as where is the strain at the initial yield stress and is the plastic strain. The viscous behaviour of the material is assumed to be governed by the Norton Law, also known as the Norton-Hoff law where is the viscous stress in the viscoelastic network and A and 2 are Norton constants. It is assumed that the mechanisms are independent and can be written as: Therefore, the elastic strain is defined as: The total strain comprising elastic, plastic and viscous strains can be expressed as follows: where = is the elastic strain in the elastic-plastic network and = is the elastic strain in the elastic-viscous network. In the ABAQUS input file, a discrete set of points is required to represent the inelastic stress-strain behaviour, which is calculated based on Eq.

Optimization Model
A non-linear optimization technique is developed using the MATLAB optimization toolbox, which can provide an interface to FE codes such as ABAQUS, through various programming languages such as C and Python. The optimization technique is used to determine the mechanical properties for a given set of target indentation data using an iterative procedure based on a MATLAB nonlinear least-squares routine to produce the best fit between the given indentation data and the optimized indentation data, produced by FE analysis. This non-linear least-squares optimization function (called LSQNONLIN) is a subspace trust-region-reflective algorithm and is based on the interiorreflective Newton model [25]. This optimization procedure minimizes the objective function, and iterations are performed until convergence is reached. The optimization model can be written as follows: where F( ) is the objective function, is the optimization variable set (a vector in the n-dimensional space, ), which for this specific case contains the full set of the material constants in the model, as follows: { 0 < < 300; 0 < < 2; 0 < 1 < 5; 0 < < 10; 0 < 2 < 10; 0 < < 10; Since the indenter is load-controlled, D( ) i pre and D i exp are the predicted total displacement and the (experimental) displacement from target data, respectively, at a specific position i, within the loops.
N is the total number of points used to represent the experimental (measured) load-displacement curves. Arbitrary values of ( , , 1 , , 2 , ) have been chosen as initial values and the proposed optimization algorithm has been used to find the optimised values of these parameters from which the best fit between the experimental and predicted load-displacement loops can be achieved.

Optimization Procedure
The general optimization algorithm used in this work is illustrated in Fig 5. Since the initial guess values for ( , , 1 , , 2 , ) are provided, the optimization procedure is carried out in several steps using MATLAB, which controls the C language EXE file to automatically generate an ABAQUS input file, running ABAQUS and a Python script to extract the load history from the resulting ABAQUS output file. In terms of pre-processing the FE analysis, the material properties in the ABAQUS input file are replaced by new two-layer viscoplasticity material properties. In ABAQUS, these are Young's modulus, Poisson's ratio, and discrete points on the post yielding true stress-true strain curve.
In the ABAQUS input file, a discrete set of points is required to represent the uniaxial stress-strain data, rather than specifying the work-hardening exponent 1 [8]. Therefore, a fixed set of plastic strain values of 0.005, 0.01, and 0.0115 are used in order to specify the plastic stress-strain data in ABAQUS. The coefficient R can be calculated by using Eq. (3) and the updated stress data related to these strains can be obtained by Eq. (2). The viscous behaviour of the material is governed by the Norton power-law with creep parameters, A and 2 and the fraction, f, is defined in the viscous section of ABAQUS input file. This procedure can be performed by a C language code or a similar computing language to replace the current material properties in the ABAQUS input file by the new calculated material properties. In terms of post-processing, the load history results, extracted by a Python script, are read by a MATLAB program and the objective function calculated. All the procedures are processed automatically until convergence is reached.  Based on the procedures in Fig. 4 Yes No

Finite Element Indentation modelling
The FE analysis of the bulk material indentation system is based on axisymmetric elements using ABAQUS. Contacts between three different types of rigid indenters and an isotropic two-layer viscoplasticity specimen are modelled. During each iteration, FE analysis is performed with the updated two-layer viscoplasticity properties determined from the optimization processes. The tip radius of the spherical indenter is 0.1 mm. For shallow indentation depths, the size effects in the reallife experimental indentation tests can affect the accuracy of the simulations [27]. It should be noted that the FE simulations do not model the indentation size effects and are therefore limited to simulating macro indentations. The friction coefficients at the contact surfaces between the indenter and the top surface of the bulk material are assumed to be zero, since friction has a negligible effect on the indentation process [4]. It is assumed that there is no temperature variation of the bulk material during FE analysis. In the case of a conical indenter, a perfectly sharp indenter tip is used.
A "master-slave" contact scheme in the FE procedure is applied on the rigid indenter and the specimen surfaces. The depth of the bulk material is 2 mm and the maximum load on the indenter is 4.62 N under load-control conditions. The entire processes have been performed by a PC running Window XP with Intel Core 2 Duo CPU E8300 processor.
For convenience, the indented specimen is modelled as an axisymmetric geometry with four node bilinear axisymmetric quadrilateral continuum elements (CAX4 in ABAQUS). For the indenter, an axisymmetric analytical rigid shell/body is used. The region of interest is in the vicinity of the indenter surface and a high element density has been used there due to the expected high stress gradients immediately beneath the indenter tip, whereas a gradually coarser mesh further from the contact region is used, as shown in Fig. 6. All nodes at the base of the specimen are constrained to prevent them from moving in the x and y directions. The simulation is carried out in three distinct steps: loading, holding and unloading. In the first step, a total indentation load 4.62 N is applied.
During the loading step, the rigid indenter moves downwards along the y-direction and penetrates the foundation up to the maximum specified force. In the second step, the indenter is held at the maximum specified force with a dwell time (3 s) to induce viscoelastic deformation. In the third step, the load is reduced to zero. In the first unloading step, significant nonlinearity occurs which requires very small load/time increments. In the second and third steps, contact between the indenter and substrate is maintained and removed, respectively and larger load increments can be used.

Optimization using a target curve obtained from a FE simulation
The optimization scheme has been applied to determine the material properties of two-layer viscoplasticity model using a spherical indenter. Firstly, a set of material properties are chosen as the target values and FE analysis is performed to obtain a simulated target loading-unloading curve. A set of initial 'guess' material parameters are then selected and implemented in the MATLAB optimization algorithm which automatically performs a new ABAQUS run for each iteration. The target material properties are shown in Table 1 for two materials, XN40F (a high nickel-chromium material) at 900°C and P91 steel at 600°C. The given material properties have been obtained by uniaxial tensile creep testing. Previous work, see. e.g. [28,29], has shown that creep parameters obtained from impression creep agree well with those obtained from conventional uniaxial creep testing.
In order to check the sensitivity of the optimization algorithms, each parameter is changed in turn, while all other parameters are fixed at their target values. Generally, the optimized results are obtained in about 8-10 iterations with a deviation of less than about 1% from the target values.
Moreover, optimizations based on a combination of parameters involving two or four parameters are performed, as shown in Table 2. The parameter 'Errnorm' is the sum of the squares of the differences between the target and optimised curves. The results show that all the variables converge from their initial guess parameters to their target values to within 1% and 'Errnorm' is nearly zero, while it takes more iterations to reach the target values when more parameters are added, except in Test 4 in Table 2. It is noted that convergence is faster, and with improved accuracy, when the initial guess values are chosen closer to the target values. These results demonstrate that the proposed optimization algorithms are capable of obtaining the material properties accurately. Fig. 7 shows the convergence history of the material properties during the iteration process. It clearly demonstrates that convergence to the target values can be achieved despite a large variation in the initial guess values. As can be seen from Fig 7 (a), convergence starts after about 6 iterations for Young's modulus and creep parameter n 2 , whereas the creep parameter A goes up and then steadily decreases until the target value is reached. The four parameters converge to their target values after 9 iterations in Fig 7 (b), which is much faster than in Fig 7 (a).
(a) Test 3 in Table 2 (b) Test 4 in Table 2

Fig. 7 Optimized parameter values versus iterations for a spherical indenter (a) Test 3 and (b) Test 4
After checking the sensitivity of the optimization algorithms, the optimizations of the full combination of six parameters of the XN40F material are investigated. results, whereas the yield stress in Test 2 and creep parameter A are generally higher than the other parameters. Fig 8 shows the convergence trends for all three tests for XN40F. Although the initial guess values are different, it is interesting to see that the trends shown in Fig 8 (a), (b) and (c) for the six parameters are similar, and gradually reach their target values.  Table 3 for XN40F material

Fig. 8. Optimized parameter values versus iterations for a spherical indenter (a) Test 1, (b) Test 2 and (C) Test 3 in
The combination sets of six parameters optimization for another material, P91 steel, are also investigated. Table 4 shows the details of the initial values and final optimised values for the P91 steel material. Most of the final optimised parameters in each case converge to within 10%, but higher errors occur in the creep parameter A. Fig 9 shows the convergence history of the material properties for each iteration which clearly illustrates that convergence to the target values can be achieved despite a large variation in the initial values. In Test 1, where the differences between the target and initial values are small, the convergence trend is similar for all material parameters, whereas the convergence trends are different in Tests 2 and 3. It is interesting to note that the convergence trend of the work hardening exponent in Fig 9(c), goes down and then steadily increases until the target value is reached. As before, the convergence rate and accuracy depend on  the initial guess values. Since this is a non-linear material behaviour, there is no guarantee that the optimization algorithm will always converge to the 'target' parameters (whether obtained by FE analysis or an experiment). The optimization approach produces impressive accuracy in Table 2,   whereas less accurate optimised results with six parameters are shown in both Tables 3 and 4.
Further studies should be undertaken to analyse the parameter correlation in terms of six parameters.
It is clearly illustrated that the convergence accuracy of creep parameter A in Tests 2 and 3 in Table   4 is not as good as the other parameters. Therefore, the sensitivity of the loading-unloading curves to changes in the creep parameter A is investigated further in Fig 10. Fig 10 (a) shows the loadingunloading curves based on the optimised and target values of the creep parameter A, based on the results of Test 3 in Table 4. Figure 10 Table 4. It is interesting to note that changes of the creep parameter A of up to 9% have a very small influence on the loading-unloading curves.  Table 4 (b) Test 2 in Table 4 (c) Test 3 in Table 4 Fig. 9 Optimized parameter values versus iterations for the P91 material using a spherical indenter

Optimization approach using a conical indenter
A previous study [8] has shown that Berkovich and Vickers indenters displace more volume and thereby produce greater local stresses due to fact that the contact areas between the indenters and the bulk materials are larger than in the case of conical indenters. Despite these differences, conical indenters have the advantage of possessing axial symmetry and equivalent projected areas of contact can be used between conical and pyramid-shaped indenters such as Berkovich and Vickers indenters.
In order to check the feasibility and the sensitivity of the optimization approach, it is appropriate to consider a numerical experimental load-unloading curve from a conical indenter to determine the time-dependent material properties. Table 5 shows the optimised results for axisymmetric conical indenters for the P91 steel material.
The differences between the target and optimised parameters with a conical indenter are within 10%.

Material parameters
Spherical indenter test 1 in Table 4 Conical indenter test 1 in Table 5 2

Young's modulus
Yield stress

Conclusions
In this study, a combined FE analysis based on a two-layer viscoplasticity model, and optimization approach is presented to determine six time-dependent material properties ( , , , , ) of unknown materials from a given loading-unloading indentation curve. Two different materials are investigated, XN40F at 900°C and P91 steel at 600 o C. The optimization algorithm automatically provides input data for the material section in the ABAQUS FE input file and automatically runs FE simulations until the optimised loading-unloading curve reaches the given simulated target loadingunloading curve.
Previous studies have shown that the determination of material properties from time-dependent material behaviour based on conventional indentation test methods does not provide an accurate estimation of the material properties of the indented specimen. The proposed approach can be used to investigate the visco-elastic-plastic material behaviour based on the two-layer viscoplasticity model and determine the time-dependent material properties from the target simulated loading-unloading curve, to within 1-10% error, using results from a spherical indenter, despite using various initial guess values and using different materials. Moreover, there are good agreements between the target and optimised values based on using a conical indenter, although the convergence rate and accuracy depend on the initial input values. Additionally, the optimisation results with six parameters are less accurate compared with that with less parameter due to a non-linear material behaviour. Further research may be targeted at using experimental target loading-unloading indentation curves for a wider range of materials.