Sensitivity based virtual fields for identifying
hyperelastic constitutive parameters

Adel Tayeb, Jean-Benoit Le Cam, Eric Robin, F. Canévet, M. Grédiac, E.

Toussaint, X. Balandraud

► To cite this version:

Adel Tayeb, Jean-Benoit Le Cam, Eric Robin, F. Canévet, M. Grédiac, et al. Sensitivity based virtual fields for identifying hyperelastic constitutive parameters. 11th European Conference on Constitutive Models for Rubber, Jun 2019, Nantes, France. pp.163-168, 10.1201/9780429324710-29 . hal-02499976

HAL Id: hal-02499976

https://univ-rennes.hal.science/hal-02499976

Submitted on 7 May 2020

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sensitivity based virtual fields for identifying hyperelastic constitutive parameters

A. Tayeb, J.B. Le Cam & E. Robin

Univ Rennes, CNRS, IPR (Institut de Physique de Rennes) -UMR 6251, Rennes, France

LC-DRIME Joint Research Laboratory, Cooper Standard -CNRS -Univ Rennes

Campus de Beaulieu, Bât. 10B, Rennes Cedex, France

F. Canévet

Cooper Standard France, 194 route de Lorient, Rennes -France

LC-DRIME Joint Research Laboratory, Cooper Standard -CNRS -Univ Rennes

Campus de Beaulieu, Bât. 10B, Rennes Cedex, France

M. Grédiac, E. Toussaint & X. Balandraud

Université Clermont Auvergne, Institut Pascal, BP 10448, Clermont-Ferrand, France

CNRS, UMR 6602, IP, Aubière, France

**ABSTRACT:** 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.

# 1 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 4-branch 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 Ogden 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.

# 2 THEORETICAL BACKGROUND

Assuming a plane stress state and large strains, the principle of virtual work can be expressed as follows in the Lagrangian configuration:

$$ -e_0 \int_{s_0} \Pi(X, t) : \frac{\partial U}{\partial X}(X, t) dS_0 + e_0 \int_{\partial S_0} (\Pi \cdot N) \cdot U*(X, t) dL_0 = 0, \quad (1) $$
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.

## 2.1 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:

$$ \Pi = -pF^{-t} + \frac{\partial W}{\partial F} \quad (2) $$

where $p$ is an indeterminate coefficient due to incompressibility, $F$ is the deformation gradient tensor and $\bullet^t$ designates the transpose of a second-order tensor. For the Mooney model Mooney (1940), the strain energy density can be written as follows:

$$ W = c_1(I_1 - 3) + c_2(I_2 - 3), \quad (3) $$

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:

$$ c_1 \int_{S_0} \Theta : \frac{\partial U}{\partial \mathbf{X}} dS_0 + c_2 \int_{S_0} \Lambda : \frac{\partial U}{\partial \mathbf{X}} dS_0 = \int_{\partial S_0} (\Pi \cdot \mathbf{N}) \cdot \mathbf{U} * dL_0, \quad (4) $$

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 numplacememt virtual fields leads to the following linear system

$$ Ac = B $$

with

$$ A : \begin{bmatrix} \int_{S_0} \Theta : \frac{\partial U_{*}^{(1)}}{\partial X} dS_0 & \int_{S_0} \Lambda : \frac{\partial U_{*}^{(1)}}{\partial X} dS_0 \\ \int_{S_0} \Theta : \frac{\partial U_{*}^{(2)}}{\partial X} dS_0 & \int_{S_0} \Lambda : \frac{\partial U_{*}^{(2)}}{\partial X} dS_0 \end{bmatrix} \quad (5) $$

and

$$ c = \begin{cases} c_1 \\ c_2 \end{cases}, \quad B = \begin{cases} \int_{\partial S_0} (\Pi \cdot N) \cdot U_{*}^{(1)} dL_0 \\ \int_{\partial S_0} (\Pi \cdot N) \cdot U_{*}^{(2)} dL_0 \end{cases} $$

The two constitutive parameters $c_1$ and $c_2$ 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

$$ W = \sum_{i=1}^{N} \frac{2\mu_i}{\alpha_i^2} (\lambda_1^{\alpha_i} + \lambda_2^{\alpha_i} + \lambda_3^{\alpha_i} - 3), \quad (6) $$

where $\mu_i, \alpha_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

$$ \Pi_i = \frac{\partial W}{\partial \lambda_i} - \lambda_i^{-1} p. \quad (7) $$

The principle of virtual work of eq. (1) becomes in this case

$$ - \int_{S_0} (\Pi \cdot U_{u,u}^* + \Pi_2 \cdot U_{v,v}^*) dS_0 + \int_{\partial S_0} (\Pi \cdot N) \cdot U^*(X, t) dL_0 = 0, \quad (8) $$

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

$$ f(\chi) = \sum_{j=1}^{nV} \sum_{t=1}^{nTime} \sum_{i=1}^{nPts} \left[ \Pi_1(\chi) \cdot U_{u,u}^{*i(j)} + \Pi_2(\chi) \cdot U_{v,v}^{*i(j)} \right] S^i - W_{ext}^* \quad (9) $$

$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$.

## 2.2 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).

### 2.2.1 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.

### 2.2.2 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

$$ \delta\Pi^{(i)}(\chi, \mathbf{t}) = \Pi(\chi + \delta\chi_i, \mathbf{t}) - \Pi(\chi, \mathbf{t}), \quad (10) $$

where $0.1 \le \delta\chi_i \le 0.2 \chi_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

$$ \delta\Pi^{(i)}(\chi, \mathbf{t}) = B_{\text{glob}} U^{*(i)}, \quad (11) $$

where $B_{\text{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_{\text{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 $\bar{B}_{\text{glob}}$ is obtained from the original matrix $B_{\text{glob}}$. The virtual displacement field is then given by

$$ U^{*(i)} = \text{pinv}(\bar{B}_{\text{glob}}) \delta\Pi^{(i)}(\chi, t), \quad (12) $$

where pinv 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

$$ \frac{\partial U^{*(i)}}{\partial \mathbf{X}} = B_{\text{glob}} U^{*(i)}. \quad (13) $$

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)).

## 3 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><img src="image_2.png" /><figcaption>Figure 1. Sample geometry</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2. Experimental setup</figcaption></figure>

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 $\lambda_{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

<figure><img src="image_4.png" /><figcaption>Figure 3. Region of interest with 10 by 10 pixels Z.O.I</figcaption></figure>

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.

# 4 RESULTS

## 4.1 Experimental kinematic fields

The displacement field obtained from the SeptD software is smoothed using a mean centered filter. The values in Z.O.I.s 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.

## 4.2 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.

<figure><img src="image_5.png" /><figcaption>Figure 4. Experimental displacement gradient fields</figcaption></figure>
<table><caption>Table 1. Identified hyperelastic constitutive parameters using random virtual displacement fields</caption><thead><tr><th>Model</th><th>Parameters</th><th>Identified value</th></tr></thead><tbody><tr><td>Mooney</td><td>c1</td><td>0.229 MPa</td></tr><tr><td>—</td><td>c2</td><td>9.4 10<sup>-3</sup> MPa</td></tr><tr><td>Ogden 1</td><td>μ<sub>1</sub></td><td>0.68 MPa</td></tr><tr><td>—</td><td>α<sub>1</sub></td><td>1:4</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 5. Randomly generated virtual fields for the identification of the Mooney model</figcaption></figure>

## 4.3 Identification from sensitivity-based virtual fields

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 are obtained for several simulations with different sensitivity parameters, i.e. different virtual fields. Furthermore, the mean values for the parameters do not affect the

<table><caption>Table 2. Reference values for the constitutive parameters</caption><thead><tr><th>Model</th><th>Parameters</th><th>Identified value</th></tr></thead><tbody><tr><td>Mooney</td><td>c1</td><td>0.4 MPa</td></tr><tr><td>—</td><td>c2</td><td>0.04 MPa</td></tr><tr><td>Ogden 1</td><td>μ<sub>1</sub></td><td>0.768 MPa</td></tr><tr><td>—</td><td>α<sub>1</sub></td><td>1.26</td></tr></tbody></table>

final result of the identification. Indeed, the mean values for each parameter could change within the reference range without affecting the final result of the identification. For the Ogden model, the least square error (the value of the objective function at the end of the identification) is about 1.5 10⁻⁵.

<figure><img src="image_7.png" /><figcaption>Figure 6. Sensitivity based virtual fields for the identification of the Mooney model</figcaption></figure>

<table><caption>Table 3. Parameters identified using sensitivity-based virtual fields</caption><thead><tr><th>Model</th><th>Parameters</th><th>Identified value</th></tr></thead><tbody><tr><td>Mooney</td><td>c1</td><td>0.22 MPa</td></tr><tr><td>—</td><td>c2</td><td>1.9 10<sup>-2</sup> MPa</td></tr><tr><td>Ogden 1</td><td>μ<sub>1</sub></td><td>0.46 MPa</td></tr><tr><td>—</td><td>α<sub>1</sub></td><td>2.11</td></tr></tbody></table>

## 4.4 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,
<figure><img src="image_8.png" /><figcaption>Figure 7. Force obtained from finite element simulations compared to experimental force</figcaption></figure>

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.

## 5 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.

## 6 ACKNOWLEDGEMENTS

The authors thank M. Marek and Prof. Pierron from Southampton University for providing algorithms dealing with the sensitivity-based virtual fields approach. The authors thank the National Center for

Scientific Research (MRCT-CNRS and MI-CNRS), Rennes Metropole and Region Bretagne for financially supporting this work, as well as the Cooper Standard company for providing the specimens.

## REFERENCES

Grédiac, M. & F. Pierron (1998). At-shaped specimen for the direct characterization of orthotropic materials. *International Journal for Numerical Methods in Engineering* 41(2), 293-309.

Grédiac, M., E. Toussaint, & F. Pierron (2002). Special virtual fields for the direct determination of material parameters with the virtual fields method. I—principle and definition. *International Journal of Solids and Structures* 39 (10), 2691-2705.

Guélon, T., E. Toussaint, J.-B. Le Cam, N. Promma, & M. Grédiac (2009). A new characterization method for rubbers. *Polymer Testing* 28, 715-723.

Johlitz, M. & S. Diebels (2011). Characterisation of a polymer using biaxial tension tests. part i: Hyperelasticity. *Arch Appl. Mech.* 81, 1333-1349.

Marek, A., F. M. Davis, & F. Pierron (2017). Sensitivity-based virtual fields for the non-linear virtual fields method. *Computational Mechanics* 60(3), 409-431.

Marek, A., F. M. Davis, M. Rossi, & F. Pierron (2018). Extension of the sensitivity-based virtual fields to large deformation anisotropic plasticity. *International Journal of Material Forming*, 1-20.

Mooney, M. (1940). A theory of large elastic deformation. *Journal of applied physics* 11(9), 582-592.

Ogden, R. W. (1972). Large deformation isotropic elasticity—on the correlation of theory and experiment for incompressible rubberlike solids. *Proc. R. Soc. Lond. A* 326(1567), 565-584.

Pierron, F. & M. Grédiac (2012). *The virtual fields method: extracting constitutive mechanical parameters from full-field deformation measurements*. Springer Science & Business Media.

Promma, N., B. Raka, M. Grediac, E. Toussaint, J.-B. Le Cam, X. Balandraud, & F. Hild (2009). Application of the virtual fields method to mechanical characterization of elastomeric materials. *International Journal of Solids and Structures* 46(3-4), 698-715.

Seibert, H., T. Scheffer, & S. Diebels (2014). Biaxial testing of elastomers—experimental setup, measurement and experimental optimisation of specimen?*s shape. *Technische Mechanik*, 81 72-89.

Toussaint, E., M. Grédiac, & F. Pierron (2006). The virtual fields method with piecewise virtual fields. *International Journal of Mechanical Sciences* 48(3), 256-264.

Vacher, P., S. Dumoulin, F. Morestin, & S. Mguil-Touchal (1999). Bidimensional strain measurement using digital images. *Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering* 213(8), 811-817.

Yeoh, O. H. (1993). Some forms of the strain energy function for rubber. *Rubber Chemistry and Technology* 66 (5), 754-771.

Zienkiewicz, O. C., R. L. Taylor, O. C. Zienkiewicz, & R. L. Taylor (1977). *The finite element method*, Volume 36. McGraw Hill London.