Sensitivity based virtual fields for identifying hyperelastic constitutive parameters

In this study, the Virtual Fields Method (VFM) is applied to identify constitutive parameters of hyperelastic models from a heterogeneous test. Digital image correlation (DIC) was used to estimate both displacement and strain fields required by the identification procedure. Two different hyperelastic models were considered: the Mooney model and the Ogden model. Applying the VFM to the Mooney model leads to a linear system which involves the hyperelastic parameters due to the linearity of the stress with respect to these parameters. In the case of the Ogden model, the stress is a nonlinear function of the hyperelastic parameters and a suitable procedure shall be used to determine virtual fields leading to the best identification. This complicates the identification procedure and affects its robustness. This is the reason why the sensitivity based virtual field approach recently proposed in case of anisotropic plasticity by Marek et al. (2017) has been successfully implemented to be applied in case of hyperelasticity. Results obtained clearly highlight the benefits of such an inverse identification approach in case of non-linear systems.


INTRODUCTION
The constitutive parameters of hyperelastic models are generally identified from several homogeneous tests, typically uniaxial tension (UT), pure shear (PS) and equibiaxial tension (EQT). An alternative methodology consists in performing only one heterogeneous test. This is typically the case when a multiaxial loading is applied to a 3-branch, see Guélon et al. (2009) or a 4branch cruciform specimen, see Promma et al. (2009), Johlitz and Diebels (2011), Seibert et al. (2014). Indeed this induces a large number of mechanical states at the specimen's surface. The resulting heterogeneous strain fields are generally measured by Digital Image Correlation (DIC) technique. Among the different identification methodologies, the Virtual Field Method (VFM) has been successfully applied to hyperelasticity in Promma et al. (2009). In this work, linear systems were obtained due to the linearity of the stress with respect to the constitutive parameters (see the Mooney (1940) and the Yeoh (1993) models). When the system becomes nonlinear, typically for the Odgen model Ogden (1972), a statistical analysis can be carried out for optimizing the choice of the virtual fields. This complicates the identification procedure and affects its robustness. This is the reason why we apply here the sensitivity-based virtual field approach recently proposed in case of anisotropic plasticity by Marek et al. (2017). In the next section the theoretical background of VFM is presented and the sensitivity-based virtual field approach is more precisely detailed. Then, the experimental setup is briefly described. The results obtained for the two approaches are presented and discussed. Concluding remarks close the paper.

THEORETICAL BACKGROUND
Assuming a plane stress state and large strains, the principle of virtual work can be expressed as follows in the Lagrangian configuration: where Π is the first Piola-Kirchhoff stress tensor, e0 is thickness of the solid, S 0 is the surface of the solid in the normal direction to the thin dimension and ∂S 0 its boundary. X are the coordinates and N denotes the normal vector to the edge. S 0 and e 0 are measured in the reference configuration.

Hyperelasticity
For hyperelastic materials, the mechanical behavior is described by the strain energy density W relating the stress to the strain through the principle stretches (λ1,λ2,λ3) or the first two principal invariants of the left Cauchy-Green strain tensor (I1 and I2). Assuming that the material is incompressible, the first Piola-Kirchhoff stress tensor for such material reads: where p is an indeterminate coefficient due to incompressibility, F is the deformation gradient tensor and • t designates the transpose of a secondorder tensor. For the Mooney model Mooney (1940), the strain energy density can be written as follows: where c 1 and c 2 are the constitutive parameters to be identified. Combining eqs.
(2) and (3) and replacing Π by its expression in eq. (1) lead to the expression of the principle of virtual work: where Θ and Λ are two functions of the principle stretches. Using eq. (4) with two independent dis-nV F, nT ime and nP ts denote respectively the numplacement virtual fields leads to the following linear system The two constitutive parameters c1 and c2 are obtained by inversion of system. The second model considered in the present study is the Ogden (1972) one. In this case, the strain energy density reads where µ i ,α i ; i =1..N are the constitutive parameters governing the model. From eqs. (2) and (6) the eigenvalues of the Piola-Kirchhoff stress tensor are given by The principle of virtual work of eq. (1) becomes in this case where (u, v) is the principle basis of the strain tensor. In this basis, the cost function to be minimized to find the constitutive parameters can be written as follows nV F, nTime and nPts denote respectively the number of independent virtual fields, the time steps and the number of Z.O.I. In this paper, only the first order (N =1) of this model is considered. The identification of the constitutive parameters is performed by minimizing the cost function f.

Choice of the virtual displacement fields
Since an infinite number of kinematically admissible virtual fields U* satisfy the principle of virtual work in Eq. (1), the choice of a set of independent virtual fields remains a topical issue. A first attempt to fulfill this aim was to consider analytic forms of the virtual fields which can work very well for elastic materials, see Grédiac et al. (2002) and composite materials, see Pierron and Grédiac (2012) and Grédiac and Pierron (1998). Another authors considered piecewise displacement virtual fields such as in Toussaint et al. (2006).

Random virtual displacement fields
To deal with hyperelastic materials, different approaches in the generation of independent virtual displacement fields should be applied. To the authors' knowledge, the virtual fields method was first applied to hyperelastic materials in Promma et al. (2009). Motivated by a noise-sensitivity study, a set of random virtual displacement fields was generated. The procedure relies on the division of the region of interest in the sample into 12 quadrangular sub-domains over which piecewise virtual fields are defined. The displacement at the nodes are then generated and the set of virtual fields leading to the best identification is chosen. The displacement is approximated in each sub-domain by four-noded quadrangular finite elements Zienkiewicz et al. (1977). In Promma et al. (2009), these random virtual displacement fields were used for two hyperelastic models for which the application of the virtual fields method led to a linear system. In this case, to ensure the independence of the virtual fields, the criterion was a good conditioning of the system in eq. (5). However, for models for which the Virtual Field Method do not lead to a linear system, such as the Ogden model, a statistical study should be done to generate independent virtual fields, which makes identification more complicated and less robust, and requires the development of alternative strategies.

Sensitivity-based virtual displacement fields
In a recent work Marek et al. (2017), a new procedure for generating independent virtual displacement fields was employed for the identification of the constitutive parameters of an anisotropic plastic material in the small strain domain. The method is based on the sensitivity of the stress to small changes of constitutive parameters, typically between 10 and 20 %. The virtual displacement fields are then generated proportionally to the stress sensitivity fields through a finite element-like approach. The method was then extended to finite strain for anisotropic plasticity in Marek et al. (2018). In this case, the stress sensitivity field is defined by where 0.1 χ i ≤ δχ i ≤ 0.2 χ i is the sensitivity of the i th parameter with numerical values picked from literature. Note that the stress sensitivity field in Eq. (10) gives the influence of each constitutive parameter in the global response of the material at each point since the experiment used is heterogeneous. Therefore, the virtual displacement fields were generated proportionally by setting the stress sensitivity fields with the following expression where B glob is the global strain-displacement matrix from a virtual mesh generated a priori. This matrix is obtained by assembling the elementary strain-displacement matrix obtained directly from the derivation of the shape functions with respect to the coordinates. U *(i) in Eq. (11) designates the virtual displacement field corresponding to the i th constitutive parameter. Note that this virtual displacement field is a test function and has no physical meaning. In practice, matrix B glob should be modified to account for the boundary conditions of the region of interest (R.O. I). Typically, for edges where external loading is unknown a null displacement should be imposed. Therefore, a new matrix B ̅ glob is obtained from the original matrix B glob . The virtual displacement field is then given by where pinυ designates the pseudo inverse operator. Once the virtual displacement field is obtained, its gradient needed in the principle of virtual work is computed using the classic equation obtained with the finite elements method The contribution of each constitutive parameter to the response of the material is very different and unique. Therefore, a scaling in the cost function should be added (see Marek et al. (2017) and Marek et al. (2018)).

EXPERIMENTS
The material used in this study is a carbon black filled natural rubber. The sample is shown in Figure 1. It is a 105 mm long and 2 mm thick cruciform sample. This shape of the specimen gives, when an equibiaxial load is applied, several strain states, especially UT, PS and EQT, and a various states in between. Therefore, the single heterogeneous test used in this study can replace all the homogeneous tests classically used to identify hyperelastic constitutive parameters.
The experimental setup is presented in Figure 2. It is composed by a home-made biaxial testing machine and a digital camera. The four independent actuators were linked to have the same movement Figure 1. Sample geometry such that the specimen center was motionless during the test. Hence, a reference point is obtained in the center of the sample with respect to the correlation procedure. A displacement of 70 mm was applied to each branch at a loading rate of 150 mm/min, which corresponds to value of λ max of 3.4.
During the mechanical test, images of the specimen surface were stored at a frequency of 5 Hz using an IDS camera equipped with a 55 mm telecentric objective. The charge-coupled device (CCD) of the camera has 1920 × 1200 joined pixels. The displacement field at the surface of the specimen was determined using the digital image correlation (DIC) technique. The correlation process is achieved thanks to the SeptD software, see Vacher et al. (1999). The spatial resolution, defined as the smallest distance between two independent points, was equal to 10 pixels. A rectangular region defined in one branch of the specimen is sufficient to apply the identification procedure described in Section 2. The rectangular R.O.I is represented in Figure 3. Due to the large displacement applied during the experiment, the correlation could not be achieved in some Z.O.I (less than 4% of the total number of Z.O.I).
The displacement in these Z.O.I is approximated with a linear interpolation of the values obtained in the surrounding Z.O.I.s. The displacement fields are then smoothed using a mean filter in order to eliminate noise, especially where significant gradient occurs. This filter is applied before and after the differentiation step needed to obtain the displacement gradient tensor.

Experimental kinematic fields
The displacement field obtained from the SeptD software is smoothed using a mean centered filter. The values in Z.O.Is where the correlation could not be achieved were interpolated beforehand, as discussed above. The displacement gradient fields are presented in Figure 4 for the rectangular R.O.I. These data were smoothed using the same filter. The data obtained experimentally were used in the identification of the constitutive parameters for random and sensitivity based virtual displacement fields.

Identification from random virtual fields
The identification procedure described in 2.2.1 is applied herein for the determination of the hyperelastic constitutive parameters. For the Mooney model, the virtual displacement fields were chosen in such a way that the conditioning of the matrix A of Eq. (5) was greater than 0.3. For the Ogden model no criteria is found in the choice of the virtual displacement fields, hence, 100 virtual fields were generated and used in the identification procedure. The parameters identified using this approach are reported in table 1. The parameters found herein were then tested, at the end of this section, to reproduce the experiments.  The identification procedure presented in 2.2.2 is used to obtain the hyperelastic constitutive parameters. First, a virtual mesh has to be generated. It can be different from the correlation grid. Then, the virtual fields are generated proportionally to the stress sensitivity fields. The reference values for the parameters used in this work are reported in Table 2. The values of the parameters identified are reported in table 3. The parameters for all the models considered here are in agreement with the reference values. Note that, the parameters in Table 3

Results comparison
To evaluate the accuracy of the identified parameters, the biaxial experiment used in this work was simulated by using the Abaqus FE package. A plane stress problem was considered with the parameters given in Tables 1 and 3. The mesh used in the model was CPS4. This element is a four node bilinear plane stress quadrilateral. The four edges of the cruciform sample were subjected to a traction displacement of 70 mm each. The displacement was blocked in the shear direction along the edges of the sample. For each set of parameters, the resulting force was recorded in each branch of the sample and compared to the experimental load obtained in the experiment. The results of this comparison are shown in figure 7, where SBVF and RVF refer to sensitivity based virtual fields and random virtual fields, respectively. For the RVF method, the Mooney model appears to have a good result for a maximum stretch up to 2.7, which is the usual range for the Mooney model. However, the Ogden model overestimates the force in the branch for the whole experiment. This is due to the choice of the virtual displacement fields, which was done randomly. For the SBVF method, the Mooney model has a good prediction for the experimental force for a maximal stretch up to 2.6. The Ogden model has a better result for wider strain range corresponding to a maximal stretch up to 3. The results of these two models are very satisfactory given that they do not take into account the stress hardening phenomenon. The second order Ogden model predicts very well the experimental force for the whole experiment. Hence, the capacity of the SBVF method in the generation of the virtual displacement fields is illustrated here in the case of hyperelastic behavior.

CONCLUSION
In this study, the Virtual Fields Method (VFM) was applied to identify constitutive parameters of hyperelastic models from a heterogeneous test in the cases of linear and non-linear relationships between the stress and the constitutive parameters to be identified. In the former case, the Mooney model was considered and virtual field were randomly generated. For the latter case, the Ogden model was used and a sensitivity-based virtual fields approach inspired from a recent work due to Marek et al. (2017) for anisotropic plasticity was applied to choose the virtual fields. Results obtained with the two approaches clearly highlight the benefits of using the sensitivity-based virtual fields approach for identifying the constitutive parameters in case of non-linear systems.