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

Laurent Stainier

► To cite this version:

Laurent Stainier. A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors. Advances in Applied Mechanics, Academic press, pp.69-126, 2013, 10.1016/B978-0-12-396522-6.00002-5. hal-01007469

HAL Id: hal-01007469

https://hal.science/hal-01007469

Submitted on 19 Jun 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors

## Laurent Stainier

Research Institute in Civil and Mechanical Engineering (GeM), Ecole Centrale Nantes, 1 rue de la Noë, F-44321 Nantes, France

## Contents

1. Introduction

2. General Modeling Framework

2.1 Local thermodynamic model

2.2 Balance equations

3. Variational Formulation of Coupled Thermo-Mechanical Boundary-Value Problems

3.1 Variational updates

3.2 Variational boundary-value problem

3.3 Dynamics

3.4 Linearization

4. Thermo-Visco-Elasticity

4.1 Linearized kinematics

4.2 Finite thermo-visco-elasticity

5. Thermo-Elasto-Visco-Plasticity

5.1 Crystal plasticity

5.2 Macroscopic plasticity

6. Numerical Approximation Methods

6.1 Variational finite element approximations

6.2 Alternative variational Ritz-Galerkin approximations

7. Examples of Applications

7.1 Variational multiscale models

7.2 Variational adaptive methods

8. Conclusions

Acknowledgements

References
Abstract

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.

# 1. INTRODUCTION

Variational principles have played an important role in mechanics for several decades, if not more than a century (see for example Lanczos, 1986 or 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 (1956, 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, 1993; Simo & 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 thermo-mechanical boundary-value problems is initially presented within a quasi-stationary 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 thermo-elasto-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.
## 2. GENERAL MODELING FRAMEWORK

### 2.1. 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 \in \mathbb{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 $\mathbf{X}$ to spatial points $\mathbf{x}$:

$$ \varphi: B_0 \times [0, t_f] \to \mathbb{R}^3, \mathbf{x} = \varphi(\mathbf{X}, t) \quad (1) $$

and the absolute temperature field:

$$ T: B_0 \times [0, t_f] \to \mathbb{R}^+, T = T(\mathbf{X}, t) > 0, \quad (2) $$

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 $\mathcal{B}$ when there is no ambiguity):

$$ B_t = \{\mathbf{x} = \varphi(\mathbf{X}, t), \forall \mathbf{X} \in B_0\}. \quad (3) $$

We adopt a local state approach, where the local material state is described by the deformation gradient $\mathbf{F}(\mathbf{X}, t)$:

$$ \mathbf{F}: B_0 \times [0, t_f] \to GL^+(3, \mathbb{R}), \mathbf{F} = \nabla_0 \varphi(\mathbf{X}, t) \quad (4) $$

where $GL^+(\mathbb{R}, 3)$ is the space of second-order tensors on $\mathbb{R}^3$ with positive determinant, and temperature $T(\mathbf{X}, t)$. Since we want to consider dissipative behaviors, a set of internal variables $\mathbf{Z}$ is added to external state variables $\mathbf{F}$ and $T$. The exact nature of $\mathbf{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(\mathbf{F}, T, \mathbf{Z})$ and a dissipation pseudo-potential $\Delta(\dot{\mathbf{F}}, \dot{\mathbf{Z}}; \mathbf{F}, T, \mathbf{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
from the free energy, and a dissipative part, derived from the dissipation pseudo-potential (Ziegler, 1977):

$$ \mathbf{P} = \mathbf{P}^e + \mathbf{P}^d \qquad (5) $$

$$ \mathbf{P}^e = \frac{\partial W}{\partial \mathbf{F}}(\mathbf{F}, T, \mathbf{Z}) \qquad (6) $$

$$ \mathbf{P}^d = \frac{\partial \Delta}{\partial \dot{\mathbf{F}}} (\dot{\mathbf{F}}, \dot{\mathbf{Z}}; \mathbf{F}, T, \mathbf{Z}) \qquad (7) $$

and

$$ \mathbf{Y} = \mathbf{Y}^e + \mathbf{Y}^d \qquad (8) $$

$$ \mathbf{Y}^e = \frac{\partial W}{\partial \mathbf{Z}}(\mathbf{F}, T, \mathbf{Z}) \qquad (9) $$

$$ \mathbf{Y}^d = \frac{\partial \Delta}{\partial \dot{\mathbf{Z}}} (\dot{\mathbf{F}}, \dot{\mathbf{Z}}; \mathbf{F}, T, \mathbf{Z}) \qquad (10) $$

where $\mathbf{P}$ is the first Piola–Kirchhoff (or Piola) stress tensor, conjugate to $\mathbf{F}$, and $\mathbf{Y}$ the forces conjugate to internal variables $\mathbf{Z}$. What differentiates external and internal variables is that the latter should produce no net work, i.e.

$$ \mathbf{Y} \cdot \dot{\mathbf{Z}} = (\mathbf{Y}^e + \mathbf{Y}^d) \cdot \dot{\mathbf{Z}} = 0, \quad \forall \dot{\mathbf{Z}}. \qquad (11) $$

As a consequence, one must have that

$$ \mathbf{Y}^d = -\mathbf{Y}^e \Leftrightarrow \frac{\partial W}{\partial \mathbf{Z}}(\mathbf{F}, T, \mathbf{Z}) + \frac{\partial \Delta}{\partial \dot{\mathbf{Z}}} (\dot{\mathbf{F}}, \dot{\mathbf{Z}}; \mathbf{F}, T, \mathbf{Z}) = 0. \qquad (12) $$

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 $\eta$ is given by

$$ \rho_0 \eta = - \frac{\partial W}{\partial T}(\mathbf{F}, T, \mathbf{Z}) \qquad (13) $$

where $\rho_0$ is the density in the reference configuration. In the following, it will prove useful to introduce the internal energy density potential $U(\mathbf{F}, \eta, \mathbf{Z})$:

$$ U(\mathbf{F}, \eta, \mathbf{Z}) = \sup_T [\rho_0 \eta T + W(\mathbf{F}, T, \mathbf{Z})] \qquad (14) $$

such that

$$ T = \frac{1}{\rho_0} \frac{\partial U}{\partial \eta}(\mathbf{F}, \eta, \mathbf{Z}). \qquad (15) $$
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:

$$W(QF, T, Z) = W(F, T, Z) \quad \forall Q \in SO(3, \mathbb{R}) \quad (16)$$

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

$$P = FS \quad \text{with} \quad S = 2 \frac{\partial W}{\partial C} \quad (17)$$

where $S$ is the second Piola–Kirchhoff (symmetric) stress tensor.

## 2.2. Balance equations

We now briefly recall local balance equations corresponding to conservation principles (mass, momentum, and energy). In Lagrangian description, mass conservation yields

$$\rho \det[F] = \rho_0 \quad (18)$$

where $\rho$ is the mass density in the deformed configuration. Linear momentum conservation yields

$$\rho_0 \ddot{\boldsymbol{\varphi}} = \nabla_0 \cdot \boldsymbol{P}^T + \rho_0 \boldsymbol{b} \quad (19)$$

where $\boldsymbol{b}$ denotes applied bulk forces (per unit mass), while angular momentum conservation yields

$$\boldsymbol{P}\boldsymbol{F}^T = \boldsymbol{F}\boldsymbol{P}^T. \quad (20)$$

Note that, given Eq. (17), this will be automatically verified for objective (frame invariant) constitutive models. Conservation of energy yields

$$\rho_0 T \dot{\eta} = \boldsymbol{P}^d \cdot \dot{\boldsymbol{F}} + \boldsymbol{Y}^d \cdot \dot{\boldsymbol{Z}} - \nabla_0 \cdot \boldsymbol{H} + \rho_0 Q \quad (21)$$

where $\boldsymbol{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 $\mathcal{D}_{\text{int}} = \boldsymbol{P}^d \cdot \dot{\boldsymbol{F}} + \boldsymbol{Y}^d \cdot \dot{\boldsymbol{Z}}$. Then, using the definition (13) of entropy, we can obtain the heat equation:

$$C \dot{T} = T \frac{\partial^2 W}{\partial \boldsymbol{F} \partial T} \cdot \dot{\boldsymbol{F}} + T \frac{\partial^2 W}{\partial \boldsymbol{Z} \partial T} \cdot \dot{\boldsymbol{Z}} + \mathcal{D}_{\text{int}} - \nabla_0 \cdot \boldsymbol{H} + \rho_0 Q \quad (22)$$
where the heat capacity C, defined by

$$C = -T \frac{\partial^2 W}{\partial T^2}(\mathbf{F}, T, \mathbf{Z}) \qquad (23)$$

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_{\text{int}}$. The Clausius–Duhem inequality (second principle of thermodynamics) then writes

$$T \dot{\Gamma} = D_{\text{int}} - \frac{1}{T} \mathbf{H} \cdot \nabla_0 T \geq 0 \qquad (24)$$

where $\dot{\Gamma}$ denotes the net entropy production rate. Choosing a dissipation pseudo-potential $\Delta(\dot{\mathbf{F}}, \dot{\mathbf{Z}}; \mathbf{F}, T, \mathbf{Z})$ which is convex with respect to $\dot{\mathbf{F}}$ and $\dot{\mathbf{Z}}$, nonnegative and such that $\Delta(\mathbf{0}, \mathbf{0}; \mathbf{F}, T, \mathbf{Z}) = 0$ will ensure that $D_{\text{int}} \geq 0$. Following Biot (1958), we then introduce a conduction potential $\chi(\mathbf{G}; \mathbf{F}, T, \mathbf{Z})$, such that

$$\mathbf{H} = \frac{\partial \chi}{\partial \mathbf{G}}(\mathbf{G}; \mathbf{F}, T, \mathbf{Z}) \quad \text{where} \quad \mathbf{G} = - \frac{\nabla_0 T}{T}. \qquad (25)$$

Assuming convexity and nonnegativeness of $\chi$, together with the condition $\chi(\mathbf{0}; \mathbf{F}, T, \mathbf{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 $\mathbf{L} = \dot{\mathbf{F}} \mathbf{F}^{-1}$, and its symmetric part defines the Eulerian strain rate

$$\mathbf{D} = \frac{1}{2} (\mathbf{L} + \mathbf{L}^T) = \frac{1}{2} \mathbf{F}^{-T} \dot{\mathbf{C}} \mathbf{F}^{-1}. \qquad (26)$$

The Cauchy stress tensor $\boldsymbol{\sigma}$ is given by

$$\boldsymbol{\sigma} = \frac{1}{J} \mathbf{P} \mathbf{F}^T = \frac{1}{J} \mathbf{F} \mathbf{S} \mathbf{F}^T \qquad (27)$$

where $J = \det[\mathbf{F}]$. As a consequence, we have that

$$\mathbf{P} \cdot \dot{\mathbf{F}} = J \boldsymbol{\sigma} \cdot \mathbf{D}. \qquad (28)$$
The linear momentum conservation equation then becomes

$$ \rho \ddot{\boldsymbol{\varphi}} = \nabla \cdot \boldsymbol{\sigma} + \rho \boldsymbol{b} \qquad (29) $$

where we have the used mass conservation equation, while conservation of energy becomes

$$ \rho T \dot{\eta} = \boldsymbol{\sigma}^d \cdot \boldsymbol{D} + \boldsymbol{Y}^d \cdot \dot{\boldsymbol{Z}} - \nabla \cdot \boldsymbol{h} + \rho Q \qquad (30) $$

with $\boldsymbol{h} = J^{-1} F \boldsymbol{H}$ is the (Eulerian) heat flux vector.

### 3. VARIATIONAL FORMULATION OF COUPLED THERMO-MECHANICAL BOUNDARY-VALUE PROBLEMS

#### 3.1. 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.

##### 3.1.1 Local evolution problem

At a given material point, evolution of internal variables is ruled by Eq. (12). Taking into account the convexity of $\Delta$ w.r.t. $\dot{\boldsymbol{Z}}$, this equation can be interpreted as the stationarity condition associated to the following variational principle:

$$ \dot{\boldsymbol{Z}} = \arg \inf_{\dot{\boldsymbol{Z}}} D(\dot{\boldsymbol{F}}, \dot{\eta}, \dot{\boldsymbol{Z}}, T; \boldsymbol{F}, \eta, \boldsymbol{Z}) \qquad (31a) $$

with
$$D(\dot{\mathbf{F}}, \dot{\eta}, \dot{\mathbf{Z}}, T; \mathbf{F}, \eta, \mathbf{Z}) = \frac{d}{dt}[U(\mathbf{F}, \eta, \mathbf{Z})] - \rho_0 \dot{\eta} T + \Delta \left( \frac{T}{\Theta} \dot{\mathbf{F}}, \frac{T}{\Theta} \dot{\mathbf{Z}}; \mathbf{F}, \Theta(\mathbf{F}, \eta, \mathbf{Z}), \mathbf{Z} \right). \quad (31b)$$

In the above expression, $\Theta(\mathbf{F}, \eta, \mathbf{Z})$ is the equilibrium or internal temperature, as defined by Eq. (15):

$$\Theta(\mathbf{F}, \eta, \mathbf{Z}) \equiv \frac{1}{\rho_0} \frac{\partial U}{\partial \eta}(\mathbf{F}, \eta, \mathbf{Z}) \quad (32)$$

while $T$ will be called external temperature. Local thermal equilibrium imposes that $\Theta = T$, which corresponds to stationarity of $D$ w.r.t. $\dot{\eta}$:

$$\underset{\dot{\eta}}{\operatorname{stat}} D(\dot{\mathbf{F}}, \dot{\eta}, \dot{\mathbf{Z}}, T; \mathbf{F}, \eta, \mathbf{Z}) \iff T = \Theta(\mathbf{F}, \eta, \mathbf{Z}). \quad (33)$$

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 $\Delta(\dot{\mathbf{F}}, \dot{\eta}, \dot{\mathbf{Z}})$, but we will not pursue this possibility here.

The apparent complexity of the function $D(\dot{\mathbf{F}}, \dot{\eta}, \dot{\mathbf{Z}}, T; \mathbf{F}, \eta, \mathbf{Z})$ is in practice motivated by the following properties. Derivatives of the function w.r.t. $\dot{\mathbf{F}}$ and $T$ are, respectively, given by:

$$\frac{\partial D}{\partial \dot{\mathbf{F}}} = \frac{\partial U}{\partial \mathbf{F}}(\mathbf{F}, \eta, \mathbf{Z}) + \frac{T}{\Theta} \frac{\partial \Delta}{\partial \dot{\mathbf{F}}} \left( \frac{T}{\Theta} \dot{\mathbf{F}}, \frac{T}{\Theta} \dot{\mathbf{Z}}; \mathbf{F}, \Theta(\mathbf{F}, \eta, \mathbf{Z}), \dot{\mathbf{Z}} \right) \quad (34)$$

and

$$\frac{\partial D}{\partial T} = -\rho_0 \dot{\eta} + \frac{1}{\Theta} \left[ \frac{\partial \Delta}{\partial \dot{\mathbf{F}}} \left( \frac{T}{\Theta} \dot{\mathbf{F}}, \frac{T}{\Theta} \dot{\mathbf{Z}} \right) \cdot \dot{\mathbf{F}} + \frac{\partial \Delta}{\partial \dot{\mathbf{Z}}} \left( \frac{T}{\Theta} \dot{\mathbf{F}}, \frac{T}{\Theta} \dot{\mathbf{Z}} \right) \cdot \dot{\mathbf{Z}} \right]. \quad (35)$$

Then, if we define the effective rate potential $D_{\text{eff}}$ as

$$D_{\text{eff}}(\dot{\mathbf{F}}, \dot{\eta}, T; \mathbf{F}, \eta, \mathbf{Z}) = \inf_{\dot{\mathbf{Z}}} D(\dot{\mathbf{F}}, \dot{\eta}, \dot{\mathbf{Z}}, T; \mathbf{F}, \eta, \mathbf{Z}), \quad (36)$$

we can write that

$$\frac{\partial D_{\text{eff}}}{\partial \dot{\mathbf{F}}} = \mathbf{P}^e(\mathbf{F}, \Theta(\mathbf{F}, \eta, \mathbf{Z}), \mathbf{Z}) + \frac{T}{\Theta} \mathbf{P}^d \left( \frac{T}{\Theta} \dot{\mathbf{F}}, \frac{T}{\Theta} \dot{\mathbf{Z}}; \mathbf{F}, \Theta(\mathbf{F}, \eta, \mathbf{Z}), \mathbf{Z} \right) \quad (37)$$

$$\frac{\partial D_{\text{eff}}}{\partial T} = -\rho_0 \dot{\eta} + \frac{1}{T} \underbrace{\left[ \frac{\partial \Delta}{\partial \dot{\mathbf{F}}} \left( \frac{T}{\Theta} \dot{\mathbf{F}}, \frac{T}{\Theta} \dot{\mathbf{Z}} \right) \cdot \frac{T}{\Theta} \dot{\mathbf{F}} + \frac{\partial \Delta}{\partial \dot{\mathbf{Z}}} \left( \frac{T}{\Theta} \dot{\mathbf{F}}, \frac{T}{\Theta} \dot{\mathbf{Z}} \right) \cdot \frac{T}{\Theta} \dot{\mathbf{Z}} \right]}_{=D_{\text{int}}} \quad (38)$$
Thus, provided that $T = \Theta(F, \eta, Z)$, it is then seen that function $D_{\text{eff}}(\dot{F}, \dot{\eta}, 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_{\text{int}}$). For example, looking at heat equation (21), we see that local adiabatic behavior would then be characterized by

$$ \frac{\partial D_{\text{eff}}}{\partial T}(\dot{F}, \dot{\eta}, T; F, \eta, Z) = 0 \\ \Leftrightarrow D_{\text{ad}}(\dot{F}; F, \eta, Z) = \underset{\dot{\eta}}{\text{stat}} \underset{T}{\text{stat}} D_{\text{eff}}(\dot{F}, \dot{\eta}, T; F, \eta, Z) \quad (39) $$

and

$$ \mathbf{P} = \frac{\partial D_{\text{ad}}}{\partial \dot{F}} (\dot{F}; F, \eta, Z) \quad (40) $$

i.e. $D_{\text{ad}}$ acts as a rate potential for the stress tensor $P$ under local adiabatic conditions.

### 3.1.2 Local time-discrete constitutive problem

Under most circumstances, nonlinear problems involving history-dependent 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, \eta_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, \eta\}$) 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 $\dot{U} - \rho_0 \dot{\eta} T = \dot{W} + \rho_0 \eta \dot{T}$, we write

$$ I(F, T, Z; F_0, T_0, Z_0) \approx \int_{t_0}^{t} D(\dot{F}, \dot{\eta}, \dot{Z}, T(\tau); F(\tau), \eta(\tau), Z(\tau)) d\tau \\ = W(F, T, Z) - W_0 + \rho_0 \eta_0 \delta T + \delta t \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t}; F(\tau), T(\tau), Z(\tau) \right) \right\rangle \quad (41) $$

where $W_0 = W(F_0, T_0, Z_0)$, $\delta t = t - t_0$, $\delta F = F - F_0$, $\delta T = T - T_0$, $\delta Z = Z - Z_0$. 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/\Theta$ have been replaced by $T/T_0$.
The incremental variational update takes the form of the following minimization problem:

$$ W(F, T; F_0, T_0, Z_0) = \inf_{Z} I(F, T, Z; F_0, T_0, Z_0). \quad (42) $$

The stationarity equation corresponding to this variational problem is given by

$$ \frac{\partial W}{\partial Z}(F, T, Z) + \delta t \frac{\partial}{\partial Z} \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t}; F(\tau), T(\tau), Z(\tau) \right) \right\rangle = 0. \quad (43) $$

In addition, we will have that

$$ \frac{\partial W}{\partial F} = P^e(F, T, Z) + \delta t \frac{\partial}{\partial F} \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t}; F(\tau), T(\tau), Z(\tau) \right) \right\rangle \quad (44) $$

$$ \frac{\partial W}{\partial T} = -\rho_0(\eta - \eta_0) + \delta t \frac{\partial}{\partial T} \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t}; F(\tau), T(\tau), Z(\tau) \right) \right\rangle. \quad (45) $$

Consistency of the incremental update thus requires that

$$ \lim_{\delta t \to 0} \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t}; F_\alpha, T_\alpha, Z_\alpha \right) \right\rangle = \Delta(\dot{F}, \dot{Z}; F, T, Z) \quad (46) $$

but also that

$$ \lim_{\delta t \to 0} \delta t \frac{\partial}{\partial F} \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t}; F(\tau), T(\tau), Z(\tau) \right) \right\rangle = \frac{\partial \Delta}{\partial \dot{F}} (\dot{F}, \dot{Z}; F, T, Z) \quad (47) $$

$$ \lim_{\delta t \to 0} \delta t \frac{\partial}{\partial Z} \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t}; F(\tau), T(\tau), Z(\tau) \right) \right\rangle = \frac{\partial \Delta}{\partial \dot{Z}} (\dot{F}, \dot{Z}; F, T, Z) \quad (48) $$

and

$$ \lim_{\delta t \to 0} \frac{\partial}{\partial T} \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t}; F(\tau), T(\tau), Z(\tau) \right) \right\rangle = \frac{1}{T} \left( \frac{\partial \Delta}{\partial \dot{F}} \cdot \dot{F} + \frac{\partial \Delta}{\partial \dot{Z}} \cdot \dot{Z} \right). \quad (49) $$

As discussed in Stainier (2011), the most intuitive expression
$$ \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t} ; F(\tau), T(\tau), Z(\tau) \right) \right\rangle = \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t} ; F_{\alpha}, T_{\alpha}, Z_{\alpha} \right) \quad (50) $$

where {$F, T, Z$}$_{α}$ = (1 - α) {$F_0, T_0, Z_0$} + α {$F, T, Z$} (with algorithmic parameter α ∈ [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):

$$ \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t} ; F_{\alpha}, T_{\alpha}, Z_{\alpha} \right) \right\rangle = \frac{T_0}{T} \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t} ; F_{\alpha}, T_{\alpha}, Z_{\alpha} \right) + \frac{\delta T}{T} \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t} ; F_{\alpha}, T_{\alpha}, Z_{\alpha} \right). \quad (51) $$

One can see that in the case α = 0, or in the case of no parametric dependence of Δ on the current state, the above expression reduces to

$$ \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t} ; F_{\alpha}, T_{\alpha}, Z_{\alpha} \right) \right\rangle_{\alpha=0} = \Delta \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t} \right). \quad (52) $$

In that specific case, the internal variable update becomes

$$ \frac{\partial W}{\partial Z}(F, T, Z) + \frac{T}{T_0} \frac{\partial \Delta}{\partial \dot{Z}} \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t} \right) = 0 \quad (53) $$

and we also have

$$ \frac{\partial W}{\partial F} = P^e(F, T, Z) + \frac{T}{T_0} \frac{\partial \Delta}{\partial \dot{F}} \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t} \right) \quad (54) $$

$$ \frac{\partial W}{\partial T} = -\rho_0(\eta - \eta_0) + \frac{1}{T_0} \left[ \frac{\partial \Delta}{\partial \dot{F}} \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t} \right) \cdot \delta F + \frac{\partial \Delta}{\partial \dot{Z}} \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t} \right) \cdot \delta Z \right]. \quad (55) $$

We can interpret these results as providing approximate incremental expressions for the viscous stress, the thermodynamic forces conjugate to internal variables:

$$ P^d \approx \frac{T}{T_0} \frac{\partial \Delta}{\partial \dot{F}} \left( \frac{T}{T_0} \frac{\delta F}{\delta t}, \frac{T}{T_0} \frac{\delta Z}{\delta t} \right) \quad (56) $$
$$ \mathbf{Y}^d \approx \frac{T}{T_0} \frac{\partial \Delta}{\partial \dot{\mathbf{Z}}} \left( \frac{T}{T_0} \frac{\delta \mathbf{F}}{\delta t}, \frac{T}{T_0} \frac{\delta \mathbf{Z}}{\delta t} \right) \qquad (57) $$

and the dissipation (average value over the time step)

$$ \mathcal{D}_{\text{int}} \approx \frac{1}{T} \left( \mathbf{P}^d \cdot \frac{\delta \mathbf{F}}{\delta t} + \mathbf{Y}^d \cdot \frac{\delta \mathbf{Z}}{\delta t} \right). \qquad (58) $$

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

$$ \frac{\partial W}{\partial \mathbf{F}} = \mathbf{P}^e + \mathbf{P}^d = \mathbf{P} \qquad (59) $$

$$ \frac{\partial W}{\partial T} = -\rho_0(\eta - \eta_0) + \delta t \mathcal{D}_{\text{int}} = \rho_0 \delta \eta^{\text{eff}}. \qquad (60) $$

In the more general case ($\alpha \neq 0$) additional terms will appear in the incremental approximations of $\mathbf{P}^d$, $\mathbf{Y}^d$, and $\mathcal{D}_{\text{int}}$, but the resulting expressions will remain consistent, i.e. will tend toward continuous values as $\delta t \to 0$.

A detailed study of the influence of algorithmic parameter $\alpha$ on precision and convergence properties of the variational update in the case of isothermal elasto-visco-plasticity can be found in Brassart and Stainier (2012). Among other results, it is shown that it is always possible to find a value of parameter $\alpha$ 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 $\alpha$ vanishes as the time increment size decreases.

## 3.2. 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.
### 3.2.1 Rate problem

Let us now consider the quasi-static boundary value problem consisting in determining fields $\dot{\varphi}$, $\dot{\eta}$, $\dot{Z}$, and $T$ on $B_0$ at a given time $t$ (the current state $\{F, \eta, Z\}$ is assumed to be known), verifying the differential equations

$$ \left\{ \begin{array}{l} \nabla_0 \cdot \mathbf{P}^T + \rho_0 \mathbf{b}(t) = 0 \\ \rho_0 T \dot{\eta} = \underbrace{\mathbf{P}^d(\dot{\mathbf{F}}, \dot{\mathbf{Z}}; T) \cdot \dot{\mathbf{F}} + \mathbf{Y}^d(\dot{\mathbf{F}}, \dot{\mathbf{Z}}; T) \cdot \dot{\mathbf{Z}}}_{=\mathcal{D}_{\text{int}}} - \nabla_0 \cdot \mathbf{H}(\mathbf{G}; T) + \rho_0 \mathbf{Q}(t) \end{array} \right. \qquad (61) $$

together with local constitutive equation (12) and boundary conditions: $\varphi = \bar{\varphi}(t)$ on $\partial_u B_0$, $\mathbf{P} \cdot \mathbf{N} = \bar{\mathbf{t}}(t)$ on $\partial_\sigma B_0$ ($\partial_u B_0 \cup \partial_\sigma B_0 = \emptyset$), $T = \bar{T}(t)$ on $\partial_T B_0$, and $\mathbf{H} \cdot \mathbf{N} = \bar{\mathbf{H}}(t)$ on $\partial_H B_0$ ($\partial_T B_0 \cup \partial_H B_0 = \partial B_0$ and $\partial_T B_0 \cap \partial_H B_0 = \emptyset$), where $\bar{\varphi}(t)$, $\bar{\mathbf{t}}(t)$, $\bar{T}(t)$, and $\bar{\mathbf{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

$$ \Phi(\dot{\varphi}, \dot{\eta}, T) = \int_{B_0} \left[ D_{\text{eff}}(\nabla_0 \dot{\varphi}, \dot{\eta}, T) - \chi \left( -\frac{\nabla_0 T}{T}; F, \Theta(F, \eta, Z), Z \right) \right] dV \\ - \int_{B_0} \rho_0 \mathbf{b}(t) \cdot \dot{\varphi} dV - \int_{\partial_\sigma B_0} \bar{\mathbf{t}}(t) \cdot \dot{\varphi} dS \\ + \int_{B_0} \rho_0 \mathbf{Q}(t) \frac{T}{\Theta} dV - \int_{\partial_H B_0} \bar{\mathbf{H}}(t) \frac{T}{\Theta} dS \qquad (62) $$

where we have omitted to explicitly indicate the sufficiently clear dependence of $D$ and $\Phi$ on the current state for brevity. Then balance equations (61) correspond to stationarity conditions for $\Phi$:

$$ \langle D_{\dot{\varphi}} \Phi, \tilde{\varphi} \rangle = \int_{B_0} \frac{\partial \mathcal{D}_{\text{eff}}}{\partial \dot{\mathbf{F}}} \cdot \nabla_0 \tilde{\varphi} dV - \int_{B_0} \rho_0 \mathbf{b}(t) \cdot \tilde{\varphi} dV - \int_{\partial_\sigma B_0} \bar{\mathbf{t}}(t) \cdot \tilde{\varphi} dS = 0 \quad \forall \tilde{\varphi} \text{ adm.} \qquad (63) $$

$$ \langle D_T \Phi, \tilde{T} \rangle = \int_{B_0} \left[ \frac{\partial \mathcal{D}_{\text{eff}}}{\partial T} \tilde{T} + \frac{\partial \chi}{\partial \mathbf{G}} \cdot \nabla_0 \left( \frac{\tilde{T}}{T} \right) \right] dV + \int_{B_0} \rho_0 \mathbf{Q}(t) \frac{\tilde{T}}{\Theta} dV \\ - \int_{\partial_H B_0} \bar{\mathbf{H}}(t) \frac{\tilde{T}}{\Theta} dS = 0 \quad \forall \tilde{T} \text{ adm.} \qquad (64) $$

while stationarity condition w.r.t. $\dot{\eta}$ yields $T = \Theta$, as seen previously. Applying Green–Ostrogradsky’s theorem to the above stationarity equations, and
replacing $\Theta$ 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 $\Phi$, 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:

$$ \{\dot{\varphi}, \dot{\eta}, T\} = \arg \inf_{\dot{\varphi}, \dot{\eta}} \sup_{T} \Phi(\dot{\varphi}, \dot{\eta}, T). \quad (65) $$

In the case of thermo-elasticity (i.e. $Z = \emptyset$, $\Delta = 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

$$ \{\dot{\varphi}, \dot{\eta}, T\} = \arg \inf_{\dot{\varphi}} \sup_{\dot{\eta}} \Phi_{\text{thel}}(\dot{\varphi}, \dot{\eta}, T). \quad (66) $$

If we exclude viscous stresses [i.e. $\Delta = \Delta(\dot{Z}; F, T, Z)$], then the effective rate potential $D_{\text{eff}}$ is concave in $T$ [this results from the implicitly assumed convexity of $\Delta(\dot{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.

### 3.2.2 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 $\{\varphi_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 = \nabla_0 \varphi, T\}$ are given. External fields $\{\varphi, T\}$ can themselves be computed as optimizers of the following functional:

$$ \begin{aligned} \mathcal{I}(\varphi, T) = & \int_{B_0} \left[ W(\nabla_0 \varphi, T; F_0, T_0, Z_0) - \delta t \left\langle \chi \left( -\frac{\nabla_0 T}{T}; F_\alpha, T_\alpha, Z_\alpha \right) \right\rangle \right] dV \\ & - \int_{B_0} \rho_0 b(t) \cdot \delta \varphi dV - \int_{\partial_\alpha B_0} \bar{t}(t) \cdot \delta \varphi dS \\ & + \int_{B_0} \delta t \rho_0 Q(t) \log \frac{T}{T_0} dV - \int_{\partial_H B_0} \delta t \bar{H}(t) \log \frac{T}{T_0} dS \end{aligned} \quad (67) $$

where the average conduction dissipation function $\langle \chi \rangle$ can be treated in a similar fashion as dissipation potential $\Delta$ occurring in the local incremental
problem. The incremental boundary-value problem thus takes a variational form:

$$ \{\boldsymbol{\varphi}, T\} = \arg \underset{\boldsymbol{\varphi}, T}{\mathrm{stat}} \mathcal{I}(\boldsymbol{\varphi}, T). \qquad (68) $$

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:

$$ \{\boldsymbol{\varphi}, T\} = \arg \inf_{\boldsymbol{\varphi}} \sup_{T} \mathcal{I}(\boldsymbol{\varphi}, T). \qquad (69) $$

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 $\mathcal{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.

### 3.2.3 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:

$$ h = h(T - T_{\mathrm{ext}})n \quad \text{for } x \in \partial_c \mathcal{B} \qquad (70) $$
where $\partial_c B$ is the part of the (deformed) boundary where heat exchange conditions are effective, $T_{\text{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 $\Phi_c$ to functional (62). This term is given by

$$ \Phi_c(T; F, \Theta) = \int_{\partial_c B_0} \frac{h}{2} \frac{(T - T_{\text{ext}}(t))^2}{\Theta} J_S \, dS \quad (71) $$

where $J_S$ is the ratio of a surface element in the deformed configuration to the same surface element in the reference configuration: $J_S = n \cdot \text{cot} F \cdot N = Jn \cdot F^{-T} \cdot N$. The variational principle then becomes

$$ \{\dot{\boldsymbol{\varphi}}, \dot{\eta}, T\} = \arg \underset{\dot{\boldsymbol{\varphi}}, \dot{\eta}, T}{\text{stat}} [\Phi(\dot{\boldsymbol{\varphi}}, \dot{\eta}, T) - \Phi_c(T)]. \quad (72) $$

In this modified variational principle, it is of course understood that $\partial_T B_0 \cup \partial_H B_0 \cup \partial_c B_0 = \partial B_0$ and $\partial_T B_0 \cap \partial_c B_0 = \emptyset$, $\partial_H B_0 \cap \partial_c B_0 = \emptyset$. 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:

$$ \mathcal{I}_c(T; F) = \int_{\partial_c B_0} \delta t \frac{h}{2} \frac{(T - T_{\text{ext}}(t))^2}{T_0} J_{S\alpha} \, dS \quad (73) $$

where $J_{S\alpha}$ denotes the surface Jacobian evaluated at time $t_\alpha = (1 - \alpha)t_0 + \alpha t$. The variational principle is then given by

$$ \{\boldsymbol{\varphi}, T\} = \arg \underset{\boldsymbol{\varphi}, T}{\text{stat}} [\mathcal{I}(\boldsymbol{\varphi}, T) - \mathcal{I}_c(T; \nabla_0 \boldsymbol{\varphi})]. \quad (74) $$

It is interesting to note that the inclusion of the term $\mathcal{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 $\Phi - \Phi_c$ with respect to $\boldsymbol{\varphi}$). This term is purely numerical, but ensures the symmetry of the mathematical formulation. It will vanish when $\delta t \to 0$, as consistency requires. In practice, this effect can be avoided by choosing $\alpha = 0$ in expression (73). The effect of deformation on the exchanged heat flux will then be delayed by one time step.
### 3.3. Dynamics

The variational formulation (68) accounts for rate-dependent behavior, including transient thermal effects. It formally corresponds to a quasi-stationary 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:

$$ \mathcal{I}(\boldsymbol{\varphi}, T) = \int_{\mathcal{B}_0} \left[ \frac{\rho_0}{2\beta\delta t^2} \delta\boldsymbol{\varphi} \cdot \delta\boldsymbol{\varphi} + W(\mathbf{F}, T) - \delta t \left\langle \chi\left(-\frac{\nabla_0 T}{T}; \mathbf{F}_\alpha, T_\alpha, \mathbf{Z}_\alpha\right) \right\rangle \right] dV \\ - \int_{\mathcal{B}_0} \rho_0 \boldsymbol{b}^*(t) \cdot \delta\boldsymbol{\varphi} dV - \int_{\partial_\sigma \mathcal{B}_0} \bar{\mathbf{t}}(t) \cdot \delta\boldsymbol{\varphi} dS \\ + \int_{\mathcal{B}_0} \Delta t \rho_0 Q(t) \log \frac{T}{T_0} dV - \int_{\partial_H \mathcal{B}_0} \delta t \bar{H}(t) \log \frac{T}{T_0} dS \quad (75) $$

with $\boldsymbol{b}^*(t) = \boldsymbol{b}(t) - \frac{1}{\beta\delta t}(\dot{\boldsymbol{\varphi}}_0 + (\frac{1}{2} - \beta)\delta t \ddot{\boldsymbol{\varphi}}_0)$. The incremental boundary-value problem then takes the variational form

$$ \{\boldsymbol{\varphi}, T\} = \arg\inf_{\boldsymbol{\varphi}} \sup_{T} \mathcal{I}(\boldsymbol{\varphi}, T). \quad (76) $$

When combined with the following update rule for accelerations and velocities:

$$ \ddot{\boldsymbol{\varphi}} = \frac{1}{\beta} - \beta \ddot{\boldsymbol{\varphi}}_0 + \frac{1}{\beta\delta t} \dot{\boldsymbol{\varphi}}_0 + \frac{\delta\boldsymbol{\varphi}}{\beta\delta t^2} \quad (77) $$

$$ \dot{\boldsymbol{\varphi}} = \dot{\boldsymbol{\varphi}}_0 + (1 - \gamma)\delta t \ddot{\boldsymbol{\varphi}}_0 + \gamma \delta t \ddot{\dot{\boldsymbol{\varphi}}} \quad (78) $$

the stationarity condition of $\mathcal{I}$ with respect to $\boldsymbol{\varphi}$ 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 ($\beta = 0.25$, $\gamma = 0.5$) in implicit dynamics (ensuring unconditional stability in linear elasto-dynamics), or with parameters ($\beta = 0.0$, $\gamma = 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, 2008), since they allow to work with an effective incremental potential formally identical to hyper-elasticity.

### 3.4. 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:

$$ \boldsymbol{\epsilon} = \frac{1}{2} (\nabla_0 \boldsymbol{\varphi} + (\nabla_0 \boldsymbol{\varphi})^T) - \mathbf{I}. \qquad (79) $$

where $\mathbf{I}$ is the identity second-order tensor. All the above variational principles still hold, provided that the deformation gradient is replaced by

$$ \mathbf{F} \approx \mathbf{I} + \boldsymbol{\epsilon} \qquad (80) $$

and Piola stress tensor $\mathbf{P}$ by Cauchy stress tensor $\boldsymbol{\sigma}$.

Similarly, if temperature variations are small (about a reference temperature $T_r$), i.e.:

$$ T = T_r + \theta, \quad \frac{\theta}{T_r} \ll 1 \qquad (81) $$

the problem can be simplified. The generalized gradient is then given by

$$ \mathbf{G} \approx -\frac{\nabla \theta}{T_r} \qquad (82) $$

while coefficients occurring in the time discrete variational formulation become

$$ \frac{T}{T_0} \approx 1 + \frac{\delta\theta}{T_r}. \qquad (83) $$

A consistent approximation to the average dissipation function (eq. 3.21) is then obtained by replacing factors as follows:

$$ \frac{T_0}{T} \approx \frac{T_r}{T_r + \delta\theta} \quad \text{and} \quad \frac{\delta T}{T} \approx \frac{\delta\theta}{T_r + \delta\theta}. \qquad (84) $$
Finally, the “linearized” balance equations become:

$$ \begin{cases} \nabla \cdot \boldsymbol{\sigma} + \rho_0 \boldsymbol{b}(t) = 0, \\ \rho_0 T_r \dot{\eta} = \mathcal{D}_{\text{int}} - \nabla \cdot \boldsymbol{h}(\nabla \theta) + \rho_0 Q(t). \end{cases} \qquad (85) $$

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.

# 4. 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.

## 4.1. Linearized kinematics

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

### 4.1.1 Kelvin-Voigt model

The simplest visco-elasticity model can be described by combining a linear (thermo-) elasticity free-energy potential:

$$ W(\boldsymbol{\varepsilon}, \theta) = \frac{1}{2} \boldsymbol{\varepsilon} \cdot \mathcal{C}(\theta) \cdot \boldsymbol{\varepsilon} - \theta \boldsymbol{\alpha} \cdot \mathcal{C}(\theta) \cdot \boldsymbol{\varepsilon} - \frac{1}{2} C \frac{\theta^2}{T_r} \qquad (86) $$

where $\mathcal{C}(\theta)$ is a (temperature-dependent) fourth-order elasticity tensor, $\boldsymbol{\alpha}$ a second-order tensor describing (possibly anisotropic) thermal expansion, $C$ the heat capacity (per unit undeformed volume), and a linear viscosity dissipation function:
$$ \Delta(\dot{\boldsymbol{\varepsilon}}; \theta) = \frac{1}{2} \dot{\boldsymbol{\varepsilon}} \cdot \mathbb{C}_v(\theta) \cdot \dot{\boldsymbol{\varepsilon}} \qquad (87) $$

where $\mathbb{C}_v(\theta)$ is a (temperature-dependent) fourth-order viscosity tensor. The entropy is then given by the following expression:

$$ \rho_0 \eta(\boldsymbol{\varepsilon}, \theta) = -\frac{\partial W}{\partial \theta} = C \frac{\theta}{T_r} + \boldsymbol{\alpha} \cdot \mathbb{C} \cdot \boldsymbol{\varepsilon} - \frac{1}{2} \boldsymbol{\varepsilon} \cdot \mathbb{C}'(\theta) \cdot \boldsymbol{\varepsilon} + \theta \boldsymbol{\alpha} \cdot \mathbb{C}'(\theta) \cdot \boldsymbol{\varepsilon} \quad (88) $$

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 thermo-mechanical 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

$$ \mathcal{W}(\boldsymbol{\varepsilon}, \theta; \boldsymbol{\varepsilon}_0, \theta_0) = W(\boldsymbol{\varepsilon}, \theta) - W(\boldsymbol{\varepsilon}_0, \theta_0) + \rho_0 \eta_0(\theta - \theta_0) \\ + \delta t \left\langle \Delta \left( \left( 1 + \frac{\delta \theta}{T_r} \right) \frac{\delta \boldsymbol{\varepsilon}}{\delta t} ; \theta(\tau) \right) \right\rangle. \qquad (89) $$

Using the consistent approximation of incremental dissipation functions proposed in Stainier (2011), the last term becomes

$$ \left\langle \Delta \left( \left( 1 + \frac{\delta \theta}{T_r} \right) \frac{\delta \boldsymbol{\varepsilon}}{\delta t} ; \theta \right) \right\rangle = \frac{T_r}{T_r + \delta \theta} \Delta \left( \left( 1 + \frac{\delta \theta}{T_r} \right) \frac{\delta \boldsymbol{\varepsilon}}{\delta t} ; \theta_0 \right) \\ + \frac{\delta \theta}{T_r + \delta \theta} \Delta \left( \left( 1 + \frac{\delta \theta}{T_r} \right) \frac{\delta \boldsymbol{\varepsilon}}{\delta t} ; \theta_\alpha \right) \quad (90) $$

Considering the specific quadratic form we chose for $\Delta$ in this setting, this expression reduces to

$$ \left\langle \Delta \left( \left( 1 + \frac{\delta \theta}{T_r} \right) \frac{\delta \boldsymbol{\varepsilon}}{\delta t} ; \theta \right) \right\rangle = \frac{1}{2} \left( 1 + \frac{\delta \theta}{T_r} \right) \frac{\delta \boldsymbol{\varepsilon}}{\delta t} \cdot \left[ \mathbb{C}_v(\theta_0) + \frac{\delta \theta}{T_r} \mathbb{C}_v(\theta_\alpha) \right] \cdot \frac{\delta \boldsymbol{\varepsilon}}{\delta t}. \qquad (91) $$

The resulting incremental expressions for stress and net entropy increase are then given by

$$ \boldsymbol{\sigma} = \frac{\partial W}{\partial \boldsymbol{\varepsilon}} = \mathbb{C}(\theta) \cdot (\boldsymbol{\varepsilon} - \theta \boldsymbol{\alpha}) + \left( 1 + \frac{\delta \theta}{T_r} \right) \left[ \mathbb{C}_v(\theta_0) + \frac{\delta \theta}{T_r} \mathbb{C}_v(\theta_\alpha) \right] \cdot \frac{\delta \boldsymbol{\varepsilon}}{\delta t} \quad (92) $$
$$ -\rho_0 \delta \eta^{\text{eff}} = \frac{\partial W}{\partial \theta} = -\rho_0 \delta \eta + \frac{\delta \boldsymbol{\epsilon}}{2T_r} \cdot \left[ \mathbb{C}_v(\theta_0) + \frac{\delta \theta}{T_r} \mathbb{C}_v(\theta_\alpha) \right. \\ \left. + \left(1 + \frac{\delta \theta}{T_r}\right) \left( \mathbb{C}_v(\theta_\alpha) + \alpha \delta \theta \mathbb{C}'_v(\theta_\alpha) \right) \right] \cdot \frac{\delta \boldsymbol{\epsilon}}{\delta t} \quad (93) $$

which are consistent with continuous expressions when $\delta t \to 0$.

### 4.1.2 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:

$$ W(\boldsymbol{\epsilon}, \theta, \boldsymbol{\epsilon}^v) = W^{(0)}(\boldsymbol{\epsilon}, \theta) + W^{(1)}(\boldsymbol{\epsilon}^e, \theta) \quad (94) $$

where $\boldsymbol{\epsilon}^e = \boldsymbol{\epsilon} - \boldsymbol{\epsilon}^v$, and

$$ \Delta(\dot{\boldsymbol{\epsilon}}^v; \theta) = \frac{1}{2} \dot{\boldsymbol{\epsilon}}^v \cdot \mathbb{C}_v(\theta) \cdot \dot{\boldsymbol{\epsilon}}^v. \quad (95) $$

Typically, $W^{(0)}$ can take a form similar to Eq. (86), while $W^{(1)}$ can take a simple quadratic form:

$$ W^{(1)}(\boldsymbol{\epsilon}^e, \theta) = \frac{1}{2} \boldsymbol{\epsilon}^e \cdot \mathbb{C}^{(1)}(\theta) \cdot \boldsymbol{\epsilon}^e \quad (96) $$

where $\mathbb{C}^{(1)}(\theta)$ 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:

$$ W(\boldsymbol{\epsilon}, \theta; \boldsymbol{\epsilon}_0, \theta_0) = \inf_{\boldsymbol{\epsilon}^v} \left[ W^{(0)}(\boldsymbol{\epsilon}, \theta) + W^{(1)}(\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^v, \theta) - W(\boldsymbol{\epsilon}_0, \theta_0, \boldsymbol{\epsilon}_0^v) \right. \\ \left. + \rho_0 \eta_0 \delta \theta + \delta t \Delta \left( \left( 1 + \frac{\delta \theta}{T_r} \right) \frac{\delta \boldsymbol{\epsilon}^v}{\delta t} \right) \right]. \quad (97) $$
Stationarity condition for the viscous strain yields:

$$-\mathbb{C}^{(1)} \cdot (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{\nu}) + \left(1 + \frac{\delta\theta}{T_r}\right)^2 \mathbb{C}_v \cdot \frac{\delta\boldsymbol{\varepsilon}^{\nu}}{\delta t} = 0 \text{ or} \\ \boldsymbol{\varepsilon}^{\nu} = \left[\mathbb{C}^{(1)} + \frac{1}{\delta t} \left(1 + \frac{\delta\theta}{T_r}\right)^2 \mathbb{C}_v\right]^{-1} \cdot \left[\mathbb{C}^{(1)} \cdot \boldsymbol{\varepsilon} + \frac{1}{\delta t} \left(1 + \frac{\delta\theta}{T_r}\right)^2 \mathbb{C}_v \cdot \boldsymbol{\varepsilon}_0^{\nu}\right]. \qquad (98)$$

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:

$$\boldsymbol{\sigma} = \frac{\partial W}{\partial \boldsymbol{\varepsilon}} = \boldsymbol{\sigma}^{(0)}(\boldsymbol{\varepsilon}, \theta) + \mathbb{C}^{(1)} \cdot (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{\nu}) \qquad (99) \\ -\rho_0 \delta \eta^{\text{eff}} = \frac{\partial W}{\partial \theta} = \rho_0 \delta \eta + \frac{1}{T_r} \left(1 + \frac{\delta \theta}{T_r}\right) \delta \boldsymbol{\varepsilon}^{\nu} \cdot \mathbb{C}_v \cdot \frac{\delta \boldsymbol{\varepsilon}^{\nu}}{\delta t} \qquad (100)$$

where $\boldsymbol{\sigma}^{(0)}(\boldsymbol{\varepsilon}, \theta) = \partial_{\boldsymbol{\varepsilon}} W^{(0)}(\boldsymbol{\varepsilon}, \theta)$ and $\boldsymbol{\varepsilon}^{\nu}$ is given by Eq. (98). Once again, these results are consistent with continuous expressions when $\delta t \to 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.

## 4.2. 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.

### 4.2.1 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:

$$W(\mathbf{F}, T) = f(J, T) + W^e(\hat{\mathbf{C}}, T) + W^t(T) \qquad (101)$$

where $J = \det[\mathbf{F}]$, $\hat{\mathbf{C}} = \hat{\mathbf{F}}^T \hat{\mathbf{F}}$, $\hat{\mathbf{F}} = J^{-\frac{1}{3}} \mathbf{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

$$f(J, T) = K[J - 1 - \log J - 3\alpha(T - T_r)\log J] \quad (102)$$

where $K$ is the bulk modulus and $\alpha$ the coefficient of thermal expansion, while for heat capacity, a typical example is

$$W^t(T) = C \left[ T - T_r - T \log \frac{T}{T_r} \right] \quad (103)$$

where $C$ is the heat capacity (per unit undeformed volume). For the strain-energy 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

$$\Delta(\dot{\mathbf{F}}; \mathbf{F}, T) = \phi(\mathbf{D}; T). \quad (104)$$

This specific form implies that no viscous stress is generated by pure rotations (objectivity). The simplest case would then be given by

$$\phi(\mathbf{D}; T) = \frac{1}{2} \mathbf{D} \cdot \mathcal{C}_v(T) \cdot \mathbf{D}. \quad (105)$$

Note that in an incremental context, one must be attentive to preserving objectivity. For example, if approximating $\mathbf{L}$ by $\delta\mathbf{F}\mathbf{F}_\alpha^{-1}$, rigid body rotations will generate spurious strain rate. Suitable approximations can be built by noting that $\mathbf{D} = \mathbf{F}^{-T}\dot{\mathbf{E}}\mathbf{F}^{-1}$, yielding

$$\mathbf{D} \approx \frac{1}{2\delta t} \mathbf{F}_\alpha^{-T} \delta \mathbf{C} \mathbf{F}_\alpha^{-1} \quad (106)$$

which clearly becomes equal to zero in case of rigid-body motions. An alternative approximation, which preserves isochoricity (i.e. if $\det \mathbf{F} = \det \mathbf{F}_0$, then $\mathrm{tr} \mathbf{D} = 0$), is given by:

$$\mathbf{D} \approx \frac{1}{2\delta t} \log(\mathbf{F}_0^{-T} \mathbf{C} \mathbf{F}_0^{-1}) \quad (107)$$

The incremental potential energy then takes the form

$$\mathcal{W}(\mathbf{F}, T; \mathbf{F}_0, T_0) = W(\mathbf{F}, T) - W(\mathbf{F}_0, T_0) + \rho_0 \eta_0 \delta T \\ + \delta t \left\langle \phi \left( \frac{T}{T_0} \mathbf{D}; T(\tau) \right) \right\rangle \quad (108)$$
where a consistent approximation to the average dissipation function is provided by

$$ \left\langle \phi \left( \frac{T}{T_0} \mathbf{D}; T(\tau) \right) \right\rangle \approx \frac{T_0}{T} \phi \left( \frac{T}{T_0} \mathbf{D}; T_0 \right) + \frac{\delta T}{T} \phi \left( \frac{T}{T_0} \mathbf{D}; T_\alpha \right). \quad (109) $$

Considering the simple quadratic form earlier, it yields

$$ \left\langle \phi \left( \frac{T}{T_0} \mathbf{D}; T(\tau) \right) \right\rangle = \frac{1}{2} \frac{T}{T_0} \mathbf{D} \cdot \left( \mathbb{C}_v(T_0) + \frac{\delta T}{T_0} \mathbb{C}_v(T_\alpha) \right) \cdot \mathbf{D} \quad (110) $$

where **D** is given by one of the above incremental approximations.

The stress tensor is then given by

$$ \mathbf{P} = \frac{\partial W}{\partial \mathbf{F}} = \frac{\partial W}{\partial \mathbf{F}} + \mathbf{F} \left[ \frac{T}{T_0} \mathbf{D} \cdot \left( \mathbb{C}_v(T_0) + \frac{\delta T}{T_0} \mathbb{C}_v(T_\alpha) \right) \cdot \left( 2\delta t \frac{\partial \mathbf{D}}{\partial \mathbf{C}} \right) \right] \quad (111) $$

while the net effective entropy variation is given by

$$ -\rho_0 \delta \eta^{\text{eff}} = \frac{\partial W}{\partial T} = -\rho_0 \delta \eta + \frac{\delta t}{2T_0} \mathbf{D} \cdot \left[ \mathbb{C}_v(T_0) + \frac{\delta T}{T_0} \mathbb{C}_v(T_\alpha) \right. \\ \left. + \frac{T}{T_0} \left( \mathbb{C}_v(T_\alpha) + \alpha \delta T \mathbb{C}'_v(T_\alpha) \right) \right] \cdot \mathbf{D}. \quad (112) $$

It is relatively straightforward to verify that these expressions are consistent with the continuous formulation when $\delta t \to 0$.

### 4.2.2 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:

$$ W(\mathbf{F}, T, \mathbf{F}^v) = W^{(0)}(\mathbf{C}, T) + W^{(1)}(\mathbf{C}^e, T) \quad (113) $$

where we adopt the multiplicative decomposition initially proposed by Sidoroff (1974) for finite visco-elasticity:

$$ \mathbf{F} = \mathbf{F}^e \mathbf{F}^v. \quad (114) $$
We recall that free-energy potentials should depend on $C = F^T F$ and $C^e = F^e F$ in order to ensure objectivity. Viscous deformations are assumed to produce no rotation:

$$ \dot{F}^v = D^v F^v \qquad (115) $$

where $D^v$ is a viscous strain-rate symmetric tensor. Viscous strains are typically assumed to be isochoric, in which case tr[$D^v$] = 0. In the simplest case (one Maxwell branch in parallel with an elastic branch), the dissipation function then takes the form:

$$ \Delta(\dot{F}^v; F^v, T) = \phi^{(1)}(D^v; T). \qquad (116) $$

In order to avoid lengthy developments, we will assume in the following that there is no parametric dependence of $\phi^{(1)}$ on $T$.

The incremental energy potential then takes the form:

$$ W(F, T; F_0, F_0^v, T_0) = \inf_{D^v} \left[ W^{(0)}(C, T) + W^{(1)}(C^e, T) - W(F_0, T_0, F_0^v) \right. \\ \left. + \rho_0 \eta_0 \delta T + \delta t \phi \left( \frac{T}{T_0} D^v \right) \right] \qquad (117) $$

where one must provide an incremental update rule for $F^v$. The exponential update formula of Weber and Anand (1990), initially proposed in the framework of finite elasto-plasticity:

$$ F^v = \exp(\delta t D^v) F_0^v \qquad (118) $$

provides the advantage of preserving the isochoric nature of viscous deformation (if enforced) and will thus be used here. We can then write

$$ C^e = \exp(-\delta t D^v) C_{tr}^e \exp(-\delta t D^v) \qquad (119) $$

where $C_{tr}^e = F_0^{v-T} C_0^e F_0^{v-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:

$$ W^{(1)}(C^e, T) = W^{(1)}(c_1^e, c_2^e, c_3^e, T) \qquad (120) $$

$$ \phi^{(1)}(D^v) = \phi^{(1)}(d_1^v, d_2^v, d_3^v) \qquad (121) $$

where $\{c_1^e, c_2^e, c_3^e\}$ (resp. $\{d_1^v, d_2^v, d_3^v\}$) are the principal values of tensor $C^e$ (resp. $D^v$), the optimization problem of the constitutive update can be reduced to a
simpler form. Indeed the ansatz that $D^ν$, $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(\lambda_I^e) = \log(\lambda_I^{tr}) - \delta t d_I^ν \quad (I = 1, 2, 3) \quad (122) $$

where $\lambda_I^e = \sqrt{c_I^e}$ and $\lambda_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 $\varepsilon_I^e = \log(\lambda_I^e)$ and $\varepsilon_I^{tr} = \log(\lambda_I^{tr})$. The optimization problem now reduces to

$$ \inf_{d_1^ν, d_2^ν, d_3^ν} \left[ W^{(1)}(\varepsilon_1^e, \varepsilon_2^e, \varepsilon_3^e, T) + \delta t \phi^{(1)} \left( \frac{T}{T_0} d_1^ν, \frac{T}{T_0} d_2^ν, \frac{T}{T_0} d_3^ν \right) \right] \quad (123) $$

where $\varepsilon_I^e = \varepsilon_I^{tr} - \delta t d_I^ν$, and where an isochoricity constraint can be added: $d_1^ν + d_2^ν + d_3^ν = 0$. 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

$$ \rho_0 \delta \eta^{\text{eff}} = \frac{\partial W}{\partial T} = -\rho_0 \delta \eta + \frac{\delta t}{T_0} \frac{\partial \phi^{(1)}}{\partial D^ν} \left( \frac{T}{T_0} D^ν \right) \cdot D^ν \quad (124) $$

which is consistent with the continuous formulation when $\delta t \to 0$.

### 4.2.3 Viscous fluids

Newtonian fluids are described by Navier-Poisson constitutive equations:

$$ \boldsymbol{\sigma} = -p(\rho, T)\boldsymbol{I} + \kappa \text{tr}[\boldsymbol{D}]\boldsymbol{I} + 2\mu_ν \text{dev} \boldsymbol{D} \quad (125) $$

where $\boldsymbol{\sigma}$ is Cauchy stress tensor, $\boldsymbol{D} = (\boldsymbol{L} + \boldsymbol{L}^T)/2$ is the strain rate tensor, and $p(J, T)$ is the hydrostatic pressure, related to the density $\rho$ and temperature $T$ through the equation of state (EOS). It is usually assumed that, either the fluid flow is incompressible ($\text{tr}[\boldsymbol{D}] = 0$), either $\kappa = 0$ (Stokes condition), such that the pressure is always equal to the hydrostatic pressure: $\text{tr}[\boldsymbol{\sigma}] = -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:

$$ W(\boldsymbol{F}, T) = f(J, T) \quad (126) $$
where $J = \det F$, and a viscous dissipation potential of the form:

$$ \Delta(\dot{F}; F, T) = \mu_v J \mathrm{dev}[D] \cdot \mathrm{dev}[D]. \qquad (127) $$

It is easily verified that this yields Navier-Poisson equations, with $\kappa = 0$ and

$$ p(\rho, T) = -\frac{\partial f}{\partial J}(J, T) = -\frac{\partial f}{\partial J}\left(\frac{\rho_0}{\rho}, T\right) \qquad (128) $$

where we recall that $J\rho = \rho_0$.

Finally, note that most non-Newtonian viscous fluid models could also be formulated in the current variational framework.

# 5. 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.

## 5.1. Crystal plasticity

### 5.1.1 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

$$ L^p = \dot{F}^p F^{p-1} = \sum_{k=1}^{N} \dot{\gamma}^{(k)} s^{(k)} \otimes m^{(k)} \qquad (129) $$
where $\gamma^{(k)}$ is the slip strain, and $\mathbf{s}^{(k)}$, and $\mathbf{m}^{(k)}$ are orthogonal unit vectors defining the slip direction and slip-plane normal corresponding to slip system $k$. The collection $\gamma$ of slip strains may be regarded as a subset of the internal variable set $\mathbf{Z}$. A zero value of a slip rate $\dot{\gamma}^{(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 $\{\mathbf{s}^{(k)}, \mathbf{m}^{(k)}\}$ remain constant throughout the deformation and are determined by crystallography. In order to account for non-monotonous loading paths, it is common to consider dislocation glide in $+\mathbf{s}^{(k)}$ and $-\mathbf{s}^{(k)}$ directions as occurring in separate systems, and adding the constraint $\dot{\gamma}^{(k)} \ge 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 $\mathbf{F}^e$, or elastic deformation, may also be expected in general. One therefore has, locally,

$$ \mathbf{F} = \mathbf{F}^e \mathbf{F}^p \qquad (130) $$

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:

$$ W(\mathbf{F}, T, \mathbf{F}^p, \gamma) = W^e(\mathbf{C}^e, T) + W^p(\gamma, T) + W^t(T) \qquad (131) $$

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 $\mathbf{C}^e$ (instead of $\mathbf{F}^e$) ensures objectivity, as explained earlier. The rate of free energy can then be written as

$$ \dot{W} = \mathbf{P} \cdot \dot{\mathbf{F}} - (\boldsymbol{\tau} - \mathbf{g}) \cdot \dot{\boldsymbol{\gamma}} - \rho_0 \eta \dot{T} \qquad (132) $$

where $\boldsymbol{\tau} = \{\tau^{(k)}\}$ is the collection of resolved shear stresses:

$$ \tau^{(k)} = (\mathbf{F}^e T \mathbf{P} \mathbf{F}^p T) \cdot (\mathbf{s}^{(k)} \otimes \mathbf{m}^{(k)}) \qquad (133) $$

while $\mathbf{g} = \{g^{(k)}\}$ is the collection of yield resolved shear stresses linked to plastic energy storage mechanisms $g^{(k)} = \partial_{\gamma^{(k)}} W^p$.
The kinematic of plastic slip can be modeled through a dissipation function $\Delta(\dot{\gamma}; \gamma, T)$. Typical expressions include power-law type formulas such as

$$ \Delta(\dot{\gamma}; \gamma, T) = \sum_{k=1}^{N} \phi_{(k)}(\dot{\gamma}^{(k)}; \gamma, T) $$

$$ \text{with } \phi_{(k)}(\dot{\gamma}^{(k)}; \gamma, T) = \begin{cases} \frac{m}{m+1} \gamma^{(k)}(\gamma, T) \dot{\gamma}_{0}^{(k)} \left( \frac{\dot{\gamma}^{(k)}}{\dot{\gamma}_{0}^{(k)}} \right)^{\frac{m+1}{m}} & \text{if } \dot{\gamma}^{(k)} \ge 0 \\ +\infty & \text{otherwise} \end{cases} \quad (134) $$

where $\mathbf{Y} = \{Y^{(k)}\}$ is a collection of dissipative yield resolved shear stresses and $\dot{\gamma}_{0}^{(k)}$ are reference slip rates. The exponent $m \in [1, +\infty]$ controls rate-dependency effects. In particular, rate-independent behavior can be recovered at the limit when $m \to +\infty$. 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)}(\gamma, T)$ (see for example Stainier, Cuitiño, & Ortiz, 2002). More complex rate-dependency 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

$$ D(\dot{\mathbf{F}}, \dot{\eta}, T, \dot{\gamma}, \mathbf{F}, \eta, \mathbf{F}^{\mathrm{p}}, \gamma) = \mathbf{P}\dot{\mathbf{F}} - (\boldsymbol{\tau} - \mathbf{g}) \cdot \dot{\gamma} - \rho_{0}\dot{\eta}(T - \boldsymbol{\Theta}) + \Delta\left(\frac{T}{\boldsymbol{\Theta}}\dot{\gamma}; \gamma, \boldsymbol{\Theta}\right) \quad (135) $$

and the minimization problem

$$ D_{\mathrm{eff}}(\dot{\mathbf{F}}, \dot{\eta}, T) = \inf_{\dot{\gamma}} D(\dot{\mathbf{F}}, \dot{\eta}, T, \dot{\gamma}; \mathbf{F}, \eta, \mathbf{F}^{\mathrm{p}}, \gamma). \quad (136) $$

Considering for example the case of the power-law dissipation function given in Eq. (134), stationarity conditions yield

$$ \tau^{(k)} = g^{(k)} + Y^{(k)}(\gamma, T) \left( \frac{\dot{\gamma}^{(k)}}{\dot{\gamma}_{0}^{(k)}} \right)^{\frac{1}{m}} \quad (137) $$

where we have accounted for the fact that posterior stationarity condition with respect to $\dot{\eta}$ will yield $T = \boldsymbol{\Theta}$.
### 5.1.2 Incremental update

In the time-discrete setting, the local variational constitutive update for crystal plasticity takes the following form:

$$ W(F, T) = \inf_{\gamma} \left[ W(F, T, F^p, \gamma) + \delta t \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta \gamma}{\delta t} ; \gamma(\tau), T(\tau) \right) \right\rangle \right. \\ \left. + \rho_0 \eta_0 \delta T - W(F_0, T_0, F_0^p, \gamma_0). \right] \quad (138) $$

In the above expression, a relation must be provided between $F^p$, $\delta\gamma$, 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):

$$ F^p = \exp \left( \sum_{k=1}^{N} \delta \gamma^{(k)} s^{(k)} \otimes m^{(k)} \right) F_0^p \quad (139) $$

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:

$$ \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta \gamma}{\delta t} ; \gamma(\tau), T(\tau) \right) \right\rangle \approx \Delta \left( \frac{T}{T_0} \frac{\delta \gamma}{\delta t} ; \gamma_{\alpha} \right). \quad (140) $$

The constitutive update problem then reduces to the following constrained optimization problem:

$$ \gamma = \arg \inf_{\gamma} \left[ W(F, T, F^p, \gamma) + \delta t \Delta \left( \frac{T}{T_0} \frac{\delta \gamma}{\delta t} ; \gamma_{\alpha} \right) \right] \quad \text{with } \delta \gamma^{(k)} \ge 0 (k = 1, \dots, N) \quad (141) $$

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.
## 5.2. 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).

### 5.2.1 Linear kinematics

Under small strains assumptions, the total deformation is additively decomposed into an elastic strain $\boldsymbol{\epsilon}^e$ and a plastic strain $\boldsymbol{\epsilon}^p$:

$$ \boldsymbol{\epsilon} = \boldsymbol{\epsilon}^e + \boldsymbol{\epsilon}^p. \qquad (142) $$

A plastic flow rule corresponding to von Mises-type plasticity can be written as follows:

$$ \dot{\boldsymbol{\epsilon}}^{\mathrm{p}}=\dot{\bar{\boldsymbol{\epsilon}}}^{\mathrm{p}} \boldsymbol{M} \qquad (143) $$

where $\dot{\bar{\boldsymbol{\epsilon}}}^p$ is the amplitude of the plastic strain rate and $\boldsymbol{M}$ its direction. The symmetric tensor $\boldsymbol{M}$ must of course be normalized, and the particular choice

$$ \boldsymbol{M} \cdot \boldsymbol{M}=\frac{3}{2} \qquad (144) $$

ensures that $\bar{\boldsymbol{\epsilon}}^p$ corresponds to the cumulated equivalent plastic strain. In addition, we will require that plastic deformation be isochoric:

$$ \mathrm{tr}[\boldsymbol{M}] = 0. \qquad (145) $$

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:

$$ W(\boldsymbol{\epsilon}, T, \boldsymbol{\epsilon}^{\mathrm{p}}, \bar{\boldsymbol{\epsilon}}^{\mathrm{p}}) = W^{\mathrm{e}}(\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^{\mathrm{p}}, T) + W^{\mathrm{p}}(\boldsymbol{\epsilon}^{\mathrm{p}}, \bar{\boldsymbol{\epsilon}}^{\mathrm{p}}, T) + W^{\mathrm{t}}(T) \qquad (146) $$

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

$$ \dot{W}=\boldsymbol{\sigma} \cdot \dot{\boldsymbol{\epsilon}}-\underbrace{\left[\left(\boldsymbol{\sigma}-\boldsymbol{\sigma}^{\mathrm{c}}\right) \cdot \boldsymbol{M}-\tau\right]}_{\equiv \gamma} \dot{\bar{\boldsymbol{\epsilon}}}^{\mathrm{p}}-\rho_{0} \eta \dot{T} \qquad (147) $$
where $\boldsymbol{\sigma} = \partial_{\boldsymbol{\epsilon}} W$ is the (Cauchy) stress, $\boldsymbol{\sigma}^c = \partial_{\boldsymbol{\epsilon}^p} W^p$ is the back-stress (kinematic hardening), and $g = \partial_{\bar{\boldsymbol{\epsilon}}^p} W^p$ is a yield stress associated to energy storage mechanisms (such as dislocation microstructures). From the above, we see that quantity $\gamma$ can be considered as conjugate to the cumulated equivalent plastic strain $\bar{\boldsymbol{\epsilon}}^p$, when accounting for the flow rule. Following the formalism of GSM, we will then describe the relation between $\gamma$ and $\dot{\bar{\boldsymbol{\epsilon}}}^p$ through the definition of a convex dissipation potential $\Delta(\dot{\bar{\boldsymbol{\epsilon}}}^p; \bar{\boldsymbol{\epsilon}}^p, T)$, or its conjugate $\Delta^*(\gamma; \bar{\boldsymbol{\epsilon}}^p, T)$.

### 5.2.1.1 Rate problem of visco-plasticity

Functional $D(\dot{\boldsymbol{\epsilon}}, \dot{\eta}, T, \dot{\bar{\boldsymbol{\epsilon}}}^p, \boldsymbol{M})$ is then given by

$$D(\dot{\boldsymbol{\epsilon}}, \dot{\eta}, T, \dot{\bar{\boldsymbol{\epsilon}}}^p, \boldsymbol{M}; \boldsymbol{\epsilon}, \eta, \boldsymbol{\epsilon}^p, \bar{\boldsymbol{\epsilon}}^p) = \boldsymbol{\sigma} \cdot \dot{\boldsymbol{\epsilon}} - \gamma \dot{\bar{\boldsymbol{\epsilon}}}^p - \rho_0 \dot{\eta} (T - \boldsymbol{\Theta}) + \Delta \left( \frac{T}{\boldsymbol{\Theta}} \dot{\bar{\boldsymbol{\epsilon}}}^p; \bar{\boldsymbol{\epsilon}}^p, \boldsymbol{\Theta} \right) \quad (148)$$

and the effective rate potential is defined as

$$D_{\text{eff}}(\dot{\boldsymbol{\epsilon}}, \dot{\eta}, T) = \inf_{\substack{\dot{\bar{\boldsymbol{\epsilon}}}^p \\ \boldsymbol{M}}} D(\dot{\boldsymbol{\epsilon}}, \dot{\eta}, T, \dot{\bar{\boldsymbol{\epsilon}}}^p, \boldsymbol{M}; \boldsymbol{\epsilon}, \eta, \boldsymbol{\epsilon}^p, \bar{\boldsymbol{\epsilon}}^p). \quad (149)$$

Note that the minimization with respect to $\boldsymbol{M}$ can be related to the principle of maximal dissipation:

$$\inf_{\boldsymbol{M}} D \Leftrightarrow \sup_{\boldsymbol{M}} \gamma \dot{\bar{\boldsymbol{\epsilon}}}^p. \quad (150)$$

Accounting for constraints on $\boldsymbol{M}$, the optimal plastic flow direction is given by

$$\boldsymbol{M} = \sqrt{\frac{3}{2} \frac{\text{dev}[\boldsymbol{\sigma} - \boldsymbol{\sigma}^c]}{\|\text{dev}[\boldsymbol{\sigma} - \boldsymbol{\sigma}^c]\|}} = \sqrt{\frac{3}{2} \frac{\boldsymbol{s}}{\|\boldsymbol{s}\|}} \quad (151)$$

where $\boldsymbol{s} = \text{dev}[\boldsymbol{\sigma} - \boldsymbol{\sigma}^c]$, which corresponds to the normal direction to von Mises yield criterion. Using this result, we can rewrite conjugate force $\gamma$ as

$$\gamma = \sqrt{\frac{3}{2}} \boldsymbol{s} \cdot \boldsymbol{s} - g = \boldsymbol{\sigma}^{\text{eq}} - g \quad (152)$$

where $\boldsymbol{\sigma}^{\text{eq}}$ is von Mises equivalent stress (accounting for backstress $\boldsymbol{\sigma}^c$). Dissipation potential $\Delta(\dot{\bar{\boldsymbol{\epsilon}}}^p; \bar{\boldsymbol{\epsilon}}^p, T)$ takes the general form

$$\Delta(\dot{\bar{\boldsymbol{\epsilon}}}^p; \bar{\boldsymbol{\epsilon}}^p, T) = \begin{cases} \phi(\dot{\bar{\boldsymbol{\epsilon}}}^p; \bar{\boldsymbol{\epsilon}}^p, T) & \text{if } \dot{\bar{\boldsymbol{\epsilon}}}^p \ge 0 \\ +\infty & \text{otherwise} \end{cases} \quad (153)$$
with $\phi$ a convex function such that $\phi(0; \bar{\epsilon}^p, T) = 0$ and $\partial_{\dot{\bar{\epsilon}}^p} \phi(0; \bar{\epsilon}^p, T) = \sigma_y(\bar{\epsilon}^p, T) \ge 0$. Minimization of $D$ with respect to $\dot{\bar{\epsilon}}^p$ then yields the following result:

$$ \arg\inf_{\dot{\bar{\epsilon}}^p} D(\dot{\boldsymbol{\epsilon}}, \dot{\eta}, T, \dot{\bar{\epsilon}}^p, \boldsymbol{M}; \boldsymbol{\epsilon}, \eta, \boldsymbol{\epsilon}^p, \bar{\epsilon}^p) = \begin{cases} 0 & \text{if } \sigma^{\text{eq}} - g < \sigma_y \\ \dot{\bar{\epsilon}}^{p*} : \sigma^{\text{eq}} - g = \partial_{\dot{\bar{\epsilon}}^p} \phi(\dot{\bar{\epsilon}}^{p*}; \bar{\epsilon}^p, T) & \text{if } \sigma^{\text{eq}} - g \ge \sigma_y. \end{cases} \quad (154) $$

where we have accounted for the fact that stationarity with respect to $\dot{\eta}$ will yield $T = \Theta$. Note that in the rate-independent case, corresponding to a function $\phi$ linear in $\dot{\bar{\epsilon}}^p$, we have $\partial_{\dot{\bar{\epsilon}}^p} \phi(\dot{\bar{\epsilon}}^p; \bar{\epsilon}^p, T) = \sigma_y(\bar{\epsilon}^p, T)$, and thus necessarily $\sigma^{\text{eq}} - g = \sigma_y$. We thus recover the classical results of von Mises elasto-visco-plasticity.

### 5.2.1.2 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:

$$ W(\boldsymbol{\epsilon}, \theta) = \inf_{\bar{\epsilon}^p, M} \left[ W(\boldsymbol{\epsilon}, \theta, \boldsymbol{\epsilon}^p, \bar{\epsilon}^p) + \delta t \left\langle \Delta \left( \left( 1 + \frac{\delta \theta}{T_r} \right) \frac{\delta \bar{\epsilon}^p}{\delta t} ; \bar{\epsilon}^p(\tau), \theta(\tau) \right) \right\rangle \right] \\ + \rho_0 \eta_0 \delta \theta - W(\boldsymbol{\epsilon}_0, \theta_0, \boldsymbol{\epsilon}_0^p, \bar{\epsilon}_0^p) \quad (155) $$

where

$$ \boldsymbol{\epsilon}^p = \boldsymbol{\epsilon}_0^p + \delta \bar{\epsilon}^p \boldsymbol{M}. \quad (156) $$

Consistent expressions for the average dissipation function $\langle \Delta \rangle$ 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 $\Delta$ on temperature, in which case we can take

$$ \left\langle \Delta \left( \left( 1 + \frac{\delta \theta}{T_r} \right) \frac{\delta \bar{\epsilon}^p}{\delta t} ; \bar{\epsilon}^p(\tau) \right) \right\rangle \approx \Delta \left( \left( 1 + \frac{\delta \theta}{T_r} \right) \frac{\delta \bar{\epsilon}^p}{\delta t} ; \bar{\epsilon}_\alpha^p \right). \quad (157) $$

Stationarity condition for $\boldsymbol{M}$ (including constraints) yields relation (151) again, where $\boldsymbol{s} = \text{dev}[\boldsymbol{\sigma} - \boldsymbol{\sigma}_c]$ depends on $\bar{\epsilon}^p$. Since dissipation function
$\Delta(\dot{\bar{\epsilon}}^p)$ is not regular, a practical approach to finding the infimum with respect to $\bar{\epsilon}^p$ consists at first evaluating the gradient of the incremental energy functional at the singularity point:

$$ \left( \frac{\partial}{\partial \delta \bar{\epsilon}^p} \left[ W(\boldsymbol{\epsilon}, \theta, \boldsymbol{\epsilon}^p, \bar{\epsilon}^p) + \delta t \Delta \left( \left( 1 + \frac{\delta \theta}{T_r} \right) \frac{\delta \bar{\epsilon}^p}{\delta t} ; \bar{\epsilon}_\alpha^p \right) \right] \right)_{\delta \bar{\epsilon}^p = 0^+} \\ = -\sigma_{\text{tr}}^{\text{eq}} + g_0 + \underbrace{(1 + \frac{\delta \theta}{T_r}) \partial_{\dot{\bar{\epsilon}}^p} \phi(0; \bar{\epsilon}_0^p) + \alpha \delta t \partial_{\bar{\epsilon}^p} \phi(0; \bar{\epsilon}_0^p)}_{\tilde{\sigma}_{y, \text{tr}}} \quad (158) $$

where $\sigma_{\text{tr}}^{\text{eq}}$ is the (trial) equivalent von Mises stress computed at material state $\{\boldsymbol{\epsilon}, \theta, \boldsymbol{\epsilon}_0^p, \bar{\epsilon}_0^p\}$, and where we have used the stationarity condition on $M$. Considering the convexity of the incremental energy (w.r.t. $\bar{\epsilon}^p$), we know that if this gradient is positive (i.e. $\sigma_{\text{tr}}^{\text{eq}} \le g_0 + \tilde{\sigma}_{y, \text{tr}}$), then the optimum is at $\delta \bar{\epsilon}^p = 0$, and the increment is elastic, while otherwise the optimum is $\delta \bar{\epsilon}^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

$$ W^e(\boldsymbol{\epsilon}^e, \theta) = f(\text{tr}[\boldsymbol{\epsilon}^e], \theta) + \mu(\theta) (\text{dev}[\boldsymbol{\epsilon}^e])^2 \quad (159) $$

where $f$ can be a quadratic function of $\text{tr}[\boldsymbol{\epsilon}^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 ($\boldsymbol{\sigma}_c \equiv 0$), we can then write

$$ \boldsymbol{\epsilon}^e = \boldsymbol{\epsilon}_{\text{tr}}^e - \delta \bar{\epsilon}^p \boldsymbol{M} \Rightarrow \boldsymbol{s} = \boldsymbol{s}_{\text{tr}} - 2\mu(\theta) \delta \bar{\epsilon}^p \boldsymbol{M} \quad (160) $$

where $\boldsymbol{\epsilon}_{\text{tr}}^e = \boldsymbol{\epsilon} - \boldsymbol{\epsilon}_0^p$ and $\boldsymbol{s}_{\text{tr}} = 2\mu \text{dev}[\boldsymbol{\epsilon}_{\text{tr}}^e]$. Considering stationarity condition (151), an immediate consequence is that $\boldsymbol{s}$, $\boldsymbol{s}_{\text{tr}}$, and $\boldsymbol{M}$ are aligned:

$$ \boldsymbol{M} = \sqrt{\frac{3}{2} \frac{\boldsymbol{s}}{\|\boldsymbol{s}\|}} = \sqrt{\frac{3}{2} \frac{\boldsymbol{s}_{\text{tr}}}{\|\boldsymbol{s}_{\text{tr}}\|}} \quad (161) $$

Thus, in the case of purely isotropic elasto-visco-plasticity, the variational formulation is equivalent to the classical radial return algorithm (Wilkins, 1964):

$$ \boldsymbol{s} = \left( 1 - \sqrt{\frac{3}{2} \frac{\delta \bar{\epsilon}^p}{\| \text{dev}[\boldsymbol{\epsilon}_{\text{tr}}^e] \|}} \right) \boldsymbol{s}_{\text{tr}} \quad (162) $$
The equivalent von Mises stress is then given by

$$ \sigma^{\text{eq}} = \sigma_{\text{tr}}^{\text{eq}} - 3\mu(\theta)\delta\bar{\epsilon}^{\text{p}} \quad (163) $$

and the stationarity condition for $\bar{\epsilon}^{\text{p}}$ is

$$ -\sigma_{\text{tr}}^{\text{eq}} + 3\mu(\theta)\delta\bar{\epsilon}^{\text{p}} + g(\bar{\epsilon}^{\text{p}}, \theta) + \bar{\sigma}_{\gamma}(\bar{\epsilon}^{\text{p}}, \theta) = 0 \quad (164) $$

where

$$ \bar{\sigma}_{\gamma} = \left(1 + \frac{\delta\theta}{T_r}\right) \frac{\partial\phi}{\partial\dot{\bar{\epsilon}}^{\text{p}}} \left(\left(1 + \frac{\delta\theta}{T_r}\right) \frac{\delta\bar{\epsilon}^{\text{p}}}{\delta t}; \bar{\epsilon}_{\alpha}^{\text{p}}\right) + \alpha\delta t \frac{\partial\phi}{\partial\bar{\epsilon}^{\text{p}}} \left(\left(1 + \frac{\delta\theta}{T_r}\right) \frac{\delta\bar{\epsilon}^{\text{p}}}{\delta t}; \bar{\epsilon}_{\alpha}^{\text{p}}\right) \quad (165) $$

is a consistent incremental approximation of the dissipative part of the yield stress. A study of the effect of the choice of algorithmic parameter $\alpha$ in the context of elasto-visco-plasticity can be found in Brassart and Stainier (2012).

### 5.2.2 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)]:

$$ F = F^{\text{e}} F^{\text{p}}. \quad (166) $$

A plastic flow rule corresponding to von Mises-type plasticity can then be written as follows:

$$ L^{\text{p}} = \dot{F}^{\text{p}} F^{p-1} = D^{\text{p}} = \dot{\bar{\epsilon}}^{\text{p}} M \quad (167) $$

with $M$ a symmetric tensor such that

$$ M \cdot M = \frac{3}{2} \quad \text{and} \quad \text{tr}[M] = 0 \quad (168) $$

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

$$ W(F, T, F^{\text{p}}, \bar{\epsilon}^{\text{p}}) = W^{\text{e}}(C^{\text{e}}, T) + W^{\text{p}}(F^{\text{p}}, \bar{\epsilon}^{\text{p}}, T) + W^{\text{t}}(T) \quad (169) $$

where $W^{\text{e}}$ is the elastically stored energy (recoverable), for which we have accounted for objectivity through the use of $C^{\text{e}}$, $W^{\text{p}}$ is the plastically stored
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

$$ \dot{W} = \mathbf{P} \cdot \dot{\mathbf{F}} - \underbrace{[(\mathbf{F}^e \mathbf{T} \mathbf{P} - \mathbf{Y}^c) \cdot (\mathbf{M} \mathbf{F}^p) - g]}_{\equiv Y} \frac{\dot{\mathbf{e}}^p}{\mathbf{e}^p} - \rho_0 \eta \dot{T} \quad (170) $$

where $\mathbf{P} = \partial_F W$ is the Piola stress tensor, $\mathbf{Y}^c = \partial_{F^p} W^p$ is a back-stress tensor (linked to kinematic hardening), and $g = \partial_{\bar{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 $\bar{e}^p$, when accounting for the specific flow rule chosen. This quantity can be rewritten as

$$ Y = (\mathbf{T} - \mathbf{T}^c) \cdot \mathbf{M} - g \quad (171) $$

where $\mathbf{T} = \mathbf{F}^e \mathbf{T} \mathbf{P} \mathbf{F}^p$ is Mandel stress tensor, and $\mathbf{T}^c = \mathbf{Y}^c \mathbf{F}^p$ is the associated backstress tensor.

### 5.2.2.1 Rate problem of finite visco-plasticity

Functional $D(\dot{\mathbf{F}}, \dot{\eta}, T, \dot{\bar{e}}^p, \mathbf{M})$ is then given by

$$ D(\dot{\mathbf{F}}, \dot{\eta}, T, \dot{\bar{e}}^p, \mathbf{M}; \mathbf{F}, \eta, \mathbf{F}^p, \bar{\mathbf{e}}^p) = \mathbf{P} \cdot \dot{\mathbf{F}} - Y \dot{\bar{e}}^p - \rho_0 \dot{\eta} (T - \mathbf{\Theta}) \\ + \Delta \left( \frac{T}{\mathbf{\Theta}} \dot{\bar{e}}^p; \bar{\mathbf{e}}^p, \mathbf{\Theta} \right) \quad (172) $$

and the effective rate potential is defined as

$$ D_{\text{eff}}(\dot{\mathbf{F}}, \dot{\eta}, T) = \inf_{\dot{\bar{e}}^p, \mathbf{M}} D(\dot{\mathbf{F}}, \dot{\eta}, T, \dot{\bar{e}}^p, \mathbf{M}; \mathbf{F}, \eta, \mathbf{F}^p, \bar{\mathbf{e}}^p). \quad (173) $$

Accounting for constraints on $\mathbf{M}$, the optimal plastic flow direction is given by

$$ \mathbf{M} = \sqrt{\frac{3}{2} \frac{\text{dev}[\mathbf{T} - \mathbf{T}^c]}{\|\text{dev}[\mathbf{T} - \mathbf{T}^c]\|}} = \sqrt{\frac{3}{2} \frac{\mathbf{s}}{\|\mathbf{s}\|}} \quad (174) $$

where $\mathbf{s} = \text{dev}[\mathbf{T} - \mathbf{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

$$ Y = \sqrt{\frac{3}{2} \mathbf{s} \cdot \mathbf{s}} - g = \sigma^{\text{eq}} - g \quad (175) $$
where $\sigma^{\text{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.

### 5.2.2.2 Exponential mapping

In the time-discrete setting, the local variational constitutive update problem takes the following form:

$$ W(F, T) = \inf_{\bar{\epsilon}^p, M} \left[ W(F, T, F^p, \bar{\epsilon}^p) + \delta t \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta \bar{\epsilon}^p}{\delta t} ; \bar{\epsilon}^p(\tau), T(\tau) \right) \right\rangle \right] \\ + \rho_0 \eta_0 \delta T - W(F_0, T_0, F_0^p, \bar{\epsilon}_0^p) \quad (176) $$

where one must provide a time-discrete update rule for plastic deformation $F^p$. The exponential update formula (Weber & Anand, 1990):

$$ F^p = \exp(\delta \bar{\epsilon}^p M) F_0^p \quad (177) $$

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 $\Delta$, in which case a consistent approximation to the average dissipation function over the time increment is given by:

$$ \left\langle \Delta \left( \frac{T}{T_0} \frac{\delta \bar{\epsilon}^p}{\delta t} ; \bar{\epsilon}^p(\tau) \right) \right\rangle \approx \Delta \left( \frac{T}{T_0} \frac{\delta \bar{\epsilon}^p}{\delta t} ; \bar{\epsilon}_\alpha^p \right). \quad (178) $$

### 5.2.2.3 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:

$$W^e(\mathbf{C}^e, T) = f(\operatorname{tr}[\boldsymbol{\varepsilon}^e], T) + \mu(T)(\operatorname{dev}[\boldsymbol{\varepsilon}^e])^2 \quad (179)$$

with

$$\boldsymbol{\varepsilon}^e = \log \mathbf{U}^e = \frac{1}{2} \log \mathbf{C}^e = \frac{1}{2} \log [\mathbf{F}^{p-T} \mathbf{C} \mathbf{F}^{p-1}]. \quad (180)$$

Considering the exponential update formula (177) and the ansatz that $\mathbf{C}_{\text{tr}}^e = \mathbf{F}_0^{p-T} \mathbf{C} \mathbf{F}_0^{p-1}$ and $\mathbf{M}$ are colinear, we can write

$$\boldsymbol{\varepsilon}^e = \frac{1}{2} \log \mathbf{C}_{\text{tr}}^e - \delta \bar{\mathbf{e}}^p \mathbf{M} = \boldsymbol{\varepsilon}_{\text{tr}}^e - \delta \bar{\mathbf{e}}^p \mathbf{M}. \quad (181)$$

Provided that there is no kinematic hardening ($\mathbf{Y}^c \equiv 0$), we thus formally recover the same variational problem as in linearized kinematics (with logarithmic elastic strains), which can also be interpreted as a radial return algorithm:

$$\mathbf{M} = \sqrt{\frac{3}{2} \frac{\operatorname{dev}[\boldsymbol{\varepsilon}_{\text{tr}}^e]}{\|\operatorname{dev}[\boldsymbol{\varepsilon}_{\text{tr}}^e]\|}} = \sqrt{\frac{3}{2} \frac{\operatorname{dev}[\boldsymbol{\varepsilon}^e]}{\|\operatorname{dev}[\boldsymbol{\varepsilon}^e]\|}} \quad (182)$$

and

$$-\sigma_{\text{tr}}^{\text{eq}} + 3\mu(T)\delta\bar{\mathbf{e}}^p + g(\bar{\mathbf{e}}^p, T) + \sigma_{\gamma}(\bar{\mathbf{e}}^p, T) = 0 \quad (183)$$

where

$$\sigma_{\text{tr}}^{\text{eq}} = \sqrt{6\mu(T)} \|\operatorname{dev}[\boldsymbol{\varepsilon}_{\text{tr}}^{e}]\| \quad (184)$$

while

$$\bar{\sigma}_{\gamma} = \frac{T}{T_0} \frac{\partial \phi}{\partial \dot{\bar{\mathbf{e}}}^p} \left( \frac{T}{T_0} \frac{\delta \bar{\mathbf{e}}^p}{\delta t} ; \bar{\mathbf{e}}_{\alpha}^p \right) + \alpha \delta t \frac{\partial \phi}{\partial \bar{\mathbf{e}}^p} \left( \frac{T}{T_0} \frac{\delta \bar{\mathbf{e}}^p}{\delta t} ; \bar{\mathbf{e}}_{\alpha}^p \right) \quad (185)$$

is a consistent incremental approximation of the dissipative part of the yield stress.

### 5.2.2.4 General isotropic hyperelasticity

In polycrystalline metals, elastic strains typically remain small, and the choice of the elastic free energy is mostly a matter of convenience (in the context of isotropic materials). Yet, visco-plasticity models can be useful in other contexts, such as thermoplastic polymers, or adhesive glues for example (Fancello, Goglio, Stainier, & Vassoler, 2008a). These materials may exhibit
larger elastic strains, and it may thus be required to use more general elastic free energy expressions to be able to reproduce their actual elastic behavior.

Considering isotropic materials, we can adopt the spectral approach described in Fancello et al. (2008b). Indeed, for such materials, thermodynamic potentials can be represented by functions of the principal values of their tensorial arguments, i.e.

$$W^e(\mathbf{C}^e, T) = W^e(c_1^e, c_2^e, c_3^e, T) \qquad (186)$$

where {$c_1^e, c_2^e, c_3^e$} are the principal values of tensor $\mathbf{C}^e$. In the same spirit, we can parameterize the plastic flow direction tensor by its principal values and vectors:

$$\mathbf{M} = \sum_{I=1}^{3} q_{I} \mathbf{N}_{I} \qquad (187)$$

where {$q_1, q_2, q_3$} are the principal values of $\mathbf{M}$, and $\mathbf{N}_I = \mathbf{n}_I \otimes \mathbf{n}_I$, with $\mathbf{n}_I$ the normalized principal directions of $\mathbf{M}$ ($\mathbf{n}_I \cdot \mathbf{n}_J = \delta_{IJ}$). Constraints on $\mathbf{M}$ translate as

$$\sum_{I=1}^{3} q_{I}^{2} = \frac{3}{2} \quad \text{and} \quad \sum_{I=1}^{3} q_{I} = 0. \qquad (188)$$

Keeping with the ansatz that $\mathbf{M}$ and $\mathbf{C}_{\text{tr}}^e$ are colinear (i.e. they share the same principal directions), we can then write

$$\varepsilon_I^e = \varepsilon^{\text{tr}} - \delta \bar{\varepsilon}^p q_I \quad (I = 1, 2, 3) \qquad (189)$$

where $\varepsilon_I^e = \log(\lambda_I^e) = \frac{1}{2}\log(c_I^e)$ are the principal logarithmic strains, while $\varepsilon^{\text{tr}} = \frac{1}{2}\log(c_I^{\text{tr}})$, with $c_I^{\text{tr}}$ the principal values of the predictor (trial) Cauchy-Green tensor. Restricting ourselves to the case without kinematic hardening ($\mathbf{Y}^c \equiv 0$), the variational constitutive update then becomes

$$\inf_{\bar{\mathbf{e}}^p, q_1, q_2, q_3} \left[ W^e(\varepsilon_1^e, \varepsilon_2^e, \varepsilon_3^e, T) + W^p(\bar{\mathbf{e}}^p, T) + \delta t \Delta \left( \frac{T}{T_0} \frac{\delta \bar{\mathbf{e}}^p}{\delta t} ; \bar{\mathbf{e}}_\alpha^p \right) \right] \qquad (190)$$

with relations (189) and constraints (188). In practice, we first determine the elastic or visco-plastic nature of the increment by looking at gradients of the objective function at singularity points corresponding to $\delta \bar{\mathbf{e}}^p = 0$. More precisely, as explained in details in Fancello et al. (2008b), we look for the minimal gradient in the $\bar{\mathbf{e}}^p$ direction, for all admissible values of $q_I$. It turns out
that this gradient can be expressed analytically and the condition for an elastic step takes the form:

$$ -\sigma_M^{\text{tr}} + g(\bar{\varepsilon}_0^p, T) + \frac{T}{T_0} \frac{\partial \Delta}{\partial \dot{\bar{\varepsilon}}^p}(0^+; \bar{\varepsilon}_0^p) > 0 \quad (191) $$

where

$$ \sigma_M(\varepsilon_1^e, \varepsilon_2^e, \varepsilon_3^e, T) = \sqrt{\frac{3}{2} \sum_{I=1}^{3} \left( \frac{\partial W^e}{\partial \varepsilon_I^e} - \frac{1}{3} \sum_{J=1}^{3} \frac{\partial W^e}{\partial \varepsilon_J^e} \right)^2} \quad (192) $$

and $\sigma_M^{\text{tr}} = \sigma_M(\varepsilon_1^{\text{tr}}, \varepsilon_2^{\text{tr}}, \varepsilon_3^{\text{tr}}, T)$. If the above gradient is negative, then there is an admissible descent direction, and the solution can be found by solving the four scalar stationarity equations with respect to $\bar{\varepsilon}^p$, $q_1$, $q_2$, $q_3$ [with constraints (188)].

It is also interesting to note that the equivalent stress $\sigma_M$ defined earlier also corresponds to a classical equivalent stress computed from Mandel stress tensor:

$$ \sigma_M = \sqrt{\frac{3}{2} \text{dev}[\boldsymbol{T}] \cdot \text{dev}[\boldsymbol{T}]} \quad (193) $$

The Mandel stress tensor is defined in the intermediate configuration and is symmetric in the case of isotropic elasticity. Thus, the suggested minimization strategy amounts to computing an elastic trial Mandel stress, evaluate the classical von Mises criterion in terms of this stress tensor, and, if the criterion is verified, apply a plastic corrector step. Once again, this process parallels the classical predictor-corrector scheme used in small strains. In general, it may not correspond to a radial return, though, since the predictor of $q_I$ (corresponding to the trial state) may be different from the overall optimal ones.

### 5.2.3 Alternative flow rules

The framework of variational constitutive update is of course not limited to von Mises plasticity. It has been, for example, directly extended to flow rules corresponding to Tresca plasticity (Ortiz & Stainier, 1999), Cam-Clay plasticity (Ortiz & Pandolfi, 2004), and porous plasticity (Weinberg, Mota, & Ortiz, 2006).

Mosler et al. have proposed a very generic approach allowing to reset a very wide class of classical visco-plasticity models into the variational
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).

### 5.2.4 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:

$$ \frac{\partial D_{\text{eff}}}{\partial \dot{\mathbf{F}}} = \mathbf{P} \qquad (194) $$

$$ \frac{\partial D_{\text{eff}}}{\partial T} = -\rho_0 \dot{\eta} + \frac{1}{\Theta} \left( \partial_{\dot{\bar{\epsilon}}^p} \Delta \left( \frac{T}{\Theta} \dot{\bar{\epsilon}}^p ; \bar{\epsilon}^p, \Theta \right) \right) \dot{\bar{\epsilon}}^p. \qquad (195) $$

If we take into account stationarity of $D_{\text{eff}}$ with respect to $\dot{\eta}$: $T = \Theta$ and stationarity of functional $D: Y = \partial_{\dot{\bar{\epsilon}}^p} \Delta(\dot{\bar{\epsilon}}^p, \bar{\epsilon}^p, T)$, we obtain that

$$ \frac{\partial D_{\text{eff}}}{\partial T} = -\rho_0 \dot{\eta} + \frac{1}{\Theta} Y \dot{\bar{\epsilon}}^p \qquad (196) $$

and comparing with Eq. (38), we obtain the following expression for the internal dissipation

$$ \mathcal{D}_{\text{int}} = Y \dot{\bar{\epsilon}}^p = (\mathbf{T} - \mathbf{T}^c) \cdot \mathbf{D}^p - g \dot{\bar{\epsilon}}^p \qquad (197) $$

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 $\mathbf{T} \cdot \mathbf{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:

$$ \mathcal{D}_{\text{int}} = \beta \mathbf{T} \cdot \mathbf{D}^p \qquad (198) $$
with $\beta$ considered as a constant material parameter (typically chosen as $\beta = 0.9$). This formula remains widely used today, despite more recent experimental work (see for example Hodowany, Ravichandran, Rosakis, & Rosakis, 2000; 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 Rosakis, Rosakis, Ravichandran, & Hodowany, 2000; Zehnder, 1991), aiming at providing an expression describing evolution of coefficient $\beta$. In the variational framework described here, the coefficient $\beta$ does not appear explicitly (although it can be computed a posteriori). Instead the ratio of dissipation to total visco-plastic 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 $\beta$, 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.

# 6. 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

$$ \mathcal{V} = \{ \boldsymbol{\varphi} : \mathcal{B}_0 \to \mathbb{R}^3 \mid \begin{array}{l} \boldsymbol{\varphi} = \bar{\boldsymbol{\varphi}}(t) \text{ on } \partial_u \mathcal{B}_0, \nabla_0 \boldsymbol{\varphi} \in \text{GL}_+(3, \mathbb{R}); \\ T : \mathcal{B}_0 \to \mathbb{R}_+ \mid T = \bar{T}(t) \text{ on } \partial_T \mathcal{B}_0 \end{array} \}. \quad (199) $$

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 $\mathcal{B}_0$:

$$ \mathcal{V}_h = \left\{ \begin{array}{l} \boldsymbol{\varphi}_h : \mathcal{B}_0 \to \mathbb{R}^3 \quad \left| \begin{array}{l} \boldsymbol{\varphi}_h(\mathbf{X}) = \displaystyle\sum_{a=1}^{N_{\text{nodes}}} N_a(\mathbf{X}) \mathbf{x}_a; \\ T_h : \mathcal{B}_0 \to \mathbb{R}_+ \quad \left| \begin{array}{l} T_h(\mathbf{X}) = \displaystyle\sum_{a=1}^{N_{\text{nodes}}} N_a(\mathbf{X}) T_a \end{array} \right. \end{array} \right. \right\} \quad (200) $$
where $\mathbf{x}_a = \boldsymbol{\varphi}(\mathbf{X}_a, t)$ are the positions of the $N_{\text{nodes}}$ nodes in the current (deformed) configuration and $T_a = T(\mathbf{X}_a, t)$ are the temperatures at these same nodes. As usual, the material shape functions $N_a(\mathbf{X})$ are defined on elements connected to node a and set to zero elsewhere. Index h refers to the mesh $\mathcal{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

$$F_h(\mathbf{X}) = \sum_{a=1}^{N_{\text{nodes}}} \mathbf{x}_a \otimes \nabla_0 N_a(\mathbf{X}) \qquad (201a)$$

$$G_h(\mathbf{X}) = - \frac{\sum_{a=1}^{N_{\text{nodes}}} T_a \nabla_0 N_a(\mathbf{X})}{\sum_{a=1}^{N_{\text{nodes}}} T_a N_a(\mathbf{X})} \qquad (201b)$$

and the variational principle (69) becomes

$$\inf_{\boldsymbol{\varphi}_h \in \mathcal{V}_h} \sup_{T_h \in \mathcal{V}_h} \mathcal{I}_h(\boldsymbol{\varphi}_h, T_h; \boldsymbol{\varphi}_0, T_0, \mathbf{Z}_0) \qquad (202)$$

where $\boldsymbol{\varphi}_0$ and $T_0$ are computed by interpolation from nodal values of $\mathbf{x}_a(t_0)$ and $T_a(t_0)$, while the internal variables $\mathbf{Z}_0$ are typically stored at integration points. Dirichlet boundary conditions are of course enforced by setting nodal variables $\mathbf{x}_a$ and $T_a$ to appropriate values on $\partial_u \mathcal{B}_0$ and $\partial_T \mathcal{B}_0$. Subscript h in $\mathcal{I}_h$ denotes the fact that volume integrals will generally be computed by numerical quadrature, based on the chosen discretization $\mathcal{T}_h$. Variations of incremental potential $\mathcal{I}_h$ are now taken with respect to nodal unknowns, and stationarity conditions are written

$$\sum_{a=1}^{N_{\text{nodes}}} f_a \cdot \delta \mathbf{x}_a = \sum_{a=1}^{N_{\text{nodes}}} \left\{ \int_{\mathcal{B}_0} [\mathbf{P}_h \cdot \nabla_0 N_a - \rho_0 \mathbf{b} N_a] dV - \int_{\partial_{\sigma} \mathcal{B}_0} \bar{t} N_a dS \right\} \cdot \delta \mathbf{x}_a \\ = 0 \quad \forall \delta \mathbf{x}_a \text{ adm.} \qquad (203a)$$

$$\sum_{a=1}^{N_{\text{nodes}}} Q_a \delta T_a = \sum_{a=1}^{N_{\text{nodes}}} \left\{ \int_{\mathcal{B}_0} \left[ -\rho_0 \delta \eta_h^{\text{eff}} N_a + \delta t \mathbf{H}_h \cdot \left( \frac{\nabla_0 N_a}{T_h} + \frac{\mathbf{G}_h}{T_h} N_a \right) + \delta t \rho_0 Q \frac{N_a}{T_h} \right] dV \right. \\ \left. - \int_{\partial_H \mathcal{B}_0} \delta t \bar{H} \frac{N_a}{T_h} dS \right\} \cdot \delta T_a = 0 \quad \forall \delta T_a \text{ adm.} \qquad (203b)$$
where $\mathbf{P}_h$, $\delta\eta_h^{\text{eff}}$ and $\mathbf{H}_h$ are, respectively, given by

$$ \mathbf{P}_h = \frac{\partial W}{\partial \mathbf{F}}(\mathbf{F}_h, T_h; \mathbf{F}_0, T_0, \mathbf{Z}_0) \quad (204) $$

$$ \rho_0 \delta \eta_h^{\text{eff}} = - \frac{\partial W}{\partial T}(\mathbf{F}_h, T_h; \mathbf{F}_0, T_0, \mathbf{Z}_0) \quad (205) $$

$$ \mathbf{H}_h = \frac{\partial}{\partial \mathbf{G}} \langle \chi(\mathbf{G}_h; \mathbf{F}_\alpha, T_\alpha, \mathbf{Z}_\alpha) \rangle. \quad (206) $$

The stationarity conditions (203a, b) yield discrete mechanical and thermal balance equations, which can alternatively be written

$$ \left\{ \begin{array}{l} \mathbf{f}_a^{\text{int}} - \mathbf{f}_a^{\text{ext}} = 0 \\ Q_a^{\text{int}} - Q_a^{\text{ext}} = 0 \end{array} \right. \quad (a=1, \dots, N_{\text{nodes}}) \quad (207) $$

where internal and external nodal forces and fluxes are given by

$$ \mathbf{f}_a^{\text{int}} = \bigcup_{e=1}^{N_{\text{lems}}} \int_{\Omega_0^e} \mathbf{P}_h \cdot \nabla_0 N_a \, dV \quad (208a) $$

$$ \mathbf{f}_a^{\text{ext}} = \bigcup_{e=1}^{N_{\text{lems}}} \int_{\Omega_0^e} \rho_0 \mathbf{b} N_a \, dV + \bigcup_{e=1}^{N_{\text{lems}}^{\sigma}} \int_{\partial \Omega_0^e} \bar{\mathbf{t}} N_a \, dS \quad (208b) $$

$$ Q_a^{\text{int}} = \bigcup_{e=1}^{N_{\text{lems}}} \int_{\Omega_0^e} \left[ \rho_0 \delta \eta_h^{\text{eff}} N_a - \delta t \mathbf{H}_h \cdot \left( \frac{\nabla_0 N_a}{T_h} + \frac{\mathbf{G}_h}{T_h} N_a \right) \right] dV \quad (208c) $$

$$ Q_a^{\text{ext}} = \bigcup_{e=1}^{N_{\text{lems}}} \int_{\Omega_0^e} \delta t \frac{\rho_0 Q}{T_h} N_a \, dV - \bigcup_{e=1}^{N_{\text{lems}}^H} \int_{\partial \Omega_0^e} \delta t \frac{\bar{H}}{T_h} N_a \, dS \quad (208d) $$

In the above expressions, $\cup$ denotes the element assembly operator, while $\Omega_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 $\delta\eta^{\text{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.
### 6.1.2 Mixed formulations

At the core of the variational formulation (69) lies the energy functional

$$I(\boldsymbol{\varphi}, T) = \int_{B_0} [\mathcal{W}(\mathbf{F}, T) - \delta t \chi(\mathbf{G})] \, dV \quad (209)$$

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

$$I(\boldsymbol{\varphi}, T) \approx \bigcup_{e=1}^{N_{\text{elements}}} I_h^e(\boldsymbol{\varphi}_h, T_h) \quad (210)$$

with elementary contributions given by

$$I_h^e(\boldsymbol{\varphi}_h, T_h) = \int_{\Omega_0^e} [\mathcal{W}(\mathbf{F}_h, T_h) - \delta t \chi(\mathbf{G}_h)] dV. \quad (211)$$

In the case of incompressible material behavior, or even quasi-incompressible 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):

$$\hat{I}_h^e(\boldsymbol{\varphi}_h, \theta_h, p_h, T_h) = \int_{\Omega_0^e} [\mathcal{W}(\hat{\mathbf{F}}_h^e, T_h) + p_h^e(J_h - \theta_h^e) - \delta t \chi(\mathbf{G}_h)] dV \quad (212)$$

where $\theta_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 $\Omega_0^e$). The modified deformation gradient $\hat{\mathbf{F}}_h^e$ is defined as

$$\hat{\mathbf{F}}_h^e = \left(\frac{\theta_h^e}{J_h}\right)^{\frac{1}{3}} \nabla_0 \boldsymbol{\varphi}_h \quad \text{with} \quad J_h = \det[\nabla_0 \boldsymbol{\varphi}_h]. \quad (213)$$

The variational principle then becomes

$$\inf_{\boldsymbol{\varphi}_h, \theta_h} \sup_{p_h, T_h} \hat{\mathcal{I}}(\boldsymbol{\varphi}_h, \theta_h, p_h, T_h) \quad \text{with} \quad \boldsymbol{\varphi}_h, \theta_h, p_h \text{ and } T_h \text{ admissible} \quad (214)$$
with

$$
\begin{equation}
\begin{split}
\hat{\mathcal{I}}(\boldsymbol{\varphi}_h, \theta_h, p_h, T_h) ={}& \hat{\mathcal{I}}(\boldsymbol{\varphi}_h, \theta_h, p_h, T_h) - \int_{\mathcal{B}_0} \rho_0 \boldsymbol{b}(t) \cdot \delta\boldsymbol{\varphi}_h dV \\
& - \int_{\partial_{\sigma}\mathcal{B}_0} \bar{\boldsymbol{t}}(t) \cdot \delta\boldsymbol{\varphi}_h dS + \int_{\mathcal{B}_0} \delta t \rho_0 Q(t) \log \frac{T_h}{T_0} dV \\
& - \int_{\partial_H \mathcal{B}_0} \delta t \bar{H}(t) \log \frac{T_h}{T_0} dS.
\end{split}
\tag{215}
\end{equation}
$$

Euler–Lagrange equations pertaining to $\theta_h$ and $p_h$ are given by, respectively:

$$
\begin{equation}
\begin{split}
\langle D_{\theta} \hat{\mathcal{I}}(\boldsymbol{\varphi}_h, \theta_h, p_h, T_h), \tilde{\theta}_h \rangle &= \bigcup_{e=1}^{N_{\text{elem}}} \int_{\Omega_0^e} \left[ \langle D_{\theta} W(\hat{\mathbf{F}}_h^e, T_h), \tilde{\theta}_h^e \rangle - p_h^e \tilde{\theta}_h^e \right] dV = 0 \\
&\qquad \forall \tilde{\theta}_h^e \quad \text{adm.}
\end{split}
\tag{216}
\end{equation}
$$

$$
\langle D_p \hat{\mathcal{I}}(\boldsymbol{\varphi}_h, \theta_h, p_h, T_h), \tilde{p}_h \rangle = \bigcup_{e=1}^{N_{\text{elem}}} \left[ \int_{\Omega_0^e} (J_h - \theta_h^e) dV \right] \tilde{p}_h^e = 0 \quad \forall \tilde{p}_h^e \text{ adm.} \tag{217}
$$

Equation (217) immediately yields the expression for the piecewise constant
volumic deformation

$$
\theta_h^e = \frac{1}{|\Omega_0^e|} \int_{\Omega_0^e} \det[\nabla_0 \boldsymbol{\varphi}_h] \, dV. \quad (218)
$$

Considering that

$$
\langle D_{\theta} W(\hat{\mathbf{F}}_{h}^{e}, T_{h}), \tilde{\theta}_{h}^{e} \rangle = \frac{\partial W}{\partial \mathbf{F}}(\hat{\mathbf{F}}_{h}^{e}, T_{h}) \cdot \left(\frac{1}{3} \theta_{h}^{e-1} \tilde{\theta}_{h}^{e} \hat{\mathbf{F}}_{h}^{e}\right) = \frac{1}{3} \theta_{h}^{e-1} \text{tr}[\boldsymbol{\tau}(\hat{\mathbf{F}}_{h}^{e}, T_{h})] \tilde{\theta}_{h}^{e} \tag{219}
$$

where we have introduced the notation:

$$
\boldsymbol{\tau}(\hat{\mathbf{F}}_h^e, T_h) \equiv \frac{\partial W}{\partial \mathbf{F}}(\hat{\mathbf{F}}_h^e, T_h) \hat{\mathbf{F}}_h^{eT}, \quad (220)
$$

Eq. (216) then yields

$$
p_h^e = \frac{1}{|\Omega_0^e|} \int_{\Omega_0^e} p(\hat{\mathbf{F}}_h^e, T_h) dV \quad (221)
$$
with

$$p(\hat{F}_h^e, T_h) = \frac{1}{3} \theta_h^{e-1} \mathrm{tr}[\boldsymbol{\tau}(\hat{F}_h^e, T_h)]. \qquad (222)$$

We can now proceed to compute elementary contributions to nodal forces and fluxes. Since we have

$$\langle D_{\varphi}, J_h, \tilde{\boldsymbol{\varphi}}_h \rangle = J_h \nabla_0 \tilde{\boldsymbol{\varphi}}_h \cdot \nabla_0 \boldsymbol{\varphi}_h^T = J_h \mathrm{tr}[\nabla \tilde{\boldsymbol{\varphi}}_h] \qquad (223)$$

we can write

$$\begin{aligned}
\langle D_{\varphi} W(\hat{F}_h^e, T_h), \tilde{\boldsymbol{\varphi}}_h \rangle &= \frac{\partial W}{\partial \boldsymbol{F}}(\hat{F}_h^e, T_h) \cdot \theta_h^{e\frac{1}{3}} \left( J_h^{-\frac{1}{3}} \nabla_0 \tilde{\boldsymbol{\varphi}}_h - \frac{1}{3} J_h^{-\frac{4}{3}} \langle D_{\varphi} J_h, \tilde{\boldsymbol{\varphi}}_h \rangle \nabla_0 \boldsymbol{\varphi}_h \right) \\
&= \frac{\partial W}{\partial \boldsymbol{F}}(\hat{F}_h^e, T_h) \cdot \theta_h^{e\frac{1}{3}} J_h^{-\frac{1}{3}} \left( \nabla_0 \tilde{\boldsymbol{\varphi}}_h - \frac{1}{3} \mathrm{tr}[\nabla \tilde{\boldsymbol{\varphi}}_h] \nabla_0 \boldsymbol{\varphi}_h \right) \\
&= \tau(\hat{F}_h^e, T_h) \cdot \mathrm{dev}[\nabla \tilde{\boldsymbol{\varphi}}_h].
\end{aligned} \qquad (224)$$

We then have

$$\langle D_{\varphi} I_h^e(\boldsymbol{\varphi}_h, T_h), \tilde{\boldsymbol{\varphi}}_h \rangle = \int_{\Omega_0^e} \left( \boldsymbol{\tau}(\boldsymbol{F}_h^e, T_h) \cdot \mathrm{dev}[\nabla \tilde{\boldsymbol{\varphi}}_h] + p_h^e J_h \mathrm{tr}[\nabla \tilde{\boldsymbol{\varphi}}_h] \right) dV \qquad (225)$$

which can alternatively be written

$$\langle D_{\varphi} I_h^e(\boldsymbol{\varphi}_h, T_h), \tilde{\boldsymbol{\varphi}}_h \rangle = \int_{\Omega_0^e} \left( \mathrm{dev}[\boldsymbol{\tau}(\boldsymbol{F}_h^e, T_h)] + J_h p_h^e \boldsymbol{I} \right) \cdot \nabla \tilde{\boldsymbol{\varphi}}_h dV. \qquad (226)$$

The latter expression defines a recovery method for Kirchhoff stress field:

$$\boldsymbol{\tau}_h^e \equiv \mathrm{dev}[\boldsymbol{\tau}(\boldsymbol{F}_h^e, T_h)] + J_h p_h^e \boldsymbol{I} \qquad (227)$$

or, alternatively, an effective Piola–Kirchhoff stress tensor:

$$\boldsymbol{P}_h^e \equiv \left(\frac{\theta_h^e}{J_h}\right)^{\frac{1}{3}} \left[ \frac{\partial W}{\partial \boldsymbol{F}}(\hat{F}_h^e, T_h) - \frac{1}{3} \left(\frac{\partial W}{\partial \boldsymbol{F}} \cdot \hat{F}_h^e\right) \hat{F}_h^{e-T} + J_h p_h^e \hat{F}_h^{e-T} \right]. \qquad (228)$$

The nodal internal forces array is then given by
$$f_a^{\text{int}} = \bigcup_{e=1}^{N_{\text{lems}}} \int_{\Omega_0^e} \mathbf{P}_h^e \cdot \nabla_0 N_a \, dV. \qquad (229)$$

Nodal internal fluxes are identical to those obtained through the standard formulation (208c), where the net entropy variation is computed by

$$\rho_0 \delta \eta_h^{\text{eff}} = \frac{\partial W}{\partial T}(F_h^e, T_h). \qquad (230)$$

## 6.2. 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).
Regarding the modeling of discontinuities and/or interfaces, an alternative approach consists in directly formulating interface models, for example with plasticity (Jérusalem, Stainier, & Radovitzky, 2007) or cohesive damage (Mosler & Scheider, 2011).

# 7. 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.

## 7.1. 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:

$$ \bar{W}(\bar{F}) = \inf_{\varphi \in K(\bar{F})} \frac{1}{|\mathcal{B}_0|} \int_{\mathcal{B}_0} \bar{W}(\nabla_0 \varphi) dV \quad (231) $$

where $\mathcal{B}_0$ is the domain occupied by a RVE, and $K(\bar{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 $\bar{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-visco-plastic phases (small strains).

Note that all these contributions have, for the moment, been restricted to the isothermal case.

## 7.2. 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, 2009).

## 8. 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 thermo-mechanics and diffusion-mechanics problems (or thermo-diffusion-mechanics 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).

## ACKNOWLEDGEMENTS

Part of this work was funded by the Région des Pays de la Loire (under project VAMPS).

## REFERENCES

* Armero, F., & Simo, J. C. (1992). A new unconditionally stable fractional step method for nonlinear coupled thermomechanical problems. *International Journal for Numerical Methods in Engineering, 35*(4), 737–766.

* Armero, F., & Simo, J. C. (1993). A priori stability estimates and unconditionally stable product formula algorithms for nonlinear coupled thermoplasticity. *International Journal of Plasticity, 9*, 749–782.
Aubry, S., & Ortiz, M. (2003). The mechanics of deformation-induced subgrain dislocation structures in metallic crystals at large strains. In *Proceedings of the Royal Society of London A*, 459, 3131–3158.

Balzani, D., & Ortiz, M. (2012). Relaxed incremental variational formulation for damage at large strains with application to fiber-reinforced materials and materials with truss-like microstructures. *International Journal for Numerical Methods in Engineering, 92*, 551–570.

Batra, G. (1989). On a principle of virtual work for thermo-elastic bodies. *Journal of Elasticity*, 21, 131–146.

Ben-Amoz, M. (1965). On a variational theorem in coupled thermoelasticity. *Journal of Applied Mechanics*, 32, 943–945.

Bertsekas, D. P. (1995). *Nonlinear programming*. Belmont, MA: Athena Scientific.

Biot, M. (1958). Linear thermodynamics and the mechanics of solids. In *Proceedings of the Third US National Congress of Applied Mechanics* (pp. 1–18). Providence, RI: ASME.

Biot, M. A. (1956). Thermoelasticity and irreversible thermodynamics. *Journal of Statistical Physics*, 27, 250–253.

Bleier, N., & Mosler, J. (2012). Efficient variational constitutive updates by means of a novel parameterization of the flow rule. *International Journal for Numerical Methods in Engineering, 89*(9), 1120–1143.

Bleier, N., & Mosler, J. (2013). A hybrid variationally consistent homogenization approach based on Ritz's method. *International Journal for Numerical Methods in Engineering, 94*(7), 625–647.

Bourdin, B., Francfort, G. A., & Marigo, J. J. (2008). *The variational approach to fracture*. Berlin: Springer.

Brassart, L., & Stainier, L. (2012). On convergence properties of variational constitutive updates for elasto-visco-plasticity. *GAMM-Mitteilungen*, 35(1), 26–42.

Brassart, L., Stainier, L., Doghri, I., & Delannay, L. (2011). A variational formulation for the incremental homogenization of elasto-plastic composites. *Journal of the Mechanics and Physics of Solids*, 59(12), 2455–2475.

Brassart, L., Stainier, L., Doghri, I., & Delannay, L. (2012). Homogenization of elasto-(visco) plastic composites based on an incremental variational principle. *International Journal of Plasticity*, 36, 86–112.

Canadija, M., & Mosler, J. (2011). On the thermomechanical coupling in finite strain plasticity theory with non-linear kinematic hardening by means of incremental energy minimization. *International Journal of Solids and Structures*, 48(7–8), 1120–1129.

Carini, A. (1996). Colonnetti's minimum principle extension to generally non-linear materials. *International Journal of Solids and Structures*, 33(1), 121–144.

Comi, C., Corigliano, A., & Maier, G. (1991). Extremum properties of finite-step solutions in elasto-plasticity with nonlinear hardening. *International Journal of Solids and Structures*, 27(8), 965–981.

Conti, S., & Ortiz, M. (2005). Dislocation microstructures and the effective behavior of single crystals. *Archive for Rational Mechanics and Analysis*, 176(1), 103–147.

Conti, S., & Ortiz, M. (2008). Minimum principles for the trajectories of systems governed by rate problems. *Journal of the Mechanics and Physics of Solids*, 56, 1885–1904.

Fancello, E., Goglio, L., Stainier, L., & Vassoler, J. M. (2008a). Identification of the strain rate parameters for structural adhesives. *Journal of Adhesion Science and Technology*, 22(13), 1523–1540.

Fancello, E., Ponthot, J.-P., & Stainier, L. (2006). A variational formulation of constitutive models and updates in non-linear finite viscoelasticity. *International Journal for Numerical Methods in Engineering, 65*(11), 1831–1864.

Fancello, E., Vassoler, J., & Stainier, L. (2008b). A variational constitutive update algorithm for a set of isotropic hyperelastic-viscoplastic models. *Computer Methods in Applied Mechanics and Engineering*, 197(49–50), 4132–4148.
Francfort, G. A., & Marigo, J. J. (1998). Revisiting brittle fracture as an energy minimization problem. *Journal of the Mechanics and Physics of Solids*, **46**, 1319–1342.

François-Lavet, V., Henrotte, F., Stainier, L., Noels, L., & Geuzaine, C. (2013). An energy-based variational model of ferromagnetic hysteresis for finite element computations. *Journal of Computational and Applied Mathematics*, **246**, 243–250.

Hackl, K. (1997). Generalized standard media and variational principles in classical and finite strain elastoplasticity. *Journal of the Mechanics and Physics of Solids*, **45**(5), 667–688.

Halphen, B., & Nguyen, Q. S. (1975). Sur les matériaux standard généralisés. *Journal de Mécanique*, **14**(1), 39–63.

Han, W., Jensen, S., & Reddy, B. D. (1997). Numerical approximations of problems in plasticity: Error analysis and solution algorithms. *Numerical Linear Algebra with Applications*, **4**, 191–204.

Herrmann, G. (1963). On variational principles in thermoelasticity and heat conduction. *Quarterly of Applied Mathematics*, **22**, 151–155.

Hodowany, J., Ravichandran, G., Rosakis, A., & Rosakis, P. (2000). Partition of plastic work into heat and stored energy in metals. *Experimental Mechanics*, **40**(2), 113–123.

Holzapfel, G. A. (2000). *Nonlinear solid mechanics. A continuum approach for engineering*. Chichester: Wiley.

Hughes, T. J. R. (2000). *The finite element method: Linear static and dynamic finite element analysis*. Mineola, NY: Dover.

Jérusalem, A., Stainier, L., & Radovitzky, R. (2007). A continuum model describing the reverse grain-size dependence of the strength of nanocrystalline metals. *Philosophical Magazine*, **87**(17), 2541–2559.

Johnson, G., Ortiz, M., & Leyendecker, S. (2012). A linear programming-based algorithm for the signed separation of (non-smooth) convex bodies. *Computer Methods in Applied Mechanics and Engineering*, **233–236**, 49–67.

Kane, C., Repetto, E. A., Ortiz, M., & Marsden, J. E. (1999). Finite element analysis of nonsmooth contact. *Computer Methods in Applied Mechanics and Engineering*, **180**(1–2), 1–26.

Kintzel, O., & Mosler, J. (2010). A coupled isotropic elasto-plastic damage model based on incremental minimization principles. *Technische Mechanik*, **30**(1–3), 177–184.

Kintzel, O., & Mosler, J. (2011). An incremental minimization principle suitable for the analysis of low cycle fatigue in metals: A coupled ductile-brittle damage model. *Computer Methods in Applied Mechanics and Engineering*, **200**(45–46), 3127–3138.

Lanczos, C. (1986). *The variational principles of mechanics*. Mineola, NY: Dover.

Lee, E. (1969). Elastic-plastic deformation at finite strains. *Journal of Applied Mechanics*, **36**, 1–6.

Leroy, Y. M., & Molinari, A. (1992). Stability of steady states in shear zones. *Journal of the Mechanics and Physics of Solids*, **40**(I), 181–212.

Lippmann, H. (1978). Extremum and variational principles in mechanics. In *Proceedings of the CISM* (2nd ed.). Austria: Springer-Verlag.

Lubliner, J. (1990). *Plasticity theory*. New York: Macmillan.

Martin, J., Kaunda, M., & Isted, R. (1996). Internal variable formulations of elastic-plastic dynamic problems. *International Journal of Impact Engineering*, **18**(7–8), 849–858.

Mialon, P. (1986). Eléments d'analyse et de résolution numérique des relations de l'élasto-plasticité. *EDF Bulletin de la Direction des Etudes et Recherches Série C*, **3**, 57–89.

Miehe, C. (1996). Exponential map algorithm for stress updates in anisotropic multiplicative elasto-plasticity for single crystals. *International Journal for Numerical Methods in Engineering*, **39**, 3367–3390.

Miehe, C. (2002). Strain-driven homogenization of inelastic microstructures and composites based on an incremental variational formulation. *International Journal for Numerical Methods in Engineering*, **55**(11), 1285–1322.
Miehe, C., Rosato, D., & Kiefer, B. (2011). Variational principles in dissipative electro-magneto-mechanics: A framework for the macro-modeling of functional materials. *International Journal for Numerical Methods in Engineering, 86*(10), 1225-1276.

Mielke, A., & Ortiz, M. (2008). A class of minimum principles for characterizing the trajectories and the relaxation of dissipative systems. *ESAIM: Control, Optimisation and Calculus of Variations, 14*, 494-516.

Molinari, A., & Ortiz, M. (1987). Global viscoelastic behavior of heterogeneous thermoelastic materials. *International Journal of Solids and Structures, 23*(9), 1285-1300.

Mosler, J. (2010). Variationally consistent modeling of finite strain plasticity theory with non-linear kinematic hardening. *Computer Methods in Applied Mechanics and Engineering, 199*(45-48), 2753-2764.

Mosler, J., & Bruhns, O. (2009). Towards variational constitutive updates for non-associative plasticity models at finite strain: Models based on a volumetric-deviatoric split. *International Journal of Solids and Structures, 46*(7-8), 1676-1684.

Mosler, J., & Ortiz, M. (2006). On the numerical implementation of variational arbitrary Lagrangian-Eulerian (VALE) formulations. *International Journal for Numerical Methods in Engineering, 67*, 1272-1289.

Mosler, J., & Ortiz, M. (2007). Variational h-adaption in finite deformation elasticity and plasticity. *International Journal for Numerical Methods in Engineering, 72*, 505-523.

Mosler, J., & Ortiz, M. (2009). An error-estimate-free and remapping-free variational mesh refinement and coarsening method for dissipative solids at finite strains. *International Journal for Numerical Methods in Engineering, 77*(3), 437-450.

Mosler, J., & Scheider, I. (2011). A thermodynamically and variationally consistent class of damage-type cohesive models. *Journal of the Mechanics and Physics of Solids, 59*(8), 1647-1668.

Noels, L., Stainier, L., & Ponthot, J. P. (2006). An energy momentum conserving algorithm using the variational formulation of visco-plastic updates. *International Journal for Numerical Methods in Engineering, 65*(6), 904-942.

Noels, L., Stainier, L., & Ponthot, J. P. (2008). A first-order energy-dissipative momentum-conserving scheme for elasto-plasticity using the variational updates formulation. *Computer Methods in Applied Mechanics and Engineering, 197*(6-8), 706-726.

Oden, T. J., & Reddy, J. N. (1976). *Variational methods in theoretical mechanics*. Berlin: Springer.

Ortiz, M., & Pandolfi, A. (2004). A variational cam-clay theory of plasticity. *Computer Methods in Applied Mechanics and Engineering, 193*, 2645-2666.

Ortiz, M., Repetto, E. A., & Stainier, L. (2000). A theory of subgrain dislocation structures. *Journal of the Mechanics and Physics of Solids, 48*(10), 2077-2114.

Ortiz, M., & Stainier, L. (1999). The variational formulation of viscoplastic constitutive updates. *Computer Methods in Applied Mechanics and Engineering, 171*(3-4), 419-444.

Pandolfi, A., Kane, C., Marsden, J. E., & Ortiz, M. (2002). Time-discretized variational formulation of non-smooth frictional contact. *International Journal for Numerical Methods in Engineering, 53*(8), 1801-1829.

Radovitzky, R., & Ortiz, M. (1999). Error estimation and adaptive meshing in strongly nonlinear dynamic problems. *Computer Methods in Applied Mechanics and Engineering, 172*(1-4), 203-240.

Rice, J. (1971). Inelastic constitutive relations for solids: an internal-variable theory and its application to metal plasticity. *Journal of the Mechanics and Physics of Solids, 19*(6), 443-455.

Rittel, D., Bhattacharyya, A., Poon, B., Zhao, J., & Ravichandran, G. (2007). Thermo-mechanical characterization of pure polycrystalline tantalum. *Materials Science and Engineering A, 447*(1-2), 65-70.

Rosakis, P., Rosakis, A., Ravichandran, G., & Hodowany, J. (2000). A thermodynamic internal variable model for the partition of plastic work into heat and stored energy in metals. *Journal of the Mechanics and Physics of Solids, 48*(3), 581-607.
Sidoroff, F. (1974). *Un modèle viscoélastique non linéaire avec configuration intermédiaire.* *Journal de Mécanique*, **13**, 679-713.

Simo, J., & Hughes, T. (1998). *Computational Inelasticity*. No. 7 in *Interdisciplinary applied mathematics*. New York: Springer.

Simo, J., & Miehe, C. (1992). Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation. *Computer Methods in Applied Mechanics and Engineering*, **98**, 41–104.

Simo, J., & Taylor, R. (1991). Quasi-incompressible finite elasticity in principal stretches—continuum basis and numerical algorithms. *Computer Methods in Applied Mechanics and Engineering*, **85**(3), 273–310.

Stainier, L. (2011). Consistent incremental approximation of dissipation pseudo-potentials in the variational formulation of thermo-mechanical constitutive updates. *Mechanics Research Communications*, **38**(4), 315–319.

Stainier, L., Cuitiño, A. M., & Ortiz, M. (2002). A micromechanical model of hardening, rate sensitivity and thermal softening in bcc single crystals. *Journal of the Mechanics and Physics of Solids*, **50**(7), 1511–1545.

Stainier, L., & Ortiz, M. (2010). Study and validation of a variational theory of thermomechanical coupling in finite viscoplasticity. *International Journal of Solids and Structures*, **47**(5), 705–715.

Su, S. (2012). Energy-based variational variational modeling of adiabatic shear band structure. Ph.D. thesis, Ecole Centrale de Nantes, http://tel.archives-ouvertes.fr/tel-00797178.

Taylor, G., & Quinney, H. (1937). The latent heat remaining in a metal after cold working. In *Proceedings of the Royal Society of London A*, **163**, 157–181.

Thomas, J., & Triantafyllidis, N. (2009). On electromagnetic forming processes in finitely strained solids: Theory and examples. *Journal of the Mechanics and Physics of Solids*, **57**(8), 1391–1416.

Titchener, A., & Bever, M. (1958). The stored energy of cold work. *Progress in Metal Physics*, **7**, 247–338.

Weber, G., & Anand, L. (1990). Finite deformation constitutive-equations and a time integration procedure for isotropic, hyperelastic viscoplastic solids. *Computer Methods in Applied Mechanics and Engineering*, **79**(2), 173–202.

Weinberg, K., Mota, A., & Ortiz, M. (2006). A variational constitutive model for porous metal plasticity. *Computational Mechanics*, **37**, 142–152.

Wilkins, M. (1964). Calculation of elastic-plastic flow. *Methods of computational physics*: Vol. 3. New York: Academic Press.

Yang, Q., Mota, A., & Ortiz, M. (2005). A class of variational strain-localization finite elements. *International Journal for Numerical Methods in Engineering*, **62**, 1013–1037.

Yang, Q., Stainier, L., & Ortiz, M. (2006). A variational formulation of the coupled thermomechanical boundary-value problem for general dissipative solids. *Journal of the Mechanics and Physics of Solids*, **54**(2), 401–424.

Zehnder, A. (1991). A model for the heating due to plastic work. *Mechanics Research Communications*, **18**(1), 23–28.

Ziegler, H. (1977). *An introduction to thermomechanics*. Amsterdam: Elsevier/North-Holland.