Léo Thiercelin

Francis Praud

Fodil Meraghni
email: fodil.meraghni@ensam.eu

Eric Fleury

Thermally-activated hardening recovery in viscoplastic materials with kinematic hardening at high temperatures

In this work, a new constitutive model is proposed to describe the thermally-activated hardening recovery mechanism in metallic materials. This model takes up the concept of hardening recovery variable, which is extended to the case of kinematic hardening within a thermo-viscoplastic formulation including high temperature dependencies for both elastic and viscoplastic properties. The model is identified for AISI 316L austenitic stainless steel using experimental data from uni-axial tests conducted over a wide range of temperature (from room temperature to 1273 K) and at two different strain rates (2.5 × 10 -4 and 2.5 × 10 -3 s -1 ). Validation is further achieved with rather good agreements by comparing the simulated responses with additional experimental data where the material is subjected to complex thermomechanical loading paths while other examples are presented to provide a better insight of the model and to illustrate its predictive capabilities.

Introduction

It is well-established that the properties of metals and alloys strongly depend on temperature. It is commonly established that the solicitations at high temperatures can have detrimental effects due to the reduction of strength properties, or the development of residual stresses, possibly resulting in material failure [START_REF] Withers | Residual stress and its role in failure[END_REF]. However, in some other cases like in the annealing process or in hot metal forming, thermal exposure may be beneficial by helping the material to relax stresses (Totten et al., 2002;[START_REF] Angkurarach | Effects of high-temperature deep rolling on fatigue, work hardening, and residual stress relaxation of martensitic stainless steel AISI 420[END_REF] and by promoting ductility [START_REF] Humpreys | Recrystallization and Related Annealing Phenomena[END_REF]. Whatever the purpose, to meet industrial challenges and needs, the thermomechanical behaviour of metallic materials must be well understood to accurately predict their response.

A wide variety of phenomenological and constitutive models have been developed over the past decades to address plastic and viscoplastic response of materials upon monotonic and cyclic loadings [START_REF] Chaboche | A unified constitutive model for cyclic viscoplasticity and its applications to various stainless steels[END_REF][START_REF] Lemaitre | Mechanics of Solid Materials[END_REF][START_REF] Chaboche | A review of some plasticity and viscoplasticity constitutive theories[END_REF][START_REF] Besson | Non-Linear Mechanics of Materials[END_REF]Karvan and Varvani-Farahani, 2019b,a;[START_REF] Shekarian | Ratcheting prediction at the notch root of steel samples over asymmetric loading cycles[END_REF][START_REF] Karvan | Viscoplastic ratcheting response of materials under step-loading conditions at various cyclic stress levels[END_REF]. These models introduced different types of hardening, like isotropic, kinematic or combined with temperature-dependent parameters to capture the rate-and temperature-dependent non-linear behaviour of metallic materials.

The temperature dependence of the material parameters reflects microstructural changes and the activation of various physical phenomena taking place at the microscopic scale. Among them, the hardening recovery is a particular phenomenon that enables the material to partially or fully restore its initial properties upon thermal exposures, as the result of dislocations rearrangement and annihilation, grain recrystallization and their subsequent growth [START_REF] Humphreys | A unified theory of recovery, recrystallization and grain growth, based on the stability and growth of cellular microstructures-I. The basic model[END_REF][START_REF] Humpreys | Recrystallization and Related Annealing Phenomena[END_REF], as illustrated in Fig. 1.

Several phenomenological models have been proposed to integrate the annealing under temperature exposures. Commonly used models consider the existence of an activation temperature from which the hardening recovery is complete or initiated. In this case, the hardening variables are artificially reset to zero (Muránsky et al., 2012a,b;[START_REF] Deng | Influence of material model on prediction accuracy of welding residual stress in an austenitic stainless steel multi-pass buttwelded joint[END_REF]. In other models, two annealing stages are considered: a starting temperature for partial recovery and another one, beyond which the recovery is total and instantaneous. Despite predictive results [START_REF] Hendili | HSNV140 -thermo-plasticité avec restauration d'écrouissage : essai de dilatométrie bloquée[END_REF][START_REF] Depradeux | Influence of viscoplasticity, hardening, and annealing effects during the welding of a three-pass slot weld (NET-TG4 round robin)[END_REF], these approaches display a lack of physical consistency, since they only consider the recovery mechanism when a certain temperature is reached without accounting for the time of exposure.

Recrystallization models based on Johnson-Mehl-Avrami-Kolmogrov laws have also been developed [START_REF] Avrami | Kinetics of phase change. I General theory[END_REF][START_REF] Avrami | Kinetics of phase change. II Transformation-time relations for random distribution of nuclei[END_REF][START_REF] Fanfoni | The Johnson-Mehl-Avrami-Kolmogorov model: A brief review[END_REF]. These models are often limited to simple heat treatments like isothermal conditions or anisothermal with constant heating rates, and are hence not applicable to practical real cases. Nevertheless, [START_REF] Suo | Influence of strain hardening and annealing effect on the prediction of welding residual stresses in a thick-wall 316 stainless steel butt-welded pipe joint[END_REF] have recently proposed a time-and temperature-dependent corrective factor directly derived from the Avrami's equation to correct the kinematic hardening variable of their model.

One can also mention the dislocation-based models that represent the competition between annihilation and creation of dislocations under the combined effects of stress and temperature [START_REF] Mecking | Kinetics of flow and strain-hardening[END_REF][START_REF] Estrin | A unified phenomenological description of work hardening and creep based on one-parameter models[END_REF][START_REF] Kocks | Physics and phenomenology of strain hardening: The FCC case[END_REF][START_REF] Blaizot | Constitutive model for nickel alloy 690 (Inconel 690) at various strain rates and temperatures[END_REF][START_REF] Lin | A unified constitutive model based on dislocation density for an Al-Zn-Mg-Cu alloy at time-variant hot deformation conditions[END_REF][START_REF] Yuan | Thermomechanical behavior of laser metal deposited inconel 718 superalloy over a wide range of temperature and strain rate: testing and constitutive modeling[END_REF][START_REF] Voyiadjis | Effect of dislocation density evolution on the thermomechanical response of metals with different crystal structures at low and high strain rates and temperatures[END_REF]Voyiadjis et al., 2019). Since these models involve dislocation physics, they rather well describe the recovery mechanism. However, the identification of those models remains quite cumbersome and requires non-conventional experimental techniques to quantify the dislocation density and to follow their kinetics during the tests, e.g., Transmission Electron Microscopy. Moreover, the hardening recovery is not only related to the annihilation of dislocations, particularly in metals with low stacking fault energy, such as austenitic stainless steel, since the recrystallization and the dislocations annihilation seem to occur simultaneously [START_REF] Humpreys | Recrystallization and Related Annealing Phenomena[END_REF].

In an alternative approach, some of the authors of this paper has previously introduced the concept of hardening recovery variable within a classical elastoplastic formulation with isotropic hardening [START_REF] Mouelle | Thermally-activated hardening recovery of thermo-elasto-plastic metals during annealing: constitutive modeling for the simulation of welding process[END_REF]. The main idea behind this concept is to consider an internal state variable, referred to as the recovery variable, that gradually counteracts the effect of the hardening through a temperature-dependent evolution law. Although this model was able to capture some interesting tendencies, some features are still clearly missing to properly reproduce the material responses at high temperatures [START_REF] Mouelle | Caractérisation et modélisation d'assemblages de tôles en acier inoxydable austénitique brasées et soudées pour des applications d'échangeurs de chaleur[END_REF]. This was mainly caused by the fact that, the hardening recovery was only formulated in the case of isotropic hardening without any other time-or temperature-dependencies with regards to the elasto-plastic properties. Indeed, isotropic hardening is not always appropriate to well describe the material behaviour upon alternate tension-compression. Such loading conditions are common in metal-forming applications like in multi-pass welding. Moreover, accounting for time-and temperature-dependencies appears to be necessary to well describe the material behaviour at moderate and high temperatures, even when the hardening recovery is not yet active.

Therefore, to address these issues, a new constitutive model is proposed in this work. The latter takes up the concept of the hardening recovery variable, but extends it to the case of kinematic hardening. The model is formulated considering an elasto-viscoplastic formulation since it is well established that metallic materials exhibit viscous effects at high temperatures. Furthermore, temperature dependencies are set in the elastic and viscoplastic properties to account for the prevailing thermal influences on the material response. The model is then identified for AISI 316L austenitic stainless steel using experimental data from [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF] and [START_REF] Hendili | HSNV140 -thermo-plasticité avec restauration d'écrouissage : essai de dilatométrie bloquée[END_REF], who conducted uni-axial tension tests over a wide range of temperature (from room temperature to 1273 K) and at two different strain rates (2.5 × 10 -4 and 2.5 × 10 -3 s -1 ). Validation is further achieved with rather good agreements by comparing the simulated responses with additional experimental data [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF] where the material is subjected to complex thermomechanical loading paths including tension-relaxation and fully constrained dilatation tests. Other examples are presented to provide a better insight of the model, as well as to illustrate the kinetic of the hardening recovery in relation with the temperature and the time of exposure.

In this work, the following notation is adopted: bold and blackboard symbols denote second and fourth order tensors, respectively, whereas other symbols are scalar quantities. The product operators ∶ and ⊗ depict the twice contracted and the dyadic products, respectively, such as:

[𝑨 ∶ 𝑩] = 𝐴 𝑖𝑗 𝐵 𝑖𝑗 , [A ∶ 𝑩] 𝑖𝑗 = 𝐴 𝑖𝑗𝑘𝑙 𝐵 𝑘𝑙 , [𝑨 ⊗ 𝑩] 𝑖𝑗𝑘𝑙 = 𝐴 𝑖𝑗 𝐵 𝑘𝑙 .
Moreover, all the second order tensors are symmetric (𝐴 𝑖𝑗 = 𝐴 𝑗𝑖 ) and all the fourth order tensors have at least the minor symmetries (𝐴 𝑖𝑗𝑘𝑙 = 𝐴 𝑗𝑖𝑘𝑙 = 𝐴 𝑖𝑗𝑙𝑘 ). Consequently, they can be respectively reduced to 6 × 1 and 6 × 6 matrices according to the Voigt notation. 𝑰 and I stand for the second and fourth order identity tensors, respectively. The positive part of a scalar quantity 𝑎 is denoted by

⟨ 𝑎 ⟩ + , such that ⟨ 𝑎 ⟩ + = 𝑎 if 𝑎 ≥ 0 or ⟨ 𝑎 ⟩ + = 0 if 𝑎 < 0.
The operators hyd(𝝈) and 𝐃𝐞𝐯(𝝈) designate the hydrostatic pressure and the deviatoric part of a stress tensor 𝝈, respectively:

hyd(𝝈) = 1 3 tr(𝝈), 𝐃𝐞𝐯(𝝈) = 𝝈 -hyd(𝝈)𝑰,
while eq(𝝈) is the equivalent von Mises stress:

eq(𝝈) = √ 3 2 𝐃𝐞𝐯(𝝈) ∶ 𝐃𝐞𝐯(𝝈).

Constitutive equations and thermodynamic framework

State laws

The proposed constitutive equations are based on a viscoplastic model with kinematic hardening [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF]. Accordingly, the energetic state of the material is described by the following form of the Helmholtz free energy potential:

𝜌𝜓(𝜺, 𝑇 , 𝜺 𝑝 , 𝝃 𝑝 , 𝝃 𝑟 ) = 1 2 ( 𝜺 -𝜺 𝑝 -𝜺 𝑡ℎ (𝑇 ) ) ∶ C(𝑇 ) ∶ ( 𝜺 -𝜺 𝑝 -𝜺 𝑡ℎ (𝑇 ) ) + 𝑐(𝑇 ) 2 (𝝃 𝑝 -𝝃 𝑟 ) ∶ (𝝃 𝑝 -𝝃 𝑟 ) + 𝜌𝜓 𝑡ℎ (𝑇 ), (1)
which involves the total strain 𝜺 and the absolute temperature 𝑇 as observable state variables. Furthermore, the Helmholtz free energy also depends on a set of internal state variables, among which are the plastic strain 𝜺 𝑝 , the back strain 𝝃 𝑝 and the back strain recovery variable 𝝃 𝑟 , simply referred to as recovery variable.

In (1), 𝜌 denotes the mass density considered as constant under the small strains assumption. C(𝑇 ) is a temperature-dependent fourth order material stiffness tensor defined as:

C(𝑇 ) = 𝑔 𝑒 (𝑇 ) C 0 , ( 2
)
where C 0 is the stiffness tensor at the reference temperature 𝑇 0 , classically defined for bulk isotropic materials by 𝐸 0 and 𝜈 0 , the Young modulus and the Poisson ratio, respectively. Note that, in this work, the reference temperature is the room temperature such that 𝑇 0 = 293 K. As written in (2), the elastic stiffness at a given temperature is connected to the one at the reference temperature through a shift function 𝑔 𝑒 (𝑇 ) that will be defined later in Section 3. This function is equal to 1 at the reference temperature and tends to 0 for high temperatures. Similarly to the elastic stiffness, 𝑐(𝑇 ) is a temperature-dependent hardening parameter defined as:

𝑐(𝑇 ) = 𝑔 𝑝 (𝑇 ) 𝑐 0 , ( 3
)
where 𝑐 0 is the hardening parameter given at the reference temperature, while 𝑔 𝑝 (𝑇 ) is another temperature-dependent shift function associated with plasticity. Like 𝑔 𝑒 (𝑇 ), 𝑔 𝑝 (𝑇 ) is equal to 1 at the reference temperature whereas it tends to 0 for high temperatures. Its expression will also be defined later in Section 3. 𝜺 𝑡ℎ (𝑇 ) stands for the thermal strain, usually expressed for an isotropic material as:

𝜺 𝑡ℎ (𝑇 ) = 𝛼𝑰(𝑇 -𝑇 0 ), (4)
where 𝛼 denotes the coefficient of thermal expansion. The last term of (1), namely: 𝜓 𝑡ℎ (𝑇 ), represents the calorific energy, which depends on the temperature. Its expression is given by:

𝜓 𝑡ℎ (𝑇 ) = 𝑐 𝑝 [ (𝑇 -𝑇 0 ) -𝑇 ln ( 𝑇 𝑇 0 )] -𝑠 0 𝑇 + 𝑒 0 , ( 5
)
where 𝑠 0 and 𝑒 0 are the entropy and the internal energy at the reference temperature 𝑇 0 , whereas 𝑐 𝑝 denotes the specific heat capacity at constant pressure. By following the standard thermodynamic derivation, it yields for the stress:

𝝈 = 𝜌 𝜕𝜓 𝜕𝜺 = C(𝑇 ) ∶ ( 𝜺 -𝜺 𝑝 -𝜺 𝑡ℎ (𝑇 ) ) , -𝝈 = 𝜌 𝜕𝜓 𝜕𝜺 𝑝 , ( 6
)
and for the back stress:

𝑿 = 𝜌 𝜕𝜓 𝜕𝝃 𝑝 = 𝑐(𝑇 )(𝝃 𝑝 -𝝃 𝑟 ), -𝑿 = 𝜌 𝜕𝜓 𝜕𝝃 𝑟 . ( 7
)
Furthermore, the specific entropy per mass unit 𝑠 is associated to the temperature such that:

𝑠 = - 𝜕𝜓 𝜕𝑇 . ( 8
)
The state and associated variables are summarized in Table 1. Note that, unlike classical kinematic hardening formulations, the back stress 𝑿 does not only depend on 𝝃 𝑝 , but rather on the difference between the 𝝃 𝑝 and 𝝃 𝑟 , as given in Eq. ( 7). Thus, through its evolution, the recovery variable 𝝃 𝑟 will play an antagonistic role with regard to the back strain 𝝃 𝑝 , by gradually cancelling the effect of the hardening when 𝝃 𝑟 gets close to 𝝃 𝑝 . This process will be further detailed later in Section 2.2.2.

Evolution laws

To guaranty the thermodynamic admissibility, it is important that the evolution laws, expressing the kinetic of the internal state variables, are well in line with the Clausius-Duhem inequality, stating that the intrinsic dissipation must be positive or null. For the proposed formulation, this gives:

Φ = 𝝈 ∶ ̇𝜺 -𝜌( ψ + 𝑠 Ṫ ) ≥ 0 = 𝝈 ∶ ̇𝜺 -𝜌 ( 𝜕𝜓 𝜕𝜺 ∶ ̇𝜺 + 𝜕𝜓 𝜕𝑇 Ṫ + 𝜕𝜓 𝜕𝜺 𝑝 ∶ ̇𝜺𝑝 + 𝜕𝜓 𝜕𝝃 𝑝 ∶ ξ𝑝 + 𝜕𝜓 𝜕𝝃 𝑟 ∶ ξ𝑟 + 𝑠 Ṫ ) ≥ 0 = 𝝈 ∶ ̇𝜺𝑝 -𝑿 ∶ ξ𝑝 + 𝑿 ∶ ξ𝑟 ≥ 0. (9)
In the above equation, two contributions can be identified: the dissipation associated with viscoplasticity, i.e., 𝝈 ∶ ̇𝜺𝑝 -𝑿 ∶ ξ𝑝 , and the one associated with the hardening recovery, i.e., 𝑿 ∶ ξ𝑟 . Since, the viscoplasticity and the hardening recovery mechanism may occur independently from one another, each of these two contributions must be positive or null. For these reasons, two potential functions will be used to derive the evolution laws, one for the viscoplasticity, and another for the hardening recovery.

Viscoplasticity with kinematic hardening

Non-linear kinematic hardening [START_REF] Armstrong | A mathematical representation of the multiaxial bauschinger effect[END_REF]) is commonly employed to represent the asymmetrical behaviour of metallic material, i.e., the Bauschinger effect, upon alternate tension/compression loading stages. Such type of hardening necessitates to be formulated in the context of non-associative plasticity. Accordingly, the evolution laws derive from the normality rule of a convex indicative potential function defined as:

𝐹 (𝝈, 𝑿; 𝑇 ) = 𝑓 (𝝈, 𝑿, 𝑇 ) + 𝛾 2 𝑐(𝑇 ) 𝑿 ∶ 𝑿 ( 10
)
where 𝛾 is an additional hardening parameter, which is assumed to be temperature-independent, while 𝑓 (𝝈, 𝑿, 𝑇 ) represents the yield function given by:

𝑓 (𝝈, 𝑿, 𝑇 ) = eq(𝝈 -𝑿) -𝜎 𝑦 (𝑇 ). (11)
One can notice that the yield function depends on the temperature through the yield threshold 𝜎 𝑦 (𝑇 ), which is defined as:

𝜎 𝑦 (𝑇 ) = 𝑔 𝑝 (𝑇 ) 𝜎 𝑦 0 . ( 12
)
Like the expression of the hardening parameter 𝑐(𝑇 ) in ( 3), the temperature dependency of 𝜎 𝑦 (𝑇 ) is expressed by multiplying the yield threshold at the reference temperature 𝜎 𝑦 0 by the same shift function 𝑔 𝑝 (𝑇 ).

Following the normality rule of 𝐹 (𝝈, 𝑿; 𝑇 ), it yields for the evolution laws:

̇𝜺𝑝 = 𝜕𝐹 𝜕𝝈 ṗ = 𝜦(𝝈, 𝑿) ṗ, ξ𝑝 = - 𝜕𝐹 𝜕𝑿 ṗ = 𝜴(𝝈, 𝑿, 𝑇 ) ṗ, ( 13
)
where 𝑝 is a multiplier ( ṗ ≥ 0), whereas 𝜦(𝝈, 𝑿) and 𝜴(𝝈, 𝑿, 𝑇 ) respectively denote the plastic strain and back strain flows:

𝜦(𝝈, 𝑿) = 3 2 𝐃𝐞𝐯(𝝈 -𝑿) eq(𝝈 -𝑿) , 𝜴(𝝈, 𝑿, 𝑇 ) = 𝜦(𝝈, 𝑿) - 𝛾 𝑐(𝑇 ) 𝑿. ( 14
)
Note that the convexity of the indicative function 𝐹 with regard to the associated variables guarantees that the plasticity-related dissipation is always positive or null, i.e., 𝝈 ∶ ̇𝜺𝑝 -𝑿 ∶ ξ𝑝 ≥ 0.

To account for a viscoplastic behaviour, the evolution of the multiplier 𝑝 is connected to the yield function 𝑓 (𝝈, 𝑿, 𝑇 ) through a Nortonlike law [START_REF] Norton | The Creep of Steel at High Temperatures[END_REF]:

ṗ = 1 𝜂 [ ⟨ 𝑓 (𝝈, 𝑿, 𝑇 ) ⟩ + 𝜎 𝑦 (𝑇 )
] 𝑁 [START_REF] Humphreys | A unified theory of recovery, recrystallization and grain growth, based on the stability and growth of cellular microstructures-I. The basic model[END_REF] where 𝜂, 𝑁 are viscosity parameters. It is worth noticing that in the above equation, the yield function is normalized by the yield threshold 𝜎 𝑦 (𝑇 ), which itself depends on the temperature through the shift function 𝑔 𝑝 (𝑇 ), as defined in Eq. ( 12). In this manner, 𝑔 𝑝 (𝑇 ) acts as a unique temperature-dependent function for the whole viscoplasticity mechanism.

Thermally-activated hardening recovery

The recovery kinetic is mainly governed by the temperature 𝑇 and the dislocation density [START_REF] Petkovic | Recovery and recrystallization of polycrystalline copper after hot working[END_REF][START_REF] Farzadi | Modeling of isothermal recovery and recrystallization kinetics by means of hardness measurements[END_REF][START_REF] Mouelle | Thermally-activated hardening recovery of thermo-elasto-plastic metals during annealing: constitutive modeling for the simulation of welding process[END_REF]. From a phenomenological standpoint, the hardening, through the back stress 𝑿, can be considered as a representative quantity of the actual dislocation density. For these reasons, 𝑇 and 𝑿 are considered as driving forces of the hardening recovery process. To this end, the recovery kinetic is based on a dual dissipation potential defined as a convex function of the back stress and in which the temperature acts as parameter. This gives:

𝜑 * 𝑟 (𝑿; 𝑇 ) = 𝑔 𝑟 (𝑇 ) 𝐴 𝑋 2 𝑿 ∶ 𝑿, ( 16
)
where 𝐴 𝑋 is an amplitude parameter and 𝑔 𝑟 (𝑇 ) is an activation function of the recovery mechanism. The latter is formulated using a Johnson-Mehl-Avrami-Kolmogorov-like law (JMAK) [START_REF] Avrami | Kinetics of phase change. I General theory[END_REF][START_REF] Avrami | Kinetics of phase change. II Transformation-time relations for random distribution of nuclei[END_REF][START_REF] Fanfoni | The Johnson-Mehl-Avrami-Kolmogorov model: A brief review[END_REF]:

𝑔 𝑟 (𝑇 ) = exp ( - ( 𝑇 𝑟 𝑇 ) 𝑛 𝑟 ) . ( 17
)
This form is motivated by the fact that the hardening covers not only the annihilation of dislocations but also the grain recrystallization whose effects are well predicted by the JMAK laws. Fig. 2 shows that, in 𝑔 𝑟 (𝑇 ), 𝑇 𝑟 represents a transition temperature from which the recovery process becomes predominant, whereas 𝑛 𝑟 is a sensitivity parameter regarding this transition.

The evolution law of the recovery variable is then achieved after derivation of this potential with respect to its associated variable, i.e., -𝑿. accordingly, it yields:

ξ𝑟 = 𝜕𝜑 * 𝑟 𝜕𝑿 = 𝑔 𝑟 (𝑇 )𝐴 𝑋 𝑿. ( 18
)
With this formalism, the hardening recovery variable rate is nearly null when the temperature is far below the transition temperature 𝑇 𝑟 , as 𝑔 𝑟 (𝑇 ≪ 𝑇 𝑟 ) → 0, making the recovery mechanism quasi-inactive at low temperatures. However, the recovery mechanism starts becoming significantly active when the temperature approaches 𝑇 𝑟 and eventually becomes fully active, when the temperature is far beyond 𝑇 𝑟 , as 𝑔 𝑟 (𝑇 ≫ 𝑇 𝑟 ) → 1. The hardening recovery variable rate is further amplified by the actual value of 𝑿 through the parameter 𝐴 𝑋 . Note that the convexity of 𝜑 * 𝑟 (𝑿; 𝑇 ) necessarily implies that the hardening recovery-related dissipation is positive or null, i.e., 𝑿 ∶ ξ𝑟 ≥ 0, making the proposed constitutive relation well thermodynamically admissible.

Description of the hardening recovery mechanism

During the hardening recovery, 𝝃 𝑟 runs towards 𝝃 𝑝 as a result of a temperature exposure. This process can be emphasized by combining Eqs. ( 7) and ( 18) together such that:

𝝃 𝑟 + 𝜏 𝑟 (𝑇 ) ξ𝑟 = 𝝃 𝑝 . ( 19
)
This equation appears to be a first order differential equation highlighting a temperature-dependent characteristic time 𝜏 𝑟 (𝑇 ) that can be identified as:

𝜏 𝑟 (𝑇 ) = 1 𝑐(𝑇 )𝑔 𝑟 (𝑇 )𝐴 𝑋 , ( 20
)
where 𝑐(𝑇 ) and 𝑔 𝑟 (𝑇 ) are given in Eqs. ( 3) and ( 17), respectively. 𝜏 𝑟 (𝑇 ) provides a useful information about which range of time the hardening recovery is expected to occur, for a given temperature. One can notice that when the recovery mechanism is quasi-inactive at low temperatures, 𝜏 𝑟 tends to very long times whereas it becomes much shorter when the hardening recovery becomes active at higher temperatures.

Let us consider the case of an already yielded material, undergoing hardening recovery without any active yielding, i.e., 𝝃 𝑝 ≠ 𝟎 and ξ𝑝 = 𝟎, as a result of a sudden temperature exposure held constant, i.e., Ṫ = 0, from an instant 𝑡 0 . Under these conditions, Fig. 3(a) then shows that 𝝃 𝑟 asymptotically tends to 𝝃 𝑝 through an exponential growth:

𝝃 𝑟 (𝑡; 𝑇 ) = 𝝃 𝑝 ( 1 -exp ( - 𝑡 -𝑡 0 𝜏 𝑟 (𝑇 ) ) ) . ( 21
)
According to Eq. ( 7), since the back stress 𝑿, for a given temperature, is proportional to the difference between 𝝃 𝑝 and 𝝃 𝑟 ,Fig. 3(b) shows that 𝑿 relaxes to zero by following an exponential decay:

𝑿(𝑡; 𝑇 ) = 𝑿 0 (𝑇 ) × exp ( - 𝑡 -𝑡 0 𝜏 𝑟 (𝑇 ) ) , with 𝑿 0 (𝑇 ) = 𝑐(𝑇 )𝝃 𝑝 , ( 22
)
gradually shifting the yield function, in Eq. ( 11), back to the stress space centre, as illustrated in Fig. 3(c). Note that according to Eqs. ( 21) or ( 22), 3 × 𝜏 𝑟 represents to time needed to recover 95 % of the hardening. It is worth noticing that a similar trend is obtained with the Zener-Wert-Avrami equation that is commonly utilized to predict the thermal relaxation of residual stresses (Totten et al., 2002;[START_REF] Angkurarach | Effects of high-temperature deep rolling on fatigue, work hardening, and residual stress relaxation of martensitic stainless steel AISI 420[END_REF].

Assumption regarding the temperature field equation

The energy balance equation governs the temperature field within the material domain. By combining together the first law of thermodynamics in its local form with the Fourier law, one obtains for an isotropic medium:

𝑟 + 𝜔 + 𝑘𝛁 2 𝒙 𝑇 = 0, (23)
where 𝜔 stands for an eventual applied heat source, 𝑘 is the material conductivity and 𝛁 2 𝒙 the Laplace operator. Furthermore, 𝑟 represents the differential energy density rate between the mechanical power and the rate of internal energy. Knowing that the latter is connected to the Helmholtz free energy by 𝑒 = 𝜓 + 𝑇 𝑠, it yields for 𝑟:

𝑟 = 𝝈 ∶ ̇𝜺 -𝜌 ė with ė = ψ + 𝑠 Ṫ + 𝑇 ṡ = 𝝈 ∶ ̇𝜺 -𝜌( ψ + 𝑠 Ṫ ) -𝜌𝑇 ṡ, = Φ -𝜌𝑇 ṡ ( 24
)
where the intrinsic dissipation Φ is given in ( 9) for the proposed model formulation.

According to ( 23) and ( 24), It is clear that the deformation as well as the hardening recovery processes obviously play a role in the energy balance through the energy dissipated by the material and the variation of entropy, resulting in temperature changes within the material. This makes the mechanical and thermal problems interdependent such that they should be solved through a fully-coupled thermomechanical analysis, where the equilibrium equations (mechanical problem) and the energy balance equation (thermal problem) are simultaneously solved, while integrating all the coupling terms as detailed in [START_REF] Chatzigeorgiou | Periodic homogenization for fully coupled thermomechanical modeling of dissipative generalized standard materials[END_REF][START_REF] Chatzigeorgiou | Thermomechanical Behavior of Dissipative Composite Materials[END_REF]. However, since the temperature sensitivity of metallic materials is usually quite low compared to the temperature elevations that can be caused by the intrinsic dissipation, its effect can be neglected in the expression of 𝑟, i.e., Φ ≈ 0. Moreover, it can also be assumed that most of the entropy results from the differentiation of the caloric energy 𝜓 𝑡ℎ , given in (5), with respect to the temperature. Therefore, since 𝜓 𝑡ℎ only depends on the temperature, it yields:

𝑠 = - 𝜕𝜓 𝜕𝑇 ≈ - 𝜕𝜓 𝑡ℎ 𝜕𝑇 ⇒ ṡ ≈ - 𝜕 2 𝜓 𝑡ℎ 𝜕𝑇 2 Ṫ = 𝑐 𝑝 Ṫ 𝑇 . ( 25
)
With these two assumptions in mind, the differential energy density rate reduces to:

𝑟 ≈ -𝜌𝑐 𝑝 Ṫ , ( 26
)
which renders the thermal problem independent from the mechanical one. In this context, the whole problem can be conveniently treated by means of an uncoupled thermomechanical analysis, where the thermal problem is first solved prior to the mechanical one. The temperature field computed from the thermal analysis is then utilized as input data for the subsequent mechanical analysis. In that sense, the temperature can be regarded as a controllable quantity, as considered in the remainder of this paper.

Numerical implementation

It is important to remind that a FE solver, like ABAQUS/Standard, employs a backward Euler (or time-implicit) integration scheme. Accordingly, a constitutive law is implemented through an integrator, e.g., a UMAT (User MATerial subroutine) in ABAQUS/Standard, that computes:

(i) The stress as well as history-dependent quantities (internal state variables) over a time increment [𝑡 (𝑛) ; 𝑡 (𝑛+1) ] of length 𝛥𝑡 = 𝑡 (𝑛+1) -𝑡 (𝑛) , while all variables known at 𝑡 (𝑛) . (ii) The tangent operators d𝝈 d𝜺 and d𝝈 d𝑇 that are necessary for the FE solver to compute the global predictor at the whole FE system's level. This is usually done through a return mapping algorithm [START_REF] Simo | Computational Inelasticity[END_REF]Praud et al., 2017a,b;[START_REF] Praud | Multi-Scale Modelling of Thermoplastic-Based Woven Composites[END_REF][START_REF] Chatzigeorgiou | Thermomechanical Behavior of Dissipative Composite Materials[END_REF], which is detailed in Appendix for the proposed model.

Parameters identification for AISI 316L austenitic stainless steel

Flow stress

The proposed model is identified for the AISI 316L austenitic stainless steel following a specific identification procedure. This is done from monotonic tensile tests taken from literature [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF][START_REF] Hendili | HSNV140 -thermo-plasticité avec restauration d'écrouissage : essai de dilatométrie bloquée[END_REF]. These experiments have been performed at six temperatures ranging from 293 K to 1273 K, and for two strain rates, namely: 2.5 × 10 -4 s -1 and 2.5 × 10 -3 s -1 . For each test, the axial plastic strain 𝜀 𝑝 is extracted from the experimentally recorded axial stress 𝜎 and total strain 𝜀:

𝜀 𝑝 = 𝜀 - 𝜎 𝐸(𝑇 ) , where 𝐸(𝑇 ) = 𝑔 𝑒 (𝑇 )𝐸 0 , ( 27
)
while the plastic strain rate is assumed to be constant and equal to the strain rate of the test, such that ε𝑝 ≈ ε. Furthermore, under monotonic uni-axial tension, the plastic strain 𝜀 𝑝 equals the multiplier 𝑝.

The identification procedure then relies on the evaluation of the flow stress 𝜎 as a function of 𝜀 𝑝 , where ε𝑝 and 𝑇 can be regarded as parameters. Based on Eqs. ( 11) and ( 15), it yields under monotonic uni-axial conditions, the following relationship between 𝜎, 𝑋 and ε𝑝 :

𝜎(𝑋; ε𝑝 , 𝑇 ) = 3 2 𝑋 + 𝜎 𝑦 (𝑇 ) [ 1 + (𝜂 ε𝑝 ) 1∕𝑁 ] , ( 28
)
whereas the axial back stress 𝑋 is expressed through Eqs. ( 7), ( 13) and ( 18) and is accordingly governed by the following differential equation:

Ẋ + [ 𝛾 ε𝑝 + 𝑐(𝑇 ) 𝑔 𝑟 (𝑇 ) 𝐴 𝑋 - 𝑔 ′ 𝑝 𝑔 𝑝 Ṫ ] 𝑋 = 𝑐(𝑇 ) ε𝑝 , ( 29
)
which can be solved under monotonic conditions at constant temperature ( Ṫ = 0). This leads to the following analytical expression:

𝑋(𝜀 𝑝 ; ε𝑝 , 𝑇 ) = 𝑐(𝑇 ) ε𝑝 𝛾 ε𝑝 + 𝑐(𝑇 ) 𝑔 𝑟 (𝑇 ) 𝐴 𝑋 [ 1 -exp ( -𝛾𝜀 𝑝 -𝑐(𝑇 ) 𝑔 𝑟 (𝑇 ) 𝐴 𝑋 𝜀 𝑝 ε𝑝 ) ] . ( 30
)
Therefore, by inserting ( 30) into (28), it yields the following expression for the flow stress:

𝜎(𝜀 𝑝 ; ε𝑝 , 𝑇 ) = 3 2 𝑐(𝑇 ) ε𝑝 𝛾 ε𝑝 + 𝑐(𝑇 ) 𝑔 𝑟 (𝑇 )𝐴 𝑋 [ 1 -exp ( -𝛾𝜀 𝑝 -𝑐(𝑇 ) 𝑔 𝑟 (𝑇 ) 𝐴 𝑋 𝜀 𝑝 ε𝑝 ) ] + 𝜎 𝑦 (𝑇 ) [ 1 + (𝜂 ε𝑝 ) 1∕𝑁 ] ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 𝜎 * 𝑦 (𝑇 , ε𝑝 ) , ( 31
)
where it is recalled that 𝑐(𝑇 ) = 𝑔 𝑝 (𝑇 )𝑐 0 and 𝜎 𝑦 (𝑇 ) = 𝑔 𝑝 (𝑇 )𝜎 𝑦 0 . In (31), it is important to point out that the last term, namely:

𝜎 * 𝑦 (𝑇 , ε𝑝 ) = 𝜎 𝑦 (𝑇 ) [ 1 + (𝜂 ε𝑝 ) 1∕𝑁 ] = 𝑔 𝑝 (𝑇 )𝜎 𝑦 0 [ 1 + (𝜂 ε𝑝 ) 1∕𝑁 ] , (32)
can be interpreted as an apparent yield stress, which can be conveniently read on a tensile curve for a given temperature 𝑇 and plastic strain rate ε𝑝 . It is furthermore worth noticing that if the hardening recovery mechanism is not active, which is likely to occur at low temperature, then Eq. ( 31) is reduced to:

𝜎(𝜀 𝑝 ; ε𝑝 , 𝑇 ) = 𝑔 𝑝 (𝑇 ) [ 3 2 𝑐 0 𝛾 [ 1 -exp ( -𝛾𝜀 𝑝 ) ] + 𝜎 𝑦 0 [ 1 + (𝜂 ε𝑝 ) 1∕𝑁 ] ] . ( 33
)
In this condition, one can remark that the 𝑔 𝑝 (𝑇 ) acts as a unique shift function to the whole apparent hardening at the reference temperature.

With or without hardening recovery, the assessment of the flow stress conveniently enables to identify the parameters of the model. The identification procedure is then performed according to the three following steps:

(i) Identification of the shift functions 𝑔 𝑒 (𝑇 ) and 𝑔 𝑝 (𝑇 ) associated with the temperature-dependency of the elastic and viscoplastic properties, (ii) Identification of the viscoplastic parameters, (iii) Identification of the hardening recovery parameters, which are detailed in the next subsections.

Identification of the shift functions

Prior to identifying the model from the flow stress, it is important to calibrate the shift function 𝑔 𝑒 (𝑇 ) and 𝑔 𝑝 (𝑇 ) associated with the elastic and viscoplastic properties, respectively, through Eqs. ( 2), ( 3) and ( 12). Although temperature dependencies are commonly formulated using Arrhenius-like law [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF], no specific forms have been defined at this stage for these two functions.

In this work, given the experimental data [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF][START_REF] Hendili | HSNV140 -thermo-plasticité avec restauration d'écrouissage : essai de dilatométrie bloquée[END_REF], a form based on the Langevin-like formulation, as proposed by [START_REF] Kittel | Introduction to Solid State Physics[END_REF], has been considered for Fig. 4. Shift function 𝑔 𝑒 (𝑇 ) and 𝑔 𝑝 (𝑇 ) associated with the temperature dependencies of the elastic and viscoplastic properties, respectively. These two functions 𝑔 𝑒 (𝑇 ) are calibrated from experimental data taken from [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF] and [START_REF] Hendili | HSNV140 -thermo-plasticité avec restauration d'écrouissage : essai de dilatométrie bloquée[END_REF].

𝑔 𝑒 (𝑇 ). This function employs a multiplicative decomposition of a second order polynomial and exponential function. Accordingly, it reads:

𝑔 𝑒 (𝑇 ) = ( 𝑎 𝑒 1 + 𝑎 𝑒 2 𝑇 + 𝑎 𝑒 3 𝑇 2 ) × exp ( - ( 𝑇 𝑇 𝑒 ) 𝑛 𝑒 ) , ( 34
)
where the parameters 𝑎 𝑒 1 , 𝑎 𝑒 2 , 𝑎 𝑒 3 , 𝑛 𝑒 , 𝑇 𝑒 are calibrated from the experimental data by calculating for each test the ratio between the Young's modulus 𝐸(𝑇 ) and the one at the reference temperature 𝐸 0 . Fig. 4(a) shows the obtained curve for 𝑔 𝑒 (𝑇 ) with the identified parameters as well as the experimental data from [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF] and [START_REF] Hendili | HSNV140 -thermo-plasticité avec restauration d'écrouissage : essai de dilatométrie bloquée[END_REF].

Similarly to the elastic properties, the viscoplastic properties are very sensitive to the temperature. For 𝑔 𝑝 (𝑇 ), a formulation similar to 𝑔 𝑒 (𝑇 ) has been considered to capture the temperature-dependency of the viscoplastic properties. This function introduces another set of material parameters 𝑎 𝑝 1 , 𝑎 𝑝 2 , 𝑎 𝑝 3 , 𝑛 𝑝 , 𝑇 𝑝 , such that:

𝑔 𝑝 (𝑇 ) = ( 𝑎 𝑝 1 + 𝑎 𝑝 2 𝑇 + 𝑎 𝑝 3 𝑇 2 ) × exp ( - ( 𝑇 𝑇 𝑝 ) 𝑛 𝑝 ) , ( 35
)
where the calibration is done by calculating for each test the ratio between the apparent yield stress 𝜎 * 𝑦 (𝑇 , ε𝑝 ) and the one at room temperature considering the same strain rate, 𝜎 * 𝑦 (𝑇 0 , ε𝑝 ), according to Eq. ( 32). The identified function 𝑔 𝑝 (𝑇 ) agrees well with the experimental data from [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF] and [START_REF] Hendili | HSNV140 -thermo-plasticité avec restauration d'écrouissage : essai de dilatométrie bloquée[END_REF] as shown in Fig. 4(b). As expected, the apparent yield stress is much more sensitive at low temperatures than the Young's modulus. It can be noticed that 𝑔 𝑝 (𝑇 ) exhibits three regimes. At low temperatures (below 673 K), the apparent yield stress decreases about 40 % compared to the reference temperature. Then, at intermediate temperatures between 673 K and 1073 K, the apparent yield stress is halved compared to the reference temperature but it decreases less quickly than for the first regime. Finally, for temperatures above 1073 K, the apparent yield stress drastically drops and nearly vanishes reaching a value close to 0 for temperatures above 1473 K. Although the experimental data are few for the apparent yield stress, it is important to point out that the identified temperature dependency is well consistent with known tendencies for other metallic materials of the same type [START_REF] Gardner | Fire testing and design of stainless steel structures[END_REF][START_REF] Gardner | Elevated temperature material properties of stainless steel alloys[END_REF].

Identification of the viscoplastic parameters

At low temperatures, typically below 873 K, the recovery mechanism can be considered as inactive [START_REF] Kerrouault | Fissuration à chaud en soudage d'un acier inoxydable austénitique[END_REF][START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF]. Under this condition, as explained in Section 3.1, the expression of the flow stress (33) only involves the shift function 𝑔 𝑝 (𝑇 ) and the parameters related to viscoplasticity, namely: 𝜎 𝑦 0 𝑐 0 , 𝛾, 𝜂 and 𝑁. Since 𝑔 𝑝 (𝑇 ) is known at this stage, one can identify these parameters by minimizing the least squares between the experimentally recorded stresses at 𝑇 < 873 K and the ones calculated using Eq. (33).

Identification of the hardening recovery parameters

At high temperatures, typically above 873 K, the hardening recovery mechanism can be considered as active [START_REF] Kerrouault | Fissuration à chaud en soudage d'un acier inoxydable austénitique[END_REF][START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF]. In this condition, the expression of the flow stress (31) now involves the parameters related to viscoplasticity as well as the one related to the hardening recovery, namely: 𝐴 𝑋 , 𝑛 𝑟 and 𝑇 𝑟 . However, since the shift function 𝑔 𝑝 (𝑇 ) and the viscoplasticity parameters are known at this stage, one can identify the hardening recovery parameters by minimizing the least squares between the experimentally recorded stresses at 𝑇 > 873 K and the ones calculated using Eq. ( 31). The obtained parameters are listed in Table 2 while Figs. 5(a) and 5(b) show good agreements between the stress-strain responses computed with these parameters and the experimental data over the whole temperature range and for the two considered strain-rates. Fig. 5. Comparison between the stress-strain computed responses and the experimental data over the whole temperature range and for the two considered strain-rates, namely: ε = 2.5 × 10 -4 s -1 and ε = 2.5 × 10 -3 s -1 . The experimental data, taken from [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF] and [START_REF] Hendili | HSNV140 -thermo-plasticité avec restauration d'écrouissage : essai de dilatométrie bloquée[END_REF], are represented by dot markers whereas the computed responses are depicted by solid lines.

Fig. 6(a) displays the identified activation function 𝑔 𝑟 (𝑇 ), which is part of the evolution laws of the hardening recovery mechanism through the parameters 𝑇 𝑟 and 𝑛 𝑟 (see Sub-Section 2.2.2). It shows that the activation of the recovery mechanism is about 973 K, which is well in line with previous works dealing with this material [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF][START_REF] Hendili | HSNV140 -thermo-plasticité avec restauration d'écrouissage : essai de dilatométrie bloquée[END_REF][START_REF] Mouelle | Thermally-activated hardening recovery of thermo-elasto-plastic metals during annealing: constitutive modeling for the simulation of welding process[END_REF]. Fig. 6(b) also shows for reference the characteristic time 𝜏 𝑟 (𝑇 ) defined in Eq. ( 20) as a function of the temperature. Below the activation temperature of the recovery mechanism, about 973 K, it appears that 𝜏 𝑟 is greater than 10 6 s, which represents about 11 days. Such amount of time is way beyond the time-scale on which thermomechanical loads are usually applied so that the recovery mechanism can be considered as quasi-inactive. The characteristic time 𝜏 𝑟 then decreases for higher levels of temperature. Around 1123 K, one can notice that 𝜏 𝑟 goes below a second. Such a drastic variation in the recovery kinetics has already been reported for other type of austenitic stainless steels [START_REF] Taylor | The post-deformation recrystallization behaviour of 304 stainless steel following high strain rate deformation[END_REF].

Note that, the proposed identification procedure only requires monotonic tension tests, which is sufficient as long as the sole presence of kinematic hardening is assumed. However, it is worth mentioning that, if a combined isotropic-kinematic hardening was considered, a much more complex procedure based on alternate tension-compression tests would certainly be necessary for a proper identification.

Validation examples

In this section, three examples of uni-axial simulations are presented to validate the model while providing a better understanding regarding its behaviour. In these examples, the material undergoes hardening recovery under complex thermomechanical loading paths which are of particular interest towards potential applications. The simulations were carried out with the parameters identified in Section 3 for AISI 316L austenitic stainless steel.

Double hit tension tests

In this example, the so-called double hit tension test is applied to the material. Fig. 7 shows that a stress of 300 MPa is first applied in 100 s before being released to zero in the next 100 s, during this first loading stage, the material is kept at room temperature (reference temperature). In a second stage, the material is kept free of stress for a certain period of time during which the temperature is raised up to temperature 𝑇 max at the heating rate of 10 K s -1 . Once reached this temperature is held constant for a time 𝑡 hold after which the material is cooled down to room temperature at a cooling rate of -10 K s -1 . Now comes the third and last stage in which a stress of 350 MPa is now applied in 100 s and next released during the following 100 s, while the material remains at room temperature.

The double hit tension test was first proposed by [START_REF] Mouelle | Thermally-activated hardening recovery of thermo-elasto-plastic metals during annealing: constitutive modeling for the simulation of welding process[END_REF] to dissociate the hardening recovery from viscoplasticity as these mechanisms are separately activated during such a test. Indeed, during the first stage, the material yields without any recovery as the temperature is too low. During the second stage, the material is kept at zero stress so that no more plastic strain and hardening can be generated. However, hardening recovery occurs as the temperature increases beyond the activation temperature. During this stage, the amount of recovered hardening depends on the reached temperature 𝑇 max as well as the time of temperature exposure 𝑡 hold . The effect of the hardening recovery can then be observed during the third stage, from the level of stress from which the material starts yielding again. If yielding occurs from the stress reached at the first stage, then there is no hardening recovery during the second stage. Contrariwise, if yielding occurs at a lower stress level, then there is hardening recovery, knowing that the initial yielding point can be restored in case of full recovery.

It is recalled that the characteristic time 𝜏 𝑟 (𝑇 ), given in Eq. ( 20) and displayed in Fig. 6(b) for the identified set of parameters, is directly related to time that is necessary to fully recover the hardening, i.e., 95 % of the hardening is recovered in 3 × 𝜏 𝑟 .

Fig. 8 illustrates a case of a partial recovery. In this simulation, the material is heated slightly above the activation temperature, at 𝑇 max = 1023 K, which is held for 𝑡 hold = 100 s. Fig. 8(b) shows that, during the second stage, the hardening was barely recovered as 𝜉 𝑟 slightly increased towards 𝜉 𝑝 . Indeed, in this case, the thermal exposure was too short and too low to achieve full recovery. From an external standpoint, this can be observed at the third stage as the material starts yielding at a stress level nearly equal to the one reached during the first stage, as shown in Fig. 8(a). One can notice that since 𝜏 𝑟 (𝑇 max = 1023 K) = 583 s is much greater than 𝑡 hold = 100 s. This is well consistent since the hardening was barely recovered for this simulation.

Fig. 9 now illustrates a case of full recovery by increasing the time of temperature exposure to 𝑡 hold = 3600 s. Fig. 9(b) shows that, in this simulation, the hardening is fully recovered as 𝜉 𝑟 had the time to asymptotically reach 𝜉 𝑝 during the second stage. Therefore during the third stage, the material starts yielding from its initial yielding point, as shown in Fig. 9(a). In this simulation, since 𝜏 𝑟 (𝑇 max = 1023 K) = 583 s is now much shorter than 𝑡 hold = 3600 s, a full recovery was well expected.

Fig. 10 illustrates another case of full recovery. In this simulation, the material is heated at a higher temperature, namely: 𝑇 max = 1073 K for 𝑡 hold = 100 s. In spite of such a short time of exposure, Fig. 10( shows that the higher temperature level results in a much faster recovery, which enables the material to get through to a full recovery and an entire restoration of its yielding point, as depicted in Fig. 10(a). In this simulation, since 𝜏 𝑟 (𝑇 max = 1073 K) = 11 s is now much shorter than 𝑡 hold = 100 s, a full recovery was well expected too.

Stress relaxation tests

In this example, it is proposed to reproduce with the model the stress relaxation tests carried out by [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF] on AISI 316L austenitic stainless steel. As depicted in Fig. 11, in these experiments, the material is first heated from room temperature up to certain temperature 𝑇 max under zero stress conditions so that it is free to expand. This temperature is then held constant over the whole remainder of the test. Once at temperature, a strain 𝜀 𝑎 is applied at the rate of 2.5 × 10 -3 s -1 before being held constant for a time 𝑡 hold . Once this time has elapsed, another 𝜀 𝑎 is applied to the material at the same rate.

According the data of [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF], two test configurations are considered. In the first one, the material is heated up to 𝑇 max = 473 K, while an applied strain of 𝜀 𝑎 = 0.01 is considered, for an holding time 𝑡 hold = 500 s. In the second configuration, the material is heated up to 𝑇 max = 1073 K with an applied strain of 𝜀 𝑎 = 0.015, for an holding time 𝑡 hold = 1000 s.

To further emphasize the role played by the hardening recovery, both configurations are simulated with the proposed model with and without hardening recovery. In this latter case, the hardening recovery is turned off by virtually setting the parameter 𝐴 𝑋 to zero. Fig. 12 shows the results for the first configuration, for which the temperature remains far below the activation temperature of the hardening recovery mechanism, about 973 K. Therefore, when the Fig. 12. Stress relaxation test for 𝑇 max = 473 K, 𝜀 𝑎 = 0.01 and 𝑡 hold = 500 s. Note that, in this case, the applied temperature, i.e., 473 K, is far below the activation temperature of the hardening recovery that is about 973 K.

strain is held constant, the material exhibits stress relaxation under the sole effect of viscoplasticity. Fig. 12(b) then shows that the recovery variable does not evolve in both simulations with and without recovery, resulting in identical stress-strain responses in Fig. 12(a), which are in relative good agreement with the experimental data.

Fig. 13 now shows the results of the second configuration for which the temperature goes beyond the activation temperature. Thus, the material exhibits stress relaxation resulting from the combined effects of viscoplasticity and the hardening recovery, promoting the development of plastic strain and therefore increasing the amount of relaxed stress. Fig. 13(b) shows that, during the relaxation, the recovery variable 𝜉 𝑟 runs towards the back strain 𝜉 𝑝 while assisting its development. This proceeds up to a point where both variables eventually take the same value. Comparatively with the response without recovery, it can be noticed in Fig. 13(a) that the response with recovery is in much better agreement with the experimental data.

Satoh test

For this last example, the Satoh test carried out by [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : Validation sur cas tests de complexité croissante[END_REF], [START_REF] Depradeux | Influence of viscoplasticity, hardening, and annealing effects during the welding of a three-pass slot weld (NET-TG4 round robin)[END_REF] on AISI 316L austenitic stainless steel is reproduced with the model. The Satoh test is an experiment that is typically representative of the multi-pass welding conditions [START_REF] Satoh | Thermal stresses developed in high-strength steels subjected to thermal cycles simulating weld heat-affected zone[END_REF][START_REF] Satoh | Transient thermal stresses ofweld heat-affected zone by both-ends fixed bar analogy[END_REF][START_REF] Mouelle | Caractérisation et modélisation d'assemblages de tôles en acier inoxydable austénitique brasées et soudées pour des applications d'échangeurs de chaleur[END_REF]. It consists in applying thermal cycles to the material while its axial strain is continuously to zero. In the present case, Fig. 14 shows that four successive cycles are considered. The material is heated at the rate of 90 K s -1 from room temperature (𝑇 0 = 293 K) up to 1348 K, 1128 K, 908 K and 685 K, for cycle 1, 2, 3 and 4, respectively, and subsequently cooled down to room temperature (𝑇 0 = 293 K) at the rate of 0.3 K s -1 . As for the Stress relaxation test, the Satoh test is simulated with the proposed model with and without hardening recovery. It is recalled that, in this latter case, the hardening recovery is turned off by virtually setting the parameter 𝐴 𝑋 to zero. Fig. 15(a) shows the comparison between the experimental data for the Satoh test and the simulated stress-temperature response without hardening recovery. During each thermal cycle, as it is prevented to expand, the material undergoes alternate compression and tension states with yielding over the course of the heating and cooling stages, respectively. According to the model, Fig. 15(b) shows that the back strain variable 𝜉 𝑝 makes the yield function shift in compression and in tension during the heating and cooling stages, respectively. Note that the back strain recovery does not evolve here since the hardening recovery is turned off in this simulation. In spite of that, it can be noticed that the model predicts a similar tendency compared to the experiment. During the heating stage of the first thermal cycle (from 293 K to 1348 K), the stress is overestimated of about 50 MPa. Then, during the next cooling stage, (from 1348 K to 293 K), the stress is quite close from the experimental one. For the next cycles, the heating stages are well predicted. However, during the cooling stages, the yielding point occurs sooner than what is experimentally observed. Fig. 16(a) now shows the comparison between the experimental data and the simulated stress-temperature response with hardening recovery. In this case, the contribution of the hardening recovery allows to improve the prediction at the end of the first heating stage, while the yielding point is slightly shifted at the beginning of the next cooling stage. This greatly improves the prediction during the second cycle whereas the following cycles keep nearly unchanged compared to the simulation without hardening recovery. As depicted in Fig. 16(b), the recovery variable 𝜉 𝑟 evolves towards 𝜉 𝑝 only during the two first cycles for which the temperature exceeds the activation temperature, i.e., about 973 K. During this time, the recovery mechanism tends to bring the yield function back to its state without hardening, counteracting the effect of the latter as it progresses.

Conclusions and perspectives

In this work, a new phenomenological model has been proposed to describe the thermally-activated hardening recovery mechanisms in metallic materials. Such a feature constitutes an original modelling approach to capture the effects caused by different phenomena occurring at the microstructure level, such as dislocation annihilations, recrystallization and grain growth. Comparatively to previous works, where the hardening recovery was only considered in a classical plasticity formalism with isotropic hardening [START_REF] Mouelle | Thermally-activated hardening recovery of thermo-elasto-plastic metals during annealing: constitutive modeling for the simulation of welding process[END_REF], the proposed model rather extends this concept in a viscoplastic formulation with kinematic hardening. Besides, the model also introduces temperaturedependent shift functions to better represent the thermomechanical response of the material over a wide temperature range whether the hardening recovery is active or not. The model's parameters have been successfully identified for AISI 316L austenitic stainless steel through a stepwise identification procedure based on the evaluation of the flow stress, while validation is further achieved on other representative tests, namely double hit tension tests, stress relaxation tests and Satoh tests for which the material undergoes hardening recovery under complex thermomechanical loading paths.

As future prospects, it is intended to extend the present formulation to combined isotropic-kinematic hardening. Although this might require additional thermomechanical experiments (Thiercelin et al., 2022b) as well as a revised identification procedure, a combined hardening might be helpful to better capture the material's response upon alternate tension-compression like it is the case for the Satoh test. Extensions to strain gradient plasticity (Voyiadjis and Song, 2019;[START_REF] Martínez-Pañeda | The role of plastic strain gradients in the crack growth resistance of metals[END_REF][START_REF] Jebahi | Scalar-based strain gradient plasticity theory to model sizedependent kinematic hardening effects[END_REF] and/or nonlocal damage [START_REF] Kiefer | A gradient-enhanced damage model coupled to plasticity-multi-surface formulation and algorithmic concepts[END_REF][START_REF] Praud | Phase field fracture models to predict crack initiation and propagation in antireflective coatings[END_REF][START_REF] Satouri | A gradient enhanced constitutive framework for the investigation of ductile damage localization within semicrystalline polymers[END_REF] are also under consideration to account for strain localization and other size effects. Finally, establishing correlations between microstructural changes and the evolution of the internal state variables [START_REF] Thiercelin | Thermomechanical modelling of the tribological surface transformations in the railroad network (White eching layer)[END_REF](Thiercelin et al., , 2022a) ) describing the material behaviour at the upper scale would be a valuable information towards enriching the proposed constitutive model.

CRediT authorship contribution statement

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Data will be made available on request.

Appendix. Return mapping algorithm

The numerical implementation of the proposed model is based on the convex cutting plane form of the return mapping algorithm, which is detailed in this appendix. The reader can refer to several references among them [START_REF] Simo | Computational Inelasticity[END_REF]Praud et al., 2017a,b;[START_REF] Praud | Multi-Scale Modelling of Thermoplastic-Based Woven Composites[END_REF][START_REF] Chatzigeorgiou | Thermomechanical Behavior of Dissipative Composite Materials[END_REF].

A.1. Residuals and linearization of the constitutive equations

The evolution laws of the viscoplastic multiplier 𝑝 and the recovery variable 𝝃 𝑟 , given in ( 15) and ( 18), respectively, are written under the form of residuals that must both satisfy a nullity condition. Accordingly, this gives:

𝑅 𝑝 ( ṗ, 𝝈, 𝑿, 𝑇 ) = ṗ - 1 𝜂 [ ⟨ 𝑓 (𝝈, 𝑿, 𝑇 ) ⟩ + 𝜎 𝑦 (𝑇 ) ] 𝑁 , (A.1a) 𝑹 𝝃 𝑟 ( ξ𝑟 , 𝑿, 𝑇 ) = ξ𝑟 -𝑔 𝑟 (𝑇 )𝐴 𝑋 𝑿. (A.1b)
Taken in their linearized forms, it yields:

𝛿𝑅 𝑝 = 𝜕𝑅 𝑝 𝜕 ṗ 𝜕 ṗ 𝜕𝑝 𝛿𝑝 + 𝜕𝑅 𝑝 𝜕𝝈 ∶ 𝛿𝝈 + 𝜕𝑅 𝑝 𝜕𝑿 ∶ 𝛿𝑿 + 𝜕𝑅 𝑝 𝜕𝑇 𝛿𝑇 = 𝐴 𝑝𝑝 𝛿𝑝 + 𝑨 𝑝𝝈 ∶ 𝛿𝝈 + 𝑨 𝑝𝑿 ∶ 𝛿𝑿 + 𝐴 𝑝𝑇 𝛿𝑇 , (A.2a) 𝛿𝑹 𝝃 𝑟 = 𝜕𝑹 𝝃 𝑟 𝜕 ξ𝑟 ∶ 𝜕 ξ𝑟 𝜕𝝃 𝑟 ∶ 𝛿𝝃 𝑟 + 𝜕𝑹 𝝃 𝑟 𝜕𝑿 ∶ 𝛿𝑿 + 𝜕𝑹 𝝃 𝑟 𝜕𝑇 ∶ 𝛿𝑇 = A 𝝃 𝑟 𝝃 𝑟 ∶ 𝛿𝝃 𝑟 + A 𝝃 𝑟 𝑿 ∶ 𝛿𝑿 + 𝑨 𝝃 𝑟 𝑇 𝛿𝑇 . (A.2b)
Let us recall that the convex cutting plane method [START_REF] Simo | Computational Inelasticity[END_REF]Praud et al., 2017a,b;[START_REF] Praud | Multi-Scale Modelling of Thermoplastic-Based Woven Composites[END_REF][START_REF] Chatzigeorgiou | Thermomechanical Behavior of Dissipative Composite Materials[END_REF] considers simplified linearized forms of the flow Eqs. ( 13

. Elastic prediction with hardening recovery correction

To assess if the material is actively yielding or not, an elastic prediction is necessary. During this step, the internal state variables, 𝜺 𝑝 , and 𝝃 𝑝 as well as the viscoplastic multiplier 𝑝 are kept constant such that: 𝜺 (𝑛+1) 𝑝 = 𝜺 (𝑛) 𝑝 , 𝝃 (𝑛+1)

𝑝

= 𝝃 (𝑛) 𝑝 and 𝑝 (𝑛+1) = 𝑝 (𝑛) . However, since the recovery variable 𝝃 𝑟 defines the yielding point while having no influence on the elastic behaviour, it is important to let it evolve during the elastic prediction. Accordingly, after an initialization step in which 𝝃 (𝑛+1)(𝑘=0) 𝑟 = 𝝃 (𝑛) 𝑟 , the recovery variable is computed by means of iterative corrections such that: Therefore, the correction loop is repeated until reaching convergence, namely, when the residual 𝑹 𝝃 𝑟 is sufficiently close to zero. Once converged, the yield function 𝑓 can be checked to assess whether the material is actively yielding or not:

• If 𝑓 < 0, then the material is not yielding and no more corrections are needed. • if 𝑓 > 0, then the material is actively yielding and another corrective step is necessary to account for the viscoplasticity in addition to the hardening recovery.

A.2.2. Viscoplastic and hardening recovery correction

If the material is actively yielding, i.e., 𝑓 > 0, then all the internal state variables, 𝜺 𝑝 , 𝝃 𝑝 and 𝝃 𝑟 as well as the viscoplastic multiplier 𝑝 evolve. Accordingly, after an initialization step in which 𝑝 (𝑛+1)(𝑘=0) = 𝑝 (𝑛) , 𝜺 (𝑛+1)(𝑘=0) 𝑝 = 𝜺 (𝑛) 𝑝 , 𝝃 (𝑛+1)(𝑘=0)

𝑝

= 𝝃 (𝑛) 𝑝 and 𝝃 (𝑛+1)(𝑘=0)

𝑟

= 𝝃 (𝑛) 𝑟 and, all these variables are iteratively computed by: 𝑝 (𝑛+1)(𝑘+1) = 𝑝 (𝑛+1)(𝑘) + 𝛿𝑝 (𝑛+1)(𝑘) , (A.12a)

Fig. 1 .

 1 Fig. 1. Representation of the microstructural mechanisms behind the thermally-activated hardening recovery phenomenon. Namely: dislocations rearrangement and annihilation, grain recrystallization and their subsequent growth. Source: Figure adapated from Humpreys et al. (2017).

Fig. 2 .

 2 Fig. 2. Activation function of the recovery mechanism 𝑔 𝑟 (𝑇 ). This function involves two parameters, namely: a transition temperature 𝑇 𝑟 and 𝑛 𝑟 , a sensitivity parameter regarding this transition.

Fig. 3 .

 3 Fig. 3.Illustrative case of an already yielded material, undergoing hardening recovery without any active yielding, i.e., 𝝃 𝑝 ≠ 𝟎 and ξ𝑝 = 𝟎, as a result of a sudden temperature exposure held constant, i.e., Ṫ = 0, from an instant 𝑡 0 .

Fig. 6 .

 6 Fig. 6. Identified function 𝑔 𝑟 (𝑇 ) and its associated characteristic time 𝜏 𝑟 (𝑇 ).

Fig. 7 .

 7 Fig. 7. Thermomechanical loading path for the double hit tension tests.

Fig. 8 .

 8 Fig. 8. Double hit tension test for 𝑇 max = 1023 K and 𝑡 hold = 100 s.

Fig. 9 .

 9 Fig. 9. Double hit tension test for 𝑇 max = 1023 K and 𝑡 hold = 3600 s.

Fig. 10 .

 10 Fig. 10. Double hit tension test for 𝑇 max = 1073 K and 𝑡 hold = 100 s.

Fig. 11 .

 11 Fig. 11. Thermomechanical loading path for the stress relaxation tests.

Fig. 13 .

 13 Fig.13. Stress relaxation test for 𝑇 max = 1073 K, 𝜀 𝑎 = 0.015 and 𝑡 hold = 1000 s. Note that, in this case, the applied temperature, i.e., 1073 K, is greater than the activation temperature of the hardening recovery that is about 973 K.

Fig. 14 .

 14 Fig. 14. Thermomechanical loading path for the Satoh test.

Fig. 15 .

 15 Fig. 15. Satoh test without hardening recovery.

Fig. 16 .

 16 Fig. 16. Satoh test with hardening recovery.

  𝑘 refers to the correction loop. Then, 𝛿𝝃 𝑟 is obtained from the nullity condition of the residual 𝑹 𝝃 𝑟 in its linearized form (A.6b), while considering 𝛿𝑝 = 0, 𝛿𝜺 𝑝 = 𝟎 and 𝛿𝝃 𝑝 = 𝟎, along with 𝛿𝜺 = 𝟎 and 𝛿𝑇 = 0 (all the quantities are taken at (𝑛+1)(𝑘) ):𝑹 𝝃 𝑟 + 𝛿𝑹 𝝃 𝑟 = 𝟎 ⇔ 𝛿𝝃 𝑟 = -K -1 𝝃 𝑟 𝝃 𝑟 ∶ 𝑹 𝝃 𝑟 . (A.11)

  𝜺 (𝑛+1)(𝑘+1) 𝑝 = 𝜺 (𝑛+1)(𝑘) 𝑝 + 𝛿𝜺 (𝑛+1)(𝑘) 𝑝 , (A.12b) 𝝃 (𝑛+1)(𝑘+1) 𝑝 = 𝝃 (𝑛+1)(𝑘) 𝑝 + 𝛿𝝃 (𝑛+1)(𝑘)

Table 1

 1 State and associated variables. The present formulation involves the total strain 𝜺 and the absolute temperature 𝑇 as observable state variables, to which the stress 𝝈 and the specific entropy 𝑠, respectively, are associated. Furthermore, the plastic strain 𝜺 𝑝 , the back strain 𝝃 𝑝 and the recovery variable 𝝃 𝑟 are involved as internal state variables, to which -𝝈, the back stress 𝑿 and -𝑿, respectively, are associated.

	State variables		Associated variables
	Observable	Internal	-
	𝜺		𝝈
	𝑇		𝑠
		𝜺 𝑝	-𝝈
		𝝃 𝑝	𝑿
		𝝃 𝑟	-𝑿

Table 2

 2 Identified parameters for AISI 316L austenitic stainless steel.

	Feature	Parameter	Value	Unit
	Thermo-elasticity	𝐸 0	195 600	MPa
		𝜈 0	0.3	-
		𝛼	17.1 × 10 -6	K -1
		𝑎 𝑒 1 𝑎 𝑒 2 𝑎 𝑒 3	1.00 6.98 × 10 -5 -3.57 × 10 -7	-K -1 K -2
		𝑛 𝑒	22.4	-
		𝑇 𝑒	1629	K
	Thermo-visco-plasticity	𝑐 0	4.0 × 10 4	MPa
		𝛾	358	-
		𝑁	10	-
		𝜂	1.45 × 10 4	s
		𝜎 𝑦 0 𝑎 𝑝 1 𝑎 𝑝 2 𝑎 𝑝 3	100 1.61 -2.52 × 10 -3 1.54 × 10 -6	MPa -K -1 K -2
		𝑛 𝑝	6.83	-
		𝑇 𝑝	1234	K
	Hardening recovery	𝐴 𝑋	9.76 × 10 -4	MPa -1 s -1
		𝑛 𝑟	11.9	-
		𝑇 𝑟	1234	K

Léo Thiercelin :

 Thiercelin Writing -original draft, Formal analysis, Investigation, Methodology, Validation. Francis Praud: Conceptualization, Writing -original draft, Formal analysis, Investigation. Fodil Meraghni: Methodology, Formal analysis, Review & editing, Conceptualization. Eric Fleury: Methodology, Formal analysis, Review & editing, Conceptualization.

  𝑩 𝑿𝑇 𝛿𝑇 + 𝑩 𝑿𝑝 𝛿𝑝 + B 𝑿𝝃 𝑟 ∶ 𝛿𝝃 𝑟 .Thus, inserting (A.4) and (A.5) into (A.2a) and (A.2b) eventually leads to the full linearization of the residuals:𝛿𝑅 𝑝 = 𝐾 𝑝𝑝 𝛿𝑝 + 𝑲 𝑝𝝃 𝑟 ∶ 𝛿𝝃 𝑟 + 𝑲 𝑝𝜺 ∶ 𝛿𝜺 + 𝐾 𝑝𝑇 𝛿𝑇 , (A.6a) 𝛿𝑹 𝝃 𝑟 = 𝑲 𝝃 𝑟 𝑝 𝛿𝑝 + K 𝝃 𝑟 𝝃 𝑟 ∶ 𝛿𝝃 𝑟 + K 𝝃 𝑟 𝜺 ∶ 𝛿𝜺 + 𝑲 𝝃 𝑟 𝑇 𝛿𝑇 , 𝑲 𝝃 𝑟 𝑝 K 𝝃 𝑟 𝝃 𝑟 𝝃 𝑟 𝝃 𝑟 = A 𝝃 𝑟 𝝃 𝑟 + A 𝝃 𝑟 𝑿 ∶ B 𝑿𝝃 𝑟 , (A.8f) 𝑲 𝝃 𝑟 𝑇 = A 𝝃 𝑟 𝑿 ∶ 𝑩 𝑿𝑇 + 𝑨 𝝃 𝑟 𝑇 ,

	A 𝝃 𝑟 𝑿 =	𝜕𝑹 𝝃 𝑟 𝜕𝑿	= -𝑔 𝑟 (𝑇 )𝐴 𝑋 I,		(A.9f)													), which
	𝑨 𝝃 𝑟 𝑇 =	𝜕𝑹 𝝃 𝑟 𝜕𝑇	= -𝑔 ′ 𝑟 (𝑇 )𝐴 𝑋 𝑿,		(A.9g)	become: 𝛿𝜺 𝑝 = 𝜦(𝝈, 𝑿)𝛿𝑝,			𝛿𝝃 𝑝 = 𝜴(𝝈, 𝑿, 𝑇 )𝛿𝑝,	(A.3)
	B 𝝈𝜺 =	𝜕𝝈 𝜕𝜺	= 𝑔 𝑒 (𝑇 ) C 0 ,			(A.9h)	so that the linearized stress and back stress, according to Eqs. (6) and (7), are given by:
	𝑩 𝝈𝑇 =		𝜕𝝈 𝜕𝑇		= 𝑔 ′ 𝑒 (𝑇 ) C 0 ∶	( 𝜺 -𝜺 𝑝 -𝜺 𝑡ℎ (𝑇 ) )	-𝑔 𝑒 (𝑇 ) C 0 ∶ 𝛼𝑰,	(A.9i)	𝛿𝝈 =	𝜕𝝈 𝜕𝜺	∶ 𝛿𝜺 +	𝜕𝝈 𝜕𝑇	𝛿𝑇 +	𝜕𝝈 𝜕𝜺 𝑝	∶ 𝛿𝜺 𝑝
	𝑩 𝝈𝑝 =	𝜕𝝈 𝜕𝜺 𝑝		∶ 𝜦(𝝈, 𝑿) = -𝑔 𝑒 (𝑇 ) C 0 ∶ 𝜦(𝝈, 𝑿),	(A.9j)		=	𝜕𝝈 𝜕𝜺	∶ 𝛿𝜺 +	𝜕𝝈 𝜕𝑇	𝛿𝑇 +	𝜕𝝈 𝜕𝜺 𝑝	∶ 𝜦(𝝈, 𝑿)𝛿𝑝	(A.4)
	𝑩 𝑿𝑇 = 𝑩 𝑿𝑝 = B 𝑿𝝃 𝑟 =	𝜕𝑿 𝜕𝑇 𝜕𝑿 𝜕𝝃 𝑝 𝜕𝑿 𝜕𝝃 𝑟	= 𝑔 ′ 𝑝 (𝑇 ) 𝑐 0 (𝝃 𝑝 -𝝃 𝑟 ), ∶ 𝜴(𝝈, 𝑿, 𝑇 ) = 𝑔 𝑝 (𝑇 ) 𝑐 0 𝜴(𝝈, 𝑿, 𝑇 ), = -𝑔 𝑝 (𝑇 ) 𝑐 0 I.	(A.9k) (A.9l) (A.9m)	= B 𝝈𝜺 ∶ 𝛿𝜺 + 𝑩 𝝈𝑇 𝛿𝑇 + 𝑩 𝝈𝑝 𝛿𝑝, 𝛿𝑿 = 𝜕𝑿 𝜕𝑇 𝛿𝑇 + 𝜕𝑿 𝜕𝝃 𝑝 ∶ 𝛿𝝃 𝑝 + 𝜕𝑿 𝜕𝝃 𝑟 ∶ 𝛿𝝃 𝑟 = 𝜕𝑿 𝜕𝑇 𝛿𝑇 + 𝜕𝑿 𝜕𝝃 𝑝 ∶ 𝜴(𝝈, 𝑿, 𝑇 )𝛿𝑝 + = (A.5) 𝜕𝑿 𝜕𝝃 𝑟 ∶ 𝛿𝝃 𝑟
	A.2. Computation of the stress
	A.2.1
																						(A.6b)
										or, equivalently:
										{	𝛿𝑅 𝑝 𝛿𝑹 𝝃 𝑟	}	=	(	𝐾 𝑝𝑝		𝑲 𝑝𝝃 𝑟	)	×	{	𝛿𝑝 𝛿𝝃 𝑟	}	+	(	𝑲 𝑝𝜺 K 𝝃 𝑟 𝜺 𝑲 𝝃 𝑟 𝑇 𝐾 𝑝𝑇	)	×	{	𝛿𝜺 𝛿𝑇	}	, (A.7)
										K 𝝃 𝑟 𝜺 = 0,									(A.8g)
																						(A.8h)
										with
										𝐴 𝑝𝑝 =		𝜕𝑅 𝑝 𝜕 ṗ 𝜕 𝜕𝑝 ṗ	=	1 𝛥𝑡	,			(A.9a)
										𝑨 𝑝𝝈 =	𝜕𝑅 𝑝 𝜕𝝈	= -	𝑁 𝜂 𝜎 𝑦 (𝑇 )	[ ⟨	⟩ 𝑓 (𝝈, 𝑿, 𝑇 ) 𝜎 𝑦 (𝑇 )	+	] 𝑁-1	𝜦(𝝈, 𝑿),	(A.9b)
										𝑨 𝑝𝑿 =	𝜕𝑅 𝑝 𝜕𝑿	=	𝑁 𝜂 𝜎 𝑦 (𝑇 )	[ ⟨	⟩ 𝑓 (𝝈, 𝑿, 𝑇 ) 𝜎 𝑦 (𝑇 )	+	] 𝑁-1	𝜦(𝝈, 𝑿),	(A.9c)
										𝐴 𝑝𝑇 =	𝜕𝑅 𝑝 𝜕𝑇	=	𝜂 𝑁	[ ⟨	𝑓 (𝝈, 𝑿, 𝑇 ) ⟩
													𝜕𝑹 𝝃 𝑟 𝜕 ξ𝑟	∶	𝜕 ξ𝑟 𝜕𝝃 𝑟	=	I 𝛥𝑡	,	(A.9e)

where

𝐾 𝑝𝑝 = 𝐴 𝑝𝑝 + 𝑨 𝑝𝝈 ∶ 𝑩 𝝈𝑝 + 𝑨 𝑝𝑿 ∶ 𝑩 𝑿𝑝 , (A.8a

)

𝑲 𝑝𝝃 𝑟 = 𝑨 𝑝𝑿 ∶ B 𝑿𝝃 𝑟 , (A.8b

)

𝑲 𝑝𝜺 = 𝑨 𝑝𝝈 ∶ B 𝝈𝜺 , (A.8c) 𝐾 𝑝𝑇 = 𝑨 𝑝𝝈 ∶ 𝑩 𝝈𝑇 + 𝑨 𝑝𝑿 ∶ 𝑩 𝑿𝑇 + 𝐴 𝑝𝑇 , (A.8d) 𝑲 𝝃 𝑟 𝑝 = A 𝝃 𝑟 𝑿 ∶ 𝑩 𝑿𝑝 , (A.8e) K + 𝜎 𝑦 (𝑇 ) ] 𝑁-1 𝜎 ′ 𝑦 (𝑇 ) [ 𝜎 𝑦 (𝑇 ) + 𝑓 (𝝈, 𝑿, 𝑇 ) ]

𝜎 𝑦 (𝑇 ) 2 , (A.9d)

A 𝝃 𝑟 𝝃 𝑟 =

Then 𝛿𝑝 and 𝛿𝝃 𝑟 are obtained from the nullity condition of the two residuals 𝑅 𝑝 and 𝑹 𝝃 𝑟 in their linearized form (A.7), while considering 𝛿𝜺 = 𝟎 and 𝛿𝑇 = 0 (all the quantities are taken at (𝑛+1)(𝑘) ):

where

while, according to the convex cutting plane method, 𝛿𝜺 𝑝 and 𝛿𝝃 𝑝 are calculated using Eqs. (A.3). This correction loop is therefore repeated until reaching convergence, namely, when both residuals 𝑅 𝑝 and 𝑹 𝝃 𝑟 are sufficiently close to zero.

A.3. Tangent operators

The tangent operators d𝝈 d𝜺 and d𝝈 d𝑇 are also required for the FE solver to compute the global predictor at the whole FE system's level. Since, the proposed implementation is based on the convex cutting plane form of the return mapping algorithm [START_REF] Simo | Computational Inelasticity[END_REF]Praud et al., 2017a,b;[START_REF] Praud | Multi-Scale Modelling of Thermoplastic-Based Woven Composites[END_REF][START_REF] Chatzigeorgiou | Thermomechanical Behavior of Dissipative Composite Materials[END_REF], it utilizes continuous tangent operators. To compute them, one needs to take back the linearized stress:

If the material is not actively yielding, i.e., the computation of the stress stopped at the elastic prediction with hardening recovery correction, then d𝑝 = 0 and the tangent operators are simply given by: d𝝈 d𝜺 = B 𝝈𝜺 and d𝝈 d𝑇 = 𝑩 𝝈𝑇 . However, if the material is actively yielding, i.e., the computation of the stress went through the viscoplastic and hardening recovery correction, then, one must consider the nullity condition of the two residuals 𝑅 𝑝 and 𝑹 𝝃 𝑟 in their linearized forms (A.7), which eventually leads to: .16) or, equivalently:

where


WICKY Lucie

PhD Student in Sociology

France lucie@wicky.fr

Background

Research on Gender-based violence provides detailed understanding of exposure, contexts and consequences of violence experienced by women. Power relations define violence, and in particular gender, race, class and sexuality. Knowledge about violence against men, and more specifically sexual violence, need to be developed. This is particularly the case in France: there is no in-depth research that describes sexual violence experienced by men during their lives.

Moreover, available research is scattered in different fields of investigation: child sexual abuse, violence against bisexual and gay men, intimate partner violence, domestic violence, etc.

Results

Different contexts of violence before the age of 18, depending to age of exposure Childhood: gender and age domination Adult men (of family) perpetrators: gender and age relations are intertwined, often with family hierarchies Qualification of violence as such during adulthood facilitated by status of "child" Teenage years: gender, sexuality and age domination Situations of male socialization that are difficult to qualify as violence (hétérosexuals' men) First experiences of sexuality between men, not qualify neither

Innovation

This study is the first research in France to look in detail at sexual violence suffered by men during their life.

The choice to consider all sexual violence together, irrespective of age, to understand the evolution of men's exposure to violence according to the context of violence and the power relations at work against them.

As men are not affected by violence throughout their lives, but at younger ages, the question of power dynamics between younger and older individuals is investigated. There is only few research in the field of age studies or aging in the case of childhood and in the context of violence suffered by men.

Objective

To In France, 3.14% of men declared sexual violence during their lives, mainly during childhood and adolescence and committed by other men. What are the conditions of exposure to violence for men? What are the different power relation that work against men? How are these situations of domination illustrated in the discourses of men themselves?

12-17 years old: Violence during teenage years

Same-sex experiences

Beggining of same-sex experiences Principally gay men "Confusion" or "non-consent situation": difficulty to speak of "violence" Protect the whole LBGT community Limiting impact on the already subordinated masculinity Sexualisation of the age difference

Spaces of socialization between men

  describe the contexts of exposure and the power relations at work according to the age at sexual violence experienced by men Qualitative data: biographical interviews 50 men who declared sexual violence in VIRAGE survey (follow-up interviews by phone) Random selection from violence within family and close social circle and/or within other spaces (incl. public spaces) A lager sample of male victims, diversity of experiences of sexual violence and diversity of sociodemographic profiles Life history calendar: trajectories of violence, childhood, conjugality, sexuality, education and work, health, norms to violence, gender and sexuality Age distribution of sexual violence experienced by men (%) Percentage Age at (first) violencePerpetrators by their gender, adulthood and family relationship, by age at experience of violence to male victims (%)

11 years old 12-17 years old Adult men perpetrators 57.7 27.8

																	0-11 years old: Violence during Childhood
	15
	10																When things happened, I wasn't able to understand what was going on. You can't
																	force to have sex with a child
																	who doesn't know what it is
	5
	0
	3	5	7	9	1 1	1 3	1 5	1 7	1 9	2 1 2 3	2 5	2 7	2 9	3 5	4 5	5 3	6 0
	0-Adult men within the family	46.3					6.9
	Adult men outside the family				11.4					20.9
																	What happened?
	Other perpetrators						42.3					72.2
	Other perpetrators within the family			24.7					14.7
	Other perpetrators outside the family		17.7					57.5
	405 men report sexual violence
	during their lives on VIRAGE

Older men perpetrators Age domination linked to family hierarchies Identify easily "violence" as such during adulthood: difference between the "child" victim he was and man he is today

	aged 20-69.
	Life spaces: within school or university, within the workplace, or committed by a partner or ex-partner,
	within the family or close social circle and, finally, within others spaces (including public spaces).
	"Violence" or "assault" never used: description purely factual
	Questions about violence declared: frequency, age, perpetrators, etc.

Sexual Violence against Men and Power RElations: Understanding contexts of exposure 4. Methods Quantitative data: Violence and Gender Relations survey ("VIRAGE") Conducted by the French Institute of Demographic Studies (INED) in 2015. A sample of almost 16,000 women and 12,000 men representative of the French population

School, hobbies, rarely sexual experience Heterosexuals' (and bisexuals') men Without mentioning "sexual violence" A teenager cannot be a victim To refuse violence limiting impact on masculinity Temporality of violence structuring the contexts Power relation related to age and gender depending to age of exposure

  I would rather say... where things, well... we feel a little bit in a kind of... how can I put it? in a corridor, it has to come, here. […] A little bit of that... that is… a confusion .

	Obligations to do things to an adult… […] You have to masturbate your grandfather at six	He didn't force it, but I still had, I don't know, some kind of moral obligation to listen to my big brother. I mean, I was still too young to understand
		and realize what was
		going on
	Philippe, 59 years			Bernard, 64 years old,
	old, assaulted by			assaulted from 3 to 13
	his grandfather			years old by a cousin (14
	from 6 to 11 years old	Bastien, 39 years old, raped from 11 to 15 years old by his elder brother	years old), a priest and a teacher
	Have you experienced psychological,
	physical or sexual violence?
	Sexual, no. Well, yes, well… it wasn't really violence, but... You see, in boarding school, it	Have you ever accepted sex without wanting to?
	is a contest between guys,
	and we do it to have peace, you know what I mean.	Without wanting to? Not at all, no. I don't know, maybe it's not very precise, but I think of things…
			and in particular gay encounters, clearly, where
			it's implied that you don't dare to say that you're
			not sure. […] But not... I have trouble saying "non-
	It was big willy competition. Masturbation competition. Oh, wait yes! About violence, we also had abs exercise forced sessions.	consent", Jordan, 31 years old, sexual experiences with other men
				(16-25 years old)


B Sauvage

F Miralles

S Wornom

Bruno Koobus

Alain Dervieux

About mesh adaptation for hybrid flow calculation

Keywords: turbulence, RANS, DDES, hybrid RANS/LES, Variational Multi-scale

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

We propose to address the increasing complexity of industrial flows by extending the predictivity of hybrid models by the introduction of intermittency together with making easier the research of efficient and converged space and time approximations by introducing adaptation algorithms.

Modelization:

The base ingredients of our hybrid turbulence strategies are:

• Baseline RANS components: two low Reynolds RANS models are used in our hybrid models, more specifically they are the k -ε model proposed in [START_REF] Goldberg | A wall-distance-free k -ε model with Enhanced Near-Wall Treatment[END_REF] and the k -R model recently introduced by [START_REF] Zhang | Development of k-R turbulence model for wall-bounded flows[END_REF]. They have been chosen for their abilities to properly predict separated flows with adverse pressure gradients.

• Transition model for boundary layer: Let P k = τ : ∇u and D k = ρε be the turbulent kinetic energy production and destruction terms. Let us define the laminar and turbulent flow regions in the boundary layer:

L = {x ∈ Ω f | Re θ (x) < Re θ,S (x)}, T u = complement(L),
where Re θ denotes the Reynolds number based on the boundary layer thickness which is defined by Re θ = 0.664 Re|x|.

Re

θ,S = 163 + exp F λ - 100F λ 6.91 2 3 k .
with

F λ = 6.91 + 2.48λ -12.27λ 2 , λ > 0, 6.91 + 2.48λ + 63.64λ 2 , λ < 0. and λ = 0.664 2 1 Re ∂u x ∂x (x -x 0 )
where x 0 denotes the abscissa of the front body point. Following the work of Menter et al. (2015) and the works of [START_REF] Akhter | Development of Prediction Method of Boundary Layer Bypass Transition using Intermittency Transport Equation[END_REF][START_REF] Akhter | Numerical Simula-tion of Bypass Transition by the Approach of Intermittency Transport Equation[END_REF][START_REF] Akhter | Numerical simulation of heat transfer coefficient on turbine blade using intermittency factor equation[END_REF], the present transition model is defined as follows :

∂ρk ∂t + ∇ • (ρuk) = P k -D k + ∇ • [(µ + µ t σ k ) ∇k] ∂ρε ∂t + ∇ • (ρuε) = (c ε1 P k -c ε2 D k + E) T -1 t + ∇ • [(µ + µ t σ ε ) ∇ε] ∂ργ ∂t + ∇ • (ρuγ) = c g1 γ(1 -γ) P k k + ρc g2 k 2 ε ∇γ • ∇γ+ ∇ • [σ γ (µ + µ t )∇γ] .
The intermittency model interacts with the turbulence model by modifying the turbulent kinetic energy equation. The new production and destruction terms are defined by : where the model constants are defined by γ 1 = 0 and γ 2 = 0.1. From the above equations, one can notice that the baseline RANS model is recovered for an intermittency value γ = 1 (fully turbulent mode). A zero normal flux is also imposed on γ at the wall.

P k = 0 if x ∈ L, max (γ, γ 1 ) P k otherwise. D k = 0 if x ∈ L, max (γ, γ 2 ) D k otherwise.
• DDES component: in this work, the classical DDES approach , [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF], is based either on the Spalart-Allmaras model or the k -ε model of [START_REF] Goldberg | A wall-distance-free k -ε model with Enhanced Near-Wall Treatment[END_REF].

• LES-like component: the DVMS approach proposed in Moussaed et al. (2014) is used as the LES part of our hybrid models. In this approach, the variational multiscale (VMS) model, aiming to limit the effects of the subgrid-scale (SGS) model to the smallest resolved scales, is combined with the dynamic procedure which provides a tuning of the SGS dissipation in space and time, so that the resulting DVMS model enjoys synergistic effects.

• Hybrid models Our hybrid strategies blend either a RANS or DDES model with the DVMS approach.

3 Mesh adaptative discretization :

The baseline discretization is an upwind, vertex centered, finite-volume approximation on tetrahedra. Diffusion terms are accounted for with the P 1 -Galerkin approximation. The numerical diffusion involved in upwinding is made of sixth-order derivatives of primitive variables.

The capture of high Reynolds number flows with a mesh adaptation algorithm demands an algorithm able to find the boundary layer when starting from a uniform flow. This is obtained by using (1) an adequate adaptation criterion and (2) a stable and convergent fixed point for the coupling of flow and mesh. Main principles for a successful adaptation can be found in [START_REF] Dervieux | Mesh adaptation for Computational Fluid Dynamics[END_REF]. The transient fixed point (TFP) was introduced and discussed in [START_REF] Alauzet | Adaptation de maillages non structurés pour des problèmes instationnaires[END_REF] and [START_REF] Alauzet | 3D transient fixed point mesh adaptation for timedependent problems: Application to CFD simulations[END_REF]. The simulation time frame [0, T ] is split into several subintervals (Figure 1):

[0, T ] = [0 = t 0 , t 1 ]∪...∪[t i , t i+1 ]∪...∪[t n-1 , t n = T ].
The TFP is extended as in [START_REF] Sauvage | A spacetime mesh Transient Fixed Point for mesh adaptation[END_REF] by defining a space-time continuous mesh (M, τ, nstep) as the knowledge of the following ingredients:

-a number nstep of time intervals.

-a timestep function:

t ∈]0, T [ → τ (t) ∈]0, T [ valid in the sense that T 0 (τ (t)) -1 dt = nstep.
-for every t ∈]0, T [ a spatial metric M(t), constant on each subinterval, of spatial complexity n(t) = C sp (M(t)).The space-time complexity C(M, τ, nstep) of a space-time continuous mesh (M, τ, nstep) is:

C(M, τ, nstep) = T 0 C sp (M(t))(τ (t)) -1 dt. (1)
To any space-time continuous mesh (M, τ, nstep) corresponds an estimate of the resulting error of a calculation working on the space-time mesh build with it. With a BDF2 time advancing scheme and a secondorder spatial approximation, this estimate can be (in practice we choose p = 4):

E 0 (M h , τ ) = T 0 Ω (τ 2 | ∂ 3 u ∂t 3 |+∆x H u ∆x p dtdx.
(2) The sensor u is in our examples taken as the local Mach number of the computed flow. Minimizing the error E 0 (M h , τ ) with respect to (M h , τ ) under the constraint of a fixed space-time complexity C(M, τ, nstep) = N prescribed can be done analytically, allowing the building of a new space-time mesh. The process is repeated in the TFP loop.

In the approach chosen here, we restrict to LES flows for which a global shedding period T s can be estimated, and we use T s as the size of the subintervals. Due to quasi-periodicity, the different subinterval meshes are very similar, and we reduce them to an unique mesh in which quasi-steady features like boundary layers are accurately followed by the mesh, while unsteady vortices travel in a rather uniformly refined region.

Impact of the intermittency model:

The new intermittency model has been tested on the flow past a cylinder at Reynolds number around the drag crisis. At Reynolds 3900, the physical boundary layer is laminar, the intermittency model inhibits the statistical model and numerically predicts the laminar boundary layer (Figure 2). The quality of the prediction is of same level as with a pure LES calculation, e.g. by Moussaed et al. (2014). At Reynolds 380K, the physical boundary layer presents a dissymmetry between top and bottom, one layer being laminar while the other one is turbulent. This is well predicted with the help of the intermittency model, Figure 3. At Reynolds 1M the physical boundary layers are both turbulent and are again well predicted by the proposed model, Figure 4. The improvement carried in drag prediction by the intermittency equation is put in evidence in Figure 5.

The new intermittency model is now applied to the flow around a 3D wing based on the NACA0018 air- foil. At zero angle of attack, the slight vortex shedding is captured, Figure 6. For the intermediate α = 6 degrees case, the unstable extrados boundary layer is also predicted (Figure 7). The α = 15 degrees case is depicted in Figure 8 and global behavior of the prediction is in accordance with experiments of [START_REF] Nakano | Experimental study on flow and noise characteristics of NACA0018 Aifoil[END_REF], of [START_REF] Du | Numerical and experimental investigations of Darrieus wind turbine start-up and operation[END_REF], and of [START_REF] Boutilier | Effects of end plates and blockage on low-reynolds-number flows over airfoils[END_REF].

Mesh adaptative computations

We essentially present experiments performed with 2D test cases.

• Flow around a cylinder

We present 2D computations of a flow around a cylinder at Reynolds number 3900 with the Spalart-Allmaras turbulence model. Mesh adaptation options are : -only one time interval, therefore only one adapted spatial mesh.

-the Space-Time complexity N st is prescribed succes- sively to 10M, 20M, 100M, 200M. An example of adapted mesh is given in Figure 9. Space-time statistics are presented in Table 1. We ob- Table 1: Space-time statistics for circular cylinder case at Reynolds number 3900. k is the number of fixed points, Csp holds for the spatial complexity (≈ number of space nodes), # ∆t for the number of time steps, Esp and Etime for the resulting space and time error integrals.

N st k C sp # ∆t E sp E time 10M 1 27K 361
serve that space error and time error tend to equilibrate, but yet imperfectly.

• Flow around an airfoil

The second case is the flow around a NACA0021 at Reynolds number 270K and an angle of attack 60 degrees, again with the Spalart-Allmaras turbulence model. Mesh adaptation options are : -only one adapted spatial mesh, -Space-time complexity is prescribed to 10M, 20M, 100M, 200M. An example of adapted mesh is presented in Figure 10. Space-time statistics are presented in Table 2.

In Figure 11, the adaptive time step (function of physical time) produced by the space-time optimization is depicted, showing the increase of time-step, reducing computational cost while preserving the global accuracy.

N st k C sp # ∆t E sp E time 2M
1 1766 1132 1.7 10 -1 2.4 10 In Figure 12 we present the evolution (with Transient Fixed point Iterations) of the two components of the error, converging to each other. The novel mesh adaptation approach has been tested on the 3D unsteady flow around a NACA0018 at an angle of attack of 8 degrees. The VMS model is applied and produces a rather unsteady flow for which the adaptation is able to both concentrate anisotropicly on boundary layers and manage a rather uniform for large structures propagation, Figures 13,14.

The method is also being applied to a MRF flow [START_REF] Chargy | A mesh adaptative method for rotating machines[END_REF]

Figure 1 :

 1 Figure 1: Time splitting of the Transient Fixed Point mesh adaptation algorithm. Sub-intervals (in green) used for the transient process and timesteps (in red).

Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 2345 Figure 2: Flow around a cylinder: impact of the intermittency model on the drag crisis prediction, vorticity at Rey=3900

Figure 6 :

 6 Figure 6: Flow around NACA0018 wing, α = 0 degree.

Figure 7 :

 7 Figure 7: Flow around NACA0018 wing, α = 6 degrees.

Figure 8 :

 8 Figure 8: Flow around NACA0018 wing, α = 15 degrees.

Figure 9 :

 9 Figure 9: Cylinder at Reynolds number 3900. Example of adapted mesh.

Figure 10 :

 10 Figure 10: Naca0021 at high angle of attack. Example of the adapted mesh.

Figure 11 :

 11 Figure 11: Bottom: Timestep lengths of initial flow at CFL=50, 1132 timesteps on 1766 vertices, and top: first timestep lengths proposed by the adaptation algorithm, 123 timesteps on 16K vertices.

Figure 13 :

 13 Figure 13: Flow around NACA0018 at 15 degrees angle of attack. LES-adapted mesh with an instantaneous flow obtained with it (velocity magnitude).

Figure 14 :

 14 Figure 14: Flow around NACA0018 at 15 degrees angle of attack. LES-adapted mesh with an instantaneous flow obtained with it (velocity magnitude).

Figure 12 :

 12 Figure 12: Evolution of theorical space error Esp (top) and time error Etime (bottom) with TFP iterations.

  around a Caradonna-Tung wing rotating above a Robin helicopter shape, Figures15,16

  .

Figure 15 :

 15 Figure 15: Flow around a Caradonna-Tung wing rotating above a Robin helicopter. LES-adapted mesh.

Figure 16 :

 16 Figure 16: Flow around a Caradonna-Tung wing rotating above a Robin helicopter. An instantaneous flow obtained with the adapted mesh (velocity magnitude).

Table 2 :

 2 Space-time statistics for NACA0021 at Reynolds number 270K.

Acknowledgments

This work was supported by the ANR NORMA