B Chatzigeorgiou

F Kiefer

Meraghni

A Benaarbia
email: adil.benaarbia@ensam.eu

G Chatzigeorgiou
email: georges.chatzigeorgiou@ensam.eu

B Kiefer

F Meraghni
email: fodil.meraghni@ensam.eu

A fully

Keywords: Thermoplastic polymers, energy rate balance, Taylor-Quinney, cyclic viscoelasticity-viscoplasticity, ductile damage

A fully coupled thermo-viscoelastic-viscoplastic-damage framework to study the cyclic variability of the Taylor-Quinney coefficient for semi-crystalline polymers

able viscoelastic effects occurring on different time scales, irreversible strains observed above the stress threshold and energy responses arising from the anelastic deformation process (storage of energy, dissipation of energy, thermomechanical coupling, etc.). One of the aims of this paper is to evaluate the partitioning between stored and heat energy without assuming any

Taylor-Quinney values. The proposed approach is incorporated into the framework of Thermodynamics of Irreversible Processes and Generalized Standard Materials formalism to offer the thermodynamic consistency of all the constitutive equations. The numerical algorithm for the proposed model is implemented in the well-known finite element code Abaqus via the user material subroutine UMAT using an implicit formulation of the constitutive equations coupled with a radial return mapping algorithm. The model is then calibrated and validated across monotonic tensile and cyclic tensile-tensile tests by comparing predicted and experimentally observed energy responses. This comparison shows a good level of accordance between experimental findings and model predictions in terms of stress-strain responses for both monotonic and cyclic thermomechanical loading conditions. The model can also numerically capture the

nomenclature

Introduction

Semi-crystalline polymers (SCPs) are widely used in a large range of industrial sectors (e.g. civil engineering, aeronautics, automotive, etc.). The great interest in SCPs is mainly motivated by their ease of manufacture, low density and long lifespan. The modeling process of such SCPs has been widely investigated in the literature, and several constitutive models have been developed to capture the complex mechanical behavior of SCPs at both short-(impact) and longterm (fatigue, creep, etc.) deformation conditions [START_REF] Dusunceli | Modelling effects of degree of crystallinity on mechanical behavior of semicrystalline polymers[END_REF][START_REF] Launay | Cyclic behaviour of short glass fibre reinforced polyamide: Experimental study and constitutive equations[END_REF][START_REF] Uchida | Micro-, meso-to macroscopic modeling of deformation behavior of semi-crystalline polymer[END_REF][START_REF] Krairi | A thermodynamically-based constitutive model for thermoplastic polymers coupling viscoelasticity, viscoplasticity and ductile damage[END_REF][START_REF] Maurel-Pantel | A thermo-mechanical large deformation constitutive model for polymers based on material network description: Application to a semi-crystalline polyamide 66[END_REF][START_REF] Gudimetla | A finite strain thermodynamically-based constitutive framework coupling viscoelasticity and viscoplasticity with application to glassy polymers[END_REF][START_REF] Praud | Phenomenological multimechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation[END_REF]. However, most of the models generated are still inconclusive and limited to mechanical realm, which makes their constitutive formulations thermomechanically inconsistent. The numerical implementation of such models in commercial simulation software is also one of the major limiting factors that constrains their rational engineering use in real world structures. Models correctly describing the material behavior should properly estimate the thermomechanical damage and account for the energy behavior. The main benefit in developing fully-coupled thermomechanical models and their numerical implementation in a commercial simulation software is the ability to treat any kind of thermomechanical loadings consistently without making rough hypotheses, which offers the potential for using SCPs materials in structural applications in more intelligent ways.

Numerous experimental studies have been conducted in the literature to substantially understand the leading failure mechanisms common in SCPs [START_REF] Benaarbia | Influence of relative humidity and loading frequency on the pa6.6 thermomechanical cyclic behavior: Part ii. energy aspects[END_REF][START_REF] Federico | Large strain/time dependent mechanical behaviour of pmmas of different chain architectures. application of time-temperature superposition principle[END_REF]. Through examination of past research, strong interaction between time-dependent (viscous effects), time-independent (plasticity effects), and temperature-activated (thermoelasticity, entropic elasticity, etc.) deformation mechanisms has been emphasized [START_REF] Benaarbia | Thermomechanical analysis of the onset of strain concentration zones in wet polyamide 6.6 subjected to cyclic loading[END_REF], which makes the modelling of the material behavior and the prediction of lifetime of SCPs more complex. The modelling of viscoplasticity can be performed based on unified and non-unified viscoplasticity theories [START_REF] Chaboche | On the plastic and viscoplastic constitutive equations -part i: Rules developed with internal variable concept[END_REF][START_REF] Contesti | Description of creep-plasticity interaction with non-unified constitutive equations: Application to an austenitic stainless steel. 6th International seminar on inelastic analysis and life prediction in high temperature environment[END_REF]. In the unified viscoplasticity framework, the inelastic strain is confined to one single measure and represents both rate-dependent viscosity and rate-independent plasticity [START_REF] Chaboche | A review of some plasticity and viscoplasticity constitutive theories[END_REF][START_REF] Al-Rub | A thermodynamic framework for constitutive modeling of time-and rate-dependent materials. part i: Theory[END_REF]Benaarbia et al., 2018a,b). In non-unified viscoplasticity theories, both mechanisms are calculated separately and controlled by different flow rules [START_REF] Cailletaud | Study of plastic/viscoplastic models with various inelastic mechanisms[END_REF], which makes their numerical implementation difficult. Non-unified models lack the ability of predicting some of the material behaviors such as ratcheting and interaction between viscous and plastic effects [START_REF] Krempl | Viscoplastic models for high temperature applications[END_REF]. Numerous viscoplastic models for polymers have been developed based on both unified and non-unified viscoplasticity theories [START_REF] Anand | A thermo-mechanically coupled theory for large deformations of amorphous polymers. part i: Formulation[END_REF][START_REF] Ames | A thermo-mechanically coupled theory for large deformations of amorphous polymers. part ii: Applications[END_REF][START_REF] Zairi | A physically-based constitutive model for anisotropic damage in rubber-toughened glassy polymers during finite deformation[END_REF][START_REF] Poulain | Finite-strain elasto-viscoplastic behavior of an epoxy resin: Experiments and modeling in the glassy regime[END_REF][START_REF] Krairi | A thermodynamically-based constitutive model for thermoplastic polymers coupling viscoelasticity, viscoplasticity and ductile damage[END_REF][START_REF] Gudimetla | A finite strain thermodynamically-based constitutive framework coupling viscoelasticity and viscoplasticity with application to glassy polymers[END_REF][START_REF] Achour | Implicit implementation and consistent tangent modulus of a viscoplastic model for polymers[END_REF][START_REF] Krairi | Thermo-mechanical coupling of a viscoelastic-viscoplastic model for thermoplastic polymers: thermodynamical derivation and experimental assessment[END_REF], which include theories of elasto-viscoplasticity, coupled viscoelasticity and viscoplasticity, coupled viscoplasticity and damage, etc. The pioneering viscoplastic model developed in the literature is the Chaboche and Rousselier model [START_REF] Chaboche | On the plastic and viscoplastic constitutive equations -part i: Rules developed with internal variable concept[END_REF], initially designed for high temperature metals. Since then, a variety of authors proposed their versions based on different forms of non-linear hardening rules and viscosity functions [START_REF] Kocks | Thermodynamics and kinetics of slip[END_REF][START_REF] Miller | An inelastic constitutive model for monotonic, cyclic, and creep deformation: Part i -equations development and analytical procedures[END_REF][START_REF] Delobelle | Sur les lois de comportement viscoplastique à variables internes[END_REF][START_REF] Krempl | Rate (time) -dependent deformation behavior: an overview of some properties of metals and solid polymers[END_REF][START_REF] Kang | Constitutive modeling of strain range dependent cyclic hardening[END_REF][START_REF] Chaboche | A review of some plasticity and viscoplasticity constitutive theories[END_REF][START_REF] Abdel-Karim | Cyclic multiaxial and shear finite deformation responses of ofhc cu. part ii: An extension to the khl model and simulations[END_REF][START_REF] Zhang | Interaction of cyclic softening and stress relaxation of 9-12% cr steel under strain-controlled fatigue-creep condition: Experimental and modeling[END_REF].

From a mechanics of materials standpoint, there have been many attempts to develop constitutive models in the context of the thermodynamics of irreversbile processes [START_REF] Callen | Thermodynamics: an introduction to the physical theories of equilibrium thermostatics and irreversible thermodynamics[END_REF][START_REF] Germain | Thermodynamics continuum mechanics[END_REF]. The idea behind this is to define internal state variables that regulate the storage and heat energies through material structure and deformation mechanisms [START_REF] Ranc | Calorimetric consequences of thermal softening in johnsoncook's model[END_REF]. These internal variables are generally introduced in order to describe the thermodynamic/microstructural state of materials under consideration. Historically, The thermodynamics of Irreversible Processes (TIP) frame was initially developed as a linear theory and then coupled with the Generalised Standard Material (GSM) formalism to expand from the limitations of linearity [START_REF] Nguyen | Matériaux élastoviscoplastiques a potentiel généralisé[END_REF][START_REF] Halphen | On the generalized standards materials (in french)[END_REF]. The GSM formalism is widely employed in the literature [START_REF] Lemaitre | Coupled elasto-plasticity and damage constitutive equations[END_REF][START_REF] Nguyen | Plastic bifurcation and postbifurcation analysis for generalized standard continua[END_REF][START_REF] Lemaitre | Mechanics of solid materials[END_REF][START_REF] Lubliner | Plasticity theory[END_REF], and has been used to ensure the satisfaction of the dissipation inequality by deriving constitutive equations from both the thermodynamic and the dissipation potential. Constitutive models are comprised of evolution laws and equations of state. The form of the state laws and the dissipation inequality are specified by the choice of the internal state variables and the particular form of the thermodynamic potential, whilst evolution formulations are derived from gradients of dissipation potentials. The difficulty in developing reliable constitutive models directly relies on the choice of both thermodynamic and dissipation potentials as well as internal state variables. SCPs undergo various processes in the fatigue regime; these include prior thermoelasticity, viscoelasticity at different time scales (highlighting the time dependent response of the material to the mechanical loading) coupled phenomenologically with ductile damage (represented by a material stiffness reduction variable)

and non-linear viscoplasticity causing the irreversible strain acccumulation generally observed above the stress threshold (Benaarbia et al., 2014b[START_REF] Benaarbia | Influence of relative humidity and loading frequency on the pa6.6 thermomechanical cyclic behavior: Part ii. energy aspects[END_REF][START_REF] Benaarbia | Thermomechanical analysis of the onset of strain concentration zones in wet polyamide 6.6 subjected to cyclic loading[END_REF].

Under cyclic loading conditions, the self-heating is mainly induced by two different types of heat sources: dissipation and coupling sources. The former source is attributed to the material irreversibility while the latter is associated with either standard thermoelasticity or rubber elasticity, also termed entropic elasticity. The standard thermoelasticity reflects the thermodilatability of the material while the entropic elasticity is related to the stress-induced orientation/disorientation of macromolecular chains. The possible competition between these two coupling sources is associated with the famous thermoelastic inversion effect first observed by [START_REF] Joule | On some thermodynamic properties of solids[END_REF]. Several attempts have been reported to model this thermoelastic inversion [START_REF] Anthony | Equations of state for naturals and synthetic rubber like materials: unaccelerated natural soft rubber[END_REF][START_REF] Chadwick | Modified entropic elasticity of rubberlike materials[END_REF][START_REF] Ogden | Aspects of the phenomenological theory of rubber thermoelasticity[END_REF]. The modelling framework is called modified entropic elasticity and belongs to the nonlinear hyperelasticity domain. In previous works (Benaarbia et al., 2014a[START_REF] Benaarbia | Thermomechanical analysis of the onset of strain concentration zones in wet polyamide 6.6 subjected to cyclic loading[END_REF], the competition of both coupling sources have been modelled using the Generalized Standard Material formalism.

In this paper a continuum thermodynamics model for SCPs is proposed to capture the governing deformation mechanisms. Constitutive equations of the model are derived from the TIP framework coupled with the GSM formalism. The model combines linear viscoelasticities coupled with a power flow rule (in the viscoplastic contribution) and non-linear isothermal cyclic evolution of isotropic hardening to capture the irreversible strain evolution.

The highly coupled nature of these mechanisms makes material parameter identification very challenging to perform. Unfortunately, constitutive equations in such situations depend on a great number of parameters that must be identified using a limited number of experiments. Thus, a step by step optimisation procedure is ultimately required to fine-tune all material parameters to properly represent these highly coupled phenomena. Another major difficulty inherent to these highly coupled models is the development of a reliable numerical implementation strategy that enables their multiaxial use for real-life engineering components. This strategy represents often a laborious task and requires constitutive tangent tensors to be implemented in a user-defined material routine. In the authors' previous work [START_REF] Praud | Phenomenological multimechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation[END_REF], a time implicit implementation of an isothermal viscoelastic-viscoplastic-damage model has been developed based on the return mapping algorithm using the convex cutting plane form [START_REF] Ortiz | An analysis of a new class of integration algorithms for elastoplastic constitutive relations[END_REF][START_REF] Simo | Computational Inelasticity[END_REF][START_REF] Qidwai | Numerical implementation of a shape memory alloy thermomechanical constitutive model using return mapping algorithms[END_REF]. This research work is intending to follow the same numerical implementation procedure by extending the existing viscoelastic-viscoplastic-damage model to cover dissipative effects (self-heating) and account for themomechanical couplings. To the author's knowledge, the numerical implementation of fully-coupled thermo-viscoelastic-viscoplastic-damage models under both the TIP and GSM frameworks has never been addressed, and this research work by developing the details of the associated computational tools is therefore an original contribution. In addition to this aspect, the principal originality and contribution of this work, compared to existing thermomechanical studies [START_REF] Yu | A thermo-mechanically coupled nonlinear viscoelasticviscoplastic cyclic constitutive model for polymeric materials[END_REF][START_REF] Krairi | Thermo-mechanical coupling of a viscoelastic-viscoplastic model for thermoplastic polymers: thermodynamical derivation and experimental assessment[END_REF], is a) the integration of the thermomechanical damage, which is a very important factor in the behavior of the SCPs and b) the investigation, both experimentally and numerically, of the dissipation and the stored energies. These two energetic terms are very important for fatigue testing of polymeric materials, since they can assist in developing proper thermomechanically-driven lifetime criteria. In addition, this paper proposes for the first time a numerical modelling framework for evaluating the cyclic thermomechanical evolution of the Taylor-Quinney coefficient for polymeric materials.

The present paper is organized as follows: in section 3 a brief introduction to the theoretical background used to derive the heat diffusion equation is presented followed by a short experimental description of the energy balance; in section 4 the constitutive equations of the thermo-viscoelastic-viscoplastic-damage model are derived using both the TPI and GSM formalisms; in section 5 the numerical implementation of constitutive laws based on the return mapping and the convex cutting plane schemes is presented; in section 6, the first part is devoted to the model validation by comparing numerical computations with experimental results.The second part presents structural FE simulations conducting a parametric study on the kinetics of the Taylor-Quinney coefficient and damage; finally in section 7 concluding remarks are presented and some prospects are suggested for further work.

Preliminaries

Heat diffusion equation

Polymer materials can undergo large deformations at high moisture conditions. In the sequel, the framework of small deformations and rotations is adopted for two reasons: i) the main goal of developing this new constitutive law is to integrate it in multiscale approaches for polymer based composites, which usually do not undergo large deformations, ii) coupling of three dissipative mechanisms (viscoelasticity, viscoplasticity, damage) in a large deformations framework is a hard task, especially due to the multiplicative decomposition requirement of the total deformation gradient.

Under the small deformation framework, the first principle of thermodynamics can be given as ρ ė = σ : ˙ σ : ˙ σ : ˙ + r e -div q q q, (1)

where ρ, ė and r e stand for the material density, the rate of specific internal energy and the external heat supply, respectively, while σ σ σ, ˙ ˙ ˙ and q q q are the Cauchy stress tensor, the total strain rate tensor, the heat influx vector, respectively. In Eq. (1), the symbol : : : denotes the double contraction product.

Irreversible processes are defined from the second principle of thermodynamics (Clausius-Duhem inequality). The expression for the second law of thermodynamics can be written as ρ ṡ ≥ ϑ -1 r e -div q q q + 1 ϑ q q q.∇ ∇ ∇ ϑ ,

where ϑ stands for the absolute temperature, defined as ϑ = ∂ s e, and s represents the material specific entropy. By combining the expressions of the first and second principles, the heat diffusion equation can be derived.

The thermodynamic state of a material can be described by a finite set of state variables [START_REF] Germain | Thermodynamics continuum mechanics[END_REF]. These can include the absolute temperature, ϑ, the overall strain tensor, and additional internal state variables, denoted here by the set ג{ ג ג k } N k=1 . For the simplicity of the expressions presented henceforth, the following operator is introduced

∂ ג ג ג k {•} ג ג ג k = ∑ p ∂ x p {•} ẋp + ∑ q ∂ x x x q {•} . . . ẋ x x q + ∑ n ∂ X X X n {•} : : : Ẋ X X n , (3)
where x p can be scalars, x x x q vectors and X X X n second-order tensors. ∂ y y y x x x denotes the partial derivative of x x x with respect to y y y.

In the GSM framework, two potentials are used to describe the material behavior: the thermodynamic potential, ψ, and the dissipation potential, ϕ. The thermodynamic potential is classically considered as the Helmholtz free energy. On the other hand, irreversibilities are introduced through the dissipation potential, a convex function of the state variable rates ( ˙ , { ג ג ג k } N k=1 and q q q). Using the definition of the dissipation potential and the Clausius-Duhem inequality, the total dissipation Ω can be written as

Ω = ∂ ∂ ∂ ˙ ϕ : : : ˙ + ∂ ∂ ∂ ג ג ג k ϕ ג ג ג k + ∂ ∂ ∂ q q q ϕ . . . q q q, = [σ σ σ -ρ ∂ ∂ ∂ ψ] : : : ˙ ˙ ˙ -ρ ∂ ∂ ∂ ג ג ג k ψ ג ג ג k -ϑ -1 q q q. . .∇ ∇ ∇ϑ ≥ 0. ( 4
)
The total dissipation is generally split into two parts which are assumed to be separately positive: the intrinsic dissipation

Ω int = [σ σ σ -ρ ∂ ∂ ∂ ψ] : : : ˙ ˙ ˙ -ρ ∂ ∂ ∂ ג ג ג k ψ ג ג ג k ,

and thermal dissipation

Ω the = -ϑ -1 q q q. . .∇ ∇ ∇ϑ. The intrinsic dissipation is associated with the material degradation while the thermal dissipation is induced by heat diffusion.

By defining the specific heat as c = -ϑ ∂ , ϑϑ ψ at constant and ג ג ג k , the local expression of the heat diffusion equation can be written as

ρc θ -Q = [σ σ σ -ρ ∂ ∂ ∂ ψ] : : : ˙ ˙ ˙ -ρ ∂ ∂ ∂ ג ג ג k ψ ג ג ג k Ω int + ρϑ ∂ ∂ ∂ , ϑ ψ : : : ˙ ˙ ˙ + ρϑ ∂ ∂ ∂ , ϑ ג ג ג k ψ ג ג ג k C tm , ( 5
)
where the thermal inertia ρc θ represents the stored or released heat rate, while Q = r e -div q q q stands for the volume heat losses by conduction and radiation. The intrinsic mechanical dissip-ation Ω int and the thermomechanical coupling sources C tm are pooled on the right hand side.

The thermomechanical coupling sources are the caloric signatures of the possible interactions between the temperature and the other mechanical or microstructural state variables.

Taylor-Quinney coefficient

When a material is inelastically deformed, a part of the mechanical energy, W de f , expended in the deformation process is converted into heat, W d , but the remainder is stored in the material, W s , thereby raising its internal energy. Values of the dissipation energy rate can be obtained as the difference between the deformation energy rate, W de f , spent during the material transformation and the sum of the elastic, W e , and stored, W s , energy rates. Following this approach, it follows that

Ω int = σ σ σ : : : ˙ W de f -ρ ∂ ∂ ∂ ψ : : : ˙ + ρ ∂ ∂ ∂ ג ג ג k ψ ג ג ג k W e +W s . ( 6
)
The Taylor-Quinney coefficient can be introduced to quantify the stored energy fraction (in the integral β int and differential β di f f forms)

β int := W s W in ; β di f f := W s W in = 1 - W d W de f -W e , ( 7
)
where W in represents the inelastic work and the non-standard â ȂIJ â ȂŸ â Ȃ˙I notation specifies that the energy rate is path-dependent (i.e. the energy term is not necessarily a state function).

The Taylor-Quinney coefficients, β int and β di f f , represent measures of how efficiently inelastic work or power are converted into heat during the deformation process. β int expresses the fraction of stored to inelastic energy and is related to the microstructural transformations during the deformation process. It reflects how much of inelastic strain energy is stored in the microstructure and how much is converted into heat. β int relates to stored and inelastic powers instead of energies and has a similar physical meaning.

In the literature, this coefficient is commonly supposed to be a phenomenological constant. It is used to allow estimating the dissipation-induced self-heating using a simple heat equation, without introducing a thermodynamic framework. The interested reader can find a wealth of references in the review of [START_REF] Bever | The stored energy of cold work[END_REF] on the Taylor-Quinney coefficient, also termed the stored energy ratio of cold work. The available Taylor-Quinney values in the literature are mostly limited to metallic materials and investigations on polymer materials are still scarce and confined to monotonous tests [START_REF] Shao | An internal state variable thermodynamic model for determining the taylor-quinney coefficient of glassy polymers[END_REF]. This research work aims at increasing our knowledge regarding the storage of energy arising from the inelastic cyclic deformation of semi-crystalline polyamides. This will enable us to ensure the thermomechanical consistency of constitutive models and to properly integrate the energy terms into classical simulation software.

Experimental observations

Many experimental results dealing with the cyclic thermomechanical behavior of polyamide 6.6 material have been thoroughly presented and summarized in the author's previous papers (Benaarbia et al., 2014a(Benaarbia et al., ,b, 2015[START_REF] Benaarbia | Thermomechanical analysis of the onset of strain concentration zones in wet polyamide 6.6 subjected to cyclic loading[END_REF]. Different relative humidity levels (dry and wet) and loading rates have been considered to investigate both time and moisture effects. Hysteretic, thermal and energy responses involved during the fatigue process of PA6.6 samples have been closely documented and analyzed. Through the examination of past research, many major features were observed under fatigue conditions:

• From a mechanical standpoint, it was found that the deformation energy (hysteresis area) decreased markedly at high relative humidity and low strain rates. The slope of hysteresis loops was also observed to decrease markedly at the very beginning of loading and seemed to stabilize after a few hundred cycles. A ratcheting behavior, characterized by a drift of the mean strain, was found to be significant for wet samples.

• From a thermal point of view, it was found that the cyclic tests cannot be considered as isothermal mechanical tests (the self-heating of dry samples was greatly lower than that of moist ones). It was also demonstrated that temperature variations were induced by two independent types of heat sources: dissipation and coupling sources. It was considered that the intrinsic dissipation was attributed to the material irreversibility. Two coupling effects were emphasized: standard thermoelasticity reflecting the thermodilatability of the material and rubber elasticity, also termed entropic elasticity, highlighting the stressinduced orientation/disorientation of macromolecular chains.

• From a general energy standpoint, it was found that both relative humidity and loading rate conditions affected the form and kinetics of the energy rate balance. Investigations on energy stored during deformation indicated that the stored ratio (i.e. the differential Taylor-Quinney coefficient) was significantly lower at low loading rates but remained high at high loading rates (see. Figure 1). These considerations are of great importance regarding the inclusion of energy balance kinetics into PA6.6 constitutive models.

Assuming a fixed value of the Taylor-Quinney coefficient has been and is still widely used for predicting mechanical behaviour in engineering contexts [START_REF] Zaera | On the taylor -quinney coefficient in dynamically phase transforming materials. application to 304 stainless steel[END_REF][START_REF] Maurel-Pantel | A thermo-mechanical large deformation constitutive model for polymers based on material network description: Application to a semi-crystalline polyamide 66[END_REF][START_REF] Rittel | The dependence of the taylor -quinney coefficient on the dynamic loading mode[END_REF]. One of the aims of this research work is to assess the partitioning between stored and heat energy as a predictor of microstructure evolution, without assuming any Taylor-Quinney values a priori.

Specification of constitutive equations for the investigated material

Experimental observations extracted from the author's previous papers have clearly shown that the response of the investigated material shows both recoverable and irrecoverable processes.

Recoverable mechanisms (prior thermoelasticity and viscoelasticities at different time-scales) are hereafter modeled using a solid-like thermo-viscoelasticity model, however, the irrecoverable mechanisms are captured using a time and rate-dependent plasticity model. The changes in the material's microstructure during deformation cause the investigated materials to experience a significant amount of damage under service conditions. The degradation in the mechanical properties of the material is caused by damage and cannot be explained only by viscoelasticity and/or viscoplasticity constitutive models. A Lemaitre model based on continuum damage mechanics can thus be used to account for degradations involved within the material [START_REF] Lemaitre | Coupled elasto-plasticity and damage constitutive equations[END_REF].

In the current context, viscoelastic, viscoplastic, damage and thermally-dependent material properties are used to describe the temperature, time-and rate-dependent behavior of the investigated material. In the following, { ve i } P i=1 , vp and r are, respectively, defined as viscoelastic strain tensors, viscoplastic strain tensor and isotropic hardening variable. The parameter P stands for the number of viscoelastic strain tensors, which are linked with Kelvin-Voigt type branches that are presented in section 6. For viscoelastic materials, one Kelvin-Voigt branch is not sufficient to describe the behaviour in a large time range and additional ones are required.

Each branch is linked with a characteristic time and covers a specific range of frequencies. The vector of internal state variables ג{ ג ג k } N k=1 thus reads ג ג ג = { ve i } P i=1 , vp , r, d , with d ∈ [0, 1] being responsible for damage. The choice of an isotropic hardening, imposing hence an isotropic ductile damage, is motivated by the fact that in wet polymer materials conditioned at a relative humidity 50% or higher, which is the examined material conditions in this work, the glass transition is around the room temperature. Thus, the material is beyond the glass transition domain and no tension-compression asymmetry response is observed (Gilles Robert, private communication). Nevertheless, in case of tension-compression asymmetry, a proper modification of the current model can be adapted by incorporating a kinematic hardening function that also affects the tension-compression damage kinetics.

Free energy, state laws and heat diffusion equation

Following what has been mentioned previously, an additively decoupled form of the thermodynamic potential is introduced

ψ ϑ, , { ve i } P i=1 , vp , r, d := ψ e ϑ, , { ve i } P i=1 , vp , d + ψ cal (ϑ) + P ∑ i=1 ψ ve i ( ve i , d) + ψ vp (r) . (8)
The right-hand side of the previous expression pools all the thermodynamic potential components:

• ψ e ϑ, , { ve i } P i=1 , vp , d denotes the damage-affected elastic free energy, which is assumed to be temperature independent in terms of elastic material parameters. This elastic component reads

ψ e ( , e , d) := 1 2 e : : : [1 -d] C C C e : : : e ; with e = - P ∑ i=1 ve i -vp -α α α [ϑ -ϑ o ] , (9)
where ϑ o stands for the reference temperature, while C C C e and α α α represent the initial fourthorder stiffness 1 and second-order thermal expansion tensors, respectively.

• ψ cal (ϑ) stands for the stored or released heat induced by temperature variations. Its expression for a reference heat capacity, c o , is given as

ψ cal (ϑ) = c o [ϑ -ϑ o ] -ϑ ln ϑ ϑ o -s o ϑ -e o , (10)
where s o and e o are both reference entropy and reference internal energy, respectively.

• ψ ve i ( ve i , d) represents the damage-affected viscoelastic free energies with i = 1, ..., P. Their 1 The fourth-order stiffness tensor has minor and major symmetries

C ijkl = C jikl = C ijlk ; C ijkl = C klij .
expressions are given as follows

ψ ve i ( ve i , d) := 1 2 ve i : : : [1 -d] C C C ve i : : : ve i with i = 1, ..., P, (11)
where C C C ve i represent the initial fourth order stiffness tensors associated with the viscoelasticities.

• Finally, ψ vp (r) denotes the viscoplastic free energy. This term is chosen such that

ψ vp (r) := F (r) = h 1 n vp + 1 r n vp +1 , ( 12
)
where n vp stands for the isotropic hardening exponent, while h 1 is a linear-type hardening coefficient. It is noted that under thermomechanical loading conditions, the viscoplastic free energy usually depends also on the temperature. In this work this dependence has not been investigated experimentally, and thus it is omitted. Nevertheless, the consequent computational algorithm is presented later in a general form, assuming that F ≡ F(r, θ).

Each individual state variable is associated with a conjugate variable defined as the partial derivative of the thermodynamic potential with respect to the state variable. Following the expressions of the foregoing elementary free energies, the state laws read:

• The total entropy (conjugate variable of temperature ϑ)

s := -∂ ϑ ψ = c o ln ϑ ϑ o + s o + α α α : : : σ σ σ; (13)
• The reversible stress tensor (conjugate variable of total strain tensor )

σ σ σ r := ∂ ∂ ∂ ψ = [1 -d] C C C e : : : e ; (14)
• The conjugate variables of the viscoelastic strain tensors ve i a a a ve

i := ∂ ∂ ∂ ve i ψ = [1 -d] C C C ve i : : : ve i -σ σ σ r with i = 1, ..., P; (15)
• The conjugate variable of the viscoplastic strain tensor vp a a a vp :

= ∂ ∂ ∂ vp ψ = -[1 -d] C C C e : : : e = -σ σ σ r ; (16)
• The drag stress (conjugate variable of the isotropic hardening variable r)

R := ∂ r ψ = ∂ r F; (17)
• The strain density release rate (conjugate variable associated with the damage variable, d)

a d := ∂ d ψ = - 1 2 e : : : σ σ σ r [1 -d] - 1 2 P ∑ i=1 ve i : : : C C C ve i : : : ve i . ( 18
)
Using Eq. ( 4), both thermal and intrinsic dissipations within the current framework can be rearranged as

Ω the = -ϑ -1 [1 -d] q q q.∇ ∇ ∇ϑ ≥ 0 , ( 19
) and damage (term c) phenomena. This means that no irreversibility is associated with the rate of the total strain tensor ˙ . Thus, the reversible stress σ σ σ r is chosen as identical to the total Cauchy stress tensor, σ σ σ.

Ω int = [σ σ σ -σ σ σ
By performing the partial derivatives of the free energy with regard to the internal variables, Eq. ( 5) can be rearranged as

ρ c θ -Q = σ σ σ : : : ˙ vp - P ∑ i=1
a a a ve i : : :

˙ ve i -R ṙ -a d ḋ Ω int + ρ ϑ ḋ [1 -d] α α α : : : σ σ σ -ρ ϑ [1 -d]
[C e : : : α α α C e : : : α α α C e : : : α α α] : : : ˙ -

P ∑ i=1 ˙ ve i -˙ vp C tm . ( 21
)
This equation is split into both dissipative and coupling heat sources. Note that only coupling between temperature and prior elasticity is considered in the current work since the hardening function is assumed to be independent of temperature.

Dissipation potential and evolution laws

As has already been mentioned, state laws associated with conservation equations are not sufficient to describe the irreversibilities accompanying the deformation of the material.

Since state laws only describe the properties of the equilibrium state, it is vital to introduce complementary equations which will, in turn, describe the evolution of the system and any possible irreversibilities that accompany it.

For the case of viscoelasticity, the GSM formalism postulates the existence of a so-called dual dissipation potential, ϕ , that reads (written as the sum of the viscous sub-potentials)

ϕ ve i {t t t ve i } P i=1 , d := 1 2 t t t ve i : [[1 -d] D D D ve i ] -1 : t t t ve i , t t t ve i = -a a a ve i , (22)
where t t t ve i stands for the thermodynamic forces associated with the viscoelastic strain fluxes { ˙ ve i } P i=1 , while D D D ve i represents the fourth-order viscosity tensors2 . Thus, the following evolution law for each individual viscoelastic strain tensor is given through the partial derivative of the dual potential with respect to the associated thermodynamic force

˙ ve i := ∂ ∂ ∂ t t t ve i ϕ = [[1 -d] D D D ve i ] -1 : : : t t t ve i . ( 23
)
Combining Eq. ( 15) with Eq. ( 23), one obtains relations governing each individual viscoelastic strain tensor. The respective differential form of these relations (the effective stress as a function of the viscoelastic strain tensors and their rates) can thus be written as

σ σ σ ( ve i , ˙ ve i ) = C C C ve i : : : ve i + D D D ve i : : : ˙ ve i i = 1, ..., P; with σ σ σ = [1 -d] -1 σ σ σ. ( 24
)
In the following, the evolution of the viscoplastic and damage internal state variables are not extracted from a dual dissipation potential, since both viscoplasticity and damage are considered as coupled mechanisms. Hence, using a J 2 -theory for viscoplasticity together with a damage potential, the rule of generalized normality expressed as

ג ג ג = ∂ ∂ ∂ t t t φ γ, ג ג ג = { ˙ vp , ṙ, ḋ}, t t t = {σ σ σ, -R, -a d }, (25)
is hereafter adopted to derive the evolution laws. The choice of the J 2 -theory implies that the current model does not capture hydrostatically dependent viscoplastic strains. In Eq. ( 25), γ stands for the viscoplastic multiplier while the indicative function φ is split into both viscoplasticity coupled with damage, φ vpd , and pure damage, φ d , functions as

φ (σ σ σ, -R, -a d ; r, d) := φ vpd ( σ σ σ, -R; r) + φ d (-a d ; d) . ( 26
)
The viscoplasticity coupled with damage is chosen as

φ vpd ( σ σ σ, -R; r) := f (σ σ σ, -R; d, r) , ( 27
)
where f (σ σ σ, -R; d, r) stands for the viscoelastic domain when f < 0 and viscoplastic flow when f > 0. The von Mises yield criterion (function of the invariants of the deviatoric stresses) is considered as follows (see [START_REF] Lemaitre | Mechanics of solid materials[END_REF])

f (σ σ σ, -R; d, r) = J 2 ( σ σ σ) -R (r) -Σ y , (28)
where Σ y defines the yield stress while J 2 stands for the second invariant defined as J 2 (σ σ σ) = 3/2 σ σ σ dev : : : σ σ σ dev with σ σ σ dev denotes the deviator of the total stress tensor 3 , σ σ σ.

The damage potential φ d for ductile evolutions is chosen as [START_REF] Lemaitre | Coupled elasto-plasticity and damage constitutive equations[END_REF])

φ d (-a d ; d) := K [1 -d] [n d + 1] - a d K n d +1 , ( 29
)
where n d stands for the damage exponent, while K represents a damage constant.

Following the generalized normality rules expressed in Eq. ( 25), one obtains the following evolution laws

• The viscoplastic flow

˙ vp := ∂ ∂ ∂ σ σ σ φ γ = [1 -d] -1 ñ n n γ, ( 30
)
where ñ n n stands for the effective viscoplastic flow direction (the partial derivative of the yield function with respect to the total effective stress tensor). Its expression is given as

ñ n n := ∂ ∂ ∂ σ σ σ f = 3 2 σ σ σ dev J 2 ( σ σ σ) . (31)
• The rate of the isotropic hardening internal variable

ṙ := -∂ R φ γ = γ. (32)
• The rate of the damage variable

ḋ := -∂ a d φ γ = 1 [1 -d] - a d K n d γ. ( 33
)
For the viscoplastic model, we adopt a power law viscosity function to define the viscoplastic

b) ijkl = a ij b kl , 1 1 1 ij = δ ij and I I I ijkl = 1 2 δ ik δ jl + δ il δ jk ,
where ⊗ stands for the dyadic product.

multiplier rate such that

γ := v vp (σ σ σ, -R; d, r, ϑ) ; v vp := f Z -1 vp p vp if f > 0, v vp = 0 if f ≤ 0, ( 34
)
where the parameter Z vp represents the viscoplastic resistance stress, while p vp stands for the viscoplastic exponent. 4

Numerical implementation

The numerical scheme utilized in this section follows the well-known return mapping algorithm. The general theory is presented in [START_REF] Simo | Computational Inelasticity[END_REF] and an extended version for materials under thermomechanically coupled processes is described in [START_REF] Chatzigeorgiou | Thermomechanical Behavior of Dissipative Composite Materials[END_REF]. This implicit numerical framework is suitable for quasi-static conditions. For known total strain field, the computation of the stress level and the internal state variables is performed through a prediction-correction type scheme. In addition, four thermomechanically related tangent moduli are estimated for the needs of the finite element implementation of the constitutive law. The essential details of the algorithmic scheme for the proposed model are briefly presented in the sequel. It should be mentioned that the indices n and n + 1 denote time steps, while the indices m and m + 1 denote iteration increments.

Constitutive law implementation

During the inelastic correction, both the viscoelastic and the viscoplastic criteria are assumed null such that

Φ Φ Φ (n+1)(m) ve i + δΦ δΦ δΦ (n+1)(m) ve i = 0 0 0; Φ (n+1)(m) r + δΦ (n+1)(m) r = 0, ( 35
)
where both increments δΦ δΦ δΦ ve i and δΦ r are determined by linearization of the state and evolution equations (see Table 2). For simplicity, all exponent notations are dropped whenever the quantities are taken at the increment (n + 1)(m) to reduce the complexity of indices and 4 It is worth mentioning that the viscoplasticity activation criterion (34) considers only the part f of the total potential ϕ. Thus, a type of non-associative flow rule is adopted in the present model. From a theoretical standpoint, one could provide a bipotential, generalizing the concept of GSM to the Implicit Standard Materials (de Saxcé and Bousshine, 2002). Evolution laws:

7. The viscoelastic strain rates

˙ ve i = -[[1 -d] D D D ve i ]
-1 : : : a a a ve i with i = 1, ..., P .

The viscoplastic multiplier ˙

vp = Λ Λ Λ r ṙ , Λ Λ Λ r = 3 2 [1 -d] -1 σ σ σ dev J 2 ( σ σ σ) , σ σ σ = [1 -d] -1 σ σ σ.

The damage variable rate

ḋ = Λ d ṙ , Λ d = 1 [1 -d] - a d K n d .
10. The intrinsic dissipation

Ω int = - P ∑ i=1
a a a ve i : : : ˙ ve i -a a a vp : : :

˙ vp -R ṙ -A d ḋ ≥ 0.
11. The heat diffusion equation h s = -ϑ ṡ + Ω int where h s stands for the heat sources.

Criteria of leading mechanisms:

12. The viscoelastic criteria

Φ Φ Φ ve i = [[1 -d] D D D ve i ] -1 : : : a a a ve i + ˙ ve i = 0 0 0 with i = 1, ..., P . 13. The viscoplastic criterion Φ r = J 2 ( σ σ σ) -R -Σ y Z vp p vp + -ṙ = 0.
subscripts. Following the stress-strain state law given in Eq. ( 14) and the internal variable rate forms expressed in Eqs. ( 30)-( 33), the linearization of the stress tensor, σ σ σ, can be written as

δσ δσ δσ = G G G σ : : : δ δ δ + g g g σϑ δϑ -G G G σ : : : ∑ i δ δ δ ve i + g g g σr δr , ( 36
)
G G G σ = [1 -d] C C C e , g g g σϑ = -[1 -d] C C C : : : α α α, g g g σr = -[1 -d] -1 Λ d σ σ σ -[1 -d] C C C e : : : Λ Λ Λ r .
Combining Eq. ( 36) and Eq. ( 15) allows the linearization of the viscoelastic conjugate variables, {a a a ve i } P i=1 , as

δa δa δa ve i = -G G G σ : : : δ δ δ -g g g σϑ δϑ + ∑ j G G G ve i ve j : : : δ δ δ ve j + g g g ve i r δr , ( 37
) G G G ve i ve j = [1 -d] C C C e if i = j, G G G ve i ve i = [1 -d] (C C C e + C C C ve i ) if i = j, g g g ve i r = -[1 -d] -1 Λ d a a a ve i + [1 -d] C C C e : : : Λ Λ Λ r .
Thus, both foregoing equations (Eq. ( 36) and Eq. ( 37)) can be combined to express increments of the viscoelastic criteria, {δΦ δΦ δΦ ve i } P i=1 , as functions of the state variables. Their expressions can be written in a more compact way as follows δΦ δΦ δΦ ve i = K K K ve i : : : δ δ δ + k k k ve i ϑ δϑ + ∑ j K K K ve i ve j : : : δ δ δ ve j + k k k ve i r δr , (38)

K K K ve i = -D D D -1 ve i : : : C C C e , k k k ve i ϑ = D D D -1 ve i : : : C C C e : : : α α α, k k k ve i r = D D D -1 ve i : : : C C C e : : : Λ Λ Λ r , K K K ve i ve j = D D D -1 ve i : : : C C C e if i = j, K K K ve i ve i = D D D -1 ve i : : : (C C C e + C C C ve i ) + ∆t -1 I I I if i = j.
In a similar manner, the usage of Eq. ( 28) and Eq. ( 34) in combination with Eq. ( 36) allows to write the increment of the viscoplastic criterion, δΦ r , in the following form5 δΦ r = k k k r : : :

δ δ δ + k r ϑ δϑ + ∑ i k k k r ve i : : : δ δ δ ve i + k r r δr , ( 39
)
k k k r = [1 -d] κ 1 Λ Λ Λ r : : : C C C e , k k r ve i = -[1 -d] κ 1 Λ Λ Λ r : : : C C C e , k r ϑ = -[1 -d] κ 1 Λ Λ Λ r : : : C C C e : : : α α α -κ 1 F ,rϑ -κ 1 Σ y,ϑ -κ 2 , k r r = -[1 -d] κ 1 Λ Λ Λ r : : : C C C e : : : Λ Λ Λ r -κ 1 F ,r r -κ 1 Σ y,r -∆t -1 , κ 1 = p vp 2Z p vp vp f + 1 + f | f | , κ 2 = p vp Z p vp +1 vp f + Z vp,ϑ .
Defining a Fischer-Burmeister type of criterion for the viscoplastic part of the model as

Φ r := (Φ r ) 2 + (∆r) 2 + Φ r -∆r, (40)
its linearisation using the convex cutting plane approach allows to write δΦ r = k k k r : : : δ δ δ + k r ϑ δϑ + ∑ i k k k r ve i : : : δ δ δ ve i + k r r δr , ( 41)

k k k r = ω r 1 k k k r , k r ϑ = ω r 1 k r ϑ , k k k r ve i = ω r 1 k k k r ve i , k r r = ω r 1 k r r + ω r 2 , ω r 1 = Φ r (Φ r ) 2 + (∆r) 2 + 1, ω r 2 = ∆r (Φ r ) 2 + (∆r) 2 -1.
Both expressions of the viscoelastic criteria increment (Eq. ( 38)) and the Fisher-Burmeister type criterion increment (Eq. ( 41)), form a corrector linear system of equations from which all unknown internal variable increments can be identified. Such a corrector system is expressed in matrix notation as follows

Φ Φ Φ + δΦ δΦ δΦ = 0 0 0 ; δΦ δΦ δΦ = X X X δ δ δ + X X X ϑ δϑ + X X X ג δג δג δ,ג (42)
where all X X X-terms in Eq. ( 42) are defined in the appendix.

Hence, all inelastic strains and damage variable are updated at the increment (n + 1)(m + 1)

following the set of relations below

ve i = ve i + δ δ δ ve i ; vp = vp + Λ Λ Λ r δr ; d = d + Λ d δr. ( 43
)
The stress tensor is then adjusted iteratively until the preset viscoelastic-viscoplastic-damage tolerance is reached (i.e. Φ Φ Φ

(n+1)(m+1) ve i and Φ (n+1)(m+1) r ≤ Tolerance).
It is worth noting that the viscoelastic correction requires only one Newton Raphson iteration since all considered viscoelasticities are linear. However, a typical full viscoelastic-viscoplastic-damage step needs many corrections and requires more computational time to update all internal state variables.

Tangent modulus

In addition to the stress and internal state variables computation, the FE analysis requires the estimation of appropriate thermomechanical tangent moduli. The latter are computed at a specific time step. In this regard, any increment that appears in the calculations is in the same spirit with the virtual increments of the constitutive law algorithm. Thus, a given quantity z z z is considered to be updated as z z z → z z z + δz z z.

More specifically, four tangent moduli need to be identified (partial derivatives of the total stress tensor, σ σ σ, and heat sources, h s , with respect to both the total strain tensor, , and temperature, ϑ). To do so, it is assumed that δ δ δ = 0 0 0, δϑ = 0 and δΦ δΦ δΦ = 0 0 0. This allows to write

δג δג δג = L L L ג + L L L ג ϑ; L L L ג = -X X X -1 X X X ; L L L גϑ = -X X X -1 X X X ϑ . ( 44
)
From Eq. ( 44), increments of all internal variables (inelastic strains and damage variable) can be written, as functions of increments of the total strain and temperature δ δ δ ve i = L L L ve i : : : δ δ δ + l l l ve i ϑ δϑ, δd = (Λ d l l l r ) : : : δ δ δ + (Λ d l rϑ ) δϑ, (45) δr = l l l r : : : δ δ δ + l rϑ δϑ, δ δ δ vp = (Λ Λ Λ r ⊗ l l l r ) : : :

δ δ δ + (Λ Λ Λ r l rϑ ) δϑ. ( 46
)
Substituting Eq. ( 45) and Eq. ( 46) into Eq. ( 36) and Eq. ( 37) yields the following linearized forms of both the stress tensor and viscoelastic conjugate variables, σ σ σ and {δa δa δa ve i } P i=1 , expressed as functions of total strain and temperature increments δσ δσ δσ = M M M σ : : : δ δ δ + m m m σϑ δϑ, δa δa δa ve i = M M M ve i : : : δ δ δ + m m m ve i ϑ δϑ, (47)

M M M σ = G G G σ -G G G σ : : : ∑ i L L L ve i + g g g σr ⊗ l l l r , m m m σϑ = g g g σϑ -g g g σϑ : : : ∑ i l l l ve i ϑ + g g g σr l rϑ , ( 48
) M M M ve i = -G G G σ + ∑ j G G G ve i ve j : : : L L L ve j + g g g ve i r ⊗ l l l r , m m m ve i ϑ = -g g g σϑ + ∑ j
G G G ve i ve j : : : l l l ve j ϑ + g g g ve i r l rϑ .

Both M M M σ and m m m σϑ tensors given in Eq. ( 48) represent the first tangent moduli of the constitutive model. The last two moduli are derived from the linearized form of the heat diffusion equation.

Hence, substituting Eq. ( 47) into Eq. ( 13) and Eq. ( 18) allows to write δs = m m m s : : :

δ δ δ + m s ϑ δϑ, δa d = m m m d : : : δ δ δ + m d ϑ δϑ, ( 49
)
m m m s = α α α : : : M M M σ -F , rϑ l l l r , m s ϑ = α α α : : : m m m σ ϑ + c o ϑ -F ,ϑϑ -F , rϑ l r ϑ , m m m d = - e : : : M M M σ [1 -d] - 1 2 e : : : σ σ σ [1 -d] 2 Λ d l l l r -∑ i ve i : : : C C C ve i : : : L L L ve i , m d ϑ = - e : : : m m m σ ϑ [1 -d] - 1 2 e : : : σ σ σ [1 -d] 2 Λ d l r ϑ -∑ i ve i : : : C C C ve i : : : l l l ve i ϑ ,
Further, the usage of foregoing linearized forms and Eq. ( 20) yields to the following formulation of the intrinsic dissipation increment

δΩ int = m m m Ω : : : δ δ δ + m Ω ϑ δϑ, ( 50
)
m m m Ω = 1 ∆t o o o r + ∑ i o o o ve i + o o o d , m Ω ϑ = 1 ∆t o r ϑ + ∑ i o ve i ϑ + o d ϑ , o o o r = [σ : Λ σ : Λ σ : Λ r ] l l l r + ∆ ∆ ∆ vp : : : M M M σ -[F , r + ∆rF , rr ] l l l r , o o o ve i = -a a a ve i : : : L L L ve i -∆ ∆ ∆ ve i : : : M M M ve i , o o o d = -a d Λ d l l l r -∆d m m m d , o rϑ = [σ : Λ σ : Λ σ : Λ r ] l rϑ + ∆ ∆ ∆ vp : : : m m m σϑ -[F , r + ∆rF , rr ] l rϑ -F , rϑ ∆r,
o ve i ϑ = -a a a ve i : : : l l l ve i ϑ -∆ ∆ ∆ ve i : : : m m m ve i ϑ , o d ϑ = -a d Λ d l r ϑ -∆d m d ϑ .
Finally, both tangent moduli (denoted below by m m m h and m hϑ ) are derived from the linearized form of the heat diffusion equation (see expression 11 of Table 2), expressed such that

δh s = m m m h : : : δ δ δ + m hϑ δϑ, ( 51
)
m m m h = - ϑ ∆t m m m s + m m m , m hϑ = - ∆s ∆t - ϑ ∆t m sϑ + m Ωϑ . ( 52
)

Model validation

The presented thermo-viscoelastic-viscoplastic-damage model is hereafter calibrated using a set of quasi-static experimental data derived from tests conducted on a wet PA6.6. In addition to existing tests performed in the author's previous work [START_REF] Praud | Phenomenological multimechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation[END_REF], monotonic tensile and stress-controlled tensile-tensile tests have been carried out at room temperature using a hydro-controlled MTS testing machine. The tensile tests have been conducted at 3 different displacement rates to highlight the viscoelastic effects by investigating the stiffness sensitivity to the loading rate, while the stress-controlled cyclic tests were conducted to emphasize the damage kinetics throughout the loading and also to highlight the dissipative and thermomechanical effects arising from the inelastic deformation. During cyclic tests, the surface temperature of the sample was recorded using an infrared camera (Flir Titanium with InSb detector). The corresponding mechanical, thermal and energy responses are presented in Figure 2. In Fig. 2, the deformation energy rate per cycle, W de f , was computed using the following expression

W de f = Γ f l v gp F du s , Γ = L o /2 -L o /2 l -1 s (y)dy/ L/2 -L/2 l -1 s (y)dy, (53)
where f l stands for the loading frequency, while F, u s and v gp are the applied force, the current length and volume of the sample gage part, respectively. In Eq. ( 53), the term Γ represents a geometrical parameter, whilst L o and l s stand for the length of the gage section and the width of the sample.

The mean intrinsic dissipation was computed using the following heat diffusion equation

W d = ρc f l k f l (k-1) f l θ + ϑ τ 0D hl dt, ( 54
)
where k stands for the cycle number and τ 0D hl the time constant of the heat exchanges between the sample and surroundings. For deep details on the energy behaviour of the PA6.6 material, readers are directed to [START_REF] Benaarbia | Influence of relative humidity and loading frequency on the pa6.6 thermomechanical cyclic behavior: Part ii. energy aspects[END_REF].

A rheological representation of the one-dimensional scheme of the proposed model is illustrated in Figure 3. The model is split into many components in order to consider all recoverable and irrecoverable mechanisms (e.g. coupled damage elasticities, thermomechanical coupling, rate-dependent viscoelasticities, rate-dependent viscoplasticity, etc.). The thermoelastic element is hereafter used to account for the thermomechanical coupling sources arising during the cyclic deformation of the PA6.6 material (see Figure 2b). The Kelvin-Voigt (KV) branches are employed to numerically capture the rate sensitivity of dissipative effects (which are more striking on wet PA6.6 materials compared to dry ones, see Benaarbia et al. (2014b)), whilst the viscoplastic element is introduced to model the irrecoverable viscoplastic strain as well as the non-linear cyclic softening observed in the PA6.6 material behavior (see. [START_REF] Praud | Phenomenological multimechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation[END_REF]). Parameters related to internal heat production and heat exchange (i.e. c, α and D, where the later stands for thermal diffusivity) were measured experimentally using a hot disk thermal analyzer and Netzsch DIL 402 dilatometer. These are set in the following as 1638 J.kg -1 . • C -1 , 98×10 -6 • C -1 and 0.343 W.m -1 . • C -1 , respectively. The mass density was determined using a pycnometer and is set as 1120 kg.m -3 , while the Poisson's ratio of elasticity and viscoelasticities is considered equal to 0.3.

It is worth noting that under the linear thermoelasticity framework and for isotropic materials, the amplitude of temperature variations, denoted hereafter by ∆ϑ the , varies linearly with the amplitude of applied stress, denoted by ∆σ [START_REF] Boulanger | Calorimetric analysis of dissipative and thermoelastic effects associated with the fatigue behavior of steels[END_REF]). Thus, one calibrates the thermal expansion coefficient using the following thermoelastic relation

∆ϑ the = α ϑ o ∆σ ρ c . ( 55
)
The constitutive model uses a set of 16 material parameters to describe the cyclic thermomechanical response of the PA6.6 material. In this work, 4 KV branches have been considered to allow capturing the viscoelastic effects at the frequencies used in the designed experimental protocol.

These parameters are determined through an optimization process by fine-tuning the model responses to the available experimental results. A Levenberg-Marquardt algorithm was used for the purposes of material parameters optimization [START_REF] Meraghni | Parameter identification of a thermodynamic model for superelastic shape memory alloys using analytical calculation of the sensitivity matrix[END_REF]. A single objective function that compares experimental and predicted responses for every data set (below referred by the index j) was applied to check the accuracy of fitting for each of the experimental data sets.

Thus, a scaled cost function was evaluated by combining all objective function values for each of the data sets, such that

F c = ∑ j ∑ i p j x exp ij -x num ij 2 / ∑ j ∑ i x exp ij 2 , ( 56
)
where the objective function, F c , is summed of j experiments (here we talk about monotonic tensile and cyclic tensile-tensile tests), while x exp ij and x num ij stand for the experimental and numerical stress (and/or strain) values. The weighting term p j is used to adjust contribution from test result sets and also to offset the magnitude difference between each data set (without this normalisation, the higher curves with higher numerical values would dominate the scoring).

To restrict the search for each parameter, penalty reasonable bounds were used. Sufficiently large bounds were placed on each of the model coefficients during optimization so that the solution was not overly constrained by the initial coarse guesses.

Determination of material parameters was achieved by optimizing the trial model parameters against the cyclic tensile-tensile mechanical results coupled with the self-heating data. Proper calibration of the model requires several types of experimental data. In authors' previous work using a material similar to the investigated PA6.6 [START_REF] Praud | Phenomenological multimechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation[END_REF], the parameter identification was performed using four types of tests: monotonic, incremental load/unload tensile, incremental load/relaxation/unload/relaxation tensile and cyclic tensile. Here, the scope is not to present a full calibration procedure, but to mainly focus on the thermomechanical coupling aspects.

Approximate orders of magnitude for parameters in the 4 Kelvin-Voigt elements were determined in [START_REF] Praud | Multi-scale modelling of thermoplastic-based woven composites, cyclic and time-dependent behaviour[END_REF]. Related viscoplastic and damage material constants were firstly fitted to the monotonic data and then recalibrated using the cyclic tensile-tensile mechanical data.

The parameter optimization procedure is considered completed once the stationary conditions are reached. In our identification study, the cost function for the final set of parameters obtains the value of 1.87.10 -6 . Optimized material parameters are presented in Table 3.

Comparative plots for the three loading rates selected for the monotonic tensile test and for the hysteretic response are presented in Fig. 4a-b. In all cases, a good level of accordance is observed between predicted and experimental results, with loading rate dependencies and viscoplastic effects appropriately represented. Particular emphasis should also be placed on the ability of the model to capture ratcheting behavior of the PA6.6 material. The novelty of the proposed model is also its ability to regenerate all energy terms involved in the energy

balance. An excellent level of agreement is observed between predicted and experimental energy results (see Fig. 5 and Fig. 2c); the model can numerically capture both dissipative and storage mechanisms, with a value of Taylor-Quinney coefficient, β di f f , estimated at 10% at the last investigated cycles (see Fig. 6).

Figure 4: The prediction of mechanical responses using the optimised material parameters (given in Table 3) under both a) monotonic tensile tests and b) cyclic tensile-tensile test. To illustrate the efficiency of the proposed numerical scheme, a three cycle stress driven loading has been studied and analyzed in terms of error estimation per algorithmic iteration.

Specifically, a material point has been subjected to adiabatic conditions and stress controlled cyclic loading with a rate of 50MPa/s. At three distinct positions, indicated in Figure 7, the algorithm converged after three iterations. The error per iteration is described by the relation

error = δV V V ∆V V V , ( 57
)
where δ and ∆ denote the iteration and time increments respectively, . stands for the norm of a vector, and V V V is a vector containing the strain, stress and temperature fields, i.e.

V V V =       σ σ σ ϑ       . ( 58
)
As shown in the tables of Figure 7, at every position the convergence rate is quadratic. This result highlights the good performance of the numerical algorithm. [-] Figure 7: Error estimation per iteration highlighting the efficiency of the proposed numerical scheme.

Finite element application of the constitutive model

The capability of the proposed model is further checked through finite element analyses conducted on a Meuwissen structure (see Figure 8) using the optimised material parameters.

The discretized equations of subsections 5.1, 5.2 and the appendix have been implemented in a homemade user material subroutine (UMAT), written in FORTRAN. The latter has been developed and integrated into the Finite Element commercial software ABAQUS/Standard [START_REF] Simulia | Abaqus/explicit user's manual[END_REF]. The UMAT takes into account all related thermo-viscoelastic-viscoplasticdamage coupling terms, as described in the model. Hence, the rate sensitivity of energy responses is investigated by performing a parametric study under tensile loading conditions conducted at different loading rates (i.e. 0.0021GN.s -1 , 0.021GN.s -1 , 0.21GN.s -1 and 2.1GN.s -1 ). The Meuwissen structure is clamped at the left side and subjected to tensile loading at the right side (the dashed red boxes in Figure 8 highlight the location of the loading machine grips.). The mesh process was performed using 4563 eight-node trilinear displacement and temperature elements, C3D8T in Abaqus notation. Figure 8 shows some features of the finite element model.

Figure 8: a) Geometry of the Meuwissen structure (All dimensions are in mm) and Mesh configuration used in the numerical simulations (three elements along the thickness). P1 is a measurement point where intensive localizations will take place.

Mechanical and self heating responses

The mechanical and self-heating responses of the Meuwissen structure at the measurement point P1 were first evaluated by checking the ability of the constitutive model to describe the material rate sensitivity (see Figure 9 where the four selected loading rates were applied). The results first show that the maximum reachable strain (resp. self heating) increases when the loading rate decreases (resp. increases). These findings are in accordance with the experimental observations conducted on the PA6.6 material where both the mean strain and self-heating have been observed to increase with increasing loading rates (the dissipative phenomena (in terms of magnitude) are hardly visible for the high loading rates, see [START_REF] Benaarbia | Influence of relative humidity and loading frequency on the pa6.6 thermomechanical cyclic behavior: Part ii. energy aspects[END_REF]). The influence of the loading rate on the dissipative effects from the very beginning of the loading is a clear hint of short-term viscoelastic effects. These effects are much more striking on SCPs [START_REF] Khan | Thermo-mechanical response of nylon 101 under uniaxial and multi-axial loadings: Part i, experimental results over wide ranges of temperatures and strain rates[END_REF][START_REF] Farrokh | A strain rate dependent yield criterion for isotropic polymers: Low to high rates of loading[END_REF], especially on wet polymers.

Taylor-Quinney and damage responses

The change in energy characteristics of the Meuwissen structure was evaluated by tracking the spatio-temporal distribution of β di f f throughout the loading. The goal is to track the possible development of energy localization zones. Figure 10 shows typical β di f f fields generated using the user material routine (UMAT) developed within the current work. All maps were computed at different loading stages (referred here by t 1 < t 2 < t 3 < t 4 ) and the values of the featured β di f f were depicted and graduated on a color bar. The analysis was performed using calorimetric and kinematic investigations based essentially on computing heat sources and inelastic energy involved within the Meuwissen structure. It should be noted that the spatial distributions of these energy maps are not uniform and vary markedly. The spatial heterogeneities start to be clearly visible as long as the loading increases. The propagation of the dissipative mechanisms seems to concentrate in the gage part, notably at the dissymmetric edges of the structure where the stored energy ratio, β di f f , is low (with a maximum value typically around 6%). the loading with a maximum damage of 25% observed in both out-plane faces (near to the measurement point P1). These spatial propagations have been observed for all loading rates selected for the current parametric study. Figure 12a-b show the damage variable evolving with respect to both the total strain and temperature. It can be noted from Figure 12a that the damage magnitude is strongly affected by the loading rate, whilst the damage trend seems not to be influenced by the time effects. Figure 12b indicates that both values and trend of the damage have strong sensitivity to temperature variations. Figure 13a shows the evolution of the differential Taylor-Quinney, β di f f , with respect to the total strain, , for the selected loading rates at the measurement point P1. It can be observed that β di f f strongly varies as deformation proceeds; it starts to increase at the very beginning of loading and then seems to reach a constant state notably for the low loading rates, and finally decreases and re-increases later for the rest of the loading. The rate sensitivity of β di f f is largely visible with higher values reported for the lowest loading rates. Moreover, the β di f f kinetics with respect to the damage variable d are depicted in Figure 13b. One can note, during the very beginning of loading, a sharp decrease in values of β di f f with increasing damage till reaching a certain minimum value after which β di f f increases with increasing damage. These observations lead us to infer that β di f f is not only confined to damage-induced mechanisms (Other mechanisms can contribute to the storage of energy).

Concluding comments

The The numerical findings have also shown a good accordance with the experimental observations conducted on the PA6.6 material where the mean strain, self-heating and Taylor-Quinney coefficient have been observed to increase with increasing loading rates.

The advantage of the validated implicit formulation is the determination of the four tangent moduli (partial derivatives of the stress tensor and heat sources with respect to both strain tensor and temperature). This offers an opportunity for the implementation of the proposed approach into micromechanics scheme that allows modeling of composite materials using an incremental multiscale strategy where the polymer matrix behaves as a thermo-viscoelastic-viscoplastic material exhibiting ductile damage.

Another outlook of the proposed work is to investigate the material for extended humidity rates and temperature conditions. Indeed, water plasticizing effects are recognized for their marked influence on the glass transition temperature (GTT) [START_REF] Benaarbia | Thermomechanical analysis of the onset of strain concentration zones in wet polyamide 6.6 subjected to cyclic loading[END_REF]; the GTT decreases as the relative humidity increases which means that the polymer shows rubber effects under standard usage temperatures. Several approaches have been developed in the past to describe the thermomechanical rubber effects and the so-called thermoelastic inversion (competition between standard and entropic thermoelasticity [START_REF] Joule | On some thermodynamic properties of solids[END_REF]). It would be particularly interesting to extend the model to cover the entropic effects above the glass transition temperature, which regards many applications for the automotive industry.

Figure 1 :

 1 Figure 1: Kinetics of the Taylor-Quinney coefficient for a PA6.6 material under different moisture contents (RH0, RH50 and RH80). The results are extracted from tensile-tensile cyclic tests conducted at both loading rates (a) 1Hz and (b) 10Hz for a load ratio set at 0.1.

  r ] : : : ˙ -∑ P i=1 a a a ve i : : : ˙ ve i a: VE -a a a vp : : : ˙ vp -R ṙ intrinsic dissipation is generated by the viscoelastic (term a), viscoplastic (term b)

3

  The spherical and deviatoric operators I I I vol and I I I dev are given by I I I vol := 1 3 1 1 1 ⊗ 1 1 1 and I I I dev = I I I -I I I vol ; with (a a a ⊗ b b

Figure 2 :

 2 Figure 2: Thermomechanical cyclic responses of a wet PA6.6 extracted from(Benaarbia et al., 2014a[START_REF] Benaarbia | Influence of relative humidity and loading frequency on the pa6.6 thermomechanical cyclic behavior: Part ii. energy aspects[END_REF]: a) hysteresis loops highlighting the ratcheting behavior, b) self-heating response showing both the thermoelasticity and the dissipation kinetics, c) energy rate balance form indicating all energy terms (deformation energy rate (in black), intrinsic dissipation (in red) and stored energy rate (in green)). All data have been derived from the cyclic tensile-tensile test under a stress ratio of 0.1, a loading frequency of 1Hz and a relative humidity of RH50.

Figure 3 :

 3 Figure 3: 1D-rheological representation of the proposed constitutive model.

Figure 5 :

 5 Figure 5: The prediction of a) strain and b) energy responses using the optimized material parameters under cyclic tensile-tensile test.

Figure 6 :

 6 Figure 6: The predicted Taylor-Quinney coefficient using the optimized parameters during the cyclic test.

Figure 9 :

 9 Figure 9: The stress strain curves of the Meuwissen structure captured at the measurement point P1 for the selected loading rates. b) The corresponding self heating responses highlighting the rate sensitivity of the material thermal behavior.

Figure 10 :

 10 Figure 10: 3D Maps of the Taylor-Quinney coefficient, β di f f , estimated at different loading stages for a loading rate of 2.1GN.s -1 and a time increment set at 0.005s.

Figure 11

 11 Figure 11 depicts the spatio-temporal distributions of the associated damage captured at the same loading conditions for the same loading stages. The damage patterns indicate a precocious and gradual development of spatial concentration zones. These spatial heterogeneities develop preferentially close to the dissymmetric edges of the structure. They spread rapidly until they cover the half of the gage width. The mean values of damage grow continuously throughout

Figure 11 :

 11 Figure 11: 3D Maps of damage captured at different loading stages for a loading rate of 2.1GN.s -1 and a time increment of 0.005s.

Figure 12 :

 12 Figure 12: Kinetics of the damage variable with respect to both a) the total strain and b) temperature variations captured at the measurement point P1.

Figure 13 :

 13 Figure 13: Kinetics of the Taylor-Quinney coefficient, β di f f , with respect to a) the total strain and b) damage.

  challenge of the current research work was to develop a modeling framework to investigate the thermomechanical anelastic behavior of wet thermoplastic polymers. The current model was mainly designed to describe the variability of the Taylor-Quinney coefficient and the material energy balance form during both monotonic and cyclic deformation. A significant amount of attention was given to both recoverable viscoelastic, irrecoverable viscoplastic and damage effects. The proposed approach was incorporated into the frame of Thermodynamics of Irreversible Processes and Generalised Standard Materials formalism to offer a thermodynamically grounded coupling of all deformation mechanisms (e.g. storage of energy, dissipation of energy, thermomechanical coupling, etc.). The constitutive model was numerically implemented in a three-dimensional frame following an implicit formulation. The computational methodology was based on the radial return mapping algorithm using the convex cutting plane form. The model parameters have been determined following an identification procedure conducted under different test conditions (monotonic tensile tests with different loading rates and stress-controlled cyclic tensile-tensile tests). A parameteric study has then been conducted on a Meuwissen structure, with the help of the commercial finite element code Abaqus, to evaluate the capability of the proposed model to describe the material rate sensitivity and the variability of both the Taylor-Quinney coefficient and the thermomechanical damage during the loading.A fair agreement, in terms of energy balance, has been demonstrated between the numerical simulations and experimental data. The proposed model was relevant for predicting the anelastic thermomechanical behavior of the wet thermoplastic polymers considered. The influence of the loading rate on the dissipative and storage energy effects has been captured by the model.

Table 2 :

 2 Summary set of the thermo-viscoelastic-viscoplastic-damage constitutive model.

	State laws:
	1. The total stress tensor	σ σ σ = [1 -d] C C C e : : : e .
	2. The total entropy	s = c o ln	ϑ ϑ o	+ s o + α α α : : : σ σ σ.
	3. The viscoelastic conjugate variables 4. The viscoplastic conjugate variable	a a a ve i = [1 -d] C C C ve i : : : ve i -σ σ σ with i = 1, ..., P. a a a vp = -σ σ σ.
	5. The drag stress	R = ∂ r F (r) .
	6. The damage conjugate variable	a d = -	1 2	e : : : σ σ σ [1 -d]	-	1 2

P ∑ i=1 ve i : : C C C ve i : : : ve i .

Table 3 :

 3 A summary of the thermo-viscoelastic-viscoplastic-damage material constants (determined by applying the optimization procedure to the monotonic tensile tests and cyclic tensile-tensile data).

	Thermomechanical feature Model parameter	Symbol Unit	Optimised value
	Thermoelasticity
	Initial Young Modulus	E o	MPa	2731.14
	Poisson's ratio	ν	-	0.3
	Thermal expansion	α	K -1	98.10 -4
	Viscoelasticities
	Viscosity -KV1	η ve 1	MPa.s	3485.98
	Stiffness modulus -KV1	E ve 1	MPa	9751.44
	Viscosity -KV2	η ve 2	MPa.s	128516.91
	Stiffness modulus -KV2	E ve 2	MPa	19125.64
	Viscosity -KV3	η ve 3	MPa.s	196773.14
	Stiffness modulus -KV3	E ve 3	MPa	30855.24
	Viscosity -KV4	η ve 4	MPa.s	870047.82
	Stiffness modulus -KV4	E ve 4	MPa	6771.25
	Viscoplasticity
	Yield stress	Σ y	MPa	47.64
	Viscoplastic resistance	Z vp	MPa.s p vp 45.86
	Viscoplastic exponent	p vp	-	0.072
	Isotropic hardening coefficient h 01	MPa	1302.71
	Isotropic hardening exponent	n vp	-	0.8
	Damage
	Damage coefficient	K	MPa	20.03
	Damage exponent	n d	-	-0.861

The inverse of a symmetric fourth-order tensor A A A is defined as the fourth-order tensor [A A A] -1 , for whichA ijmn [A mnkl ] -1 = I ijkl = 1 2 [δ ik δ jl + δ il δ jk ],where δ ij stands for the Kronecker delta.

While in the current model Σ y and Z vp are constants, the numerical algorithm is presented in a general fashion considering that both can be functions of r and θ.

Acknowledgements

The authors gratefully Acknowledge Solvay Engineering Plastics for providing material data and specimens. They also would like to warmly thank a very present partner, Dr. Gilles Robert, for his effective tips that significantly improved the current research work.

Appendix. The VE-VP corrector system

The corrector linear system from which the unknown increments of internal variables can be determined is expressed such that

where

,