# On the nonlinear viscoelastic behavior of rubber-like materials: Constitutive description and identification

Adel Tayeb, Makrem Arfaoui, Abdelmalek Zine, Adel Hamdi, Jalel

Benabdallah, Mohamed Ichchou

► To cite this version:

Adel Tayeb, Makrem Arfaoui, Abdelmalek Zine, Adel Hamdi, Jalel Benabdallah, et al. On the nonlinear viscoelastic behavior of rubber-like materials: Constitutive description and identification. International Journal of Mechanical Sciences, 2017, 130, pp.437-447. 10.1016/j.ijmecsci.2017.06.032. hal-02444061

HAL Id: hal-02444061

https://hal.science/hal-02444061

Submitted on 17 Jan 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
On the nonlinear viscoelastic behavior of rubber-like materials: Constitutive description and identification

Adel Tayeb<sup>a,b</sup>, Makrem Arfaoui<sup>a</sup>, Abdelmalek Zine<sup>c</sup>, Adel Hamdi<sup>a</sup>, Jalel Benabdallah<sup>a</sup>, Mohamed Ichchou<sup>b</sup>

<sup>a</sup>Université de Tunis El Manar, École Nationale d'Ingénieurs de Tunis, LR-11-ES19
Laboratoire de Mécanique Appliquée et Ingénierie, Tunis 1002, Tunisie

<sup>b</sup>Laboratoire de tribologie et dynamique des systèmes, École Centrale de Lyon, Ecully
69130, France

<sup>c</sup>Université de Lyon, Institut Camille Jordan, CNRSUMR5208, Département de
Mathématiques et Informatique, École Centrale de Lyon, 36 av. Guy de Collongue,
69134 Ecully Cedex, France

**Abstract**

The main concern of this paper is the development of a fully three dimen-
sional viscoelastic model at finite strain to describe nonfactorizable behavior
of rubber-like materials within the framework of rational thermodynamics
and internal state variable approach such that the second law of thermody-
namics in the form of Clausius-Duhem inequality is satisfied. The nonfactor-
izable aspect of the behavior is introduced via a strain dependent relaxation
times. The model is applied to describe the response of the isotropic Pip-
kin multi-integral viscoelastic model and the BIIR material, several parame-
ters involved are then identified using quasi-static and dynamic experiments
thanks to a least-square minimization procedure performed with Matlab soft-
ware. The proposed model is able to reproduce quasi-static experimental
response and show an ability to predict the dynamic response of nonfactor-
izable rubber-like materials in a wide range of strain.

***Keywords:*** nonfactorizable viscoelasticity, Rubber, Finite-strain, Identification

# 1. Introduction

It is well known that rubber-like materials exhibit nonlinear viscoelastic
behavior over a wide range of strain and strain rates confronted in several
engineering applications such as civil engineering, automotive and aerospace industries. This is due to their capacity to undergo high strain and strain rates without exceeding the elastic range of behavior. Further, the time dependent properties of these materials, such as shear relaxation modulus and creep compliance, are, in general, functions of the history of the strain or the stress [1]. Therefore, in a wide range of strain, a linear viscoelasticity theory is no longer applicable for such material and new constitutive equations are required to fully depict the behavior of rubber-like materials for quasi-static and dynamic configurations of huge interest in engineering applications.

In the literature, several phenomenological models have been developed to describe the nonfactorizable behavior of rubber-like materials, namely the Solid-Liquid viscoelastic model in the series of papers by Chang et al. [2], [3] and Sullivan in [4] for which a generalized measure of deformation has replaced the strain tensor in the linear Boltzmann convolution integral model and the nonlinear viscoelastic model by Schapery [5] in which the creep compliance and the shear relaxation functions were considered stress-dependent and strain-dependent functions respectively and the model of Valanis [6] in which a total thermodynamic formulation led to a constitutive equation depending on the deformation via a deformation shift function in analogy with the so-called rheologically simple materials.

In the other hand, other models based on the microstructure of the polymeric chain have been proposed such as the model by Knauss et al. [7] in which, following polymer science, time dependent functions were dependent upon the volumetric strain via a strain shift function and the model by Caruthers et al. [8] in which the strain shift function was expressed in terms of the configurational energy of the molecular structure.

In this work we shall develop a nonlinear viscoelastic model at finite strain within the framework of rational thermodynamics and the approach of internal state variables, the model is derived through a combination of approaches in [9], [10] and [11] taking into account the dependence of the time dependent functions upon the state of the strain. The model's parameters are then identified using data generated from the multi-integral viscoelastic model of Pipkin [12] and experimental data for bromobutyl (BIIR) from [13].

This paper is organized as follows: in section 2, a one dimensional nonlinear viscoelastic model is developed using a modified Maxwell rheological model. In the section 3, this model is extended to the fully nonlinear formulation using a nonlinear set of evolution equation of the internal state variables within the rational thermodynamic framework. The shear relaxation modulus is set
to be a function upon the invariants of the right Cauchy-Green strain tensor
via a strain shift function analogous to the temperature shift function for the
thermorheologically simple materials, this choice is motivated experimentally
[13] and rheologically [14]. The constitutive equation for the stress is then
obtained by resolving the set of nonlinear evolution equations. In section 4,
a systematic identification procedure of several parameters involved in the
model is highlighted. The optimization problems arising from this identifica-
tion procedure are solved by a modified least square minimization algorithm
with Matlab software. The section 5 is devoted to the results of this identifi-
cation procedure using an experimental characterization of the Bromobutyl
BIIR [13] and a theoretical data using the Pipkin model [12]. The capacity
of the model to describe the behavior of the material is then outlined.

## 2. Experimental and rheological motivations

In this **section** we develop the rheological and experimental arguments
leading to the proposed finite strain viscoelastic model. To motivate the three
dimensional model developed below, we first highlight some experimental
results leading to this model and then we consider a suitable modification
to the generalized Maxwell rheological model to build the one dimensional
nonlinear viscoelastic model.

### 2.1. *Experimental motivation*

A significant class of rubbers shows nonfactorizable behavior at low and
average range of strain. This phenomenon consists on the dependence of
the shear relaxation modulus upon strain level. Several works were dedi-
cated to deal with this class of behavior especially the series of papers by
[4] and [15]. In a recent work [13], an experimental characterization was
carried out with three rubber-like materials: the natural rubber (NR), the
Bromobutyl (BIIR) and a mixture of these materials (NR-BIIR). These ma-
terials showed a dependence of the shear relaxation modulus upon strain.
In figure 1 it is plotted the logarithm of the shear relaxation module G(t)
versus the logarithm of time for two different level of strain 10% and 50% for
BIIR material. The shear relaxation modulus shows a dependence upon the
strain level which leads according to [16] to a shift in the time with a strain
dependent function since the shear relaxation module at any level could be
obtained through a combination of a vertical and horizontal translation from
the reference curve at a strain level of 10%. Therefore a one dimensional
77 viscoelastic model, taking in consideration these results, is developed in the
78 next section through a generalization of the Maxwell rheological model.

<figure><img src="image_2.png" /><figcaption>Figure 1: Dependence of the shear relaxation modulus upon strain for BIIR rubber</figcaption></figure>

## 2.2. Rheological motivation

Before we develop the three-dimensional viscoelastic model, we shall investigate the following formulation for a standard linear solid of the figure 2. In this model, $\sigma$ denotes the total stress, $\epsilon$ denotes the total strain, $k_i$ and $\tau_i$ are the parameters of the Maxwell model. Unlike the rheological model used in [17], the relaxation times $\tau_i$ are, due to the experimental result outlined above, functions of the total strain $\epsilon$. Furthermore, the stress in the spring of each Maxwell branch is denoted by $q_i$ and its governed by the following evolution equation.

$$ \dot{q}_i + \frac{1}{\tau_i(\epsilon)} q_i = \frac{1}{\tau_i(\epsilon)} k_i \epsilon, \quad q_i|_{t=0} = 0 \qquad (1) $$

The total stress $\sigma$ derive directly from the rheological model of figure 2 as the difference between the elastic equilibrium stress and the non-equilibrium stresses $q_i$.

$$ \sigma = k_i \epsilon - \sum_i q_i \qquad (2) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Generalized Maxwell model</figcaption></figure>

The time parameters of the Maxwell model are set to be a strain dependent function; this idea follows from the description of thermorheologically simple materials behavior see [16] and [18], for which all parameters are temperature dependent via a single variable function called temperature shift-function. [5] [19] and [15] among others generalized this notion to describe thermorheologically complex materials behavior where the shift function depend upon temperature and stress or strain. Hence, the time parameters take the following form.

$$ \tau_i(\varepsilon) = a(\varepsilon) \tau_i \qquad (3) $$

a(ε) is a non-negative strain function called strain shift function. Therefore, the law of evolution of the equation 1 became a linear differential equation over the reduced time, after considering the form of the time parameters of
102 equation 3.

$$ \frac{dq_i}{d\xi} + \frac{1}{\tau_i}q_i = \frac{1}{\tau_i}k_i \epsilon \quad \text{with} \quad \xi(t) = \int_0^t \frac{dt'}{a(\epsilon)} \qquad (4) $$

103 $\xi(t)$ is an increasing function of time, integration of 4 and substitution of $q_i$ by its expression into 2 yield the expression of the total stress $\sigma$ as a Boltzmann convolution integral of the strain.
104
105

## 3. Fully nonlinear viscoelastic model

107 In this section we extend the formulation outlined above to the fully nonlinear range. Hence, the mechanical framework and the thermodynamic assumptions leading to the model are outlined. It should be noted, however, that this model is derived through an isothermal conditions.
108
109
110

### 3.1. Mechanical framework and form of the Helmholtz free energy density

111 Consider a viscoelastic material with reference placement $\Omega_0$ in the reference configuration $C_0$. It occupies at the time $t$ the placement $\Omega$ in the deformed configuration $C_t$ see figure 3. Let $\phi$ denote a macroscopic motion between the two configurations, which maps any point $X$ in the reference configuration $C_0$ to the point $x$ in the deformed configuration. Let $F(X,t) = \partial x/\partial X$ be the deformation gradient tensor. Likewise, let $J = \det(F)$ be the jacobian of the deformation gradient tensor. From the deformation gradient $F(X,t)$ the deformation tensor of Green Lagrange $E = 1/2(C - I)$, the right and left Cauchy-Green strain tensors $C = F^t F$ and $B = FF^t$ are obtained, together with their principal invariants.
112
113
114
115
116
117
118
119
120
121

$$ I_1 = \mathrm{tr}C, \quad I_2 = \frac{1}{2}[(\mathrm{tr}C)^2 - \mathrm{tr}C^2] \quad \text{and} \quad I_3 = \det(C) = J^2 \qquad (5) $$

122 which, otherwise, can be expressed in terms of principal stretches by

$$ I_1 = \lambda_1^2 + \lambda_2^2 + \lambda_3^2, \quad I_2 = \lambda_1^2\lambda_2^2 + \lambda_2^2\lambda_3^2 + \lambda_1^2\lambda_3^2 \quad \text{and} \quad I_3 = \lambda_1^2\lambda_2^2\lambda_3^2 \qquad (6) $$

123 The formulation in the nonlinear range is based on the decomposition of the gradient $F(X,t)$ into a volume-preserving and pure dilatational part as its originally proposed by Flory in [20] and used in several works such as [9] and [21] among others as follow:
124
125
126

$$ F = \bar{F} \det(F)^{1/3} I \quad \text{where} \det(\bar{F}) = 1 \qquad (7) $$
127 $\bar{F}$ is the volume-preserving gradient tensor, the advantage of such decomposition is its validity near and far away from the thermodynamic equilibrium [11]. The Cauchy-Green strain tensor associated and the Lagrangian strain tensor associated with the volume-preserving gradient are expressed as

$$ \bar{C} = \bar{F}^t \bar{F} = J^{-2/3} C, \quad \bar{E} = \frac{1}{2} (\bar{C} - I) \qquad (8) $$

131 $I$ is the metric tensor in the reference configuration. Furthermore, several applications of the chain rule lead to the following

$$ \frac{\partial \bar{E}}{\partial E} = \frac{\partial \bar{C}}{\partial C} = J^{-2/3} \left[ I - \frac{1}{3} C \otimes C^{-1} \right] \qquad (9) $$

133 $I$ is the fourth order unit tensor i.e. ($I_{ijkl} = 1$ if ($i = j = k = l$); else $I_{ijkl} = 0$), the sign $\otimes$ designates the tensorial product in the reference basis. Hence, we postulated an uncoupled free energy density as its expressed in [9] by a Taylor series in which terms higher than the second order are omitted.

$$ \Psi(\bar{C}, Q) = U^0(J) + \bar{\Psi}^0(\bar{C}) - \frac{1}{2}Q : \bar{C} + \Psi_I(Q) \qquad (10) $$

137 $Q$ is a second order overstress tensor internal variable akin to the second Piola-Kirchhoff stress tensor $S$. The first two terms of the free energy density of equation 10 are the dilatational and volume-preserving parts of the stored elastic energy density. The third and fourth terms are responsible for the time-dependent behavior of the material. Note that $\Psi_I(Q)$ is a convex function of the internal variable $Q$. This decomposition of the free energy density leads to a decomposition in the stress into a deviatoric (shear) and hydrostatic (bulk) parts.

## 3.2. Rate and constitutive equations

145 The rate equation of the internal variable $Q$ is motivated by the rate equation 4 of the rheological model in which the elastic stress is replaced by the deviatoric part of the hyperelastic Second Piola-Kirchhoff stress tensor as it's expressed in [9] for which the behavior in bulk is considered purely elastic.

$$ \frac{\partial Q}{\partial \xi} + \frac{1}{\tau} Q = \frac{\gamma}{\tau} DEV \left[ 2 \frac{\partial \Psi^0(\bar{C})}{\partial \bar{C}} \right] \quad \text{with} \quad \xi(t) = \int_0^t \frac{dt'}{a(\bar{C})} \qquad (11) $$
As in the previous section, $a(\bar{C})$ is a function of the invariants of the volume-preserving right Cauchy-Green strain tensor $\bar{C}$ and $\xi$ is referred to as the reduced time and its an increasing function of time. The second law of thermodynamic is expressed in terms of the Clausius-Duhem inequality in the reference configuration $C_0$.

$$ -\dot{\Psi} + \frac{1}{2}S : \dot{C} \geq 0 \qquad (12) $$

Standard arguments [22] and [23] using inequality 12 lead to the expression of intrinsic dissipation and the second Piola-Kirchhoff stress tensor.

$$ -\frac{\partial \Psi(C, Q)}{\partial Q} : \dot{Q} \geq 0 \quad \text{and} \quad S = \frac{1}{2} \frac{\partial \Psi(C, Q)}{\partial C} \qquad (13) $$

Let $\tau = J\sigma - FSF^t$ be the Kirchhoff stress tensor, where $\sigma$ is the Cauchy stress tensor and $S$ is the second Piola-Kirchhoff stress tensor. Considering relations 9, 11 and 13 one could simply lead to the convolution representation of the second Piola-Kirchhoff stress tensor.

$$ S = J^{-2/3} \int_{0}^{\xi} g(\xi - \xi') \frac{\partial}{\partial \xi'} DEV \left( 2 \frac{\partial \Psi^{0}(\bar{C})}{\partial \bar{C}} \right) d\xi' + JpC^{-1} \qquad (14) $$

in which $DEV(\bullet) = (\bullet) - \frac{1}{3}[C:(\bullet)]C^{-1}$ denotes the deviator operator in the reference configuration. $p = \partial U^0(J)/\partial J$ is the hydrostatic part of the stress, for an incompressible material, $p$ is an undetermined pressure to be obtained by the boundary conditions. $g$ is the normalized shear relaxation modulus and its a decaying function of time, its often expressed by a power law function or a decaying exponential function. For computational reasons its more efficient to consider the Cauchy stress tensor rather than the second Piola-Kirchhoff stress tensor. Application of an integration by parts to the expression of the second piola-Kirchhoff stress tensor of relation 14 and considering the relative distortional deformation gradient tensor $\bar{F}_t(t') = J^{-1/3}\partial\Phi(X, t')/\partial\Phi(X, t)$ the Cauchy stress tensor reads

$$ \sigma = \sigma_o^d + \frac{1}{J} dev \int_0^\xi \frac{\partial g(\xi')}{\partial \xi'} (\bar{F}_\xi^{-1}(\xi - \xi') \tau_o^d(\xi - \xi') \bar{F}_\xi^{-t}(\xi - \xi')) d\xi' + pI \qquad (15) $$

$\sigma_o$ designates the elastic instantaneous Cauchy stress tensor and may be written [23]

$$ \sigma_o = \beta_0 I + \beta_1 B + \beta_{-1} B^{-1} \qquad (16) $$
where $\beta_j = \beta_j(I_1, I_2, I_3)$ are the elastic response functions. In terms of the free energy density they are given by

$$ 
\begin{aligned}
\beta_0 (I_1, I_2, I_3) &= \frac{2}{J} [I_2 \Psi_1^0 + I_3 \Psi_3^0] \\
\beta_1 (I_1, I_2, I_3) &= \frac{2}{J} \Psi_1^0 \\
\beta_{-1} (I_1, I_2, I_3) &= -2J \Psi_2^0
\end{aligned} 
\quad (17) $$

where

$$ \Psi_k^0 = \frac{\partial \Psi^0}{\partial I_k}, \quad k = 1, 2, 3 \quad (18) $$

The free energy density has an alternative form in terms of the principle stretches given by

$$ \tilde{\Psi}^o (\lambda_1, \lambda_2, \lambda_3) = \Psi^o (I_1, I_2, I_3) \quad (19) $$

$\tau_o^d = dev(J\sigma_o)$ is the deviatoric part of the instantaneous elastic Kirchhoff stress tensor and may be expressed in lights of relations 17 and 18 by

$$ \tau_o^d = 2 \left[ \frac{1}{3} (I_2 \Psi_2^0 - I_1 \Psi_1^0) I + \Psi_1^0 B - I_3 \Psi_2^0 B^{-1} \right] \quad (20) $$

Or in lights of relation 19 in the base of principle stretches by

$$ \tau_{oi}^d = \lambda_i \tilde{\Psi}_i^o - \frac{1}{3} \sum_{j=1}^{3} \lambda_j \tilde{\Psi}_j^o \quad (21) $$

where $\tilde{\Psi}_i^o$ refers to the derivative $\tilde{\Psi}^o$ with respect to $\lambda_i$. The first term of the right hand side of 15 designates the instantaneous elastic response of the material, the second one denotes the time dependent part of the material whereas the third one is the hydrostatic pressure. For an incompressible material relation 15 holds with $J = 1$ and $\bar{F}_t(t') = F_t(t') = \partial\Phi(X, t')/\partial\Phi(X, t)$. Henceforth, the material is considered incompressible so that the Cauchy stress tensor $\sigma$ and the Kirchhoff stress tensor $\tau$ are the same quantity.

## 4. Model identification

In this section, a systematic identification of the material parameters for an incompressible nonfactorizable viscoelastic constitutive equation at finite strain is highlighted. This procedure relies on the separate identification of hyperelastic potential, viscoelastic kernel and the reduced time function.
Considering the form of the constitutive equation 14, each characteristic function identification reduces to the solution of a nonlinear optimization problem. The identification procedure is carried out considering homogenous (in space) uniaxial transformations of simple extension and pure shear. For such transformation, assuming incompressibility, in the material basis the deformation gradient tensor may be specified as

$$F(t) = diag(\lambda_1(t), \lambda_2(t), \lambda_3(t)) \qquad (22)$$

The condition of incompressibility $J = 1$ leads to the following expression of the deformation gradient tensor

$$F(t) = diag(\lambda(t), \lambda^{-1/2}(t), \lambda^{1/2}(t)) \qquad (23)$$

For simple extension and

$$F(t) = diag(\lambda(t), \lambda^{-1}(t), 1) \qquad (24)$$

For pure shear. In both cases only one component of the stress remains and the indeterminate hydrostatic pressure is eliminated i.e.

$$\sigma_2 = \sigma_3 = 0 \quad \text{and} \quad \sigma_1 = \sigma_1 - \sigma_2 \qquad (25)$$

From 15 and 25 the expression of the stress follows for these histories of deformations considered.

$$\sigma = \sigma_{01}^{d} - \sigma_{02}^{d} + dev \int_{0}^{\xi} \frac{\partial g(\xi')}{\partial \xi'} \left( \left[ \frac{\lambda^{2}(\xi)}{\lambda^{2}(\xi - \xi')} \sigma_{01}^{d} \right] - \left[ \frac{\lambda(\xi - \xi')}{\lambda(\xi)} \sigma_{02}^{d} \right] \right) d\xi' \qquad (26)$$

For the simple extension and

$$\sigma = \sigma_{01}^{d} - \sigma_{02}^{d} + dev \int_{0}^{\xi} \frac{\partial g(\xi')}{\partial \xi'} \left( \left[ \frac{\lambda^{2}(\xi)}{\lambda^{2}(\xi - \xi')} \sigma_{01}^{d} \right] - \left[ \frac{\lambda^{2}(\xi - \xi')}{\lambda^{2}(\xi)} \sigma_{02}^{d} \right] \right) d\xi' \qquad (27)$$

For pure shear. The first two terms of relations 26 and 27 refer to the principle components of the deviatoric instantaneous elastic part of the stress which can be obtained from the equilibrium deviatoric elastic stress via:

$$\sigma_{0i}^{d} = \frac{G_{0}}{G_{\infty}} \sigma_{\infty i}^{d} \qquad (28)$$
in which $G_0$ and $G_\infty$ refer to the instantaneous and equilibrium shear relaxation modulus see [24] and [25], whereas the integral term depicts the dissipative or the time dependent part of the stress. A general identification procedure could be applied separately to each component of the stress. Hence, let $\Lambda = [\Lambda_1, \Lambda_2, ..., \Lambda_m]^t$ be the vector of experimental input data and $\Theta = [\Theta_1, \Theta_2, ... \Theta_m]^t$ be the vector of corresponding experimental response. For each component of the stress the response function is written $F(\Lambda, p) : \mathbb{R} \times \mathbb{R}^n \rightarrow \mathbb{R}$ in which $p = [p_1, p_2, ..., p_n]^t$ is a vector of material parameters. The objective function is defined by the squared 2-norm

$$S_F(p) : \|F(\Lambda, p) - \Theta\|_2^2 = \sum_{i=1}^{m} [F(\Lambda_i, p) - \Theta_i]^2 \quad (29)$$

The identification procedure turns out into a minimization problem which reads as follow

$$\min_{p} S_F(p) \qquad (30)$$

### 4.1. Identification of the hyperelastic potential

The free energy density $\Psi^o$ is a function of, either, the invariants of the right Cauchy-Green strain tensor or the principles stretches. The condition of incompressibility reads

$$\lambda_1\lambda_2\lambda_3 = 1 \quad \text{or} \quad I_3 = J^2 = 1 \qquad (31)$$

The general form of Mooney and Rivlin [26] free energy density is considered which reads for an incompressible hyperelastic material as follow

$$\Psi^o(I_1, I_2) = \sum_{i,j} c_{ij}(I_1 - 3)^i (I_2 - 3)^j \qquad (32)$$

$c_{ij}$ are the material parameters of the free energy density which usually should satisfy the stability conditions which ensure an admissible response of the model for any process see [27] these conditions are expressed as follow:

$$\frac{\partial \Psi^o}{\partial I_1} > 0 \quad \text{and} \quad \frac{\partial \Psi^o}{\partial I_2} \geq 0 \qquad (33)$$

In the case of uniaxial experiment, the nominal stress which is the measured quantity, actual force over reference area, and the principle stretch are related
234 through the free energy $\Psi^o$ by the equation

$$ \Pi^o = \frac{\partial \Psi^o}{\partial \lambda} = \sum_{i,j} c_{ij} \phi(i, j, \lambda) \qquad (34) $$

235 where $\phi(i, j, \lambda)$ is a nonlinear function of $i, j$ and $\lambda$

$$ \phi(i, j, \lambda) = 2i \left(\lambda - \frac{1}{\lambda^2}\right) \left(\lambda^2 + \frac{2}{\lambda} - 3\right)^{i-1} \left(2\lambda + \frac{1}{\lambda^2} - 3\right)^j + 2j \left(1 - \frac{1}{\lambda^3}\right) \left(\lambda^2 + \frac{2}{\lambda} - 3\right)^i \left(2\lambda + \frac{1}{\lambda^2} - 3\right)^{j-1} \qquad (35) $$

236 For the simple extension and

$$ \phi(i, j, \lambda) = 2(i + j) \left( \lambda - \frac{1}{\lambda^3} \right) \left( \lambda^2 + \frac{1}{\lambda^2} - 2 \right)^{i+j-1} \qquad (36) $$

237 For pure shear. An alternative useful representation of equation 34 with respect to the identification procedure is used.

238 Let $c^t = \{c_{01}, \dots, c_{0j}, c_{10}, \dots, c_{1j}, \dots, c_{i0}, \dots, c_{ij}\}$ be the vector of material parameters, $\Phi$ be a matrix representation of the function $\phi(i, j, \lambda)$ and $\Pi^e$ the discrete vector of nominal stress. Equation 34 became

$$ \Pi^e = \Phi c \qquad (37) $$

242 The identification of the material parameters $c_{ij}$ is performed using data for simple extension and pure shear simultaneously. Therefore, a modification of the objective function 29 is adopted [28]. The new objective function reads as follows

$$ \min_{c \in \mathbb{R}^{i,j}} \left( \left\| \Phi^{se} c - \tilde{\Pi}^{se} \right\|_2^2 + \left\| \Phi^{ps} c - \tilde{\Pi}^{ps} \right\|_2^2 \right) \qquad (38) $$

246 The superscript $se$ and $ps$ refers to the simple extension and pure shear respectively. $\tilde{\Pi}$ denotes the recorded experimental nominal stress vector. A least square minimization procedure is then employed under conditions 33 using Matlab software to reach the numerical values of $c_{ij}$. The results of this identification and its efficiency are discussed in the latter section of this work.

## 4.2. Identification of the viscoelastic kernel

252 The time dependent part of the stress is characterized by the shear relaxation function $G(\xi)$ which is a decaying positive function of the reduced time
ξ. Its often expressed by, either, a sum of decaying exponential functions called Prony series function or a power law functions. This identification is performed using experimental results from relaxation tests and dynamic tests in the linear range of behavior so that the reduced time is equal to the real time ξ = t and the behavior of the material is described by the single convolution integral of Boltzmann:

$$
\sigma(t) = \int_{0}^{t} G(t-t') \dot{\varepsilon}(t') dt' \qquad (39)
$$

ϵ is the linearized strain tensor.

4.2.1. Identification from relaxation test

The relaxation test is performed in shear deformation. The strain is
suddenly increased to a value ε₀ and kept constant

$$
\epsilon(t) = H(t) \epsilon_o \quad \text{with} \quad H(t) = \begin{cases} 0, & t < 0 \\ 1, & t > 0 \end{cases} \tag{40}
$$

From equations 39 and 40 the shear relaxation modulus follows

$$
G(t) = \frac{\sigma(t)}{\epsilon_o} \tag{41}
$$

In this work we adopted the Prony series form of the shear relaxation modulus

$$
G(t) = G_{\infty} + \sum_{i=1}^{N} G_i \exp \left( -\frac{t}{\tau_i} \right) \quad (42)
$$

G∞ denotes the long term shear relaxation modulus, Gᵢ i = 1, .., N are the coefficients of the Prony series and τᵢ i = 1, .., N are the relaxation time constants. Furthermore, in order to avoid the ill-conditioning of the optimization problem the set of the relaxation times τᵢ are a-priori fixed as one time constant per decade in the logarithmic time scale for the experimental time window [29] and [30]. The optimization problem arising from the identification of the N-terms Prony coefficients is

$$
\min_{\{G\} \in \mathbb{R}^N} \left\| \Gamma\{G\} - \hat{\Gamma} \right\|_2^2 \quad (43)
$$
where $\Gamma \in \mathbb{R}^{MN}$ is the matrix representation of relation 42 taking the following form

$$ \Gamma = \begin{bmatrix} 1 & \exp(-t_1/\tau_1) & \dots & \exp(-t_1/\tau_N) \\ 1 & \exp(-t_2/\tau_1) & \dots & \exp(-t_2/\tau_N) \\ \vdots & \vdots & \ddots & \vdots \\ 1 & \exp(-t_M/\tau_1) & \dots & \exp(-t_M/\tau_N) \end{bmatrix} \quad (44) $$

$t = \{t_1, \dots, t_M\}$ are the discrete time instants and $\hat{G} = \{\hat{G}_1, \dots, \hat{G}_M\}$ are the corresponding experimental values of the shear relaxation modulus using relation 41. A linear least square algorithm is used to solve the optimization problem 43 using Matlab software.

### 4.2.2. Identification from dynamic tests

The dynamic tests are performed using a cylindrical shear sheet loaded by a sinusoidal deformation without a predeformation and with small amplitude

$$ \varepsilon(t) = \varepsilon_o \exp(j\omega t) \quad \text{with } \varepsilon_o \ll 1 \quad (45) $$

$\omega$ is the circular frequency and $j$ is the unit imaginary number. Hence, from equations 39 and 45 the stress-strain relation follows

$$ \sigma = G^* \varepsilon_o \quad (46) $$

$G^*$ is the complex dynamic shear modulus, its real and imaginary parts are denoted $G'$ and $G''$ are called storage and loss modulus respectively and may be obtained by a Fourier transform of equation 42 and given by :

$$ \begin{aligned} G' &= G_{\infty} + \sum_{i=1}^{N} G_{i} \frac{(\tau_{i}\omega)^{2}}{1+(\tau_{i}\omega)^{2}} \\ G'' &= \sum_{i=1}^{N} G_{i} \frac{\tau_{i}\omega}{1+(\tau_{i}\omega)^{2}} \end{aligned} \quad (47) $$

As mentioned in the previous section, the relaxation times $\tau_i$ are a-priori fixed as one time constant per decade in the logarithmic scale of time. Thereby, both storage and loss modulus are linear with respect to the N-terms Prony coefficients. The arising optimization problem from this identification procedure reads

$$ \min_{\{G\} \in \mathbb{R}^N} \left( \| \Gamma' \{G\} - \hat{G}' - G_{\infty} \|_2^2 + \| \Gamma'' \{G\} - \hat{G}'' \|_2^2 \right) \quad (48) $$
G∞ is directly identified from the storage modulus curve as $\omega \to 0$. $\hat{G}''$ and $\hat{G}'''$ are the experimental vectors of storage and loss modulus, as recorded by the DMA machine, respectively. $\Gamma'$ and $\Gamma''$ are two M by N matrices representing equation 47 and can be expressed through the relaxation time constants and the discrete frequency vector by

$$ 
\begin{aligned}
\Gamma' &= \begin{bmatrix} \frac{(\tau_1 \omega_1)^2}{1+(\tau_1 \omega_1)^2} & \cdots & \frac{(\tau_N \omega_1)^2}{1+(\tau_N \omega_1)^2} \\ \vdots & \ddots & \vdots \\ \frac{(\tau_1 \omega_M)^2}{1+(\tau_1 \omega_M)^2} & \cdots & \frac{(\tau_N \omega_M)^2}{1+(\tau_N \omega_M)^2} \\ \frac{\tau_1 \omega_1}{1+(\tau_1 \omega_1)^2} & \cdots & \frac{\tau_N \omega_1}{1+(\tau_N \omega_1)^2} \\ \vdots & \ddots & \vdots \\ \frac{\tau_1 \omega_M}{1+(\tau_1 \omega_M)^2} & \cdots & \frac{\tau_N \omega_M}{1+(\tau_N \omega_M)^2} \end{bmatrix} \\
\Gamma'' &= \begin{bmatrix} \frac{(\tau_1 \omega_1)^2}{1+(\tau_1 \omega_1)^2} & \cdots & \frac{(\tau_N \omega_1)^2}{1+(\tau_N \omega_1)^2} \\ \vdots & \ddots & \vdots \\ \frac{(\tau_1 \omega_M)^2}{1+(\tau_1 \omega_M)^2} & \cdots & \frac{(\tau_N \omega_M)^2}{1+(\tau_N \omega_M)^2} \\ \frac{\tau_1 \omega_1}{1+(\tau_1 \omega_1)^2} & \cdots & \frac{\tau_N \omega_1}{1+(\tau_N \omega_1)^2} \\ \vdots & \ddots & \vdots \\ \frac{\tau_1 \omega_M}{1+(\tau_1 \omega_M)^2} & \cdots & \frac{\tau_N \omega_M}{1+(\tau_N \omega_M)^2} \end{bmatrix}
\end{aligned} 
\quad (49) $$

The optimization problem of equation 48 is an ill-posed problem [31]. Therefore, a Tikhonov [32] regularization method was employed to solve this system. The results of this identification using randomly perturbed simulated and real experimental data are shown in the latter section of this paper.

### 4.3. Identification of the reduced time function

Once the hyperelastic potential and the viscoelastic kernel are identified, the problem of determining the reduced time function can be addressed. This identification relies on the discretization of the stress-strain relation (equation 15) with respect to the time. Let $t = \{t_1, ..., t_M\}$ be the discrete experimental time vector and $\xi = \{\xi_1, ..., \xi_M\}$ be the corresponding reduced time vector, $\Delta t$ is the experimental time increment and $\Delta \xi$ is the reduced time increment. The general form of this discretization formula [33] for a nonlinear viscoelastic behavior as follows is described in equation 50. The identification of the reduced time vector $\xi$ is performed thanks to a recursive dichotomy algorithm applied to the error between the discretized stress 50 and the experimental stress $\tilde{\sigma} = \{\tilde{\sigma}_1, ..., \tilde{\sigma}_M\}$.

$$ 
\begin{aligned}
\sigma(t_{n+1}) &= \sigma_o^d(t_{n+1}) - \sum_{i=1}^{N} \sigma_i^d(t_{n+1}) + pI \\
\sigma_i^d(t_n) &= \frac{g_i}{\tau_i} \int_0^\xi \bar{F}_\xi^{-1}(\xi - \xi') \sigma_o^d(\xi - \xi') \bar{F}_\xi^{-t}(\xi - \xi') \exp\left(-\frac{\xi'}{\tau_i}\right) d\xi' \\
\sigma_i^d(t_{n+1}) &= \alpha_i g_i \sigma_o^d(t_{n+1}) + \beta_i g_i \hat{\sigma}_o^d(t_n) + \gamma_i \hat{\sigma}_i^d(t_n)
\end{aligned} 
\quad (50) $$

with

$$ 
\begin{aligned}
\gamma_i &= \exp\left(-\frac{\Delta\xi}{\tau_i}\right); \quad \alpha_i = 1 - \frac{\tau_i}{\Delta\xi}(1 - \gamma_i); \quad \beta_i = \frac{\tau_i}{\Delta\xi}(1 - \gamma_i) - \gamma_i \\
\hat{\sigma}_j^d(t) &= \bar{F}_t(t + \Delta t) \sigma_j^d(t) \bar{F}_t^t(t + \Delta t); \quad j = 0, 1, \dots
\end{aligned} 
$$
Once the reduced time vector $\xi$ is obtained the identification of the reduced time function $a(C)$ can be addressed since it is the inverse of the derivative of the reduced time with respect to real time:

$$ \frac{1}{a} = \frac{d\xi}{dt} \qquad (51) $$

The derivative in equation 51 is obtained numerically since the reduced time and the real time are two discrete vectors. Hence, one leads to the numerical vector of function $a(C) : a = \{a_1, ..., a_M\}$. Furthermore, a sufficient condition on this function with respect to the second principle of thermodynamics in terms of Clausius-Duhem inequality is to adopt a positive function of the invariants of the right Cauchy-Green strain tensor [14]. A suitable form of this function is denoted by:

$$ a(C) = \exp[-c_1(I_1 - 3) - c_2(I_2 - 3)] \qquad (52) $$

where $c_1$ and $c_2$ are two positive material parameters to be fitted using a nonlinear curve fitting algorithm with Matlab software.

## 5. Identification of the model using data from the Pipkin isotropic model

In this section, the capacity of the proposed model to depict the response of other complicated viscoelastic models is presented. The main concern is to reformulate a complicated model namely the isotropic viscoelastic model by Pipkin [12] in the form of our simple model presented herein. To this end the identification procedure outlined above is applied using data generated from the isotropic viscoelastic model proposed by Pipkin [12] see equations 53 and 54. Data were generated from the stress-strain relation in the case of simple extension and pure shear experiments. Several strain histories were considered to provide a complete description of the behavior. The hyperelastic potential was identified using data of simple extension and pure shear at equilibrium, the relaxation function was obtained using a relaxation test performed in simple extension and the reduced time was calculated using monotonic test in simple extension for different strain rates. The identification procedure is validated by predicting the behavior in pure shear monotonic tests for different strain rates.
### 5.1. Pipkin isotropic model

Pipkin [12] proposed a third order development of the tensorial response function $Q$ for an isotropic incompressible material. The principle of material indifference requires that the Cauchy stress tensor takes the following form:

$$ \sigma = RQR^t + pI \quad (53) $$

$R$ is the rotation tensor obtained from the polar decomposition of the transformation gradient tensor $F$ and $p$ is the indeterminate parameter due to incompressibility. The third functional development of $Q$ reads

$$ \begin{aligned} Q(t) = & \int_0^t r_1(t-t') \dot{E}(t') dt' + \int_0^t \int_0^t r_2(t-t', t-t'') \dot{E}(t') \dot{E}(t'') dt' dt'' + \\ & \int_0^t \int_0^t \int_0^t r_3(t-t', t-t'', t-t''') \operatorname{tr}[\dot{E}(t') \dot{E}(t'')] \dot{E}(t'') dt' dt'' dt''' + \\ & \int_0^t \int_0^t \int_0^t r_4(t-t', t-t'', t-t''') \dot{E}(t') \dot{E}(t'') \dot{E}(t''') dt' dt'' dt''' \end{aligned} \quad (54) $$

$r_i$ are the relaxation kernels expressed by a decaying exponential functions and $\dot{E}(t)$ is the time derivative of the Green-Lagrange deformation tensor $E = 1/2(C-I)$. Expression of $r_i$ according to [34] is reported in equation 55, the choice of $r_2(t_1, t_2) = 0$ is motivated by thermodynamic arguments to ensure the positivity of the free energy density. Further arguments could be found in [34] and references therein.

$$ \left\{ \begin{aligned} r_1(t) &= a_1 + b_1 \exp(c_1 t) \\ r_2(t_1, t_2) &= 0 \\ r_3(t_1, t_2, t_3) &= a_3 + b_3 \exp(c_3 (t_1 + t_2 + t_3)) \\ r_4(t_1, t_2, t_3) &= b_4 \exp(c_4 (t_1 + t_2 + t_3)) \end{aligned} \right. \quad (55) $$

A crucial choice of the parameters $a_i$, $b_i$ and $c_i$ enables us to describe the behaviour of the material for any given strain history.

### 5.2. Identification results

#### 5.2.1. Hyperelastic potential

The identification of the free energy density requires data at equilibrium in the case of simple extension and pure shear experiments. Hence, data were generated by omitting the time-dependent part of the stress. Considering the incompressibility of the behaviour of equations 23 and 24 its straightforward to obtain from 54 the relations for the equilibrium stress

$$ \sigma = \left( \lambda^2 - \frac{1}{\lambda} \right) \left[ \frac{a_1}{2} + \frac{a_3}{8} \left( \lambda^4 - 2\lambda^2 - \frac{4}{\lambda} + \frac{2}{\lambda^2} + 3 \right) \right] \quad (56) $$
³⁶⁵ in the case of simple extension and

$$ \sigma = \left( \lambda^2 - \frac{1}{\lambda^2} \right) \left[ \frac{a_1}{2} + \frac{a_3}{8} \left( \lambda^4 - 2\lambda^2 - \frac{2}{\lambda^2} + \frac{1}{\lambda^4} + 2 \right) \right] \quad (57) $$
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Equilibrium stresses versus principle stretch for the Pipkin model (diamond) and the Mooney-Rivlin model (solid curve)</figcaption></figure>
<table><caption>Table 1: Prony series parameters</caption><thead><tr><th>g<sub>i</sub></th><th>τ<sub>i</sub>(s)</th></tr></thead><tbody><tr><td>6.25 10<sup>-2</sup></td><td>2.003</td></tr><tr><td>2.84 10<sup>-5</sup></td><td>14.06</td></tr><tr><td>1.12 10<sup>-4</sup></td><td>82.76</td></tr></tbody></table>

For the simple shear. Results of the identification using the generalized Mooney-Rivlin free energy density in terms of the first Piola-Kirchhoff stress are reported in figure 3 for simple extension and pure shear experiments. A second order generalized Mooney-Rivlin potential, in relation 32, was satisfactory to describe the hyperelastic part of the Pipkin model.

### 5.2.2. Viscoelastic kernel

The identification of the Prony series requires shear relaxation data at low level of strain. To this end, a Heaviside strain history of relation 40 is considered. Introduction of this strain history into 53 and 54 yields the relaxation stress-strain relationship.

$$ \sigma(t) = \frac{r_1(t)}{2} \left(\lambda^2 - \frac{1}{\lambda^2}\right) + \frac{r_3(3t)}{8} \left(\lambda^2 - \frac{1}{\lambda^2}\right) \left(\lambda^4 - 2\lambda^2 - \frac{2}{\lambda^2} + \frac{1}{\lambda^4} + 2\right) + \frac{r_4(3t)}{8} \left(\lambda^2 - \frac{1}{\lambda^2}\right) \left(\lambda^4 - 3\lambda^2 - \frac{3}{\lambda^2} + \frac{1}{\lambda^4} + 4\right) \quad (58) $$

In figure 4 are reported curves of the normalized shear relaxation modulus versus time for four different levels of strain. Its well shown that the hypothesis of separability doesn't hold for the Pipkin model since the normalized shear relaxation modulus depends upon strain level. It should be noted, however, that the identification procedure was performed using results of the 5% level of strain. Prony series parameters are reported in Table 1.
<figure><img src="image_6.png" /><figcaption>Figure 4: Normalized shear relaxation modulus of Pipkin model versus time for four different strain levels.</figcaption></figure>

### 5.2.3. Reduced time function

In this part monotonic tests of simple extension and pure shear were generated from the Pipkin model. Simple extension test was used in the identification of the reduced time function whereas pure shear test was used in the validation of the results. For computational convenience with respect to the multi-integral form involved in 54, the principle stretch corresponding to a monotonic test was set to be an increasing exponential function of time of the form:

$$ \lambda(t) = \exp(\alpha t) \qquad (59) $$
<figure><img src="image_7.png" /><figcaption>Figure 5: Simple extension Cauchy stress versus principle stretch for two different strain rates $\alpha_1 = 1.19 \ 10^{-2}$ and $\alpha_2 = 6 \ 10^{-3}$</figcaption></figure>
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: reduced time function and reduced time ratio versus principle stretch for two strain rates $\alpha_1 = 1.19 \ 10^{-2}$ and $\alpha_2 = 6 \ 10^{-3}$</figcaption></figure>
where $\alpha$ is a positive constant could be interpreted as a strain rate. Replacing the principle stretch in 54 by its expression yields the expression of the Cauchy stress. For each test, two different strain rates were considered $\alpha_1 = 1.19 \ 10^{-2}$ and $\alpha_2 = 6 \ 10^{-3}$. Data for simple extension Cauchy stress are plotted versus principle stretch in figure 5. Hence, the identification procedure highlighted above was applied to identify the reduced time function. Results are reported in figure 6 by means of the reduced time ratio $\xi(t)/t$ and the reduced time function $a(C)$ for the two strain rates considered. The reduced time function is obtained numerically via a numerical derivation of the reduced time with respect to time. It shows a significant dependence upon strain level which is consistent with the results shown in Figure 4. Furthermore, this function is independent of the strain rate which motivate the choice of the form of the reduced time function of equation 52. The capacity of the nonlinear viscoelastic model developed herein is evaluated by predicting the behaviour of the Pipkin model using the parameters identified in this section. In order to avoid a division by small value of the force when the principle stretch is near to one, a modified relative error formula was used as proposed in [35]

$$err_i = \frac{|\sigma_i - \sigma_i^P|}{\max\{0.5, \sigma_i^P\}} \qquad (60)$$
<figure><img src="image_10.png" /><figcaption>Figure 7: Relative error of the predicted Cauchy stress of the Pipkin model for pure shear experiment</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 8: Pure shear Cauchy stress for the model (solid curve) and the Pipkin model (diamond and square)</figcaption></figure>

in which $\sigma_i$ is the Cauchy stress computed using 15 and $\sigma_i^P$ is the Pipkin
Cauchy stress computed using 53 and 54. This function is plotted versus the
principle stretch in figure 7 in the case of pure shear experiment. For the
two strain rates considered, the relative error remains under . In figure 8 is
plotted the Cauchy stress versus principle stretch for the Pipkin model and
the proposed model.

# 6. Application of the identification procedure to experimental data

In this section, the identification procedure outlined in section 4 is used to
identify the parameters of the proposed model using experimental data for a
bromobutyl (BIIR) rubber material. Experimental data used here are those
from [13], in which a complete experimental characterization was performed
to obtain the response of the material for several strain history configuration
and several temperatures. In what follows, results of the identification of the
model' parameters are highlighted and discussed.
6.1. *Hyperelastic potential*

<figure><img src="image_12.png" /><figcaption>Figure 9: Equilibrium stresses versus principle stretch: Experimental (diamond) and the identified Mooney-Rivlin model (solid curve)</figcaption></figure>
<table><caption>Table 2: Prony series parameters for BIIR rubber</caption><thead><tr><td>g<sub>i</sub></td><td>τ<sub>i</sub>(s)</td></tr></thead><tbody><tr><td>4.46 10<sup>-3</sup></td><td>14.79</td></tr><tr><td>3.77 10<sup>-2</sup></td><td>125.71</td></tr><tr><td>5.69 10<sup>-2</sup></td><td>460.7</td></tr><tr><td>5.84 10<sup>-2</sup></td><td>1761.6</td></tr><tr><td>8.76 10<sup>-2</sup></td><td>9598.5</td></tr></tbody></table>

The identification of the free energy density coefficients of relation 38 is performed under stability conditions of the relation 33 using Matlab software. A second order Mooney-Rivlin potential was found to be satisfactory with respect to the value of the total relative error of the stress for simple extension and pure shear experiments. In figure 9 are plotted experimental and identified Piola-Kirchhoff stresses versus principle stretch at equilibrium for simple extension and pure shear. The relative error of the relation 60 was calculated for both experiments, its average value is 0,5 % for simple extension and 2,3 % for pure shear which are very satisfactory considering the non-linearity of the material.

## 6.2. Viscoelastic kernel

The identification of the viscoelastic kernel, as it's described in section 4, is performed using two different experimental data: shear relaxation experiment in the linear range of the behavior and dynamic tests for low level of dynamic amplitude and without pre-strain. In what follows, results of this identification procedure are discussed.

### 6.2.1. From shear relaxation experiment

The shear relaxation experiment is performed in simple shear deformation at a strain level of 10 %. Despite that the value of the strain is quite big, it's considered, however, in the linear range of the behavior since the material is highly deformable. From figure 1 data for shear relaxation experiment at 10 % are extracted and used in the identification of the Prony series parameters. These parameters are reported in table 2. The average relative error between the experimental and the identified viscoelastic kernel is in the order of 0.1 %. The identified and the experimental normalized shear relaxation modulus are plotted versus time in figure 10.
<figure><img src="image_13.png" /><figcaption>Figure 10: Normalized shear relaxation modulus of BIIR rubber versus time.</figcaption></figure>

### 6.2.2. From dynamic experiments

The dynamic experiments are performed in simple shear deformation with a small dynamic amplitude and without a pre-strain of the form of equation 45. It's recalled that the problem of the identification of the viscoelastic parameters from dynamic data of relation 48 is an ill-posed problem. Hence, a regularization procedure of Tikhonov is used. In what follows, this method is recalled and applied to theoretical dynamic data using parameters from [36] and then applied to dynamic data from [13].

* Tikhonov regularization method:

The linear system arising from the identification of the Prony series parameters from dynamic data is an ill-posed problem [31]. From the original system of equation 48 the following system arise:

$$Ax = b \qquad (61)$$

in which A is the global matrix of the system to be calculated from 49, b is the vector of experimental storage modulus and loss modulus vectors $\hat{G}'$ and $\hat{G}''$ and x is the vector of the Prony series parameters $G_i, i = 1..N$. Tikhonov regularization method replaces system 61 by:

$$(A^t A + \mu I) x = A^t b \qquad (62)$$
in which $\mu > 0$ is the regularization parameter and $I$ is the identity matrix. The regularization parameter is determined via a L-curve technique using Matlab software. It's well established that the solution of system 62 noted $x_\mu$ gives the minimum residual for the minimization problem arising from system 61 which means:

$$ \forall x \in \mathbb{R}^N \qquad (63) \\ \|Ax_\mu - b\|^2 \le \|Ax - b\|^2, \text{ such that } \|x\|^2 \le \|x_\mu\|^2 $$

The proof of 63 and further development of the convergence of the regularized Tikhonov problem are well studied in [37].

• Application of the Tikhonov method to simulated dynamic data:
<figure><img src="image_14.png" /><figcaption>(a) Storage modulus versus frequency for simulated and perturbed data</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Loss modulus versus frequency for simulated and perturbed data</figcaption></figure>

Figure 11: Storage and loss moduli versus frequency for simulated and perturbed data
<table><caption>Table 3: Prony series parameters from [36]</caption><thead><tr><th>G<sub>i</sub>(Pa)</th><th>τ<sub>i</sub>(s)</th></tr></thead><tbody><tr><td>1.94 10<sup>8</sup></td><td>2 10<sup>-2</sup></td></tr><tr><td>2.83 10<sup>8</sup></td><td>2 10<sup>-1</sup></td></tr><tr><td>5.54 10<sup>8</sup></td><td>2 10<sup>0</sup></td></tr><tr><td>6.02 10<sup>8</sup></td><td>2 10<sup>1</sup></td></tr><tr><td>3.88 10<sup>8</sup></td><td>2 10<sup>2</sup></td></tr><tr><td>1.56 10<sup>8</sup></td><td>2 10<sup>3</sup></td></tr><tr><td>4.1 10<sup>7</sup></td><td>2 10<sup>4</sup></td></tr><tr><td>1.38 10<sup>7</sup></td><td>2 10<sup>5</sup></td></tr><tr><td>3.68 10<sup>6</sup></td><td>2 10<sup>6</sup></td></tr><tr><td>7.9 10<sup>5</sup></td><td>2 10<sup>7</sup></td></tr><tr><td>9.6 10<sup>5</sup></td><td>2 10<sup>8</sup></td></tr></tbody></table>

The Tikhonov regularization procedure described above was applied to a simulated dynamic data generated using Prony series parameters from [36] and relations in (47). Further, in order to test the capacity of the method to deal with noisy experimental data, the second member of system 62 was perturbed randomly as follow:

$$ \tilde{b} = (1 \pm \varepsilon) b \qquad (64) $$

in which $\varepsilon$ takes three different values: 10%, 15% and 20%. The Prony series parameters from [36] are reported in table 3 with an equilibrium modulus $G_{\infty} = 2.24 10^6 Pa$. The results of this identification are reported in figure 11 in terms of the dynamic moduli versus frequency for the perturbed and original simulated data. The mean relative error for the three perturbed data remains under 10% and hence this procedure shows a huge capacity to predict the dynamic response functions despite the perturbation of the second member of the system 63.

* Application of the Tikhonov method to dynamic data from [13]:

The dynamic experiment are represented by the storage and loss moduli $G''$ and $G''$ as functions of the frequency for the experimental frequency window of $[10^{-5}, 10^7 Hz]$. The tikhonov regularization is applied and a regularization parameter $\mu = 2.229$ is obtained for a 12 parameters Prony series reported in table 4. The mean relative error is about
20%, this is due to the fact that the storage and loss moduli were approximated via the WLF law [38] following the time-temperature superposition principle.
<figure><img src="image_16.png" /><figcaption>(a) Storage modulus versus frequency for dynamic experimental data</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Loss modulus versus frequency for dynamic experimental data</figcaption></figure>

Figure 12: Storage and loss moduli versus frequency for dynamic experimental data for different values of the regularization parameter $\mu$
<table><caption>Table 4: Prony series parameters from experimental dynamic data</caption><thead><tr><th>G<sub>i</sub>(Pa)</th><th>τ<sub>i</sub>(s)</th></tr></thead><tbody><tr><td>3.13 10<sup>8</sup></td><td>7.1 10<sup>-8</sup></td></tr><tr><td>2.61 10<sup>7</sup></td><td>1 10<sup>-6</sup></td></tr><tr><td>9.86 10<sup>6</sup></td><td>1.4 10<sup>-5</sup></td></tr><tr><td>1.94 10<sup>6</sup></td><td>1.98 10<sup>-4</sup></td></tr><tr><td>9.42 10<sup>5</sup></td><td>2.8 10<sup>-3</sup></td></tr><tr><td>3.03 10<sup>5</sup></td><td>3.92 10<sup>-2</sup></td></tr><tr><td>1.4 10<sup>5</sup></td><td>5.52 10<sup>-1</sup></td></tr><tr><td>8.53 10<sup>4</sup></td><td>7.77 10<sup>0</sup></td></tr><tr><td>8.45 10<sup>4</sup></td><td>1.09 10<sup>2</sup></td></tr><tr><td>7.72 10<sup>4</sup></td><td>1.53 10<sup>3</sup></td></tr><tr><td>7.71 10<sup>4</sup></td><td>2.16 10<sup>4</sup></td></tr><tr><td>2.06 10<sup>4</sup></td><td>3.05 10<sup>5</sup></td></tr></tbody></table>

### 6.3. Reduced time function

The reduced time is identified using the discretization formula 50 and monotonic experiments of simple extension for two strain rates: 100% s⁻¹ and 200% s⁻¹. Cauchy stress versus principle stretch are plotted in figure 13. Results of the identification are reported in figure 14 in terms of the reduced time coefficient which is a nonlinear function of time for the strain rates considered. The pure shear experiment was predicted using the reduced time, the predicted and experimental data for this experiment are plotted in figure 15 against the principle stretch. From this result the relative error of the Cauchy stress is calculated using relation 60, its mean value remains under 2.5%. Hence, the proposed model is suitable to describe the material’s behavior at low and moderate strains.
<figure><img src="image_18.png" /><figcaption>Figure 13: Cauchy stress versus principle stretch for simple extension experiment</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 14: Reduced time coefficient for two strain rates 100% s⁻¹ and 200% s⁻¹</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 15: Cauchy stress versus principle stretch for pure shear experiment: experimental (diamond) model (solid curve)</figcaption></figure>

## 7. Conclusion

A three-dimensional viscoelastic model at finite strain that incorporate a strain dependent relaxation times has been proposed to describe nonfactorizable behavior of rubber-like materials. The model is based upon the internal state variables approach and the framework of rational thermodynamics and experimental arguments. The free energy density is decomposed into a volumetric and deviatoric parts. Furthermore, thermodynamic restrictions are fulfilled via a sufficient condition on the model's parameters resulting from the application of the Clausius-Duhem inequality for an arbitrary process. The identification of several functions involved in the model was, then, addressed. Each function's identification procedure turns out to the resolution of a linear or nonlinear system. Moreover, a regularization procedure of Tikhonov was applied in the resolution of the ill-posed problem arising from the identification of the viscoelastic kernel from dynamic data. This identification procedure was applied to a generated data from a multi-integral viscoelastic model and a static and dynamic experimental characterization of a BIIR rubber. The results of the identification has shown a huge capacity of the model to describe both experimental and multi-integral
526

527

viscoelastic model in the time domain of the behavior within and outside the
range of strain rates considered in deriving the model.

528

529

Our main object was to develop a constitutive equation compatible with the
second law of thermodynamics and suitable with the Finite-elements theory.

530

531

The implementation of the model within Abaqus software and the simulation
of real component with complex industrial application is the goal of future
work.

532

533

References

534

535

536

[1] R. A. Schapery, An engineering theory of nonlinear viscoelasticity with applications, International Journal of Solids and Structures 2 (1966) 407–425.

537

538

539

[2] W. Chang, R. Bloch, N. Tschoegl, Time-dependent response of soft polymers in moderately large deformations, Proceedings of the National Academy of Sciences 73 (1976) 981–983.

540

541

542

[3] W. Chang, R. Bloch, N. Tschoegl, Study of the viscoelastic behavior of uncrosslinked (gum) rubbers in moderately large deformations, Journal of Polymer Science: Polymer Physics Edition 15 (1977) 923–944.

543

544

545

[4] J. Sullivan, A nonlinear viscoelastic model for representing nonfactorizable time-dependent behavior in cured rubber, Journal of Rheology (1978-present) 31 (1987) 271–295.

546

547

[5] R. A. Schapery, On the characterization of nonlinear viscoelastic materials, Polymer Engineering & Science 9 (1969) 295–310.

548

549

[6] K. C. Valanis, Thermodynamics of large viscoelastic deformations, Studies in Applied Mathematics 45 (1966) 197–212.

550

551

552

[7] W. G. Knauss, I. Emri, Volume change and the nonlinearly thermo-viscoelastic constitution of polymers, Polymer Engineering & Science 27 (1987) 86–100.

553

554

555

[8] J. M. Caruthers, D. B. Adolf, R. S. Chambers, P. Shrikhande, A thermodynamically consistent, nonlinear viscoelastic approach for modeling glassy polymers, Polymer 45 (2004) 4577–4597.
556

557

558

[9] R. A. Simo, On a fully three-dimensional finite-strain viscoelastic damage model: Formulation and computational aspects, Computer Methods in Applied Mechanics and Engineering 60 (1987) 153–173.

559

560

561

[10] G. Holzapfel, R. A. Simo, A new viscoelastic constitutive model for continuous media at finite thermomechanical changes, International Journal of Solids and Structures 33 (1996) 3019–3034.

562

563

564

[11] S. Reese, S. Govindjee, A theory of viscoelasticity and numerical applications, International Journal of Solids and Structures 35 (1998) 3455–3482.

565

566

[12] A. C. Pipkin, Small finite deformations of viscoelastic solids, Reviews of Modern Physics 36 (1964) 1034–1041.

567

568

569

[13] J. Nidhal, Characterization and comparison of viscoelastic properties of three vulcanized rubber materials for a damping application: Natural rubber, bromobutyl and a natural rubber/bromobutyl mixture (????).

570

571

[14] E. Peña, J. A. Peña, M. Doblar, On modelling nonlinear viscoelastic effects in ligaments, Journal of Biomechanics 41 (2008) 2659–2666.

572

573

574

[15] P. O'connell, G. McKenna, Large deformation response of polycarbonate: Time-temperature, time-aging time, and time-strain superposition, Polymer Engineering & Science 37 (1997) 1485–1495.

575

576

577

578

[16] N. W. Tschoegl, W. G. Knauss, I. Emri, The effect of temperature and pressure on the mechanical properties of thermo-and/or piezorheologically simple polymeric materials in thermodynamic equilibrium–a critical review, Mechanics of Time-Dependent Materials 6 (2002) 53–99.

579

580

[17] J. C. Simo, T. J. R. Hughes, Computational inelasticity, volume 7, Springer Science & Business Media, 2006.

581

582

[18] N. W. Tschoegl, Time dependence in material properties: An overview, Mechanics of Time-Dependent Materials 1 (1997) 3–31.

583

584

585

[19] S. Matsuoka, C. J. Aloisio, H. E. Bair, Interpretation of shift of relaxation time with deformation in glassy polymers in terms of excess enthalpy, Journal of Applied Physics 44 (1973) 4265–4268.
586

587

[20] P. J. Flory, Thermodynamic relations for high elastic materials, Transactions of the Faraday Society 57 (1961) 829-838.

588

589

590

[21] G. A. Holzapfel, J. C. Simo, A new viscoelastic constitutive model for continuous media at finite thermomechanical changes, International Journal of Solids and Structures 33 (1996) 3019-3034.

591

592

[22] B. D. Coleman, M. E. Gurtin, Thermodynamics with internal state variables, The Journal of Chemical Physics 47 (1967) 597-613.

593

594

[23] C. Truesdell, W. Noll, The non-linear field theories of mechanics, volume 3, Springer, 2004.

595

596

597

[24] V. Shim, L. Yang, C. Lim, P. Law, A visco-hyperelastic constitutive model to characterize both tensile and compressive behavior of rubber, Journal of Applied Polymer Science 92 (2004) 523-531.

598

599

600

[25] L. Yang, V. Shim, A visco-hyperelastic constitutive description of elastic foam, International Journal of Impact Engineering 30 (2004) 1099-1110.

601

602

603

604

[26] R. Rivlin, Large elastic deformations of isotropic materials. iv. further developments of the general theory, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 241 (1948) 379-397.

605

606

[27] E. Pucci, G. Saccomandi, A note on the gent model for rubber-like materials, Rubber Chemistry and Technology 75 (2002) 839-852.

607

608

[28] R. Ogden, G. Saccomandi, I. Sgura, Fitting hyperelastic models to experimental data, Computational Mechanics 34 (2004) 484-502.

609

610

611

[29] N. Tschoegl, I. Emri, Generating line spectra from experimental responses. part ii: Storage and loss functions, Rheologica Acta 32 (1993) 322-327.

612

613

614

[30] W. G. Knauss, J. Zhao, Improved relaxation time coverage in ramp-strain histories, Mechanics of Time-Dependent Materials 11 (2007) 199-216.
615

616

617

[31] C. Elster, J. Honerkamp, J. Weese, Using regularization methods for the determination of relaxation and retardation spectra of polymeric liquids, Rheologica acta 31 (1992) 161–174.

618

619

[32] M. Nashed, The theory of tikhonov regularization for fredholm equations of the first kind (cw groetsch), SIAM Review 28 (1986) 116–118.

620

621

[33] D. Hibbit, B. Karlsson, P. Sorensen, Abaqus/theory manual, Rhode Island 6 (2007).

622

623

624

[34] S. Hassani, A. A. Soulimani, A. Ehrlacher, A nonlinear viscoelastic model: the pseudo-linear model, European Journal of Mechanics-A/Solids 17 (1998) 567–598.

625

626

[35] R. W. Ogden, G. Saccomandi, I. Sgura, Fitting hyperelastic models to experimental data, Computational Mechanics 34 (2004) 484–502.

627

628

629

630

[36] S. W. Park, R. A. Schapery, Methods of interconversion between linear viscoelastic material functions. part ia numerical method based on prony series, International Journal of Solids and Structures 36 (1999) 1653–1675.

631

632

633

[37] D. Calvetti, S. Morigi, L. Reichel, F. Sgallari, Tikhonov regularization and the l-curve for large discrete ill-posed problems, Journal of Computational and Applied Mathematics 123 (2000) 423–446.

634

635

636

637

[38] M. L. Williams, R. F. Landel, J. D. Ferry, The temperature dependence of relaxation mechanisms in amorphous polymers and other glass-forming liquids, Journal of the American Chemical society 77 (1955) 3701–3707.