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

HAL
open science

# Inverse identification of constitutive parameters from heat source fields: A local approach applied to hyperelasticity

Sylvain Charles, Jean-Benoit Le Cam

► To cite this version:

Sylvain Charles, Jean-Benoit Le Cam. Inverse identification of constitutive parameters from heat source fields: A local approach applied to hyperelasticity. Strain, 2020, 56 (2), 10.1111/str.12334. hal-02469675

HAL Id: hal-02469675

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

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Inverse identification of constitutive parameters from heat source fields: A local approach applied to hyperelasticity

Sylvain Charlès, Jean-Benoît Le Cam

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

**Funding information**

PCM Technologies S.A.S; Région Bretagne; Rennes Métropole

In this paper, a new inverse identification method of constitutive parameters is developed from full kinematic and thermal field measurements. It consists in reconstructing the heat source field from two different approaches by using the heat diffusion equation. The first one requires the temperature field measure-ment and the value of the thermophysical parameters. The second one is based on the kinematic field measurement and the choice of a thermo-hyperelastic model that contains the parameters to be identified. The identification is carried out at the local scale, ie, at any point of the heat source field, without using the boundary conditions. In the present work, the method is applied to the challeng-ing case of hyperelasticity from a heterogeneous test. Due to large deformations undergone by the rubber specimen tested, a motion compensation technique is developed to plot the kinematic and the thermal fields at the same points before reconstructing the heterogeneous heat source field. In the present case, the constitutive parameter of the Neo-Hookean model has been identified, and its distribution has been characterized with respect to the strain state at the surface of a cross-shaped specimen.

**KEYWORDS**

digital image correlation, elastomer, heat source reconstruction, hyperelasticity, infrared thermog-raphy, inverse identification

## 1 | INTRODUCTION

Several methods have been recently developed for identifying constitutive parameters from field measurements. The measurement techniques employed generally deal with digital image correlation (DIC), Moiré and speckle interferometry, and grid methods. As these techniques provide a large number of experimental data, they are coupled with suitable computational strategies for analyzing the data, typically the finite element model updating method, the constitutive equation gap method, the virtual fields method, the equilibrium gap method, and the reciprocity gap method. They have been fully reviewed in Avril et al.[1] In many of these approaches, the boundary conditions are required to solve the identification problem. In the cases dealing with hyperelasticity, the boundary conditions and the measurement of the global quantities were needed in the inverse identification methodology (see for instance Promma et al[2] for the virtual fields method and Giton et al[3] for the finite element model updating method). The present study aims at developing a methodology for inverse identification from local quantities only, ie from a zone at the specimen surface and without knowing the boundary conditions of the problem. This implies that the local quantities explicitly depend on the constitutive parameters, ie, on the strain-stress relationship. In this work, we propose to identify the constitutive parameters by reconstructing the
heat source* field from two different ways: (a) the kinematic field and a given thermomechanical model that contains the parameters to be identified, and (b) the temperature field and the thermophysical parameters. This technique is well suitable for materials with strong thermoelastic couplings and low thermal conductivity such as elastomers. In Section 2, the mechanical behavior of rubber-like materials and mechanical tests classically performed for identifying the constitutive parameters are briefly recalled. Then, the theoretical framework used for reconstructing the bidimensional heat source field from thermal and kinematic field measurements is detailed in Section 3. The experimental setup is presented in Section 4, where the full thermal and kinematic fields measurement, the motion compensation technique developed to track the temperature variation during the large displacements of the material points in the IR images, and the determination of the surface emissivity of the material are precisely detailed. Section 5 gives the results obtained in terms of the kinematic and thermal fields, the loading cases map, and the heat source field reconstructed from the two field measurement techniques. The constitutive parameters are identified according to different identification strategies. Concluding remarks close the paper.

## 2 | MECHANICAL BEHAVIOR OF RUBBER-LIKE MATERIALS AND CLASSICAL TESTS FOR IDENTIFICATION

The mechanical response of rubber-like materials is governed by numerous physical phenomena: rate-dependent effects,[4] permanent set,[5] stress softening,[6–10] mechanical hysteresis,[6,9] strain-induced crystallization (SIC),[11–14] energy storage,[15–17] anisotropic effects,[5,9,18–20] change in volume,[21] and the competition between entropy and isentropy at low strains.[22] It should be noted that most of these phenomena are induced or amplified by adding fillers in the rubber matrix. Despite this, the framework of hyperelasticity is generally used as a first approximation to predict the mechanical response of these materials.[23,24] Numerous constitutive relations are available in the literature and have been reviewed and compared in Marckmann and Verron.[25] As these hyperelastic models do not account for the above-mentioned phenomena, the values of the constitutive parameters strongly depend on the strain state.[26] This is the reason why constitutive parameters are classically identified from several homogeneous tests, namely, the uniaxial tensile (UT), the pure shear (PS), and the equibiaxial tensile (EQT), or uniaxial compression. These three tests completely describe the domain of possible loading paths.[22,27–29] A trade-off between the sets of values obtained with the different tests has therefore to be found to obtain parameters that can reasonably be considered as intrinsic to the mechanical behavior of the material. Such identification approach has several disadvantages: (a) Each of the tests carried out is assumed to induce a homogeneous strain state, which is a strong assumption in case of PS and EQT tests. (b) Several specimen geometries, ie molds, are required. (c) Several testing machines are needed to apply these loading conditions. (d) Dispersion obtained for each test requires testing several specimens for each loading condition. (e) The comparison between the constitutive parameters identified from different loadings is a matter of debate. For instance, should the identification be made for a same maximum stretch applied, the same energy, etc? (f) The elaboration process may differ from one specimen geometry to another (typically compression molding vs. injection molding). This can have a significant effect on the value of the identified constitutive parameters.

Recently, an alternative approach has been proposed. It is based on a heterogeneous test that induces the three tests mentioned above and a wide range of intermediary loadings. This approach is further detailed in Promma et al[2] and Guélon et al.[30] Such a heterogeneous biaxial test was also explored in Johlitz and Diebels[31] and Seibert et al.[32] In these studies, several strategies have been used for the inverse identification of the constitutive parameters, either from the global quantities (force and/or displacement) or by coupling these global quantities with local ones, typically the kinematic field measured by using the DIC technique. In every case, the boundary conditions were required.

## 3 | HEAT SOURCE APPROACH

In this section, the suitable formulation of the heat diffusion equation is first presented. Then, the methodology applied for the heat source field reconstruction from two experimental measurements, the temperature and the kinematic fields, is explained.

<footnote>*The term “heat source” is used in this paper to mean the heat power density that is produced or absorbed by the material.</footnote>
## 3.1 | Heat diffusion equation

Considering that the constitutive state equations are derived from the Helmholtz free energy function and that the heat conduction followed Fourier's law, the local heat diffusion equation writes in the reference (Lagrangian) configuration[33]:

$$ \rho_0 C \dot{T} - \text{Div}(\mathbf{K}_0 \text{ Grad } T) - R = D_{\text{int}} + T \frac{\partial \mathbf{P}}{\partial T} : \dot{\mathbf{F}} + T \sum_{\beta=1}^{m} \frac{\partial \mathbf{A}_{\beta}}{\partial T} : \dot{\xi}_{\beta}, \quad (1) $$

where $\rho_0$ is the density, C is the heat capacity, $\mathbf{K}_0$ is the thermal conductivity tensor, T is the absolute temperature, R is the external heat source (radiation, for instance), and $D_{\text{int}}$ is the intrinsic dissipation. It is also named mechanical dissipation. P is the first Piola–Kirchhoff stress tensor (or nominal stress tensor). F is the deformation gradient tensor. S denoted the heat source. The superposed dot indicates the material time derivative.

In the right member, the term $T \frac{\partial \mathbf{P}}{\partial T} : \dot{\mathbf{F}}$ corresponds to the heat sources due to the thermo-elastic couplings.† The term $T \frac{\partial \mathbf{A}}{\partial T} : \dot{\xi}_{\beta}$ corresponds to the other thermomechanical couplings. $\xi_{\beta}$ are the internal tensorial variables. $\mathbf{A}_{\beta}$ are the thermodynamic forces associated with the internal tensorial variables.

## 3.2 | Heat source field reconstruction from IR measurements

IR thermography provides full surface temperature fields. Therefore, reconstructing the heat source field requires the formulation of a two-dimensional version of the heat diffusion equation.

For that purpose, several assumptions are used. First, the heat conduction is considered as isotropic, which leads to $K_{0ij} = k_0 \delta_{ij}$, where $k_0$ is the thermal conductivity coefficient. Consequently, the heat diffusion equation (1) can be written as follows:

$$ \rho_0 C \dot{T} - k_0 \Delta T - R = S, \quad (2) $$

where $\Delta$ is the Laplacian operator in the Lagrangian configuration.

Second, by considering the temperature as homogeneous through the specimen thickness, the surface temperature is very close to the mean temperature through the thickness and Equation (2) can be integrated over the thickness. This means that the specimen must be thin and that the plane stress state can be generally assumed. This leads to the two-dimensional formulation of the heat diffusion equation:

$$ \rho_0 C \left( \dot{T} + \frac{T - T_{amb}}{\tau_{2D}} \right) - k_0 \Delta_{2D} T - R = S, \quad (3) $$

where $T_{amb}$ is the ambient temperature and $\Delta_{2D}$ is the Laplacian operator in the specimen plane in the Lagrangian configuration. $\tau_{2D}$ is a time characterizing the heat exchanges along the Z-direction by convection with the air at the specimen's surface. It should be noted that the relevancy of this assumption for a given situation can be evaluated by calculating the Biot number $B_i$. For instance, in the case of a rectangular thin specimen, the Biot number is calculated as follows:

$$ B_i = \frac{h}{k_0} \cdot \frac{e}{2}, \quad (4) $$

where $h$ is the convection coefficient and $e$ is the thickness. Below 0.1, the temperature variation in the thickness is less than 5%.[34] which is typically the case for a 2-mm thick rectangular rubber specimen.

Third, the external radiations (R) are difficult to measure. Nevertheless, if they remain constant during the deformation process, and by using the temperature variation $\theta(X, Y, t) = T(X, Y, t) - T_{ref}(X, Y)$ instead of the temperature itself, the external radiations can be removed from the heat diffusion equation. This implies specific preparation of the thermal measurement and the specimen's environment. It should be noted that $T_{ref}$ is the ambient temperature if it is constant

<footnote>†The term thermo-elastic couplings means here the couplings between temperature and strain, whatever the origin of the elasticity is entropic, non-entropic or the both.</footnote>
over the test. In this case, Equation (3) writes

$$ \rho_0 C \left( \dot{\theta} + \frac{\theta}{\tau_{2D}} \right) - k_0 \Delta_{2D} \theta = S. \qquad (5) $$

This two-dimensional formulation of the heat diffusion equation has been recently applied to the thermomechanical characterization of the influence zone of crack tips in filled rubber specimens.[35,36] Note that elastomeric materials are generally assumed to be incompressible, ie, $\rho_0 = \rho$. Finally, if the heat conduction is negligible in the specimen plane (see Remark 1 hereafter) or if the heat exchanges with the grips of the testing machine meet certain conditions (see Remark 2 hereafter), the two-dimensional heat diffusion equation can be simplified[37,38] and can correspond to a “0D” formulation of the heat diffusion equation:

$$ \rho C \left( \dot{\theta} + \frac{\theta}{\tau} \right) = S. \qquad (6) $$

Such simplification has previously been done for metallic materials. The reader can refer to previous works[39–43] for further information. This equation was more recently applied for characterizing the calorimetric signature of the different mechanisms in play in the homogeneous deformation of rubber.[44–47] In the present case, we consider that $\tau_{2D} = \tau_{0D}$.[4] This is a reasonable hypothesis since the material is quasi-incompressible and undergoes large strains and that most of the points in the measurement area are far from the grip. Therefore, the change in thickness is the preponderant effect on $\tau$ (see Remark 2 and Equation 17).

**Remark 1.** The effect of the in-plane heat conduction has to be discussed according to the spatial resolution of the temperature measurement. A thermal diffusion length can be calculated from the thermal diffusivity $D$ and the frequency $f_L$ as $\sqrt{\frac{D}{\pi f_L}}$. In the case of unfilled rubber, the thermal diffusion length for a test which one cycle duration is 60 s is about 1.3 mm. With an IR pixel with a size of 300 µm, this means that the heat diffusion affects the temperature measurement if the stress varies significantly within a windows of about 5 × 5 pixels. In the case where the material has a bad thermal conductivity, typically an elastomer, and that the stress gradient is low, the in-plane heat conduction is negligible.

**Remark 2.** If the heat exchanges with the grips of the testing machine are proportional to the average temperature change of the specimen, which is nearly the case when heat sources are homogeneous (see section 6.1.2 in Chrysochoos and Louche[39]), resulting in a “total” heat exchange (including convention with air) that is also proportional to the average temperature change of the specimen (see Jongchansitto et al[48] for further details).

## 3.3 | Heat source prediction from the kinematic field

The prediction of the heat sources evolution during the deformation process requires the choice of a free energy density $\Psi$. Here, the material is assumed to behave as a hyperelastic material that is mechanically incompressible and isotropic. At low strain levels, typically inferior to 250%, the neo-Hookean model can be chosen to predict the mechanical behavior.[49] The free energy function is then given by the following strain energy density:

$$ \rho \Psi(\mathbf{F}, T) = W(\mathbf{F}, T) = \frac{1}{2} N k T \quad (I_1 - 3), \qquad (7) $$

where $N$ is the number of network chains per unit volume, $k$ is the Boltzmann constant, $T$ is the temperature in kelvin, and $I_1$ is the first invariant of the left Cauchy–Green deformation tensor $FF^T$.

If the material is mechanically incompressible, the deformation gradient tensor writes in the principal directions $e_I$, $e_{II}$, and $e_{III}$:

$$ \mathbf{F} = \lambda \quad e_I \otimes e_I + \lambda^B e_{II} \otimes e_{II} + \lambda^{-(B+1)} e_{III} \otimes e_{III}, \qquad (8) $$

where $\lambda$ is the stretch ratio in the direction $e_I$, defined as the ratio of the current length to the initial length in this direction, the operator $\otimes$ between two vectors is such that $[a \otimes b]_{ij} = a_i b_j$, and $B$ is the biaxiality ratio.[50,51]

<footnote>‡For sake of simplicity, $\tau_{0D}$ is denoted $\tau$.</footnote>
According to Equation (8), the strain energy density reads

$$W(\lambda, T) = \frac{1}{2} N k T [\lambda^2 + \lambda^{2B} + \lambda^{-2(B+1)} - 3]. \qquad (9)$$

In the case where the material does not produce any intrinsic dissipation and no other thermomechanical couplings come into play (due to strain-induced crystallization for instance), the heat source is given by

$$S = T \frac{\partial \mathbf{P}}{\partial T} : \dot{\mathbf{F}}. \qquad (10)$$

Using Equation (9), the heat source finally writes

$$\rho C \left(\dot{\theta} + \frac{\theta}{\tau}\right) = NkT \left[ (\lambda + B\lambda^{2B-1} - (B+1)\lambda^{-2B-3}) \frac{d\lambda}{dt} + \dot{B}ln(\lambda)(\lambda^{2B} - \lambda^{-2(B+1)}) \right]. \qquad (11)$$

Here, we recall that T is the value of the temperature in kelvin. For temperature variations that do not exceed a few degrees, NkT remains therefore nearly equal to NkT^ref, where T^ref is the initial temperature and does not affect significantly the heat source (0.34% for a temperature increase equal to 1 K from 293 K). This is the reason why NkT^ref is considered instead of NkT in the present study. For the sake of simplicity, the constant 1/2 NkT^ref is denoted Cnh in the following.

If the biaxiality coefficient remains constant during the test, ie $\dot{B} = 0$, the equation can be simplified to

$$\rho C \left(\dot{\theta} + \frac{\theta}{\tau}\right) = NkT \left(\lambda + B\lambda^{2B-1} - (B+1)\lambda^{-2B-3}\right) \frac{d\lambda}{dt}. \qquad (12)$$

# 4 | EXPERIMENTAL SECTION

## 4.1 | Material and specimen geometry

In the present study, an unfilled nitrile rubber was considered. Such a material does not produce detectable intrinsic dissipation and does not crystallize under tension. Therefore, the heat source writes only from the coupling between

<figure><img src="image_3.png" /><figcaption>FIGURE 1 Specimen geometry</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 2 Overview of the experimental setup</figcaption></figure>

temperature and strain, and Equation (11) applies. The specimen geometry is shown in Figure 1. Its total length is 125 mm and its thickness is 2 mm.

## 4.2 | Loading conditions

Figure 2 presents an overview of the experimental setup composed of an optical camera and an infrared camera, on both side of a home-made biaxial testing machine. The machine is composed of four independent RCP4-RA6C-I-56P-4-300-P3-M (IAI) electrical actuators. They were driven by a PCON-CA-56P-I-PLP-2-0 controller and four PCON-CA (IAI) position controllers. The actuators were piloted by an in-house LabVIEW program. Two load cells measure the force in the two perpendicular directions. In the present study, an equibiaxial loading was carried out. The four independent actuators were controlled in a way that the specimen's center was motionless. The displacement and loading rate were set at 70 mm and 150 mm/min, respectively, for each of them. It should be noted that the two cameras are triggered, and the images, the displacements of the grips, and the forces applied in the two perpendicular directions are stored at the same rate.

## 4.3 | Full kinematic field measurement

Images of the specimen surface at increasing stretches were stored at a frequency equal to 5 Hz with an IDS camera equipped with a 55-mm telecentric objective. The charge-coupled device (CCD) of the camera has 1,920 × 1,200 joined pixels. The displacement field at the specimen surface was determined by using the DIC technique. It consists in correlating the grey levels between two different images of a given zone, each image corresponding to different strain levels. In order to improve the image contrast, a white paint is sprayed on the surface before testing specimens. This leads to a black and white random grey field. A uniform cold lighting at the specimen surface is ensured by a home-made LED lamp. The camera was fixed on a multidirectional adjustable support. In this configuration, a region of interest (ROI) of 186 × 116 mm² is observed by the digital camera. The software used for the correlation process was SeptD.[52] The spatial resolution, defined as the smallest distance between two independent points, was equal to 10 pixels, which corresponds to 0.97 mm.

## 4.4 | Deformation gradient tensor computation

The $F_{ij}$ components of the deformation gradient tensor are determined in the centre of the square elements§ defined in the underformed DIC grid. Within each square element $i$, the horizontal and vertical displacements, $iU$ and $iV$ respectively, are assumed to be bilinear functions of the coordinates $(X, Y)$:

$$ \left\{ \begin{array}{l} iU(X, Y) = a + bX + cY + dXY \\ iV(X, Y) = e + fX + gY + hXY \end{array} \right. , \qquad (13) $$

<footnote>§One square element is formed from four ZOI, each corner of the square element corresponds to a ZOI centre</footnote>
where a, b, c, d, e, f, g, and h are constants that can be obtained from the ZOI values. The displacement is interpolated at the centre of each square element. The calculation is carried out by considering the element centre is the origin of the coordinate system. Finally, the deformation gradient tensor F components are calculated using Equation (14).

$$F = \begin{pmatrix} \frac{\partial x}{\partial X} & \frac{\partial x}{\partial Y} \\ \frac{\partial y}{\partial X} & \frac{\partial y}{\partial Y} \end{pmatrix} = \begin{pmatrix} \frac{\partial X + iU}{\partial X} & \frac{\partial X + iU}{\partial Y} \\ \frac{\partial Y + iV}{\partial Y} & \frac{\partial Y + iV}{\partial Y} \end{pmatrix} = \begin{pmatrix} 1 + b & c \\ f & 1 + g \end{pmatrix}. \quad (14)$$

The three principal stretches ($\lambda_1 > \lambda_2 > \lambda_3$) are defined as the square roots of the eigenvalues of the left Cauchy-Green tensor $B$ ($B = F F^t$). For a 2D measurement where only $\lambda_1$ and $\lambda_2$ are determined, $\lambda_3$ is deduced by assuming that the material is incompressible, ie, $J = \det F = \lambda_1 \lambda_2 \lambda_3 = 1$. The biaxiality ratio $B$, defined as the ratio between the logarithm of $\lambda_2$ and the logarithm of $\lambda_1$, is then calculated. It is equal to 0 for PS, 1 for EQT, and -0.5 for UT. The two invariants $I_1$ and $I_2$ of $B$ are written as a function of the principal stretch ratios $\lambda_i$ as follows:

$$\left\{ \begin{aligned} I_1 &= \lambda_1^2 + \lambda_2^2 + \lambda_3^2 \\ I_2 &= \lambda_1^2 \lambda_2^2 + \lambda_2^2 \lambda_3^2 + \lambda_3^2 \lambda_1^2 . \end{aligned} \right. \qquad (15)$$

These invariants can be expressed with respect to the maximum principal stretch ($\lambda_1 = \lambda$) and $B$:

$$\left\{ \begin{aligned} I_1 &= \lambda^2 + \lambda^{2B} + \lambda^{-2(B+1)} \\ I_2 &= \lambda^{2(1+B)} + \lambda^{-2} + \lambda^{-2B} . \end{aligned} \right. \qquad (16)$$

These quantities are used for characterizing the heterogeneity of the stretch states and the distribution in the maximum principal stretch value in the ROI observed by the optical camera.

## 4.5 | Full thermal field measurement

Temperature measurements were performed by using a FLIR infrared camera equipped with a focal plane array of 640 × 512 pixels and detectors operating in wavelengths between 1.5 and 5.1 µm. The integration time was equal to 2,700 µs, and the acquisition frequency was equal to procedure 25 Hz. The calibration of camera detectors was performed with a black body using a one-point NUC procedure at this acquisition frequency. The thermal resolution or noise equivalent temperature difference (NETD) is equal to 20 mK for a temperature range between 5 and 40°C. The spatial resolution of the thermal field was equal to 300 µm/px. The infrared camera is switched on several hours before testing in order to ensure its internal temperature to be stabilized. The infrared camera was trigged at the same frequency as the optical camera, ie 5 Hz. To measure the temperature field with an IR camera, the emissivity of the material must be evaluated.

<figure><img src="image_5.png" /><figcaption>**FIGURE 3** Scheme and image of the calibration pattern with the optical camera (on the left-hand side) and with the IR camera (on the right-hand side)</figcaption></figure>
In the present case, the specimen has been placed in an oven and heated at 50°C. The emissivity has been evaluated by comparing the temperature of the specimen with the one of a black body regulated at the same temperature. The emissivity obtained was equal to 0.94. This is in a good agreement with the emissivity measured or assumed for rubber surface in the literature.

**4.6 | Movement compensation technique**

Due to large deformations undergone by the material, the material points observed by the IR camera move from pixel to pixel in the IR images. The temperature variation at a given material point has therefore to be processed by compensating its movement within the IR images. This requires first describing the kinematic and the thermal fields in the same coordinate system. Several movement compensation techniques have been proposed in the case of elastomers: from infrared images only in Pottier et al.[53] by coupling full infrared and displacement fields in Samaca Martinez et al.[35,36] It should be noted that such technique has also been developed for other materials; see for instance Chrysochoos et al.[54] and Bodelot et al.[55] for steel, or Louche et al.[56] for shape memory alloys. In this paper, coupled full thermal and kinematic fields measurements have been carried out on both sides of the specimen. A calibration test pattern is placed in one of the grips (see Figure 3). It has two sets of four holes that form a large (ABCD) and a small rectangle. The rectangle used for the calibration is chosen with respect to the image resolution. In the present case, the calibration is done with holes A, B, C, and D. In the optical and IR images, the coordinates of each hole are determined with a mark tracking algorithm after image thresholding. In this figure, the centre of the holes is highlighted by a white cross symbol. Then, one of the holes is chosen as the reference one in order to plot optical and IR images in the same coordinate system.

Thanks to this common coordinate system, each point of the DIC grid where the deformation gradient tensor is calculated is plotted in the infrared image. As the resolution of the optical and IR images are not the same and the displacement

<figure><img src="image_6.png" /><figcaption>**FIGURE 4** Sample used to validate the motion compensation technique with the 3 points followed</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**FIGURE 5** Temperature evolution at the three points in Figure 4 during a 70-mm upwards and downwards specimen translation</figcaption></figure>
is not an entire number of pixels, the temperature of this point is interpolated from the four closer IR pixels. This step is reconducted for each displacement level applied to the specimen branches. This methodology has been validated with a test where the displacement of a rigid body has been studied. In this case, the temperature of points at the surface should not change during the displacement. Three points have been followed as shown in Figure 4. The specimen used is a cross-shaped specimen where two branches were cut. A translation of 70 mm of the grip is imposed at a loading rate equal to 150 mm/min. The results obtained are presented in Figure 5. This figure shows that the temperature tracked at each point does not change. The standard deviation is about 0.019°C, which corresponds to the IR camera noise equiva- lent temperature difference. This validates the motion compensation technique proposed. It should be noted that optical effects such as narcissus effect or halo observed using a Cedip Jade III-MWIR camera for instance in Pottier et al[53] does not affect significantly the temperature measurement in the present case.

5 | RESULTS

In this section, the full kinematic and thermal fields are first presented. Then, the heat source field reconstruction in the DIC grid in the undeformed configuration is explained. Finally, the different identification strategies applied are detailed.

**5.1 | Full kinematic field and heterogeneity induced by the mechanical test**

The interpolated displacement fields are plotted on the DIC grid. Figures 6a,b illustrate the displacement field obtained for the maximum loading applied. In these figures, a 3D representation of the displacement fields is given.

<figure><img src="image_8.png" /><figcaption>(a) Horizontal displacement field</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Vertical displacement field</figcaption></figure>

<caption>**FIGURE 6** Displacement field for the maximum loading applied (U = 70 mm in the four branches)</caption>
The deformation gradient tensor components are determined by using the methodology described in Section 4.3, which enables us to investigate the heterogeneity induced by the mechanical test. To do so, different representations are used, as explained in Promma et al[2]:

* Each ZOI in the ROI can be plotted as a point in the ($I_1, I_2$) plane (see Figure 7a). $I_1$ and $I_2$ are calculated by using Equation (16). Note that in this representation, the EQT, PS, and UT states correspond to three curves meeting at point (3, 3) corresponding to the undeformed state. The diagram has been established for the maximum displacement applied. Thus, all the points obtained for the different displacement levels applied are located below the dotted line in red color and between the UT and EQT curves. This diagram highlights that the three strain states are gathered in the strain field and that the strain range (represented by the distance between the origin (3,3) and the dotted line) is larger for UT than for EQT. This is logical as the EQT is obtained by stretching the four branches. It should be noted that specimens with shorter branches reduce this difference.[57]

<figure><img src="image_10.png" /><figcaption>**FIGURE 7** Characterization of the strain state heterogeneity</figcaption></figure>
* The stretch states are mapped over the ROI (see Figure 7b). A color code is employed. The EQT, PS, and UT states are represented by the blue, green, and red colors, respectively. Intermediate states between two pure states are represented by intermediate colors using a weighting coefficient between two colors defined by using the lowest difference between the two corresponding values of the biaxiality ratio. This figure highlights the spatial distribution of the strain states: EQT in a large zone at the specimen centre, UT in the branches, and PS in a ring-shaped zone surrounding the EQT zone. The latter zone is the smallest one.

* The maximum principal stretch is mapped over the ROI (see Figure 7c). This representation is complementary to the first one dealing with the $(I_1, I_2)$ plane.

It should be noted that the biaxiality coefficient does not vary during the test as soon as the global displacement level is sufficient to avoid buckling (at the end of the unloading) and a bad signal to noise ratio (at the beginning of the loading). This is highlighted in Figure 8, where the biaxiality coefficient is plotted with respect to time for five different points of the ROI. For such an equibiaxial test, $\dot{B} = 0$ and Equation (12) can therefore be applied.

## 5.2 | Full thermal field

The motion compensation technique has been applied to track the temperature variation of each point of the DIC grid where the deformation gradient tensor is determined. The first column in Figure 9 presents the thermal images where the grid of the interpolated displacement has been projected for the undeformed state and two global displacement levels (35 and 70 mm). The temperature field is then plotted in the undeformed grid of the interpolated displacement in the visible image as shown in the second column.

<figure><img src="image_11.png" /><figcaption>**FIGURE 8** Time evolution of biaxiality coefficient during the biaxial loading</figcaption></figure>
## 5.3 | Heat source reconstruction

The reconstruction of the heat source field requires the characterization of $\tau$ in Equation (6) according to the stretch and the biaxiality coefficient. $\tau$ is evaluated by heating the specimen and measuring the temperature field during the return to thermal equilibrium, then by fitting the curve by an exponential function, as shown in Figure 10 for the undeformed state. In this case, the average value of $\tau$ is denoted $\tau_0$ and is equal to 198 s. For the different displacement applied, $\tau$ spatially varies according to the change in the thickness, ie, to the stretch and to the biaxiality ratio ($\tau = \tau(\lambda, B)$). Either $\tau(\lambda, B)$ is evaluated for different loading levels or by considering that the link between the thickness changes and the biaxiality ratio and by assuming the material to be incompressible:

$$ \tau(\lambda, B) = \tau_0 \lambda^{-B-1}. \qquad (17) $$

The accurate evaluation of $\tau(\lambda, B)$ is all the more important that the test is performed at a low loading rate, leading to significant heat diffusion by convection and therefore to low temperature variations. Figure 11 presents the surface obtained by plotting $\tau(\lambda, B)$ from the measurement of $\tau_0$. The black crosses correspond to a measurement performed at the maximum displacement applied to the four specimen branches. This figure highlights the good estimation of $\tau(\lambda, B)$ by Equation (17).

Once $\tau(\lambda, B)$ is characterized, the heat source field can be reconstructed by applying Equation (5). As an illustration, Figure 12 provides the heat source field obtained for the maximum displacement applied, plotted in the undeformed DIC grid. The heat sources produced are the lowest in the EQT zone. This is explained by the fact that the maximum principal stretch and therefore the corresponding stretch rate reached in this zone are the lowest.

## 5.4 | Identification of the constitutive parameter

The constitutive parameter can be determined by the least square method. Two different approaches have been considered:

<figure><img src="image_12.png" /><figcaption>FIGURE 9 Results of the motion compensation technique developed</figcaption></figure>
* at the global scale, by assuming a unique value of the constitutive parameter for the whole ROI. In that case, the objective function to be minimized is defined by $\Phi(C_{nh}) = \sum_{i} \sum_{t} (S(i, t) - S_{nh}(i, t, C_{nh}))^2$, where $i$ is the number of the square elements, $t$ is the time, $C_{nh}$ is the neo-Hookean constitutive parameter to be identified, $S_{nh}$ is the heat source obtained using the neo-Hookean model, and $S$ is the heat source obtained from the temperature field by applying Equation (6). With this approach, $C_{nh}$ is equal to 0.204 MPa for the deformation state corresponding to the maximum loading applied (70 mm per branch). Calculations were also made for the whole load-unload cycle. The average value of $C_{nh}$ was

<figure><img src="image_13.png" /><figcaption>**FIGURE 10** Return to thermal equilibrium for determining $\tau$. In this figure, temperature is measured at a pixel in the specimen center</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**FIGURE 11** $\tau$ versus the stretch $\lambda$ and the biaxiality coefficient $B$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**FIGURE 12** Heat field in the undeformed state for the maximum displacement applied</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIGURE 13 Field of neo-Hookean parameter computed</figcaption></figure>

equal to 0.213 MPa, with a standard deviation equal to 0.0109 MPa. Due to buckling and bad signal-to-noise ratio, this calculation excluded the measurements performed for a global displacement inferior to 17 mm.

* at the local scale, by considering each square element separately and by minimizing the difference between the heat source obtained by the thermal and the kinematic fields measured at any point. This leads to an objective function defined by $\Phi(i, C_{nh}) = \sum_{t} (S(i, t) - S_{nh}(i, t, C_{nh}))^2$. Figure 13 shows results obtain using this approach. This figure highlights how the constitutive parameter value depends on the biaxiality ratio.

Such variations in the constitutive parameter, typically between uniaxial and equibiaxial loading conditions, are a well-known result since the pioneering work carried out by Treloar in 1944.[26] This clearly highlights the interest in gathering all the possible loading conditions in only one heterogeneous test for identifying constitutive parameters of hyperelastic models such as the neo-Hookean one. It should be noted that in the EQT zone, the constitutive parameter is noisier. This is due to the fact that the maximum level of the principal stretch reached is lower than for the other strain states and that it might be more affected therefore by nonentropic effects that are not taken into account in the model. It should be noted that the temperature variations during the deformation is lower than 1 K whatever the considered point in the field (see Figure 9). This confirms that in this range of temperature, $C_{nh}$ can be considered as temperature independent and its variations are therefore only due to the strain state.

# 6 | CONCLUSION

In this paper, a new inverse identification method is developed by reconstructing the heat source field from two approaches. The first one requires the measurement of the temperature field and the value of the thermophysical parameters. The second one requires the measurement of the kinematics field and the choice of a thermo-hyperelastic model that contains the parameters to be identified. In the present work, the method is applied to the challenging case of large deformation of rubber, where the movement of the points observed at the specimen surface has to be compensated. For that purpose, a motion compensation technique has been developed and validated. Finally, the identification has been performed at the local and global scales.

# ACKNOWLEDGEMENTS

The authors thank the National Center for Scientific Research (MRCT-CNRS and MI-CNRS), Rennes Metropole, and Region Bretagne for financially supporting this work and the PCM Technologies S.A.S company for providing the specimens. The authors also thank Dr. Mathieu Miroir, M Vincent Burgaud, and M Mickael Le Fur for having designed the biaxial tensile machine. Dr. Eric Robin is acknowledged for providing the mark tracking algorithm.

# ORCID

Jean-Benoît Le Cam https://orcid.org/0000-0002-2366-2512
REFERENCES

[1] S. Avril, M. Bonnet, A. S. Bretelle, M. Grédiac, F. Hild, P. Ienny, F. Latourte, D. Lemosse, S. Pagano, S. Pagnacco, F. Pierron, *Exp. Mech.* **2008**, *4*, 381.

[2] N. Promma, B. Raka, M. Grédiac, E. Toussaint, J.-B. Le Cam, X. Balandraud, F. Hild, *Int. J. Solids Struct.* **2009**, *46*, 698.

[3] M. Giton, A.-S. Caro-Bretelle, P. Ienny, *Strain* **2006**, *42*.

[4] J. S. Bergström, M. C. Boyce, *J. Mech. Phys. Solids* **1998**, *46*(5), 931.

[5] L. Mullins, *Rubber Chem. Technol.* **1949**, *22*(4), 1036.

[6] H. Bouasse, Z. Carrière, *Annales de la Faculté des Sciences de Toulouse* **1903**, *5*, 257.

[7] G. Chagnon, E. Verron, L. Gornet, G. Marckmann, P. Charrier, *J. Mech. Phys. Solids* **2004**, *52*, 1627.

[8] J. Diani, B. Fayolle, P. Gilormini, *Eur. Polym. J.* **2009**, *45*, 601.

[9] L. Mullins, *Rubber Chem. Technol.* **1948**, *21*, 281.

[10] L. Mullins, *Rubber Chem. Technol.* **1969**, *42*, 339.

[11] C. W. Bunn, *Proc. R. Soc. London, Ser. I* **1942**, *180*, 40.

[12] J. R. Katz, *Naturw* **1925**, *4*, 169.

[13] S. Toki, I. Sics, B. S. Hsiao, S. Murakami, M. Tosaka, S. Poompradub, S. Kohjiya, Y. Ikeda, *J. Polym. Sci. Part B: Polym. Phys.* **2004**, *42*, 956.

[14] S. Trabelsi, P.-A. Albouy, J. Rault, *Macromolecules* **2003**, *36*, 9093.

[15] A. Lachhab, E. Robin, J.-B. Le Cam, F. Mortier, Y. Tirel, F. Canevet, *Strain* **2018**, e12271.

[16] J.-B. Le Cam, *Polym. 2017*, 127, 166.

[17] M. T. Loukil, G. Corvec, E. Robin, M. Miroir, J.-B. Le Cam, P. Garnier, *Eur. Polym. J.* **2018**, *98*, 448.

[18] J. Diani, M. Brieu, P. Gilormini, *Int. J. Solids Struct.* **2006**, *43*, 3044.

[19] F. Laraba-Abbes, P. Ienny, R. Piques, *Polym. 2003*, *44*, 807.

[20] A. Robisson. **2000**, Comportement mécanique d'un élastomère chargé en silice. étude de l'influence des charges et modélisation par une loi visco-hyperélastique endommageable, Ph.D. Thesis.

[21] J.-B. Le Cam, *Rubber Chem. Technol.* **2010**, *83*, 247.

[22] L. R. G. Treloar, *The Physics of Rubber Elasticity*, Oxford University Press, Oxford **1975**.

[23] J. Bonet, R. D. Wood, *Nonlinear continuum mechanics for finite element analysis*, Cambridge University Press **1997**.

[24] G.-A. Holzapfel, *Non linear solid mechanics - a continuum approach for engineering Edited by* Wiley, Wiley **2000**.

[25] G. Marckmann, E. Verron, *Rubber Chem. Technol.* **2007**, *79*, 835.

[26] L. R. G. Treloar, *Trans. Faraday Soc.* **1944**, *40*, 59.

[27] C. G'Sell, A. Coupard, *Génie Mécanique des Caoutchoucs*, Appollor and INPL, Ecole des Mines de Nancy **1997**.

[28] M. Sasso, G. Palmieri, G. Chiappini, D. Amodio, *Polym. Test.* **2008**, *27*, 995.

[29] I. M. Ward, D. W. Hadley, *An Introduction to the Mechanical Properties of Solid Polymers*, Second, John Wiley and Sons Ltd, New-York **1993**.

[30] T. Guélon, E. Toussaint, J.-B. Le Cam, N. Promma, M. Grédiac, *Polym. Test.* **2009**, *28*, 715.

[31] M. Johlitz, S. Diebels, *Arch Appl. Mech.* **2011**, *81*, 1333.

[32] H. Seibert, T. Scheffer, S. Diebels, *Tech. Mech.* **2014**, *81*, 72.

[33] X. Balandraud, J. B. Le Cam, *Arch. Appl. Mech.* **2014**, *84*(6), 773.

[34] H. Louche. **2009**, Etudes de certains phénomènes de localisation à partir de champs thermomécaniques, Habilitation thesis, Savoie University.

[35] J. R. Samaca Martinez, E. Toussaint, X. Balandraud, J. B. Le Cam, D. Berghezan, *Mech. Mater.* **2015**, *81*, 62.

[36] J. R. Samaca Martinez, X. Balandraud, E. Toussaint, J. B. Le Cam, D. Berghezan, *Polym. 2014*, 55(24), 6345.

[37] B. Berthel, A. Chrysochoos, B.Wattrisse, A. Galtier, *Exp. Mech.* **2008**, *48*, 79.

[38] A. Chrysochoos, in *Colloque photomécanique*, vol. 95 **1995**, 201.

[39] A. Chrysochoos, H. Louche, *Int J Eng Sci 2000*, *38*, 1759.

[40] A. Chrysochoos, H. Louche, *Mat Sci Eng A-struct* **2001**, *307*, 15.

[41] X. Balandraud, A. Chrysochoos, S. Leclercq, R. Peyroux, *C R Acad Sci - Ser. IIB - Mech.* **2001**, *329*, 621.

[42] A. Chrysochoos, V. Huon, F. Jourdan, J.-M. Muracciole, R. Peyroux, B. Wattrisse, *Strain* **2010**, *46*, 117.

[43] F. Maquin, F. Pierron, *Mech. Mater.* **2009**, *41*, 928.

[44] J. B. Le Cam, J. R. Samaca Martinez, X. Balandraud, E. Toussaint, J. Caillard, *Exp. Mech.* **2015**, *55*(4, SI), 771.

[45] J. R. Samaca Martinez, J.-B. Le Cam, X. Balandraud, E. Toussaint, J. Caillard, *Polym. Test.* **2013**, *32*, 835-841.

[46] J. R. Samaca Martinez, J.-B. Le Cam, X. Balandraud, E. Toussaint, J. Caillard, *Polym. 2013*, 54, 2727-2736.

[47] J. R. Samaca Martinez, J.-B. Le Cam, X. Balandraud, E. Toussaint, J. Caillard, *Eur. Polym. J.* **2014**, *55*, 98.

[48] P. Jongchansitto, C. Douellou, I. Preechawuttipong, X. Balandraud, *Strain* **2019**, *55*(3), e12307.

[49] G. Marckmann. **2004**, Ph.D. Thesis, École Centrale de Nantes.

[50] H. Baaser, C. Hopmann, A. Schobel, *Arch. Appl. Mech.* **2013**, *83*, 273.

[51] E. Verron, A. Andriyana, *J. Mech. Phys. Solids* **2008**, *56*(2), 417. http://www.sciencedirect.com/science/article/pii/S0022509607001081

[52] P. Vacher, S. Dumoulin, F. Morestin, S. Mguil-Touchal, *Proc. IME. C. J. Mech. Eng. Part C: J. Mech. Eng. Sci.* **1999**, *213*.

[53] T. Pottier, M.-P. Moutrille, J.-B. Le Cam, X. Balandraud, M. Grédiac, *Exp. Mech.* **2009**, *49*, 561.

[54] A. Chrysochoos, B. Wattrisse, J.-M. Muracciole, Y. El Kaim, *J. Mech. Mater. Struct.* **2009**, *4*, 245.
[55] L. Bodelot, L. Sabatier, E. Charkaluk, P. Dufrenoy, *Mater. Sci. Eng. A-Structural Mater. prop. Microstruct.Process.* **2009**, *501*, 52.

[56] H. Louche, P. Schlosser, D. Favier, L Orgeas, *Exp. Mech.* **2012**, *52*, 1313.

[57] A. Tayeb, J.-B. Le Cam, M. Grédiac, E. Toussaint, E. Robin, X. Balandraud, F. Canévet, *Submitted Inter. J. Eng. Sci.* **2019**, *25*.