Elimination of Pseudo Initial Stress for Proposed Invariants based Hyper elastic Strain Energy

Although fitting of experimental data is important in selection of the best hyper elastic strain energy function, the zero stress condition in the initial (unstressed) configuration also should be taken into account. Therefore, in the current research a modification for an invariants based model was developed to eliminate the pseudo initial stress. To this goal, additional terms were added to hyper elastic strain energy as a function of the right Cauchy-Green deformation tensor. To justify the proposed algorithm, uniaxial loading of liver tissue and unfilled silicone rubber were studied via constitutive modeling in VUMAT user subroutine of ABAQUS software. Excellent agreement between the experimental and numerically predicted stress-stretch curves showed that the modified strain energy function enables to characterize the mechanical behavior of the liver tissue and unfilled silicone rubber. As the deformed shapes of soft tissues are important especially in the surgery operations, the effect of modification on the final deformed geometry was investigated. Results were shown the deformed profile of the modified and unmodified hyper elastic strain energy models were different.


Introduction
In the last few years, many studies have concentrated on the accurate prediction of nonlinear behavior of materials such as soft tissues and rubber like materials.Modeling of the hyper elastic behavior involves selection of an appropriate strain energy function, determination of material constants, and constitutive modeling to investigate the response of model in real loading conditions [1,2].In this context, a large number of hyper elastic strain energy functions have been proposed by researchers [3].In general, the hyper elastic strain energy functions have been based on invariants or principal values of the stretch (or Cauchy-Green) tensor [4].Due to the nonlinear mechanical behavior of soft tissues and rubber like materials, constitutive equations that calibrated with one loading condition can not necessarily evaluate the mechanical response under other loading modes [5][6][7][8].Therefore, in order to predict accurately the mechanical properties of biological tissues, experimental tests in various loading conditions are preferred [9][10][11].
Meunier et al. [12] performed experimental and numerical tests to describe the mechanical behavior of unfilled silicone rubber.Results indicated that the Haines-Wilson, Gent, and Ogden models can present proper behavior of the silicone rubber and the Mooney-Rivlin model prediction is not satisfactory.Jamshidi and Ahmadianused the Mooney-Rivlin, polynomial, Yeoh, and Ogden models to represent the mechanical properties of heart tissues and found the Yeoh material model as the best one [13].Kaster et al. [14] employed the polynomial, Yeoh, Arruda-Boyce and Ogden models for accurate characterization of the mechanical behavior of brain tissue.Although they pointed out the Yeoh model is more efficient because of its accuracy, fewer material parameters, and shorter computational time, but this model depends only on the first invariant of stretch tensor and in pure shear could not predict appropriate behavior [4].Gokgol et al. [15] used the five terms Mooney-Rivlin model to validate the experimental results on the bovine liver tissue for estimation of fracture toughness.They used the inverse solution to estimate the material parameters.the combined mathematical functions for better characterization of the nonlinear behavior of materials.Gao et al. [17] proposed a new constitutive model by adding the logarithmic function to Ogden hyper elastic model.The proposed model represented the liver tissue behavior under tension, compression, and pure shear deformation modes in large strains.In order to describe the significance of accurate zero-strain (or stress) in estimation of nonlinear material properties, Chui et al. [18] used the full cycle of compression and tension on the liver tissue.They proposed a new constitutive equation based on the combination of polynomial and logarithmic functions.In order to simulate the indentation test of porcine liver tissue, Fu and Chui [19] used the compression and tension experimental data for estimation of material parameters and showed the logarithmic-polynomial model is partly better than the five terms Mooney-Rivlin model.Darijani and Naghdabadi [20] proposed hyper elastic strain energy models based on polynomial, power, logarithmic and exponential function of principal stretches to capture the large deformation behavior.Khajehsaeid et al. [11] proposed a strain energy function for large deformation of isotropic incompressible rubber-like materials with combination of exponential and logarithmic functions and found the proposed model presents good agreement with experimental data.Hosseinzadeh et al. [21] applied the Mooney-Rivlin and Ogden models as well as the general exponential-exponential and general exponential-power law models for investigation of the demineralized and deproteinized bovine cortical femur bones.
Their results illustrated that among the employed models the general exponential-exponential and general exponential-power law models have a good agreement with the experimental results and material stability.Elyasi et al. [22] used the general exponential and invariants models to find the best hyper elastic function that fit the experimental data of extensor apparatus.They investigated the material instability and accordingly obtained the stress distribution in flexion movement.Mansouri and Darijani [4] proposed a new framework for strain energy functions of the isotropic hyper elastic materials in terms of the first and second invariants of the right Cauchy-Green deformation tensor.According to their results, the in variants model illustrates precisely the stress-strain behavior of elastomers and soft tissues in special test conditions and its performance for general cases should be investigated in FE-based modeling techniques.
In the current research, the invariants model proposed by Mansouri and Darijani is modified to consider the initial free stress configuration and elimination of spurious mathematical stress.To justify the significance of initial free stress configuration, the deformation behavior of the porcine liver tissue and silicon rubber is simulated.To this end, in the next sections firstly the modification of the hyper elastic model and accordingly the constitutive modeling are presented.Finally, the results of simulations will be presented and discussed.

Hyper Elastic Model and Unstressed Configuration Modification
As it was mentioned, Darijani and Naghdabadi [20] proposed hyper elastic strain energy models based on a combination of polynomial, power law, logarithmic and exponential functions.Hyper elastic strain energy functions not only should be able to fit the experimental data, but also should be as simple as possible and having a small number of parameters.Meanwhile, Mansouri and Darijani [4] generalized the models of Darijani and Naghdabadi [20] and they called the hyper elastic functions as general power and general exponential models.Their hyper elastic functions were based on the first and second invariants of the right Cauchy-Green deformation tensor and they showed the robustness of models in describing the mechanical behavior of isotropic materials.Thereby, in this study the model in terms of the first ( 1I ) and second ( 2 ) invariants of the right Cauchy-Green deformation tensor ( C ) was used.In the following, for simplicity, this model is called the invariants model and it can be expressed as [4]: Where 1 B , 1 m , and 1 n are the unknown material parameters and W represents the hyper elastic strain energy function.The second Piola-Kirchhoff stress components ( S) can be derived from the hyper elastic strain energy function as [21,22]: Derivatives of the invariants of the symmetric tensor C with respect to itself can be mentioned as: Where I is identity tensor.Assuming of Eqs. 1, 2 and 3, the second Piola-Kirchhoff stress can be derived as: Substituting of the initial (undeformed) configuration condition ) in Eq. 4 the following initial values can be obtained: Where initial S represents the initial stress values and as it can be seen it only includes the normal components.It should be noted the initial stress does not any physical interpretation like the physical initial stresses in biomaterials or residual stresses in some elastomers.Therefore, it has the mathematical sources (spurious or pseudo stresses) and it should be eliminated.Indeed, the constant "-1" in Eq. 1 was added to consider the zero value for the hyper elastic strain energy function in the initial configuration as a basic requirement [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20].However, applying of Eq. 2 results to mathematical initial pseudo stress presented in Eq. 5.In order to consider the zero initial stress condition, the constraint term was added to Eq. 1 and the modified strain energy equation was derived as: Where, Modify W is the modified strain energy function that consists of the invariants model terms and the zero initial stress constraint.Moreover, the initial energy yet is zero and as Ogden et al has mentioned it is one of the methods to impose a constraint on the hyper elastic strain energy [23][24][25][26].Regarding Eqs. 2, 3 and 6the modified second Piola-Kirchhoff stress can be stated as: ] Now, by inspecting of Eqs. 6 and 7 it could be ensured the stress and energy are zero in the initial configuration.

Constitutive Modeling
To use a hyper elastic strain energy function for an incompressible material, another constraint due to incompressibility should be imposed on the hyper elastic energy function.In order to apply the incompressibility, Eq. 6 should be modified as [27]: Here, k represents a small number and the incompressibility condition is imposed by the penalty method.Considering of Eqs.
2 and 8, the second Piola-Kirchhoff stress components in the presence of the incompressibility condition could be stated as:

m I n I S A m e I B ne I I C A m B n I C I I C k
Now by inspecting of Eqs. 8 and 9 the initial stress and energy are zero and the incompressibility condition is also imposed.In order to implement the above formulas in VUMAT user subroutine of the ABAQUS software, the second Piola-Kirchhoff stress should be converted to the co rotational Cauchy stress ( σ )by the following relation: Here U represents the right stretch tensor and S the second Piola-Kirchhoff stress tensor.

Experimental Data
In order to evaluate the modified model, two datasets including tension of porcine liver tissue and compression of unfilled silicone rubber have been used from literatures.The nonlinear biomechanical property of the liver tissue was extracted from experimental study of Gao et al. [17].Mansouri and Darijani [4] utilized this test data for proofing of robustness of the invariants model and showed the invariants model fitted both of the pure shear and uniaxial tension loading.The silicone rubber is common material in biomechanical applications because of its unique properties and excellent biocompatibility [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27].The material constants of the porcine liver tissue and unfilled silicone rubber can be seen in Table 1.

Table 1:
Material constants of the porcine liver tissue and unfilled silicone rubber [4].

Simulations
Among several types of virtual surgery simulators, the finite element method (FEM) has been commonly used.This method is useful for solving of problems with complex geometry and loading as well as complex nonlinear materials response.This method can be an appropriate method for analysis of the materials stress and strain [28][29][30] because of realistic materials behavior.
A 3D deformable model was built to simulate the uniaxial tension of porcine liver sample.The length, width, and thickness of the sample were considered as 25, 15, and 17mm respectively similar to the experimental conditions.The 8-node linear brick elements with reduced integration and hourglass control (C3D8R element) were used for simulation.Reduced integrations and enhanced hourglass control were utilized to prevent the volumetric locking and elimination of numerical instability [31].This model consists of 7488 nodes and 6375 elements.The boundary conditions were assigned according to experimental (Figure 1) and geometry symmetry conditions such that the bottom of the sample was fixed in the Y direction, the upper surface was stretched in the Y direction by and the X and Z directions were fixed.The material parameters of Table 1 were assigned to the invariants model as the material properties.The silicon rubber three dimensional (3D) model was created for the compression test simulation.Using symmetry constraints, only quarter of the specimen was considered to increase the numerical stability and reduce the computational time.The height of 20 mm and diameter of 37 mm similar to experimental sample was considered.The compression platen was modelled as a rigid body with diameter of 60 mm.The contact between the platen surface and the rubber specimen were applied.Because in experimental the silicone grease was used as lubricant, the Coulomb friction coefficient between the platen and rubber was considered as 0.005.The simulation was performed by applying a displacement of to the platen in the vertical direction to produce the maximum stretch of 0.48.The applied boundary conditions to the sample were shown in Figure 2. The reaction force and displacement were calculated along the loading direction in the reference point of platen.By dividing of the reaction force to the initial surface of sample, the nominal stress was calculated.The simulations details of two tests were summarized in Table 2 and (Figures 2 ).

Results and Discussion
In order to examine the proposed modification on the hyper elastic model, in this section the results of porcine liver tension and silicon rubber compression are presented and discussed.

Porcine Liver Tension
For simulation of the porcine liver tissue behavior in uniaxial tension, the invariants model was used.As mentioned by previous researchers [5][6][7][8][9][10] a hyper elastic model that enables to simultaneously fit all loading conditions and represents material stability could be considered as a suitable model for simulation of finite deformation in complex geometry and loadings.In other hand, a satisfactory prediction significantly depends on the accuracy of geometrical model and constitutive equation [32].The results of simulations with both of the modified and unmodified invariants model were illustrated in Figure 3.As it can be seen from (Figure 3), there are significant differences between the results of the unmodified model and the experimental data.According to definition of isotropic stress free initial configuration, the reference state for an isotropic model should be undistorted and stress-free (zero-strain state) [33].As it can be seen in Figure 3, by applying the unstressed initial configuration condition on the modified model, the results are very close to the experimental curve.The von Mises stress contours of the porcine liver tension simulation for the modified and unmodified models are provided in Figure 4. To be clear the difference of the deformed geometries resulted from the modified and unmodified models, in Figure 5 the normalized coordinate with respect to the initial coordinate of nodes (i.e. x X , where X is the initial coordinate and x is the current coordinate) were plotted against to the distance along the free edge (from start to end node).As it can be recognized the error due to use of unmodified formulation can be significant and it is known the deformed profile in the biomechanics application has great importance in the interpretation of loading results.In addition, the mesh convergence study was revealed increasing of elements do not impact the results.

Silicon Rubber Compression
For further investigation of the modification of hyper elastic model, the compression of the unfiled silicon rubber was studied.This problem is operational in the simulation of materials behavior especially the soft tissues as in plenty of the surgery operations, the soft tissue is in contact with the surgery devise.Therefore, understanding the soft tissue behavior in this kind of tests can be helpful in the surgical simulations [14][15][16][17][18][19][31][32][33][34].The stress-stretch curves obtained from the modified and unmodified invariants model for compression of unfiled silicon rubber were shown in Figure 6.As it can be seen when the unmodified formulation was used, it was not possible to capture the nonlinear behavior of the silicon rubber.On the other hand, the modified model showed a good agreement with the experimental data.
Figure 7 illustrates the von Mises stress contours and Different stress distribution can be seen for the modified and unmodified formulations.The normalized z coordinate with respect to the initial coordinate ( z Z , where Z is the initial coordinate and z is the current one) were plotted against to the distance along the free edge in Figure 8.As in the porcine liver tension, there is discrepancy between two formulations.

Conclusion
It was observed by previous researches the invariants model is enable to fit the experimental data of tension, compression, and balanced biaxial tests.As this model does not provide the zero initial stress, in this study, the role of unstressed configuration on the finite element results was investigated.To this goal, the initial stress was removed from the hyper elastic strain energy via a function of the third invariant of the right Cauchy-Green deformation tensor.Inspection of the new hyper elastic strain energy revealed the zero energy and zero stress could be captured in the unstressed initial configuration.To justify the effects of the modification, the deformation behavior of porcine liver tissue and unfilled silicone rubber were considered.Results indicated that there is a significant difference between the predicted and experimental stress-stretch curves when the modification is not considered.In addition, it was seen the deformed geometry were different for the modified and unmodified models.

Figure 1 :
Figure 1: Geometry and boundary conditions in uniaxial tension test of porcine liver.

Table 2 :Figure 2 :
Figure 2: Geometry and boundary conditions in uniaxial compression test of silicon rubber.

Figure 3 :
Figure 3: Effect of modification on stress-stretch curve of tension test of porcine liver (experimental data from Gao et al. [17]).

Figure 5 :
Figure 5: (a) Normalized x coordinate of the free edge versus distance along this edge and (b)Path of free edge for tension test.

Figure 6 :
Figure 6: Effect of modification on the stress-stretch curve of compression test of unfiled silicon rubber (experimental data from Meunier et al. [12]).

Figure 7 :
Figure 7: Finite element simulations of compression test: (a) modified and (b) unmodified invariants model.

Figure 8 :
Figure 8: (a) Normalized z coordinate of the free edge and (b) path along free edge for compression test.