# A thermo-viscoelastic model for elastomeric behaviour and its numerical application

Adnane Boukamel, Stéphane Méo, Olivier Débordes, Marc Jaeger

► To cite this version:

Adnane Boukamel, Stéphane Méo, Olivier Débordes, Marc Jaeger. A thermo-viscoelastic model for elastomeric behaviour and its numerical application. Archive of Applied Mechanics, 2001, 71 (12), pp.785-801. 10.1007/s004190100191. hal-01236417

HAL Id: hal-01236417

https://univ-tours.hal.science/hal-01236417

Submitted on 3 Feb 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
# A thermo-viscoelastic model for elastomeric behaviour and its numerical application

A. Boukamel, S. Méo, O. Débordes, M. Jaeger

**Summary** This paper presents a model of thermo-mechanical behaviour of viscoelastic elastomers under large strain. A formulation is proposed with a generalisation to large strain of the Poynting-Thomson rheological model. A finite element formulation is then exposed taking the incompressibility constraint for mechanical equilibrium into account. On the thermo-mechanical coupling aspect, an algorithm of time discretisation is proposed with two time scales corresponding respectively to mechanical and thermal behaviours. Finally, an application for the simulation of a double-shearing test is presented with an analysis of parameters' influence and a comparison between numerical and experimental results.

**Keywords** Viscoelasticity, Large Strain, Thermomechanical Coupling, Elastomeric Behaviour

## 1

### Introduction

Elastomers are frequently employed in many industrial sectors such as the automobile and aeronautics industries. Their use has been extended to parts that undergo strong mechanical and thermal loadings. Since their mechanical properties highly depend on the temperature, the prediction of the behaviour and the assessment of the fatigue strength of elastomeric pieces, which are often closely linked to safety, require a local analysis based on a formulation of a thermo-mechanical model.

In literature, among the phenomena considered in order to analyse the local behaviour of elastomers under dynamic and thermal loadings, several approaches can be distinguished:

* hyperelastic modelling to characterise the static behaviour of the material;

* nonlinear viscoelasticity to simulate dissipative response and relaxation phenomena;

* micro-mechanical analysis for modelling the softening behaviour under deformation (Mullins' effect) or damage evolution;

* thermo-mechanical coupling to describe the changes of temperature due to the mechanical dissipation.

For hyperelasticity, several laws have been proposed which are based on the expression of strain energy density in isotropic, incompressible materials. Among these laws, statistical models were carried out with entropic considerations of the molecular chains configurations

A. Boukamel (✉)
Ecole Supérieure d'Ingénieurs de Marseille & Laboratoire de Mécanique et D'Acoustique, Technopôle de Château-Gombert,
13451 Marseille Cedex 20, France
e-mail: boukamel@esim.imt-mrs.fr

S. Méo, O. Débordes
Laboratoire de Mécanique et d'Acoustique,
CNRS, 31, Chemin Joseph Aiguier,
13402 Marseille Cedex 20, France

M. Jaeger
Institut Universitaire des Systèmes Thermiques Industriels,
Technopôle de Château-Gombert,
13453 Marseille Cedex 13, France
[1]. The other class of models reflects a phenomenological approach deduced from the isotropy and incompressibility hypothesis, which must be adjusted to experiments, [2-5].

The nonlinear viscoelastic aspects are treated in essentially two phenomenological approaches:

1. *The integral approach* using restrictive hypothesis such as the finite linear viscoelasticity (FLV) theory. This method is based on a modified linear Boltzmann superposition principle with a nonlinear strain measure. It uses an asymptotic expansion of relaxation functions. The separability principle between time and strain is expressed by the first term of this expansion, while the higher order terms constitute a correction, [6, 7]. More recently, several authors have proposed a modified FLV model introducing a generalised deformation measure, [8, 9]. The theoretical frame of this approach was essentially developed for nonlinear materials with fading memory and based on integral formulation of the stress tensor with respect of deformations history, [10-12].

2. *The differential approach* is based on the concept of intermediate state commonly used in describing finite elastic-plastic deformations, [13], and consists in a generalisation to large strains of classical rheological models, [14-16]. This theory introduces a viscoelastic internal variable. Its evolution is governed by constitutive laws deduced from thermodynamics of irreversible processes. This approach appears more adapted to numerical developments than the integral formulations.

Another important aspect of elastomeric behaviour is the softening effect, experimentally observed and called Mullins' effect, [17, 18]. This loss of stiffness in filled rubbers is micro-mechanically described by the rupturing of the bonds between the rubber matrix and embedded particles. A phenomenological model of damage and numerical applications are proposed in [19]. A statistical approach is used for modelling this phenomenon in [20], relating the damage mechanism to the maximum stretch of the deformation history. In [21], a transition to a phenomenogical model is proposed for numerical implementation. More recently, an application of continuum damage mechanics to Ogden-type material was proposed in [22]. The description of damage accumulation is handled similarly to that in [21] but another variable is added, which characterises the evolution of free energy during the process.

The dissipative nature of the viscoelastic behaviour involves an internal heat production, and so a temperature evolution takes place in the material. This thermo-viscoelastic coupling phenomenon was investigated recently in [23] in order to establish a general framework of nonlinear thermo-viscoelasticity. In the isothermal case, it is shown that a nonlinear generalisation of the Maxwell rheological model is obtained. Some guidelines were given to handle non-isothermal problems, taking the thermo-mechanical coupling into account.

Many works have contributed to the development of mathematical formulations and numerical methods, allowing precise simulations of hyperelastic or viscoelastic materials. In particular, some finite element (FE) formulations have been largely proposed in literature. They were adapted to large strain and took the incompressibility constraint into account, using the penalty method, [24, 25], or an augmented lagrangian method, [26], in order to solve the equilibrium problem in hyperelasticity. An equivalent procedure was also adopted for nonlinear viscoelasticity with a local treatment of internal variables, [15].

In this paper, we present a thermo-viscoelastic model for large deformation and finite variations of temperature. Our approach is fully phenomenological and based on the principles of irreversible thermodynamics. Micro-structural aspects, like interaction between the rubber components, stress softening phenomenon (Mullins' effect) or other damaging mechanism are not considered.

An FE implementation of this model was performed in [27] but restrained to problems with a uniform dissipation field and no temperature dependence of the mechanical behaviour. In the present work, we propose a numerical model including, on one hand, sequential and local coupling algorithm, and on the other hand, a temperature dependence of the mechanical properties. To our knowledge, no numerical developments for the case of large viscoelastic strains and large temperature changes have been proposed for this subject.

The physical problem is described firstly with regard to the basic principles of thermodynamics of continuous media adapted to the large deformation theory. A generalisation of a rheological model is then presented in order to obtain a law governing viscoelasticity, and finally the thermal equations are written in the Lagrangian configuration.

On the mechanical aspect, a perturbed Lagrangian formulation allows to handle the incompressibility constraint. For the thermal problem, a classical variational equation is given. These two variational forms are associated to FE approximations in order to obtain an explicit
coupling system which will be solved using a two-time scales algorithm according to the viscoelastic and thermal behaviour respectively.

Finally, this approach is validated by a simulation of a shearing test on a specimen composed of a two-layer test piece (elastomer/steel). A comparison between the numerical results and experimental measurements is finally offered and interpreted in terms of numerical and physical aspects.

## 2
Description of the physical problem

### 2.1
Description of the motion

Consider a body occupying the domain $\Omega$ in the undeformed reference configuration $C_0$ and the domain $\omega$ in the current configuration $C_t$. The motion of the body from $C_0$ to $C_t$ is locally described by the gradient tensor defined as

$$ \overline{\overline{F}}(\vec{X}, t) = \frac{\overline{\overline{\partial_x}}}{\partial \vec{X}}(\vec{X}, t), \qquad (1) $$

where $\vec{X}$ and $\vec{x}$ define, respectively, the position vector in the reference configuration and in the current one.

The dilatation tensor or right Cauchy-Green tensor is built from the gradient tensor, such that

$$ \overline{\overline{C}}(\vec{X}, t) = \overline{\overline{F}}^T \cdot \overline{\overline{F}}. \qquad (2) $$

To describe elastoplastic or viscoelastic behaviour in large deformations, an intermediate state $C_i$ is usually introduced which could be interpreted as the locally relaxed configuration, [28]. Thus, we consider the multiplicative decomposition

$$ \overline{\overline{F}} = \overline{\overline{F}}_e \cdot \overline{\overline{F}}_v, \qquad (3) $$

where $\overline{\overline{F}}_e$ defines a pseudo-gradient of the elastic motion, and $\overline{\overline{F}}_v$ a pseudo-gradient of the viscous motion. A viscous and an elastic Cauchy-Green tensors are associated

$$ \overline{\overline{C}}_v = \overline{\overline{F}}_v^T \cdot \overline{\overline{F}}_v, \quad \overline{\overline{C}}_e = \overline{\overline{F}}_e^T \cdot \overline{\overline{F}}_e. \qquad (4) $$

The total Cauchy-Green tensor is then given by

$$ \overline{\overline{C}} = \overline{\overline{F}}^T \cdot \overline{\overline{F}} = \overline{\overline{F}}_v^T \cdot \overline{\overline{C}}_e \cdot \overline{\overline{F}}_v. \qquad (5) $$

### 2.2
Thermodynamic formulation

The three conservation laws (mass, momentum and energy) of the classical thermodynamics, which are locally written in the Lagrangian description, govern the transformation of the system

$$ \rho \det(\overline{\overline{F}}) = \rho_0, \quad \rho_0 \ddot{u} = \operatorname{div} \overline{\overline{\Pi}} + \rho_0 \vec{f}, \quad \rho_0 \dot{e} = \overline{\overline{\Pi}} : \dot{\overline{\overline{F}}} - \operatorname{div} \overline{\overline{Q}} + \rho_0 r, \qquad (6) $$

where $\rho_0$ and $\rho$ are, respectively, the local mass density in $C_0$ and $C_t$, $\overline{\overline{\Pi}}$ is the first Piola-Kirchoff stress tensor, $\vec{f}$ corresponds to the specific (i.e. per mass unit) body force, $\overline{\overline{Q}}$ is the specific internal energy, $\overline{\overline{Q}}$ is the Piola-Kirchoff heat flux and $r$ is the rate of heat production per unit mass.

In order to be admissible, the thermodynamic process must satisfy the Clausius-Duhem inequality which states that the total dissipation is positive or equal zero. In the case of a Lagrangian description this inequality is written as

$$ \Phi = \Phi^m + \Phi^\theta \geq 0, \qquad (7) $$
where

$$
\begin{align}
\Phi^m &= \bar{\bar{\Pi}} : \overset{\bullet}{\bar{\bar{F}}} - \rho_0 (\dot{\Psi} + s \dot{T}) \quad \text{is the mechanical dissipation} \tag{8} \\
\Phi^\theta &= - \vec{Q} \cdot \frac{\bar{\bar{\nabla}}_L T}{T} \quad \text{is the thermal dissipation}
\end{align}
$$

and

$$
\Psi = e - Ts, \tag{9}
$$

$\Psi$ is the specific free energy, $s$ the specific entropy and $T$ the absolute temperature.

The first step in the development of the thermodynamic model consists in the choice of a set of state variables adapted to the problem. For the present study, we add to the usual state variables ($\bar{F}$, $T$), [28], another one associated to the viscous behaviour of material: the viscous dilatation tensor

$$
\bar{\bar{C}}_v = \bar{\bar{F}}_v^T \cdot \bar{\bar{F}}_v .
$$

So the rate of specific free energy can be written as

$$
\dot{\Psi} = \frac{\partial \Psi}{\partial \bar{\bar{F}}} : \overset{\bullet}{\bar{\bar{F}}} + \frac{\partial \Psi}{\partial \bar{\bar{C}}_v} : \overset{\bullet}{\bar{\bar{C}}}_v + \frac{\partial \Psi}{\partial T} \dot{T} . \quad (10)
$$

The mechanical dissipation has the following form:

$$
\Phi^m = \left( \bar{\bar{\Pi}} - \rho_0 \frac{\partial \Psi}{\partial \bar{\bar{F}}} \right) : \overset{\bullet}{\bar{\bar{F}}} - \rho_0 \frac{\partial \Psi}{\partial \bar{\bar{C}}_v} : \overset{\bullet}{\bar{\bar{C}}}_v - \rho_0 \left( s + \frac{\partial \Psi}{\partial T} \right) \dot{T} . \quad (11)
$$

On the other hand, the thermodynamics of irreversible processes states the existence of a
pseudo-potential of dissipation $\varphi$ such that

$$
\Phi^m = \frac{\partial \varphi}{\partial \bar{\bar{C}}_v} : \overset{\bullet}{\bar{\bar{C}}}_v . \qquad (12)
$$

3

**Constitutive equations**

3.1

**Mechanical behaviour**

The mechanical constitutive equations are determined by the extension to the case of large
deformation of the well-known Poynting–Thomson rheological model Fig. 1.

Let $\Psi_e$ and $\Psi_v$ define the specific free energies associated, respectively, to the instantaneous and delayed responses of the material. For these two quantities, the forms of strain energy corresponding, respectively, to the hyperelastic Hart–Smith model and the neoHookian model are adopted, [29]

$$
\rho_0 \Psi_e = \int c_1 \exp(c_3(I_1^e - 3)^2) dI_1^e + c_2 \ln \frac{I_2^e}{3}, \quad (13)
$$

$$
\rho_0 \Psi_v = a_1 (I_1^v - 3) . \tag{14}
$$

Here, $I_1^e, I_2^e$ denote the first and second principal invariants of $\bar{\bar{C}}_e$ and $I_1^v$ is the first one of $\bar{\bar{C}}_v$. The coefficients $c_1, c_2, c_3$ and $a_1$ are material characteristics that must be identified experimentally.

These strain energies are currently used to describe the mechanical behaviour of the hy-
perelastic incompressible media. According to the multiplicative decomposition (3) and the
rheological model, the incompressibility constraint applies both to observable variable $\bar{\bar{F}}$
and internal variable $\bar{\bar{C}}_v$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Poynting–Thomson model</figcaption></figure>

$$ \det(\bar{\bar{F}}) = 1, \quad \mathrm{tr}\left(\bar{\bar{C}}_v \cdot \bar{\bar{C}}_v^{-1}\right) = 0. \qquad (15) $$

The total free energy density of the system is then given by

$$ \rho_0 \Psi(\bar{\bar{F}}, \bar{\bar{C}}_v, T) = \begin{cases} \rho_0 (\Psi_e(\bar{\bar{F}}, T) + \Psi_v(\bar{\bar{C}}_v, T)) & \text{if } \det(\bar{\bar{F}}) = 1 \\ \cdot & \text{otherwise} \end{cases} \qquad (16) $$

The pseudo-potential of dissipation $\varphi$ takes the following quadratic form and only depends on $\bar{\bar{C}}_v$, [28]

$$ \varphi(\bar{\bar{C}}_v) = \begin{cases} \frac{1}{2}\eta(\bar{\bar{C}}_v : \bar{\bar{C}}_v) & \text{if } \mathrm{tr}(\bar{\bar{C}}_v \cdot \bar{\bar{C}}_v^{-1}) = 0 \\ \cdot & \text{otherwise} \end{cases} \qquad (17) $$

where $\eta$ is an experimentally identified viscosity parameter.

Substituting (16) and (17) in (11) and (12), and for processes, respectively, isothermally, reversible ($\dot{T} = 0$ and $\bar{\bar{C}}_v = 0$) and only isothermal, this choice of thermodynamic formulation leads to the following constitutive laws:

$$ \begin{aligned} \bar{\bar{\Pi}} &= 2\rho_0 \bar{\bar{F}}_e \cdot \left[ \left( \frac{\partial \Psi_e}{\partial I_1^e} + \frac{\partial \Psi_e}{\partial I_2^e} I_1^e \right) \bar{\bar{I}} - \frac{\partial \Psi_e}{\partial I_2^e} \bar{\bar{C}}_e \right] \cdot \bar{\bar{F}}_v^{-T} + p \det(\bar{\bar{F}}) \bar{\bar{F}}^{-T}, \\ -\rho_0 \bar{\bar{C}}_v^{-1} \cdot \bar{\bar{C}}_v \cdot \left[ \left( \frac{\partial \Psi_e}{\partial I_1^e} + \frac{\partial \Psi_e}{\partial I_2^e} I_1^e \right) \bar{\bar{I}} - \frac{\partial \Psi_e}{\partial I_2^e} \bar{\bar{C}}_v^{-1} \cdot \bar{\bar{C}}_v \right] \cdot \bar{\bar{C}}_v^{-1} + \rho_0 \frac{\partial \Psi_v}{\partial I_1^v} \bar{\bar{I}} + \eta \bar{\bar{C}}_v + q \bar{\bar{C}}_v^{-1} = 0, \\ \det(\bar{\bar{F}}) = 1, \quad \mathrm{tr}(\bar{\bar{C}}_v \cdot \bar{\bar{C}}_v^{-1}) = 0. \end{aligned} \qquad (18) $$

In these expressions, the incompressibility constraints (18c) are imposed through the Lagrange multipliers $p$ and $q$.

Left-multiplying equation (18b) by $\bar{\bar{C}}_v$ and considering only its deviatoric part, the Lagrange multiplier $q$ disappears and the complementary law becomes

$$ \begin{aligned} & -\rho_0 \left[ \bar{\bar{C}}_v \cdot \left( \left( \frac{\partial \Psi_e}{\partial I_1^e} + \frac{\partial \Psi_e}{\partial I_2^e} I_1^e \right) \bar{\bar{I}} - \frac{\partial \Psi_e}{\partial I_2^e} \bar{\bar{C}}_v^{-1} \cdot \bar{\bar{C}}_v \right) \cdot \bar{\bar{C}}_v^{-1} \right]^D \\ & + \rho_0 \left[ \frac{\partial \Psi_v}{\partial I_1^v} \bar{\bar{C}}_v \right]^D + \left[ \eta \bar{\bar{C}}_v \cdot \bar{\bar{C}}_v \right]^D = 0, \quad \mathrm{tr} \left( \bar{\bar{C}}_v \cdot \bar{\bar{C}}_v^{-1} \right) = 0, \end{aligned} \qquad (19) $$
which will be subsequently summarised by

$$ \dot{\overline{\overline{C}}}_v(t) = \overline{\overline{\chi}}(t, \overline{\overline{F}}(t), \overline{\overline{C}}_v(t)) . \qquad (20) $$

## 3.2

### Thermal behaviour

The Fourier law in the reference configuration gives the Piola–Kirchhoff heat flux

$$ \overline{\overline{Q}} = -\overline{\overline{K}}_L \cdot \overline{\overline{\nabla}}_L T(\overline{\overline{X}}, t) . \qquad (21) $$

In this relation, the conductivity tensor $\overline{\overline{K}}_L$ is expressed in the initial configuration and is given in [27]

$$ \overline{\overline{K}}_L = \overline{\overline{F}}^{-1} \cdot \overline{\overline{K}} \cdot \overline{\overline{F}}^{-T}, \qquad (22) $$

while $\overline{\overline{K}}$ is the usual Eulerian conductivity tensor. It is reduced, in the isotropic case, to $\overline{\overline{K}} = K \overline{\overline{I}}$, where K is the material conductivity factor.

# 4

## Thermo-mechanical model

### 4.1

#### Physical problem

As mentioned above, the state of the system is given by the three conservation laws of the classical thermodynamics.

Under the assumption of a quasi-static evolution, the mass and momentum conservation laws, added to the hypothesis of incompressibility medium, yield the following equations which govern the mechanical evolution:

$$ \frac{\rho}{\rho_0} = \det(\overline{\overline{F}}) = 1, \quad \vec{0} = \operatorname{div} \overline{\overline{\Pi}} + \rho_0 \vec{f} \quad \forall \vec{X} \in \Omega. \qquad (23) $$

This formulation has to be supplemented with the constitutive laws (18), (19).

The heat transfers are governed by the heat equation. It can be derived from the energy conservation law (6), by introducing the expression of the dissipation (7), (8) and the corresponding expression of $\overline{\overline{\Pi}}$ Eq. (18). We finally obtain, in Lagrangian description

$$ \rho_0 C_{\varepsilon} \dot{T} = -\operatorname{div} \overline{\overline{Q}} + \rho_0 r + D_S \quad \forall \vec{X} \in \Omega, \qquad (24) $$

where $C_{\varepsilon}$ is the heat capacity, $\overline{\overline{Q}}$ is the heat flux of Piola–Kirchhoff given by the Fourier law (22), $D_S$ is the internal source term per volume unit in the initial configuration, given by

$$ D_S = \eta \overline{\overline{C}}_v : \overline{\overline{C}}_v + T \left( \frac{\partial \overline{\overline{\Pi}}}{\partial T} : \overline{\overline{F}} - \frac{\partial \eta}{\partial T} \overline{\overline{C}}_v : \overline{\overline{C}}_v \right) . \qquad (25) $$

In this expression, the second term gives the power dissipation due to the variation of the mechanical properties with respect to the temperature, whereas the first term is the power dissipation due to the viscous behaviour of the material.

The local equations (23) to (25) are associated to the classical set of boundary conditions.

### 4.2

#### Variational formulation

The variational formulation of the quasi-static equilibrium problem is written under a perturbed Lagrangian form, [30]. Thus, the displacement $\vec{u}$ and the hydrostatic pressure p have to cancel the following integral form (where $\alpha$ is a positive perturbation coefficient):

$$ \mathcal{Q}^{\text{mec}}(\vec{u}, p) = \int_{\Omega} \overline{\overline{\Pi}} : \delta \overline{\overline{F}} dV - \int_{\partial \Omega_F} \overline{\overline{F}} \cdot \delta \vec{u} dS - \int_{\Omega} \rho_0 \vec{f} \cdot \delta \vec{u} dV + \int_{\Omega} \delta p (J - 1 - \alpha p) dV, \qquad (26) $$
for all the test functions $\delta\vec{u}$, $\delta p$.

The solution $(\vec{u}, p)$ must also locally satisfy the constitutive equations (18, 19) and the incompressibility constraint (15), when the parameter $\alpha$ tends to zero. However, a choice of $\alpha$ strictly positive but small, allows to obtain a solution satisfying approximately the incompressibility condition.

On the other hand, the temperature $T$ has to cancel the following integral form:

$$ \Im \text{ther}(T) = \int_{\Omega} \left( \rho_0 C_{\varepsilon} \dot{T} \delta T + \vec{\nabla}_L(\delta T) \cdot \overline{\vec{K}}_L \cdot \dot{\vec{\nabla}}_L(T) \right) dV - \int_{\partial \Omega_q} \vec{Q} \cdot \vec{N} \delta T \, dS \\ - \int_{\Omega} (D_s + \rho_0 r) \delta T \, dV, \qquad (27) $$

for all the test functions $\delta T$.

The FE approximation is classically introduced into the integral forms and yields a nonlinear system of equations. The unknowns are the nodal displacements and the nodal temperatures. Then, the mechanical and thermal problems are solved separately.

## 4.3

### Finite element approximation

#### Mechanical problem

For the mechanical problem, a classical FE is built using the commonly called T6/P1 triangular element, with a quadratic interpolation for displacement and a constant pressure. This element satisfies the discrete LBB condition and yields stable pressure approximations, [31].

The hydrostatic pressure may to be discontinuous between the elements. It is considered as an internal degree of freedom and thus eliminated by a static condensation at element level.

#### Thermal problem

Using the same geometric discretization, the thermal problem approximation is based on a linear interpolation for temperature. The global system so obtained is then solved by an explicit Eulerian time-integration scheme.

## 4.4

### Local solving of complementary law

For the mechanical problem, the construction of the elementary matrices and vectors requires the evaluation, at each integration point, of the internal variable $\overline{C}_v$ and the derivative term $\partial\overline{C}_v/\partial\overline{F}$. They are given by

$$ \overline{\overline{C}}_v(t) = \overline{\overline{\chi}}_v t, \overline{\overline{F}}(t), \overline{\overline{C}}_v(t) \quad \text{on } [t_0, t_0 + \Delta t], \quad \overline{\overline{C}}_v(t_0) = \overline{\overline{C}}_v^0, \qquad (28) $$

and

$$ \frac{\partial \overline{\overline{C}}_v}{\partial \overline{\overline{F}}} = \frac{\partial \overline{\overline{\chi}}_v}{\partial \overline{\overline{C}}_v} : \frac{\partial \overline{\overline{C}}_v}{\partial \overline{\overline{F}}} + \frac{\partial \overline{\overline{\chi}}_v}{\partial \overline{\overline{F}}} \quad \text{on } [t_0, t_0 + \Delta t], \quad \left. \frac{\partial \overline{\overline{C}}_v}{\partial \overline{\overline{F}}} \right|_{t=t_0} = 0. \qquad (29) $$

By means of an implicit Eulerian scheme, [32], the system of differential equations (28) gives on $[t_n, t_n + \delta t]$ a sub-interval of $[t_0, t_0 + \Delta t]$

$$ \overline{\overline{C}}_v(t_{n+1}) = \overline{\overline{C}}_v(t_n) + \delta t \overline{\overline{\chi}}_v (t_{n+1}, \overline{\overline{F}}(t_{n+1}), \overline{\overline{C}}_v(t_{n+1})). \qquad (30) $$

Then, with a Newton-Raphson scheme, the obtained nonlinear system is solved.

The second system (29) is linear, so the solution could be approached by a classical Crank-Nicholson scheme, with the same time step as used in (28), because the approximate solution of this system is taken into account to solve (29).
## 4.5

### Coupling algorithm

Since the time scales of the mechanical and thermal problems are very different, they can be solved separately, [33]. The coupling between the two problems Fig. 2 is then taken into account through the algorithm shown in Fig. 3.

The mechanical problem needs smaller time-increments than the thermal one. In other words, it is possible to realise many mechanical increments without an actualisation of the mechanical characteristics according to the temperature.

Moreover, during a mechanical step (less than one second), the evolution of the temperature remains imperceptible. Therefore, in the numerical implementation, the coupling terms of (25),

$$T \frac{\partial \overline{\Pi}}{\partial T} : \overline{\ddot{F}} \quad \text{and} \quad T \frac{\partial \eta}{\partial T} : \overline{\ddot{C}_v} : \overline{\ddot{C}_v},$$

can be neglected. Thus, the source term $D_s$ is reduced to the mechanical dissipation $\Phi^m$. The mechanical properties are updated according to the temperature obtained after a thermal cycle. This means that the thermo-mechanical coupling algorithm comes to a fixed-point method on temperature, as illustrated Fig. 2.

## 5

### Validation and applications

#### 5.1

##### Test description

In order to check if the proposed model can describe the thermo-mechanical behaviour of elastomers, a two-layer elastomer-steel test piece has been instrumented, Fig. 4. The elastomer part is made of dimethyl-vinyl-siloxan vulcanised by peroxide.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Thermo-mechanical coupling</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Coupling algorithm</figcaption></figure>
The test piece is put into a thermal enclosure at a temperature of 27 °C, and subjected to cyclic shearing deformation

$$ \gamma(t) = \Gamma \sin(2\pi ft) , $$

with a frequency $f = 4.5$ Hz. The amplitude $\Gamma = 0.5$ defines the relative shearing displacement applied to the external armatures, while the central armature is fixed.

Two thermocouples, Fig. 5, measure the temperature on one of the elastomeric layers.

## 5.2

### Material characteristics

The mechanical properties of the steel armatures are taken as classical values: Young modulus $E = 210000$ MPa, Poisson ratio $\nu = 0.3$. For the elastomeric material, the parameters of the Poynting-Thomson model are assumed to be temperature-dependent. These parameters are determined from experimental measurements, like relaxation under axial traction and cyclic shearing tests, using the following identification procedure:

(i) By identifying the average curve of a cyclic-shearing response, see Fig. 6a, to a shearing hyperelastic curve (Hart-Smith model), a first set of values ($c_1, c_2, c_3$) is estimated.

(ii) Then, through an adjustment to the relaxation test and by considering ($c_1, c_2, c_3$) as given, $a_1$ and $\eta$ are determined.

(iii) Finally, a correction and the complete determination of all parameters is provided according to a least-square fitting on the hysteresis curve of the shearing test. To do so, the curve is decomposed into an average curve and a difference one, obtained by the stress difference at constant strain, see Fig. 6b.

These three steps are realised for different temperatures and followed by an Arrhenius-law interpolation, [34], in order to obtain the evolution according to temperature as shown in Fig. 7. The thermal properties of the materials are given in Table 1, [27].

## 5.3

### Validation

#### Numerical simulation

In this section, the first series of computations is performed to determine the mesh-refinement influence. Then, the results are analysed according to the computation time of each problem, mechanical and thermal.

Due to the symmetry of the problem and the plane-strain hypothesis, the computation of the shearing test has been carried out on a half cross section of the test piece.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Two-layer test piece</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Half cross section of the test piece. Position of the thermocouples</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Identification procedure</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Elastomer parameters in dependence of the temperature</figcaption></figure>

For the thermal problem (see Fig. 5 and Table 3), linear triangular, bilinear quadrilateral and linear flux elements are used to mesh, respectively, the elastomeric layer, the metallic armatures and the boundaries. A zero-flux condition is applied on the symmetry axis, and a linear convective heat-transfer condition is imposed on the other boundaries.

For the mechanical problem using the plane strain hypothesis, the mesh is the same as the thermal one (without the linear flux elements), but a quadratic interpolation is used. The boundary conditions are given in Table 3.

## Mesh influence

Three spatial discretizations are selected (Fig. 8. and Table 2). The characteristics of each mesh are given in the Table 2. Each computation were decomposed into three times (mechanical cycle $\Delta t = 0.89$ s followed by a thermal cycle $\Delta t = 200$ s.

Some mechanical and thermal results are presented for the purpose of analysing the mesh refinement influence, Fig. 9. On Fig. 9a, the quantity $(p - p_0)/p_0$ is evaluated, where p is the
Lagrange multiplier in the center of the elastomeric layer and $p_0$ is its analytical value in the case of a uniform shearing. We see that for the mesh B, we have an error of only 3%. Concerning the amplitude of the reaction force, the meshes B and C show a difference of less than 3%, Fig. 9b.

The ratio of the mechanical dissipation and the external power gives also less than 1% of variation, Fig. 9c. Finally, in Fig. 9d, the temperatures in the center (thermocouple 1, Fig. 5) of

<table><caption>**Table 1.** Thermo-physical properties of the materials</caption><thead><tr><th></th><th>Elastomer</th><th>Steel</th></tr></thead><tbody><tr><td>Conductivity factor K (W·m<sup>-1</sup>·K<sup>-1</sup>)</td><td>0.127</td><td>45.0</td></tr><tr><td>Specific heat capacity ρ<sub>0</sub>C<sub>e</sub> (J·m<sup>-3</sup>·K<sup>-1</sup>)</td><td>0.74E+06</td><td>3.5E+06</td></tr><tr><td>Transfer coefficient h (W·m<sup>-2</sup>·K<sup>-1</sup>)</td><td>17.0</td><td>30.0</td></tr></tbody></table>

<table><caption>**Table 2.** Numbers of elements and degrees of freedom according to the problem and the refinement</caption><thead><tr><th>Problem</th><th>Refinement</th><th>A</th><th>B</th><th>C</th></tr></thead><tbody><tr><td rowspan="3">Mechanical problem</td><td>Elements</td><td>130</td><td>400</td><td>1420</td></tr><tr><td>d.o.f. (Disp.)</td><td>742</td><td>1938</td><td>6250</td></tr><tr><td>d.o.f. (Pres.)</td><td>80</td><td>320</td><td>1280</td></tr><tr><td rowspan="2">Thermal problem</td><td>Elements</td><td>171</td><td>459</td><td>1515</td></tr><tr><td>d.o.f. (Temp.)</td><td>121</td><td>285</td><td>853</td></tr></tbody></table>

<table><caption>**Table 3.** Boundary conditions</caption><thead><tr><th>Segments</th><th>Mechanical boundary conditions</th><th>Thermal boundary conditions</th></tr></thead><tbody><tr><td>AB</td><td>Stress free</td><td>Convective flux</td></tr><tr><td>BC</td><td>u<sub>x</sub> = U, u<sub>y</sub> = 0</td><td>Imposed temperature</td></tr><tr><td>CD</td><td>Stress free</td><td>Convective flux</td></tr><tr><td>DE</td><td>Stress free</td><td>Convective flux</td></tr><tr><td>EF</td><td>u<sub>x</sub> = 0, u<sub>y</sub> = 0</td><td>Imposed temperature</td></tr><tr><td>FA</td><td>u<sub>y</sub> = 0 (symmetry)</td><td>Zero flux (symmetry)</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>coarse mesh (A)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>medium mesh (B)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>accurate mesh (C)</figcaption></figure>

<caption>**Fig. 8.** Mesh refinement</caption>
the elastomeric layer are compared at the time 50 s for the three meshes. The results for the models B and C are quite the same.

These investigations point out that the more accurate mesh C gives the best results; but it also demands more CPU time. A good compromise between the results validity and the computational time seems to be the mesh B. For these reasons, all the next simulations presented in this paper have been done with the mesh B.

## Time influence

In order to analyse the influence of the time partition for the sequential thermo-mechanical coupling, two strategies are checked.

Strategy (i): three times a mechanical cycle $\Delta t = 0.89$ s followed by a thermal cycle $\Delta t = 200$ s;

Strategy (ii): four times a mechanical cycle $\Delta t = 0.89$ s followed by a thermal cycle $\Delta t = 50$ s and

two times a mechanical cycle $\Delta t = 0.89$ s followed by a thermal cycle $\Delta t = 200$ s.

Figure 10 presents the evolution of the relative difference between the temperature in the piece centre obtained with these two strategies. The curve grows to a relative value of about 7% and then decreases to a stable value of 3% when the temperature tends to the stationary solution.

<figure><img src="image_12.png" /><figcaption>Fig. 9. Mesh influence</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Relative difference of the temperature in the center of the parts according to the time partition</figcaption></figure>
The time discretization (ii) seems to be more consistent with the important variation of the temperature at the beginning, since it takes better into account the quick variation of the mechanical characteristics. In the following simulations, the configuration (ii) will be adopted.

## 5.4
### Numerical results and comparison with experiment

#### Mechanical results

Figure 11 shows the map of the time-average of the mechanical dissipation power on the deformed shape of the specimen, and this for a few cycles. It can be noticed that the maximums

<figure><img src="image_14.png" /><figcaption>Fig. 11. Map of the mean dissipation over time power (in MW m⁻³)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Reaction force in dependence of the displacement (stabilised cycles)</figcaption></figure>
occur in the corners of the elastomer bloc. The maps of the last two cycles are nearly the same. The same phenomenon is established in the evolution of the reaction force with the displacement, Fig. 12. The three curves indicate a loss of the stiffness. It can be concluded that the spatial and the time stabilisations are reached after the fourth cycle. Figure 13 presents the evolution over time of the spatial average of the dissipated power (continuous line) and its time average (dashed line) during few first mechanical cycles. It reaches a stabilised value after a few mechanical periods. This seems to prove that it is sufficient to consider only three periods for a mechanical cycle.

## Thermal results and comparison with measurements

It can be noticed that the temperature initially increases near the maximum dissipation zone (in the corners), and then spreads toward the center of the elastomer, Fig. 14.

The curves displayed in Figs. 15, 16 enable us to make two remarks. If we consider that the entire mechanical dissipation is transformed into heat, we can see a wide difference between the computed results and the experimental results on Fig. 17.

<figure><img src="image_16.png" /><figcaption>Fig. 13. Variation of the dissipation power (continuous line) and its mean over time (dashed line)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Evolution of the temperature field over time</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 15. Comparison of the computed temperature with the measured temperature in the center of the elastomer</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. Comparison of the computed temperature with the measured temperature after a correction of the rate of transmitted power in the center and the corner of the elastomer</figcaption></figure>

To fit the computed results with the experimental ones, we must assume that only a part of the mechanical dissipation is transformed into heat. In the Fig. 16, the evolution of the numerical temperature in the center and in the corner (thermocouple 3, Fig. 5) of the elastomer layer is plotted, with only 70% of mechanical dissipation taken in the thermal computation into account. We state that the two curves have a good agreement with the experiment.

This result was already invoked in previous works concerning the thermo-mechanical behaviour of other classes of materials, like metallic material, [35], and more recently, in [36], for filled rubber vulcanizates like styrene-butadiene rubber (SBR) and natural rubber. Thus, we can suggest that the rest of the energy is used by the material to reorganise its molecular structure. This possibility requires more investigations and provides a perspective for future works.

# 6

## Conclusion

A new model has been developed to deal with mechanical and thermal problems, allowing for an interaction of these two phenomena.

This explicit coupling has a solving algorithm adaptable to the two distinct time-scales. Moreover, the developments of these two models are totally independent. As an immediate consequence, they can be easily and separately upgraded. For the mechanical model, a hyperviscoelastic behaviour is proposed, initially based on the Hart–Smith hyperelastic law according to the incompressibility constraint and on a simple form of the pseudo-potential of dissipation. The thermal formulation is carried out under finite transformation hypothesis in the lagrangian configuration in order to fit it to the mechanical formulation. The thermo-
viscoelastic model has allowed the simulation of the shearing test with a specimen subjected to cyclic loading. The simulation gives an estimation of the dissipation field and consequently the evolution of temperature in the elastomeric layer.

The examination of the dissipation field allows to notice in the initial steps the prominence of two high dissipation zones. They are located at the corners of the elastomeric layer where the shear is the more important. Over time, the field tends to become uniform.

The computed temperature field is similar to the experimental one, but the internal source term seems to be overestimated. This could be explained by the fact that a part of this source is used for internal restructuring.

This phenomenon could provide an object for future research. A formulation including a new internal variable, could help to characterise the microscopic restructuring.

## References

1. Treolar, L.R.G.: The elasticity of a network of long chain molecules. Trans Faraday Soc 39 (1943) 36-64

2. Mooney, M.: A theory of large elastic deformation. J Appl Phys 11 (1940) 582-592

3. Rivlin, R.S.: Some topics in finite elasticity. 1st Symp Naval Struct Mech, 11-14 August, 1958, Stanford

4. Hart-Smith, L.J.: Elasticity parameters for finite deformations of rubber-like materials. J Appl Math Phys 17 (1966) 608-626

5. Ogden, R.W.: Large deformation isotropic elasticity, on the correlation of theory and experiment for incompressible rubber-like solids. Proc Roy Soc A326 (1972) 565-584

6. Coleman, B.D.; Noll, W.: Foundations of linear viscoelasticity. Rev Mod Phys 33 (1961) 239-249

7. Lianis, G.: Constitutive equations of viscoelastic solids under finite deformation. Purdue University Report, A&ES (1963) 63-11

8. Chang, W.V.; Bloch, R.; Tschoegl, N.W.: The behaviour of rubber-like materials in moderately large deformations. J Rheol 22 (1978) 1-32

9. Morman, K.N.: Original contributions. An adaptation of finite linear viscoelasticity theory for rubber-like viscoelasticity by use of the generalised strain measure. Rheol Acta 27 (1988) 3-14

10. Green, A.E.; Rivlin, R.S.: The mechanics of non-linear materials with memory. Arch Rational Mech Anal 1 (1957) 1-21

11. Coleman, B.D.: Thermodynamics of materials with memory. Arch Rational Mech Anal 17 (1964) 1-46

12. Christensen, R.M.: Theory of Viscoelasticity, an Introduction. New York: Academic Press 1971

13. Sidoroff, F.: Variables internes en viscoélasticité, II. Milieux avec configuration intermédiaire. J Méc 14 (1975) 571-595

14. Le Tallec, P.: Numerical Analysis of Viscoelastic Problems. Paris: Masson 1990

15. Le Tallec, P.; Rahier, C.: Numerical models of steady rolling for non-linear viscoelastic structures in finite deformations. Int J Numer Methods Eng 37 (1994) 1159-1186

16. Leonov, A.I.: On thermodynamics and stability of general Maxwell-like viscoelastic constitutive equations. Theor Appl Rheol. Proc XIth Int Cong Rheol, Brussels, 1992

17. Harwood, J.A.C.; Mullins, L.; Payne, A.R.: Stress softening in rubbers: a review. J IRI (1967) 17-27

18. Mullins, L.; Tobin, N.R.: Stress softening in rubber vulcanizates. Part I: use of a strain amplification factor to describe the elastic behavior of filler-reinforced vulcanized rubber. J Appl Polymer Sci 9 (1965) 2933-3009

19. Simo, J.C.: On a fully three-dimensional finite-strain viscoelastic damage model: formulation and computational aspects. Comput Methods Appl Mech Eng 60 (1987) 153-173

20. Govindjee, S.; Simo, J.: A micro-mechanically based continuum damage model for carbon black-filled rubbers incorporating Mullins' effect. J Mech Phys Solids 39 (1991) 87-112

21. Govindjee, S.; Simo, J.: Transition from micro-mechanics to computationally efficient phenomenology: carbon black filled rubbers incorporating Mullins' effect. J Mech Phys Solids 40 (1992) 213-233

22. Miehe, C.: Discontinuous and continuous damage evolution in Ogden-type large-strain elastic materials. Eur J Mech A/Solids 14 (1995) 697-720

23. Holzapfel, G.A.; Simo, J.: A new viscoelastic constitutive model for continuous media at finite thermomechanical changes. Int J Solids Struct 33 (1996) 3019-3034

24. Malkus, D.S.: Finite element with penalties in nonlinear elasticity. Int J Num Meth Eng 16 (1980) 121-136

25. Oden, J.T.: A theory of penalty methods for finite element approximations of highly nonlinear problems in continuum mechanics. Comput Struct 8 (1978) 445-449

26. Glowinski, R.; Le Tallec, P.: Numerical solution of problems in incompressible finite elasticity by augmented Lagrangian methods; II. three-dimensional problems. SIAM J Appl Math 44 (1984) 710-733

27. Bérardi, G.; Jaeger, M.; Martin, R.; Carpentier, C.: Modelling of a thermo-viscoelastic coupling for large deformations through finite element analysis. Int J Heat Mass Transfer 39 (1996) 3911-3924

28. Sidoroff, F.: Un modèle viscoélastique non linéaire avec configuration intermédiaire. J Méc 13 (1974) 679-713

29. Oden, J.T.: Finite Element on Nonlinear Continua. New York: McGraw-Hill 1972

30. Carey, G.F.; Oden, J.T.: Finite Elements, vol 2, pp 96-166. New Jersey: Prentice-Hall 1986
31. Oden, J.T.; Kikuchi, N.: Finite element methods for constrained problems in elasticity. Int J Num Meth Eng 18 (1982) 701-725

32. Bonelli, S.; Golay, F.; Débordes, O.: Résolution autoadaptative par éléments finis de problèmes de diffusion fortement non-linéaires. Rev Eur Eléments Finis 2 (1993) 495-515

33. Bérardi, G.: Modélisation numérique du comportement thermo-viscoélastique d'un élastomère en grandes déformations. Doctorat de l'Université de Provence 1995

34. Carpentier-Gabrieli, C.: Modélisation théorique et numérique du comportement viscoélastique d'élastomères sous sollicitations harmoniques. Doctorat de l'Université de la Méditerranée 1995

35. Chysochoos, A.: Energy balance for elastic plastic deformation at finite strain. J Theo Appl Mech 4 (1985) 589-614

36. Kar, K.K.; Bhowmick, A.K.: Hysteresis loss in filled rubber vulcanizates and its relationship with heat generation. J Appl Polym Sci 64 (1997) 1541-1555