A Variational Framework for the Thermomechanics of Gradient-Extended Dissipative Solids -- with Applications to Diffusion, Damage and Plasticity

The paper presents a versatile framework for solids which undergo nonisothermal processes with irreversibly changing microstructure at large strains. It outlines rate-type and incremental variational principles for the full thermomechanical coupling in gradient-extended dissipative materials. It is shown that these principles yield as Euler equations essentially the macro- and micro-balances as well as the energy equation. Starting point is the incorporation of the entropy and entropy rate as canonical arguments into constitutive energy and dissipation functions, which additionally depend on the gradient-extended mechanical state and its rate, respectively. By means of (generalized) Legendre transformations, extended variational principles with thermal as well as mechanical driving forces can be constructed. On the thermal side, a rigorous distinction between the quantity conjugate to the entropy and the quantity conjugate to the entropy rate is essential here. Formulations with mechanical driving forces are especially suitable when considering possibly temperature-dependent threshold mechanisms. With regard to variationally consistent incrementations, we suggest an update scheme which renders the exact form of the intrinsic dissipation and is highly suitable when considering adiabatic processes. It is shown that this proposed numerical algorithm has the structure of an operator split. To underline the broad applicability of the proposed framework, we set up three model problems as applications: Cahn-Hilliard diffusion coupled with temperature evolution, where we propose a new variational principle in terms of the species flux vector, as well as thermomechanics of gradient damage and gradient plasticity. In a numerical example we study the formation of a cross shear band.


Introduction
In order to model dissipative size effects in solid materials that are, for example, related to the width of shear bands or the grain size in polycrystals, nonstandard continuum theories have to be elaborated that are based on characteristic length scale parameters. The idea is to incorporate at least first-order spatial gradients of micro-structural variables that describe (possibly in a homogenized sense) the irreversibly evolving microstructure. It is of great importance here to account for thermomechanical coupling effects such as thermal softening and temperature increase due to intrinsic dissipation. Practical applications include technological processes like sheet-metal forming or extrusion of metals. Within this work, we present a unified framework for the fully coupled thermomechanics of gradient-extended dissipative solids that is applicable to a wide range of model problems such as diffusion, (crystal) plasticity and damage. The formulation is embedded into the concept of standard dissipative solids that is characterized by energy and dissipation functions, see Biot (1965), Ziegler (1963) and Halphen & Nguyen (1975). In particular, the strongly coupled multifield problem will exhibit an incremental variational structure which is an extension of the framework of gradient-extended dissipative solids presented in Miehe (2011Miehe ( , 2014 towards nonisothermal processes. Variational Thermomechanics of Gradient-Extended Continua 2 Historically, the consideration of long-range effects via length scales can be traced back to the work of the brothers Cosserat & Cosserat (1896). They investigated a microstructure with rigid particles by introducing an additional microrotation governed by three additional degrees of freedom. Since then, a lot of research has been done on so called micropolar, micromorphic and microstrain continua, i.e., by Eringen (2012), see also Leismann & Mahnken (2015) for a comparative study. The book of Capriz (2013) provides a general setting for materials with microstructure where order parameters are considered as components of elements of an abstract manifold. Here, we also refer to the work of Mariano (2002). In Maugin (1990) the standard thermodynamic theory of local internal variables is extended by an additional dependency of the energy function on the first-order spatial gradient of a (not necessarily scalar) internal variable, see also Maugin & Muschik (1994) and Frémond (2013). An important ingredient of such theories is a proper energetic treatment of the microstructural processes yielding additional balance-type equations that are coupled with the standard macroscopic balance equations (mass, linear/angular momentum and energy). For an embedding into the theory of micro-forces we refer to Gurtin (2003).
Typical applications of above mentioned nonlocal theories including additional balance equations are theories of phase transformation, gradient plasticity and gradient damage. Theories of gradient plasticity are dealt with on a phenomenological level in Aifantis (1987), Fleck & Hutchinson (2001), Gurtin (2003), Forest & Sievert (2003), Gurtin & Anand (2005) just to name a few. Gradient damage theories that are basically in the same spirit are considered in Peerlings et al. (1996), Dimitrijevic & Hackl (2008), Pham et al. (2011), Frémond (2013 among many others. An important role in the modeling of dissipative processes in solids is played by incremental variational principles. They offer an elegant way to couple different physical fields, possess inherent symmetry and require only Jacobian and Hessian matrices of the underlying potential density for an implementation into typical finite element codes. In addition, if a minimization structure is at hand, structural and material stability analysis is possible, i.e., the formation of complex microstructures can be described. For contributions on the variational theory of local plasticity we refer to Hill (1950), Simo & Honein (1990), Hackl (1997), Ortiz & Stainier (1999), Ortiz & Repetto (1999), Carstensen et al. (2002) and Mosler & Bruhns (2010). The works of Miehe (2002), Petryk (2003) and Hackl & Fischer (2007) deal with general inelastic behavior. Extensions towards incorporation of (at least) the first gradient of the internal variable into constitutive functions can be found in Mühlhaus & Aifantis (1991), Fleck & Willis (2009), Mainik & Mielke (2009), Nguyen (2011, Miehe (2014) and Lancioni et al. (2015) for gradient plasticity and in Lorentz & Andrieux (1999), Mielke & Roubíček (2006), Dimitrijevic & Hackl (2008) and Pham et al. (2011) for gradient damage. In addition, we mention the general formulations of gradient-extended dissipative solids by Svendsen (2004), Francfort & Mielke (2006) and Miehe (2011). A gradient plasticity theory that incorporates a fractional derivative of the plastic strain is presented in Dahlberg & Ortiz (2019). Alternatively, in damage mechanics one can account for nonlocality by considering the gradient of the total strain measure (instead of the damage variable) in the constitutive functions. Non-smooth (hemi)variational formulations of such an approach are presented in, e.g., , Placidi et al. (2019) and Timofeev et al. (2021). Studies comparing these formulations with the gradient damage approach can be found in Le et al. (2018) and . Note that all works cited so far in this paragraph deal with dissipative processes under isothermal conditions except Hackl (1997) in the last section. In the latter work, an original idea of Simo & Miehe (1992) is taken up, namely to introduce a plastic configurational entropy as an additional internal variable that, as part of the total entropy, leaves the internal energy unchanged. From the viewpoint of variational principles however, the main difficulty of fully coupled thermomechanics is to account for the internal dissipation in the energy equation. The construction of a single unified potential from which all coupled thermomechanical field equations follow is outlined in the seminal work of Yang et al. (2006). They introduce a specific integrating factor that makes the linearized integral expressions, showing up in the fully coupled weak form, symmetric. This factor is based on distinguishing the equilibrium temperature from a so-called external temperature. Both temperatures coincide only at equilibrium. With this variational principle at hand, the fully coupled thermomechanical boundary value problem can be solved monolithically by a sequence of symmetric algebraic systems. Alternatively, staggered numerical algorithms were mainly used before which lead to symmetric finite element stiffness matrices for each subproblem, see Simo & Miehe (1992) and . In the recent years, the work of Yang et al. (2006) has been the basis for further model developments in thermoplasticity, see Stainier & Ortiz (2010), Canadija & Mosler (2011), Su et al. (2014), Bartels et al. (2015), Mielke (2016) and Fohrmeister et al. (2018). In the latter, a specific phenomenological model of gradient plasticity together with a micromorphic extension in the sense of Forest (2009) is considered. Besides incorporating the first gradient of a generic internal variable into the energy function, Mielke (2016) also discusses in detail the gradient structure of thermoplasticity. In Bartels et al. (2015) special focus is put on the thermodynamic (in)consistency of the Taylor-Quinney factor. This factor, usually introduced in an ad-hoc way, takes into account that only a certain amount of the plastic power is transformed into heat as observed in the experiments of Taylor & Quinney (1934), see also Rosakis et al. (2000) and Oliferuk & Raniecki (2018).
According to the authors' knowledge, a universal variational framework for the thermomechanics of solids with energetic as well as dissipative gradient extensions is still missing in literature. The present work outlines a generalization of the versatile framework of Miehe (2011Miehe ( , 2014 for gradient-extended dissipative solids towards nonisothermal processes. Point of departure is the definition of two constitutive functions, namely (i) the internal energy function that depends on the (gradient-extended) mechanical state and the entropy and (ii) a (nonsmooth) dissipation potential function that depends on the rates of the (gradient-extended) mechanical state and the entropy. Then, at current time the rates of the macro-and micro-motion and the entropy rate are governed by a canonical saddle point principle in case of an adiabatic process and a canonical minimization principle in case of a process including heat conduction. However, since specific forms of both constitutive functions are difficult to access, (generalized) Legendre transformations of those functions are considered. Here, the concept of dual variables is rigorously pursued, i.e., the quantity conjugate to the entropy (the temperature) is distinguished from the quantity conjugate to the entropy rate (the thermal driving force). Then, the necessary condition of the (generalized) Legendre transformation of the dissipation potential function must render the structure of the energy equation in form of an entropy evolution equation. The obtained rate-type variational principle is in line with Yang et al. (2006) but derived in an alternative way and extended by an additional long-range micro-motion that accounts for effects related to length scales. In addition, extended rate-type vari-Variational Thermomechanics of Gradient-Extended Continua 4 ational principles are constructed that incorporate dissipative mechanical driving forces and a temperature-dependent scalar threshold function via a Lagrange multiplier. By construction of consistent numerical time integration algorithms, variational principles valid for the time increment under consideration are obtained. Here, we especially point out a new semi-explicit numerical algorithm that, in contrast to a fully implicit scheme, evaluates the integration factor of Yang et al. (2006) to the value one and, hence, gives the algorithmically consistent form of the intrinsic dissipation.
In the subsequent sections, the above aspects are carried out in detail and the novel key contributions of the presented work for the thermomechanics of dissipative multifield problems are • a canonical variational principle based on a dissipation potential function that depends on the entropy rate, • the construction of rate-type and incremental mixed variational principles via (generalized) Legendre transformations, where we propose a new semi-explicit numerical update scheme that has the structure of an operator split and is highly suitable when considering adiabatic processes, • a general and versatile framework for the full thermomechanical coupling in gradientextended dissipative continua and • applications to gradient plasticity, gradient damage and Cahn-Hilliard diffusion, where for the latter we propose a new minimization structure with respect to the species flux vector.
The paper is organized as follows: In Section 2 we consider as motivating example a very simple adiabatic thermomechanical material element and derive rate-type as well as incremental variational principles based on the entropy and entropy rate as canonical variables in the energy and dissipation potential functions. Two update schemes, namely the known implicit and a new semi-explicit numerical algorithm are presented that govern the time-discrete evolution of the thermomechanical state. In Section 3 we set up general forms of governing equations for the thermomechanics of gradient-extended dissipative solids in a three-dimensional large-strain continuum setting. Section 4 shows that this fully coupled system is related to a variational statement. In addition, we construct extended rate-type and incremental variational principles that incorporate dissipative mechanical driving forces and threshold mechanisms that depend on the thermal state. The formulations outlined in Section 4 are general and can be applied to a wide spectrum of problems in thermomechanics. As examples, we specify in Section 5 the developed setting to Cahn-Hilliard diffusion, gradient damage and (additive) gradient plasticity with strong coupling to temperature evolution. The former is given a new variational structure in terms of the species flux vector. Finally, to demonstrate the capability of the newly proposed semi-explicit variational update scheme, we show a numerical example that is concerned with adiabatic shear band localization in softening plasticity. PSfrag replacements Figure 1: Rheological material element with internal thermoelastic and dissipative (viscous) mechanisms, loaded by an external stress σ ext (t). The time-discrete evolution of the material's state (ε, q, η) is described by a finite sequence of incremental variational principles. We consider a very simple rheological material element depicted in Figure 1 describing a smooth linear thermo-visco-elastic behavior. The material element is understood to be embedded in a thermal environment characterized by the absolute temperature θ > 0. The spring exhibits thermoelastic behavior by Hooke's law coupled with thermal expansion. The dashpot device characterizes viscous response via Newton's law. In addition, the material element is able to store heat. We have the fundamental constitutive relationships Here, σ denotes as usual a stress, ε a strain and e an internal energy. The four involved material parameters are Young's modulus E, the thermal expansion coefficient α T , the viscosity H and the heat capacity C at constant internal and external deformation. Without loss of generality, these material parameters are assumed to be constant, i.e., they do not depend on the temperature. Finally, θ 0 stands for a reference temperature. Such a simple model is able to describe classical thermomechanical coupling effects, such as the Gough-Joule effect or thermal expansion, as well as heating due to dissipation in the dashpot. The mechanical state of the rheological device is described by the total strain ε and the internal variable q characterizing the strain of the inner dashpot. We summarize these two quantities in the mechanical state array: c = ( ε, q ) .
Identifying the strain in the spring device by the simple kinematic relationship ε e = ε − q, we define the well known thermoelastic free energy function Variational Thermomechanics of Gradient-Extended Continua 6 that is not defined for nonpositive temperatures and is concave in θ. In addition, we define the dissipation potential function associated with the two dashpot devices (2)

Governing System of Equations
The thermomechanical response of the material element within a time interval T = (0, t e ) ⊂ R + is fully governed by the two constitutive functions (1) and (2) and characterized by four equations (3)-(6). First, on the mechanical side, the equilibrium equation is a relationship between internal and external stresses. Second, the evolving internal strain state is governed by Biot's equation and characterizes an internal dissipative mechanism, see Biot (1965). Third, on the thermal side, the state equation determines the current entropy η or in an inverse manner the current temperature θ. Finally, the evolution of the entropy is governed by a nonnegative dissipatioṅ Equation (6) 1 together with the constitutive definition of the dissipation D is a form of the first law of thermodynamics, i.e., the balance of energy for an adiabatic process under consideration. The inequality (6) 2 is the second law of thermodynamics that is ensured a priori by a dissipation potential function φ(·) that is (i) nonnegative, (ii) zero in the origin φ(0 ) = 0 and (iii) convex, see, e.g., Frémond (2013). In case of a dissipation potential function that is positively homogeneous of degree ℓ, i.e., φ(γċ) = γ ℓ φ(ċ) for all γ > 0, the dissipation can be written alternatively as 1 D = ℓ φ(ċ). Aim of the subsequent treatment is the construction of variational principles that account for the governing equations (3)-(6) of the thermomechanical device. Starting point is the definition of a (maybe unusual) dissipation potential function that depends on the entropy rate.

Canonical Variational Principle with Entropy Variable
The main difficulty in the formulation of rate-type variational principles in thermomechanics is to account for the full update of the thermal state via (6). For the nonisothermal case, we start by modeling the internal energy e whose corresponding function naturally depends on the entropy η. Hence, our investigation of thermomechanics in generalized standard materials is based on assuming the constitutive functions e(t) = e(c, η) and v(t) = v(ċ,η; c, η) .
The internal energy function e determines the current internal energy stored in the material element. The newly introduced dissipation potential function v is assumed to be a function of the rates (ċ,η) and might additionally depend on the current thermomechanical state (c, η). Alternatively, instead of the current entropy η the dissipation potential function v can depend on the current temperature θ by making use of the constitutive relationship (5). Comparing v with the dissipation potential function φ defined in (2), we observe an additional thermal slot with the dependence onη. This assumed dependency is the key ingredient of the subsequent construction of variational principles in dissipative thermomechanics. Based on the two constitutive functions (7), we construct the rate-type potential at given thermomechanical state (c, η), that accounts for energetic as well as dissipative mechanisms. Next, we define an external load function p ext (ε; t) = σ ext (t)ε depending on the given external stress σ ext , that loads the material element, see Figure 1. Then, the rate of the thermomechanical state at current time t under adiabatic conditions is determined by the canonical saddle point principle Note that the saddle point structure of this variational principle is governed by an assumed concavity of the dissipation potential function v with respect toη, see Section 2.3.2 and Footnote 2. Taking the first derivative of the potential (8) yields as necessary conditions of the variational principle (9) three Euler equations for the rate of the thermomechanical state of the material element at current time t, namely 1. Evolving external state ∂ε p ≡ ∂ ε e + ∂ε v = σ ext , 2. Evolving internal state ∂q p ≡ ∂ q e + ∂q v = 0 , 3. Evolving thermal state ∂η p ≡ ∂ η e + ∂η v = 0 .
The first equation determines the rateε of the external strain state and is a form of the stress equilibrium (3) without inertia terms. The second equation governs the ratė q of the internal variable and is identical to (4). Finally, the third equation is a form of the balance of energy, i.e., the first law of thermodynamics (6) combined with an inverse representation of (5), which determines the rateη of the entropy. However, this is a quite unusual representation. The main difficulty in this canonical setting is the formulation of the two constitutive functions e and v in terms of the entropy η and entropy rateη, respectively 2 . Hence, a formulation in terms of the temperature is more convenient.

Mixed Variational Principle with Temperature Variable
For practical applications, the above canonical variational principle is transformed into a mixed saddle point principle that incorporates the temperature as a primary variable. This transformation is achieved within two steps.

Variable dual to Entropy
First, we define the internal energy function occuring in the rate-type potential (8) by a partial Legendre transformation where ψ is the free energy function depending on the thermal variable θ. The necessary optimality condition of this Legendre transformation defines the entropy in terms of the thermal variable θ. The latter is the dual quantity to the entropy η, i.e., the temperature of the material element.

Variable dual to Entropy Rate
In a second step, we define the dissipation potential function v occuring in the rate-type potential (8) by a partial Legendre transformation v(ċ,η; c, η) = inf in terms of the dissipation potential function φ, that depends on the thermal variable T and, by making use of (11), is defined at current state (c, θ). The necessary condition of this Legendre transformation defines the evolution of the entropẏ in terms of the thermal variable T . We identify T as the quantity dual to the entropy rateη and call it the thermal driving force that, up to this point, is rigorously distinguished from the temperature θ of the material element. For the adiabatic case under consideration, (13) must have the form of the balance of energy (6) and we identify This relationship restricts the form of the functional dependence of the dissipation potential function φ and is fulfilled if a multiplicative dependence on the thermal driving force T is assumed To arrive at (14), we use the result that at equilibrium the thermal driving force T can be identified with the temperature θ, which is an outcome of the mixed variational principle (17) below. The scaling factor T /θ first arised via a symmetry argument as integrating factor in the seminal work of Yang et al. (2006).

Mixed Variational Principle
Starting from the rate-type potential (8), the two partial Legendre transformations (10) and (12) motivate the definition of a mixed rate-type potential 3 π(ċ,η,θ, T ) : at given thermomechanical state (c, η, θ) in terms of the two constitutive functions ψ and φ. Inserting the necessary condition (11) on the given thermomechanical state yields the reduced mixed rate-type potential at given thermomechanical state (c, η, θ). Based on this definition, we obtain the mixed saddle point principle that determines at current time t the rates of the external strain, internal variable and entropy as well as the thermal driving force. The Euler equations of this saddle point principle are 1. Evolving external state ∂ε π red ≡ ∂ ε ψ + ∂ε φ = σ ext , 2. Evolving internal state ∂q π red ≡ ∂ q ψ + ∂q φ = 0 , 3. Thermal driving force ∂η π red ≡ T − θ = 0 , 4. Evolving thermal state ∂ T π red ≡ −η + ∂ T φ = 0 .
Again, the first equation represents the stress equilibrium and the second one Biot's equation. The third relationship identifies the thermal driving force with the given temperature as reported in Yang et al. (2006). The last equation represents the first law of thermodynamics in the form of an evolution equation for the entropy. Clearly, within this mixed setting the entropy rateη plays the role of a Lagrange multiplier. With known rates (ċ,η) of mechanical state and entropy, the temperature rate at current time t can be computed by taking the time derivative of (5). A constrained minimization formulation dual to the variational principle (9) is obtained by commuting the order of the infimum and the supremum in (17) (which we assume to be possible) such that For a general discussion of constructing mixed and dual variational principles we refer to the books Maugin (1992), Nguyen (2000) and Berdichevsky (2009).

Incremental Variational Principles
We consider a finite time interval [t n , t n+1 ] ⊂ T with step length τ := t n+1 − t n > 0 and assume all thermomechanical state quantities at time t n to be known. The goal is then to determine all the approximate state quantities at time t n+1 based on variational principles valid for the time increment under consideration. Subsequently, all variables without subscript are understood to be evaluated at time t n+1 .

Mixed Variational Principle
A mixed variational principle associated with the finite time interval [t n , t n+1 ] is based on the incremental potential in terms of the continuous rate-type potential π defined in (16). Here, Algo stands for a numerical integration algorithm in time within the interval [t n , t n+1 ]. The incremental potential π τ is understood to be defined at given thermomechanical state (c n , η n , θ n ) at time t n . Then, the finite-step sized incremental mixed variational principle determines the mechanical state, the entropy, the temperature and the thermal driving force at current time t n+1 . The Algo operator is constructed such that the incremental variational principle (20) yields as Euler equations consistent algorithmic counterparts of the Euler equations (18) stemming from the rate-type variational principle (17). In what follows two different forms of the Algo operator are discussed.

Implicit Variational Update
As a first approach, we consider a numerical algorithm that performs an exact incrementation of the internal energy and an incrementation of the dissipative term according to a backward Euler scheme see Yang et al. (2006). Here,ċ τ = (c − c n )/τ denotes an algorithmic expression of the rate of the mechanical state. The incremental potential (19) takes the form Comparing this set of equations with (23), we observe essential modifications in (27) 3 and (27) 4 . These equations identify the current thermal driving force T with the given temperature θ n at time t n and define the current temperature θ in terms of the current mechanical state c and the given entropy η n at time t n . The key equations (27) 4 and (27) 5 on the thermal side can be recast into where T has been eliminated by (27) 3 and again a positively homogeneous dissipation potential function φ of degree ℓ is assumed. One observes, that in sharp contrast to the governing equation (24) 2 stemming from the implicit variational approach, equation (28) 2 does not contain the scaling factor. Hence, by the explicit variational update the algorithmically consistent form ℓ φ(ċ τ ) of the dissipation (6) 2 is obtained. Furthermore, also the dissipative terms in the stress equilibrium (27) 1 and Biot's equation (27) 2 do not have the scaling factor. Finally, we observe that as a consequence of the temperature update (28) 1 for given entropy η n , a decoupled incremental formulation is at hand.

Operator Split
The two update equations (28) characterize a variational-based staggered algorithmic treatment of dissipative thermomechanics based on the potential π τ defined in (26): (i) update the mechanical state c as well as the temperature θ at frozen entropy η n and given (predicted) thermal driving force T = θ n , and (ii) update the entropy η as well as the thermal driving force T which turns out to be identical to the predicted one. Hence, the semi-explicit variational update in the time interval [t n , t n+1 ] can be seen as a composition of two fractional steps Algo = Algo η,T • Algo c,θ , which are both variational. This numerical algorithm falls under the category of incrementally isentropic operator splits, a specific form of which is considered in . Within the isentropic predictor step we first solve the variational sub-problem that gives the mechanical state and the temperature at time t n+1 . The resulting optimality conditions are the equations (27) 1−2 and (27) 4 . Within a second step, the entropy corrector, the entropy as well as the thermal driving force at time t n+1 are obtained via the variational sub-problem for given mechanical state c * and temperature θ * stemming from the isentropic predictor (29). The resulting optimality conditions are the equations (27) 3 and (27) 5 .

Thermomechanics of Gradient-Extended Dissipative Solids
In this section, we shortly summarize the basic constitutive framework for the thermomechanics of gradient-extended dissipative solids in an abstract setting and adopt the notion of Miehe (2011Miehe ( , 2014. Let B ⊂ R d with d ∈ {2, 3} be the reference placement (d-manifold) of a material body B into the Euclidean space with smooth boundary ∂B.
We study the thermomechanical behavior of the body under mechanical as well as thermal loadings in a time interval T = (0, t e ) ⊂ R + . To this end, we focus on a multiscale viewpoint in the sense that we relate dissipative effects at (X, t) ∈ B × T to changes in the microstructure (besides the dissipative effect stemming from heat conduction). In the phenomenological context, we account for these microstructural mechanisms by micromotion fields that generalize the classical concept of internal variables. The evolution of the related dissipative system is considered to be nonsmooth and therefore the governing equations are written as differential inclusions. Without loss of generality, we assume within our treatment homogeneous solids only, i.e., material properties do not depend on the position X ∈ B. In what follows, we denote by( ·) = ∂ ∂t (·)(X, t) the material time derivative and by ∇(·)(X, t) the gradient on the reference manifold B. In addition, to keep notation short, we understand the operation a · b either as duality product between a vector and a one-form a · b = a c b c or as inner product between two vectors a · b = a c b d δ cd and two one-forms a · b = a c b d δ cd , respectively. The Kronecker deltas δ ab and δ ab represent the coefficients of the metric and inverse metric tensor in a Cartesian coordinate system, respectively. Finally, we give the usual definition of a double contraction between two second-order tensors, e.g., for two mixed-variant tensors A :

Macro-Motion of a Body
Within the geometrically nonlinear theory, the macroscopic motion of the body is described by the macro-motion field ϕ : which maps at time t ∈ T points X ∈ B of the reference configuration B to points x ∈ ϕ t (B) of the current configuration ϕ t (B). Let G = δ AB E A ⊗E B and g = δ ab e a ⊗e b denote the Euclidean metric tensors associated with the reference and current configuration, where the Kronecker symbols refer to Cartesian coordinate systems on both manifolds. The deformation gradient F = Dϕ(X, t) is defined as the tangent corresponding to the macro-motion (30). Next, the convected current metric C = F T gF is introduced, that is often denoted as the right Cauchy-Green tensor. The macro-motion (30) is constrained by a positive Jacobian J = √ det C > 0 that rules out penetration of matter. With respect to mechanical loading of the solid, the boundary of the reference configuration is decomposed into nonoverlapping parts ∂B ϕ and ∂B t such that ∂B = ∂B ϕ ∪ ∂B t . We prescribe the macro-motion (Dirichlet boundary condition) and the macro-tractions (Neumann boundary condition) on ∂B t × T as specified in Section 3.3 below.

PSfrag replacements
macro-motion field micro-motion field thermal driving force Figure 2: Primary fields in thermomechanics of gradient-extended dissipative solids. At time t, the macro-motion field ϕ(X, t) is constrained by Dirichlet and Neumann boundary conditions ϕ =φ on ∂B ϕ and P n 0 ∋t on ∂B t . The long-range micro-motion field q l (X, t) is at time t restricted by the conditions q l =q l (X) on ∂B q and H l n 0 ∋ 0 on ∂B H . Finally, for a heat conducting solid the thermal driving force field T (X, t) is at time t constrained by the conditions T =T on ∂B T and q · n 0 =q on ∂B q .

Micro-Motion of a Body
Microstructural dissipative changes of the body are described by additional fields related to the concept of internal variables. These variables are assembled in the micro-motion field of the solid q : The array q with in total δ scalar valued entries may contain internal variables of any tensorial rank that describe in a homogenized sense the micro-motion of the material due to structural changes on lower scales. Classical examples for members of q are damage variables, plastic strains or phase fractions. In addition, we assume all tensorial elements of q to be Lagrangian objects, i.e., they are not affected by rigid body macro-motions superimposed on the current configuration ϕ t (B). Like in Miehe (2011Miehe ( , 2014, we distinguish between long-range variables q l (order parameters) that are governed by PDEs in form of additional balance equations and connected to given length scale parameters, and short-range variables q s that are determined by ODEs and represent the standard concept of local internal variables. We summarize both, long-as well as short-range variables in the array q = (q s , q l ) .
For the treatment of each long-range micro-motion field, we decompose the boundary of the reference configuration into nonoverlapping parts ∂B q and ∂B H such that ∂B = ∂B q ∪ ∂B H . We prescribe the micro-motion (Dirichlet boundary condition) and the micro-tractions (Neumann boundary condition) on ∂B H × T as specified in Section 3.2.2 below. Note that the given micro-motion (32) on the Dirichlet boundary is considered to be independent of time. This is because we assume these variables to be "passive" in the sense that the deformation process cannot be driven by externally applying them.

Thermal Driving Force
As a third primary field in the thermomechanics of dissipative materials, we introduce the (macroscopic) thermal driving force field T : As outlined in the motivating Section 2.3.2, the thermal driving force appears as the dual quantity to the entropy rate and is identified as the temperature. It is governed by the generalized Legendre transformation (48) introduced below. With respect to thermal loading of a heat conducting solid, we decompose the boundary of the reference configuration into nonoverlapping parts ∂B T and ∂B q such that ∂B = ∂B T ∪ ∂B q . We prescribe the thermal driving force (Dirichlet boundary condition) and the heat flux (Neumann boundary condition) on ∂B q × T as specified in Section 3.3 below. The primary fields, namely the macro-motion, the long-range micro-motion and, in case of heat conduction, the thermal driving force are depicted in Figure 2.

Free Energy Function
The free energy storage in continua is governed by a free energy function. We specify it by focusing on simple materials of grade one, i.e., we include as arguments the first gradients of the micro-and macro-motions 4 defining the mechanical constitutive state. Together with the temperature, this state is invariant under any rigid body macro-motion superimposed on the current configuration ϕ t (B), i.e., the function (34) 1 is objective. We define the energetic contribution to the (contravariant) first Piola-Kirchhoff stress tensor P as P e = 2F ∂ C ψ , and the well known Clausius-Planck inequality takes the form Here, P d = P − P e denotes the dissipative part of the first Piola-Kirchhoff stress tensor.

Dissipation Potential Functions
Dissipative mechanisms are described by dissipation potential functions. First, we take mechanical effects into account via an objective intrinsic dissipation potential function at given thermomechanical state (c, θ). For generality, we assume in this abstract setting φ int to be nonsmooth with respect toċ, e.g., positively homogeneous of degree one, φ int (γċ; c, θ) = γ φ int (ċ; c, θ) for all γ > 0, which characterizes in the adiabatic case a rate-independent evolution, see also Footnote 7. With the intrinsic dissipation potential function at hand, we define the dissipative stress that arises in the reduced Clausius-Planck inequality (35). As suggested in Miehe (2014), we assume the evolution equations 5 for the micro-motions together with the Neumann boundary conditions In view of the nonsmoothness of the dissipation potential function (36), we write the microand macroscopic balance equations (37) and (40) as differential inclusions, i.e., ∂ċ(·) is understood as subdifferential. Considering the global form of (35), doing integration by parts two times and taking into account the homogeneous rate-type Dirichlet boundary conditionq = 0 on ∂B q according to (32), yields the alternative representation for the intrinsic dissipation, see again Miehe (2014). Note that the micro-force balance (37) 1 and the homogeneous Neumann boundary condition (37) 2 are outcomes of the variational principle (59) set up below. 6 As a second contribution to entropy production, we take into account heat conduction via an objective dissipation potential φ con (X, t) = φ con ( ; c, θ) with = −∇θ/θ at given thermomechanical state (c, θ). With such a function at hand, we define the material heat flux vector q = ∂ φ con , and the well known Fourier inequality takes the form The inequalities (38) and (39) serve as fundamental physically-based constraints on the dissipation potential functions φ int and φ con . These two conditions are satisfied a priori for dissipation potential functions that are (i) nonnegative φ(· ; c, θ) ≥ 0, (ii) zero in the origin φ(0 ; c, θ) = 0 and (iii) convex inċ and , respectively.

Governing Equations for Thermomechanics of Gradient-Extended Solids
We now summarize the governing equations for the thermomechanics of gradient-extended dissipative solids undergoing large deformations. First, we have the balance of linear momentum whereγ(X, t) andt(X, t) are prescribed body force and nominal surface traction fields. The balance of angular momentum P F T = F P T is a priori satisfied by the objectivity of the free energy function (34) and the dissipation potential function (36). The evolving micro-motion is governed by the micro-force balance equation as already given in (37). On the thermal side, the entropy is defined by the local state equation Finally, the evolution of the entropy is governed by the energy equation wherer(X, t) andq(X, t) are prescribed heat source and material surface heat flux fields. Note that if φ int (· ; c, θ) is positively homogeneous of degree one 7 , ∂ċ φ int is a set, but the expression ∂ċ φ int ·ċ, which represents the intrinsic dissipation, evaluates to a unique value since it coincides with φ int (ċ; c, θ), see also Section 2.1.2. The equations (40)-(43) generalize the ODEs (3)-(6) for the rheological device to the large-strain continuum setting and include intrinsic gradient-type dissipative effects as well as heat conduction.

Variational Principles for Thermomechanics of Gradient-Extended Solids
We now discuss the variational structure of thermomechanics of gradient-extended dissipative solids undergoing large deformations, whose Euler equations were summarized in Section 3.3. Again, the starting point is the definition of canonical energy and dissipation potential functions. By means of (generalized) Legendre transformations, different rate-type and incremental mixed formulations are derived.

Canonical Energy and Dissipation Potential Functionals
We generalize the variational framework for the rheological model presented in Section 2 to the large-strain continuum setting of gradient-extended dissipative solids. To this end, based on the internal energy and dissipation potential functions e(X, t) = e(c, η) and v(X, t) = v(ċ,η; c, η, X, t) , we introduce the energy and dissipation potential functionals E represents the internal energy stored in the entire body B due to coupled micro-macro deformations and thermal effects. The introduced functional V is related to intrinsic dissipative mechanisms and entropy production due to heat conduction. In analogy to (7), the entropy η as well as the entropy rateη are used as canonical thermal variables.
On the mechanical side, the constitutive functions (44) are assumed to depend on the constitutive state array c defined in (34) 2 that makes those functions a priori objective. Note, that the dissipation potential function v in general depends on the current state (c, η) as well as explicitly on position and time (X, t) ∈ B × T stemming from a possible inhomogeneous and time-dependent thermal loading, see below.

Energy and Dissipation Functionals in terms of Temperature
For practical modeling, we transform the above energy and dissipation potential functionals into functionals that depend additionally on the temperature. This we do in analogy to Section 2.3.
The necessary optimality condition of (46) is the statement (42) which identifies the thermal state variable θ as the dual quantity to the entropy η, i.e., as the temperature.

Variable dual to Entropy Rate
In a second step, we define at time t the dissipation potential functional by a generalized Legendre transformation 8 in terms of the mixed dissipation potential functional governed by a dissipation potential function φ and the thermal load functional Note that by use of the local state equation (42), the mixed dissipation potential functional V + depends on the current temperature. The maximization in (48) at time t is performed under the constraint T =T on ∂B T and defines as necessary condition the evolution of the entropy along with a thermal boundary conditioṅ These equations generalize the local statement (13) for the rheological device to a largestrain continuum setting including intrinsic gradient-type dissipative effects as well as heat conduction. As in Section 2.3.2, we call T the thermal driving force that is dual to the entropy rateη. The entropy evolution (51) 1 must have the form (43) 1 and we identify The first of these conditions is fulfilled for T = θ in B if the dissipation potential function is specified as which is in line with Yang et al. (2006), but derived in an alternative way. The second condition (52) 2 is satisfied for a volumetric thermal loading function It remains to find an expression for the thermal surface load function s T . Note that −∂ ∇T φ con = 1/θ ∂ φ con for T = θ in B and we identify from (43) 2 which is fulfilled for a thermal surface load function 8 The generalized Legendre transformation is conceptually in line with Miehe et al. (2014a).

Load Functionals
Besides the mixed energy and dissipation potential functionals (47) and (49), we have an external thermomechanical load functional with decoupled mechanical and thermal contributions. On the mechanical side, we define the load functional in terms of a given body force fieldγ and nominal surface traction fieldt. The thermal load functional (50) attains with the identifications (54) and (55) the form in terms of a given heat source fieldr and material surface heat fluxq.

Fundamental Mixed Variational Principle for Thermomechanics 4.3.1 Rate-type Formulation
Based on the internal energy and dissipation potential functionals E + and V + and the thermomechanical load functional P ext , we define at current time t the rate-type potential 9 with given state (ϕ, q, η, θ). We write this potential with its internal and external contributions Π + (φ,q,η,θ, in terms of the internal potential density Recalling the mixed functions (47) 2 and (49) 2 together with (53) and inserting the necessary condition (42) on the given thermomechanical state, yields a reduced internal potential density of the form Then, the rates of the macro-and the micro-motion as well as the rate of the entropy and the thermal driving force at current time t are governed by the variational principle 10 {φ,q,η, T } = Arg{ inḟ Here, one has to account for the rate forms of the Dirichlet boundary conditions (31) and (32) for the macro-and micro-motions, i.e., ϕ =φ on ∂B ϕ andq = 0 on ∂B q as well as for the Dirichlet boundary condition (33) for the thermal driving force. The variational principle (59) is as stated in Yang et al. (2006), but extended by a long-range micro-motion field. By the first variation of the functional (57), we have the necessary optimality conditions for all admissible test functions δη and (δφ, δq, δT ) fulfilling homogeneous forms of the Dirichlet boundary conditions. We get the Euler equations 1. Evolving macro-motion δφ π + red ≡ δ ϕ ψ + δφ φ int ∋ gγ ,

Evolving thermal state
in B along with the Neumann boundary conditions on ∂B t , ∂B H and ∂B q , respectively. In contrast to (18), the equations are now exclusively governed by variational derivatives of the reduced potential density π + red defined in (58). The central three field equations are the quasi-static mechanical equilibrium that governs the rateφ of the macro-motion, the micro-force balance determining the rateq of the micro-motion and the energy equation for the evolutionη of the entropy.

Incremental Formulation
Consider a finite time interval [t n , t n+1 ] ⊂ T with step length τ = t n+1 − t n > 0 and assume all thermomechanical field variables at time t n to be known. The goal is then to determine all the approximate fields at time t n+1 based on variational principles valid for the time increment under consideration. Subsequently all variables without subscript are understood to be evaluated at time t n+1 . We may formulate the incremental potential where E +τ , V +τ and P τ ext are incremental energy, dissipation and load functionals associated with the time interval [t n , t n+1 ]. These functionals are defined at given state (ϕ n , q n , η n , θ n ) at time t n . In analogy to the rate-type formulation (57), we rewrite the incremental potential in terms of an incremental internal potential density π +τ which is defined at given state (c n , η n , θ n ). Such a function is obtained by a numerical integration algorithm that has to be constructed in such a way that the subsequent incremental variational principle gives consistent algorithmic counterparts of the Euler equations (61). In what follows, we construct an implicit as well as a semi-explicit numerical integration algorithm, compare Section 2.4. As a typical example we consider an integration using the approximations of the rates of state quantitieṡ and the incremental internal potential density (67) with k = n + 1 for an implicit numerical integration algorithm according to (21) and k = n for a semi-explicit numerical integration algorithm according to (25). In (67) we dropped terms that are associated with previous time t n . Then, defining the incremental load functional P τ ext (ϕ, T ; ϕ n , θ n , t n+1 ) = P ϕ ext (ϕ − ϕ n ; t n+1 ) + τ P T ext (T ; θ n , t n+1 ) in B along with the Neumann boundary conditions on ∂B t , ∂B H and ∂B q , respectively. As a fundamental difference to the fully implicit numerical algorithm, a semi-explicit update identifies the thermal driving force with the given temperature at time t n . Hence, the scaling factor results in T /θ n = 1 in B such that the algorithmically consistent form of the intrinsic dissipation defined in (35) and reformulated in (38) is obtained. Especially, the incremental energy equation reads Additionally, also the dissipative terms in the quasi-static equilibrium (69) 1 , micro-force balance (69) 2 and the boundary conditions (70) 1−2 do not contain the scaling factor. However, there are two issues that arise when using the proposed semi-explicit update for a heat conduction process: (i) on the thermal side it might be restricted to homogeneous Neumann boundary conditions (70) 3 on the whole boundary, i.e.,q = 0 on ∂B and (ii) the Courant-Friedrichs-Lewy (CFL) condition gives, depending on the mesh size associated with the spatial discretization, an upper bound for the time step size τ in order to obtain a stable numerical solution. 11 Note that the semi-explicit update can be seen as an incrementally isentropic operator split that consists of two fractional steps Algo = Algo η,T • Algo ϕ,q,θ .
First, in the isentropic predictor step we optimize the incremental potential (64) with respect to the macro-and micro-motions ϕ and q as well as the temperature θ, i.e., where the entropy is frozen. Then, the entropy η and the thermal driving force T are updated via the entropy corrector step (Algo η,T ) : {η * , T * } = Arg{ stat η,T Π +τ (ϕ * , q * , η, θ * , T ) } .
11 A derivation of the CFL condition for the suggested semi-explicit update scheme in case of a heat conduction process is beyond the scope of this paper. For a stability analysis of fractional step methods in thermomechanics we refer to .

Mixed Variational Principle with Mechanical Driving Forces
We now put the focus on mixed variational principles for the thermomechanics of gradientextended dissipative continua, where not only the microstructural variables, but also the corresponding local driving forces are taken into account and considered as additional variables. Following Miehe (2011Miehe ( , 2014, we consider the equivalent representation of the intrinsic dissipation (38) by the dual intrinsic dissipation potential function φ * int depending on the mechanical driving forces f = (m, d, g) conjugate to c = (C, q, ∇q) .
The Legendre transformation motivates the definition of an extended mixed dissipation potential functional in terms of the mixed dissipation potential function which governs the subsequent extended mixed variational principle. The necessary condition of (72) reads 12ċ ∈ ∂ f φ * int (f; c, θ) . 12 Perzyna-type dual dissipation potential function. An important example of a smooth intrinsic dual dissipation potential function is in terms of a function f (· ; c, θ) that differs from a gauge just by an additive constant and serves as a threshold function in the (adiabatically) rate-independent setting considered in Section 4.5 below. η f > 0 is a material parameter and · + : R → R + , x → 1 2 (|x| + x) the ramp function. Note carefully that φ * int (· ; c, θ) defined in (75) is not a homogeneous function and its dual φ int (· ; c, θ) is still nonsmooth in general, i.e.,ċ = ∂ f φ * int but δ ϕ ψ + δφ φ int ∋ gγ and δ q ψ + δq φ int ∋ 0 .
We can write (76) 1 in the particular forṁ which regularizes the rate-independent structure (90)-(91) to be discussed later in Section 4.5. This rate-dependent setting is in line with the formulation of visco-plasticity according to Perzyna (1966).

Rate-type Formulation
Based on the internal energy and dissipation potential functionals E + in (47) and V * in (73), we are in the position to formulate a mixed rate-type variational principle that accounts for the mechanical driving forces f. We define at current time t the rate-type potential 13 with given state (ϕ, q, η, θ). We write this potential with its internal and external contributions in terms of the extended internal potential density Recalling the mixed functions (47) 2 and (74) and inserting the necessary condition (42) on the given thermomechanical state, yields a reduced internal potential density of the form Then, the rates of the macro-and micro-motion as well as the rate of the entropy and the thermal and mechanical driving forces at current time t are governed by the variational principle {φ,q,η, T, f} = Arg{ inḟ Like in (59), one has to account for the Dirichlet boundary conditions (60) and (33). By the first variation of the functional (77), we have the necessary optimality conditions for all admissible test functions (δη, δf) and (δφ, δq, δT ) fulfilling homogeneous forms of the Dirichlet boundary conditions. We obtain the Euler equations 3. Thermal driving force ∂η π * red ≡ T − θ = 0 ,

Evolving thermal state
in B along with the Neumann boundary conditions on ∂B t , ∂B H and ∂B q , respectively. The central three field equations are the quasi-static equilibrium the micro-force balance and the energy equation They are complemented by the inverse definitions (79) 5 of the mechanical driving forces which split into three evolution equationṡ as already given in Miehe (2014), where the identity T = θ in B from (79) 3 has been used. Note carefully that the driving forces f are variables and the variational derivatives in (79) 1−2 with respect toφ andq are understood in the ordinary sense. Hence, (79) 1−2 as well as the boundary conditions (80) 1−2 do not represent differential inclusions.

Incremental Formulation
Within a time interval [t n , t n+1 ] ⊂ T a variational principle can be constructed by the same avenue as outlined in Section 4.3.2. Using the algorithmic approximations (66) of rates of state quantities, we get the incremental internal potential density where again k = n + 1 corresponds to a fully implicit and k = n to a semi-explicit update. In addition, we dropped terms that are associated with previous time t n . Then, with the use of the incremental load functional (68) the incremental variational principle determines all thermomechanical fields at time t n+1 . The corresponding Euler equations in B are time-discrete forms of (79) together with (69) 4 stemming from the variation with respect to the temperature θ. As before, the proposed semi-explicit integration of the rate of the internal energy yields for the scaling factor T /θ n = 1 in B and one obtains the algorithmically consistent form of the intrinsic dissipation, i.e., the incremental energy equation reads In addition, the dissipative terms in the time-discrete forms of the quasi-static equilibrium (79) 1 , the micro-force balance (79) 2 and the boundary conditions (80) 1−2 as well as the dissipative terms in the time discrete forms of the evolution equations (79) 5 do not contain the scaling factor. The isentropic operator split is modified by an additional optimization in the isentropic predictor step with respect to the mechanical driving forces f.

Mixed Variational Principle with Threshold Function
Intrinsic dissipation potential functions are often modelled by the principle of maximum dissipation. For the classical local theories of plasticity, this principle can be traced back among others to the work of Hill (1950), see also Moreau (1976), Simo (1988) and Lubliner (2008). For a general discussion of this principle and its connection to evolution laws governed by dissipation potentials we refer to Hackl & Fischer (2007) and Hackl et al. (2011b,a). The intrinsic dissipation potential function for thermomechanics is defined by the constrained maximum principle that includes at given state (c, θ) a domain of admissible mechanical driving forces which we assume to be smooth 14 Clearly, the function φ int (· ; c, θ) defined in (87) is positively homogeneous of degree one. The set (88) is governed by a threshold function f (f; c, θ) = w(f; c, θ) − c(c, θ), where c(c, θ) > 0 is a positive threshold constant that might depend on the given thermomechanical state and w a level set function that is a gauge, i.e., (i) nonnegative w(· ; c, θ) ≥ 0, (ii) zero in the origin w(0 ; c, θ) = 0, (iii) convex in f and (iv) positively homogeneous of degree one in f. As a result of (iii), the constrained optimization problem (87) has a unique solution. By the use of the Lagrange multiplier method, we can put (87) into the form whose necessary optimality condition defines the evolution of the constitutive statė where the Lagrange multiplier λ satisfies classical loading-unloading conditions in Kuhn-Tucker form This is the typical structure of flow rules associated with rate-independent material behavior. Note that even though ∂ f f (0 ; c, θ) is a subdifferential (see Footnote 7), we get for f = 0 the unique valueċ = 0 . The optimization problem (89) now motivates the definition of a modified mixed dissipation potential functional in terms of the mixed dissipation potential function that governs the subsequent modified mixed variational principle.
It should be mentioned that unlike for local theories, the Lagrange multiplier λ for f = 0 cannot be determined at current time by a local consistency condition in terms of rates of the external quantities deformation and temperature. Hence, as, e.g., mentioned in De Borst & Mühlhaus (1992) a nonlocal version has to be elaborated, see Section 4.5.1.

Rate-type Formulation
Based on the internal energy and dissipation potential functionals E + in (47) and V * λ in (92), we are in the position to formulate a mixed rate-type variational principle that accounts for dissipative threshold mechanisms. We define at current time t the rate-type potential with given state (ϕ, q, η, θ). We write this potential with its internal and external contributions in terms of the extended internal potential density π * λ (ċ,η,θ, T, ∇T, f, λ) = d dt e + (c, η, θ) + v * λ (ċ,η, T, ∇T, f, λ) .
Recalling the mixed functions (47) 2 and (93) and inserting the necessary condition (42) on the given thermomechanical state, yields a reduced internal potential density of the form Then, the rates of the macro-and micro-motion as well as the rate of the entropy, the thermal and mechanical driving forces and the Lagrange multiplier at current time t are governed by the variational principle Like in (59), one has to account for the Dirichlet boundary conditions (60) and (33). By the first variation of the functional (94) we have the necessary optimality conditions for all admissible test functions (δη, δf, δλ) with λ + δλ ≥ 0 in B and (δφ, δq, δT ) fulfilling homogeneous forms of Dirichlet boundary conditions. We obtain the Euler equations 1. Evolving macro-motion δφ π * λ,red ≡ δ ϕ ψ + δφ( T θ f ·ċ) = gγ , 2. Evolving micro-motion δ˙q π * λ,red ≡ δ q ψ + δ˙q( T θ f ·ċ) = 0 , 3. Thermal driving force ∂η π * λ,red ≡ T − θ = 0 , in B along with the Neumann boundary conditions (80). Note, that the condition f ≤ 0 in (98) 6 follows from (97) 3 if we set λ = 0, which necessarily demands δλ ≥ 0. On the other hand, by choosing λ > 0 the test function δλ can have any sign and we obtain from (97) 3 the equality f = 0, or in summary λ f = 0 as given in (98) 6 . The central three field equations are identical to (81)-(83) and are complemented by the set of evolution equations (98) 5 which reaḋ as already given in Miehe (2014), where the identity T = θ in B from (98) 3 has been used. Note that the only formal difference to (84) is that the nonsmooth evolution is not written in form of differential inclusions but is governed by the Karush-Kuhn-Tucker conditions (98) 6 containing the scalar variable λ. We also mention that the representation (99) is possible because of the assumption of a smooth elastic domain E, i.e., its boundary only consists of regular points.
In addition, we have for f = 0 the nonlocal consistency conditions which at current time t are supplemented by rate forms of the equations (42) and (98) 1−2 , the evolution equations (98) 5 and the rate forms of the Neumann boundary conditions (80) related to the macro-and micro-motion, respectively. To see conditions (100), consider at current time t a nonzero evolution of the mechanical constitutive state, i.e., λ t > 0. Then, the first variation of the dissipation potential functional (92) with respect to the Lagrange mulitplier vanishes for all δλ. Next, at time t + τ , τ ≥ 0 we consider a state with λ t+τ ≥ 0 and the first variation of the dissipation potential functional (92) with respect to the Lagrange multiplier is nonnegative for all λ t+τ + δλ ≥ 0. Subtracting (101) from (102) and dividing by τ yields for all λ t+τ + δλ ≥ 0. For τ → 0 we have λ t+τ → λ t > 0, O(τ 2 )/τ → 0 and get d dt f = 0 since δλ can have any sign. For τ small enough, we assume λ t+τ = 0 and δλ must be nonnegative yielding d dt f ≤ 0. When summarizing, we arrive at the nonlocal consistency conditions (100). From this condition the Lagrange multiplier field can be determined in terms of the rates (φ,θ) of the external fields deformation and temperature.

Incremental Formulation
Within a time interval [t n , t n+1 ] ⊂ T a variational principle can be constructed by the same avenue as outlined in Section 4.3.2. Using the algorithmic approximations (66) of rates of state quantities, we get the incremental internal potential density where like before k = n + 1 corresponds to a fully implicit and k = n to a semi-explicit update. Again, terms that are associated with previous time t n are dropped. Then, with the use of the incremental load functional (68) the incremental variational principle determines all fields at time t n+1 . The corresponding Euler equations in B are timediscrete forms of (98) together with (69) 4 stemming from the variation with respect to the temperature θ. Considering the semi-explicit integration, the incrementally isentropic operator split is modified by an additional optimization in the isentropic predictor step with respect to the Lagrange multiplier λ ≥ 0.

Representative Model Problems
We now apply the general variational setting for thermomechanics of gradient-extended dissipative solids to complex multifield problems related to evolutions of species content in elastic solids, damage and plasticity. Especially, we highlight a new minimization structure for Cahn-Hilliard diffusion with respect to the species flux. The focus is put on incremental potentials which also inherently contain the proposed operator split. The latter is applied to a numerical example which shows the formation of an adiabatic shear band.

Cahn-Hilliard Diffusion coupled with Temperature Evolution
We consider as a first model problem the Cahn-Hilliard theory of diffusive phase separation in a rigid solid coupled with temperature evolution. In the following c : B × T → [0, 1] denotes a dimensionless concentration field whose evolution is governed by the local form of conservation of species contentċ = − Div À , where À : B × T → R d is the species flux vector field. To give the concentration field the character of an order parameter q l = (c), we impose homogeneous boundary conditionṡ c = 0 on ∂B q and ∂ ∇c ψ · n 0 = 0 on ∂B H , see Figure 2. Hence, the dynamic process is only driven by an initial inhomogeneous distribution c 0 (X) of the concentration field in the domain B. In the following, we neglect the phenomenon of thermal diffusion (Soret effect) that is species flow caused by a temperature gradient, see De Groot & Mazur (2013) and the recent contribution Nateghi & Keip (2021). The free energy function decomposes into a local, nonlocal and purely thermal contribution where the last term is given in (1). As considered in Cahn & Hilliard (1958) we choose in terms of the threshold and mixing energy parameters A and B as well as the diffuse interface parameter L. Note the nonconvexity of ψ l for B > 2A which is related to phase separation. The evolution of the concentration is driven by the chemical potential µ given by It can be understood as a constitutive representation of a micro-force balance in the sense of Gurtin (1996).

Rate-type Minimization Principles in Isothermal Case
Point of departure is the definition of the energy and dissipation functionals where φ is the smooth dissipation potential function accounting for diffusion mechanisms. A minimization of the corresponding rate-type potential with respect toċ then yields as Euler equation the mass balance in the form Alternative to this setting, we now propose a new minimization formulation in terms of the species flux vector. In line with Miehe et al. (2015), we reformulate the rate of the energy functional (108) 1 at current concentration where we inserted the balance equation (105). The dissipation functional is defined as in terms of the dissipation potential function χ which has the simple quadratic form Here, M > 0 is a so-called mobility parameter. Note that χ(· ; c) is a positively homogeneous function of degree two and its image coincides with half the dissipation in a material element, see below. With the functionals (109) and (110) at hand, we define the potential Π(À ; c) = E(À ; c) + X(À ; c) with given concentration c. Its minimization with regard to homogeneous Dirichlet-type boundary conditions À · n 0 = 0 and − Div À = 0 on ∂B q determines the current species flux field. We obtain the Euler equations 1. Species flux They contain necessary compatibility conditions for the given concentration field. As a post-processing step, the current rateċ of the concentration is determined via (105). Starting from (35), we now calculate the dissipation whose unique source is diffusion where we performed integration by parts two times and inserted the balance (105) as well as the homogeneous boundary conditions. Using (112) 1 , we can express the dual to the species flux vector by the dissipation potential χ and obtain with (111) for the dissipation associated with a volume element D = 2 χ(À ; c) ≥ 0.

Incremental Formulation with Temperature Evolution
In addition to (111), we consider the conductive dissipation potential function φ con ( ; θ) = k θ( · )/2 , which is convex in = −∇θ/θ with k > 0 being the thermal conductivity. For the incremental setting, we consider the implicit update of the species concentration c = c n − τ Div À , and specify the incremental internal potential density (67) as π +τ = ψ(c n − τ Div À, ∇c n − τ ∇[ Div À ], θ) Here, χ is the dissipation potential function (111) related to diffusion mechanisms, where an additional temperature dependence M = M (θ) of the mobility parameter is taken into account. We obtain the Euler equations where we recall the definition of the chemical potential (107) together with (113).

Thermomechanics of Gradient Damage
We consider as a second application the thermomechanics of a gradient damage model with an elastic stage. The scalar micro-motion field d : B × T → [0, 1], referred to as damage variable, measures at a macroscopic point X ∈ B the ratio between an arbitrary oriented area of microcracks and a representative reference surface in which the mentioned crack surfaces are embedded, see, e.g., Lemaitre (1996). In this sense, a value d = 0 characterizes an unbroken state, whereas d = 1 represents a fully broken state. The irreversibility of micro-cracking is usually expressed by the inequality constraintḋ(X, t) ≥ 0 on the evolution of the damage variable. The mechanical constitutive state is specified as and contains the right Cauchy-Green tensor C, the damage variable as well as its first gradient. In addition, we introduce the elastic right Cauchy-Green tensor C e = F eT gF e that is based on the definition of an elastic, stress producing part F e = J −1/3 θ F of the deformation gradient in terms of a volumetric thermal expansion J θ = exp[3α T (θ − θ 0 )], see Lu & Pister (1975). One can then write C e by means of C as C e = J −2/3 θ C. A simple model of thermo-gradient-damage at large deformations may then be based on the objective free energy function where g(d) = (1 − d) 2 is a degradation function and ψ θ the purely thermal contribution as given in (1). Note that the gradient of damage does not arise in this constitutive function but will exclusively enter the dissipation potential function, see below. µ > 0 is the shear modulus and the second parameter δ > 0 models a weak compressibility. Note that only the full elastic energy storage is degraded. From (114) we obtain the tensor functions that represent constitutive relationships for driving forces.
Here, irreversibility of damage is ensured by the indicator function I + (ḋ) of the set of positive real numbers defined as with ∂ denoting the subdifferential. The parameter c(θ) > 0 is a temperature-dependent force-like threshold value for the onset of damage with d dθ c < 0 and l a length scale parameter. Note that (116) is a positively homogeneous function of degree one in (ḋ, ∇ḋ) and hence models for an adiabatic process a rate-independent evolution of damage. In addition, we have the conductive dissipation potential function where the thermal conductivity is a function of the damage variable and may take the simple form k(d) = g(d)k b in terms of the heat conduction coefficient k b > 0 of the undamaged bulk. With the constitutive functions (114), (116) and (119) at hand, we specify the internal potential density (58) as Then, the variational principle (59) determines at current time t the rates of deformation, damage and entropy as well as the thermal driving force and gives the Euler equations δφ π + red ≡ − Div gP e = gγ , in B. Observe, that the evolution of the entropy is driven by the rates of damage and gradient of damage. The nonsmooth evolution of the damage variable is governed by the differential inclusion (120) 2 . From there and the relation (120) 3 , we conclude for the determination of the rates of deformation and damage the nonlocal consistency conditions 15 whereβ e = ∂ 2 dC ψ :Ċ + ∂ 2 dd ψḋ + ∂ 2 dθ ψθ andθ follows from taking the time derivative of the state equation (42). In addition, the rate form Div[ ∂ F P e :Ḟ + ∂ d P eḋ + ∂ θ P eθ ] =γ of mechanical equilibrium (120) 1 has to be considered together with the boundary conditions [ ∂ F P e :Ḟ + ∂ d P eḋ + ∂ θ P eθ ] n 0 =ṫ and ∇ḋ · n 0 = 0 (121) on ∂B t and ∂B H , respectively. (121) are rate forms of the Neumann boundary conditions P e n 0 =t on ∂B t and ∇d·n 0 = 0 on ∂B H which are outcomes of the variational principle.
From (38) the intrinsic dissipation is found to be where we used integration by parts and the homogeneous boundary conditions. Hence, thermodynamic consistency is shown.

Incremental Formulation based on Threshold Function
We specify the array (71) of dissipative driving forces as where β is the quantity conjugate to d. The Legendre transform of the local part (117) 1 of the intrinsic dissipation potential function reads and enters the incremental internal potential density (85). Note that φ * l is the indicator function of the set (88) of admissible driving forces governed by the threshold function 16 The latter defines the local intrinsic dissipation potential function by the constrained optimization problem 16 An alternative intrinsic local dissipation potential function may be given by that in contrast to (117) 1 depends explicitly on the given damage state d. We obtain the threshold function f (β; d, θ) = β − c(θ)d , which, due to occurrence of the damage variable in the resistance term, corresponds to a model without an elastic stage, i.e., the damage starts to evolve from d = 0 at the instant of loading. Figure 3: Geometry of additive plasticity. a) Definition of the total Hencky strain tensor ε = 1 2 ln C in terms of the pull-back C of the standard current metric g. b) Definition of the plastic metric tensor G p on the reference configuration governing the plastic Hencky strain ε p = 1 2 ln G p whose evolution is given by a local flow rule. Then, ε e = ε − ε p is the elastic strain measure entering the free energy function.

Thermomechanics of Additive Gradient Plasticity
As third model problem, we consider a thermomechanically coupled formulation of additive gradient plasticity. Besides the standard metrics G and g, we introduce on the reference configuration the (covariant) plastic metric tensor G p ∈ Sym + (3) that evolves in time starting from the initial state G p (X, t 0 ) = G. Following Miehe et al. (2002) and as visualized in Figure 3, a Lagrangian elastic strain variable may be based on an explicit dependence on the right Cauchy-Green tensor C, that is the current metric pulled back to the reference configuration, and the plastic metric G p in an additive format where the total and plastic Hencky strain tensors ε = 1 2 ln C and ε p = 1 2 ln G p are introduced. Hence, instead of G p we consider in what follows the logarithmic plastic strain ε p as the local internal variable whose evolution from ε p (X, t 0 ) = 0 is governed by a standard flow rule. Note that within this framework the condition of plastic incompressibility det G p = 1 is equivalent to a standard statement of vanishing trace tr ε p = 0. We specify the mechanical constitutive state which contains a scalar hardening variable α as well as its first gradient. In the following, we focus on metal plasticity characterized by small elastic but large plastic deformations and consider the free energy function ψ(c, θ) = ψ e (ε, ε p ) + ψ p (α, θ) + 1 2 µl 2 |∇α| 2 + ψ e−θ (ε, θ) + ψ θ (θ) .
Here, ψ e is the purely elastic contribution that is assumed to have the quadratic form where κ > 0 and µ > 0 are the bulk and shear modulus, respectively. ψ p is an isotropic hardening function that also takes into account thermally induced softening. The gradient extension related to a length scale parameter l is assumed to affect the scalar hardening variable only. The coupled thermoelastic response is modelled by the constitutive function (1), where also the pure thermal contribution ψ θ is specified. Note, that the function (127), known as Hencky energy, is not polyconvex 17 with respect to the deformation gradient F = Dϕ in the sense of Ball (1976), but rank-one convex for a moderately high elastic deformation range, see Bruhns et al. (2001). Hence, it is applicable to the typical scenario of metal plasticity. With the free energy function (126) at hand, we can derive the tensor functions that represent constitutive relationships for driving forces. Note the occurrence of the Laplacian term µl 2 ∆α in (128) 3 that is in line with the approach of Aifantis, see, e.g., Aifantis (1987).

Incremental Formulation based on Threshold Function
We specify the array (71) of dissipative driving forces as where (s, β) are the quantities conjugate to (ε p , α). For von Mises-type gradient plasticity we define the yield function that restricts the set of admissible driving forces according to (88). y(θ) is a temperature dependent yield stress function with d dθ y < 0. The corresponding intrinsic dissipation potential function is defined by the constrained optimization problem The intrinsic dissipation follows as as Euler equations in B. This set of equations can again be reduced by expressing the dissipative driving forces s and β from (132) 2 and (132) 3 yielding for k = n the nonlocal (regularized) update equations ε p = ε p n − (τ /η f ) | β e | − 2/3 ( y(θ n ) + β e ) + Gβ e G/|β e | , α = α n + (τ /η f ) | β e | − 2/3 ( y(θ n ) + β e ) + 2/3 for the plastic Hencky strain as well as the hardening variable. A local finite strain thermoplasticity model that uses the same additive strain kinematics together with the plastic configurational entropy in the sense of Simo & Miehe (1992) is proposed by Ulz (2009). For a comparison of rate-independent and rate-dependent formulations in isothermal gradient-plasticity of Fleck-Willis-type we refer to Nielsen & Niordson (2019).

Numerical Example: Cross Shear Localization
For softening visco-plasticity of von Mises-type, we analyze the development of shear bands in a rectangular plate B = (0, L) × (0, 2L) with L = 50 mm subject to tensile loading under the condition of plane strain. The geometric setup is depicted in Figure 4. We use the viscous over-force formulation of the mixed large deformation setting from Section 5.3.1 above and specify the isotropic hardening function in (126) as ψ p (α, θ) = 1 2 h(θ)α 2 .
Here, h is a temperature dependent hardening function which together with the temperature dependent yield stress function is specified as PSfrag replacements

2L
Lū u X Y B Figure 4: Cross Shear Localization. Geometry and mechanical loading. The process of heat conduction is neglected. Due to the symmetry of the boundary value problem, only the top right quarter of the domain is discretized by finite elements. To trigger plasticity in the center, the initial yield stress y 0 is reduced by 3% in the dark grey element.
see Simo & Miehe (1992). The used material parameters are summarized in Table 1. To trigger plasticity in the center, the initial yield stress is reduced by 3% in the element shaded in dark grey in Figure 4. For simplicity, we neglect the effect of heat conduction which is a reasonable assumption in case of a fast formation of the shear band generated by a high loading rate. We stretch the specimen with a constant displacement ratė u = 5 mm/s within the time interval T = (0, 1) s that is divided into 1000 equal increments 18 . We use the proposed semi-explicit variational update with index k = n in the incremental internal potential density (133). Since heat conduction is not taken into account, the time step size τ is not restricted by the CFL condition. Due to the variational structure, the resulting stiffness matrix within a typical Newton-Raphson iteration step is symmetric. As (global) primary fields we take the macroscopic deformation ϕ, the scalar hardening/softening variable α and its dual driving force β. The temperature θ is calculated via the implicit local equation (132) 5 . Due to symmetry, only one quarter of the domain is discretized by 15 × 30, 20 × 40 and 25 × 50 quadrilateral finite elements. We use a Q 1 E 4 -Q 1 -Q 1 element pairing which bases on a (local) enhancement of the macroscopic displacement gradient according to . Figure 5 shows contour plots of the equivalent plastic strain α and the temperature θ at final displacementū = 5 mm for three different plastic length scale parameters l ∈ {0.05, 0.1, 0.2} mm. Clearly, the specimen experiences a rise in the temperature in the region of plastic dissipation. When increasing the plastic length scale, the equivalent plastic strain α as well as the temperature θ spread over more elements. At the same time, one observes decreasing maximum values of α and θ, see also Aldakheel & Miehe (2017) and the references cited therein, 18 Since the process of heat conduction is neglected and the viscosity is chosen very low in order to have a formulation that is close to the nonsmooth setting, the overall thermomechanical material behavior is de facto rate-independent. Hence, the specific loading rate applied on the specimen is practically irrelevant.
Variational Thermomechanics of Gradient-Extended Continua i.e., Voyiadjis & Faghihi (2012). As widely known, in case of a local theory (l = 0 mm) the plastic deformation localizes within one element width. This mesh dependency also manifests itself in the load-displacement curve of the structure as shown in Figure 6a). In contrast, the regularization provided by the gradient theory yields mesh independent results, see Figure 6b). There, one also observes the additional softening effect in the nonisothermal case due to locally decreasing yield stresses according to (134) 2 . At this point, we want to note that Wcis lo & Pamin (2017) incorporate a gradient-enhanced relative temperature field in their formulation in order to regularize adiabatic thermal softening behavior.
For the used mixed setting of gradient plasticity, alternative finite element formulations are presented in Miehe et al. (2014b). Note that in this context the construction of infsup stable finite element pairings is a difficult task which in detail has recently been investigated by Krischok & Linder (2019). Especially, our chosen element pairing results in a nonphysical oscillatory behavior of the driving force field β. However, this instability seems to have no visible influence on the macroscopic deformation field ϕ, the scalar plastic strain field α and the temperature field θ.

Conclusion
We presented a unified framework for the fully coupled thermomechanics of gradientextended dissipative solids which undergo large deformations. The key of the formulation is the consideration of the entropy and the entropy rate as canonical variables which enter besides the gradient-extended mechanical state and the rate of this state, respectively, the internal energy and dissipation functions. Here, the rate-type formulation of local thermoplasticity outlined in Yang et al. (2006) is recovered by the concept of dual variables. Especially, the quantity conjugate to the entropy is rigorously distinguished from the quantity conjugate to the entropy rate. The coupled macro-and micro-balances as well as the energy equation are outcomes of the stated saddle point principle. Within this setting, the entropy is also driven by additional gradient-type dissipative effects. Emphasis was also put on extended variational formulations which incorporate dissipative mechanical driving forces and temperature dependent threshold mechanisms. In addition, we discussed variationally consistent time incrementations and suggested a semi-explicit numerical algorithm that renders the algorithmically consistent form of the intrinsic dissipation. It was shown that this algorithm has the structure of an operator split. A special focus was put on applying the framework to complex multifield problems which are of interest in the thermomechanics of solids. Three model problems, i.e., Cahn-Hilliard diffusion, gradient damage and (additive) gradient plasticity strongly coupled with temperature evolution, showed the capability of the proposed formulation. In a numerical example we studied the formation of a cross shear band under adiabatic condition.