J M Martinez

A Boukamel
email: adnane.boukamel@ec-marseille.fr

S Méo

S Lejeunes

Statistical approach for a hyper-visco-plastic model for filled rubber: experimental characterization and numerical modeling

Keywords: Hyperviscoplastic behavior, Rheological model, Payne effect, Gent-Fletcher effect, Filled elastomer

This paper presents a campaign of experimental tests performed on a silicone elastomer filled with silica particles. These tests were conducted under controlled temperatures (ranging from -55 o C to +70 o C) and under uniaxial tension and in shearing modes. In these two classes of tests, the specimens were subjected to cyclic loading at various deformation rates and amplitudes and relaxation tests at various levels of deformation. A statistical hyper-visco-elasto-plastic model is then presented, which covers a wide loading frequency spectrum and requires indentifying only a few characteristic parameters. The method used to identify these parameters consists in performing several successive partial identifications with a view to reducing the coupling effects between the parameters. Lastly, comparisons between modeling predictions and the experimental data recorded under harmonic loading, confirm the accuracy of the model in a relatively wide frequency range and a large range of deformations.

Introduction

Elastomers belong to the high polymer family i.e. they consist of macromolecular chains of various lengths, with and without ramifications. This structure confers on these materials a low level of rigidity and a high level of deformability. In addition, the reinforcement of these materials with fillers accentuates their dissipative behavior. Because of these properties, especially their damping capacity, these materials are widely used in industry. The application on which this study focuses is that of the drag dampers for helicopters. These parts connect helicopter blades to the rotor and attenuate the drag movement. Designing these parts, which are often related to safety, imposes a guarantee of high reliability under extreme operating conditions (dynamic loading with multi-frequency and large amplitudes, thermal constraints, etc). Meeting these specifications requires good knowledge of the mechanical behavior of the constitutive materials. In addition, the behavior of an elastomer can depend heavily on the temperature, the degree of cross-linking and the type of particles incorporated (carbon or silica), etc.

During recent decades, several approaches have been used by previous authors to model various behavioral aspects of elastomers:

• To describe the static behavior of the material, a hyperelastic approach was used in: [START_REF] Treloar | The elasticity of a network of long chain molecules i[END_REF]Treloar ( , 1957) ) where statistical models were proposed; and [START_REF] Mooney | A theory of large elastic deformation[END_REF]; [START_REF] Rivlin | Some topics in finite elasticity[END_REF]; [START_REF] Hart-Smith | Elasticity parameters for finite deformations of rubber-like materials[END_REF]; [START_REF] Ogden | Large deformation isotropic elasticity, on the correlation of theory and experiment for incompressible rubberlike solids[END_REF], which involved the use of phenomenological approaches.

• Some authors have used the damage mechanics approach to describe the softening behavior occurring during the first loading cycles, which is known as the Mullins effect, see [START_REF] Mullins | Effect of strectching on the properties of rubber[END_REF]; [START_REF] Harwood | Stress softening in rubbers : a review[END_REF]. A theoretical framework was proposed by Govindjee andSimo (1991, 1992) in the case of a hyperelastic behavior. A similar approach was described in [START_REF] Simo | On a fully three-dimensional finite-strain viscoelastic damage model : fomulation and computational aspects[END_REF]; [START_REF] Miehe | Discontinuous and continuous damage evolution in ogden-type large-strain elastic materials[END_REF] in the case of a viscoelastic material.

• In [START_REF] Holzapfel | Fully coupled thermomechanical behaviour of viscoelastic solids treated with finite elements[END_REF]; Holzapfel and Simo (1996a,b); [START_REF] Lion | On the large deformation behaviour of reinforced rubber at different temperatures[END_REF], a thermomechanical coupling model was developed, which takes into account the temperature dependence of the mechanical characteristics and describes the temperature changes resulting from the mechanical dissipation.

• Furthermore, to model the viscoelastic effects of these materials, a framework of the Finite Non Linear Viscoelasticity has been proposed by several authors. These models can be classified in the following groups, depending on the type of formulation used:

Those using an integral approach, which were mainly developed for modeling non linear materials with evanescent memory. These approaches describe the behavior of the material using equations giving the stress tensor in terms of the strain history. [START_REF] Rivlin | Some topics in finite elasticity[END_REF]; [START_REF] Coleman | Thermodynamics of materials with memory[END_REF]; [START_REF] Christensen | Theory of viscoelasticity, an introduction[END_REF]; [START_REF] Coleman | Foundations of linear viscoelasticity[END_REF]; [START_REF] Lianis | Constitutive equations of viscoelastic solids under finite deformation[END_REF]; [START_REF] Chang | The behaviour of rubber-like materials in moderatly large deformations[END_REF]; [START_REF] Morman | Original contributions. an adaptation of finite linear viscoelasticity theory for rubber-like by use of the generalised strain measure[END_REF].

Those using a differential approach, based on the concept of intermediate states commonly used to describe finite elastic-plastic deformations (see [START_REF] Sidoroff | The geometrical concept of intermediate configuration and elastic finite strain[END_REF][START_REF] Sidoroff | Un modèle viscoélastique non linéaire avec configuration intermédiaire[END_REF]). Defining intermediate states provides the internal variables needed to describe the behavior. This approach can be said to be an extension of rheological models in the case of large strains: [START_REF] Sidoroff | Rhéologie non-linéaire et variables internes tensorielles[END_REF]; Le [START_REF] Tallec | Numerical analysis of viscoelastic problems[END_REF]; Le [START_REF] Tallec | Numerical models of steady rolling for non-linear viscoelastic structures in finite deformations[END_REF]; [START_REF] Leonov | On thermodynamics and stability of general maxwell-like viscoelastic constitutive equations, theoretical and applied rheology[END_REF]. The local state method, [START_REF] Lemaître | Mécanique des matériaux solides[END_REF], provides the theoretical framework of this formulation, and the internal variables are provided by the intermediate states.

And those using micro-physically motivated models for filled elastomers, which are often based on hypotheses about the interactions between the agglomerates of fillers and the gum matrix: Drozdov (2001a,b); Drozdov andDorfmann (2002, 2003); [START_REF] Drozdov | Linear thermo-viscoelasticity of isotactic polypropylene[END_REF], or about the mechanisms underlying the deformation and rearrangement of the macro-molecular network: [START_REF] Tanaka | Viscoelastic properties of physically cross-linked networks. transient network theory[END_REF]; [START_REF] Drozdov | A model for the nonlinear viscoelastic response in polymers at finite strains[END_REF][START_REF] Drozdov | A model of cooperative relaxation in finite viscoelasticity of amorphous polymers[END_REF]; [START_REF] Reese | A micromechanically motivated material model for the thermo-viscoelastic material behaviour of rubber-like polymers[END_REF].

In this paper, a meso-physically motivated approach is used to model the response of the material, in large strain and at various frequencies and temperatures. A statistical approach is then proposed to develop a model based on the generalization of an assembly of rheological models. The advantage of this statistical rheological model is that it can be used to simulate the behavior of the material in a wide frequency range while requiring only a few parameters to be identified.

First we present the results of a series of experimental tests, which were carried out on a silicone elastomer filled with silica. These tests were uniaxial tension and shear tests and were performed under controlled temperature (ranging from -55 o C to +70 o C) and under various loading conditions (Relaxation tests, quasistatic and dynamic loading at various strain rates). The results show the dependence of the behavior of the material on the temperature, as well as on the strain-rate (Fletcher-Gent effect, see [START_REF] Fletcher | Non-linearity in the dynamic properties of vulcanised rubber compounds[END_REF]) and the amplitude of the strain (Payne effect, see [START_REF] Harwood | Stress softening in rubbers : a review[END_REF]). The constitutive model is then developed on the basis of the fundamental principle of thermodynamics of continuous media, adapted to finite strain theory. Using the concept of intermediate configurations (multiplicative decomposition of the deformation gradient) and in line with the theory of thermodynamics of irreversible processes, and under the hypothesis of the normal dissipation depending only on the internal variables, the constitutive equation and the flow rules are obtained. A statistical approach is then applied, in order to extend this rheological model to a wide range of strain rates and to account for the plastic behavior of the material. In the following section, this statistical hyper-visco-plastic model is analyzed in the case of simple loads, with a view to propose a strategy for identifying its parameters. For this purpose, analytical solutions are developed to simulate the relaxation response and the hardening test, respectively. In the case of cyclic loading, a semianalytical response was obtained using a symbolic and numeric computation software. These identifications were performed at various temperatures. Lastly, using the semi-analytical solution under sinusoidal shear loading conditions at various frequencies and amplitudes, the effect of parameters such as the temperature, frequency and loading amplitude on the harmonic response of the elastomer are analyzed.

Experimental analysis

Description of the experimental tests

An experimental campaign was conducted on a silicone (dimethyl-vinyl-siloxan vulcanized by peroxide) reinforced with silica particles. The glass transition temperature of this elastomer is approximately -105 o C.

The following tests were carried out:

• Uniaxial tensile tests on specimens with a dumbbell shape (H2 according to standard NF T46-002), to determine the quasi-static behavior and the relaxation response of the material. • Shear tests on Double-Shearing specimens (DS, see Figure 1). These specimens were successively subjected to: a quasi-static loading-unloading cycle; relaxation tests at various shearing levels; triangular cyclic loading, at various strain rates (from 0.03s -1 to 10s -1 ) and various amplitudes (12.5%, 25% and 50%).

All these tests were performed under controlled temperatures (ranging from -55 o C to +70 o C) in a climatic chamber cooled by injecting nitrogen and heated with an electrical resistance and the airflow. In the tensile tests, monitoring and deformation measurements were performed with a laser extensometer.

Remark 1 (Mullins effect) To eliminate the Mullins effect (see [START_REF] Mullins | Effect of strectching on the properties of rubber[END_REF]; [START_REF] Harwood | Stress softening in rubbers : a review[END_REF]) and therefore to characterize the behavior of the stabilized material, a softening process was first induced by applying about ten cycles with an amplitude greater than the maximum strain imposed during the series of tests.

Remark 2 (Temperature stabilization) To avoid errors in the temperature measurement, a waiting period of ten minutes was fixed between each characterization test to allow the temperature to reach equilibrium inside the specimen. The characterization time was sufficiently short to avoid a too strong self-heating phenomena in the specimen.

Experimental results

Relaxation tests: In the relaxation tests, the specimen was subjected to various strain levels: 25%, 50% and 100% under tension loading; 20%, 30% and 54% under shear strain. The response of the material is described by the evolution of the normalized stress 1 versus time. The curves presented in Figure 2(a) and Figure 2(c) show that at temperatures above the ambient temperature, the evolution of these stresses during relaxation was always independent of the strain amplitude. At these temperatures, the relaxation mechanism seems independent of the strain level, under both tension and shear loading; whereas at lower temperatures, the responses doesn't show the same linearity of the stress depending on the deformation, especially in the case of uniaxial tension tests (see Figure 2(b)). The graphs in Figure 2(b) and Figure 2(d) 1 For relaxation tests, normalized stress is obtained by dividing the total stress by the instantaneous stress. show the dependence of the relaxation response on the temperature. The relaxation response was therefore more sensitive to the temperature in the [-55 o C, -25 o C] range than at higher temperatures (above 25 o C).

Quasi-static shear response: The quasi-static test was a loading-unloading test, performed at low strain rate ( γ = 0.03s -1 ) and for three shear amplitudes (γ max = 12.5%, 25% and 50%). The stress-strain curves given in Figure 3 show that even at low rates of deformation, the material shows dissipative behavior. It will therefore be necessary to take the plasticity into account when developing the constitutive model. • Shear rate ( γ): 0.03s -1 , 0.1s -1 , 0.3s -1 , 1s -1 , 3s -1 and 10s -1 .

• Shear amplitude (γ max ): 12.5%, 25% and 50%.

Figures 4(a) and 5 show the effects of the loading amplitude on the stabilized response. Qualitatively, these responses show the strong non-linearity at high amplitudes. In addition, it is worth noting the decrease in the global stiffness observed when the strain amplitude increases. These results therefore clearly confirm that this phenomenon, which is known as the Payne effect, is more pronounced at low temperatures (see Figure 5(a)).

In previous studies, this softening has often been attributed to the plastic behavior of elastomers reinforced with fillers.

As with other visco-elastic polymer materials, the influence of the strain rate was more classical (see Figures 4(b) and6) at all the temperatures tested: an increase in the global stiffness and the cyclic dissipation with the strain rate were clearly observed.

Lastly, the hysteresis loops presented in Figure 7 show the strong influence of the temperature on the behavior of the material. It can be seen in particular that: a softening and a decrease in the cyclic dissipation occurs as the temperature increases; at low temperatures the hysteresis loop shows a nonlinear behavior, which is characterized by: the angular point, the stiffening observed at the end of the cycle and the contraction of the loop at zero strain. These non-linear features, which are more pronounced at low temperatures, are consistent with plastic behavior. These results therefore confirm that plastic behavior begins to predominate when the material approaches the glass transition point.

In conclusion, the following aspects of the behavior have to be taken into account in the model:

• The geometric non-linearities due to large strains. • The dissipative behavior induced by viscous effects which should be coupled to the hyperelasticity.

• The model must be able to describe the behavior in a wide range of strain rates, and in particular, to reflect the Feltcher-Gent effect.

• The effects of plasticity on the behavior, including the Payne effect in particular.

Other phenomena, such as the Mullins effect and the self-heating of the material, were also observed during in this experiments. However, these aspects will not be integrated directly into the model developed in study, because their analysis has been widely discussed in the literature, see [START_REF] Mullins | Effect of strectching on the properties of rubber[END_REF][START_REF] Mullins | Determination of degree of crosslinking in natural rubber vulcanizates. part i[END_REF][START_REF] Mullins | Determination of degree of crosslinking in natural rubber vulcanizates. part iii[END_REF]; [START_REF] Mullins | Stress softening in rubber vulcanizates. part i, use of a strain amplification factor to describe the elastic of filler-reinforced rubber[END_REF]; Govindjee andSimo (1991, 1992).

Constitutive equations

Some generalities about rheological modeling

Using the concept of the local intermediate configuration, introduced by [START_REF] Sidoroff | Un modèle viscoélastique non linéaire avec configuration intermédiaire[END_REF][START_REF] Sidoroff | Variables internes en viscoélasticité, 2. milieux avec configuration intermédiaire[END_REF], the transformation gradient tensor F is split into a viscous and an elastic parts:

F = F e • F v . (1)
Then, assuming that the Clausius-Duhem inequality can be written in Eulerian terms as follows (neglecting the thermal effects):

φ = σ : D -J -1 ρ 0 ψ, (2)
where ρ 0 is the density in the initial configuration, φ is the intrinsic dissipation, σ is the Cauchy stress tensor and D represents the Eulerian strain rate tensor:

D = 1 2 L + L T with L = Ḟ • F -1 , (3)
J denotes the determinant of the gradient tensor F, and ψ is the free specific energy which is expressed as the sum:

ψ(B, B e ) = ψ v (B e ) + ψ 0 (B) (4)
where B = F • F T and B e = F e • F T e are sets of independent thermodynamic variables. So, one can express ψ as follows:

ψ = ∂ψ 0 ∂B : Ḃ + ∂ψ v ∂B e : Ḃe (5)
The time derivatives of the left Cauchy-Green tensor and the local changes in volume are given by:

Ḃ = L • B + B • L T (6) J = J (1 : L) (7)
where 1 is the identity tensor, and:

Ḃe = L • B e + B e • L T -2V e • D o v • V e (8)
where V e is the purely elastic strain tensor (i.e. coming from the polar decomposition F e = V e • R e ) and

D o v = R e • D v • R T e (9)
is the objective measure of the anelastic strain rate. By injecting equations ( 6), ( 7) and ( 8) in ( 5), the variation of the free energy can be written as follows:

ψ = 2B • ∂ψ 0 ∂B : D + 2B e • ∂ψ v ∂B e : D -2V e • ∂ψ v ∂B e • V e : D o v ( 10
)
with the incompressibility conditions:

D : 1 = 0, D o v : 1 = 0 (11)
Using the assumption of the normal dissipation, (we choose a quadratic pseudo-potential of dissipation, ϕ v , depending only on D o v ), the constitutive equation and the evolution law are obtained as follows 2 :

σ = σ 0 + σ v -p1 with              σ 0 = 2ρ 0 J -1 B • ∂ψ 0 ∂B D σ v = 2ρ 0 J -1 B e • ∂ψ v ∂B e D ( 12
)
∂ϕ v ∂D o v = 2ρ 0 J -1 V e • ∂ψ v ∂B e • V e D ( 13
)
where p is the hydrostatic pressure due to the local incompressibility condition. Equations ( 12) and ( 13) can be said to be a generalization of the classical Zener rheological model to the case of finite strain.

A statistical approach for a hyper-visco-plastic model

The constitutive model must first reflect the behavior of the material in a wide range of strain rates, but it also has to account for the effects of the plasticity, such as the Payne effect, and the behavior of the material at low temperatures. Previous studies have shown that plasticity gives good agreement between the 2 the symbol D stands for the deviatoric operator. 2005)), and that some rheological models are suitable for modeling the behavior under a given range of loads [START_REF] Olsson | A fitting procedure for a viscoelastic-elastoplastic material[END_REF]; [START_REF] Miehe | Surimposed finite elastic-viscoelastic-plastoelastic stress response with damage in filled rubbery polymers. experiments, modelling and algorithmic implementation[END_REF]; [START_REF] Nedjar | Frameworks for finite strain viscoelastic-plasticity based on multiplicative decomposition. part i: Continuum formulations[END_REF]). Under more complex loading conditions, rheological models with several branches (see Figure 9(a)) seem to account satisfactorily for the behavior of the material. However, the disadvantage of these models is that they require identifying a large number of parameters.

In order to overcome this difficulty, a statistical approach was developed, whereby the assembly of discrete rheological branches is extended to a continuous model with an infinite number of branches. The advantage of this method is that it covers a wide range of retardation times (or a large frequency spectrum). This approach gives the advantages of a multi-branch assembly without increasing the number of parameters in the model.

The model presented here, can be motivated in micro-physical terms by the heterogeneity of reinforced elastomers, especially in the case of a silicone elastomer slightly filled with silica particles. In order to account for this heterogeneity, it is therefore assumed that the elastomeric matrix is dense and that the inclusions, which consist of particles of silica agglomerated together with a thin rubber bond, are supposed to be slightly reticulated (see for instance [START_REF] Drozdov | A micro-mechanical model for the response of filled elastomers at finite strains[END_REF]). Based on these assumptions, the behavior of the elastomer can be defined as follows:

• The behavior of elastomeric matrix is hyperelastic ;

• The behavior of the inclusions is hyper-visco-elastic (an extension of the Maxwell model to large strains);

• The inclusion/matrix interfaces is assumed to have a hyper-elasto-plastic behavior (an extension of the Saint-Venant model to large strains).

Two statistic quantities are introduced, namely:

• ω i , which denotes the activation energy of the mechanism inclusion/matrix, i.e. the energy required to break the links at the interface [START_REF] Drozdov | A model of cooperative relaxation in finite viscoelasticity of amorphous polymers[END_REF]),

• and P i which represents the probability of that a population of inclusion corresponds to the activation energy ω i .

The discrete form of the statistical model can be written:

           ψ = ψ 0 (B) + N i=1 ψv (ω i , B e (ω i )) P i + ψ p (B ep ) ϕ = N i=1 φv (ω i , D o v (ω i )) P i + ϕ ⋆ p (σ p ) (14)
where D o p denotes the anelastic objective strain rate of the elasto-plastic branch (same definition as in equation ( 9)), ψ 0 denotes the specific free energy associated with the matrix, whereas ψv and ψ p are the free energies associated with the inclusions and the inclusions/matrix interfaces, respectively. φv and ϕ ⋆ p are the pseudo-potential of dissipation 3 corresponding to the inclusions and to the interface, respectively.

Using the continuous statistical model in Figure 9(b) to generalize this formulation, we obtain:

       ψ = ψ 0 (B) + ∞ 0 ψv (ω, B e (ω)) P (ω)dω + ψ p (B ep ) ϕ = ∞ 0 φv (ω, D o v (ω)) P (ω)dω + ϕ ⋆ p (σ p ) ( 15
)
where ω is a random variable associated with the activation energy of a relaxation micromechanism, and P (ω) is the probability that a population of fillers has the given value ω.

Substituting the potentials (15) in equations ( 12) give the constitutive equation of the statistical model 4 :

σ = σ 0 + ∞ 0 σv (ω)P (ω)dω + σ p -p1 with                                σ 0 = 2ρ 0 J -1 B • ∂ψ 0 ∂B D σv (ω) = 2ρ 0 J -1 B e (ω) • ∂ ψv (ω) ∂B e (ω) D σ p = 2ρ 0 J -1 B ep • ∂ψ p ∂B ep D ( 16
)
3 The symbol ⋆ denotes a Legendre-Fenchel transformation. 4 In the hyperelasto-plastic branch, F is split into an elastic part Fep and a plastic part Fpp, (F = Fep • Fpp). The objective measure of the plastic strain rate D o p is therefore defined in the same way as D o v (see ( 9)).

and the following evolution laws:

∂ φv (ω) ∂D o v (ω) = 2ρ 0 J -1 B e (ω) • ∂ ψv (ω) ∂B e (ω) D ( 17
)
D o p = ∂ϕ ⋆ p ∂σ p D ( 18
)
The various forms of the potentials can be chosen so that: the neo-Hookean incompressible hyperelastic model for the matrix behavior; the neo-Hookean form, for the hyperelasticity of the inclusions and the inclusions/matrix interface; a quadratic form for the pseudo-potential of viscous dissipation of the inclusions;

and a perfectly plastic pseudo-potential at the interface. These choice can be written as follows:

                                             ρ 0 ψ 0 = C 1 (I 1 (B) -3) ρ 0 ψv (ω) = G(ω) (I 1 (B e (ω)) -3) φv (ω) = η(ω) 2 D o v : D o v ρ 0 ψ p = A p (I 1 (B ep ) -3) ϕ ⋆ p =< σ p -χ > (19)
where C 1 is the coefficient of the neo-Hookean density, G(ω) and η(ω) are two functions of the random variable ω, χ and A p are the parameters involved in the hyper-elasto-plastic branch, < . > denotes the Mac-Cauley brackets and σ p = √ σ p : σ p . To define the distribution function P (ω), a classical Gaussian distribution centred at the origin and characterized by the standard deviation Ω was adopted:

P (ω) = 1 ∞ 0 P (ω)dω Exp - ω Ω 2 (20)
To choose the functions G(ω) and η(ω), which describe the variations in the hyperelastic and viscous characteristics depending on ω, several forms were tested. The functions giving the best match with the experimental data were of the form:

             G(ω) = G 0 Exp [ω] , η(ω) = η ∞ ln [ √ ω + 1] ω + 1 (21)
In fact, these expressions for the hyperelastic and viscous characteristics lead to a decreasing evolution of the retardation time depending on ω. Combining this variation with the distribution function ( 20) makes it possible to focus on the instantaneous elastic response rather than on the delayed response. Lastly, by injecting ( 19) in ( 16), we obtain the constitutive equation:

σ = σ 0 + ∞ 0 σv (ω)P (ω)dω + σ p -p1 with                      σ 0 = 2C 1 B D σv (ω) = 2G(ω)B e (ω) D σ p = 2A p B ep D ( 22
)
and by using ( 19), ( 17) and ( 18) in ( 8), we obtain the following flow rules:

Ḃe (ω) = L • B e (ω) + B e (ω) • L T -4 G(ω) η(ω) B e (ω) • B e (ω) D (23)
Ḃep = L • B ep + B ep • L T -2 < σ p -χ > σ p σ p • B ep (24)
The statistical hyper-visco-plastic model given by expressions (22) therefore includes 6 parameters which have to be identified, 5 of which are determinist parameters (C 1 , G 0 , η ∞ , A p , χ) and one of which is a statistical parameter (Ω).

Identification of the model parameters

To identify the model parameters by fitting the response of the model to the experimental data, an algorithm implemented in the MATHEMATICA software was used. This algorithm, is based on a minimization of the sum of the squared differences between the experimental data and the analytical or semi-analytical responses. The latter are obtained by simulating the uniaxial tension tests and the double shear tests, in which the responses are assumed to be homogeneous and incompressible.

Analytical forms of tension responses

Under uniaxial tension, the elastic and anelastic gradients of the transformation are written as follows, respectively:

F =      λ 0 0 0 1 √ λ 0 0 0 1 √ λ      , F a =      λ a 0 0 0 1 √ λa 0 0 0 1 √ λa      (25)
Substituting these expressions into the constitutive equations and the complementary law ( 22), the response of the material under various loading modes can be obtained using an analytical form or after a numerical solving.

Relaxation test

To obtain the relaxation response, the instantaneous and delayed stresses are written in terms of the axial component of the first Piola-Kirchoff stress tensor Π, (σ = Π 11 ):

                             σ 0 = 2 λ 3 -1 C 1 + ∞ 0 G(ω)P (ω)dω 1 λ 2 + 3 2 χ λ σ ∞ = 2 C 1 λ 2 λ 3 -1 + 3 2 χ λ σ| t=0 = - 8 3 2λ 3 -λ - 1 λ 3 ∞ 0 G(ω) 2 η(ω) P (ω)dω (26)
where σ 0 is the instantaneous stress response and σ ∞ is the infinite (long-time) stress response.

Hardening test

The axial stress is written here in the quasi-static case, in the form:

σ(λ) = 2 λ C 1 λ λ 3 -1 + A p λ 3 e -1 λ e (27)
The plastic and elastic elongations λ p and λ e are given by:

λp =< λ > H (λ -λ y ) and λ e = λ λ p ( 28
)
where H is a hardening function which is obtained from eq. ( 24) λ y is the elongation corresponding to the plastic yield:

λ y =                    3 1 2 + 1 4 - b 3 3 + 3 1 2 - 1 4 - b 3 3 if b ≤ 3 3 √ 4 2 b 3 Cos 1 3 π -Acos - 1 2 3 b 3 2 if b ≥ 3 3 √ 4 with b = χ 2A p 3 2 (29)
Lastly, the residual strain at zero stress λ 0 corresponds to the equation σ(λ 0 ) = 0 when λ < 0, which can be written as a function of (C 1 , A p , χ) and the maximum strain.

Shear responses

In the case of shear tests, the gradient tensors are taken to be as follows:

F =      1 γ 0 0 1 0 0 0 1      , F a =      λ a1 γ a 0 0 λ a2 0 0 0 1 λa 1 λa 2      (30)

Relaxation test

With the approximation, λ ai = 1, the instantaneous and delayed stress relaxation terms (τ = Π 12 ) are given by:

                             τ 0 = 2γ C 1 + ∞ 0 G(ω)P (ω)dω + √ 2 2 χ τ ∞ = 2C 1 γ + √ 2 2 χ τt=0 = -8γ ∞ 0 G(ω) 2 η(ω) P (ω)dω (31)

Hardening test

The response of the material under quasi-static loading/unloading can be approximated as follows. The shearing stress is written:

τ (γ) = 2C 1 γ + 2A p γ e (32)
The plastic and elastic shear strain γ p and γ e are given by: γp =< γ > H (γγ y ) and

γ e = γ -γ p (33)
where γ y is the shear strain corresponding to the plastic yield:

γ y = 9 + 3χ 2 2A 2 p -3 2 ≃ √ 2χ 4A p (34)

Cyclic test

More generally, based on expressions (30), the response to a cyclic shear test can be obtained by writing the constitutive equations ( 22). This leads to a system of differential equations, the solutions of which are {γ a , λ a1 , λ a2 }. These systems can be solved using a Runge-Kutta scheme.

Identification algorithm

The identification of the parameters of the model X = {C 1 , G 0 , η ∞ , A p , χ, Ω} can be reduced to the minimization of the difference between the experimental curves {(λ i , σ i ), i = 1, N T } and {(γ i , τ i ), i = 1, N S } and the theoretical responses (λ, σ(λ, X)) and (γ, τ (γ, X)). This difference is characterized by the least square distance:

E(X) = NT i=1 ξ i (σ i -σ(λ i , X)) 2 + NS i=1 η i (τ i -τ (γ i , X)) 2 , ( 35
)
where ξ i and η i are the weights associated with the tensile and shear tests, respectively.

This minimization problem is solved using Powell's iterative algorithm (see [START_REF] Fletcher | Pratical methods of optimization[END_REF]), which is a conjugate direction method without gradient calculation. This algorithm is combined with a one-dimensional minimization procedure in each direction which is based on a quadratic interpolation of the function to be minimized.

In short, the identification procedure consists of:

1. determining of the response of the material with a set of model parameters, under a given loading mode;

2. calculating the least square difference between the modeling predictions and the experimental data;

3. applying the iterative procedure to minimize the least square difference.

Identification strategy

Given the complexity of the model, the number of parameters which have to be identified and the multiplicity of the experimental data required to identify these parameters, it was necessary to develop a strategy for decoupling the various stages in the identification procedure. This identification strategy was based on the distinction between the instantaneous and delayed responses, as well as between the effects of viscosity and plasticity. Based on the analytical or semi-analytical results outlined in the previous paragraphs, the following strategy was therefore adopted:

1. Quasi-static and delayed responses are used to identify the parameters C 1 , A p and χ.

2. The other parameters, G 0 , Ω and η ∞ , can be obtained by fitting the values to the instantaneous or cyclic responses at various strain rates.

Table (1) summarizes the successive steps in the identification strategy.

Results and Discussion

Identification results

The six model parameters The table (2) also gives the relative identification error obtained for each temperature and each strain rate. A maximum of 15% of error is obtained and the figure 11 shows two examples of identification results.

(C 1 , G 0 , η ∞ , A p , χ ,
A good agreement between the predictions of the model and the experimental data is observed in the shear tests, at various temperatures and strain rates.

Lastly, Figure 12 shows the evolution of the model parameters identified (normalized) versus the temperature. These curves suggest an exponential decay of all the parameters with the increase of the temperature, except for Ω, which has been fixed at all temperatures, since it characterized only the range of retardation times to be covered by the model. This result shows that the present model seems consistent as the evolution of the parameters between the temperature is monotonous in the temperature range considered.

Relevance of the model

Comparisons between the results of the model predictions and the experimental data (which has not been used for identification), obtained under sinusoidal shear loading conditions, show that the model accurately predicts the effects observed experimentally, namely, the Payne effect (see [START_REF] Harwood | Stress softening in rubbers : a review[END_REF]), as shown in Figure 14(c), and The Gent-Fletcher effect (see [START_REF] Fletcher | Non-linearity in the dynamic properties of vulcanised rubber compounds[END_REF]) as shown in Figure 14(a).

In addition, figures 14(b) and 14(d) show the existence of good agreement between the simulated and the experimental data on the cyclic dissipation depending on the frequency and the amplitude.

Other comparisons made in multi-harmonic loading in shear tests also shows the ability of this model to accurately simulate the behavior of materials subjected to a combination of several loads at different frequencies (see figure ( 13)). These results show that the present model can successfully predict the complex behavior of a highly dissipative silicone rubber for a large range of strain amplitudes and strain rates with a few number of material parameters.

Conclusion

In this study, a statistical approach was used to develop a hyper-visco-plastic model covering a wide frequency range. This approach has the same advantages of classical multi-branch models, such as the ability to simulate the behavior of material for several decades of retardation time, but it do not show the same inconvenient as only a few number of material parameters are required (6 in the present model).

A series of experiments were conducted under various loading and temperature conditions, in order to identify the parameters involved in the model, using an algorithm developed with the MATHEMATICA software library. To optimize this identification procedure, a relevant strategy was adopted, which consisted in distinguishing between the various stages in the procedure and thus reducing the effects of coupling

Figure 1 :

 1 Figure 1: Double-shearing specimen.

  Effects of the temperature under shear loading.

Figure 2 :

 2 Figure 2: Response of the material in relaxation tests: first Piola-Kirchoff stress vs time. The vertical axis corresponds to the total stress divided by the instantaneous stress (normalized stress).

Figure 3 :

 3 Figure 3: Quasi-static responses recorded in loading-unloading shear tests: first Piola-Kirchoff stress vs shear strain (shear rate: γ = 0.03s -1 , temperature: T = 25 o C).

  Stabilized cycles (T = 25 o C, γ = 0.3s -1 ), effect of the strain amplitude. Stabilized cycles (T = 25 o C, γmax = 50%), effect of the strain rate.

Figure 4 :-

 4 Figure 4: Cyclic responses in triangular shear tests: first Piola-Kirchoff stress vs shear strain (temperature: T = 25 o C). The first cycles has been removed to keep only stabilized responses.

Figure 7 :

 7 Figure 6: Effects of the strain rate on the response of the material under triangular cyclic shearing tests at various temperatures (γmax = 50%).

Figure 8 :

 8 Figure 8: Intermediate configuration.

Figure 9 :

 9 Figure 9: Statistical hyper-visco-plastic model

  Evolution of the distribution function P(ω).

Figure 10 :

 10 Figure 10: Evolution of statistical functions: predominance of the instantaneous elasticity.

  Ω) were identified successively, at various temperatures (-55 o C, -40 o C, -25 o C, 25 o C, 40 o C and 70 o C) and at various strain rates (3s -1 , 10s-1 ). The values obtained are given in table (2).

Figure 11 :

 11 Figure 11: Comparison between stabilized hysteretic cycles at various temperatures and strain rates. The solid line show the model response, the points show the experimental results.

  Figure 12: Evolution of parameters vs temperature.

Figure 14 :

 14 Figure 14: Comparisons between modeling predictions (solid line) and experimental data (points).

Table 1 :

 1 Identification strategy: C 1 , Ap and χ are identified on quasi-static and delayed responses, G 0 , Ω and η∞ are obtained by fitting the values to the instantaneous or cyclic responses at various strain rates.

	Error (%)

Table 2 :

 2 Parameters identified at various temperatures.

		0.3
		0.2
	Stress (MPa)	-0.1 0 0.1
		-0.2
		-0.3
		-0.2	-0.1	0	0.1	0.2
				Strain

Treloar, L., 1957. The present status of the theory of large elastic deformations. In: The rheology of elastomers.

between the parameters.

The results obtained at various temperatures show the ability of this model to simulate the behavior of the material in a wide range of temperatures. In addition, the comparisons between the modeling predictions and the experimental data recorded at various frequencies and strain amplitudes have shown a good agreement.

The present model is therefore capable to reproduce the complex behavior of filled rubber in particular the Gent-Fletcher and Payne effects.


T A N'guyen

S Lejeunes
email: lejeunes@lma.cnrs-mrs.fr

D Eyheramendy

A Boukamel

A thermodynamical framework for the thermo-chemo-mechanical couplings in soft materials at finite strain

Keywords: polymer, rubber, thermodynamics, irreversible processes

This paper focus on the modeling of thermo-chemo-mechanical behaviors for soft materials (polymer, rubber, ...) that can exhibit large strain in a rigorous thermodynamical framework. As an application it is considered the case of a thermo-viscoelastic material that can undergo chemical reactions which are described by a kinectic approach. Simple idealized and homogeneous numerical tests are considered and illustrate the availability of this framework to take into account reciprocal coupling between each physics.

Introduction

Thermo-chemo-mechanical couplings in soft materials (polymers, rubbers, ...) are an important topic that is linked to many different industrial applications. The most evident one concerns the processing of these materials. A precise simulation of material processing can be of great interest to optimize the process, to predict residual strain/stress due to process, to anticipate material heterogeneities (different grade of material due to different chemical states). Another application concerns the prediction of the behavior of industrial parts submitted to severe operating conditions (in aeronautics or spatial applications for instance). In these applications: thermal and chemical aging can be coupled to mechanical damage phenomena. Concerning modelisation, the thermodynamic of irreversible processes is the main basic tool to address these problematics. The pioneer work of Prigogine and co-authors is very important as it allows the interpretation and the modelisation of chemical reactions in the context of irreversible processes (see for instance [START_REF] Prigogine | Introduction to thermodynamics of irresversible processes[END_REF]). Prigogine has considered that a thermodynamical system can be described by classical state variables (temperature, volume,...) and by supplementary ones which are the extend of reactions in the case of reactive systems. Thermodynamical energy potential (Gibbs or Helmholtz free energy) are then dependent on these new thermodynamical variables. Chemical thermodynamical forces can be defined, called affinity, and linked to earlier work of De Donder (1928). The second law of thermodynamic (entropy production) can be satisfied as these chemical thermodynamical forces derive from the thermodynamical potential of energy. The generalisation of these fundamental concepts to the case of continuum mechanic can be done in a straightforward manner by considering the concept of local state and the introduction of internal variables (see [START_REF] Germain | Continuum thermodynamics[END_REF] and reference therein). Therefore all the thermodynamical variables are function of time and position. Applying Onsager reciprocity concept (see [START_REF] Onsager | Reciprocal Relations in Irreversible Processes[END_REF]) or using the convex framework of standard generalized materials (see [START_REF] Halphen | On generalized standard materials[END_REF]) evolution equations for chemistry or mechanical irreversibilities can be established from the definition of a thermodynamical energy potential and eventually from a pseudo potential of dissipation (in the context of standard generalized materials). These general tools has been used by many authors in the literature in the context of soft material subjected to chemical reactions. [START_REF] Gigliotti | Chemo-mechanics couplings in polymer matrix materials exposed to thermo-oxidative environments[END_REF]; [START_REF] Gigliotti | Assessment of chemo-mechanical couplings in polymer matrix materials exposed to thermo-oxidative environments at high temperatures and under tensile loadings[END_REF] has proposed a chemo-mechanical coupled model in the context of elastic behavior with several reactions for a thermo-oxydation process in polymer matrix. They also consider the diffusion of species inside the materials to predict aging. In [START_REF] Lion | On the phenomenological representation of curing phenomena in continuum mechanics[END_REF]; [START_REF] Mahnken | Thermodynamic consistent modeling of polymer curing coupled to viscoelasticity at large strains[END_REF] one can find thermo-chemo-mechanical models for the curing of polymer. Theses models take into account of mechanical, thermal and chemical deformations (dilatation and shrinkage) in a thermo-viscoelastic context. The chemistry is phenomenologically described by a kinetic approach. [START_REF] André | Thermo-mechanical behaviour of rubber materials during vulcanization[END_REF] has proposed a thermo-chemo-mechanical model in small strain to simulate the vulcanization of rubber materials. The mechanical behavior is assumed to be elasto-visco-plastic and two phenomenological chemical reactions are considered to describe the vulcanization process. In [START_REF] Kannan | A thermodynamical framework for chemically reacting systems[END_REF], the authors have developed a framework in finite strain that is applied to vulcanization of rubbers in a viscoelastic context with several chemical reactions. In this contribution, we propose a rigorous thermodynamical framework that can be used as a basis of development for thermo-chemo-mechanical models. As previously mentioned the basic tools are the thermodynamics of irreversible processes and the local state hypothesis. As already done by some authors, we assume that the volume variation is decomposed in a thermal part (dilation) a chemical part (shrinkage) and a mechanical part (compressibility). Hydrostatic pressure is therefore directly related to the chemical, mechanical and thermal states. This decomposition also implies a coupling of the heat capacity and the chemical evolution upon pressure or volume variation. A kinetic (phenomenological) description of chemistry is chosen. However, due to the fact that mechanical parameters (modulus) and volume variation are dependent on the chemical state, the thermodynamical chemical force is also directly dependent on the mechanical state. The hydrostatic pressure and deviatoric strain can therefore have a favorable or a unfavorable effect on chemical evolution. The main originality of this paper lies in the following points. First, it is proposed a new form of chemical potential of energy, this potential takes into account of an initiation temperature under which no reactions take place and energy is stored as heat. This potential is defined so as to have a clear definition and admissible form of the heat capacity. Second, the chemical evolution that is associated to this potential naturally takes into account of the influence of the mechanical state. It is proposed a new evolution law for chemistry that is inspired from visco-plasticity behavior: the chemical evolution rate is assumed to be null if the thermodynamical force becomes negative (therefore chemical reversion is not allowed) and a chemical viscosity parameter is introduced. This parameter allows the balancing of mechanical influence on chemical evolution. The paper is organised as follows: the general thermodynamical framework and conservation equations are presented in a first section. As an application it is derived a phenomenological thermo-chemo-visco-elastic model in a second section. This model is general and can be applied to different problematic. In a third section, some simple idealized examples are provided. These examples illustrate both the capability of the proposed model to be applied to different problematic (aging or material processing) and the influence of some material parameters on the reciprocal couplings. In the last section some remarks are addressed to conclude.

Thermodynamical framework

Definition of a chemical internal variable

In this work, it is considered a solid body submitted to thermo-mechanical loadings. This body is viewed as a closed 1 thermodynamical continuum system and therefore no mass exchange with the outside can occur. Furthermore as classically assumed neither creation nor destruction of mass is allowed. Let us assumed that the initial configuration (at t = 0) is defined by a stress and heat flux free state. In this configuration the chemical species are in equilibrium: the reaction rates are null in this state. Using the classical continum mechanics point of view, the body can be considered as a continuum of material points that are defined by an initial position X and a current position x in an euclidean space. The material point can be viewed as an infinitesimal element of volume and it can be defined the current material density ρ(x, t), which evolves in times as the volume evolves (the initial density is denoted ρ 0 ). The closed system hypothesis leads to:

ρ + ρdiv(v) = 0 (1)
where div(v) is the eulerian divergence of the velocity of the material point and ρ is the total time derivative (so called material time derivative2 ) of the density. Assuming that this infinitesimal element of volume is a mixture of all chemical species that compose the material, one can defined the current mass concentration of the i th species relative to the current infinitesimal volume (of the mixture), denoted ρ i (x, t). In the case of null diffusion of species inside the body (i.e. no relative velocity of one species from another) the mass conservation is also defined from:

     ρi + ρ i div(v) = m i i = 1, 2, .., n n i=1 m i = 0 (2)
where m i (x, t) is the rate of production of mass of the i th species per unit of current volume of mixture.

Knowing the molar mass of each species, denoted M i equations 2 can be rewritten in terms of volume concentration (mole per unit of current volume) denoted

Y i = ρ i /M i :        Ẏi + Y i div(v) = m i M i i = 1, 2, .., n n i=1 m i = 0 (3)
As already mentioned by [START_REF] Kannan | A thermodynamical framework for chemically reacting systems[END_REF] in finite strain it is not a good idea to consider Y i (x, t) as thermodynamical variables that characterize the chemical state because volume variation leads to a change in concentration (in terms of volume). Following aforementioned authors, we prefer to consider n i the number of mole of the ith species per unit of mass of the mixture:

n i (x, t) = Y i /ρ = ρ i /(ρM i ).
Inserting this definition in 3 and using 1, the mass balance becomes:

n i=1 ṅi M i = 0 (4)
Therefore only n -1 variables are sufficient to describe fully the chemical state. If the reactive scheme is known and if this reactive scheme implies r stoichiometric reactions (with the stochiometric ratios ν ir ), the extend ζ r of the rth reaction can be defined from:

ζ r (x, t) = n r i (x, t) -n r i (X, 0) ν ir ∀i, r = 1..m (5)
where n r i (x, t) is the molar concentration of the ith species for the rth reaction per unit of mass of the mixture. Finally if both the initial concentration and the stochiometric coefficients are known the chemical state can be defined by only r independant variables that describe the advance of reactions. In this paper, it is only considered one reaction3 and the chemical state is defined by the following normalized concentration (chemical conversion):

ξ(x, t) = n 1 vu (x, t) -n 1 vu (X, 0) n 1 vumax (6)
where n 1 vu (x, t) is the current molar concentration of the product of interest (in the case of material processing this could be the vulcanized or polymerized product) and n 1 vumax is the maximum attainable concentration. The previously defined chemical conversion is chosen to be a thermodynamical internal variable in the following.

Kinematic

To define kinematic, it is assumed that the motion of all the material points of the body is described by a bijective vectorial function χ(X, t) such that: x = χ(X, t). The transformation gradient F(X, t) is therefore defined by F = ∂χ/∂X. Following [START_REF] Flory | Thermodynamic relations for highly elastic materials[END_REF], the transformation is split into volumetric and isochoric part:

F = (J 1/3 1) • F (7)
where J = detF is the volume variation and 1 the identity tensor. It is assumed that the volume variation can result form three independent terms: the mechanical compressibility J m , the thermal dilatation J Θ and the chemical shrinkage

J ξ : J = J m J Θ J ξ (8)
The volume variation terms are defined from the following relations:

J Θ = 1 + α(Θ -Θ 0 ) (9)
J ξ = 1 + βg(ξ) (10)
J m = J(J ξ J Θ ) -1 (11)
where α and β are respectively the thermal expansion and chemical shrinkage coefficients, Θ 0 is the temperature in the initial configuration, g is a shrinkage function. For the sake of simplicity, the thermal dilatation is assumed as linear (eq. 9). More complex expressions for the thermal dilatation can be considered that depend on the material composition and on the temperature range of the desired application (see for instance [START_REF] Holzapfel | Entropy elasticity of isotropic rubber-like solids at finite strains[END_REF] who have proposed an exponential form). The linearized form of the thermal dilatation is commonly assumed in the case of small thermal variations. However for some materials this hypothesis seems valuable for a large range of temperatures (far from the glass transition temperature). For instance, figure 1 show the evolution of density of a partially hydrogenated poly(acrylonitrileco-1,3-butadiene) upon temperature for various density of cross-links. To describe inelastic effects an intermediate configuration is introduced and the isochoric transformation gradient is classically decomposed with the following multiplicative split:

F = F e • F i (12)
where F i represents the inelastic transformation and F e represents elastic one. This decomposition implicitly assumes that inelastic flows are incompressible.

Clausius-Duhem inequality

It is chosen in the following to consider the Helmoltz free specific energy to characterize thermodynamical states. This potential of energy is classical adopted by mechanician whereas chemists are more familiar with Gibbs potential. It is therefore defined ψ the Helmoltz free specific energy:

ψ = e -Θs ( 13
)
where e is the specific internal energy and s is the specific entropy. Combining the first and second laws of thermodynamics together with the previous expression, one can obtain the Clausius-Duhem inequality. In an Eulerian configuration, it can be written as:

φ = σ : D -ρ ψ -ρs Θ -q • grad x Θ Θ ≥ 0 ∀ D, Θ, q (14)
where φ stands for the dissipation, σ(x, t) is the Cauchy stress, D = Ḟ•F -1 the eulerian rate of deformation, q(x, t) is the eulerian heat flux, grad x is the eulerian gradient operator. Introducing the left Cauchy-Green deformation B = FF T and Be = F e F e T and assuming that the free energy is a function of B, Θ, J, Be , ξ, it can be obtained:

ψ = ∂ψ ∂ B : Ḃ + ∂ψ ∂ Be : Ḃe + ∂ψ ∂J J + ∂ψ ∂Θ Θ + ∂ψ ∂ξ ξ ( 15
)
where:

J = J (1 : D) (16)
Ḃ = L • B + B • L T - 2 3 (1 : D) B (17)
Ḃe = L • Be + Be • L T -2 Ve • Do i • Ve - 2 3 (1 : D) Be ( 18
)
with Ve is the pure deformation coming from the polar decomposition of F e = Ve • Re , Do i is the objective rate of inelastic deformation, defined from:

Do i = Re • Di • RT e .
Inserting equations ( 16), ( 17), ( 18) in ( 15) and putting the result in ( 14) and regrouping terms, we have:

φ = σ -2ρ B • ∂ψ ∂ B + Be • ∂ψ ∂ Be D -ρJ ∂ψ ∂J 1 : D + 2ρ Be • ∂ψ Be D : Do i -ρ s + ∂ψ ∂Θ Θ -ρ ∂ψ ∂ξ ξ - grad x Θ Θ • q ≥ 0 (19)
To proceed further, it is made the following hypothesis: entropy is fully defined from the free specific energy variation and there are no dissipation for the thermodynamical force associated with the thermodynamical flux D. This leads to:

s = - ∂ψ ∂Θ (20) σ = 2ρ B • ∂ψ ∂ B D σeq + 2ρ Be • ∂ψ ∂ Be D σneq + ρJ ∂ψ ∂J 1 σ vol (21)
It remains the following terms in the dissipation:

φ = φm 2ρ Be • ∂ψ ∂ Be D : Do i + φ ξ -ρ ∂ψ ∂ξ ξ + φΘ - grad x Θ Θ • q ≥ 0 (22)
where φ m is the intrinsic dissipation, φ ξ the chemical dissipation and φ θ the thermal dissipation. It is assumed that φ m , φ ξ and φ Θ are independently positives. Interesting readers can refer to [START_REF] Germain | Continuum thermodynamics[END_REF] and references therein for a general discussion on thermodynamics of local state models. From eq. ( 22), it can be defined a mechanical thermodynamical force σ i = 2ρ( Be • ∂ψ/∂ Be ) D , a chemical force A ξ = -ρ∂ψ/∂ξ and a heat force A Θ = -grad x Θ/Θ.

Heat equation

The heat equation can be obtained from the first thermodynamical principle (energy conservation), which takes the following local form in the eulerian configuration:

ρ ė = σ : D + ρr -div x q (23)
Where r is a volumetric heating source term (defined by unit of volume). The material time derivative of eq. ( 13) leads to: ė = ψ + s Θ + ṡΘ (24) using eq (24) in eq (23), one can obtain:

ρ ṡΘ = σ : D + ρr -div x q -ρ ψ -ρs Θ (25)
The material time derivative of the free energy is given by eq (15) and using eqs ( 16), ( 17), ( 18), ( 20) in ( 25):

ρ ṡΘ = φ m + φ c + ρr -div x q (26)
The material time derivative of the entropy is obtained from eq (20):

ṡ = - ∂ 2 ψ ∂θ 2 Θ- ∂ 2 ψ ∂θ∂ξ ξ -2 B • ∂ 2 ψ ∂Θ∂ B + Be • ∂ 2 ψ ∂Θ∂ Be D : D+2 Be ∂ 2 ψ ∂Θ∂ Be D : Do i -J ∂ 2 ψ ∂Θ∂J (1 : D) (27)
Finally, the heat equation is obtained by replacing eq (27) in eq (26):

ρC Θ = φ m + φ c + l m + l c + ρr -div x q (28)
where C is the specific heat capacity which is defined from:

C = -Θ ∂ 2 ψ ∂θ 2 (29)
the coupling terms l m and l c are defined from:

l m = Θ ∂σ ∂Θ : D -Θ ∂σ neq ∂Θ : Do i (30) l c = -Θ ∂A ξ ∂Θ ξ (31)

Application: a thermo-chemo-visco phenomenological model

One can consider the specific free energy splitting into mechanical, thermal and chemical parts as:

ψ = ψ m ( B, Be , ξ, Θ, J) + ψ Θ (Θ) + ψ ξ (Θ, ξ) (32)
This decomposition is based on the following ideas: (i) the mechanical behavior is strongly dependent on chemical and thermal state even if the behavior is elastic, (ii) without mechanical deformations or mechanical stresses the chemical free energy will mainly depend on the chemical state and on the temperature (classical thermo-kinetic approach), (iii) in a purely thermal process free energy depends only on temperature. Using the previous decomposition leads to an additive splitting of the entropy and eventually to the specific heat capacity splitting:

s = - ∂ψ ∂Θ = s m + s Θ + s ξ C = -Θ ∂ 2 ψ ∂Θ 2 = C m + C Θ + C ξ (33)
PSfrag replacements

ψ vol (J, Θ, ξ) C10(Θ, ξ), C01(Θ, ξ) G(Θ, ξ) η(Θ, ξ)

mechanical part

Without loss of generality it is chosen the very simple Zener viscoelastic model for the mechanical behavior. The model is schematized on the figure 2. The mechanical free energy is decomposed into an equilibrium, a non-equilibrium and a volumetric part:

ψ m ( B, Be , ξ, Θ, J) = ψ eq ( B, ξ, Θ) + ψ neq ( Be , ξ, Θ) + ψ vol (J, ξ, Θ) (34)
The equilibrium and non-equilibrium free energy are based on Mooney-Rivlin and neo-Hookean hyperelastic model:

ρ 0 ψ eq = C 10 (Θ, ξ)(I 1 ( B) -3) + C 01 (Θ, ξ)(I 2 ( B) -3) (35) ρ 0 ψ neq = G(Θ, ξ)(I 1 ( Be ) -3) (36)
where I 1 (•) and I 2 (•) are the first and second invariant of the tensorial argument (•), defined by

I 1 (•) = tr(•) and I 2 (•) = 1/2(I 1 (•) 2 -tr(• 2 )).
The material parameter C 10 , C 01 and G, are assumed to depend both on temperature and chemical state. The volumetric part is assumed to be related only to compressibility and the compressibility modulus K v is assumed to be independent of the thermal and chemical state:

ρ 0 ψ vol = K v 2 (J m -1) 2 = K v 2 ( J J Θ J ξ -1) 2 (37)
The stresses are therefore expressed as:

σ eq = 2C 10 (Θ, ξ)J -1 BD + 2C 01 (Θ, ξ)J -1 (I 1 ( B) B -B2 ) D σ neq = 2G(Θ, ξ)J -1 BD e σ vol = K v (J m -1)J -1 Θ J -1 ξ 1 = p1 (38)
The visco-elastic behavior is assumed to be described by the following flow rule:

Do i = 2 η(Θ, ξ) Be • ∂ψ ∂ Be D = 2G(Θ, ξ) η(Θ, ξ) BD e ( 39
)
where η(Θ, ξ) is a viscosity parameter. Using (39) in [START_REF] Onsager | Reciprocal Relations in Irreversible Processes[END_REF], it is obtained the following evolution equation (classical Maxwell viscosity):

Ḃe = L • Be + Be • L T - 2 3 (1 : L) Be - 1 τ (Θ, ξ) BD e • Be ( 40
)
where τ (Θ, ξ) = η(Θ, ξ)/4G(Θ, ξ) is a characteristic time of viscosity. If the dependency of C 01 , C 10 upon temperature is assumed to be linear, the mechanical contribution to the heat capacity is only given by:

ρ 0 C m = ΘK v α 2 (2 -3J m )JJ -3 Θ J -1 ξ -Θ ∂ 2 G ∂Θ 2 (I 1 ( Be ) -3) (41)
This contribution can be either positive or negative depending on the value of J m and the sign of the second derivative of the modulus G. However, C m is usually small compared to C Θ .

Thermal part

For the thermal part of the free energy it is adopted the form proposed by Reese and Govindjee (1997); [START_REF] Behnke | An extended tube model for thermo-viscoelasticity of rubberlike materials: Parameter identification and examples[END_REF] which leads to a linear dependency of heat capacity upon temperature:

ρ 0 ψ Θ (Θ) = C 0 Θ -Θ 0 -Θ log Θ Θ 0 -C 1 (Θ -Θ 0 ) 2 2Θ 0 (42)
An isotropic Fourier transport is assumed for the heat flux in the eulerian configuration:

q = -k t grad x Θ ( 43
)
where k t is the thermal conductivity coefficient and it is neglected the dependence of this coefficient upon temperature and chemical state. One can easily verify that the previous heat flux is thermodynamically consistent, ie. the mechanical dissipation is always positive independently of the thermodynamic state.

The contribution of the thermal part to the heat capacity is therefore:

ρ 0 C Θ = C 0 + C 1 Θ Θ 0 (44)

Chemical part

The choice of a free energy for a chemical reactive system within a rigorous thermodynamical framework is not obvious. For well established reactive systems (when stochiometric equations can be written with a clear reactive scheme), the Gibbs free energy can be written as a function of chemical potentials of each species. For perfect or ideal systems the chemical potentials have simple form which can depend on temperature and hydrostatic pressure. This approach has been used by [START_REF] Kannan | A thermodynamical framework for chemically reacting systems[END_REF] for the vulcanization of natural rubber. The authors consider a complex reactive scheme in isothermal conditions. In [START_REF] Mahnken | Thermodynamic consistent modeling of polymer curing coupled to viscoelasticity at large strains[END_REF], it is proposed to formulate the chemical part of the free energy as a linear function of the chemical conversion and the amount of heat generated during curing. This amount of heat depends on the all history during processing and it is considered as a material parameter in the aforementioned study.

In this work it is proposed a different approach: at constant temperature the free energy has to be the integral of the chemical force for an increment of the chemical conversion. As a kinetic approach is assumed, the chemical force is multiplicatively split into a function of the conversion (ξ) and a function of the temperature. The conversion part is inspired by the work of [START_REF] Prime | Differential scanning calorimetry of the epoxy cure reaction[END_REF]. The temperature part is motivated by the fact that in general it exists a temperature of activation of chemical reactions. Below this temperature no reactions occur and heat is stored by the material. The temperature part also control the thermochemical coupling term in the heat transport equation and eventually the dependency of the heat capacity to the chemical conversion. This two contributions are rarely discussed in the literature but wrong choices may lead to non-physical models. It is proposed to use the following form for the chemical part of the free energy:

ρ 0 ψ ξ (Θ, ξ) = C 2 Θ i log Θ Θ i (1 -ξ) m+1 m + 1 -Θ 0 log Θ 0 Θ i ( 45
)
where C 2 and m are material parameter, Θ i is an initiation temperature: below this temperature at a null hydrostatic pressure (pressure of free stress state in the present model) no reaction can occur. The chemical force is therefore given by:

A ξ = -ρ 0 J -1 ∂ψ m ∂ξ -ρ 0 J -1 ∂ψ ξ ∂ξ ρ 0 ∂ψ ξ ∂ξ = -C 2 Θ i log Θ Θ i (1 -ξ) m ρ 0 ∂ψ m ∂ξ = ∂C 10 ∂ξ (I 1 ( B) -3) + ∂C 01 ∂ξ (I 2 ( B) -3) + ∂G ∂ξ (I 1 ( Be ) -3) -β ∂g ∂ξ J -1 ξ p (46)
The chemical force is therefore coupled to the mechanical state through the dependency of the mechanical parameter upon ξ and through the hydrostatic pressure. Furthermore the sign of this coupling determine the influence of the mechanics upon chemistry a given state of strain and stress can be favorable or unfavorable (a hydrostatic compression state is a priori favorable to a vulcanization process for rubbers or polymerization for polymers).

It is assumed the following evolution equation:

ξ = k(Θ) A ξ ( 47
)
where < . > are the Mac-Cauley brackets4 , k(Θ) is a kinetic rate term which is here assumed to follow an Arrhenius law:

k(Θ) = A exp -Ea RΘ ( 48
)
where A is a material parameter, Ea is an activation energy and R is the ideal gaz constant (R = 8.314J/mol/K). The Mac-Cauley brackets are introduced to guarantee a null or positive chemical evolution5 .

If one consider a deformation for which I1 = I2 = 3, the kinetic evolution is therefore defined by:

ξ = A J exp -Ea RΘ C 2 Θ i log Θ Θ i (1 -ξ) m + β ∂g ∂ξ J -1 ξ p (49)
It can be seen from the previous equation that if Θ < Θ i the first term inside the Mac-Cauley Brackets is negative and therefore no reactions can occur if the second term is null (null hydrostatic pressure) or negative (depending on the sign of ∂g/∂ξ and the sign of p). The relative influence of mechanic under chemical reaction is parametrized by the product C 2 × A and β × A.

Inserting (47) in the expression of φ ξ leads to a straight forward proof of the admissibility of the previous chemical flow rule:

φ ξ = k(Θ) A ξ A ξ ≥ 0 ∀ξ, Θ, F (50)
The contribution of the chemical free energy to the heat capacity is obtained as follows:

ρ 0 C ξ = -Θ ∂ 2 ψ ξ ∂Θ 2 = C 2 Θ i Θ (1 -ξ) m+1 m + 1 (51)
This chemical heat capacity is a positive decreasing function of ξ. It can be seen that the proposed form of the chemical free energy is a priori compatible with the idea of a heat energy storage to be released during chemical reactions. In the literature, it can be found that the impact of the chemical evolution upon the specific heat capacity is greatly dependent on the material and the chemical process considered. In [START_REF] Likozar | Cross-linking of polymers: Kinetics and transport phenomena[END_REF], authors show that the heat capacity of a HPAB polymer has mainly a linear dependency upon temperature and a decreasing behavior upon the formation of chemical bonds in the polymer (this effect is important: 25% of relative influence on the heat capacity from the uncured to the fully cured state).

In Dimier ( 2003) the vulcanisation of a natural rubber with carbon black fillers is considered, the authors show that the heat capacity evolution upon chemical state is very low and could be neglected (for a Natural Rubber with carbon black).

The shrinkage function g(ξ) must fulfill some restrictions: this function has to be chosen such as -1 < g(ξ) ≤ 0 ∀ξ and g(ξ = 0) = 0 to guaranty that 0 < J ξ ≤ 1. In the present model the following form has been adopted:

g(ξ) = 0.367 -exp -(1-ξ) n 0.632 (52)
where n is a material parameter that controls the chemical shrinkage function as well as the material parameter β that is the shrinkage coefficient (see eq. ( 3)). Figure 3 shows the influence of n on the shrinkage function. It is also assumed in this work that the chemical rate ξ, is identically null when ξ = 1, therefore g (ξ = 1) must be null. This condition is fulfilled if n > 1.

-

Numerical simulation for simple homogeneous tests

The material parameters are synthesized in table 1. These parameters have not been identified from a real material however the order of magnitude of most of them are closed from those of a filled rubber material. For the mechanical parameters, it has been adopted the following hypothesis:

• the chemical conversion is assumed to stiffen the mechanical behavior

• far from the glass transition temperature, an increase of the temperature leads to a stiffening effect for the elastic part (parameter C 10 ) of the model (entropic elasticity) and to a softening effect for the viscoelastic part (parameter G).

• the characteristic time of viscosity: η/(4G), is a decreasing function of the temperature.

Exemple 1: cyclic shearing at fixed temperature and fixed chemical conversion

To illustrate the influence of chemical conversion and temperature on the mechanical response, it is considered a cyclic simple shear test. The response is assumed as homogeneous and the deformation gradient is stated as following (thermal dilatation and chemical shrinkage are neglected):

F(t) =   1 γ(t) 0 0 1 0 0 0 1   γ(t) = g sin(2Πf t) (53)
where f = 1Hz is the chosen frequency and g is the amplitude of the signal. It is also assumed the following form for the viscoelastic Cauchy-Green gradient:

Be (t) =   B e11 (t) B e12 (t) 0 B e12 (t) B e22 (t) 0 0 0 1/(B e11 (t)B e22 (t))   (54) Density ρ 0 (Kg/m 3 ) 1000 Thermal part α(K -1 ) C 0 (J/m 3 /K) C 1 (J/m 3 /K) k t (W/m/K) 2.2e -4
8e 5 1e 6 0.22 Mechanical part K v (P a) C 10 (P a) C 01 (P a) 1.0e 9

1.e 6 (0.2 + 1.5e For this test, the shear stress is fully determined by the material behavior model independently of the equilibrium equation. The equations are:

σ 12 = 2γ(t)(C 10 (Θ, ξ) + C 01 ) + 2G(Θ, ξ)B e12 (t) (55) Ḃe = L • Be + Be • L T - 1 τ (Θ, ξ) BD e • Be (56)
Using ( 53) and ( 54) in ( 56), one can obtain a system of three differential equation with three unknowns (B e11 , B e12 , B e22 ). This system is solved using the NDSolve function of Mathematica (Wolfram Research ( 2014)), Θ and ξ are fixed. The results are synthetized on the figures 4. These results clearly show that the chemical conversion has a stiffening effect on the dynamical mechanical response (at a fixed temperature) and that temperature increase leads to a softer and less dissipative dynamical response (at fixed chemical conversion).

Example 2: Isothermal case with fixed hydrostatic pressure

It is considered a block of matter in isothermal conditions and submitted to an hydrostatic pressure. The thermal evolution and dilatation effect are neglected. The deformation, the temperature and the chemical state are assumed to be homogeneous. In this case, the equations of the model become:

J ξ = 1 + βg(ξ) (57) p = K v (J m -1)J -1 ξ = K v ( J J ξ -1)J -1 ξ (58) ξ = A J exp -Ea RΘ C 2 Θ i log Θ Θ i (1 -ξ) m + β ∂g ∂ξ J -1 ξ (ξ)p , ξ(t = 0) = 0 (59)
From eqs ( 58) and ( 57) it can be obtained the volume variation as a function of the hydrostatic pressure and the chemical conversion:

J = p K v (1 + βg(ξ)) 2 + (1 + βg(ξ)) (60)
Considering the case p << K v this simplify to: J = J ξ = (1 + βg(ξ)), inserting this result in eq (59) the chemical conversion is therefore defined by: By noting f (ξ) = -g (ξ)(1 + βg(ξ)) -1 the influence of chemical shrinkage function on kinetic evolution, eq (61) can be written as:

ξ = A exp -Ea RΘ 1 + βg(ξ) C 2 Θ i log Θ Θ i (1 -ξ) m + β ∂g ∂ξ (1 + βg(ξ)) -1 p , ξ(t = 0) = 0 (
ξ = A exp -Ea RΘ 1 + βg(ξ) C 2 Θ i log Θ Θ i (1 -ξ) m + βf (ξ)p , ξ(t = 0) = 0 (62)
The previous differential equation is numerically solved with the NDSolve function of Mathematica.

Figures 5 show the response of the model in isothermal condition at fixed temperature. As observed for rubber system, the chemical conversion is very sensitive to the temperature. The influence of the hydrostatic pressure can be clearly observed: hydrostatic compression leads to an acceleration of the conversion contrary to hydrostatic tension that could stop the chemical process when a limit is reached. The relative influence of mechanic on chemistry can be very sensitive to the shrinkage coefficient as can be seen on figure 5(c).

Example 3: adiabatic case with hydrostatic pressure, application to material processing

In this example it is considered the fictitious case of the vulcanization of a block of rubber. The matter is submitted to a thermal and pressure cycle which is supposed to be representative of rubber curing process. No heat exchange is considered with the outside (adiabatic case) and the temperature is assumed to be homogeneous. Thermal diffusion is neglected however thermal dilatation is taken into account. The deformation is assumed to be homogenous and purely hydrostatic. The figure 6 shows the pressure and heating cycle considered: pressure is first applied linearly at the end of this ramp heating start and is then released. Pressure is kept fixed during chemical conversion followed by a cooling phase. Finally pressure is released to zero. In the following the amplitude values of heat source are the same for each test (max/min value of 6 M J/m 3 of initial volume) and different values of hydrostatic pressure amplitude and sign are chosen to illustrate the influence of mechanical couplings. For this problem as the deformation is purely hydrostatic, the system of equations is reduced to:

p = K v J J Θ (Θ)J ξ (ξ) -1 1 J Θ (Θ)J ξ (ξ) (63) Θ = J φ ξ (Θ, ξ, J) + l m (Θ, ξ, J) + l c (Θ, ξ, J) + J -1 ρ 0 r ρ 0 C(Θ, ξ, J) , Θ(t = 0) = Θ 0 (64) ξ = A J exp -Ea RΘ C 2 Θ i log Θ Θ i (1 -ξ) m + β ∂g ∂ξ J -1 ξ (ξ)p , ξ(t = 0) = 0 (65)
The previous system of equations is resolved numerically with a standard Euler (explicit) method with a small step size (dt = 0.2s). Each physics are therefore integrated at the same times with the same scheme (Mathematica is also used for this test). The results of the numerical simulation are presented on figures 7. It can be made the following observations: -a hydrostatic tension state of stress (p = 50M pa) leads to the smallest chemical conversion and therefore reach the smaller temperature during curing; -the higher the compression the fastest is the chemical evolution and as the reaction is exothermic the higher is the temperature. The competition between thermal dilatation/chemical shrinkage can be clearly seen on the figure that shows the volume variation. The figures 8 show that at the end of curing after cooling with the same amount of volumetric energy than this of heating and after pressure release the final temperature is slightly different than the initial one which is of 20C. The final temperature is also dependent on the hydrostatic pressure applied during curing. The final shrinkage value correspond to the value given for the coefficient β, i.e. 3% of volume variation at the maximum conversion. The final value of volume variation is therefore directly related to the reached conversion value at the end. The figure 9, illustrate the dependency of heat capacity upon temperature and chemical conversion. 4.4. Example 4: adiabatic case with cyclic shearing and heating, application to chemo-thermal aging upon thermo-mechanical loadings In this example, it is considered the case of an homogeneous sinusoidal shear test with an initial volumetric heating (as previously with an imposed value of ρ 0 r). No heat exchange is considered therefore the solution 14 is assumed to be homogeneous. The total deformation gradient is defined from eq (53) and the viscoelastic Cauchy Green tensor is defined from eq (54). The thermal and the chemical shrinkage are taken into account, the imposed deformation gradient leads to J = 1 and therefore J m is explicitly defined as a function of ξ and Θ:

J m = J -1 ξ J -1 Θ = (1 + βg(ξ)) -1 (1 + α(Θ -Θ 0 )) -1 (66)
In this case hydrostatic pressure is defined by:

p = K v (J -1 ξ J -1 Θ -1)J -1 ξ J -1 Θ (67)
As for the previous example, it is imposed in the first 30 mechanical cycles a volumetric heating to initiate chemical reaction (initial condition of Θ 0 = 20C is assumed). The volumetric heating is defined by: ρ 0 r(t) = 4.2e 6 × (t/30) J/m 3 t < 30 0 t ≥ 30 (68)

This idealized example is viewed as a very simplified simulation of thermo-chemical aging of a polymer material. These phenomena are distinct from the previous virtual curing processing however it is assumed that couplings upon each physics can be represented by the same model. As we do not intend to simulate a real material but only to qualitatively represent observed phenomena, the material parameters are kept identical than previous examples.

In the particular case of simple shear, the shear stress is explicitly given from the constitute equation, it is obtained the following set of equations:

σ 12 = 2γ(t)(C 10 (Θ, ξ) + C 01 ) + 2G(Θ, ξ)B e12 (t) (69) Ḃe = L • Be + Be • L T - 1 τ (Θ, ξ) BD e • Be (70) Θ = J ρ 0 C(Θ, ξ, Be ) φ m (Θ, ξ, Be ) + φ ξ (Θ, ξ, Be ) + l m (Θ, ξ, Be ) + l c (Θ, ξ, Be ) + J -1 ρ 0 r (71) ξ = A J exp -Ea RΘ C 2 Θ i log Θ Θ i (1 -ξ) m - ∂C 10 ∂ξ γ(t) 2 - ∂G ∂ξ (tr( Be ) -3) -β K v 2 ∂g ∂ξ J -1 ξ (J -1 ξ J -1 Θ -1) 2 (72) ξ(0) = 0, Θ(0) = Θ 0 , Be (0) = 1 (73)
For this example, the mechanical evolution problem has the smaller time scale and it is adopted a staged coupling strategy to save computing time and memory. It is assumed that Θ and ξ does not evolve much during a mechanical period (T = 1s). Therefore Θ and ξ are considered as functions of the number of mechanical cycles N . It is defined the following quantity:

Θ(N ) = Θ N = 1 T cycle Θ(t)dt ξ(N ) = ξ N = 1 T cycle ξ(t)dt (74)
The system of equations ( 69) to ( 73) is replaced by:

         σ 12 (t) = 2γ(t)(C 10 (Θ N , ξ N ) + C 01 ) + 2G(Θ N , ξ N )B e12 (t) Ḃe (t) = L(t) • Be (t) + Be (t) • L(t) T - 1 τ (Θ N , ξ N ) Be (t) D • Be (t) Be (0) = 1 (75)                          δΘ N δN = cycle φ m (Θ N , ξ N , Be ) + φ ξ (Θ N , ξ N , Be ) + l m (Θ N , ξ N , Be ) + l c (Θ N , ξ N , Be ) + ρ 0 J -1 r dt cycle J -1 ρ 0 C(Θ N , ξ N , Be (t))dt δξ N δN = A exp -Ea RΘ cycle J -1 C 2 Θ i log Θ N Θ i (1 -ξ N ) m - ∂C 10 ∂ξ γ(t) 2 - ∂G ∂ξ (tr( Be ) -3) -β K v 2 ∂g ∂ξ J -1 ξ (J -1 ξ J -1 Θ -1) 2 dt ξ N =0 = 0, Θ N =0 = Θ 0 , (76)
For a given cycle N , the following resolution is adopted: eq (75) are integrated with a forward Euler method. Temperature and chemical state are frozen and all the integral quantity needed to compute system (76) are numerically computed with a rectangle rule. System ( 76) is also integrated with a rectangle rule.

Figures 10 show the results of the numerical simulation up to N = 3000 cycles. As expected, the higher is the mechanical amplitude the higher is the temperature and the faster is the chemical conversion. The chemo-mechanical coupling is controlled by the dependency of mechanical parameter upon chemical conversion and by the ratio of A × C 2 versus β × A. For the chosen values of C 2 and A, it can be seen that a small change of β can lead to strong differences for the chemical conversion. In this example and for the chosen material parameters, it is needed to bring external heating for initiating chemical reactions. The proposed model is not limited to this case and it could be considered the case of chemical reaction activated only by self heating due to mechanical loadings (however more cycles would be necessary).

Remarks and conclusion

In this paper, we have developed a rigorous thermodynamical framework for the thermo-chemo-mechanical couplings in soft materials at finite strain. A thermo-viscoelastic model with a kinetic approach for chemistry is proposed. The main originality resides in the definition of the chemical Helmoltz free energy potential and the associated chemical flow rule. The free energy potential takes into account of an initiation temperature under which energy is stored as heat and above which energy is released to initiate chemical conversion. The chemical conversion which is derived from the energy potential naturally takes into account of the hydrostatic pressure. This two points lead to a clear definition of the heat capacity (and the respective influence of temperature and chemistry on it) and to a clear influence of hydrostatic pressure on chemical process. The proposed model is a first step to analyze and simulate complex multiphysics phenomena that take place during processing or aging of soft materials submitted to large strain. In this paper it has been chosen to consider only idealized application without any characterization of a real material. Therefore the chemical (one phenomenological reaction scheme) and the mechanical (zener model) behaviors are chosen as simple as possible to focus on the influence of chemo-mechanical couplings. For a specific application, the question of material parameters identification will certainly become a central one. In the case of non-linear mechanical behavior at finite strain the identification of phenomenological model usually required several experiments to be valid. As example simple DMA experiments are in general not sufficient. Conducing these experiments at controlled (homogeneous in the sample and stable) chemical states seems challenging. Nevertheless, if one consider curing of polymers (or rubbers), chemorheology can give useful informations in the linear regime (see for instance [START_REF] Joshi | Modeling of steady and time-dependent responses in filled, uncured, and crosslinked rubbers[END_REF]). It is therefore needed to develop new testing methods and experimental procedures that could help to investigate the non-linear regime or to impose complex mechanical loadings within a rheometer (for instance a precise control of the hydrostatic pressure, see [START_REF] Mackley | Experimental observations on the pressure-dependent polymer melt rheology of linear low density polyethylene, using a multi-pass rheometer[END_REF]). We believe that numerical simulation with thermo-chemo-mechanical models will greatly help to introduce new experimental technics and will be useful for a better understanding of existing ones. The finite element implementation of the present model requires to pay a careful attention on the variational formulation and on the numerical integrations schemes (for both evolution and conservation equations). The intrinsic times of each physics can be greatly different depending on the application (curing, aging, ...), on the material (polymer, rubbers, ...) and the loading conditions. Depending on these characteristics a monolithic or a staged strategy can be more appropriate. A first attempt has been proposed in Nguyen [START_REF] Nguyen | A finite strain thermo-chemo-mechanical coupled model for filled rubber[END_REF]; [START_REF] Eyheramendy | Advances in symbolic and numerical approaches in computational mechanics[END_REF] and need further improvements. These numerical developments will permit to take into account of more complex and localized phenomena like heat conduction (thermal gradients) or chemical diffusion.

Figure 1 :

 1 Figure 1: Dependence of the density of HPAB polymer with different degrees of cure on temperature. The legend shows the concentration of cross-links. Issued from[START_REF] Likozar | Cross-linking of polymers: Kinetics and transport phenomena[END_REF]

Figure 2 :

 2 Figure 2: A visco-elastic model

Figure 3 :

 3 Figure 3: Shrinkage function g(ξ)

  Cyclic shear test at ξ = 0.5 for various temperature values. Thermal softening and hysteresis area decreasing due to temperature increase.

Figure 4 :

 4 Figure 4: Mechanical response to cyclic shear test at constant temperature and constant chemical conversion values.

  Isothermal chemical conversion at p = 0M pa, effect of the temperature. Isothermal chemical conversion at Θ = 160C and p = -100M pa, effect of the shrinkage coefficient.

Figure 5 :

 5 Figure 5: Chemical conversion at fixed temperature and fixed hydrostatic pressure.

Figure 6 :

 6 Figure 6: Hydrostatic pressure and volumetric heating cycle.

Figure 7 :

 7 Figure 7: Results of a virtual material processing test under various hydrostatic pressure.

Figure 9 :

 9 Figure 8: Results of a virtual material processing test under various hydrostatic pressure: zoom at the end of curing.

  Figure 10: Thermo-Chemical aging for cyclic shearing

Table 1 :

 1 Material parameters: mechanical and thermal parts

	-3 (Θ -273) + 0.35ξ)	1.e 4

Chemical reactions in soft material often involve matter (gaz) diffusion and therefore the problem should be considered as on open thermodynamical system, however we neglected this effect in this work.

derivative of quantity upon time holding initial position constant: ρ = (∂ρ(x, t)/∂t) X

The proposed thermodynamical framework can be easily extended to more complex reacting systems: one has to define at least the same number of internal variables that the number of reactions.

< f >= f if f > 0 and < f >= 0 if f ≤ 0

A negative chemical evolution (reversion) could be thermodynamically admissible, however the authors believe that one have to introduce supplementary mechanisms (reactions or damage) that leads to reversion to be consistent.