Determination of plastic properties using instrumented indentation test with hybrid particle swarm optimization

instrumented indentation test is a promising non-destructive method to determine mechanical properties. This paper proposes a new approach to determine the plastic properties of bulk metal materials (including yield stress, strain-hardening exponent (n) and strain-hardening rate (K)), which couples an experimental load-displacement curve with finite element method. The load–displacement curve was obtained from continuous instrumented indentation test. Then a hybrid particle swarm optimization was employed to minimize the deviation between experimental and simulated load-displacement curves. As a combination of particle swarm optimization and simulated annealing, the simulated annealing particle swarm optimization is an economical and effective algorithm to identify plastic parameters. It was observed that the maximum error of strain-hardening rate extracted from the macro indentation test was 8.2 percent contrast to that determined by the conventional tensile test, and the maximum error of strain-hardening exponent was 4.7% respectively.


Introduction
Indentation tests have been employed to estimate material characteristics extensively. The indentation test was originally used to measure material hardness(H) by dividing maximum indentation load to the contact area [1]. Then the load-displacement curve obtained from indentation test was used to determine Young's modulus(E), based on Hertzian contact theory [2] and Sneddon's axisymmetric elastic contact problem. With the development of highresolution, depth-sensing instruments, Doerner and Nix [3] determined both hardness and Young's modulus based on the indentation load-depth curve(P-h curve) of thin films. Further improvements were made by Oliver and Pharr [4] , They induced Young's modulus and hardness from the indentation load-displacement and demonstrated the possibility of using data directly from the indentation test to determine the area function and mechanical properties without images of indentation. Thereafter instrumented indentation was well developed to obtain material properties such as yield stress, ultimate tensile strength, flow characteristics [5][6][7][8][9][10][11][12], fracture toughness [13][14][15][16][17][18][19][20][21], creep property [22,23] and fatigue property [24,25] because of its convenience, non-destruction and economy.
After several decades' research, varieties of data processing with P-H curve led to varieties theories. Most studies focused on the relationship between hardness (H) and mechanical properties of materials. It has been proved that the P-h curve is sensitive to the material properties [26]. Generally, P-h curve is a function of uniaxial stressstrain, although an accurate prediction of its relationship is sophisticated due to the three-dimensional stress-strain state. Moreover, high nonlinearity due to large deformation beneath the indenter made the characterization process very difficult.
Based on experimental results of spherical indentation test, Meyer found an empirical relationship between force (P) and contact diameter (D) as equation (1): Where α β are material constant. In 1951, Tabor [1] found that the material constant in equation (1) had a relationship with work-hardening exponent (n) according to a power-law stress-strain behavior and it can be expressed as equation (2) as below: This equation is valid especially when applied load is large enough to form a plastic zone near the indenter.
Doerner and Nix proved that the unloading curve can be described with a linear equation as below: Where hc is the true contact indentation depth which takes pile-up and sink-in into account. The top one third of the unloading data should be used to obtain best-fit results.
Oliver and Pharr analyzed six materials to establish the combination of hardness and Young's modulus using Berkovich indenter. It has been shown that the unloading curve of indentation test can be nonlinear and it can be fitted in the form of exponent as below: Where h is the maximum indentation depth during one load-unloading cycle, hf is impressed depth after unloading. A and m are fitting coefficients responding to given materials.
In 2006, Luo J and Lin J [27] extracted mechanical properties from P-h curves using an optimization approach. They minimized the errors between the simulated curves and experimental curves with an objective function using Bates and Watts' optimization method. This approach could deduce material properties (E, σy, n) effectively, but it had to preset a lot of initial values. Thus it was not intelligent and efficient enough. Zhang Tairui [28] proposed a method to determine the optimal constitutive model through spherical indentation test. Dimension analysis and finite element method were employed to determine the constitutive model of materials (linear hardening or power-law) and its parameters. However, too many fitting coefficient are utilized and no specific physical meaning was given.
To deduce material properties from indentation P-h curves, presentative strain was also employed. Ahn and Kwon presented a representative stress and strain method to derive plastic stress-strain relationship from the P-h curves of ball indentation [29]. They made a new definition of strain and investigated the ratio of mean contact pressure to representative stress in indentation. However, it is hard to determine the contact depth because of pileup and sink-in processing around the indenter. With the help of finite element method (FEM), Jian Lu defined a new presentative strain which is a function of the representative stresses and reduced Young's modulus [5]. Combined with dimension analysis and inverse approach, the plastic properties of given material had been derived. The presentative strain had been defined in various ways, but this approach was based on the numerical estimations of empirical parameters.
It has been shown that the pile-up and sink-in had an obvious influence on the precision of determined material characteristics. To avoid the effect of pile-up and sink-in on the contact area edges, an intelligent optimization algorithm named a simulated annealing particle swarm optimization algorithm and finite element simulation were employed to derive mechanical properties of materials from P-h curves directly. This novel approach can obtain material characteristics regardless of any inspection of indentation morphology or various definition of representative strain. The uniqueness of identified parameters was discussed and the stress-strain curves were also verified.

Optimization Algorithm
PSO algorithm, a kind of intelligent algorithm created by Eberhart and Kennedy [30,31], is originally inspired by the social behavior of birds swarm prey. PSO algorithm is widely applied to solve global optimization problems. For this algorithm, particles which have no mass and volumes are regarded as birds. The position of particle i, denoted as x i ( ) = [x ,1 ( ), x ,2 ( ), … x , ( )] at the t-th iteration, is represented as potential solution to optimization problem in a D-dimensional space. The flight direction and speed of particle i are determined by the velocity, denoted as ( ) = [ ,1 ( ), ,2 ( ), … , ( )] . During searching for the optimum solution, particles move around randomly in the multidimensional solution space and record their own previous best position simultaneously. After each iteration, particle i changes its velocity and position according to Pi(t) and Pg, shown in equation (5) and (6). Pi(t) was the best position of particle i which fits the objective function F best after the t-th iteration. The particle with the smallest objective function value which is the best position of particles swarm is denoted by Pg. The set of each particle's best position is marked as P, and = [ 1 ( ), 2 ( ) … ( )] . Pg is an element of set P. Particles keep flying until the iteration criterion is satisfied. The iteration process is introduced in detail as below.
Step 1, n particles are randomly generated in the solution space. xi and vi represent the position and velocity of particle i (i=1,2…n). Then its fitness for the objective function F of optimization problem was evaluated and = [ 1 ( ), 2 ( ) … ( )] and Pg are obtained. Steps 2, initial velocity and position of particles get changed according to the following formulas: , ( + 1) = , ( ) + , ( + 1)，j = 1,2, … , D Where w is inertial weight that controls the impact of previous velocity of particle on current one; positive constant c1 and c2 are learning factor; r1 and r2 are random numbers within the range [0, 1].
Step 3, fitness of each particle in new position is evaluated. Then P and Pg are updated. If the fitness meets the requirements or iteration times exceed the permitted times, Pg will be recognized as the final solution, otherwise turn to step 2.
Although PSO algorithm has a strong ability to find global optimum solution, it may lead to premature convergence because Pg is the exclusive parameter used to update particles' velocity. Therefore the final solution may be affected by the initial position of particles [32,33]. Simulated annealing (SA) algorithm, introduced by Kirkpatrick S [34], has a strong ability to search the local optimum solution while it is weak at finding global optimum solution. The core of SA algorithm is the introduction of annealing temperature T. In the process of annealing, SA algorithm not only accepts the best solution, but also accepts the poor solutions in a probability [35]. And the probability of accepting poor solutions drops with the decrease of annealing temperature T.
Considering the advantage of SA algorithm, the annealing temperature T is introduced into PSO algorithm. So, a new algorithm called simulated annealing particle swarm optimization (SAPSO) algorithm is created by combining PSO algorithm with SA algorithm. The updating of particles' velocity is modified in SAPSO algorithm and described as below.
Where _ is an element of set P. Compared with the PSO algorithm, Pg is replaced by _ in SAPSO algorithm, which means not only the best solution but also a poor solution can be applied to update the new positions of particles. The determination of _ is based on the following inequality.
exp(ΔE /T) ≥ Random(0,1) Where ΔE = ( ) − ( ) is non-positive; T denotes the annealing temperature that is positive; Random (0,1) is a random number between 0 to 1. If the value of exp(ΔE /T) is larger than Random(0,1), the best position of particle i is accepted as the current _ , as shown in Fig.1. As mentioned above, SA algorithm can accept poor solutions probability including the best solution. The SAPSO algorithm inherit this typical feature of SA algorithm when the annealing temperature is introduced. Initially, at large values of T, poor solutions will be accepted high probability; as T decrease, acceptance of poor solutions will become lower; finally, when the value of T approaches 0, poor solutions will not be accepted any more. This feature means that SAPSO algorithm, compared with PSO algorithm, can help escape from local minima while it still retains an excellent ability in searching global optimum solution. Note that the probability of accepting poor solutions is implemented by comparing the value of exp(ΔE /T) with a random number generated from a uniform distribution on the interval (0, 1).
The calculation process of SAPSO algorithm was shown in Figure.1. Similar to the PSO algorithm, the SAPSO algorithm begins with initializing a group of random particles. The initial annealing temperature T 0 is determined based on the objective function F and the first best position of particles swarm Pg. The annealing time is denoted by t. New particles are obtained after transforming each particle's velocity and position according to the equation (7) and (8). If the termination condition is satisfied, Pg will be accepted as the final solution. Otherwise, the iteration will continue.

Identification Procedure
In this work, SAPSO was employed to identify the material properties directly form P-h curves of indentation test. SAPSO is an optimization algorithm integrating particle swarm optimization and simulated annealing. It has a strong ability for searching the global optimal solution of multiple solution problems. Therefore, it can be applied to find the global optimal solution of the objective function efficiently. The procedure of deriving the plastic parameters including yield stress, strain strength coefficient (K) and strain hardening exponent (n) are described in Figure 2. First, Young's modulus can be deduced from the conventional tensile test of a smooth round bar. It also can be deduced from the first unloading part of the P-h curve with Oliver and Pharr's method [26] if the material is difficult to machine into standard tensile specimen. Second, the indentation test with the same material is conducted and Ph curves of indentation are obtained. Numerical simulation of indentation test is performed by a commercial code ABAQUS (version 6.13) with Young's modulus and strain strength coefficient (K) and strain hardening exponent (n) deduced from standard tensile test to ensure that the simulation is valid. Third, preset a set of material parameters (k, n) and apply them to the simulation of indentation test. Then the simulated P-h curve is compared with the experimental curve. The error between the simulated and experimental curves is characterized by an objective function as shown in equation (9) as below.
Where F(x) is the objective function related to E, K, n. and represent the simulated force and the experimental force at the same indentation depth. N is the total number of simulated force points in the loading process. It relies on output request intervals in ABAQUS. Change the value of K and n according to SAPSO and minimize the error between the simulated and experimental P-h curves until it is equal to or less than a preset value.

Experiment and Simulation
A continuous indentation test was carried out on the pressure vessel steel 16MnR and ASTM A193 B16 steel using a DDL 20 machine whose load and displacement resolution is 0.1N and 1um. The chemical composition of 16MnR and ASTM A193 B16 steel is shown in Table 1.All the specimens were machined into a size of Φ30×15mm and were mechanically ground and polished using emery paper (2000 grit) and 1 mm alumina powder respectively. The ball indenter with 0.5 mm diameter was made of tungsten carbide. Instrumented indentation tests were performed according to ISO/TR29381 standard using Instron 5965 system of 5 kN capacity. Specimens were restrained on the test bed by a clamp. The maximum depth was 0.1 mm, and 10 unloading procedure were conducted with a fixed loading and unloading rate (0.1mm/min) by displacement control. Three indentations per specimen were used to examine repeatability and the load-displacement curves were shown in Figure 3. To ensure the accuracy and decrease time consuming, only the middle P-h curves were utilized to compare with the simulated ones. To verify the tensile properties obtained through the instrumented indentation test, conventional tensile tests have been conducted at room temperature using Instron 8800 system of 100 kN capacity. The tests were conducted at crosshead speed of 1 mm/min for specimens with gauge length of 25 mm and diameter of 6 mm according to the ASTM E-8M standard. The commercial software ABAQUS was used to simulate the process of the indentation test. Figure. 4 shows the finite element model of the indentation test. The indenter was modeled as an analytic rigid and the specimens were modeled as deformable body respectively. 4-node axisymmetric linear reduced integral was utilized to mesh the specimen. The minimal mesh size in the contact region where it was under the indenter was 0.025mm. Away from the symmetry axis, the size became larger in order to save time. Because the plastic deformation didn't generate through whole specimen, it's not necessary to establish FE model according to practice size. 10 times of diameter in length and 5 times of diameter in height is sufficient. In this finite element model, the specimen's displacement was constrained in vertical direction (Y-axis). The displacement constraint at an increment of 0.01mm each step was applied to the reference point of the indenter. Surface to surface contact with finite sliding was established. The friction coefficients between the indenter and specimen was set to 0.2. Poisson's ratio was fixed at 0.3 which is proper for metals.

Results and discussion
When the objective function was less than 0.01 in consideration of precision and time consuming. These simulation curves fit well with experiment curves, as shown in Figure 5. At the same time, a set of material properties (K, n) was determined which was close to the value derived from conventional tensile test, shown in Table 2. The maximal error between them was 8.2% which was reliable in engineering. The error introduced in this work mainly by three ways. Firstly, the processing of indentation test can't be accurately descripted by the elastic-plastic analysis because damage caused by plastic deformation was not taken into account. Secondly, power hardening model which is a simplified assumption of material constitutive relationship might introduce error especially where strain is relatively high. Thirdly, wear of experiment apparatus like indenter and precision of sensor also had effect on the accuracy of this approach.  It has been proved that the material properties can't be derived only from the P-h curves obtained from a single indentation test. Because many set of mechanical properties will result in the same P-h curve. To exclude the interference of the improper set of mechanical characteristics, continuous indentation test was utilized. Material strain hardening phenomenon which is affected by strain strength coefficient (K) and strain hardening exponent (n) can be reflect through more loading and unloading cycles. Large number of loading and unloading cycles will lead to precise identification of material properties while computation time will increase obviously. Since ten loading cycles can result in a unique and satisfying set of material parameters, continuous indentation test which has more than ten loading cycles were not carried out and identified.
In this work, power-hardening constitutive relationship was preset based on the premise that the material were coincident with the power-hardening model. However when the material's constitutive relationship is unknown, application of this approach may be restricted. More attention can be paid to the identification of unknown material constitutive relationships in further study.

Conclusions
SAPSO was an effective method to extracted plastic properties from load-displacement curve. The uniqueness and precision can be ensured. It was observed that the maximum error of strain-hardening rate extracted from the macro indentation test was 8.2 percent contrast to that determined by the conventional tensile test, and the maximum error of strain-hardening exponent was 4.7% respectively.