Extracting the isotropic uniaxial stress-strain relationship of hyperelastic soft materials based on new nonlinear indentation strain and stress measure

Instrumented indentation technique has been increasingly utilized to measure the mechanical properties of soft polymers and biological tissues. However, the indentation behaviors of these materials has not been well understood, especially the parameter identification of their hyperelastic material properties. In this paper, we developed a spherical indentation data analysis method to directly extract the isotropic uniaxial stress-strain relationship of hyperelastic soft materials from the measured spherical indentation load-displacement curves. The proposed method mainly included new measure of indentation stress and strain, which was built based on the Hertz load-displacement relationship and further revised by considering the non-Hertzian effects of neo-Hookean hyperelastic contact problems. Numerical and actual indentation experiments showed the proposed definition of indentation strain can properly evaluate the amount of nonlinear strain for neo-Hookean, Yeoh and Arruda-Boyce hyperelastic materials. Meanwhile, the proposed spherical indentation data analysis method was applicable only in certain deformation range for Yeoh and Arruda-Boyce hyperelastic materials, because their nonlinear material parameters might cause very complicated contact pressure distributions. Building a universal data processing technique for characterizing the hyperelastic mechanical properties of soft materials through indentation experiments still needed further investigations.Instrumented indentation technique has been increasingly utilized to measure the mechanical properties of soft polymers and biological tissues. However, the indentation behaviors of these materials has not been well understood, especially the parameter identification of their hyperelastic material properties. In this paper, we developed a spherical indentation data analysis method to directly extract the isotropic uniaxial stress-strain relationship of hyperelastic soft materials from the measured spherical indentation load-displacement curves. The proposed method mainly included new measure of indentation stress and strain, which was built based on the Hertz load-displacement relationship and further revised by considering the non-Hertzian effects of neo-Hookean hyperelastic contact problems. Numerical and actual indentation experiments showed the proposed definition of indentation strain can properly evaluate the amount of nonlinear strain for neo-Hookean, Yeoh and Arruda-Boyce hyperelastic materials. M...


I. INTRODUCTION
Instrumented indentation technique 1 has been an important experimental method to measure the mechanical properties of various soft materials including soft elastomers, polymeric gels and biological tissues. 2,3 These materials are "soft" and commonly exhibit nonlinear elastic stress-strain responses under external loadings, which can be well described by the hyperelastic constitutive relationship under finite deformation. [4][5][6][7] The hyperelastic indentation problems involved the coupling interaction of geometric, material and boundary nonlinearity, which make it difficult to be solved analytically. 8 Due to the lack of a proper nonlinear contact model of the indentation test, accurate identification of the hyperelastic material properties of soft materials through indentation experiments still faces great difficulties.
A large amount of work has been carried out on the identification of nonlinear material properties or the constitutive parameters of hyperelastic soft materials through indentation experiments using inverse methods. For example, Hyun et al. 9 developed a spherical indentation approach to evaluate the material property of rubber materials described by the Yeoh model. Giannakopoulos et al. 10 derived a theoretical solution for spherical indentation of incompressible Mooney-Rivlin materials, which was reliable under indentation with the ratio of the indentation depth to the indenter radius smaller than 10%. Based on the Sneddon solution, 11 Liu et al. 12 proposed a nonlinear elastic contact model by utilizing a second-order approximation of spherical function to predict the large indentation behaviors of rubberlike materials. Zhang et al. 13 developed a simple data analysis method to extract the initial shear modulus of hyperelastic soft materials. Zisis et al. 14 present a new spherical indentation method of determining the constitutive parameters of incompressible, initially isotropic, hyperelastic substrates under equal-biaxial prestretching. By using the finite element computation combined with a numerical optimization subroutine, Chen et al. 15 developed an inverse analysis method to characterize the finite viscoelasticity and to simultaneously identify the influence of adhesion in nanoindentation experiments of soft polymers. Despite these advancements, these analytical or numerical models are all built based on a given hyperelastic constitutive model, the inverse analysis process needs to firstly choose a suitable hyperelastic constitutive model. Besides, the hyperelastic parameters extracted from the indentation load-displacement curve by these data analysis methods are poorly defined, the existence, uniqueness and stability of solutions obtained through inverse analysis should be further clarified. 16,17 Compared with the inverse analysis method, the indentation stress-strain(ISS) curves offered the potential for characterizing and understanding the local stress-strain behavior of a material at the microscale. As we all know that the stress-strain curves under uniaxial loading can be easily extracted due to the uniform stress state of materials. While obtaining the complete stress-strain responses of materials under contact loading is more complex because of the significantly inhomogeneous stress state. Tabor 18,19 was the earliest to investigate the analogy between uniaxial compression and spherical indentation, and introduced the concept of indentation stress-strain curves. The indentation stress was usually defined as the mean pressure on the projected contact area. Then, the indentation strain can be obtained by transforming the Hertzian equation into a linear relationship. For spherical indentation on a flat sample surface, the indentation stress and strain can be expressed as where E * = E/(1 − υ 2 ) denoted the reduced elastic modulus, E and υ were the Young's modulus and Poisson's ratio of the indented sample, respectively. In Tabor's approach, the indentation strain was further invoked to be 0.2a/R, where the strain prefactor of 0.2 was empirically determined for spherical indentation on elastic-plastic metals. Based on the Tabor's definitions of indentation stress and strain, Lin et al. 20,21 proposed several approximate nonlinear contact models for several common hyperelastic strain energy functions, and validated these models under indentation deformation with the indentation depth to indenter radius smaller than 0.2. The other representative definitions of indentation stress and strain was proposed by Kalidindi et al., 22 Donohue et al. 23 and Pathak et al. 24 by recasting the Hertzian equation into a different linear relationship between indentation stress and strain: This definition of elastic indentation strain can be visualized by idealizing the primary zone of indentation deformation as being equivalent (in an average sense) to compressing by h a cylindrical region of radius a and height 3πa/4. Besides, Patel et al. 25 presented simple revised protocols for recovering uniaxial mechanical responses of linear/power hardening isotropic elastic-plastic materials by involving simple scaling relationships. The spherical indentation stress and strain measure above mainly originated from the Hertzian load-displacement relationship was only applicable to linear elastic indentation under small penetration or contact radius. In this paper, we proposed a data analysis method for directly extracting the nonlinear stressstrain relationship of hyperelastic soft materials from the spherical indentation load-depth curves. The proposed method mainly included novel indentation stress and strain measure, which was developed by considering the non-Hertzian effects associated with spherical indentation of hyperelastic soft materials. Moreover, the correlations between the extracted spherical indentation stress-strain curves(ISS) and the uniaxial stress-strain curves (USS) of several kinds of hyperelastic soft materials were further analyzed and discussed.

II. INDENTATION STRESS AND STRAIN MEASURE
The spherical indentation stress and strain measures in the literature mainly originated from the Hertzian load-displacement relationship. However, the indentation behaviors of hyperelastic soft materials involved coupled constitutive, geometric and contact nonlinearities, which made the indentation stain/stress fields obviously deviating from the prediction of the Hertz contact model. Fig. 1(a) showed a schematic diagram of a deformed hyperelastic half-space indented by a rigid spherical indenter. In this diagram, P, h and a were the indentation load, indentation displacement and contact radius, respectively. And R was the radius of spherical indenter, z was the axis of symmetry and x-y plane was the common tangent plane between the spherical indenter and the indented hyperelastic half-space. So, the spherical indentation stress and strain measures originated from the Hertz model were not applicable to the data analysis of spherical indentation on hyperelastic materials. Here, novel measure of indentation stress and strain were adopted and expressed as follows: where σ was the indentation(effective) stress, and ε was the indentation(effective) strain, P was the indentation load, h was the indentation displacement, a was the contact radius. And was the reduced elastic modulus, E 0 and υ were the initial elastic modulus and Poisson's ratio of the indented hyperelastic half-space, respectively. The indentation stress and strain measures given in Eq. (3) was constructed based on the classic Hertz contact model, which was not applicable to the data analysis of spherical indentation of hyperelastic materials due to the non-Hertzian effects. Although, Wu et al. 26 has showed that the Hertzian load-displacement relation still held for spherical indentation on soft materials undergoing indentation deformation with h/R as large as 0.66, and pointed out that this agreement aroused from a near cancellation of two non-Hertzian effects: the spherical (as opposed to paraboloidal) shape of the indenter, and the large deformation behavior of linear elastic system. Considering that the two non-Hertzian effects would cause the indentation stain/stress fields of soft materials obviously deviating from the prediction of the Hertz contact model, two dimensionless revised functions ψ(h/R) and χ(a/R) were introduced to reliably extract the spherical ISS curve of hyperelastic soft materials. The revised functions should be related to the non-Hertzian effects including the spherical shape of indenter, the large deformation behaviors and the nonlinear constitutive relationships. The specific illustration and determination of the two revised functions would be shown in the following section.

A. Finite element modeling
In this section, an axisymmetric finite element model was developed to simulate the spherical indentation experiments of typical hyperelastic materials by using the finite element package ABAQUS. 27 Rigid Spherical indenters with different radius were applied, their elastic modulus for all the indenters was set to 1142GPa, and Poisson ratio 0.07. The substrate could be regarded as a half-space since its height and width were both much larger than the contact radius and penetration. The bilinear axisymmetric reduced integration element was used to discrete the substrate and a refined mesh was employed in the contact zone beneath the indenter tip to ensure the accuracy of the numerical results, as shown in Fig. 1(b). The contact between the spherical indenter and the substrate was defined as frictionless. Radial displacements and rotations of elements along the axial line of symmetry were fixed, which enforced an axisymmetric boundary condition. The vertical displacement of elements at the bottom of the substrate was fixed, and the top surface outside the contact was traction-free. The rigid indenter could move only in the vertical direction. In all the simulations, the maximum indentation depth, h max , was set as h max =0.8R and the geometric nonlinearity, such as large rotation and large strain, was included.
Three kinds of widely used hyperelastic constitutive models, including neo-Hookean, Yeoh and Arruda-Boyce model (described in Appendix), were used to describe the nonlinear mechanical behaviors of the incompressible hyperelastic substrate. The neo-Hookean model was the simplest hyperelastic constitutive model with one material parameter C 10 . The Yeoh model was applicable for large deformation with simple experimental data with three material parameters C 10 , C 20 and C 30 . The initial shear modulus µ 0 is given by µ 0 = 2C 10 for both Yeoh model and neo-Hookean model. The Arruda-Boyce model, could accurately predict the nonlinear mechanical behavior of most engineering elastic materials with only two material parameters, i.e., the shear modulus µ and locking stretch λ m . All the three models were included in the material library of ABAQUS software, could be directly used. Numerical simulation of displacement controlled indentation on these three kinds of hyperelastic soft solids were performed. Also, numerical simulation of spherical indentation on linear elastic materials undergoing very large deformation with the indentation depth 0.8 times the indenter radius, was also performed as a comparison.

B. Contact radius
Accurate estimation of the contact radius is central to the data analysis of indentation experiments. Here, to better describe the shape of spherical indenter, the second-order approximation of spherical function was utilized and expressed as follows: Ignoring the adhesion interaction between the spherical indenter and the surface of the sample, the relationship between indentation depth and contact radius developed by Sneddon et al. 11 was further applied and expressed as h Due to the contact radius cannot be directly measured in most commercial instrumented indentation systems, this relationship was convenient to be used to calculate the indentation stress for spherical indentation on soft materials undergoing very large deformation. Meanwhile, when this relationship was used to process the indentation load-displacement data, the revised function for indentation strain should be intensively defined as the strain field obviously deviated from the prediction of the Hertz model which applied a parabolic to describe the shape of spherical indenter.

C. Revised function for indentation strain
The developed finite element model was used to simulate the large spherical indentation experiments of materials exhibiting both linear elastic and hyperelastic constitutive behaviors, respectively. Our previous research 28 has showed that the contact surface profiles of indented elastic and hyperelastic half-space coincided well with each other, but deviated from the prediction of Hertz model. Moreover, Fig. 2

D. Revised function for indentation stress
After the determination of indentation strain, proper evaluation of the indentation stress must be addressed in order to extract the indentation stress-strain curve of hyperelastic materials. The calculated indentation load-displacement curves of neo-Hookean materials were obtained and converted into the spherical ISS curves by applying the measure of indentation stress with no revision and the nonlinear indentation strain. Fig. 4 compared the extracted spherical ISS curve with the USS curve of the indented neo-Hookean materials. It could be seen that the two curves exhibited obvious disagreements, which was because that the effects of material nonlinearity were not considered. Considering the deviation of the indentation load-displacement curves between linear elastic and neo-Hookean materials was only originated from the material nonlinearity of neo-Hookean hyperelastic constitutive behaviors, the revised function for indentation stress could be determined by matching the extracted spherical ISS curve with the corresponding USS curve of neo-Hookean materials, which was expressed as follows: Fig. 4 also compared the final revised spherical ISS curve with the USS curve of neo-Hookean materials. Given that the dimensionless contact pressure distributions at different depth were independent on the initial shear modulus for spherical indentation on neo-Hookean materials, this revised function for indentation stress was general for this kind of hyperelastic soft materials undergoing very large indentation deformation.

IV. APPLICATION FOR OTHER HYPERELASTIC MATERIALS
The proposed definitions of indentation stress and strain considered the non-Hertzian effects associated with spherical indentation of soft materials including the spherical shape of the indenter, the large deformation behavior and the material nonlinearity of neo-Hookean hyperelastic model. Applying the proposed definitions of indentation stress and strain, the isotropic uniaxial stress-strain relationship could be extracted by processing the measured indentation load-displacement curves of neo-Hookean hyperelastic materials. Meanwhile, whether the proposed definitions of indentation stress and strain as well as the data analysis method were applicable likewise to other kind of hyperelastic materials, should be further clarified.
Here, the literature data 16 obtained from the uniaxial tensile tests and spherical indentation tests of polydimethylsiloxane (PDMS) was reanalyzed. The uniaxial tensile stress-strain curve was firstly fit using the neo-Hookean, Yeoh and Arruda-Boyce hyperelastic models to extract its hyperelastic material parameters. It was found that the neo-Hookean model could only fit it well under small deformation, other two models gave very precise predictions even under very large deformation. Then, the proposed data analysis method was utilized to extract the ISS curve by processing the experimentally measured indentation load-displacement curve. Fig. 5 compared the extracted spherical ISS curve with the USS curves predicted by the neo-Hookean, Yeoh and Arruda-Boyce hyperelastic models. It was shown that these curves agreed well with each other under indentation deformation with indentation strain smaller than 0.2, but they showed obvious disagreements as the deformation further increased. This disagreement might originate from that the nonlinear mechanical behaviors characterized by the three hyperelastic models have different influence on their spherical indentation responses. The hyperelastic material parameter was also obtained by fitting the extracted indentation stress-strain curve with indentation strain smaller than 0.2 using these three hyperelastic constitutive models. The initial shear modulus was about 0.34MPa for neo-Hookean, Yeoh and Arruda-Boyce models, which was very close to the uniaxial results.
To further illustrate the application of the proposed definitions of indentation stress and strain for Yeoh and Arruda-Boyce hyperelastic materials undergoing large indentation deformation, the FE model developed in Section II was utilized to simulate the spherical indentation of Yeoh and Arruda-Boyce hyperelastic soft materials. The key contact characteristics including the load-displacement curves, contact radius and contact pressure distributions were obtained and compared with those of neo-Hookean materials.  materials with different material parameter λ m under several different indentation deformation with those of neo-Hookean materials, respectively. It was shown that the contact surface profiles of Arruda-Boyce materials coincided well with that of neo-Hookean materials even under the deformation range with the h/R equal to 0.8. And the λ m has no obvious influence on the contact surface profile under the same h/R, which reflected that the indented materials undergo the same contact deformation. Similarly, for the Yeoh hyperelastic materials, both the material parameter C 20 and C 30 has no obvious influence on the contact surface profile under the same h/R, as shown in Fig. 7. So, the proposed definition of indentation strain could properly evaluate the amount of nonlinear strain for both Arruda-Boyce and Yeoh hyperelastic materials. Fig. 8 and Fig. 9 further compared the contact pressure distributions of Arruda-Boyce and Yeoh hyperelastic materials with those of neo-Hookean materials, respectively. It was shown that the contact pressure distributions of Arruda-Boyce hyperelastic materials showed good agreements with that of neo-Hookean materials under indentation deformation h/R smaller than 0.3. Under larger indentation deformation range, they deviated obviously from each other and the degree of deviation showed a gradually decrease tendency under the same indentation deformation h/R as the material parameter λ m increased. This might originated from that as the material parameter λ m was increased, the indentation responses exhibited a gradually softening tendency. For the Yeoh hyperelastic materials, the material parameter C 20 and C 30 also has no obvious influence on the contact pressure distributions under indentation deformation with h/R smaller than 0.3, but the nonlinear mechanical behaviors characterized by the Yeoh hyperelastic models would cause the contact pressure distributions significantly deviating from that of neo-Hookean materials under larger deformation range. Fig. 9(a) showed the contact pressures exhibited a gradually decrease tendency under the same indentation deformation h/R as the material parameter C 20 decreased; Fig. 9(b) showed that as the material parameter C 30 increased, the contact pressures showed a gradually increase tendency under the same indentation deformation h/R. This was mainly because the material parameters C 10 , C 20 and C 30 have different influences on the indentation behaviors. In the initial stage of deformation, C 10 dominated the material mechanical responses. As the deformation was increased, C 20 became a second influential variable and brought softened effect, followed by C 30 with hardening effect. Although the proposed definition of indentation strain can properly evaluate the amount of nonlinear strain for Yeoh and Arruda-Boyce hyperelastic materials, the proposed spherical indentation data analysis method was applicable only in certain deformation range for Yeoh and Arruda-Boyce hyperelastic materials, because their nonlinear material parameters might cause the contact pressure distributions significantly different from that of neo-Hookean materials under larger indentation deformation.

V. CONCLUSIONS
In the present paper, a novel spherical indentation data analysis method was proposed based on the concept of indentation stress-strain(ISS) curves. The definitions of indentation stress and strain was constructed based on the Hertz load-displacement relationship and further revised by considering the non-Hertzian effects of neo-Hookean hyperelastic contact problems. Applying the proposed method, the isotropic uniaxial stress-strain relationship of hyperelastic soft materials could be directly extracted from the spherical indentation loading curves. The validation and application was analyzed and discussed through numerical simulations and actual indentation experimental data. It was found that the proposed definition of indentation strain can properly evaluate the amount of nonlinear strain for neo-Hookean, Yeoh and Arruda-Boyce hyperelastic materials. While, due to the contact pressure distribution was dependent on the nonlinear material parameters, it was difficult to develop a general indentation stress measure for Yeoh and Arruda-Boyce hyperelastic materials. The proposed spherical indentation data analysis method was applicable only in certain deformation range for Yeoh and Arruda-Boyce hyperelastic materials. In any event, the overall approach presented in this work represented a new path for characterizing the nonlinear mechanical behaviors of soft materials, where no prior knowledge about the constitutive responses of indented materials was required.

ACKNOWLEDGMENTS
This work was supported by the National Natural Science Foundation of China (No. 11802166).

APPENDIX
In this study, three incompressible hyperelastic constitutive models, including neo-Hookean, Yeoh, and Arruda-Boyce models, which have been widely applied to model the nonlinear deformation behavior of elastomers, polymeric gels and biological soft tissues. These models are briefly given as follows:

Neo-Hookean model
The neo-Hookean model is usually considered to be the reduced polynomial form of Mooney-Rivlin model. In the neo-Hookean model, the strain energy density function can be written as where C 10 and D are material parameters, and the initial shear modulus was given by µ 0 = 2C 10 , initial bulk modulus K 0 was related to D by K 0 = D 2 , where D = 0 for incompressible materials. So, the Poisson's ratio can be expressed by υ = 3K 0 −2µ 0 6K 0 +2µ 0 .ψ is the strain energy density per unit volume in the undeformed configuration, J is the volume ratio. And the first deviatoric strain invariantsĪ 1 can be expressed asĪ 1 =λ 2 1 +λ 2 2 +λ 2 3 , whereλ i (i=1,2,3) is the deviatoric principal stretch and related to the principal stretch λ i byλ i = J − 1 3 λ i . The incompressible neo-Hookean model was the simplest hyperelastic constitutive model with one material parameter C 10 and was effective only in relatively small strain range.

Yeoh model
Yeoh hyperelastic model which is a phenomenological model and applicable for larger deformation of nearly incompressible, nonlinear elastic materials. The strain energy density function is represented by ψ = c 10 (I 1 − 3) + c 20 (I 1 − 3) 2 + c 30 (I 1 − 3) 3 (A2) where C 10 , C 20 and C 30 are material parameters, and the initial shear modulus µ 0 is related to C 10 by µ 0 = 2C 10 . Although the Yeoh model is also a reduced polynomial form of hyperelastic constitutive model, it can predict well the behavior in different deformation modes from data gained in just one deformation mode. It is noteworthy that when n = 1, the Yeoh model reduces to the neo-Hookean model for incompressible materials.

Arruda-Boyce model
Arruda-Boyce hyperelastic constitutive model, also known as the eight-chain model, is developed based on the statistical mechanics of a material with a cubic representative volume element containing eight chains along the diagonal directions. Applying the first five terms of the inverse Langevin function, the strain energy density function for the incompressible Arruda-Boyce model is given by where β = 1 λ 2 m , α 1 = 1 2 , α 2 = 1 20 , α 3 = 11 1050 , α 4 = 19 7000 , α 5 = 519 673750 . This model could accurately predict the nonlinear mechanical behavior of most engineering elastic materials with only two material parameters, i.e., the shear modulus µ and the limiting network stretch(locking stretch) λ m . In this model, the initial shear modulus µ 0 is related to the shear modulus µ by µ 0 =