Thomas Heuzé
email: thomas.heuze@ec-nantes.fr

Laurent Stainier

A variational formulation of thermomechanical constitutive update for hyperbolic conservation laws

Keywords: Variational constitutive update, Hyperbolic conservation laws, Dissipative solid media, Thermo-mechanics, Flux-difference splitting finite volume method, Thermo-hyper-elasto-viscoplasticity

In this paper, a variational framework is proposed for the constitutive update of thermomechanical constitutive models in the special case where their input results from quantities directly updated by hyperbolic conservation laws. Both a continuum and a consistent first order accurate discrete settings are derived. The originality of this work lies in that the constitutive update is driven by the rates of some strain measure and the internal energy density in the continumm setting, leading to a rate-type description of the local constitutive problem, and by the updated values at some discrete time of these strain measure and internal energy density in the discrete setting. These quantities are updated by the solution of a system of discrete conservation laws including the first principle of thermodynamics, ensuring that the right shock speeds will be computed. This point is of crucial importance when simulating impact on structures for instance. The proposed variational approach is illustrated for thermo-hyperelastic-viscoplastic solid media, especially using the parameterization of the flow rule direction based on pseudo-stresses proposed by Mosler & co-workers. The proposed discrete variational solver is then coupled with the second order accurate flux difference splitting finite volume method, which permits to solve the set of conservation laws. Comparisons are performed on a set of test cases with numerical solutions obtained with finite elements coupled to an explicit time-stepping and to a temperature-driven variational constitutive update. They allow to show the good behavior of the proposed approach.

Introduction

Many engineering applications require to perform numerical simulations of impact on structures involving dissipative solids. This is the case in a variety of important fields of application: crash, high-velocity impact such as ballistic penetration, or high-speed forming processes. In these applications, the numerical simulations often involve the solution of both mechanical and thermal effects, generally in a coupled way. Different mechanisms of thermomechanical coupling may occur, but the main one in the aforementioned applications results from large inelastic strains which are generated and contribute to both the temperature rise and the drop of the mechanical strength due to thermal softening. Traditionally, authors prefer to solve thermal phenomena using the temperature T as main thermal unknown [START_REF] Simo | Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation[END_REF], although the entropy is also sometimes used [START_REF] Armero | A new unconditionally stable fractional step method for non-linear coupled thermomechanical problems[END_REF]. This is so mainly because the temperature is a measurable quantity, and this has led to develop the modelling of thermal phenomena as a direct function of that quantity. Especially, thermomechanically coupled constitutive models are written for dissipative media in solid mechanics using Helmholtz's free energy [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF][START_REF] Ziegler | An introduction to thermomechanics[END_REF], which is a function of the temperature, while their constitutive updates are driven with some strain measure and also the temperature [START_REF] Simo | Computational inelasticity[END_REF][START_REF] Simo | Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation[END_REF][START_REF] Miehe | Entropic thermoelasticity at finite strains. aspects of the formulation and numerical implementation[END_REF], obtained from the solution of balance equations. Generally, the latter consist of the heat equation and the linear momentum balance [START_REF] Simo | Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation[END_REF], whose weak forms are solved using any discretization method for the spatial part like finite elements [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF][START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF] or particle methods [START_REF] Li | Meshfree and particle methods and their applications[END_REF], following either Lagrangian [START_REF] Camacho | Adaptive lagrangian modelling of ballistic penetration of metallic targets[END_REF][START_REF] Marusich | Modelling and simulation of high-speed machining[END_REF][START_REF] Ba | Thermomechanical total lagrangian sph formulation for solid mechanics in large deformation problems[END_REF], Arbitrary-Lagrangian-Eulerian [START_REF] Beni | Consistent arbitrary lagrangian eulerian formulation for large deformation thermo-mechanical analysis[END_REF][START_REF] Crutzen | Lagrangian and arbitrary lagrangian eulerian simulations of complex roll-forming processes[END_REF][START_REF] Feulvarch | A simple and robust moving mesh technique for the finite element simulation of friction stir welding[END_REF][START_REF] Bussetta | Two 3d thermomechanical numerical models of friction stir welding processes with a trigonal pin[END_REF] or Eulerian [START_REF] Kumar | Eulerian thermo-mechanical simulations of heterogeneous solid propellants using an approximate projection method[END_REF] descriptions, plus some explicit time discretization scheme. In very fast processes, heat conduction effects can be neglected, then the heat equation can be solved locally on the temperature increment.

However in simulations of impact phenomena, it is well known that the solution of the hyperbolic initial boundary value problem may consist of both continuous and discontinuous waves [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF]. The accurate capturing of discontinuous waves, especially shock waves, is of crucial importance to properly understand the mechanical phenomena occuring within that medium. For solid-type media, a correct capturing of plastic waves will also allow for an accurate assessment of the plastic strain field and hence that of residual stresses and distortions within the structure [START_REF] Trangenstein | A higher-order Godunov method for modeling finite deformation in elastic-plastic solids[END_REF][START_REF] Miller | A high-order eulerian godunov method for elastic-plastic flow in solids[END_REF][START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF][START_REF] Heuzé | Lax-wendroff schemes for elastic-plastic solids[END_REF][START_REF] Renaud | On loading paths followed inside plastic simple waves in two-dimensional elastic-plastic solids[END_REF]. But the heat equation traditionally solved when considering the temperature as main thermal unknown is only valid for smooth solutions. When discontinuous solutions occur, a set of conservations laws including the first principle of thermodynamics written on the total energy should rather be considered. This set of conservation laws reduces to the well-known Rankine-Hugoniot jump conditions across any discontinuity [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF], and ensures that the right shock speeds will be computed once this system will have been discretized, while the use of the heat equation does not ensure it. If the first principle of thermodynamics is to be used in place of the heat equation, the temperature is not anymore the input variable of the thermomechanical constitutive update, and the design of such constitutive update driven with quantities obtained from conservation laws should be considered.

The variational framework has emerged to be very convenient to build thermodynamically consistent and numerically efficient constitutive updates, which are also known as variational constitutive updates [START_REF] Comi | Extremum properties of finite-step solutions in elastoplasticity with nonlinear mixed hardening[END_REF][START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF]. Following the set of Standard dissipative solids [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF], hence taking advantage of the description of the constitutive response through both state and dissipation (pseudo-)potentials, the local constitutive problem was recast into an equivalent optimization problem, from which many benefits can be exploited. Especially, the existence of solutions can be analyzed by using the same tools originally designed for hyperelastic material models [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF], well-developed optimization algorithms can be reused, the Hessian matrix is now symmetric if a Newton scheme is used, and the minimum/maximum nature of the optimal point of the functional can be used as an error indicator [START_REF] Mosler | Variational h-adaptation in finite deformation elasticity and plasticity[END_REF][START_REF] Mosler | An error-estimate-free and remapping-free variational mesh refinement and coarsening method for dissipative solids at finite strains[END_REF][START_REF] Pethe | Variational h-adaption for coupled thermomechanical problems[END_REF][START_REF] Pethe | Remapping-free variational h-adaption for strongly coupled thermo-mechanical problems[END_REF] for adaptive methods. Especially, Ortiz, Stainier and co-workers [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF][START_REF] Fancello | A variational constitutive update algorithm for a set of isotropic hyperelastic-viscoplastic material models[END_REF][START_REF] Stainier | Study and validation of a variational theory of thermo-mechanical coupling in finite visco-plasticity[END_REF][START_REF] Stainier | Consistent incremental approximation of dissipation pseudo-potentials in the variational formulation of thermo-mechanical constitutive updates[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF] introduced a variational formulation of the coupled thermomechanical boundary-value problem for general dissipative solids. This principle relied on the one hand on a two-field thermal formulation (an external temperature T , appearing in the heat equation, and an internal temperature Θ obtained through a state law), whose equality is enforced as an internal constraint in the constitutive model, and relaxed in state laws. On the other hand, an integration factor T/Θ weighting rate arguments of the pseudo-dissipation potential was identified by considering a time rescaling of these quantities, that allowed to find a variational form to the general thermomechanical rate problem, and recover the requisite symmetry of the strong form. The associated incremental variational update was consistent with standard finite element codes, separating the local time-discrete constitutive update performed pointwise at the integration points, from the solution of an incremental boundary value problem giving the unknown deformation mapping and external temperature on a mesh. An extension of this work to non-associated evolutions equations, especially to nonlinear kinematic hardening, and to any yield function being positively homogeneous of degree one was then proposed by Mosler and co-workers [START_REF] Canadija | On the thermomechanical coupling in finite strain plasticity theory with non-linear kinematic hardening by means of incremental energy minimization[END_REF][START_REF] Bartels | On the thermomechanical coupling in dissipative materials: a variational approach for generalized standard materials[END_REF][START_REF] Canadija | A variational formulation for thermomechanically coupled low cycle fatigue at finite strains[END_REF]. Following former works developed in the isothermal setting by these authors [START_REF] Mosler | Towards variational constitutive updates for non-associative plasticity models at finite strain: models based on a volumetric-deviatoric split[END_REF][START_REF] Mosler | On the implementation of rate-independent standard dissipative solids at finite strain-variational constitutive updates[END_REF][START_REF] Mosler | Variationally consistent modeling of finite strain plasticity theory with non-linear kinematic hardening[END_REF][START_REF] Bleier | Efficient variational constitutive updates by means of a novel parameterization of the flow rule[END_REF], the potentially non-associative evolution equations are a priori enforced by employing a suitable parameterization of the flow rule and the evolution equations using pseudo-stresses, yielding an unconstrained optimization problem. Such parameterization was also shown to reduce the numerical complexity of variational constitutive update with respect to classical return-mapping algorithms [START_REF] Bleier | Efficient variational constitutive updates by means of a novel parameterization of the flow rule[END_REF]. Besides, being thermodynamically consistent, these approaches allow to correctly model and compute the non-constant partition of plastic work into heat and stored energy [START_REF] Hodowany | Partition of plastic work into heat and stored energy in metals[END_REF][START_REF] Rosakis | A thermodynamic internal variable model for the partition of plastic work into heat and stored energy in metals[END_REF][START_REF] Stainier | Study and validation of a variational theory of thermo-mechanical coupling in finite visco-plasticity[END_REF], in contrast with empirical approaches inspired by the pioneering work of Taylor and Quinney [START_REF] Taylor | The latent energy remaining in a metal after cold working[END_REF], and followed in some other works [START_REF] Simo | Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation[END_REF][START_REF] Wriggers | On the coupled thermomechanical treatment of necking problems via finite element methods[END_REF].

In the present work, we introduce a particular variational framework for the constitutive update of thermomechanical constitutive models which are driven by quantities directly updated by hyperbolic conservation laws. Both a continuum and a consistent first order accurate discrete settings are derived. In the continuum setting, the constitutive update is driven by the rates of some strain measure and the internal energy density, leading to a rate-type description of the local constitutive problem. More precisely, a Lagrangian functional is built from the one initially considered in [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF], and from the rate of the residual of the Legendre tranform of Helmholtz'free energy enforced to vanish through a Lagrange multiplier. In the discrete setting, a first order accurate discrete variational constitutive update is derived, whose update is driven by the increments of these strain measure and internal energy density. The parameterization of the flow rule and evolution equations is performed with the concept of pseudo-stresses, following the work of Mosler and co-workers [START_REF] Bleier | Efficient variational constitutive updates by means of a novel parameterization of the flow rule[END_REF]. The proposed variational approach is illustrated with thermo-hyperelastic-viscoplastic solid media, especially with the Johnson-Cook viscoplastic flow rule [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF]. The main advantage of this new variational formulation is that it is naturally compatible with any numerical scheme dedicated to the approximation of the solution of hyperbolic systems, like finite volumes [START_REF] Favrie | Mathematical and numerical model for nonlinear viscoplasticity[END_REF][START_REF] Maire | A nominally second-order cell-centered Lagrangian scheme for simulating elastic-plastic flows on two dimensional unstructured grids[END_REF][START_REF] Aguirre | An upwind vertex centred finite volume solver for lagrangian solid dynamics[END_REF], particle methods [START_REF] Ghavamian | An entropy-stable smooth particle hydrodynamics algorithm for large strain thermo-elasticity[END_REF], Discontinuous Galerkin approaches [START_REF] Cockburn | Runge-kutta discontinuous galerkin methods for convection-dominated problems[END_REF][START_REF] Busto | High order ader schemes for continuum mechanics[END_REF], or finite element method [START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part iii: Thermo-elasticity[END_REF], written in conservation form. Here, the second order accurate flux difference splitting finite volume method [START_REF] Leveque | Wave propagation algorithms for multidimensional hyperbolic systems[END_REF][START_REF] Heuzé | Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF] is used, which permits to solve the set of conservation laws on examples involving either small or large strains. Finally, a set of numerical tests are conducted to show the good behavior of the proposed approach, especially comparisons are performed with numerical solutions obtained with classical finite elements coupled with an explicit time-stepping and with a temperature-driven variational constitutive update [START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF][START_REF] Stainier | ZorgLib. User's manual[END_REF].

Thermo-mechanical initial boundary value problem

Conservation laws

Geometrical conservation law

We consider a continuum body Ω, whose motion is described by the mapping φ(X, t). This mapping relates the position of a material point of coordinates X ∈ Ω 0 in the initial configuration, to its current coordinates x ∈ Ω(t). In the initial configuration, the domain Ω 0 has a boundary denoted ∂Ω 0 , of outward unit normal N, while in the current configuration, the domain Ω has a boundary denoted ∂Ω and n is the associated outward unit normal. Based on this mapping, the rate of the deformation gradient two-point tensor Ḟ reads as the material gradient of the material velocity vector

v Ḟ(X, t) = GRAD v ∀X ∈ Ω 0 , (1)
which can be rewritten as a geometrical conservation law [START_REF] Plohr | A conservative eulerian formulation of the equations for elastic flow[END_REF][START_REF] Trangenstein | A higher-order Godunov method for modeling finite deformation in elastic-plastic solids[END_REF][START_REF] Lee | Development of a cell centred upwind finite volume algorithm for a new conservation law formulation in structural dynamics[END_REF]]

Ḟ -DIV(v ⊗ 1) = 0. (2)
The material divergence DIV is computed with respect to initial coordinates X, and 1 is the identity of second order. From Equation (1), it is evident that the material CURL of the deformation gradient F vanishes CURL F = 0.

(3)

Remark 1. The above system of conservation laws can be extended if required with conservation laws written on the complementary minors of the deformation, as shown in the work of Bonet, Gil and co-workers [START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part iii: Thermo-elasticity[END_REF][START_REF] Ghavamian | An entropy-stable smooth particle hydrodynamics algorithm for large strain thermo-elasticity[END_REF]. These minors consist of the deformation gradient F, its cofactor H and the jacobian determinant J. Following the notations of these authors, conservation laws associated with the two latter read

Ḣ -CURL(v × F) = 0 (4) J -DIV(H T • v) = 0. ( 5
)
Such extension is interesting when considering nearly incompressible and truly incompressible materials [START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part i: Total lagrangian isothermal elasticity[END_REF][START_REF] Gil | A first order hyperbolic framework for large strain computational solid dynamics. part ii: Total lagrangian compressible, nearly incompressible and truly incompressible elasticity[END_REF][START_REF] Haider | A first-order hyperbolic framework for large strain computational solid dynamics: an upwind cell centred total lagrangian scheme[END_REF] governed by a polyconvex constitutive law [START_REF] Bonet | A computational framework for polyconvex large strain elasticity[END_REF]. The derivation of the variational principle in this work is not written as a function of the three minors of the deformation, but only as a function of the deformation gradient. However, its extension is straightforward and the formulation could thus be applied to the schemes proposed in [START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part iii: Thermo-elasticity[END_REF][START_REF] Ghavamian | An entropy-stable smooth particle hydrodynamics algorithm for large strain thermo-elasticity[END_REF].

Other conservation laws

The conservations of the linear momentum and of the total energy (or first principle of thermodynamics) read in their material form as

∂p ∂t -DIV P = 0 (6) Ė -DIV(P T • v) = 0 (7)
where p = ρ 0 v denotes the density of linear momentum per unit undeformed volume, ρ 0 (X) = ρ(X, t = 0) is the reference mass density, P denotes the first Piola-Kirchhoff stress tensor, and E is the total energy density defined as

E = E + ρ 0 v 2 2 ( 8
)
where E refers to the internal energy density and K = ρ 0 v 2 /2 to the kinetic energy density. Equations ( 6) and ( 7) are here written without any source term. Especially no heat flux density vector is considered in Equation [START_REF] Beni | Consistent arbitrary lagrangian eulerian formulation for large deformation thermo-mechanical analysis[END_REF], which amounts to consider adiabatic transformations, and allows the whole system of equations to be hyperbolic. However, this particular case has no consequences on the main purpose of this work.

First order system of conservation laws

Gathering Equations ( 2), ( 6) and ( 7) allows to form a system of conservation laws

∂U ∂t + DIV F = 0, (9)
where U and F denote the vector of conserved quantities and the flux vector respectively, defined as

U =          F p E          ; F =          -v ⊗ 1 -P -P T • v          . ( 10
)
System ( 9) also read in castesian coordinates as

∂U ∂t + 3 α=1 ∂F α ∂X α = 0, (11)
with the flux vector F α = F • E α in the α-th material direction E α , defined as

F α =          -v ⊗ E α -P • E α -(P T • v) • E α          . ( 12
)
It is important to note that System (9) reduces to the well-known Rankine-Hugoniot jump conditions across any discontinuity of fields

S [U] = [F] • N, (13)
where the vectors U and F are defined by Equation [START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part i: Total lagrangian isothermal elasticity[END_REF], N is the material normal of the discontinuity surface moving at speed S , and [•] denotes the jump of the quantity (•) across the discontinuity, such that

[•] = (•) + -(•) -.
Clearly, the discretization of System (9) with any conservative numerical scheme [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] will ensure that the right shock speeds will be computed, since conditions (13) will be correctly approximated.

Regarding the flux vector [START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part i: Total lagrangian isothermal elasticity[END_REF], it is convenient to introduce the auxiliary vector

Q =          P v η          (14)
of which the flux vector F(Q) is a function, and η denotes here the entropy density. The closure of the above system of conservation laws ( 9) is performed by means of a set of constitutive equations, relating the vector of conserved quantities U (10) to the auxiliary vector Q (14). Especially, stresses P involved in the auxiliary vector Q are functions of both the deformation gradient F and the internal energy density E

P = G(F, E(F, η, Z)) (15)
where G denotes some functional defining the constitutive response, and Z denote a set of internal variables describing dissipative phenomena. The internal energy density E appearing as an argument of the functional G in Equation ( 15) is deduced from Equation ( 8), provided the total energy density and the linear momentum density are updated from their respective conservation laws [START_REF] Beni | Consistent arbitrary lagrangian eulerian formulation for large deformation thermo-mechanical analysis[END_REF] and [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]. Besides, Systems [START_REF] Bleier | Efficient variational constitutive updates by means of a novel parameterization of the flow rule[END_REF] or [START_REF] Bonet | A computational framework for polyconvex large strain elasticity[END_REF] must also be supplemented with appropriate initial and boundary conditions.

2.2.

Constitutive laws for thermomechanical dissipative solids 2.2.1. General framework Following the class of Generalized Standard Materials [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF], we assume the existence of a Helmholtz free energy density potential W(F, T, Z) which is a function of the deformation gradient F, the temperature T , and some internal state variables Z describing dissipative phenomena. Following [START_REF] Germain | Continuum Thermodynamics[END_REF][START_REF] Lemaitre | Mechanics of Solid Materials[END_REF], Helmholtz's free energy W is a convex function of F and Z, concave with respect to T , and contains the zero. It is also assumed the existence of a dissipation pseudo-potential φ( Ḟ, Ż; F, Z, T ), which is a convex function of its arguments, positive, and vanishes at zero.

Thermodynamic forces consist of the first Piola-Kirchhoff stress tensor P conjugate to the deformation gradient F, and of forces Y conjugate to internal variables Z. These thermodynamic forces can be additively decomposed into reversible and irreversible components [START_REF] Ziegler | An introduction to thermomechanics[END_REF]:

P = P rev + P irr Y = Y rev + Y irr ( 16
)
such that their reversible components are conjugate to state variables (F, Z) through Helmholtz's free energy, and their irreversible components are conjuguate to the rate of these state variables through the dissipation pseudo-potential:

P rev = ∂W ∂F , P irr = ∂φ ∂ Ḟ Y rev = ∂W ∂Z , Y irr = ∂φ ∂ Ż (17)
Since internal variables should not produce any work, i.e. Y • Ż = 0, ∀ Ż, it follows that

Y rev + Y irr = 0. ( 18
)
Following [START_REF] Canadija | On the thermomechanical coupling in finite strain plasticity theory with non-linear kinematic hardening by means of incremental energy minimization[END_REF], evolutions laws of internal variables are obtained as

∂W ∂Z + ∂Φ ∂ Ż = 0, (19)
and from [START_REF] Busto | High order ader schemes for continuum mechanics[END_REF], the first Piola-Kirchchoff stress tensor P reads:

P = P rev + P irr = ∂W ∂F + ∂Φ ∂ Ḟ . (20)
Next, the local thermal equilibrium is assumed to be always verified. The consequence is that the entropy density η is given by its sole reversible component:

η = - ∂W ∂T , (21)
and is conjuguate to the temperature T through Helmholtz's free energy. Finally, the second law of thermodynamics is introduced, that reads in the absence of heat conduction as

D int = P : Ḟ -(η Ṫ + Ẇ) ≥ 0, (22)
where D int denotes the mechanical dissipation that should be non-negative, which is caused by the dissipative forces P irr and Y irr , as shown when combining with Equations ( 16) and [START_REF] Camacho | Adaptive lagrangian modelling of ballistic penetration of metallic targets[END_REF]:

D int = P irr : Ḟ + Y irr • Ż ≥ 0. ( 23
)

Variational formulation of the thermomechanical local constitutive problem

Following the works of Ortiz, Stainier and co-workers [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF], we start with the following functional

D( Ḟ, η, Ż, T ; F, η, Z) = Ė(F, η, Z) -T η + φ( Ḟ, Ż; F, Z, T ) (24)
where Ė, η and φ are the rate of internal energy density, the rate of entropy density and the dissipation pseudo-potential respectively. Since the reversible power per unit volume received by the system ẇτ can be defined as

ẇτ = Ė(F, η, Z) -T η = P rev : Ḟ + Y rev : Ż, (25)
the functional [START_REF] Fancello | A variational constitutive update algorithm for a set of isotropic hyperelastic-viscoplastic material models[END_REF] can thus be interpreted as the sum of the two homogeneous terms representing a power per unit volume:

D( Ḟ, η, Ż, T ; F, η, Z) = ẇτ + φ. ( 26
)
The first one is associated with reversible processes, the other is associated with irreversible processes [START_REF] Pethe | Remapping-free variational h-adaption for strongly coupled thermo-mechanical problems[END_REF]. Another interpretation of the functional D (24) proposed by [START_REF] Canadija | On the thermomechanical coupling in finite strain plasticity theory with non-linear kinematic hardening by means of incremental energy minimization[END_REF] in the case of finite elastoplasticity is that it corresponds to the stress power P : Ḟ. One way among others to combine the above functional with a description of the reversible part of the constitutive response through Helmholtz's free energy W(F, T, Z) is to introduce the residual of the rate of the associated Legendre transform, enforced to vanish through a Lagrange multiplier. We are then led to introduce the following Lagrangian functional L( q, λ;

q) = Ė -T η + φ( Ḟ, Ż; F, Z, T ) + λ d dt (T η + W(F, T, Z) -E) (27)
where the following state vector q = {E, F, η, Z, T } has been introduced, which is assumed to be known and fixed here, and λ denotes the Lagrange multiplier. Moreover, it is assumed that the values of both the rate of deformation gradient Ḟ and the rate of internal energy density Ė are given and known from the set of conservation laws [START_REF] Bleier | Efficient variational constitutive updates by means of a novel parameterization of the flow rule[END_REF]. The optimization problem defined with this Lagrangian functional thus reads

W = stat η, Ṫ ,λ inf Ż L( q, λ; q) (28)
where the variable with respect to which the stationarity conditions are computed are the rate of entropy density η, the temperature rate Ṫ , the Lagrange multiplier λ, and the rate of internal variables Ż. Stationarity conditions of this Lagrangian functional yields

stat η L ⇔ -T + λT = 0 (29) stat Ṫ L ⇔ λ η + ∂W ∂T = 0 (30) stat λ L ⇔ d dt (T η + W(F, T, Z) -E) = 0 (31) inf Ż L ⇔ ∂φ ∂ Ż + λ ∂W ∂Z = 0 ( 32
)
The stationarity with respect to η [START_REF] Ghavamian | An entropy-stable smooth particle hydrodynamics algorithm for large strain thermo-elasticity[END_REF] gives the expected result of the Lagrange multiplier which is equal to unity λ = 1. Substitution of Equation [START_REF] Ghavamian | An entropy-stable smooth particle hydrodynamics algorithm for large strain thermo-elasticity[END_REF] into Equations ( 30) and [START_REF] Grüneisen | Theorie des festen zustandes einatomiger elemente[END_REF], namely the stationarities with respect to the temperature rate and the rate of internal variables, allows to recover the definition of the entropy [START_REF] Collela | Multidimensional upwind methods for hyperbolic conservation laws[END_REF] and that of the evolution equations of internal variables [START_REF] Canadija | A variational formulation for thermomechanically coupled low cycle fatigue at finite strains[END_REF] respectively. The optimization problem [START_REF] Germain | Continuum Thermodynamics[END_REF] is thus formulated as a rate-type problem, because only rate quantities appears in the optimization variables (apart the Lagrange multiplier), and the state vector q is considered fixed.

Remark 2. Substitution of the solved Lagrange multiplier (λ = 1) into the Lagrangian functional [START_REF] Gavrilyuk | An example of a one-parameter family of rank-one convex stored energies for isotropic compressible solids[END_REF] allows to recover the functional initially considered in the works of Ortiz, Stainier and co-workers [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF]]

L(λ = 1) = Ẇ + Ṫ η + φ. (33)
The reason why this Lagrange multiplier is kept as an unknown of the optimization problem (28) will appear more explicitly when deriving the incremental variational constitutive update in Section 4.

Remark 3. The stationarity equation with respect to the temperature rate (30) allows to recover the definition of the entropy. This is made possible because the dissipation pseudo-potential φ( Ḟ, Ż; F, Z, T ) is not a function of Ṫ , but only of the temperature T .

Remark 4. Previous variational principles of the coupled thermomechanical problem [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF][START_REF] Canadija | On the thermomechanical coupling in finite strain plasticity theory with non-linear kinematic hardening by means of incremental energy minimization[END_REF] relied on a twofield temperature formulation, an external temperature arising in the heat equation plus some internal temperature conjuguate to the entropy density through the thermodynamic potential. Equality of both was imposed as an internal constraint weakly. Observe now that such two-field temperature formulation is not required anymore to recover the sought state and evolution constitutive equations, thus it becomes useless.

Finally, the first Piola-Kirchhoff stresses are computed taking the partial derivative of the optimized functional with respect to the rate of the deformation gradient Ḟ:

∂W ∂ Ḟ = λ ∂W ∂F + ∂φ ∂ Ḟ = P rev + P irr = P (34)
which is identical to [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF].

A first order accurate discrete variational constitutive update

From the previous continuous variational principle, we are now interested in deriving an incremental variational constitutive update. Consider a discrete time increment [t n , t n+1 ], the material state vector q n = {E n , F n , η n , Z n , T n } known at time t n , and data {E, F} n+1 known and updated at time t n+1 through the solution of the discrete conservation laws. Then, an incremental functional I(q n+1 , λ n+1 ; q n ) is sought in such a way that it approximates the integral of the Lagrangian functional L (27) over the time increment ∆t:

I(q n+1 , λ n+1 ; q n ) ≈ t n+1 t n L( q(τ), λ(τ); q(τ))dτ = ∆E -T n ∆η + ∆tφ ∆F ∆t , ∆Z ∆t ; F n+α , Z n+α , T n+α + λ n+1 ∆(T η + W(F, T, Z) -E) (35)
where the operator ∆(•) = (•) n+1 -(•) n denotes the finite difference between the values of the quantity (•) at times t n+1 and t n . Parameter arguments of the dissipation pseudo-potential after the semi-colon are computed at time

t n+α = t n + α∆t such that (•) n+α = (1 -α)(•) n + α(•) n+1 , α ∈ [0, 1]
. Indeed, the most part of the numerical error performed in the approximation [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF] may come from the integral of the pseudo-potential φ. It has been shown in [START_REF] Brassart | On convergence properties of variational constitutive updates for elasto-visco-plasticity[END_REF] that the evaluation at some intermediate time of these parameter arguments may reduce this numerical error, the optimal value of α will thus depend on the chosen material model and parameters. Due to that, but also regarding the temperature evaluated at time t n in the second term of (35), the numerical scheme here derived thus appears to be semi-implicit. The incremental variational update takes thus the following form

W n+1 = stat (η,T,λ) n+1 inf Z n+1 I(q n+1 , λ n+1 ; q n ) ( 36
)
where the incremental functional [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF] is optimized with respect to the entropy density, the temperature, the Lagrange multiplier and the set of internal variables, all evaluated at time t n+1 . The stationarity conditions in the discrete case thus read: stat

η n+1 I ⇔ -T n + λ n+1 T n+1 = 0 ( 37
)
stat

T n+1 I ⇔ ∆tα ∂φ ∂T ∆F ∆t , ∆Z ∆t ; F n+α , Z n+α , T n+α + λ n+1 η n+1 + ∂W ∂T n+1 (38)
stat

λ n+1 I ⇔ ∆(T η + W(F, T, Z) -E) = 0 (39)
inf

Z n+1 I ⇔ ∂φ ∂ Ż ∆F ∆t , ∆Z ∆t ; F n+α , Z n+α , T n+α + λ n+1 ∂W ∂Z n+1 = 0 ( 40
)
The stationarity with respect to the entropy density [START_REF] Heuzé | Lax-wendroff schemes for elastic-plastic solids[END_REF] gives the updated value of the Lagrange multiplier

λ n+1 = T n T n+1 , (41)
which is not equal to unity anymore, but is close to it, especially since the absolute temperatures are considered here.

Introducing the Lagrange multiplier (41) into stationarity equations ( 38) and ( 40) yields the discrete update of the entropy density and of the discrete evolution equations of internal variables:

η n+1 = - ∂W ∂T n+1 - T n+1 T n ∆tα ∂φ ∂T ∆F ∆t , ∆Z ∆t ; F n+α , Z n+α , T n+α O(∆t) (42)
∂φ ∂ Ż ∆F ∆t , ∆Z ∆t ; F n+α , Z n+α , T n+α + T n T n+1 ∂W ∂Z n+1 = 0 (43)
Observe that the definition of the entropy density [START_REF] Collela | Multidimensional upwind methods for hyperbolic conservation laws[END_REF] is obtained in Equation (42) up to a term which is of the order of ∆t, which is the signature of the first order accuracy of this variational integrator. This is an approximation, some sort of numerical relaxation of the thermal equilibrium, but that tends to zero as ∆t goes to zero, hence ensuring the consistency of the variational integrator. Next, the discrete evolution equations [START_REF] Lee | Elastic-plastic deformation at finite strains[END_REF] are identical to these obtained in [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF]. As the time step ∆t goes to zero, this expression also tends to its continuous counterpart [START_REF] Grüneisen | Theorie des festen zustandes einatomiger elemente[END_REF].

The first Piola-Kirchhoff stresses are then computed taking the partial derivative of the optimized functional (35) with respect to the updated deformation gradient F n+1 :

∂W ∂F n+1 = T n T n+1 ∂W ∂F n+1 + ∂φ ∂ Ḟ ∆F ∆t , ∆Z ∆t ; F n+α , Z n+α , T n+α (44)
= P rev n+1 + P irr n+1 = P n+1 (45)
where the Lagrange multiplier (41) weights the partial derivative of the free energy with respect to the deformation gradient in the definition of reversible stresses in the discrete case, which is another consequence of the chosen approximation [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF] and of results of the continuous principle [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF].

Remark 5. If the Lagrange multiplier λ n+1 would have been set to unity, then the incremental functional [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF] would have reduced to that of [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF], up to a term of second order. However, this is not the case, the variational integrator derived is therefore original, and different from that of these works. More precisely, the incremental functional [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF] expressed with the optimal value of the Lagrange multiplier reads

I λ = T n T n+1 = T n T n+1 ∆W + ∆E ∆T T n+1 + η n ∆T T n T n+1 + ∆tφ. ( 46
)
Remark 6. A ratio of updated and previous values of the temperature appears here through the expression of the Lagrange multiplier [START_REF] Kumar | Eulerian thermo-mechanical simulations of heterogeneous solid propellants using an approximate projection method[END_REF], which is reminiscent of the one already appearing in the discrete integration factor of the discrete incremental potential introduced in [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF]. However, this ratio serves different purposes in this work and in that of these authors.

5. Application to thermo-hyperelastic-viscoplastic solid media

General form of Helmholtz's free energy

We are interested in this work in the particular case of thermo-hyperelastic-viscoplastic solid media. After the pioneering work of Lee [START_REF] Lee | Elastic-plastic deformation at finite strains[END_REF], the kinematics of such media follows a multiplicative decomposition of the deformation gradient into an elastic part F e and a viscoplastic part F p :

F = F e • F p , det F e > 0, det F p > 0. ( 47
)
Next, from the elastic part of the deformation gradient, the elastic right Cauchy-Green strain tensor

C e = (F e ) T • F e (48)
is generally used to compute stresses, in order to satisfy to material frame indifference. The set of internal variables Z now consists of the viscoplastic part of the deformation gradient F p , some strain-like second-order tensor α k associated with kinematic hardening, and a scalar strain-like variable α i associated with isotropic hardening, hence

Z = {F p , α k , α i }.
From this, Helmholtz's free energy is then defined by summing several contributions

W(F, T, Z) = W e (C e , T ) + W p (F p , α k , α i , T ) + W th (T ) (49)
where W e (C e , T ) is the elastically stored energy (recoverable), W p (F p , α k , α i , T ) is the plastically stored energy (not recoverable), and W th (T ) is the thermally stored energy due to heat capacity. State laws associated with strain-like hardening variables read as

Q k = - ∂W ∂α k ; Q i = - ∂W ∂α i , (50)
where Q k and Q i are the thermodynamic forces associated with the kinematic and isotropic hardenings respectively. The thermally stored energy W th (T ) is classically chosen after [START_REF] Simo | Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation[END_REF] as

W th (T ) = ρ 0 C 0 (T -T 0 ) -T ln T T 0 , (51)
where ρ 0 C 0 is some reference volumic heat capacity, and T 0 refer to some reference temperature. Moreover, in order to account for different material response of volumetric and isochoric parts, the isochoric elastic deformation gradient and right Cauchy-Green strain tensors are introduced, such that

Fe = (J e ) -1/3 F e , J e = det F e , det Fe = 1 (52) Ce = ( Fe ) T • Fe (53)
The elastically stored energy is then decomposed additively into volumetric and isochoric parts

W e (C e , T ) = W H (J e , T ) + We ( Ce , T ) (54)
where the volumetric part W H (J e , T ) depends on the elastic jacobian determinant J e and on the temperature T , while the isochoric component We ( Ce , T ) is a function of both the isochoric elastic Cauchy-Green strain tensor Ce and the temperature T .

Remark 7. In the sequel, it will be assumed for convenience that the plastic flow is isovolume, namely

J = J e , J p = det F p = 1, (55)
hence W H (J e , T ) ≡ W H (J, T ).

Hyperelasticity

The isochoric component We ( Ce , T ) describes hyperelastic distortions. A family of rank-one convex stored energies was proposed in [START_REF] Gavrilyuk | An example of a one-parameter family of rank-one convex stored energies for isotropic compressible solids[END_REF], and reused in [START_REF] Boscheri | A cell-centered lagrangian ader-mood finite volume scheme on unstructured meshes for a class of hyper-elasticity models[END_REF], which reads

We ( Ī1 , Ī2 , T ) = µ(T ) 4 -2a( Ī1 -3) + (1 + a) 3 ( Ī2 -9) ( 56
)
where a is an adjustable parameter, and invariants Ī1 , Ī2 are defined such as

Ī1 = tr Ce , Ī2 = tr [( Ce ) 2 ]. ( 57
)
The shear modulus µ(T ) may also depend on temperature, and is often described in terms of the Young modulus E and Poisson's ratio ν as

µ = E 2(1 + ν) . ( 58
)
In the range a ∈ -1, 1 2 , the authors of [START_REF] Gavrilyuk | An example of a one-parameter family of rank-one convex stored energies for isotropic compressible solids[END_REF] showed that the resulting first order system of equations is hyperbolic. In the particular case where a = -1, a more classical neo-Hookean material is obtained, whose distortional free energy reads

We ( Ī1 ,

T ) = µ(T ) 2 ( Ī1 -3). ( 59
)
Notice that polyconvex hyperelastic models (e.g. Mooney-Rivlin) can also be used. The isochoric elastic free energy We then becomes a function of the three minors of the elastic deformation, the system of conservation laws ( 9) is then extended to the other deformation minors, Equations ( 4) and ( 5), and the resulting first order system of equations is also shown to be hyperbolic [START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part iii: Thermo-elasticity[END_REF].

Equation of state

The volumetric component W H (J, T ) allows to introduce an equation of state, which we shall always choose as being convex to ensure hyperbolicity of the volumetric part of the first order system of equations. Moreover, this will generate shock waves in the simulation. A well-known equation of state is that of Mie-Gruneïsen [START_REF] Mie | Zur kinetischen theorie der einatomigen körper[END_REF][START_REF] Grüneisen | Theorie des festen zustandes einatomiger elemente[END_REF], establishing a nonlinear relationship between the pressure and the internal energy which is only a function of J:

J d p dE J=constant = -Γ 0 J q (60)
where Γ 0 is the Grüneisen constant. Once this expression integrated, and following [START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part iii: Thermo-elasticity[END_REF], the volumetric free energy

W H (J, T ) = U(J) -T η R (J) (61)
is expressed as a function of the volumetric part of the internal energy U(J) introducing thermoelastic effects

U(J) = κ 2 (J -1) 2 + ρ 0 CΓ 0 T 0 (J -1), ( 62
)
where κ is the bulk modulus, and as a function of some reference entropy depending on J written after integration of (60) as

-η R (J) = ρ 0 CΓ 0 (1 -J q ) q (63) if q ∈]0, 1],
or that reads as -η R (J) = -ρ 0 CΓ 0 ln J if q = 0, hence obtaining a response of fluid. From (61), the hydrostatic pressure defined as

P = ∂W H ∂J (64)
reads

P = κ(J -1) -ρ 0 CΓ 0 (T J q-1 -T 0 ), (65)
from which the effective bulk modulus including thermoelastic effects reads

κ = d p dJ J=1 = κ + ρ 0 CT 0 Γ 0 (1 -q). (66)
Another simple expression of equation of state is the following

W H (J, T ) = κ 2 (ln J) 2 -ρ 0 CΓ 0 (T -T 0 ) ln J, (67)
where the logarithm of the jacobian has appeared, leading to the following expression of the hydrostatic pressure

P = κ ln J J - ρ 0 CΓ 0 (T -T 0 ) J . ( 68
)
Such equation of state will be used in Section 7.3 for comparison purpose with finite element temperature-based numerical solutions in the large strain framework.

Continuous variational formulation of the constitutive model

Kinematical approach

The viscoplastic velocity gradient

L p = Ḟp • (F p ) -1
, defined in the intermediate configuration, is generally chosen to write the viscoplastic flow rule

L p = ṗM, ( 69
)
where M is the viscoplastic flow direction, and ṗ denotes the effective viscoplastic strain rate or flow intensity. From Equations ( 69), the set of internal variables Z = {F p , α k , α i } considered for viscoplastic solids can be reduced to consider only the following rate of internal variables Ż = { ṗ, M}, that is the effective viscoplastic strain rate and the flow direction. Indeed, the latter enters the flow rule defining F p by Equation ( 69), but is also part of that associated with the kinematic strain variable α k , namely H k (M). Such approach could be called a kinematical one, and was followed in [START_REF] Fancello | A variational constitutive update algorithm for a set of isotropic hyperelastic-viscoplastic material models[END_REF][START_REF] Brassart | On convergence properties of variational constitutive updates for elasto-visco-plasticity[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF].

The stationarity conditions ( 32) is thus written with these two variables. The first one reads

∂L ∂ ṗ = ∂W ∂p + ∂φ ∂ ṗ = 0. ( 70
)
Since the deformation mapping is assumed to be constant during the constitutive update Ḟ = 0, the rate of the elastic part of the deformation gradient (47

) reads Ḟe = -F e • L p , (71)
hence the differential of the elastic part of the free energy can be expressed (at fixed temperature T ) as

Ẇe Ṫ =0 = ∂W e ∂F e : Ḟe = -(F e ) T • ∂W e ∂F e : L p = -Σ : L p , (72)
where Σ = (F e ) T • ∂W e ∂F e is the Mandel stress tensor. Combining Equations ( 69) and ( 72) yields the contribution of the elastic part of the free energy to the first term of ( 70)

∂W e ∂p = -Σ : M. ( 73
)
Considering state equations [START_REF] Marusich | Modelling and simulation of high-speed machining[END_REF] for the contribution associated with the plastic part of the free energy, the stationarity condition (70) finally simplifies as

∂L ∂ ṗ = -Σ : M + Q i -Q k : H k (M) + ∂φ ∂ ṗ = 0. ( 74
)
The stationarity condition written with respect to the flow direction M then reads

∂L ∂M = ∂ Ẇ ∂M = 0. ( 75
)
However, the above minimization has to comply with physical constraints resulting on the one hand from the flow rule (and det F p > 0), and on the other hand from evolution equations relating stresses to the internal variables α k and α i . For instance, for a von Mises-type yield criterion, the flow direction must satisfy the following constraints

M = M T , tr[M] = 0, M : M = 3 2 , ( 76
)
which should be taken into account as equality-type constraints in the optimization problem inf M L. This can be done using a Lagrange multiplier approach, as it was shown analytically in the isothermal setting in [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF]. More precisely, an analytical treatment of this minimization was only found to be possible when identifying an additive structure between the updated elastic strain and its elastic predictor, which is the case for the small strain framework, and can be extended to the large strain one using Hencky hyperelasticity which relies on the logarithmic strain. Such approach is restricted to associative plasticity, i.e. H k = -M, the resulting flow rule direction reads in this case

M = 3 2 dev[Σ -Q k ] dev[Σ -Q k ] . (77)
The variational integrator was then identified to the well-known radial return algorithm [START_REF] Wilkins | Methods in computational physics[END_REF]. To account for more general hyperelastic laws, a spectral decomposition of the flow direction M in its principal values and directions was proposed in [START_REF] Fancello | A variational constitutive update algorithm for a set of isotropic hyperelastic-viscoplastic material models[END_REF]. The minimization is then performed numerically, then an integrator different from the radial return is obtained.

Parameterization of the flow rule direction based on pseudo-stresses

In the work of Mosler and co-workers [START_REF] Mosler | Towards variational constitutive updates for non-associative plasticity models at finite strain: models based on a volumetric-deviatoric split[END_REF][START_REF] Mosler | On the implementation of rate-independent standard dissipative solids at finite strain-variational constitutive updates[END_REF][START_REF] Bleier | Efficient variational constitutive updates by means of a novel parameterization of the flow rule[END_REF], another parameterization of the flow rule direction has been introduced. Pseudo-stresses Σ are introduced, which are a priori different from their physical counterparts, that is Σ (Σ -Q k ), and serve for the parameterization of the flow direction, namely

M = M( Σ). ( 78
)
Such parameterization is compatible with the equivalent stress measure Σ eq , defining the shape of the yield function, which are positively homogeneous functions of degree one, satisfying

Σ eq (cA) = cΣ eq (A), ∀A, ∀c ∈ R + . (79)
A wide range of yield criteria can then be accounted for, especially anisotropic ones. In order to define the flow direction [START_REF] Trangenstein | Numerical solution of hyperbolic partial differential equations[END_REF], it is convenient to introduce some flow function

g(Σ, Q k , Q i , T ), which is a function of thermodynamic forces Y ≡ Y irr = {Σ, Q k , Q i }
and whose purpose is to define the flow direction of the strain rate variables Ż = {L p , αk , αi }. More precisely, in the absence of irreversible stresses P irr , a dual dissipation pseudo-potential of φ( Ż) can be defined through the following Legendre transform

φ * (Y) = sup Ż (Y • Ż -φ( Ż)), (80)
such that the associated evolution equations can be written as

Ż = ∂φ * ∂Y . (81)
The dual dissipation pseudo-potential φ * (Y) is a convex function of its arguments, and defines a family of equipotential surfaces on which any point yields the same dissipation and effective viscoplastic strain rate. This potential may depend on flux variables Y through the flow function g(Y, T ), so that the evolution laws read:

Ż = ∂φ * ∂g ∂g ∂Y = ṗ ∂g ∂Y ( 82
)
where ṗ = ∂φ * ∂g ( 83
)
denotes the effective viscoplastic strain rate and ∂g/∂Y is the flow direction, normal to the flow function g. It is usual to consider such flow function g of the form [START_REF] Canadija | On the thermomechanical coupling in finite strain plasticity theory with non-linear kinematic hardening by means of incremental energy minimization[END_REF]:

g(Σ, Q k , Q i , T ) = f (Σ, Q k , Q i , T ) + f (Q k ) ( 84
)
to account for non-associativity of the kinematic hardening with respect to the yield function

f (Σ, Q k , Q i , T ), of the form f (Σ, Q k , Q i , T ) = Σ eq (Σ -Q k ) -Q i (α i , T ) -σ y (T ), (85)
where σ y is the tensile yield stress. The yield function allows to define the space of elastic stresses:

C Σ = {(Σ, Q k , Q i , T )| f (Σ, Q k , Q i , T ) ≤ 0} (86)
Adopting pseudo-stresses to parameterize the flow direction [START_REF] Trangenstein | Numerical solution of hyperbolic partial differential equations[END_REF], and accounting for the property (79), the dual evolutions laws (82) are then written as

L p = ṗ ∂g ∂Σ Σ = ṗ ∂ f ∂Σ Σ (87) αk = ṗ ∂g ∂Q k = -L p + ṗ ∂ f ∂Q k (88) αi = ṗ ∂g ∂Q i = -ṗ (89)
where the derivative ∂ f ∂Σ Σ means it is evaluated with Σ = Σ. Accordingly, the stationarity equation ( 74) is rewritten as

∂L ∂ ṗ = -(Σ -Q k ) : ∂ f ∂Σ Σ + Q i -Q k : ∂ f ∂Q k + ∂φ ∂ ṗ = 0, (90)
and simplifies accounting for the formula Σ eq = Σ : ∂ f ∂Σ [START_REF] Mosler | Towards variational constitutive updates for non-associative plasticity models at finite strain: models based on a volumetric-deviatoric split[END_REF], resulting from the positive homogeneity of the equivalent stress Σ eq , as

-f -σ y -Q k : ∂ f ∂Q k + ∂φ ∂ ṗ = 0. ( 91
)
For rate-independent elastoplasticity, the dissipation pseudo-potential φ identifies to the mechanical dissipation

D int = ṗ σ y + Q k : ∂ f ∂Q k ≥ 0 (92)
hence Equation (91) readsf = 0, and ensures that the stress state lies on the boundary of the elastic convex during plastic loading. It should then be added the consistency condition ḟ = 0, to determine the cumulated plastic strain p (or plastic multiplier). In rate-dependent plasticity, or viscoplasticity, the expression of the dissipation pseudo-potential φ allows the existence of viscous or overstresses, from which the viscoplastic flow follows.

The stationarity condition (90) should be supplemented with another one written on the pseudo-stresses Σ, which in a sense, replaces Equation [START_REF] Taylor | The latent energy remaining in a metal after cold working[END_REF], such that

∂L ∂ Σ = ∂ Ẇ ∂ Σ = -ṗ(Σ -Q k ) : ∂ 2 f ∂Σ∂Σ Σ = 0. ( 93
)
Equation ( 75) enforces the correct flow direction, which is compatible with the stresses. Hence, with such parameterization, non-associative evolution equations are a priori prescribed, and it results in an unconstrained optimization problem. However, the counterpart of this approach is that it requires to introduce the dual dissipation pseudo-potential φ * , and therefore the flow function g, in addition to the primal dissipation pseudo-potential φ. The former allows on the one hand to account for constraints put on the flow direction for a certain range of yield functions by introducing pseudo-stresses, and on the other hand ensure its compatibility with evolution equations. The latter permits to actually introduce viscous (or viscoplastic) effects in the computation of the viscoplastic strain rate ṗ.

To complete the solution process, the stationarity condition (31) associated with the thermal part of the problem should be added to Equations ( 90) and (93).

Discrete variational formulation of the constitutive model

Following the parameterization of Mosler and co-workers [START_REF] Mosler | Towards variational constitutive updates for non-associative plasticity models at finite strain: models based on a volumetric-deviatoric split[END_REF][START_REF] Mosler | On the implementation of rate-independent standard dissipative solids at finite strain-variational constitutive updates[END_REF][START_REF] Bleier | Efficient variational constitutive updates by means of a novel parameterization of the flow rule[END_REF] with pseudo-stresses, the incremental functional I (35) written with the expression of the Lagrange multiplier (41) can be expressed as

I(T, ∆p, Σ) = ∆E -T n ∆η + ∆tφ ∆p ∆t , p n+α , T n+α + T n T n+1 ∆(T η + W(C e , p, T, Σ) -E) (94)
where the increment of internal energy density ∆E, and the updated deformation gradient F n+1 are assumed fixed and known from the solution of discrete conservation laws. The unknown vector at time t n+1 is defined as

q n+1 = {T, ∆p, Σ} n+1 . ( 95
)
For the sake of simplicity, the functional ( 35) is here simplified for the sole isotropic hardening, the reader is referred to [START_REF] Mosler | Variationally consistent modeling of finite strain plasticity theory with non-linear kinematic hardening[END_REF] for the details about kinematic hardening. Taking the exponential update formula from [START_REF] Weber | Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids[END_REF] for the viscoplastic part of the deformation

F p n+1 = exp ∆p ∂ f ∂Σ Σn+1 • F p n ( 96
)
permits to preserve the isochoric nature of the viscoplastic deformation. From Equation (96), the updated elastic part of the deformation gradient reads

F e n+1 = F e trial • exp -∆p ∂ f ∂Σ Σn+1 , F e trial = F n+1 • (F p n ) -1 . ( 97
)
After Equation ( 40), the discrete stationarity equation of the incremental functional (94), with respect to the increment of cumulated viscoplastic strain ∆p is written as

∂I ∂∆p = T n T n+1 ∂W ∂∆p n+1 + ∂φ ∂ ṗ ∆p ∆t , p n+α , T n+α + α∆t ∂φ ∂p ∆p ∆t , p n+α , T n+α = 0. ( 98
)
The partial derivative of the free energy with respect to ∆p reads

∂W ∂∆p n+1 = ∂W e ∂∆p n+1 + ∂W p ∂∆p n+1 (99
)

∂W e ∂∆p n+1 = -2 C e trial • exp -∆p ∂ f ∂Σ Σn+1 • ∂W e ∂C e n+1 : D exp -∆p ∂ f ∂Σ Σn+1 : ∂ f ∂Σ Σn+1 ( 100
)
∂W p ∂∆p n+1 = Q i ( 101
)
where D exp(A) stands for the derivative of the exponential mapping of the matrix A, computed with standard procedures [START_REF] Ortiz | The computation of the exponential and logarithmic mappings and their first and second linearizations[END_REF]. The partial derivative ∂W e ∂C e is given after the splitting (54) between isochoric and volumetric contributions as

∂W e ∂C e = ∂ We ∂C e + P J 2 (C e ) -1 , (102)
where the hydrostatic pressure P is defined by Equation ( 64), and may follow the expressions given by the Mie-Grüneisen equation of state [START_REF] Renaud | On loading paths followed inside plastic simple waves in two-dimensional elastic-plastic solids[END_REF] or by Equation ( 68).

The discrete stationarity equation of (94) with respect to the pseudo-stresses Σn+1 gives

∂I ∂ Σn+1 = ∂W e ∂ Σn+1 = -2∆p C e trial • exp -∆p ∂ f ∂Σ Σn+1 • ∂W e ∂C e n+1 : D exp -∆p ∂ f ∂Σ Σn+1 : ∂ 2 f ∂Σ∂Σ Σn+1 = 0. ( 103
)
The consistency of the discrete equations ( 98) and ( 103) with respect to their continuous counterparts (90) and ( 93) are shown in [START_REF] Mosler | On the implementation of rate-independent standard dissipative solids at finite strain-variational constitutive updates[END_REF][START_REF] Bleier | Efficient variational constitutive updates by means of a novel parameterization of the flow rule[END_REF], and are not repeated here. Again, the discrete scalar equation ( 39) should be added to Equations ( 98) and ( 103) to obtain a balanced system of equations, solved on the unknown vector q n+1 (95). If a Newton method is used, the following linear system of equations should be solved at each iteration k

K (k) δq (k) = -R (k) , (104)
with δq (k) = q (k+1)q (k) , and with the Hessian matrix and the residual vector expressed as

K (k) =              ∂ 2 I ∂T 2 ∂ 2 I ∂T ∂∆p ∂ 2 I ∂T ∂ Σ ∂ 2 I ∂∆p∂T ∂ 2 I ∂∆p 2 ∂ 2 I ∂∆p∂ Σ ∂ 2 I ∂ Σ∂T ∂ 2 I ∂ Σ∂∆p ∂ 2 I ∂ Σ∂ Σ              , R (k) =            ∂I ∂T ∂I ∂∆p ∂I ∂ Σ            (k) , (105)
the Hessian matrix being symmetric. The entries of the tangent matrix K (k) associated with unknowns {∆p, Σ} are already detailed in [START_REF] Mosler | On the implementation of rate-independent standard dissipative solids at finite strain-variational constitutive updates[END_REF][START_REF] Bleier | Efficient variational constitutive updates by means of a novel parameterization of the flow rule[END_REF], and are not repeated here, while others involving partial derivative with respect to temperature can be derived in a straightforward manner. Notice also that partitioned schemes can be adopted to solve System (104), such as Gauss-Seidel, staggered or nested approaches [START_REF] Keyes | Multiphysics simulations: Challenges and opportunities[END_REF]. Finally, it is convenient to parameterize the pseudo-stresses Σ in the solution process with spherical coordinates to enforce Σ = 1, as shown in [START_REF] Bleier | Efficient variational constitutive updates by means of a novel parameterization of the flow rule[END_REF], which allows to eliminate any singularity of the Hessian matrix.

Since the dissipation pseudo-potential φ is not a smooth function for ∆p = 0, a prediction-correction scheme is classically followed. A thermoelastic prediction is first performed by solving the discrete Legendre transform [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF], with the updated trial deformation gradient F e trial (97), and ∆p = 0. Next, the slope of the incremental potential ∂I ∂∆p (98) is evaluated for ∆p = 0 + . Since I is convex with respect to ∆p, if the slope ∂I ∂∆p < 0 is negative, the optimal ∆p is positive. Actually, it amounts to check that f trial n+1 > 0. Then, the system consisting of Equations ( 39), (98) and ( 103) is solved, for instance through a Newton Method (104). If ∂I ∂∆p > 0, the optimal ∆p is zero, since negative values of ∆p are prohibited as it would lead to an infinite dissipation.

Once the unknown vector q n+1 has been updated, the plastic part of the deformation gradient is updated with the exponential formula (96), the elastic part of the deformation gradient and the elastic right Cauchy-Green strain tensor are computed with Equations ( 47) and [START_REF] Maire | A nominally second-order cell-centered Lagrangian scheme for simulating elastic-plastic flows on two dimensional unstructured grids[END_REF]. Then, following Equation ( 44), the second Piola-Kirchhoff stress tensor is computed as follows

S n+1 = 2 T n T n+1 ∂W ∂C n+1 = 2 T n T n+1 (F p n+1 ) -1 • ∂W ∂C e n+1 • (F p n+1 ) -T , (106)
while the entropy density η n+1 can be updated with Equation (42).

Fully isotropic constitutive models

If both the elastic response and the yield function are isotropic, on the one hand both the physical Mandel stresses Σ and the pseudo-stresses Σ tensors are symmetric, and on the other hand they share the same eigenbasis with the elastic right Cauchy-Green updated C e and the trial elastic strains C e trial :

Σ = 3 k=1 Σ k B trial k , Σ = 3 k=1 Σk B trial k C e = 3 k=1 λ C e k B trial k , C e trial = 3 k=1 λ C e trial k B trial k , B trial k = N trial k ⊗ N trial k ( 107
)
therefore only two unknowns are required to span all admissible flow directions if a parameterization with spherical coordinates is used [START_REF] Bleier | Efficient variational constitutive updates by means of a novel parameterization of the flow rule[END_REF]. Next, the computation of the exponential mapping of a second order tensor (96) and its derivative involved in Equations ( 100) and (103) can be avoided, which is quite convenient and computationally more efficient.

For a von Mises-type yield function, the flow rule is traceless, hence it requires that the pseudo-stresses also do, thus tr

[ Σ] = 0. ( 108
)
Following the parameterization introduced in [9]

Σk (ψ) = 2 3 sin 2 3 πk -ψ , k = 1, 2, 3 (109)
only one parameter ψ is sufficient to describe all flow directions, then the unknown vector consists now of only three scalar unknowns q n+1 = {T, ∆p, ψ} n+1 , which is numerically very efficient. More precisely, in the stationarity equation of the incremental functional with respect to ∆p (98), the partial derivative of the elastic part of the free energy (100) is now simplified as

∂W e ∂∆p n+1 = -2 3 k=1       λ C e trial k ∂W e ∂λ C e k exp -2∆p ∂ f ∂Σ k Σk (ψ) ∂ f ∂Σ k Σk (ψ)       n+1 , (110)
whereas the stationarity equation of the incremental functional with respect to pseudo-stresses ( 103) is now replaced by the following scalar one

∂I ∂ψ n+1 = ∂W e ∂ψ n+1 = -2∆p 3 k=1 3 l=1        λ C e trial k ∂W e ∂λ C e k exp -2∆p ∂ f ∂Σ k Σk (ψ) ∂ 2 f ∂Σ k ∂Σ l Σk (ψ) ∂ Σl ∂ψ        n+1 = 0.
(111)

Simplifications for the linearized geometrical framework

First, Helmholtz's free energy [START_REF] Malkus | Mixed finite element methods-reduced and selective integration techniques: a unification of concepts[END_REF] can be simplified in the linearized geometrical framework, and accounting for a sole isotropic hardening, as follows

W(ε e , p, T ) = W e (ε e , T ) + W p (p, T ) + W th (T ) (112) W e (ε e , T ) = µ(T )dev[ε e ] : dev[ε e ] + κ 2 (tr ε e ) 2 -ρCΓ 0 (T -T 0 )(tr ε e ) (113) W th (T ) = - ρC 2T 0 (T -T 0 ) 2 -η 0 (T -T 0 ) (114)
where ε e denotes the elastic strain, and the thermally stored energy W th (T ) [START_REF] Mie | Zur kinetischen theorie der einatomigen körper[END_REF] has been linearized to obtain the quadratic form (114). The Grüneisen coefficient Γ 0 = 3κα ρC [START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part iii: Thermo-elasticity[END_REF] can be linked to the dilatation coefficient α. Next, as already mentioned in Section 5.4.1, for linear isotropic elasticity in the small strains framework, the optimization problem inf M L subject to constraints (76) yields the von Mises-type flow direction, and can be solved analytically [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF][START_REF] Brassart | On convergence properties of variational constitutive updates for elasto-visco-plasticity[END_REF]. The updated flow direction is shown to be equal to the trial one

M n+1 = M trial n+1 = 3 2 dev[σ] dev[σ] trial n+1 , ( 115
)
where σ is the Cauchy stress tensor, resulting in the well-known radial return algorithm [START_REF] Wilkins | Methods in computational physics[END_REF]. Therefore, pseudostresses are useless if a von Mises-type flow rule is retained; the unknown vector consists now of only two scalar unknowns q n+1 = {T, ∆p} n+1 , namely the temperature and the increment of the cumulated viscoplastic strain, which makes the integration algorithm even more efficient. Otherwise, for more complex yield functions, the description of the flow direction through pseudo-stresses remains, and Equation ( 103) should be adapted and solved together with Equations ( 98) and [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF].

Since the small strain flow rule has now an additive structure in the linearized geometrical framework

ε p n+1 = ε p n + ∆pM trial n+1 , (116)
and replaces the exponential formula (96), where ε p stands for the plastic strain, the derivative of the elastic part of the free energy (100) is now simplified as [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF][START_REF] Brassart | On convergence properties of variational constitutive updates for elasto-visco-plasticity[END_REF]]

∂W e ∂∆p n+1 = ∂W e ∂ε p n+1 : ∂ε p ∂∆p = -σ n+1 : M trial n+1 = -σ eq n+1 , (117)
which enters the stationarity condition with respect to ∆p, namely Equation (98). Once the unknown vector q n+1 computed, the Cauchy stresses are then updated as follows

σ n+1 = ∂W ∂ε n+1 = T n T n+1 ∂W ∂ε n+1 , (118)
where ε denotes the total strain.

The extended Johnson-Cook viscoplastic solid

In the rheological equation proposed by Johnson and Cook [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF], the flow stress is assumed to be a function of the cumulated viscoplastic strain p = t 0 ṗ(τ)dτ, of the effective viscoplastic strain rate ṗ, and of the temperature T :

σ y (p, ṗ, T ) = (A + Bp m )(1 + C ln ṗ * )(1 -θ q * ) (119)
where A, B, m, C and q are constitutive parameters. The dimensionless effective viscoplastic strain rate ṗ * = ṗ ṗ0 is also introduced, involving some reference effective viscoplastic strain rate ṗ0 . In practice, the flow stress (119) is valid for ṗ > ṗ0 , otherwise ṗ * = 1 can be used. Next, θ * refers to some dimensionless temperature, defined as

θ * = T -T t T melt -T t ( 120
)
where T melt and T t refer to a melting temperature and a transition temperature respectively. If θ * ≤ 0, it is assumed that no thermal softening occurs, then f th (T ) = (1 -θ q * ) can be set to one. At the other end of the range, if the temperature is greater than the melting temperature, yielding θ * ≥ 1, the flow stress vanishes and f th (T ) = 0.

This model was first recast in a consistent thermodynamical framework in [START_REF] Ranc | Calorimetric consequences of thermal softening in johnson-cook's model[END_REF]. Next, an extension of it was also proposed in [START_REF] Su | Energy-based variational modeling of adiabatic shear bands structure evolution[END_REF] and set into a variational framework. The viscoplastic free energy density and the dissipation pseudo-potential then reads:

W p (p, T ) = A s p + B s m + 1 p m+1 (1 -θ q * ) (121) φ( ṗ; p, T ) = {[A d + B d p m ] ṗ + [A + Bp m ]C ṗ0 ( ṗ * ln( ṗ * ) -ṗ * + 1)} (1 -θ q * ) (122)
where A s , A d , B s , B d are constitutive parameters associated with stored and dissipative contributions, coupled such that

A = A s + A d B = B s + B d . (123)
The Johnson-Cook flow stress (119) is then obtained by the stationarity condition [START_REF] Grüneisen | Theorie des festen zustandes einatomiger elemente[END_REF] with respect to the effective viscoplastic strain rate ṗ:

∂L ∂ ṗ = ∂W ∂p + ∂φ ∂ ṗ = -Σ eq (Σ -Q k ) + (A s + B s p m )(1 -θ q * ) + [(A d + B d )p m ) + (A + Bp m )C ln ṗ * ](1 -θ q * ) = 0 (124)
which gives the expected flow stress (119):

⇔ Σ eq (Σ -Q k ) = (A + Bp m )(1 + C ln( ṗ * )(1 -θ q * ) = σ y (p, ṗ, T ) (125)
The interest of the split of parameters A and B into stored and dissipative components (123) lies in that it allows to weight more either the stored or the dissipative energies in the thermomechanical response. The distribution in the Johnson-Cook flow stress (119) considered in [START_REF] Ranc | Calorimetric consequences of thermal softening in johnson-cook's model[END_REF] is retrieved by setting A s = 0 and B d = 0. In the discrete setting, the stationarity equation ( 98) is computed with the potentials (121) and (122) :

∂I ∂∆p = T n T n+1 ∂W e ∂∆p + (A s + B s p m n+1 (1 -θ q * ) n+1 + (A d + B d p m n+α ) + (A + Bp m n+α )C ln ∆p ∆t ṗ0 (1 -θ q * ) n+α + α∆t mB d p m-1 n+α ∆p ∆t + mBp m-1 n+α C ṗ0 ∆p ∆t ṗ0 ln ∆p ∆t ṗ0 - ∆p ∆t ṗ0 -1 (1 -θ q * ) n+α = 0 (126)
where ∂W e ∂∆p is given by Equation (100), or by Equation (110) in the fully isotropic case. The discrete equation ( 126) is consistent with its continuous counterpart:

lim ∆t→0 ∂I ∂∆p = -Σ n+1 : ∂ f ∂Σ Σ + (A s + B s p m n+1 )(1 -θ q * ) n+1 + (A d + B d p m n+1 ) + (A + Bp m n+1 )C ln ∆p ∆t ṗ0 (1 -θ q * ) n+1 = -Σ eq n+1 + (A + Bp m n+1 )(1 + C ln( ṗ * )(1 -θ q * ) n+1 = 0 (127)
which is identical to Equation (124).

6. The Flux-difference splitting Finite Volume Method

Decomposition of interface fluxes into waves and fluctuations

The finite volume method, like other mesh-based methods, is based on the subdivision of the computational domain into elementary cells. In cell-centered versions of finite volume methods, an approximation U I of the vector of the conserved quantities U (10) is defined in each cell I by integral averaging. If we consider the quadrangular grid cell I shown in Figure 1, of area |A I |, each edge s (1 ≤ s ≤ 4) of outward unit normal n s and of length L s joins the points P s and P s+1 . The integration of the system of conservation laws (9) over the grid cell I yields the following system of ordinary differential equations:

dU dt I = - 1 |A I | N s=1 L s F s ( 128
)
where F s , 1 ≤ s ≤ N, denote the numerical fluxes defined at cell interfaces. The order of accuracy, the physical content but also the computation cost of the finite volume method essentially result from the definition of these numerical fluxes. Commonly, the approach consists in defining a Riemann problem at each cell interface, whose approximate

A + 1 ∆U 1 A - 1 ∆U 1 A + 3 ∆U 3 (1) (2)
(3) (4)

A + 2 ∆U 2 A + 4 ∆U 4 n 2 n 3 n 4 I Figure 2
: Fluctuations defined at each cell interface solution allows to compute these fluxes. For instance, when they are computed with the stationary solution (x/t = 0) of the Riemann problem, the well known Godunov's method [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF] is obtained. The latter is also retrieved by the flux-difference splitting formulation, introduced by Leveque [START_REF] Leveque | Wave propagation algorithms for multidimensional hyperbolic systems[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF], which splits the interface numerical fluxes into fluctuations, hence accounting for waves contributions. These fluctuations are denoted by the operators A ± k ∆U k , and the weighted sum of numerical fluxes expresses as a function of them as:

N s=1 L s F s = P k=1 L k A + k ∆U k + Q l=1 L l A - l ∆U l (129)
where P + Q = N, N being the number of edges of grid cell I. This summation is performed on negative fluctuations for the Q edges having an outward unit normal, and on positive fluctuations for the P edges having an inward unit normal. These fluctuations provide the contribution of first order numerical fluxes to grid cell I.

Assuming that the edge k of unit normal n k has left (L) and right (R) states known in adjacent grid cells (see Figure 1), rightward (+) and leftward (-) fluctuations defined in the local frame of edge k read

A + k ∆U k = (F(U R ) -F(U * )) • n k A - k ∆U k = (F(U * ) -F(U L )) • n k (130)
where U * denotes the stationary solution (given for x/t = 0) of the Riemann problem, which can be computed after having projected the jump of the (averaged) conserved vector

∆U k = (U R -U L ) k across the edge k onto the considered characteristic basis R (p) k ≡ R (p) (n k ) ∆U k = M w p=1 W (p) k = M w p=1 α (p) k R (p) k = R k α k , (131)
where M w is the number of waves, then determining the coefficients α (p) k the wave strengths W (p) k , 1 ≤ p ≤ M w , consist of.

High order fluxes

The class of total variation non-increasing methods [START_REF] Leveque | Wave propagation algorithms for multidimensional hyperbolic systems[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] represents one way to improve the above first order scheme, and allows to meet both a high order of accuracy in areas where the solution is regular together with a high resolution of discontinuities without spurious numerical oscillations when they occur. This class of methods can be implemented by adding high order fluxes to first order ones, which are limited so that a non-increasing total variation of the numerical solution be satisfied at each time step. These second order numerical fluxes read

FHO k = 1 2 M w p=1 |λ (p) k | 1 - ∆t ∆s k |λ (p) k | W(p) k , (132)
where ∆s k refers to the distance between barycenters of grid cells sharing edge k, as shown in Figure 3, and

W(p) k = α(p) k R (p)
k denotes the limited wave strength. Waves are limited based on an upwind ratio θ (p) k defined for the wave p at

e 1 e 2 n k L R (k) (l) if λ (p) k < 0 (l) if λ (p) k > 0 n l n l + G L + G R ∆s k
θ (p) k = W (p) l (n k ) • W (p) k W (p) k 2 ( 133
)
where l denotes the upwind edge, or more precisely either the opposed edge of grid cell L to edge k if λ (p) k > 0, or the opposed edge of grid cell R to edge k if λ (p) k < 0, see Figure 3. The upwind ratio (133) can be understood as a certain measure of the local regularity of the solution. For noncartesian quadrangles, upwind and downwind edges do not necessarily share the same normal. Following [START_REF] Heuzé | Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF], the computation of the upwind ratio (133) is performed with wave strengths recomputed in the same local reference frame of edge k. The weighting coefficients α (p) l of the wave strengths W (p) l express in the local frame of edge k as:

α l (n k ) = R -1 (n k ) • ∆U l (134)
where ∆U l is the jump across edge l of the conserved vector. The wave strengths associated with edge l are then corrected when expressed in the frame of edge k as:

W l (n k ) = diag (α l (n k )) • R(n k ) = diag ([K(n k )] -1 • ∆U l ) • R(n k ) ( 135
)
where W l (n k ) is the matrix whose columns are wave strength vectors W (p) l , 1 ≤ p ≤ M w . The wave strength W(p) k of wave p associated with edge k is limited using some classical limiting function φ(θ (p) k ) applied to wave coefficients such as:

α(p) k = φ(θ (p) k )α (p) k , (136)
many of which permit to obtain different known finite volume schemes [START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF]. Here, the classical minmod limiter is used [START_REF] Trangenstein | Numerical solution of hyperbolic partial differential equations[END_REF]: φ(θ) = max(0, min(1, θ)).

(137)

Transverse fluxes

Transverse numerical fluxes can also be added to first and second order ones to improve the stability of the numerical scheme and increase the time step, so that the Courant number can theoretically be set at unity. These fluxes permit to account for information travelling in bias with respect to the considered grid cell, through the neighbour cell only sharing a node (but not an edge) with it. This method is known after [START_REF] Collela | Multidimensional upwind methods for hyperbolic conservation laws[END_REF] as the Corner Transport Upwind (CTU) method. For elastic media, these fluxes permit to significantly improve the maximum allowable Courant number since elasticity couples strain components through Poisson's effect. Consider the patch of grid cells shown in Figure 4, which share the common edge (k) of unit normal n k . The normal fluctuations A + k ∆U k and A - k ∆U k contributing respectively to grid cells R and L can be split into transverse fluctuations contributing to neighboring cells of cell L across the edges (i) and ( j) , and across edges (l) and (m) for cell R. These normal fluctuations are first projected on the characteristic basis computed with the approximate Riemann solvers defined on the adjacent edges. The negative normal fluctuation is for instance decomposed on the characteristic basis associated with edge (i) as

e 1 e 2 n k A + k ∆U k A - k ∆U k L R T (k) (i) ( j) (l) (m) B + i A - k ∆U k B + j A - k ∆U k B + m A + k ∆U k B + l A + k ∆U k n i n j n m n l
A - k ∆U k = M w p=1 β p R (p) i = R i β (138)
where R i accounts here for the normal n i of edge (i), but also of different material properties between grid cells L and T . For plane strain linear elasticity for example, the coefficients β p can be determined analytically. The transverse fluctuations are then computed with the positive operator B + , more precisely only waves with positive characteristic speeds in the local frame of the edge will contribute to this transverse fluctuation

B + i A - k ∆U k = M w p=1 λ + p β p R (p) i , (139)
from which the additional numerical transverse flux is built as:

Ftran i = ∆t 2∆s i B + i A - k ∆U k (140)
which contributes to grid cell T .

Explicit time integration

Gathering first order fluctuations and additional numerical fluxes, and considering an explicit Euler time integration, the state of grid cell i is updated at time t n+1 with the following formula:

U n+1 I = U n I - ∆t |A I |         P k=1 L k A + k ∆U k + Q l=1 L l A - l ∆U l         - ∆t |A I |         P k=1 L k Fout k - Q l=1 L l Fin l         (141)
where Fin l and Fout k refer to inward and outward additional numerical fluxes, respectively associated with the Q and P edges having either an outward or inward normal relative to grid cell I. These additional numerical fluxes sum both second order and transverse ones Fin l = FHO l + Ftran l .

(142)

Numerical examples

Test cases at one material point

At first, elementary loading paths are investigated at the scale of one material point. Since the discrete constitutive update is driven by updated values of the internal energy density E n+1 and the deformation gradient F n+1 at time t n+1 , a certain consistency between both should be enforced to provide the material point with some physically compatible loading paths. This is performed at each loading step by solving the local discrete constitutive problem [START_REF] Heuzé | Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF] with input data satisfying the discrete balance of internal energy written between times t n and t n+1 , still considering adiabatic conditions ∆E = P : ∆F,

from which the stresses P are computed so that

P = ∂ ∂F W ∆E, ∆F ∆t n+1 . ( 144
)
The stresses (144) are computed by means of a fixed point loop, until the balance of internal energy (143) is satisfied. Notice that Equation ( 143) is only valid for smooth solutions, assumed here for the purpose of carrying out test cases on one material point, while the conservation of the total energy ( 7) is rather solved through some discrete conservative scheme when considering a whole domain, which is valid whatever the smoothness of the solution. Each loading path is computed with a constant strain rate ∆F /∆t and a constant time step ∆t.

One-dimensional strain test

Let's start with a one-dimensional strain test, whose deformation gradient is of the form

F = F 11 e 1 ⊗ e 1 + e 2 ⊗ e 2 + e 3 ⊗ e 3 , (145)
simulating the kinematics of a plane wave. The longitudinal stretch F 11 is driven so that to achieve a cycle. The Mie-Grüneisen equation of state [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF], the neo-Hookean shear elastic energy [START_REF] Ortiz | The computation of the exponential and logarithmic mappings and their first and second linearizations[END_REF], the Mises strength criterion and the Johnson-Cook flow stress with the original distribution between stored energy and dissipated heat (i.e. A s = 0, B d = 0) are used. The values of material parameters considered for this test are gathered in Table 1. Moreover, a constant stretch rate of Ḟ11 = 10 3 s -1 is considered during the loading cycle.

Elasticity parameters Figure 5 shows different quantities plotted as a function of this longitudinal stretch during the cycle of loading. The longitudinal Cauchy stress component shows the effect of the pressure, while its deviatoric part follows a classical elastic-plastic open cycle due to isotropic hardening. Both the entropy density and the cumulated viscoplastic strain follow a monotonic increasing evolution along the loading cycle, while the convex nature of both internal and free energies are shown. The temperature evolution results from both the thermoelastic effects (thermal dilatation) computed through the equation of state, and from the work generated by viscoplasticity. At last, the values reached by the incremental energy are also shown, and lie at several orders of magnitude below these of internal and free energy densities.

E = 200 GPa ν = 0.3 ρ 0 = 7800 m 3 a = -1 Thermal parameters C = 452 J/(kg.K) α = 12 × 10 -6 K -1 q = 1 T 0 = 293 K Johnson-Cook A = 400 MPa B = 770 MPa m = 0.557 C = 5 × 10 -3 q = 1 T t = 290 K T melt = 1800 K ṗ0 = 10 -3 s -1
Next, Figure 6 shows convergence curves of the solution extracted at the end of the cycle as a function of the time step ∆t, plotted for the different quantities. Relative errors are computed with respect to a numerical solution obtained with a finer time step. Very little effect of the parameter α (see Section 4) is actually observed for the Johnson-Cook dissipation pseudo-potential as the three convergence curves are superposed. As expected from the first order accurate incremental variational update [START_REF] Heuzé | Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF], a convergence rate of about one is observed for all plotted quantities.

Pure shear test

A pure shear test is then considered, with a deformation gradient of the form

F = 1 + F 12 e 1 ⊗ e 2 . (146)
The slip F 12 is also driven so that to achieve a cycle. The same equation of state, shear elastic energy, strength criterion and energy distribution in the Johnson-Cook flow stress are still used, as well as the corresponding set of material parameters shown in Table 1. Moreover, a constant slip rate of Ḟ12 = 10 3 s -1 is also imposed during the loading cycle. - monotonic increasing evolution along the loading cycle. The temperature now rises along the cycle due to the work generated by viscoplasticity. Figure 6 shows similar convergence curves of some relative errors of different quantities for the shear test, computed with the solution extracted at the end of the cycle. The same comments apply where a convergence rate of about one is observed for all plotted quantities. The parameter α now shows some small differences on convergence curves, although very limited.

⋅ m -3 ) 1e6 -4 -3 -2 -1 0 1 2 3 4 F12 1e-1 0 2 4 6 8 Cumulated plastic strain 1e-1 -4 -3 -2 -1 0 1 2 3 4 F12 1e-1 -2 -1 0 1 2 F p 12 1e-1 -4 -3 -2 -1 0 1 2 3 4 F12 1e-1 -6 -4 -2 0 2 4 6 inc (J ⋅ m -3 ) 1e5 -4 -3 -2 -1 0 1 2 3 4
The partition of viscoplastic work into dissipated heat and stored energy is also naturally computed with the variational framework [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF][START_REF] Stainier | Study and validation of a variational theory of thermo-mechanical coupling in finite visco-plasticity[END_REF][START_REF] Canadija | On the thermomechanical coupling in finite strain plasticity theory with non-linear kinematic hardening by means of incremental energy minimization[END_REF] through the definition of the plastic part of the free energy (121) and the dissipation pseudo-potential (122). Such partition is known to be non-constant [START_REF] Hodowany | Partition of plastic work into heat and stored energy in metals[END_REF][START_REF] Rosakis | A thermodynamic internal variable model for the partition of plastic work into heat and stored energy in metals[END_REF]. It can be computed for a given constitutive model. More precisely, the mechanical dissipation (23) can be expressed for the Johnson-Cook viscoplastic constitutive model as

D int = - ∂W ∂p ṗ = σ y - ∂W p ∂p ṗ, (147)
with the flow stress σ y given by Equation (119). According to the expression of the plastic part of the free energy (121), it simplifies as

D int = σ y -(A s + B s p m )(1 -θ q * ) ṗ ≡ βσ y ṗ, (148)
where the Taylor-Quinney coefficient β appears as the fraction of viscoplastic power density converted into heat, which reads as

β = 1 - (A s + B s p m )(1 -θ q * ) σ y . (149)
Its expression actually depends on the effective viscoplastic strain rate ṗ, the cumulated viscoplastic strain p and the temperature T . Figure 9 shows several plots of its evolution along the shear loading cycle for different distributions between stored energy and dissipated heat. More precisely, the first case considered is the original distribution of Johnson-Cook, A s = 0 leads to that the density of plastic power associated with the initial and isothermal yield thresh- old A contributes completely to the dissipated heat, while B d = 0 assigns the whole effect of the isotropic hardening to the stored energy. Since on the one hand the cumulated viscoplastic strain p is a monotonic increasing function of the loading (see Figure 7), and on the other hand the effective viscoplastic strain rate ṗ is almost constant (apart during the elastic loading-unloading-reloading parts of the cycle) because the slip rate Ḟ12 is imposed constant, then the Taylor-Quinney coefficient is a monotonic decreasing function during the loading cycle (apart the aforementioned areas). The second case considers that a given fraction (arbitrarily set so that A s /A d = 0.5 here) of the density of plastic power associated with the initial and isothermal yield threshold A contributes to the stored energy. As could be expected, the Taylor-Quinney coefficient is lower than in the previous case, because the density of plastic power contributes less to the dissipated heat. The evolution of the coefficient appears as translated to the lower values, consistently with Equation (149). The third case considers that the effect of isotropic hardening partly contributes to the dissipated heat (arbitrarily set so that B d /B s = 0.5 here). Consequently, the dissipated heat is given more contribution with respect to the first case, and the Taylor-Quinney coefficient becomes greater. The fourth and last case combines the second and the third ones.

Sudden velocity loading and unloading of a heterogeneous volume

Generalizing the test case introduced in [START_REF] Heuzé | Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF], we consider a two-dimensional numerical simulation carried out on a heterogeneous volume, that consists of a circular inclusion of radius R, which is associated a thermo-elasticviscoplastic constitutive response, centered in a square matrix of side length 2a whose constitutive response is thermoelastic. This volume is initially in a natural state. It is suddenly loaded on its left side at time t = 0 with a constant rightward velocity v. After time t u , the applied velocity is unloaded to zero. Symmetry conditions are set at the top and bottom sides of the volume, while transmissive boundary conditions have been set on the right side. Due to the symmetry condition, the computational domain only consists of one half of the whole volume, as shown in Figure 10.

The analysis is carried out in the two-dimensional plane strain case, and within the linearized geometrical framework. The elastic and thermal potentials (113) and (114) are used. In the inclusion, the Mises criterion, the Johnson-Cook plastic potential (121) and the dissipation pseudo-potential (122) are also employed with the original distribution between stored energy and dissipated heat (i.e. A s = 0, B d = 0). The numerical values of parameters considered for tion is conducted with the finite volume method which has been recalled in Section 6, and using a CFL number set at 0.8. It is compared with a numerical solution computed with the finite element method [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF] with Q1 finite elements, coupled with an explicit central difference time integrator, and to a variational constitutive update driven in both strain ε and temperature T [START_REF] Stainier | ZorgLib. User's manual[END_REF]. The details of such formulation can be found in [START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF]. The same elastic, plastic and thermal potentials as well as the dissipation pseudo-potential are also used. Figures 11(a) and 11(b) show some comparison between finite element and finite volume numerical solutions at two different times of computation. The temperature map is shown on the left part of each subfigure, while graphs on their right parts show the superposed plots along the symmetry line of the temperature and of the normal stress component σ 11 . The solution basically consists of a compression slot travelling rightward in the volume. Wave interactions are generated at the matrix/inclusion interface due to their mismatch of elastic impedance. Finite element and finite volume numerical solutions are globally in agreement, althought the prescribed discontinuity of velocity imposed by the loading introduces some spurious numerical oscillations in the finite element one. These oscillations do not appear in the finite volume solution thanks to the use of limiters which prevents from their appearance, enforcing here a TVD criterion [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF].

Figures 12(a) and 12(b) show the maps and the superposed plot along the symmetry line of the cumulated viscoplastic strain in the inclusion at two different instants, computed with the two methods. Right after the compressive wave has gone through the matrix/inclusion interface, the cumulated viscoplastic strain rises in the inclusion. It is slightly overestimated in the finite element solution due to the appearance of spurious oscillations, while both solutions reach the same level once the plastic wave has travelled over a sufficient length in the inclusion. The circular geometry of the inclusion generates a convergent profile at its rear side that the compression slot follows, which generates a net increase of cumulated viscoplastic strain and temperature. Especially, the cumulated viscoplastic strain rises in a quite narrow band, whose distribution is computed a little differently in the two numerical solutions. Figure 13 shows the time evolutions of the kinetic energy, the internal energy, and of the total energy of the volume during the numerical simulation conducted with the finite volume and finite element solvers. After a first rise of energies due to external energy introduced in the system through the non-zero velocity prescribed at the left side of the volume during time t u , the total energy is conserved. Its level reached in the two solutions slightly differs. The finite element solution appears to be slightly more diffusive that the finite volume one as shown in Figure 14, through the time evolution of the global entropy. The latter starts to really increase as soon as the compression slot reaches the matrix/inclusion interface, and that viscoplasticity starts to develop within the inclusion. Its time slope decreases as the compression slot has almost crossed the inclusion.

Impacted bidimensional rectangular specimen

We consider now a bidimensional rectangular specimen, infinite in the third direction, whose half square computational domain of side length l = 3 m is sketched in Figure 15. The specimen is initially animated of a uniform leftward velocity -ve 1 , v = 200 m.s -1 , and hits a rigid plane (without friction) on its left side, while a symmetry condition is accounted for on the bottom side of the square domain. The problem is treated in the plane strain case, and within the large strain framework. The specimen is assigned a thermo-hyperelastic-viscoplastic constitutive response. As for the heterogenous volume, two numerical simulations are computed in this example with the finite volume and the finite element solvers. Since large plastic incompressible strains (due to Mises criterion) occur in this example, some treatment is required to avoid volumetric locking in the two numerical solutions. On the one hand, the conservation law written on the jacobian determinant ( 5) is included in the first order system of conservation laws [START_REF] Bleier | Efficient variational constitutive updates by means of a novel parameterization of the flow rule[END_REF] when computing the finite volume numerical solution. The volumetric part of the free energy [START_REF] Malkus | Mixed finite element methods-reduced and selective integration techniques: a unification of concepts[END_REF], and more precisely here Equation [START_REF] Simo | Computational inelasticity[END_REF], is then computed with the jacobian updated through the associated conservation law. On the other hand, finite elements with selective reduced integration of the volumetric part of the constitutive response (constant pressure in the finite element) are employed [START_REF] Malkus | Mixed finite element methods-reduced and selective integration techniques: a unification of concepts[END_REF]. The CFL number is set at 0.4 for the computation of the two numerical solutions. The variational constitutive update coupled with finite elements is driven with the deformation gradient F and the temperature T [START_REF] Stainier | ZorgLib. User's manual[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF]. The two numerical solutions are computed with the same volumetric (67) and thermal energies [START_REF] Mie | Zur kinetischen theorie der einatomigen körper[END_REF], and with the same Johnson-Cook plastic potential (121) and dissipation pseudo-potential (122), still employed with the same distribution between stored and dissipated heat (i.e. A s = 0, B d = 0). The numerical values employed for the material parameters are the same than these defined in Table 1, except the definition of q which does not exist here since the pressure is here defined with Equation (68).

with e = 1 2 ln[C e ], that may give slightly different results than these obtained with the Neo-Hookean energy (59) employed for the finite volume solution. As noted in Section 5.4.1, the Hencky hyperelastic shear energy (150) was at first implemented in that variational constitutive update [START_REF] Stainier | ZorgLib. User's manual[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF] since the logarithmic elastic strain allows to recover an additive structure between the updated elastic strain and its elastic predictor as in the case of small strain framework, which permits to solve the stationarity condition [START_REF] Taylor | The latent energy remaining in a metal after cold working[END_REF] analytically. Similarly to the previous example, Figure 16 shows the time evolutions of the kinetic energy, the internal energy, and of the total energy of the computational domain during the numerical simulation conducted with the finite volume and finite element solvers. The total energy is perfectly conserved in the finite volume solution here since the loading comes from the initial velocity. That computed with the finite elements is very close to the finite volume one although the heat equation is solved locally rather than the conservation of the total energy, and may be due to the variational solver. In this test case, the finite volume solution appears slightly more diffusive than the finite element one (i.e. the converse of the previous example), as shown by the time evolutions of the internal and kinetic energies, oscillations being associated with wave round trips in the specimen during the time interval. This dissipative behaviour is also shown in Figure 17 by the time evolutions of the global entropy of the system, which is a monotonic increasing function of time as expected. The difference between finite volume and finite element solutions here essentially comes from the limiters used in the former one. Only pressure waves are here limited in the finite volume solver, which allows to keep stability in large strains at the price of more numerical viscosity. Once the first wave round trip performed, the limiters are not active anymore and the differences between time evolutions of the internal and kinetic energies and the entropy do not continue to grow.

Figure 18 shows the maps of the stress component P 11 , the cumulated plastic strain p and the temperature T , here plotted at time t = 1 × 10 -3 seconds in the initial configuration of the bidimensional specimen, which is the one where the two solutions are computed. Row outputs consist of cell data for the finite volume solver and nodal data for the finite element one. The two solutions are in qualitative agreement, although some differences appear at the top left corner on the temperature and the cumulated plastic strain fields. Figure 19 shows superposed graphs of the two solutions at the same time plotted along the bottom line of the square computational domain. The profiles shown of the stress component P 11 , the cumulated plastic strain p and the temperature T computed with the two methods are globally in agreement.

At last, Figures 20 and21 show the maps of the cumulated plastic strain p and the temperature at two different times of the deformation process, here plotted in the current configuration of the bidimensional specimen. For visualization purpose and comparison with the finite element solution, cell data extracted from the cell-centered finite volume solution are averaged at nodes to compute the deformed shape. The two solutions are qualitatively in agreement, although some differences appear here and there, especially close to the bead formed at the top left corner of the specimen. The temperature, the cumulated plastic strain and the deformed shape slightly differ in that area, which may come from the different spatial discretizations (having different numerical viscosities), the different variational constitutive updates, and the different expressions of the shear part of the hyperelastic free energy.

Multi-holed elementary cell

A last example involving large strains consider a multi-holed elementary cell submitted to some impact loading. The considered computational domain with associated boundary conditions is sketched in Figure 22. It is drilled with several circular holes where free boundary conditions are considered. The domain is initially in a natural state. A constant leftward velocity -ve 1 is suddenly applied at its right side at time t=0, until it is released to zero at time t u . As in Section 7.3, the problem is treated in the plane strain case within the large strain framework. The same formulations are adopted for the finite volume and the finite element solvers, as well as for the constitutive response whose values of material parameters are gathered in Table 1. Table 3 summarizes the numerical values of geometrical and loading parameters of this test case. Figure 23 shows the maps of the jacobian, the cumulated plastic strain and the temperature plotted in the initial configuration of the multi-holed cell at time t = 4 × 10 -6 seconds. The two solutions are in qualitative agreement, although one can observe that the finite volume solution predicts sligthly higher temperatures and lower cumulated plastic strains that these computed with finite elements. The jacobian J is computed from the solution of its associated conservation law [START_REF] Bartels | On the thermomechanical coupling in dissipative materials: a variational approach for generalized standard materials[END_REF] in the finite volume solution, whereas it is computed as the determinant of the deformation gradient J F = det F in the finite element one.

Next, Figures 24 and25 show the maps of the cumulated plastic strain p and the temperature at two different times of the deformation process, here plotted in the current configuration of the multi-holed elementary cell. The two solutions are qualitatively in agreement, although the same tendency of the finite volume solution to predict sligthly higher temperatures and lower cumulated plastic strains that these computed with finite elements still appear. Figure 24 shows the multi-holed cell after a continuous compression on its right side during the first 1 × 10 -5 seconds, just before the prescribed velocity is set at zero. The contour of the deformed holes computed with the two solutions look in agreement. After the unloading of the applied velocity, the two rightmost holes show a slightly re-increased volume in Figure 25, whereas the two leftmost ones have been compressed.

Finally, Figure 26 shows the time evolutions of the different energies of the system during the numerical simulations conducted with the finite volume and finite element solvers. The prescribed velocity at the right side of the cell makes the total energy evolving in time as for the inclusion test case (either increasing or decreasing according to wave round trips in the computational domain), while it becomes constant as soon as the velocity is set at zero. Figure 27 shows the associated time evolutions of the global entropy of the cell, which still follows a monotonic increasing trend in time.

Conclusion

A variational framework has been proposed in this work for the constitutive update of thermomechanical constitutive models when coupled with a set of hyperbolic conservation laws. Both a continuum and a consistent first order accurate discrete settings were derived. More precisely, the constitutive update is driven by the rates values of some strain measure and the internal energy density in the continumm setting, and by the updated values at some discrete time of these strain measure and internal energy density in the discrete setting. This is achieved by defining a Lagrangian functional enforcing the rate of the residual of the Legendre transform of Helmholtz's free energy to vanish through a Lagrange multiplier. Contrary to other thermomechanical formulations, the temperature is not part of the input data of the constitutive update, although it appears as one of the unknown of the solution process. Doing so, the two-field temperature formulation introduced in [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF] becomes now useless.

The main interest of such formulation lies in that it is naturally compatible with the writing of a set of conservations laws, which reduce to the well-known Rankine-Hugoniot jump conditions across discontinuities of the fields. Hence in the discrete setting, the use of any conservative numerical scheme will ensure that the right shock speeds will be computed, while the proposed variational constitutive update will permit to account for dissipative and thermomechanically coupled constitutive reponse. In addition, the total energy will be conserved in the discrete setting since the associated conservation law is solved.

The proposed variational approach is illustrated in this work with a thermo-hyperelastic-viscoplastic constitutive response. The elastic energy consists of a volumetric component related to some equation of state, and of some hyperelastic model for the shear component. The flow rule direction is parameterized with pseudo-stresses as proposed by Mosler & co-workers, resulting in the solution of an unconstrained optimization problem. For a fully isotropic medium, the discrete stationarity conditions consist in solving a system of three scalar nonlinear equations. The proposed discrete variational solver is coupled with the second order accurate flux difference splitting finite volume method for the solution of the set of conservation laws. A set of numerical test then show on the one hand the convergence of the discrete variational constitutive update on test cases conducted at one material point, and on the other hand some numerical solutions which are in agreement with these obtained with with finite elements coupled to an explicit time-stepping and to a temperature-driven variational constitutive update.

Figure 1 : 1 e 2 n 1

 1121 Figure 1: Quadrangular finite volume

Figure 3 :

 3 Figure 3: Illustration of upwind edges of edge k for wave comparison purpose

Figure 4 :

 4 Figure 4: Normal and transverse fluctuations defined from edge k.

22 0Figure 5 :

 225 Figure 5: Plots of the evolution of different quantities as a function of the stretch F 11 over a loading cycle performed at one material point in the case of a one-dimensional strain test.

Figure 7 Figure 6 :

 76 Figure 6: Plots of convergence curves of different quantities in the case of a one-dimensional strain test, with the solution extracted at the end of the loading cycle performed at one material point. Convergence curves are identical for all values of α in this test case. Only those obtained for α = 1 are traced here.

F12

F12Figure 7 :Figure 8 :

 78 Figure 7: Plots of different quantities as a function of the slip F 12 over a loading cycle performed at one material point in the case of a shear test.

25 - 5 Figure 9 :

 2559 Figure 9: Taylor-Quinney coefficient computed with the Johnson-Cook viscoplastic constitutive model along the shear loading cycle, with different distributions between stored energy and dissipated heat.

Figure 10 :

 10 Figure10: Half of a heterogeneous volume that consists of a circular inclusion centered in a matrix, suddenly loaded then unloaded after a time t u on its left side by a prescribed velocity v.

  (a) Time t = 1.63838 × 10 -7 seconds. (b) Time t = 2.88355 × 10 -7 seconds.

Figure 11 :

 11 Figure 11: Comparison of the numerical solutions computed with the finite volume method (FVM) and the finite element method (FEM). Maps on the left part of each subfigure show the temperature fields computed with the two methods. Graphs on their right parts show the superposed plots for the two solutions of the temperature and of the longitudinal stress component σ 11 , both plotted along the symmetry line.

  (a) Time t = 1.63838 × 10 -7 seconds. (b) Time t = 2.88355 × 10 -7 seconds.

Figure 12 :

 12 Figure 12: Comparison of the numerical solutions computed with the finite volume method (FVM) and the finite element method (FEM). Maps on the left part of each subfigure show the cumulated plastic strain fields in the circular inclusion computed with the two methods. Graphs on their right parts show its superposed plots along the symmetry line for the two solutions.

Figure 13 :Figure 14 :

 1314 Figure 13: Time evolutions of the global kinetic, internal and total energies integrated on the half of the heterogeneous volume during the finite volume (FVM) and the finite element (FEM) numerical simulations.

1 e 2 Figure 15 :

 1215 Figure 15: Sketch of the bidimensional square computational domain.

Figure 16 :Figure 17 :

 1617 Figure 16: Time evolutions of the global kinetic, internal and total energies integrated on the impacted bidimensional specimen during the finite volume (FVM) and the finite element (FEM) numerical simulations.

Figure 18 :Figure 19 :

 1819 Figure 18: Maps of the stress component P 11 , the cumulated plastic strain p and the temperature T plotted in the reference configuration of the bidimensional specimen, computed with the finite volume method (FVM) and the finite element method (FEM) at time t = 1 × 10 -3 seconds.

Figure 20 :

 20 Figure 20: Maps of the temperature and the cumulated plastic strains plotted in the current configuration of the bidimensional specimen at time t = 1 × 10 -3 seconds, computed with the finite volume method (FVM) and the finite element method (FEM).

Figure 21 :

 21 Figure 21: Maps of the temperature and the cumulated plastic strains plotted in the current configuration of the bidimensional specimen at time t = 3 × 10 -3 seconds, computed with the finite volume method (FVM) and the finite element method (FEM).

1 Figure 22 :

 122 Figure22: Sketch of the computational domain of the multi-holed elementary cell. Holes 1, 2, 3 of radii r = r 1 = r 2 = r 3 are centered at locations (w, 3w), (3w, w) and 4w,7w 2 respectively, where w ∈ R 3 , min 2 7 (Hr), (L-r) 4 is a parameter.

Geometry R = 8 ×H = 9w 2 Table 3 :

 823 10 -3 m r = 4 × 10 -3 m w = 5 × 10 -3 m L = 11w 2 Loading v = 200 m.s -1t u = 1 × 10 -5 s Geometrical and loading parameters for the multi-holed medium test case.

Figure 23 :

 23 Figure 23: Maps of the jacobian, the cumulated plastic strain and the temperature plotted in the initial configuration of the multi-holed elementary cell at time t = 4 × 10 -6 seconds, computed with the finite volume method (FVM) and the finite element method (FEM).

Figure 24 :

 24 Figure24: Maps of the cumulated plastic strain and the temperature plotted in the current configuration of the multi-holed elementary cell at time t = 1 × 10 -5 seconds (just before the unloading), computed with the finite volume method (FVM) and the finite element method (FEM).

Figure 25 :

 25 Figure 25: Maps of the cumulated plastic strain and the temperature plotted in the current configuration of the multi-holed elementary cell at time t = 1.4 × 10 -5 seconds, computed with the finite volume method (FVM) and the finite element method (FEM).

Figure 26 :Figure 27 :

 2627 Figure 26: Time evolutions of the global kinetic, internal and total energies integrated on the multi-holed cell during the finite volume (FVM) and the finite element (FEM) numerical simulations.

Table 1 :

 1 Material parameters for the convergence tests at one material point.

Table 2 :

 2 Parameters for the heterogeneous volume test case.