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

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.


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 con®gurations [1]. The other class of models re¯ects 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 ®nite linear viscoelasticity (FLV) theory. This method is based on a modi®ed 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 ®rst term of this expansion, while the higher order terms constitute a correction, [6,7]. More recently, several authors have proposed a modi®ed 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 ®nite 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 ®lled rubbers is micromechanically 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 ®nite 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 ®nite 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 ®eld 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 ®rstly 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 ®nally the thermal equations are written in the Lagrangian con®guration.
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 ®nally offered and interpreted in terms of numerical and physical aspects. Consider a body occupying the domain X in the undeformed reference con®guration C 0 and the domain x in the current con®guration C t . The motion of the body from C 0 to C t is locally described by the gradient tensor de®ned as whereX andx de®ne, respectively, the position vector in the reference con®guration and in the current one. The dilatation tensor or right Cauchy±Green tensor is built from the gradient tensor, such that 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 con®guration, [28]. Thus, we consider the multiplicative decomposition where F e de®nes a pseudo-gradient of the elastic motion, and F v a pseudo-gradient of the viscous motion. A viscous and an elastic Cauchy±Green tensors are associated The total Cauchy±Green tensor is then given by

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 where q 0 and q are, respectively, the local mass density in C 0 and C t , P is the ®rst Piola± Kirchoff stress tensor,f corresponds to the speci®c (i.e. per mass unit) body force, e is the speci®c internal energy,Q is the Piola±Kirchoff heat¯ux 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 and W e À Ts ; 9 W is the speci®c free energy, s the speci®c entropy and T the absolute temperature. The ®rst 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 (F; T), [28], another one associated to the viscous behaviour of material: the viscous dilatation tensor So the rate of speci®c free energy can be written as The mechanical dissipation has the following form: On the other hand, the thermodynamics of irreversible processes states the existence of a pseudo-potential of dissipation u such that Constitutive equations

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 W e and W v de®ne the speci®c 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] ; 13 Here, I e 1 , I e 2 denote the ®rst and second principal invariants of C e and I v 1 is the ®rst one of C v . The coef®cients c 1 , c 2 , c 3 and a 1 are material characteristics that must be identi®ed experimentally.
These strain energies are currently used to describe the mechanical behaviour of the hyperelastic incompressible media. According to the multiplicative decomposition (3) and the rheological model, the incompressibility constraint applies both to observable variable F and internal variable C v The total free energy density of the system is then given by The pseudo-potential of dissipation u takes the following quadratic form and only depends on where g is an experimentally identi®ed viscosity parameter.
Substituting (16) and (17) in (11) and (12), and for processes, respectively, isothermally, reversible ( _ T 0 and C v 0) and only isothermal, this choice of thermodynamic formulation leads to the following constitutive laws: In these expressions, the incompressibility constraints (18c) are imposed through the Lagrange multipliers p and q.
Left-multiplying equation (18b) by C v and considering only its deviatoric part, the Lagrange multiplier q disappears and the complementary law becomes

Thermal behaviour
The Fourier law in the reference con®guration gives the Piola±Kirchhoff heat¯ux Q ÀK L Á r L TX; t : 21 In this relation, the conductivity tensor K L is expressed in the initial con®guration and is given in [27] while K is the usual Eulerian conductivity tensor. It is reduced, in the isotropic case, to K K1, where K is the material conductivity factor.

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: 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 P Eq. (18). We ®nally obtain, in Lagrangian description where C e is the heat capacity,Q is the heat¯ux of Piola±Kirchhoff given by the Fourier law (22), D S is the internal source term per volume unit in the initial con®guration, given by 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 ®rst 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.

Variational formulation
The variational formulation of the quasi-static equilibrium problem is written under a perturbed Lagrangian form, [30]. Thus, the displacementũ and the hydrostatic pressure p have to cancel the following integral form (where a is a positive perturbation coef®cient): for all the test functions dũ, dp. The solutionũ; p must also locally satisfy the constitutive equations (18,19) and the incompressibility constraint (15), when the parameter a tends to zero. However, a choice of a 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: for all the test functions dT. 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.

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 satis®es 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.

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 C v and the derivative term oC v =oF. They are given by By means of an implicit Eulerian scheme, [32], the system of differential equations (28) gives on t n ; t n dt a sub-interval of t 0 ; t 0 Dt 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).

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), can be neglected. Thus, the source term D s is reduced to the mechanical dissipation U 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 ®xed-point method on temperature, as illustrated Fig. 2.

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.  The test piece is put into a thermal enclosure at a temperature of 27 C, and subjected to cyclic shearing deformation with a frequency f 4:5 Hz. The amplitude C 0:5 de®nes the relative shearing displacement applied to the external armatures, while the central armature is ®xed.
Two thermocouples, Fig. 5, measure the temperature on one of the elastomeric layers.

Material characteristics
The mechanical properties of the steel armatures are taken as classical values: Young modulus E 210000 MPa, Poisson ratio m 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 identi®cation procedure: (i) By identifying the average curve of a cyclic-shearing response, see Fig. 6a, to a shearing hyperelastic curve (Hart±Smith model), a ®rst 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 g are determined. (iii) Finally, a correction and the complete determination of all parameters is provided according to a least-square ®tting on the hysteresis curve of the shearing test. To do so, the curve is decomposed in to 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].

Numerical simulation
In this section, the ®rst series of computations is performed to determine the mesh-re®nement in¯uence. 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. For the thermal problem (see Fig. 5 and Table 3), linear triangular, bilinear quadrilateral and linear¯ux elements are used to mesh, respectively, the elastomeric layer, the metallic armatures and the boundaries. A zero-¯ux 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¯ux elements), but a quadratic interpolation is used. The boundary conditions are given in Table 3.

Mesh in¯uence
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 Dt 0:89 s followed by a thermal cycle Dt 200 s. Some mechanical and thermal results are presented for the purpose of analysing the mesh re®nement in¯uence, 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 an 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 17.0 30.0 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 in¯uence
In order to analyse the in¯uence of the time partition for the sequential thermo-mechanical coupling, two strategies are checked.
Strategy (i): three times a mechanical cycle Dt 0:89 s followed by a thermal cycle Dt 200 s; Strategy (ii): four times a mechanical cycle Dt 0:89 s followed by a thermal cycle Dt 50 s and two times a mechanical cycle Dt 0:89 s followed by a thermal cycle Dt 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.  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 con®guration (ii) will be adopted. 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 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 to 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 ®rst mechanical cycles. It reaches a stabilised value after a few mechanical periods. This seems to prove that it is suf®cient to consider only three periods for a mechanical cycle.

Mechanical results
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.  To ®t 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 ®lled 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.

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 ®nite transformation hypothesis in the lagrangian con®guration in order to ®t it to the mechanical formulation. The thermo- 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 Fig. 15. Comparison of the computed temperature with the measured temperature in the center of the elastomer 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 ®eld and consequently the evolution of temperature in the elastomeric layer.
The examination of the dissipation ®eld 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 ®eld tends to become uniform.
The computed temperature ®eld 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.