A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors

This chapter provides a general and self-contained overview of the variational approach to nonlinear dissipative thermo-mechanical problems initially proposed in Ortiz and Stainier (1999) and Yang, Stainier, and Ortiz (2006). This approach allows to reformulate boundary-value problems of coupled thermo-mechanics as an optimization problem of an energy-like functional. The formulation includes heat transfer and general dissipative behaviors described in the thermodynamic framework of Generalized Standard Materials. A particular focus is taken on thermo-visco-elasticity and thermo-visco-plasticity. Various families of models are considered (Kelvin–Voigt, Maxwell, crystal plasticity, von Mises plasticity), both in small and large strains. Time-continuous and time-discrete (incremental) formulations are derived. A particular attention is dedicated to numerical algorithms which can be constructed from the variational formulation: for a broad class of isotropic material models, efficient predictor–corrector schemes can be derived, in the spirit of the radial return algorithm of computational plasticity. Variational approximation methods based on Ritz–Galerkin approach (including standard finite elements) are also described for the solution of the coupled boundary-value problem. Some pointers toward typical applications for which the variational formulation proved advantageous and useful are finally given.


INTRODUCTION
Variational principles have played an important role in mechanics for several decades, if not more than a century (see for example Lanczos, 1986or Lippmann, 1978. They have been mostly developed, and widely used, for conservative systems: the most eminent examples are Hamilton's principle in dynamics and the principle of minimum potential energy in statics. Some variational principles with application to dissipative systems have been around for a long time as well, such as principles of maximum dissipation for limit analysis (notably in plasticity).
Variational approaches present many attractive features, especially regarding the possibilities that they offer for mathematical analysis, but also for numerical approaches. They open an easier way to unicity, convergence, and stability analysis of mathematical formulations and associated numerical methods. This has motivated a very large quantity of published work and an exhaustive review is thus out of the scope of this chapter. To directly focus on the category of variational approaches envisioned here, let us then simply say that, following the pioneering work of Biot (1956Biot ( , 1958, the variational form of the coupled thermo-elastic and thermo-visco-elastic problems has been extensively investigated (see for example Batra, 1989;Ben-Amoz, 1965;Herrmann, 1963;Molinari & Ortiz, 1987;Oden & Reddy, 1976). On the other hand, several authors have proposed variational principles for the equilibrium problem of general dissipative solids in the isothermal setting: see for example Carini (1996), Comi, Corigliano, & Maier (1991), Hackl (1997), Han, Jensen, & Reddy (1997), Martin, Kaunda, & Isted (1996), Mialon (1986), Ortiz & Stainier (1999), in elasto-visco-plasticity, and also Balzani & Ortiz (2012), Bourdin, Francfort, & Marigo (2008), Francfort & Marigo (1998), Kintzel & Mosler (2010, 2011, in brittle and ductile damage. By contrast, the case of thermo-mechanical coupling (i.e. with conduction) in these latter classes of dissipative materials has received comparatively less attention (cf. Armero & Simo, 1992, 1993Simo & Miehe, 1992, for notable exceptions).
This chapter is intended to provide an overview of recent and less recent work by the author and colleagues on a specific variational approach (initially described in Yang et al., 2006) to coupled thermo-mechanical problems involving nonlinear dissipative behaviors, such as thermo-visco-elasticity and thermo-elasto-visco-plasticity. It will also be the occasion to fill a few gaps between previously published material, in particular by providing a more detailed account of thermal coupling aspects for a variety of constitutive models written under variational form. Links toward closely related work by other researchers are also provided.
We start in Section 2 by setting the general thermodynamic modeling framework serving as a foundation for the proposed variational formulation of coupled thermo-mechanical boundary-value problems. This framework follows closely that of Generalized Standard Materials (GSM) (Halphen & Nguyen, 1975), with a local state description based on internal variables. We will also recall some elements of finite transformations kinematics, as well as balance equations in Lagrangian and Eulerian formulations (although we will mostly work in a Lagrangian setting). This part is quite standard, with a few departures from the mainstream approach (e.g. the use of a Biot conduction potential). We then proceed (Section 3) to reset the constitutive and balance equations defining thermo-mechanical boundary-value problems under a variational form. By variational, we here understand that the problem is formulated as an optimization (or at least a stationarity) problem, with respect to fields of state variables. Since we work with a local state description based on internal variables, we will split the presentation in two parts: first, the local constitutive problem (determination of internal variables or their rate), and, second, the boundary-value problem (determination of fields of external variables). Each part is itself structured in two subparts: we first present the time-continuous evolution problem and its variational formulation, followed by the time-discrete (or incremental) variational formulation. This structure of presentation, which somewhat differs from that adopted in Yang et al. (2006), allows to show that the variational boundary-value problem is formally identical to a thermo-elastic problem, internal variables being handled locally through a nested constitutive variational problem. This is probably the most interesting result presented in that specific part of the paper. After presenting the variational formulation for the incremental boundary-value problem, we show how to add more complex thermo-mechanical boundary conditions, such as mixed thermal conditions (e.g. convective exchange). The variational formulation of coupled thermomechanical boundary-value problems is initially presented within a quasistationary context (yet including combined heat capacity and conduction effects), but we show in subsection 3.3 how it can be extended to account for inertia effects in the time-discrete framework. We conclude this part by describing linearization procedure in the case of infinitesimal (small) displacements and temperature variations. Note that nonlinearities can remain within this "linearized" context, due to the presence of thermo-mechanical coupling terms.
In Sections 4 and 5, we look in more details at the variational formulation of (continuous and incremental) constitutive equations for some specific models in thermo-visco-elasticity and thermo-elasto-visco-plasticity. We start with the simplest thermo-visco-elastic model possible: linearized kinematics Kelvin-Voigt model with linear elasticity and viscosity. Given its relative simplicity of formulation, the variational update for this constitutive model is treated in details, including all possible temperature dependence effects (thermo-elasticity, thermal softening of elastic and viscous moduli). Note that this is the only model for which a complete treatment is provided here, some simplifying hypotheses being taken for later, more complex, models, for the sake of clarity in the presentation. We then move on to generalized Maxwell models, which introduce internal variables (viscous strains). The analysis of Kelvin and Maxwell models is then repeated for finite strains kinematics. Maybe the most interesting point in that part is the fact that, for isotropic materials at least, constitutive updates can be reduced to solving a reduced number of scalar equations by adopting a spectral approach [as given in Fancello, Ponthot, and Stainier (2006)], independently of the complexity of elastic and viscous potentials adopted in the models. Section 5 deals with thermo-elasto-visco-plasticity. We start with crystal plasticity, which describes fine scale behavior of crystalline materials (mostly metals for our purpose, but also some organic materials and rocks). For this class of models, the variational formulation offers the power of optimization algorithms to solve the complex problem of determining (incremental) plastic slip activity, a problem which can become quite acute when considering complex latent hardening models. The section continues by considering plasticity models at the macroscopic scale, chiefly von Mises ( J 2 ) plasticity. We look at both the time-continuous and time-discrete variational formulations of linearized kinematics and finite strain thermoelasto-visco-plasticity. Special emphasis is given to the fact that, in the isotropic case, the variational update can be interpreted as a radial return algorithm. Another important result is that, admitting a specific choice of elastic free energy potential, finite strain kinematics can be uncoupled from the constitutive update itself, which can then be treated as in the linearized kinematics case. This result, initially shown in Ortiz and Stainier (1999) in the isothermal case, is here described in a coupled thermo-mechanical setting. We complete this section on thermo-elasto-visco-plasticity by briefly illustrating how to consider general elastic free energy potentials (Fancello, Vassoler, & Stainier, 2008b) or more general plastic flow rules. The section closes by a discussion on the problem of partition of plastic work into stored and dissipated energy, and how it is naturally and implicitly treated within the present formulation (Stainier & Ortiz, 2010).
In Section 6, we complete the description of our variational framework by considering numerical approximation methods for solving coupled thermo-mechanical boundary-value problems. The variational formulation presented before is particularly well suited to finite element approaches, and we first describe what results from a standard Galerkin formulation. We also look at mixed formulations for handling (nearly) incompressible behaviors (occurring in plasticity, for example). The variational formulation can also be exploited by using more general Ritz-Galerkin approaches, and we describe in particular such approaches applied to the case of general discontinuities/interfaces such as shear bands or cohesive zones.
Section 7 provides some pointers toward typical applications for which the variational formulation proved advantageous and useful. We focus in particular on multiscale and adaptive variational approaches. Finally we conclude by summarizing the main features and advantages of the variational approach presented in this chapter and discuss some related open problems.

Local thermodynamic model
We consider thermo-mechanical problems and place ourselves in a continuum modeling framework. Let us then consider a thermo-mechanical problem defined on a body occupying a domain B 0 2 R 3 in its reference configuration (which, for simplicity, we will identify with the configuration occupied at the initial time t ¼ 0). The solution of the thermo-mechanical problem is described by the displacement mapping material points X to spatial points x: and the absolute temperature field: where [0, t f ] is the time interval under consideration. The domain occupied by the body in the deformed configuration at time t will be denoted by B t (or more simply B when there is no ambiguity): We adopt a local state approach, where the local material state is described by the deformation gradient F(X, t): where GL þ (R, 3) is the space of second-order tensors on R 3 with positive determinant, and temperature T(X, t). Since we want to consider dissipative behaviors, a set of internal variables Z is added to external state variables F and T. The exact nature of Z depends on the type of constitutive behavior, as will be detailed below (Sections 4 and 5). Following the framework of GSM (Halphen & Nguyen, 1975), we assume the existence of a Helmholtz free-energy density potential W(F, T, Z ) and a dissipation pseudo-potential D( _ F, _ Z; F, T, Z) (both defined per unit undeformed volume). The last three arguments of the dissipation pseudo-potential, separated from the main ones by a semi-colon, denote the possible dependence of this function on the current state of the material point. We additionally assume that forces conjugate to state variables are additively decomposed in an equilibrium part, derived 6 from the free energy, and a dissipative part, derived from the dissipation pseudo-potential (Ziegler, 1977): and where P is the first Piola-Kirchhoff (or Piola) stress tensor, conjugate to F, and Y the forces conjugate to internal variables Z. What differentiates external and internal variables is that the latter should produce no net work, i.e.
As a consequence, one must have that This equation provides the evolution law for internal variables. Finally, we will assume in this work that the local thermal equilibrium is always verified, such that the specific entropy is given by where r 0 is the density in the reference configuration. In the following, it will prove useful to introduce the internal energy density potential U (F, , Z ): such that Note that objectivity (material frame indifference) imposes that state functions such as the free or internal energy be invariant under rotations in the spatial configuration: where we have assumed that all internal variables were of Lagrangian nature. This condition can be verified by enforcing that the free (resp. internal) energy depends on F only through the right Cauchy-Green stretch tensor C ¼ F T F. The Piola stress is then given by where S is the second Piola-Kirchhoff (symmetric) stress tensor.

Balance equations
We now briefly recall local balance equations corresponding to conservation principles (mass, momentum, and energy). In Lagrangian description, mass conservation yields where r is the mass density in the deformed configuration. Linear momentum conservation yields where b denotes applied bulk forces (per unit mass), while angular momentum conservation yields Note that, given Eq. (17), this will be automatically verified for objective (frame invariant) constitutive models. Conservation of energy yields where H is the nominal (Lagrangian) heat flux vector and Q the applied bulk heat source (per unit mass). We will define the internal dissipation as Then, using the definition (13) of entropy, we can obtain the heat equation: 8 where the heat capacity C, defined by is in general a function of the current state. In the following, we will refer to the first two terms in the right-hand side of Eq. (22) as entropic heat source terms, in opposition to the dissipative heat source term D int . The Clausius-Duhem inequality (second principle of thermodynamics) then writes where _ G denotes the net entropy production rate. Choosing a dissipation pseudo-potential D( _ F, _ Z; F, T, Z) which is convex with respect to _ F and _ Z, nonnegative and such that D(0, 0; F, T, Z) ¼ 0 will ensure that D int ! 0. Following Biot (1958), we then introduce a conduction potential w(G; F, T, Z ), such that Assuming convexity and nonnegativeness of w, together with the condition w(0; F, T, Z ) ¼ 0, then ensures that Clausius-Duhem inequality will always be verified. In order to provide Eulerian versions of the above balance equations, let us first recall some results linked to finite strain kinematics. The velocity gradient is defined as L ¼ _ F F À1 , and its symmetric part defines the Eulerian strain rate The Cauchy stress tensor s is given by The linear momentum conservation equation then becomes where we have the used mass conservation equation, while conservation of energy becomes with h ¼ J À1 F H is the (Eulerian) heat flux vector.

Variational updates
In the isothermal setting, Ortiz and Stainier (1999) showed how a wide class of constitutive models for dissipative solids could be recast under the form of variational principles. But when the equilibrium and heat conduction problems for general dissipative solids are combined, the resulting coupled problem lacks an obvious variational structure. This lack of variational structure reveals itself upon linearization of the coupled problem, which results in a nonsymmetric operator. This essential difficulty accounts for the lack of variational formulations of the coupled thermo-mechanical problem for general dissipative solids. However, Yang et al. (2006) showed that an integrating factor exists which delivers the sought-for variational structure. This integrating factor hinges critically on a careful distinction between two types of temperature: an equilibrium (or internal) temperature, which follows as a state variable; and an external temperature, which equals the equilibrium temperature at equilibrium.

Local evolution problem
At a given material point, evolution of internal variables is ruled by Eq. (12). Taking into account the convexity of D w.r.t. _ Z, this equation can be interpreted as the stationarity condition associated to the following variational principle: In the above expression, Y(F, , Z ) is the equilibrium or internal temperature, as defined by Eq. (15): while T will be called external temperature. Local thermal equilibrium imposes that Y ¼ T, which corresponds to stationarity of D w.r.t. _ : Note that we could have introduced local thermal dissipation (e.g. transient heat transfer at a finer scale than that of the representative volume element (RVE) associated to each material point) through a more general dissipation function D( _ F, _ , _ Z), but we will not pursue this possibility here. The apparent complexity of the function D( _ F, _ , _ Z, T; F, , Z ) is in practice motivated by the following properties. Derivatives of the function w.r.t. _ F and T are, respectively, given by: and Then, if we define the effective rate potential D eff as we can write that 11 Thus, provided that T ¼ Y(F, , Z), it is then seen that function D eff ( _ F, _ , T ) plays the role of a rate potential for the stress P and the effective local entropy rate (i.e. local entropy rate corrected by internal dissipation D int ). For example, looking at heat equation (21), we see that local adiabatic behavior would then be characterized by and i.e. D ad acts as a rate potential for the stress tensor P under local adiabatic conditions.

Local time-discrete constitutive problem
Under most circumstances, nonlinear problems involving historydependent behavior are solved numerically by incremental methods. We thus proceed by considering a discrete time increment [t 0 , t]: the local material state at time t 0 ({F 0 , T 0 , Z 0 } or equivalently {F 0 , 0 , Z 0 }) is assumed to be completely known and, in a first step, we would like to compute the internal material state Z associated to a given external material state {F, T} (or equivalently {F, }) at time t. Starting from the variational formulation of the rate problem (36), we thus seek to define an incremental function I(F, T, Z) which approximates the integral of function D over the time increment. Noting that _ where The last term on the right-hand side, between brackets, denotes an average value of the dissipation function over the time increment, the expression of which is discussed below. Note already that rates have been approximated by a first-order finite difference (alternative expressions are possible), while factors T/Y have been replaced by T/T 0 .
The incremental variational update takes the form of the following minimization problem: The stationarity equation corresponding to this variational problem is given by In addition, we will have that Consistency of the incremental update thus requires that but also that and lim dt!0 As discussed in Stainier (2011), the most intuitive expression where {F, T, Z} a ¼ (1 À a){F 0 , T 0 , Z 0 } þ a{F, T, Z} (with algorithmic parameter a 2 [0, 1]) does not verify all the above consistency conditions. An alternative expression which does verify all the above conditions was also provided in Stainier (2011): One can see that in the case a ¼ 0, or in the case of no parametric dependence of D on the current state, the above expression reduces to In that specific case, the internal variable update becomes and we also have We can interpret these results as providing approximate incremental expressions for the viscous stress, the thermodynamic forces conjugate to internal variables: and the dissipation (average value over the time step) These quantities are not state functions, and thus can only be approximated in an incremental setting. Incremental approximations to the stress and effective entropy increase are then given by In the more general case (a 6 ¼ 0) additional terms will appear in the incremental approximations of P d , Y d , and D int , but the resulting expressions will remain consistent, i.e. will tend toward continuous values as dt ! 0. A detailed study of the influence of algorithmic parameter a on precision and convergence properties of the variational update in the case of isothermal elasto-visco-plasticity can be found in . Among other results, it is shown that it is always possible to find a value of parameter a which yields incremental approximations identical to what would be obtained using a standard backward Euler scheme on continuous rate equations. This value strongly depends on the model used, of course, and will also vary with actual loading conditions (e.g. strain rate). Note that, by consistency, the effect of a vanishes as the time increment size decreases.

Variational boundary-value problem
The variational treatment of the local constitutive problem, as detailed in the previous section, allows to construct effective potentials of strain (or strain rate) and temperature (and entropy rate) describing a formally thermo-elastic behavior. The resulting constitutive relations nonetheless include the effect of internal variables (and the associated internal dissipation). The coupled thermo-mechanical boundary-value problem can then in turn be formulated variationally using these effective potentials.

Rate problem
Let us now consider the quasi-static boundary value problem consisting in determining fields _ w, _ , _ Z, and T on B 0 at a given time t (the current state {F, , Z} is assumed to be known), verifying the differential equations together with local constitutive equation (12) and boundary conditions: , and H(t) are imposed motion, traction vector, temperature, and normal heat flux, respectively.
As shown in Yang et al. (2006), this boundary value problem can be restated in variational form. To this end, let us define the functional where we have omitted to explicitly indicate the sufficiently clear dependence of D and F on the current state for brevity. Then balance equations (61) correspond to stationarity conditions for F: while stationarity condition w.r.t. _ yields T ¼ Y, as seen previously. Applying Green-Ostrogradsky's theorem to the above stationarity equations, and replacing Y by T, yields local equations (61) and associated boundary conditions.
As discussed in Yang et al. (2006), we may focus our attention on the problem of determining extremal points of F, expecting these to correspond to stable solutions to the boundary-value problem. In a number of cases, extremal points will actually correspond to a saddle point of the functional: In the case of thermo-elasticity (i.e. Z ¼ ;, D 0), this was demonstrated under assumptions of pure Dirichlet boundary conditions and sufficient regularity of solution fields (see Yang et al., 2006). More precisely, in that case, we can write If we exclude viscous stresses [i.e. D ¼ D( _ Z; F, T, Z )], then the effective rate potential D eff is concave in T [this results from the implicitly assumed convexity of D( _ Z )]. The boundary-value problem is then formally identical to that of thermo-elasticity, and the previous result still holds. In the presence of viscous stresses (for example, Kelvin-Voigt thermo-visco-elasticity), the analysis is more subtle, and additional studies still have to be conducted for that general case.

Incremental boundary-value problem
We now come back to the incremental problem. As explained earlier, we consider a discrete time increment [t 0 , t], for which we consider the fields {w 0 , T 0 , Z 0 } as known. We have seen in Section 3.1.2 how we can compute updated internal variables Z at each material point, provided external variables {F ¼ -0 w, T} are given. External fields {w, T} can themselves be computed as optimizers of the following functional: where the average conduction dissipation function hwi can be treated in a similar fashion as dissipation potential D occurring in the local incremental problem. The incremental boundary-value problem thus takes a variational form: In many cases (we will consider specific examples in the following), the incremental potential W(F, T) shows up to be convex in F and concave in T (i.e. unless limits of material stability have been attained). We can then characterize the solution fields as a saddle point of the incremental functional: By repeating this optimization problem, taking results of the previous time step as initial conditions for the current increment, one can thus compute the evolution of dissipative thermo-mechanical systems. It is important to note that the functional I which is to be optimized changes at each time step. Note that Canadija and Mosler (2011) have proposed an alternative incremental variational formulation, which is written in terms of entropy and internal temperature, in keeping closer to the time-continuous formulation presented earlier.

Mixed thermal boundary conditions
The variational formulation of coupled thermo-mechanical boundary-value problems described earlier only includes pure Dirichlet or Neumann boundary conditions. More complex boundary conditions, such as contact or heat convection, frequently occur. We will not discuss contact here, since it is a topic in itself. Variational approaches to contact have for example been proposed in Johnson, Ortiz, and Leyendecker (2012), Kane, Repetto, Ortiz, and Marsden (1999), and Pandolfi, Kane, Marsden, and Ortiz (2002), and the reader is referred to these works and references therein. With respect to the heat transfer problem, we can easily include mixed boundary conditions, where the imposed heat flux depends on the temperature at the boundary. This kind of boundary condition can for example represent heat exchange with a surrounding fluid by (forced or natural) convection. It can also be used to model heat transfer to another contacting solid.
For illustration purposes, we will consider the simplest heat exchange model, where the heat flux at a boundary is proportional to the temperature difference between the boundary and its immediate environment: where @ c B is the part of the (deformed) boundary where heat exchange conditions are effective, T ext is the temperature of the external fluid or solid, h is a given heat exchange coefficient, and n is the outward normal to the (deformed) boundary. Such boundary conditions can be included in the variational formulation by adding an additional term F c to functional (62). This term is given by where J S is the ratio of a surface element in the deformed configuration to the same surface element in the reference configuration: In this modified variational principle, it is of course understood that In the time-discrete setting, the variational principle describing the incremental coupled thermo-mechanical boundary-value problem can similarly be modified by the following function: where J S a denotes the surface Jacobian evaluated at time t a ¼ (1 À a)t 0 þ at. The variational principle is then given by It is interesting to note that the inclusion of the term I c introduces an additional source of coupling between displacement and temperature. Indeed the modification of the boundary area caused by deformation has a direct effect on the total heat flux exchanged with the environment. Because of the variational nature of the formulation, this also implies that an additional term will appear in the mechanical balance equations (variations of F À F c with respect to w). This term is purely numerical, but ensures the symmetry of the mathematical formulation. It will vanish when dt ! 0, as consistency requires. In practice, this effect can be avoided by choosing a ¼ 0 in expression (73). The effect of deformation on the exchanged heat flux will then be delayed by one time step.

Dynamics
The variational formulation (68) accounts for rate-dependent behavior, including transient thermal effects. It formally corresponds to a quasistationary problem, where heat capacity terms are treated as rate-dependent thermo-mechanical behavior. Inertia terms can be accounted for in the time-discrete setting, but at the cost of directly embedding a specific time discretization within the formulation. Indeed, we can extend the approach initially proposed by Radovitzky and Ortiz (1999) in the isothermal context, leading to the following modified functional: The incremental boundaryvalue problem then takes the variational form When combined with the following update rule for accelerations and velocities: the stationarity condition of I with respect to w indeed yields a discrete conservation of momentum equation corresponding to the classical Newmark integration scheme of dynamics (Hughes, 2000). The classical Newmark scheme is typically used with parameters (b ¼ 0.25, g ¼ 0.5) in implicit dynamics (ensuring unconditional stability in linear elasto-dynamics), or with parameters (b ¼ 0.0, g ¼ 0.5) in explicit dynamics (with only a trivial matrix inversion when using lumped mass matrix). Clearly, the latter case cannot be implemented within the proposed variational framework, which thus appears to be intrinsically implicit.
Note that variational constitutive updates also proved very useful in deriving energy and momentum-conserving time integration schemes (Noels, Stainier, & Ponthot, 2006, since they allow to work with an effective incremental potential formally identical to hyper-elasticity.

Linearization
For a wide range of engineering cases, the problem can be simplified by considering that displacements and displacement gradients are small. The difference between the initial and deformed configurations can then be neglected, and the relevant strain measure is the engineering strain: where I is the identity second-order tensor. All the above variational principles still hold, provided that the deformation gradient is replaced by and Piola stress tensor P by Cauchy stress tensor s. Similarly, if temperature variations are small (about a reference temperature T r ), i.e.: the problem can be simplified. The generalized gradient is then given by while coefficients occurring in the time discrete variational formulation become A consistent approximation to the average dissipation function (eq. 3.21) is then obtained by replacing factors as follows: 21 Finally, the "linearized" balance equations become: Note that the assumptions of small displacements and small temperature variations are independent and not necessarily linked. Although the above approximations lead to linear balance equations in the uncoupled case, this is not true in general in the coupled case. Hence, the term linearization should be considered only as meaning linearization of strain and temperature gradient measures.

THERMO-VISCO-ELASTICITY
In this section, we look in more detail at variational formulations for thermo-visco-elastic constitutive models. We start by the simplest model: small-strains Kelvin-Voigt visco-elasticity. This model does not require internal variables, and its relative simplicity allows to develop the details of the variational treatment of thermo-mechanical coupling. Afterwards, we look at more general Maxwell models, involving internal variables, and repeat the previous analysis in the context of finite strains. Note that in this section, we mostly consider time-discrete variational constitutive updates, although time-continuous variational formulations are also available.

Linearized kinematics
We will first briefly review the basic models of thermo-visco-elasticity under linearized kinematics assumptions, and their variational formulation.

Kelvin-Voigt model
The simplest visco-elasticity model can be described by combining a linear (thermo-) elasticity free-energy potential: where C(y) is a (temperature-dependent) fourth-order elasticity tensor, a a second-order tensor describing (possibly anisotropic) thermal expansion, C the heat capacity (per unit undeformed volume), and a linear viscosity dissipation function: where C u y ð Þ is a (temperature-dependent) fourth-order viscosity tensor. The entropy is then given by the following expression: where we have a heat capacity term, a thermo-elastic term, and two terms linked to thermal softening of elastic moduli. Note that, although the individual mechanical and thermal behaviors are linear, the coupled thermomechanical model will be nonlinear because of the presence of dissipative terms (which will show in the effective entropy increase below). In this case, there are no internal variables, and the incremental energy potentials are given by Using the consistent approximation of incremental dissipation functions proposed in Stainier (2011), the last term becomes Considering the specific quadratic form we chose for D in this setting, this expression reduces to The resulting incremental expressions for stress and net entropy increase are then given by which are consistent with continuous expressions when dt ! 0.

Generalized Maxwell model
We can now consider slightly more complex visco-elasticity models, using internal variables. In particular, generalized Maxwell visco-elasticity models consist of an elastic branch in parallel with one or more visco-elastic branches. For the purpose of illustrating the variational formulation in that case, we will limit ourselves to a single visco-elastic branch. Such a model is then described by the following thermodynamic potentials: where « e ¼ « À « u , and Typically, W (0) can take a form similar to Eq. (86), while W (1) can take a simple quadratic form: where C (1) (y) is another (temperature-dependent) fourth-order elasticity tensor. Note that one could also add a thermo-elastic term to W (1) . In the following, we will restrict ourselves to temperature-independent viscous moduli, for the sake of clarity and conciseness. The treatment of such temperature dependence follows the same lines as shown in the Kelvin-Voigt case. The incremental energy potential now takes the following form: 24 Stationarity condition for the viscous strain yields: Note that in the case of multiple Maxwell branches, the viscous strain in each of these branches can be optimized independently. Thermodynamic forces conjugate to strain and temperature through the incremental potential are then given by: Àr where s (0) («, y) ¼ @ « W (0) («, y) and « u is given by Eq. (98). Once again, these results are consistent with continuous expressions when dt ! 0. The computation of material tangents is straightforward from the previous expressions and will not be detailed here. More sophisticated models can be constructed by multiplying the number of Maxwell branches and/or adding a Kelvin branch in parallel.

Finite thermo-visco-elasticity
We are now ready to repeat the previous analysis in the context of finite transformations, following the presentation of Fancello et al. (2006). It will appear that a key aspect for an efficient treatment of constitutive updates in this context, at least for isotropic materials, is the use of a spectral description.

Kelvin-Voigt model
The above thermo-visco-elasticity models and their variational formulation can be extended to finite strains. For example, one can consider the following general isotropic elastic free energy potential: where J ¼ det[F],Ĉ ¼F TF ,F ¼ J À 1 3 F, f( J, T) an energy function linked to volume changes (possibly including thermo-elastic effects), W e a strain-energy function with temperature-dependent elastic moduli, and W t a thermal stored energy (heat capacity) function. A typical example of volumetric function is where K is the bulk modulus and a the coefficient of thermal expansion, while for heat capacity, a typical example is where C is the heat capacity (per unit undeformed volume). For the strainenergy potential W e , many possibilities are available, among which neo-Hookean, Mooney-Rivlin, or the general Ogden potential (see for example Holzapfel (2000) for other possibilities).
The above free-energy potential can be combined with a general dissipation potential of the form This specific form implies that no viscous stress is generated by pure rotations (objectivity). The simplest case would then be given by Note that in an incremental context, one must be attentive to preserving objectivity. For example, if approximating L by dF F a À1 , rigid body rotations will generate spurious strain rate. Suitable approximations can be built by which clearly becomes equal to zero in case of rigid-body motions. An alternative approximation, which preserves isochoricity (i.e. if det F ¼ det F 0 , then tr D ¼ 0), is given by: The incremental potential energy then takes the form where a consistent approximation to the average dissipation function is provided by Considering the simple quadratic form earlier, it yields where D is given by one of the above incremental approximations. The stress tensor is then given by while the net effective entropy variation is given by It is relatively straightforward to verify that these expressions are consistent with the continuous formulation when dt ! 0.

Generalized Maxwell model
Following closely the presentation adopted in the linearized kinematics section, we now introduce a finite strains version of the generalized Maxwell model considered before. A presentation of general finite strains visco-elastic models, including a Kelvin-Voigt branch and several Maxwell branches can be found in Fancello et al. (2006), yet we will limit ourselves to a single branch here. We thus consider the following free-energy thermodynamic potential: where we adopt the multiplicative decomposition initially proposed by Sidoroff (1974) for finite visco-elasticity: We recall that free-energy potentials should depend on C ¼ F T F and C e ¼ F eT F e in order to ensure objectivity. Viscous deformations are assumed to produce no rotation: where D u is a viscous strain-rate symmetric tensor. Viscous strains are typically assumed to be isochoric, in which case tr[D u ] ¼ 0. In the simplest case (one Maxwell branch in parallel with an elastic branch), the dissipation function then takes the form: In order to avoid lengthy developments, we will assume in the following that there is no parametric dependence of f (1) on T. The incremental energy potential then takes the form: where one must provide an incremental update rule for F u . The exponential update formula of Weber and Anand (1990), initially proposed in the framework of finite elasto-plasticity: provides the advantage of preserving the isochoric nature of viscous deformation (if enforced) and will thus be used here. We can then write where C tr e ¼ F 0 uÀT C 0 e F 0 uÀ1 is a predictor (trial) elastic Cauchy-Green tensor.
Considering the specific case of isotropic materials, for which thermodynamic potentials can be represented as functions of principal values of their tensorial arguments: where {c 1 e , c 2 e , c 3 e } (resp. {d 1 u , d 2 u , d 3 u }) are the principal values of tensor C e (resp. D u ), the optimization problem of the constitutive update can be reduced to a simpler form. Indeed the ansatz that D u , C tr e , and thus C e , are co-linear (i.e. share the same principal directions) can then be verified a posteriori (see Fancello et al., 2006), yielding log l e I À Á ¼ log l tr I À Á À dtd u where l e I ¼ ffiffiffi ffi c e I p and l I tr are the principal stretches associated to the elastic and trial Cauchy-Green tensors. It proves convenient in practice to use the logarithmic principal strains e I e ¼ log(l I e ) and e I tr ¼ log(l I tr ). The optimization problem now reduces to , e e 2 , e e 3 , T where e I e ¼ e I tr À dt d I u , and where an isochoricity constraint can be added: The problem is then reduced to solving three (or four in the isochoric case) nonlinear scalar stationarity equations. The net effective entropy variation is then given by which is consistent with the continuous formulation when dt ! 0.

Viscous fluids
Newtonian fluids are described by Navier-Poisson constitutive equations: where s is Cauchy stress tensor, D ¼ (L þ L T )/2 is the strain rate tensor, and p(J, T) is the hydrostatic pressure, related to the density r and temperature T through the equation of state (EOS). It is usually assumed that, either the fluid flow is incompressible (tr[D] ¼ 0), either k ¼ 0 (Stokes condition), such that the pressure is always equal to the hydrostatic pressure: tr[s] ¼ Àp.
The Navier-Poisson constitutive equations can be seen as a particular case of a finite-strain Kelvin-Voigt visco-elasticity model and can thus easily be put under variational form. This is obtained by considering a purely volumic Helmholtz free energy: 29 where J ¼ det F, and a viscous dissipation potential of the form: It is easily verified that this yields Navier-Poisson equations, with k ¼ 0 and where we recall that Jr ¼ r 0 . Finally, note that most non-Newtonian viscous fluid models could also be formulated in the current variational framework.

THERMO-ELASTO-VISCO-PLASTICITY
Elasto-(visco-)plastic behavior is characterized by the existence of a domain in stress space within which the material behaves elastically. On the boundary of this domain (and outside in the case of visco-plasticity), plastic deformation can occur. Traditionally, plasticity models are thus described by the definition of a function of stress defining the elastic domain (Lubliner, 1990). In the case of nonassociated plasticity, another function of stresses can be defined, from which the plastic strain rate is computed. In the current setting, we will not make explicit use of such yield and flow functions, but instead introduce a flow rule and a dissipation function (which, in the rate-independent case, is actually the dual to the indicator function of the admissible domain). The latter approach is more kinematic in nature, but it is in the end equivalent to the more classical stress-based approach. It nonetheless offers the advantage of being well adapted to the variational formulation presented in this paper.

Constitutive modeling
In (poly-)crystalline materials, plastic deformation is due to the motion of dislocations along certain slip directions on specific slip planes (the combination of a particular slip direction and slip plane will be referred to as a slip system). Following Rice (1971), we thus adopt a flow rule of the form where g (k) is the slip strain, and s (k) , and m (k) are orthogonal unit vectors defining the slip direction and slip-plane normal corresponding to slip system k. The collection g of slip strains may be regarded as a subset of the internal variable set Z. A zero value of a slip rate _ g k ð Þ signifies that the corresponding slip system is inactive. The flow rule (129) allows for multiple slip, i.e. for simultaneous activity on more than one system over a region of the crystal. The vectors {s (k) , m (k) } remain constant throughout the deformation and are determined by crystallography. In order to account for nonmonotonous loading paths, it is common to consider dislocation glide in þs (k) and Às (k) directions as occurring in separate systems, and adding the constraint _ g k ð Þ ! 0 (for each of the duplicate systems) in order to preserve consistency of the formulation. Plastic deformations leave the crystal lattice undistorted and unrotated, and, consequently, induce no long-range stresses. Some degree of lattice distortion F e , or elastic deformation, may also be expected in general. One therefore has, locally, This multiplicative elastic-plastic kinematics was first suggested by Lee (1969), and further developed and used by many others. A classical assumption is to consider that the elastic behavior is unaffected by other internal processes (dislocations in this case), yielding the following expression for the free energy: where W e is the elastically stored energy (recoverable), W p is the plastically stored energy (not directly recoverable), for example under the form of dislocation microstructures, and W t the thermally stored energy (heat capacity). The dependence of the elastic energy on C e (instead of F e ) ensures objectivity, as explained earlier. The rate of free energy can then be written as where t ¼ {t (k) } is the collection of resolved shear stresses: while g ¼ {g (k) } is the collection of yield resolved shear stresses linked to plastic energy storage mechanisms g k ð Þ ¼ @ g k ð Þ W p .

31
The kinematic of plastic slip can be modeled through a dissipation function D _ g; g, T ð Þ. Typical expressions include power-law type formulas such as where Y ¼ {Y (k) } is a collection of dissipative yield resolved shear stresses and _ g k ð Þ 0 are reference slip rates. The exponent m 2 [1, þ1] controls ratedependency effects. In particular, rate-independent behavior can be recovered at the limit when m ! þ1. This power-law expression is indeed often used a way to regularize rate-independent models. Complex hardening models (e.g. including latent hardening) can be included in the generic expression Y (k) (g, T ) (see for example Stainier, Cuitiño, & Ortiz, 2002). More complex ratedependency models, e.g. based on thermal activation theories, can also be formulated under the form of a dissipation function (again, see Stainier et al., 2002). The rate problem of crystal plasticity can then be formulated variationally by introducing functional and the minimization problem Considering for example the case of the power-law dissipation function given in Eq. (134), stationarity conditions yield where we have accounted for the fact that posterior stationarity condition with respect to _ will yield T ¼ Y.

Incremental update
In the time-discrete setting, the local variational constitutive update for crystal plasticity takes the following form: In the above expression, a relation must be provided between F p , dg, and F 0 p , which is consistent with the continuous flow rule (129). The exponential update (Weber & Anand, 1990), first used in the context of crystal plasticity by Miehe (1996): provides such an expression, with the additional advantage that it preserves isochoricity of the plastic flow. Consistent approximation of the average dissipation function can be constructed following the general formula (51), but for the remainder of this section, we will consider the simpler case without parametric dependence of the dissipation function on temperature: The constitutive update problem then reduces to the following constrained optimization problem: and F p given by Eq. (139). Obtaining a numerical solution to this optimization problem may not be a trivial task, especially in the presence of complex (latent) hardening phenomena. Without going into the details, let us simply say that, from our experience, constrained optimization methods such as proposed by Bertsekas (1995), especially two-metric projection methods, proved quite robust in this context.

Macroscopic plasticity
In this section, we will mainly focus on J 2 (von Mises) plasticity, since it is well representative of macroscopical plasticity models, and although without doubt the most widely used. In addition, its specific structure allows for pushing analytic expressions farther, including in the large deformation setting, yielding efficient numerical algorithms, which we will relate to the classical radial return of Wilkins (1964).

Linear kinematics
Under small strains assumptions, the total deformation is additively decomposed into an elastic strain « e and a plastic strain « p : A plastic flow rule corresponding to von Mises-type plasticity can be written as follows: where _ e p is the amplitude of the plastic strain rate and M its direction. The symmetric tensor M must of course be normalized, and the particular choice ensures that e p corresponds to the cumulated equivalent plastic strain. In addition, we will require that plastic deformation be isochoric: Relation (143) can be seen as a reparametrization of plastic strain rate. For most polycrystalline metals, it is generally accepted that work hardening does not modify the elastic behavior of the material, leading to consider the following expression for free energy: where W e is the elastically stored energy (recoverable), W p is the plastically stored energy (not recoverable), and W t the thermally stored energy (heat capacity). Accounting for Eq. (143), we can write the rate of free energy as where s ¼ @ « W is the (Cauchy) stress, s c ¼ @ « p W p is the back-stress (kinematic hardening), and g ¼ @ e p W p is a yield stress associated to energy storage mechanisms (such as dislocation microstructures). From the above, we see that quantity y can be considered as conjugate to the cumulated equivalent plastic strain e p , when accounting for the flow rule. Following the formalism of GSM, we will then describe the relation between y and _ e p through the definition of a convex dissipation potential D _ e p ; e p , T À Á , or its conjugate D * y; e p , T ð Þ.

Rate problem of visco-plasticity
and the effective rate potential is defined as Note that the minimization with respect to M can be related to the principle of maximal dissipation: Accounting for constraints on M, the optimal plastic flow direction is given by where s ¼ dev[s À s c ], which corresponds to the normal direction to von Mises yield criterion. Using this result, we can rewrite conjugate force y as where s eq is von Mises equivalent stress (accounting for backstress s c ). Dissipation potential D _ e p ; e p , T À Á takes the general form with f a convex function such that f 0; e p , T ð Þ ¼ 0 and @ _ e p f 0; e p , T ð Þ ¼ s y e p , T ð Þ ! 0. Minimization of D with respect to _ e p then yields the following result: where we have accounted for the fact that stationarity with respect to _ will yield T ¼ Y. Note that in the rate-independent case, corresponding to a function f linear in _ e p , we have @ _ e p f _ e p ; e p , T À Á ¼ s y e p , T ð Þ, and thus necessarily s eq À g ¼ s y . We thus recover the classical results of von Mises elasto-visco-plasticity.

Constitutive updates
In most applications, the problem of interest is to compute the time evolution of stress and plastic deformations. The local time-discrete constitutive problem then takes the following form: where Consistent expressions for the average dissipation function hDi have been proposed in Stainier (2011) and will not be detailed here. In order to avoid cluttering the presentation, we will instead consider the special case where there is no parametric dependence of D on temperature, in which case we can take Stationarity condition for M (including constraints) yields relation (151) again, where s ¼ dev[s À s c ] depends on e p . Since dissipation function D _ e p À Á is not regular, a practical approach to finding the infimum with respect to e p consists at first evaluating the gradient of the incremental energy functional at the singularity point: where s tr eq is the (trial) equivalent von Mises stress computed at material state «, y, « p 0 , e p 0 f g, and where we have used the stationarity condition on M. Considering the convexity of the incremental energy (w.r.t. e p ), we know that if this gradient is positive (i.e. s eq tr g 0 þ e s y, tr ), then the optimum is at d e p ¼ 0, and the increment is elastic, while otherwise the optimum is d e p > 0. We thus recover a predictor-corrector scheme, such as classically used in computational plasticity.
The von Mises criterion mostly makes sense in the context of isotropic elasticity, in which case we can write the elastic free energy as where f can be a quadratic function of tr[« e ] in the simplest case, or include a linear term if accounting for (isotropic) thermo-elasticity. If, in addition, we consider the case of pure isotropic hardening (s c 0), we can then write where « tr e ¼ « À « 0 p and s tr ¼ 2m dev[« tr e ]. Considering stationarity condition (151), an immediate consequence is that s, s tr , and M are aligned: Thus, in the case of purely isotropic elasto-visco-plasticity, the variational formulation is equivalent to the classical radial return algorithm (Wilkins, 1964): The equivalent von Mises stress is then given by s eq ¼ s eq tr À 3m y ð Þd e p ð163Þ and the stationarity condition for e p is Às eq tr þ 3m y ð Þd e p þ g e p , y ð Þ þ s y e p , y ð Þ ¼ 0 where is a consistent incremental approximation of the dissipative part of the yield stress. A study of the effect of the choice of algorithmic parameter a in the context of elasto-visco-plasticity can be found in .

Finite strains
Under the general finite strains regime, several options are possible to introduce the notions of elastic and plastic strains. Here, we will adopt the multiplicative decomposition, first suggested by Lee (1969) and commonly used by many other authors [see for example Simo and Hughes (1998)]: A plastic flow rule corresponding to von Mises-type plasticity can then be written as follows: with M a symmetric tensor such that as before (linear kinematics), i.e. we assume that plastic deformation is isochoric and that it does not generate any rotation (thus defining a unique intermediate configuration). Just as in previous section, we will consider that work hardening does not modify the (hyper-)elastic behavior of the material, yielding where W e is the elastically stored energy (recoverable), for which we have accounted for objectivity through the use of C e , W p is the plastically stored 38 energy (not recoverable), and W t the thermally stored energy (heat capacity). Accounting for Eq. (167), the rate of free energy can be expressed as where P ¼ @ F W is the Piola stress tensor, Y c ¼ @ F p W p is a back-stress tensor (linked to kinematic hardening), and g ¼ @ e p W p is a yield stress associated to energy storage mechanisms, as introduced in the linearized kinematics section. From this relation, we can see that quantity Y can be considered as conjugate to the cumulated plastic strain e p , when accounting for the specific flow rule chosen. This quantity can be rewritten as where T ¼ F eT PF pT is Mandel stress tensor, and T c ¼ Y c F pT is the associated backstress tensor.

Rate problem of finite visco-plasticity
and the effective rate potential is defined as Accounting for constraints on M, the optimal plastic flow direction is given by where s ¼ dev[T À T c ], which corresponds to the normal direction to von Mises yield criterion, written in terms of Mandel stress tensor. Using this result, we can rewrite conjugate force Y as where s eq is von Mises equivalent stress (computed in the intermediate configuration from Mandel stress tensor, and accounting for backstress T c ). We thus recover a set of expressions formally similar to those obtained in the setting of linearized kinematics, where Mandel stress tensor (related to the intermediate configuration) is used instead of Cauchy stress tensor. This results naturally from the variational formulation, and this is where it mainly differs from most other formulations, which often will work with Cauchy stress in the deformed (spatial) configuration.

Exponential mapping
In the time-discrete setting, the local variational constitutive update problem takes the following form: where one must provide a time-discrete update rule for plastic deformation F p . The exponential update formula (Weber & Anand, 1990): provides the advantage of preserving the isochoric nature of plastic deformation and will thus be used here. For reasons similar to those provided in the linearized kinematics case, we will consider no direct parametric dependence on temperature in the dissipation function D, in which case a consistent approximation to the average dissipation function over the time increment is given by:

Hencky hyperelasticity
In finite strains, many choices are possible for the elastic free energy, even under the assumption of an isotropic material. As was shown in Ortiz and Stainier (1999), the specific choice of a Hencky hyperelastic free energy allows to perform analytically the optimization with respect to M in the variational principle above. Hencky hyperelastic free energy is obtained by considering a quadratic potential in logarithmic (natural) strains: with Considering the exponential update formula (177) and the ansatz that C tr e ¼ F 0 pÀT C F 0 pÀ1 and M are colinear, we can write framework (Bleier & Mosler, 2012;Mosler, 2010;Mosler & Bruhns, 2009). Without going into the details, let us simply say that this approach is based on a reparametrization of the plastic flow rule in terms of a normalized stress tensor: flow rules can then be obtained by deriving classical yield functions (provided that they are homogeneous in stress). The only downside of this generic approach is that it may sometimes be more difficult to integrate with a detailed description of thermo-mechanical coupling effects (partition of yield stress into stored and dissipative parts, see next section).

Heat generated by visco-plastic dissipation
Let us come back to the rate problem of finite thermo-visco-plasticity. If we look at derivatives of the effective rate potential, we obtain: If we take into account stationarity of D eff with respect to _ : and comparing with Eq. (38), we obtain the following expression for the internal dissipation where we have used the definition (171) and flow rule (167). A long-standing issue in metal plasticity has been to estimate the internal dissipation from the measure of plastic power TÁD p . In their pioneering work, Taylor and Quinney (1937) experimentally measured dissipation amounting to somewhere between 90% of plastic work (they actually used quantities integrated over time) [see also the later review by Titchener and Bever (1958)]. On the basis of these results, many (most) contemporary authors use the following formula to evaluate visco-plastic dissipation in their numerical simulations: with b considered as a constant material parameter (typically chosen as b ¼ 0.9). This formula remains widely used today, despite more recent experimental work (see for example Rittel, Bhattacharyya, Poon, Zhao, & Ravichandran, 2007) clearly showing that the ratio of dissipation to visco-plastic power varies strongly with plastic strain and strain rate, and with temperature. Some models have been proposed (see for example Zehnder, 1991), aiming at providing an expression describing evolution of coefficient b. In the variational framework described here, the coefficient b does not appear explicitly (although it can be computed a posteriori). Instead the ratio of dissipation to total viscoplastic power directly derives from the choice of free energy and dissipation potentials. In this way, it is possible to take into account complex evolutions of b, such as result from combined rate-dependent visco-plastic behavior and general temperature dependence (including of stored energy, for example associated to recrystallization). As illustrated in Stainier and Ortiz (2010), such cases cannot be modeled by approaches such as that of Rosakis et al. (2000), yet are naturally and implicitly accounted for in the variational approach.

NUMERICAL APPROXIMATION METHODS
6.1. Variational finite element approximations 6.1.1 Standard galerkin formulation Variational boundary-value problems such as Eqs. (69) or (74) quite naturally call for finite element numerical solution methods. For example, the space V of admissible thermo-mechanical configurations (at time t) is given by Then, approximate finite element solutions to problems (69) or (74) can be derived by a Ritz-Galerkin approach. To this end, consider the admissible subspace V h , built on a given discretization T h of B 0 : where x a ¼ w(X a , t) are the positions of the N nodes nodes in the current (deformed) configuration and T a ¼ T(X a , t) are the temperatures at these same nodes. As usual, the material shape functions N a (X) are defined on elements connected to node a and set to zero elsewhere. Index h refers to the mesh T h supporting the shape functions. Note that we have chosen here to use the same shape functions for displacement and temperature fields, for the sake of simplicity, but that it would also be possible to choose different shape functions and different meshes for these two unknown fields. The discrete deformation and temperature gradients are given by and the variational principle (69) becomes where w 0 and T 0 are computed by interpolation from nodal values of x a (t 0 ) and T a (t 0 ), while the internal variables Z 0 are typically stored at integration points. Dirichlet boundary conditions are of course enforced by setting nodal variables x a and T a to appropriate values on @ u B 0 and @ T B 0 . Subscript h in I h denotes the fact that volume integrals will generally be computed by numerical quadrature, based on the chosen discretization T h . Variations of incremental potential I h are now taken with respect to nodal unknowns, and stationarity conditions are written where P h , d h eff and H h are, respectively, given by The stationarity conditions (203a, b) yield discrete mechanical and thermal balance equations, which can alternatively be written where internal and external nodal forces and fluxes are given by In the above expressions, S denotes the element assembly operator, while O h e are the elementary domains. Note that balance equations (207) suggest a (nonlinear) quasi-stationary problem, even in the presence of heat capacity terms in the constitutive model. These terms, traditionally treated as transient terms, are here included in the effective entropy variation d eff , and thermal balance (both in its continuum or discrete form) is actually treated through an entropy balance equation written over the time increment. It is thus a different, yet equivalent, treatment from the classical approach, which consists in writing a balance equation for instantaneous heat fluxes. Our approach instead yields a nonlinear, rate-dependent, quasi-stationary set of balance equations.

Mixed formulations
At the core of the variational formulation (69) lies the energy functional i.e. the incremental functional without terms coming from external applied thermo-mechanical loads, and where we dropped the arguments denoting parametric dependence of potentials, in order to lighten notations for the rest of the section. In the context of a finite element approach, this integral is approximated by with elementary contributions given by In the case of incompressible material behavior, or even quasiincompressible behavior such as in cases dominated by plastic strains, it is well known [e.g. see Hughes (2000) or Simo and Hughes (1998)] that a standard finite element formulation will lead to locking (overestimated stiffness) for linear triangles and tetrahedrons, as well as for bilinear quadrangles and trilinear hexahedrons. In order to overcome this difficulty, we adopt a hybrid (or mixed) formulation, and, following the approach proposed in Simo and Taylor (1991), add a penalty term to integral (211): where y h e and p h e are piecewise-constant volumic deformations (not to be confused with temperature increments in this setting) and pressures, respectively (constant over each element O 0 e ). The modified deformation gradient The variational principle then becomes Euler-Lagrange equations pertaining to y h and p h are given by, respectively: Equation (217) immediately yields the expression for the piecewise constant volumic deformation Considering that where we have introduced the notation: Eq. (216) then yields We can now proceed to compute elementary contributions to nodal forces and fluxes. Since we have we can write We then have which can alternatively be written The latter expression defines a recovery method for Kirchhoff stress field: The nodal internal forces array is then given by Nodal internal fluxes are identical to those obtained through the standard formulation (208c), where the net entropy variation is computed by

Alternative variational Ritz-Galerkin approximations
Galerkin approximations are not limited to the piecewise polynomial functions of finite elements. Various Ritz-Galerkin approximations can be derived by using different test functions in the time-continuous and/or time-discrete variational principles.
A particular example that we detail here is that of adiabatic shear bands. Adiabatic shear bands are localized regions of intense shear (which take the form of a band in 2D, of a layer in 3D), where thermo-mechanical coupling effects play a prominent role. The thickness of these bands or layers can be very small compared to other characteristic lengths of the problem, in which case explicit resolving of the displacement (strain) and temperature fields with the band becomes prohibitive. As illustrated in Yang, Mota, and Ortiz (2005), shear bands can then be represented as discontinuities (in displacement) to which one associates a specific behavior. The macroscopic behavior of the shear band can be obtained by assuming specific (families of ) profiles of (plastic) strain and temperature: for example, Yang et al. (2005) considered uniform strain rate, directly computed from the macroscopic velocity jump and the thickness of the band, together with a Gaussian temperature profile. The main inconvenient of this specific choice is that the band thickness must be chosen a priori and becomes a material parameter, although it is actually an evolving quantity depending on global strain, strain rate, and temperature. This inconvenience was partly overcome in Su (2012), where the analytical solutions of Leroy and Molinari (1992) are used as basis functions. These functions allow to represent arbitrary (and independent) thicknesses of localization for strain and temperature. This latter approach was applied both to the stationary case (adiabatic shear bands in the established regime) and to the transient or evolutionary case (in an incremental time-discrete context).

EXAMPLES OF APPLICATIONS
This section provides a (very) brief illustration of some applications for which a variational formulation offers significant advantages: multiscale and adaptive approaches.

Variational multiscale models
Variational formulations are well adapted to multiscale approaches. Indeed they tend to provide systematic rules to operate scale transitions: microstructures and/or micro-scale fields should minimize (optimize in general) some overall macroscopic energy functional. For example, postulating a specific family of recursive microstructures (laminates), it is possible to construct macroscopic (variational) constitutive models accounting for the formation of evolving subgrain microstructures in polycrystalline materials submitted to (very) large plastic deformations (Aubry & Ortiz, 2003;Conti & Ortiz, 2005;Ortiz, Repetto, & Stainier, 2000).
We can also look at the more general problem of constructing constitutive models able to describe the overall (homogenized) behavior of heterogeneous materials. Under isothermal conditions, the homogenization problem takes the following variational form: where B 0 is the domain occupied by a RVE, and K( F) is the set of admissible displacement mappings, depending on the kind of boundary conditions imposed on the RVE in order to enforce the average gradient of deformation F. Thanks to the use of the effective incremental potential W, this expression remains valid for all heterogeneous materials which phases' behavior can be described in the framework of GSM. This homogenization problem can then be solved numerically (Bleier and Mosler, 2013;Miehe, 2002) or semianalytically in some specific cases (Brassart, Stainier, Doghri, & Delannay, 2011, 2012. In this latter example, a variational mean-field approach was derived for composites materials with elasto-viscoplastic phases (small strains).
Note that all these contributions have, for the moment, been restricted to the isothermal case.

Variational adaptive methods
Variational approaches are also well suited to adaptive methods. Once again, the main key is the use of the incremental potential W, which allows to extend methods initially designed for (hyper-)elasticity to general dissipative behaviors.
Optimizing with respect to mesh position, in addition to the optimization with respect to the displacement field, leads to the variational Arbitrary Lagrangian-Eulerian formulation of Mosler and Ortiz (2006). In a variation to that method, the same authors also proposed a technique to locally refine mesh discretization on variational (minimization of effective incremental energy) approach (Mosler & Ortiz, 2007).

CONCLUSIONS
In this chapter, we presented an overview of a variational approach to coupled thermo-mechanical boundary-value problems involving nonlinear dissipative behaviors. Starting from a thermodynamic description in the framework of GSM, time-continuous and time-discrete variational principles can be established by writing energy or energy-rate functionals to be optimized with respect to material state fields. Having adopted a local approach, optimization with respect to internal variables can be done at the local (i.e. material point) level, the remaining optimization problem then corresponding to an effective thermo-elastic boundary-value problem. The case of thermo-visco-elastic and thermo-elasto-visco-plastic behaviors was treated in more details. First noteworthy results are that, for isotropic materials, a spectral description of thermodynamic potentials leads to efficient update algorithms. Another one is that the predictor-corrector algorithmic structure of the classical radial return of computational plasticity can be recovered in most cases. The case of damage, either elastic damage or combined with one of the previous models, was not treated here, although variational approaches can also be developed in that direction (e.g. Balzani & Ortiz, 2012;Kintzel & Mosler, 2010, 2011. We then looked at variational approximation methods for the boundary-value problem, either by finite elements or by other Ritz-Galerkin approaches, and briefly illustrated applications in the context of multiscale problems and adaptive techniques. Another potential class of applications of the variational approach to coupled thermo-mechanics is in the field of algorithmic solution techniques: partitioned solvers, staggered schemes, could be derived by taking advantage of the variational structure (e.g. symmetry, but also the fact that the solution corresponds to an optimum).
An interesting perspective is the extension of this variational framework to other types of coupling. Indeed there is a clear parallel between thermomechanics and diffusion-mechanics problems (or thermo-diffusionmechanics since temperature typically plays an important role in this type of problems). Beyond that, coupling between mechanics and electromagnetism also leans itself to a variational treatment (see for example François-Lavet, Henrotte, Stainier, Noels, & Geuzaine, 2013;Miehe, Rosato, & Kiefer, 2011;Thomas & Triantafyllidis, 2009).
As it was previously mentioned, another interest of variational formulations is that they are well adapted to mathematical analysis. From this point of view, questions such as existence and unicity of solutions to problems combining dissipative, time-dependent mechanical behavior and heat transfer can benefit from the proposed variational treatment. A particularly interesting case in that regard seems to be the problem of thermo-visco-elasticity, to the solution of which no extremum property can be clearly associated.
Finally a challenging open question is that of a variational characterization of entire trajectories in the space of state fields. Indeed, the current approach relies on the construction of a separate variational problem at each time increment. The possibility to define optimal trajectories for rate-dependent dissipative systems, as can be done through Hamilton's principle in the dynamics of conservative systems is an appealing objective, not yet reached despite some first results (Conti & Ortiz, 2008;Mielke & Ortiz, 2008).