Structure-preserving space-time discretization of large-strain thermo-viscoelasticity in the framework of GENERIC

Large-strain thermo-viscoelasticity is described in the framework of GENERIC. To this end, a new material representation of the inelastic part of the dissipative bracket is proposed. The bracket form of GENERIC generates the governing equations for large-strain thermo-viscoelasticity including the nonlinear evolution law for the internal variables associated with inelastic deformations. The GENERIC formalism facilitates the free choice of the thermodynamic variable. In particular, one may choose (i) the internal energy density, (ii) the entropy density, or (iii) the absolute temperature as the thermodynamic variable. A mixed finite element method is proposed for the discretization in space which preserves the GENERIC form of the resulting semi-discrete evolution equations. The GENERIC-consistent space discretization makes possible the design of structure-preserving time-stepping schemes. The mid-point type discretization in time yields three alternative schemes. Depending on the specific choice of the thermodynamic variable, these schemes are shown to be partially structure-preserving. In addition to that, it is shown that a slight modification of the mid-point type schemes yields fully structure-preserving schemes. In particular, three alternative energy-momentum-entropy consistent schemes are devised associated with the specific choice of the thermodynamic variable. Numerical investigations are presented which confirm the theoretical findings and shed light on the numerical stability of the newly developed schemes.

The rest of this article is organized as follows. In Section 2, the bracket form of the GENERIC formalism is introduced for thermo-viscoelasticity with heat conduction. In particular, in Section 3, a new material version of the inelastic dissipative bracket is proposed, which makes possible the free choice of the thermodynamic variable. In Section 4, a mixed finite element approach is proposed which yields a GENERIC-consistent semi-discrete form of the evolution equations. The mid-point type discretization in time is dealt with in Section 5 leading to alternative (partially) structure-preserving schemes. After the numerical investigations in Section 6, conclusions are drawn in Section 7.

GENERIC-BASED FORMULATION OF LARGE STRAIN THERMO-VISCOELASTICITY
The present work relies on a material description of large strain thermo-viscoelasticity within the framework of GENERIC. The developments presented herein extend our previous work 39 on thermo-elasticity to the realm of inelastic deformations. For more background on the GENERIC framework the reader is referred to Betsch and Schiebl 39 and the references therein. In the GENERIC framework, the time-evolution of any functional  is governed by The evolution equation (1) is valid for closed systems in which the boundaries are disregarded. We consider a continuum body with material points X = X A e A in the reference configuration  ⊂ R 3 ( Figure 1). The material coordinates X A refer to canonical base vectors e A ∈ R 3 . Here and in the sequel, the summation convention applies to repeated indices. Within the Lagrangian description of continuum mechanics the deformed configuration of the body at time t is characterized by the deformation map ∶  ×   → R 3 , where  = [0, T] is the time interval of interest. Accordingly, the placement at time t of the material point X ∈  is given by x = (X, t). The corresponding velocity field v ∶  ×   → R 3 is defined by v =̇, where a superposed dot denotes the material time derivative. The deformation gradient is given by or, in components, (F) i⋅ ⋅A = x i ∕ X A . Furthermore, the right Cauchy-Green tensor reads In the GENERIC (1), we consider functionals of the form where p = v is the linear momentum density and ∶   → R + is the mass density in the reference configuration. Moreover, ∶  ×   → R is a generalized thermodynamic field which may be chosen from among three alternatives, F I G U R E 1 Reference configuration  and deformed configuration (, t) at time t. For now the focus is on the motion of isolated thermomechanically coupled solids. That is, external tractions acting on the boundary as well heat fluxes across the boundary are disregarded ∈ { , , u}. The three alternative fields are the absolute temperature , the entropy density , and the internal energy density u.
To describe inelastic deformations, we make use of the internal variable C −1 p ∶  ×   → R 3×3 . Since we focus on isotropic inelastic deformations, C −1 p is assumed to be symmetric, that is, (C −1 p ) T = C −1 p . Therefore the density function a in (4) is subject to the same isotropic restrictions as the Helmholtz free energy function, to be introduced in Section 3.
GENERIC (1) decouples reversible and irreversible processes. In particular, the Poisson bracket { ⋅ , ⋅ } is responsible for reversible processes, while the dissipative (or friction) bracket [⋅ , ⋅] embodies irreversible processes. Accordingly, these two brackets constitute fundamental entities of the GENERIC framework and shall be specified next. Since inelastic deformations are purely irreversible, the Poisson bracket retains its thermo-elastic form (see Betsch and Schiebl 39 ) where  and  are arbitrary functionals of the form (4). Furthermore, the functional derivatives in (5) are given by a = a − Div F a, p a = p a, a = a, In (5), denotes the entropy density giving rise to the total entropy of the system The dissipative bracket for the thermo-viscoelastic problem at hand can be additively decomposed into a part due to heat conduction and a part taking into account inelastic deformations: In analogy to thermo-elasticity with heat conduction, the dissipative bracket accounting for heat conduction is given by (see Betsch and Schiebl 39 ) Here, K = K(C, ) is the material conductivity tensor, for which we assume that K T = K and a : K : a ≥ 0 holds for all a ∈ R 3 . An important element of the GENERIC framework is that the absolute temperature takes the form In (9) and (10), u denotes the internal energy density which, together with the kinetic energy and the potential of dead loads, constitutes the total energy of the system Here, b ∶   → R 3 represent prescribed body forces which, for simplicity, are assumed to be dead loads.
Concerning the contribution of inelastic deformations, we introduce the following form of the dissipative bracket (see Section 3 for further details): The fourth-order inelastic material flow tensor  has the properties  T =  (major symmetry) and M ∶  ∶ M ≥ 0 (positive semi-definiteness) for all M given in (17). It can be easily verified that the Poisson bracket (5) is skew-symmetric ({, } = −{, }), the dissipative bracket (12) is symmetric ([, ] = [, ]) and further satisfies the non-negativity condition [, ] ≥ 0. Moreover, the two brackets satisfy the fundamental degeneracy (or non-interaction) conditions {, } = 0 and [, ] = 0. In conjunction with these properties of the two brackets, GENERIC (1) ensures for closed systems (i) the conservation of total energy (d∕dt = 0), and (ii) a non-decreasing total entropy (d∕dt ≥ 0). Due to these structural properties, the GENERIC-based formulation offers an ideal starting point for the design of structure-preserving numerical methods.

Local evolution equations
We provide a short outline of the local evolution equations emanating from GENERIC (1). With regard to (1) the total energy  generates the reversible contribution through the Poisson bracket {, }. Using expression (11) for the total energy along with formulas (6) for the variational derivatives, Poisson bracket (5) yields where the first Piola-Kirchhoff stress tensor P = P(F, , C −1 p ) takes the form The irreversible contribution to GENERIC (1) is generated by the total entropy  through the dissipative bracket [, ]. Inserting expression (7) for the total entropy into the dissipative bracket (8) leads to where the material heat flux vector Q = Q(F, , C −1 p ) is given by and M = M(F, , C −1 p ) denotes the material representation of the Mandel stress tensor, which takes the form It is important to realize that in (14), (16) and (17), Θ is the temperature field which has been introduced in (10). These relationships are representative for the GENERIC-based formulation.
With regard to the left-hand side of GENERIC (1) and expression (4) for functional , the time derivative of  can be written as In the last equationĊ −1 p stands for C −1 p ∕ t. Substituting (18), (13), and (15) into GENERIC (1), we arrive at This equation has to hold for arbitrary functionals  of the form (4). Standard arguments now lead to the local evolution equations emanating from GENERIC (1): We stress again that in the above formulas the first Piola-Kirchhoff stress tensor P, the Mandel stress tensor M and the heat flux vector Q are given by formulas (14), (17), and (16), respectively.

GENERIC-based weak form of the IBVP
Starting from the GENERIC-based local evolution equations (20), we deduce the weak form of the initial boundary value problem (IBVP) pertaining to finite-strain thermo-viscoelasticity. To this end, we decompose the boundary  of the continuum body into a displacement boundary , on which = , and a traction boundary , on which PN = t, where and t are prescribed functions for t ≥ 0 ( Figure 2). Moreover,  ∪  =  and  ∩  = ∅. Similarly, for the thermal part, we consider the subsets of the boundary  and q , with the properties  ∪ q  =  and  ∩ q  = ∅ ( Figure 3). The thermodynamic variable is prescribed on , that is, = , whereas the heat flux is prescribed on q , that is, Q ⋅ N = q. Both and q are assumed to be given for t ≥ 0. We aim at the determination of the motion The unknown fields are subject to initial conditions of the form (⋅ , 0) = X, p(⋅, 0) = V 0 , (⋅, 0) = ini and C −1 p (⋅, 0) = I in . Here, V 0 is a prescribed material velocity field, and ini is a prescribed field of the thermodynamic variable ∈ { , , u}. The unknown fields are determined by applying a space-time discretization to the weak form of the IBVP at hand.
The weak form of the IBVP can be obtained by scalar multiplying the local evolution Eqs. 20 by suitable test functions and subsequently integrating over domain . The standard procedure yields F I G U R E 2 Mechanical part of the IBVP. Note that t = PN denotes prescribed external Piola tractions acting on the current boundary expressed per unit area of the reference boundary  F I G U R E 3 Thermal part of the IBVP.
Note that q = Q ⋅ N is the prescribed rate of heat transfer across the current boundary expressed per unit area of the reference boundary q  where the velocity field v = −1 p has been introduced. In weak form (21), w , w p ∶   → R 3 , w ∶   → R and w C −1 p ∶   → R 3×3 are test functions that have to satisfy the boundary conditions w = 0 and w p = 0 on , and w = 0 on .

Balance laws
We verify the pertinent balance laws of the IBVP at hand. For that purpose it suffices to consider the pure Neumann problem (i.e.,  = q  = ). We start with the total linear momentum of the continuum body defined by L = ∫  p dV.
Due to the arbitrariness of ∈ R 3 , (22) coincides with the balance law for linear momentum. Note that the parentheses on the right-hand side of (22) contain the resultant external force exerted on the continuum body (see also Figure 2).
The total angular momentum relative to the origin of the inertial frame is defined by J = ∫  × p dV. Choosing (w , w p , w , w C −1 p ) = (p × , × , 0, 0), weak form (21) yields Note that the symmetry condition FP T = PF T has been employed. Since ∈ R 3 is arbitrary, (23) corresponds to the balance of angular momentum. The parentheses on the right-hand side of (23) contain the resultant external torque about the origin (see also Figure 2).
Next, we substitute (w , w p , w , w C −1 p ) = (−b, −1 p, u, C −1 p u) into weak form (21). A straightforward calculation leads to the balance law for total energy We choose (w , w p , w , w C −1 p ) = (0, −1 p, , C −1 p ) concerning the balance of entropy in weak form (21) to obtain Here, use has been made of formula (10) for the temperature along with expressions (16) for the material heat flux vector and (17) for the Mandel stress tensor, respectively. The above equation can be rewritten as where  inel is the local production of entropy due to inelastic deformations and  cond is the local production of entropy due to heat conduction. The last equation corresponds to the Clausius-Duhem form of the second law of thermodynamics (see, for example, Gonzalez and Stuart, 42 Sec. 5).

INELASTIC PART OF THE DISSIPATIVE BRACKET
The present model for large strain thermo-viscoelasticity relies on the introduction of the internal variable C −1 p ∶  ×   → R 3×3 , whose time-evolution accounts for local inelastic deformations. Following Reese and Govindjee, 43 C −1 p can be associated with the multiplicative decomposition of the deformation gradient 44 into an elastic part F e and an inelastic (or viscous) part F p . Decomposition (26) gives rise to the relationship The restriction to the isotropic case implies that the free energy takes the separable form (see, for example Reese and Govindjee 45 ) where Ψ eq is the equilibrium part and Ψ neq is the non-equilibrium part.
Since inelastic deformations are purely irreversible in nature, they lead to a contribution to the dissipative bracket in GENERIC (1), cf. (8). To derive the inelastic dissipative bracket (12), we resort to the dissipative bracket derived in Hütter and Svendsen 24 within a temperature-based framework for GENERIC. In particular, Hütter and Svendsen 24 consider functionals of the form with associated density function a( , F, p, , F p ). The inelastic dissipative bracket from Hütter and Svendsen 24 is given by where u(F, , F p ) is the internal energy density and N is a fourth-order inelastic flow tensor which has the properties N T = N (major symmetry) and A : N : A ≥ 0 for all A ∈ R 3×3 . In components, these properties read As has been shown in Hütter and Svendsen, 24 inelastic dissipative bracket (30) comes along with the flow rulė In the last equation, η(F, , F p ) is the entropy density of the temperature-based formulation.

Change of variables
We perform a change of variables to transform the inelastic bracket (30) to the present setting which is based on functionals of the form (4). In particular, to link the current density functions a( , F, p, , C −1 p ) to those in (29), we express the generalized thermodynamic variable ∈ { , , u} in terms of the temperature by inverting relation (10) to get Moreover, relation (27) implies Now, the two density functions under consideration can be connected through where relationships (32) and (33) are employed on the right-hand side of the last equation. A straightforward application of the chain rule to the last equation yields Note that a = a, and F p a = F p a. Furthermore, with regard to (10) and (32) we have = Θ(F, (F, , C −1 p ), C −1 p ), from which follows that We thus obtain Substituting from (35) into (34) yields a = a Θ , Note that a = a and F p a = F p a. Making use of (36), the terms in the parenthesis of (30) can be rewritten as Taking into account the relationship the inelastic dissipative bracket (30) can be recast in the form This bracket coincides with the one introduced in (12). Note that in the above formula, C −1 p a = C −1 p a. Moreover,  is the material form of the fourth-order flow tensor N in (30). In components, Accordingly, the material flow tensor inherits symmetry and positive semi-definiteness. That is, in components, for all M introduced in (17). As has been shown in Section 2.1, the above inelastic dissipative bracket gives rise to flow rules of the form (cf. (20) It is worth noting that this flow rule can be viewed as material version of the viscoelastic evolution equation derived in Reese and Govindjee 43,45 and applied in, for example, Budday et al. 46 In the last equation is an isotropic, positive definite fourth-order tensor. Using these relationships, flow rule (42) can be recast in the forṁ where the relation has been employed. Note that in the last equation, definition (17) of the Mandel stress tensor has been taken into account.
Comparing (43) with present flow rule (41), leads to the conclusion that the respective fourth-order flow tensors are related by We further remark that the present viscoelastic evolution equation (41) can also be brought into the forṁ This version of the viscoelastic evolution equation has been used in Krüger et al.; 29 see also Groß et al. 16 Thus we conclude that the newly proposed inelastic dissipative bracket (12) gives rise to evolution equations for the internal variables that have been previously developed in the context of finite deformation thermo-viscoelasticity.

DISCRETIZATION IN SPACE
Concerning the space discretization of the present GENERIC-based formulation, we essentially apply an isoparametric finite element approach (see, for example, Hughes 47 ). Our main goal is to achieve a GENERIC-consistent space discretization in the sense of Betsch and Schiebl. 40 In particular, a GENERIC-consistent space discretization ensures that the discrete formulation inherits the fundamental balance laws for both the energy and the entropy from the underlying continuous formulation (cf. Section 2.3).
We first restate the governing equations of the IBVP to be discretized in space and time. With regard to the GENERIC-based weak form (21), we consider the following set of equations: Here, the first equation represents the local form of the kinematic relationshiṗ= v, while the second equation is the viscoelastic evolution equation. The third and fourth equation serve the purpose to introduce the new fields u and . This procedure facilitates a mixed finite element approach which turns out to be crucial for a GENERIC-consistent discrete formulation. Due to the arbitrariness of the test functions w u and w , (44) 3 and (44) 4 impose the conditions u = u and = , respectively. Moreover, are the GENERIC-specific representations of the first Piola-Kirchhoff stress tensor and the material Mandel stress tensor previously introduced in (14) and (17), respectively. Similarly, the material heat flux vector Q has been introduced in (16). We emphasize again that the GENERIC-based formulation is based on expression (10) for the temperature field. In the present mixed formulation this implies The finite element method is based on finite-dimensional approximations of the following quantities and As before, the summation convention applies, where a = 1, … , N, and N denotes the total number of nodes in the finite element mesh. Moreover, N a ∶  → R are the nodal shape functions with associated nodal values q a , v a ∈ R 3 and a , (u ) a , ( ) a ∈ R. Analogous approximations are used for the test functions w u , w , w p , and w , denoted by w h u , w h , w h p , and w h . In what follows, we summarize the space-discrete version of (44). Nodal collocation of kinematic equation (44) for a = 1, … , N. Viscoelastic evolution equation (44) 2 is collocated at the integration points X g ∈  used for the numerical evaluation of the volume integrals in (44). Accordingly, for g = 1, … , G, where G denotes the total number of integration points. Here and in the sequel, index g indicates evaluation at the integration point X g ∈ . For example, where are the internal energy and entropy densities at point X g . Furthermore, the discrete deformation gradient at X g , F h g , and the discrete generalized thermal variable, h g , follow from (47) 1,3 and thus take the form Similarly, the discrete temperature at X g follows from (46) and is given by where, in view of interpolations (48), In the discrete setting, the fields u h and h are determined through (44) 3,4 . In particular, inserting the approximations (48) along with the corresponding formulas for a = 1, … , N. To calculate the spatial integrals, appropriate numerical quadrature formulas of the form have been applied to obtain (55). * Here, w g play the role of generalized weighting coefficients resulting from the specific quadrature rule along with the isoparametric description of reference domain . Now, (55) 1 can be solved for the nodal quantities (u ) a , a = 1, … , N. To this end, we introduce the components H ab of the positive definite Gram matrix [H ab ], so that (55) 1 can be rewritten as The where c a denotes the Kronecker delta. Now, interpolation formula (48) 1 along with (58) lead to the result which will be utilized below. Similarly, interpolation formula (48) 2 in conjunction with (55) 2 yield the result Next, we turn to the discretization of (44) 5 , which can be done in a straightforward way to obtain For simplicity, we have neglected the contribution of the external tractions which could be easily added to the above equation. In the above equation, first Piola-Kirchhoff stress tensor P g at point X g ∈  is given by Moreover, in (62), the components M ab of the mass matrix are given by We further introduce nodal momentum vectors p a conjugate to nodal position vectors q a through the standard relation Eventually, we consider the space-discrete version of (44) 6 . Straight-forward application of our approach yields where the material heat flux vector Q g at point X g ∈  is given by .
For simplicity, in the space-discrete evolution equation (66) for the generalized nodal thermal variable b , the term accounting for heat transfer across the boundary has been neglected. To summarize, the resulting evolution equations for the space-discrete system at hand can be written in the forṁ for 1, … , N and 1, … , G. In (67) 1 , M ab stands for the components of the inverse mass matrix satisfying M ab M bc = c a . The set of Equation (67) constitutes nonlinear first-order ordinary differential equations for the determination of the unknowns which can be collected in the state vector In the sequel, state vector (68) will be viewed as column vector. In particular, this implies that the six independent components (C −1 p ) AB g , g = 1, … , G, of the internal variable (C −1 p ) g (at quadrature point X g ) are arranged in a column vector. The set of evolution equations (67) fits into the GENERIC framework for discrete systems, as shown next.

GENERIC-consistent space discretization
Our discretization approach presented above is GENERIC-consistent in the sense of Betsch and Schiebl 40 and thus can be framed in the context of GENERIC for discrete systems. To see this, the set of evolution equation (67) needs to be put into the forṁz Here, the focus is again on closed systems in which the boundary contributions are disregarded. In analogy to (1), the time-evolution of the state vector is decomposed additively into a reversible part generated by the total energy  and an irreversible part generated by the total entropy . Poisson matrix L needs be skew-symmetric, while friction matrix M needs be symmetric positive semi-definite.
In the above equation, the total energy of the discrete system under consideration is given by This is the space-discrete version of total energy (11). Similarly, the space-discrete version of total entropy (7) takes the form To get the gradient of the total energy, ∇(z), we consider the derivative of (70) with respect to time. Accordingly, the left-hand side of (70) yields while the right-hand side of (70) gives Comparing (72) with (73) yields the following expressions for the respective derivatives of the total energy In particular, inserting from (74) 3 into (60), we obtain the important relationship Similarly, the derivatives of total entropy (71) take the form Inserting from (76) 3 into (61), we obtain Now, guided by discrete GENERIC (69), the evolution equations in (67) can be recast. In particular, taking into account (74) 2 , kinematic relation (67) 1 can be rewritten asq Next, evolution equation (67) 2 for the nodal momentum vectors together with (63) and (54) yieldṡ Taking into account (74) 1 and (75), we obtaiṅ Concerning the evolution of the nodal thermal variable a , (67) 3 along with (51) 2 , (54), and (77) lead tȯ Eventually, evolution equation (67) 4 for the internal variable, together with (51) 2 , (54), and (77) result in Next, we aim at introducing relation (76) 4 into the last term of (81). To this end, we introduce the fourth-order tensor  with components Note that  enjoys major symmetry, () ABIJ = () IJAB , due to the major symmetry of the material flow tensor  and the symmetry of C −1 p . Now, (81) can be recast in index form Altogether, evolution equation (67) pertaining to the state variables of the discrete system at hand can be recast in the formq where For simplicity of exposition, in (85), N a g = N a (X g ). Note that the properties m ab = m ba and m ABIJ g = m IJAB g hold. Now, evolution equations (84) give rise to specific forms of Poisson matrix L and friction matrix M in GENERIC (69). In particular, the Poisson matrix takes the form where I is the identity matrix (with appropriate dimension corresponding to the partitioning of the state vector (68)), and matrix [l a⋅ ⋅b ] consists of vectors l a⋅ ⋅b defined in (85) 1 . Specifically, we have It can be observed that Poisson matrix (86) is skew-symmetric. Furthermore, the friction matrix is given by Note that the last matrix is block-diagonal and symmetric. It can be easily observed that friction matrix (88) is symmetric and positive semi-definite.

Conservation properties
We next verify the conservation properties of the semi-discrete formulation of the closed system dealt with in Section 4.1.
Since the evolution equations pertaining to the semi-discrete formulation can be brought into GENERIC form (69), (i) conservation of total energy, and (ii) non-decreasing total entropy are automatically satisfied. To see this, we first consider the time-derivative of the total energy, where (69) It can be verified by a straight-forward calculation that the second non-interaction condition is satisfied. In addition to that, the positive semi-definiteness of friction matrix (88) implies the result d∕dt ≥ 0. In particular, this result represents the semi-discrete version of (25) (apart from the boundary term in (25) originating from heat flux across the boundary). Thus, the total entropy of the closed system at hand is non-decreasing, due to the irreversible nature of heat conduction and visco-elastic deformations.
Since the material response is assumed to be frame-indifferent (or objective), specific symmetry properties are inherent to the discrete system under consideration. In particular, invariance under rigid rotations implies satisfaction of the following conditions (see Appendix A for further details): for all b = 1, … , N. We tacitly assume that no resultant external torque is acting on the system (i.e., the right-hand side of (23) is assumed to vanish). The semi-discrete version of the angular momentum relative to the origin is given by where formulas (47) 1,2 along with definition (65) of the nodal momentum vectors p a have been used. The time-derivative of (94) reads where (67)

Choice of the thermodynamic variable
We shortly outline the impact of the specific choice of the thermodynamic variable, ∈ { , , u}, on the structure of GENERIC (69). For simplicity, in this section we neglect body forces, that is, b = 0 and still focus on closed systems. Choosing the internal energy density as thermodynamic variable, that is, = u, the total energy (70) takes a particularly simple form given by As before, u h g stands for u h (X g , t). In particular, interpolation formula (53) 2 gives rise to u h g = N a (X g )u a (t). The derivatives of the total energy in (74) simplify to q a = 0, Consequently, Moreover, for the choice = u friction matrix (88) attains a particularly simple block-diagonal form, since m AB g,b = 0, and coefficients m ab only contain contributions due to heat conduction (cf. (85)).
Choosing the internal entropy density as thermodynamic variable, that is, = , yields a particularly simple form of the total entropy (71) given by where interpolation formula (53) 2 gives rise to h g = N a (X g ) a (t). Consequently, the derivatives of the total entropy in (76) simplify to q a = 0, Thus d dt (z) = ∇(z)̇z Moreover, the choice = leads to l a b = 0 (see (85)), so that Poisson matrix (86) yields a particularly simple form.
In contrast to the above considerations, selecting the total temperature as thermodynamic variable, that is, = , essentially does not lead to any simplifications. We eventually remark that these conclusions also affect the discretization in time, which will be treated next.

DISCRETIZATION IN TIME
We now turn to the discretization in time of the semi-discrete GENERIC-consistent evolution equations derived in Section 4.1. To this end, we focus on a representative time interval [t n , t n + 1 ] with corresponding time-step size Δt = t n+1 − t n . We aim at second-order accurate, implicit time-stepping schemes based on the mid-point rule. Application of the mid-point rule to (69) yields Here, z n stands for the discrete vector of state variables at time t n , and Provided that the state variables z n are given, the state variables z n + 1 can be determined by solving (102).

Partially structure-preserving schemes
Next, we check whether, or under what conditions, structure-preserving properties hold in the discrete setting. In this connection, we shall see that the specific choice of the thermodynamic variable ∈ { , , u} plays a crucial role. It can be easily verified that non-interaction conditions (90) and (92)  where (102) has been used. In analogy to the time-continuous case the right-hand side of the above equation vanishes due to the skew-symmetry of L(z n+ 1 2 ) and non-interaction condition (103) 1

. Thus
On the other side, in general. This inequality complies with the well-known fact that the mid-point rule is not capable to conserve nonlinear first integrals in general. However, there exists the exceptional case related to the choice = u, for which (104) turns into an equality. This is due to the fact that for = u the total energy takes the form (96) and thus(z) is merely quadratic. Since the mid-point rule preserves quadratic first integrals (see Leimkuhler,1 Sec. 4.4.2), the choice = u yields a structure-preserving scheme which is capable to conserve total energy.
Concerning the evolution of total entropy, guided by (91), we consider where again (102) has been used. Employing non-interaction condition (103) 2 , we obtain where the positive semi-definiteness of friction matrix M(z n+ 1 2 ) has been taken as a basis. On the other hand, in analogy to (104), we have This implies that, despite the encouraging result (105), the mid-point scheme in general does not guarantee a non-decreasing entropy. However, there again is an exception related to the choice = . For this particular case, the total entropy takes the form (99) and thus(z) is merely linear. Accordingly, the choice = turns inequality (106) into an equality and the entropy-based mid-point scheme is therefore capable to correctly reproduce the second law of thermodynamics in the discrete setting.
We eventually verify that all mid-point-based schemes under consideration are capable to conserve angular momentum. The incremental change of angular momentum (94) can be written in the form Mid-point scheme (102) gives rise to q a n+1 − q a n = ΔtM ab p b ) . (108) Inserting from (108) into (107) yields Symmetry conditions (93) imply Inserting from (110) into (109) and taking into account the symmetry of M ab together with the skew-symmetry of the vector product leads to the result J h n+1 = J h n . To summarize, depending on the choice for the thermodynamic variable we get three alternative mid-point schemes which are partially structure-preserving. Correspondingly, the resulting schemes are denoted by (EM) u (EM scheme related to = u), (ME) (momentum-entropy scheme associated with = ), and M (momentum scheme related to = ).

Fully structure-preserving schemes
In this section, we show that one specific modification of the three alternative mid-point-based schemes considered above turns all of them into EME schemes. That is, independent of the choice of ∈ { , , u}, all scheme are (i) thermodynamically consistent in the sense that they obey discrete versions of the two fundamental laws of thermodynamics, and (ii) capable to conserve angular momentum.
To reach this goal, the aforementioned modification of the mid-point integrator should (i) maintain the structure-preserving properties (103), and (ii) ensure that inequalities (104) and (106) are turned into equalities for general nonlinear functions (z) and (z). For that purpose, we resort to the notion of discrete derivative introduced in Gonzalez. 41 In particular, we replace the mid-point derivatives ∇(z n+ 1 2 ) and ∇(z n+ 1 2 ) by discrete derivatives d(z n , z n+1 ) and d(z n , z n+1 ), respectively. Accordingly, starting with the discrete derivatives of the total energy, based on (74), we introduce the discrete derivatives where the discrete derivatives of the internal energy density function u(F, , C −1 p ) are denoted by du∕d(C −1 p ) AB , d u, and respectively. In the last equation, the frame-indifferent representation of the internal energy density has been accounted for (see Appendix A). We refer to Appendix B for the specific definitions of the discrete derivatives d C u, d u, and du∕d(C −1 p ) AB , respectively. Similarly, the discrete derivatives of the total entropy rely on the application of the discrete derivative to the internal entropy density. That is, based on (76), we introduce where the discrete derivatives of the internal entropy density function (F, , C −1 p ) are denoted by d ∕d(C −1 p ) AB , d , and respectively (see Appendix B for further details). In addition to (111) and (113), the derivatives of u and contained in Poisson matrix (86) and friction matrix (88) need to be replaced by the corresponding discrete derivatives. The thus obtained discrete versions of the Poisson and friction matrix are denoted by L(z n , z n + 1 ) and M(z n , z n + 1 ), respectively. These modifications to mid-point integrator (102) yield fully structure-preserving EME schemes of the form z n+1 − z n = ΔtL(z n , z n+1 )d(z n , z n+1 ) + ΔtM(z n , z n+1 )d(z n , z n+1 ).
It is important to note that the above-described modifications of the mid-point rule retain the crucial non-interaction and symmetry conditions. In particular, non-interaction conditions (103) are now replaced by M(z n , z n+1 )d(z n , z n+1 ) = 0, while symmetry conditions (110) are replaced by 0 = q a n+ 1 2 × d q a (z n , z n+1 ), Note that l a b has been introduced in (85) 1 . Accordingly, following the present procedure, the discrete derivatives d F and d are used in (85) 1 , to get l a b (z n , z n+1 ). The so-called directionality property of the (partitioned) discrete derivatives (cf. Gonzalez 41 ) ensures that inequalities (104) and (106) are now replaced by the equalities Note that directionality properties (118) are also verified in Appendix B.
The new schemes are indeed fully structure-preserving, independent of the choice for ∈ { , , u}. This can be shown in a straightforward manner by following the steps in Section 5.1. For obvious reasons, the new EME consistent schemes are abbreviated with (EME) , (EME) , and (EME) u .

NUMERICAL INVESTIGATIONS
In this section, the alternative mid-point type schemes newly developed in the present work are applied to representative numerical examples dealing with finite strain thermo-viscoelastodynamics. Depending on the specific choice for the thermodynamic variable ∈ { , , u}, the following methods are applied: Section 5.2 (EME) u (EME) (EME) In the numerical investigations, we shall focus on momentum maps associated with symmetries of the mechanical system at hand, and the balance laws associated with the two fundamental laws of thermodynamics. In this connection, we also consider the functional According to Gurtin,48 for certain types of environments,  is a natural Lyapunov function and thus qualifies as estimate for the numerical stability of the schemes under consideration.
In each time step, the schemes emanating from (102) and (115) generate a system of nonlinear algebraic equations for the determination of the state variables z n + 1 . To this end, we apply a Multilevel-Newton algorithm; † see Hartmann 49 and references therein for more details.

Material model
In order to particularize the Helmholtz free energy density (28) used in the numerical examples, we start from a temperature-based description. In particular, we consider where Here, J = √ det C is the determinant of the deformation gradient, J e = √ det(CC −1 p ) and , e , , and e are prescribed parameters, c > 0 is the specific heat capacity at constant deformation, is the coefficient of thermal expansion, and 0 is the reference temperature. We refer to Betsch and Schiebl 39 and the references therein for a detailed investigation of the thermoelastic part of the specific Helmholtz free energy density (120). Further, for the viscoelastic part, we refer to Groß 15 and the references therein. For simplicity, we assume incompressible material behavior. Quasi-incompressible material models for finite strain thermo-viscoelasticy are considered in, for example, References 50 and 51. It is now a straightforward exercise to calculate further quantities emanating from (120), depending on the specific choice for the thermodynamic variable ∈ {u, , } (see also Betsch and Schiebl 39 for additional details). In particular, the temperature-based formulation yields The formulation based on the entropy density gives Moreover, the formulation based on the internal energy density leads tô Concerning the constitutive equation for the material heat flux vector, we assume thermally isotropic material, with material conductivity tensor given by Here, k is a prescribed coefficient of thermal conductivity and, as before, J = √ det(C). Finally, the constitutive equation for the isotropic fourth-order material inelastic flow tensor is given by (see Groß and Betsch 12 or Reese and Govindjee 43 and references therein for the spatial representation of the isotropic fourth-order inelastic flow tensor)

Flying L-shaped block
The first numerical example deals with the L-shaped block depicted in Figure 4. The spatial discretization of the block relies on 117 tri-linear finite elements leading to 224 nodes. The initial temperature field is varying linearly over the height (x 3 direction) of the block. In particular, at x 3 = 0, the initial temperature is prescribed as a , while at x 3 = h, the temperature is prescribed as b . The whole block is assumed to be thermally insulated (q = 0 on q  = ). Starting at rest, Piola traction vectors t a and t b are acting on two parts of the boundary surface of the block (Figure 4). The external loads are applied in the form of a hat function over time. In particular, the traction vectors are given by Table 1 provides a summary of the data used in the simulations. During the loading phase (t ≤ 4s) the time step size for all simulations is Δt = 0.05, such that all systems start from the same energy-and entropy level directly after vanishing external loads. No Dirichlet boundary conditions are applied. Since in the initial loading phase the external forces are equilibrated, the total linear momentum of the block is a conserved quantity. In addition to that, after the loading phase (t ≥ 4s) no external torque is acting on the block.
Consequently, the total angular momentum is conserved as well. All of the integrators under consideration are capable to exactly conserve both momentum maps (up to numerical round-off), independent of the chosen time step. This can be observed from Figure 5, where representative numerical results of the EM integrator (EM) u are shown. After the loading phase, the total energy must be a conserved quantity. As expected, the (EM) u scheme does exactly reproduce this conservation law (up to numerical round-off); see Figure 6. In contrary, the schemes (M) and (ME) are not capable of conserving the total energy for larger time step sizes. Depending on the time step size, both schemes tend to increase the energy which can be observed from Figure 7. Typically, such energy blow-ups eventually lead to a failure of the iterative (Newton-Raphson) solution procedure. In the diagrams, the break down of the simulation is indicated by vertical lines.
Regardless of the capability of the (EM) u scheme to conserve the total energy, the simulation still breaks down at about the same point in time as the break down of (M) and (ME) occurs. The numerical instability of (EM) u is accompanied by a nonphysical decrease of the total entropy as can be observed from Figure 8. Although not as pronounced as (EM) u , (M) occasionally yields an incremental decrease of the total entropy ( Figure 8). However, the total entropy ought to be a non-decreasing function of time. In contrast to (EM) u and (M) , (ME) is capable to correctly adhere to the second law of thermodynamics, independent of the time step ( Figure 9). Indeed, in each time step, the total entropy does increase, as can be observed from Figure 9. The incremental change of total entropy of the closed system at hand is related to heat conduction and inelastic deformations. For the mid-point rule, it can be directly calculated from (105). A straightforward calculation taking into account friction matrix (88) yields where Note that (125) can be viewed as discrete version of (25). Accordingly, for sufficiently small time steps, the incremental change of total entropy can be calculated for the mid-point schemes from (125). Similarly, for the EME schemes, the incremental change of total entropy can be calculated from (118) 2 . The two contributions to the discrete evolution of the incremental change of total entropy are visualized in Figures 10 and 11.
After the loading phase, the environment of the present example can be characterized as thermally perfect in the sense of Gurtin. 48 Thus  defined in (119) plays the role of a Lyapunov function that has to decrease with time ( Figure 12). However, the partially structure-preserving schemes (EM) u , (ME) , and (M) do not correctly reproduce this behavior, as can be seen from Figures 13 and 14. That is, depending on the time step and the duration of the simulation, all of the schemes inevitably exhibit numerical instabilities characterized by increasing values of .
Only the EME schemes developed in Section 5.2 completely reproduce the required behavior and thus can prevent numerical instabilities. This can be concluded from Figures 12 and 14, where the results of (EME) are shown. The corresponding results of (EME) and (EME) u are practically indistinguishable from those of (EME) .
Eventually, the motion of the L-shaped block is illustrated in Figure 15 with snapshots at successive points in time. In addition to that, the distribution of the temperature over the block is shown.

Rotating disk
The second example deals with a rotating disk subjected to prescribed heat flow over part of the boundary surface ( Figure 16). The spatial discretization of the disk is based on 200 tri-linear finite elements leading to a total of 360 nodes. The initial velocity distribution over the disk results from a prescribed angular velocity 0 ∈ R 3 and is given by The initial temperature of the disk is homogeneously distributed and equal to the reference temperature 0 . In an initial period of time, t ∈ [0, 4]s, heat flow is prescribed over one quarter of the lateral boundary surface ( Figure 16). In particular, the heat flow into the disk is described by A plot of function f (t) can be found in Figure 16. The rest of the boundary surface of the disk is assumed to be thermally insulated (q = 0). Note that the prescribed heat flow vanishes after t = 4s. For t > 4s, the complete boundary surface of the disk is assumed to be thermally insulated (q = 0 on q  = ). Then the environment of the disk is thermally perfect in the sense of Gurtin. 48 A summary of the data used in the simulations of the rotating disk can be found in Table 2.
During the loading phase (t ≤ 4s) the time step size for all simulations is Δt = 0.04, such that all systems start from the same energy and entropy level directly after vanishing external loads. Due to the fact that neither external loads act on the disk, nor displacement boundary conditions are imposed, the mechanical system at hand has translational and rotational symmetry. Consequently, the corresponding momentum maps are first integrals of the motion. All integrators under consideration are capable to conserve the respective momentum map. Representative numerical results are shown in Figure 17. Since heat flow into the system is prescribed in the initial time period [0, 4]s, the total energy is expected to increase. For t > 4s, the system is closed and the total energy should stay constant. Again the (EM) u scheme is capable to correctly reproduce the first law ( Figure 18).
However, despite the algorithmic energy conservation (for t > 4s), the (EM) u scheme is not devoid of numerical instabilities, depending on the time step. The corresponding point in time of the break down of the simulation is indicated with a vertical line in the diagrams. At about the same points in time, (ME) and (M) break down as well ( Figure 19). For these schemes the break down is accompanied by a sudden increase of the total energy leading to the divergence of the Newton-Raphson iterations. For the considered duration of the simulation (t ∈ [0, 30]s), a time step of Δt = 0.04s is small enough to retain numerical stability of the three partially structure-preserving schemes at hand. The EME schemes developed in Section 5.2 do not lead to any numerical instabilities even for the large time step Δt = 0.1s. In particular, all of the EME schemes are capable to conserve the total energy for t > 0.4s. This is shown exemplarily for (EME) in Figure 24.
Due to the prescribed heat flow into the disk, the total entropy of the disk is expected to increase in the initial time period [0, 4]s. For t > 4s, the system is closed and the total entropy should be a non-decreasing function of time. The (EM) u scheme does not correctly reproduce the second law as can be seen from Figure 20. Accordingly, the divergence of the iterative solution procedure is accompanied by a nonphysical decrease of the total entropy. The (M) closely adheres to the second law as can be observed from Figure 20. As expected, the (ME) scheme is capable to exactly satisfy the second law, independent of the time step. This can be observed from Figure 21. In particular, Figure 21(right) confirms that the change per time step of the total entropy is always positive. As expected all of the EME schemes are capable to correctly reproduce the non-decreasing evolution of the total entropy. This is shown exemplarily for (EME) in Figure 24.
In addition, the two contributions to the discrete evolution of the incremental change of total entropy, see (125), are visualized in Figures 22 and 23 for the mid-point-based schemes. Most of the contribution is due to conduction of heat, only about 5.9% is due to inelastic deformations in the present example ( Figure 24).
To shed further light on the numerical stability of the present schemes, we consider the Lyapunov function  defined in (119). For t > 4s the system is closed and the function  should decrease with time. As expected, the partially structure-preserving schemes (EM) u , (ME) , and (M) in general do not correctly reproduce this behavior, depending on  (Figures 25 and 26). In particular, it can be seen that the smallest time step, Δt = 0.04s, yields a stable numerical simulation, at least in the considered time interval [0, 30]s. However, for larger time steps Δt = 0.08s and Δt = 0.1s, the numerical instability of each scheme becomes visible through the increase of. Again the fully structure-preserving EME schemes developed in Section 5.2 provide enhanced numerical stability. All of them lead to a steadily decreasing , independent of the time step size. This is shown exemplarily for EME in Figure 26.
Eventually, the motion of the disk is illustrated in Figure 27 with snapshots at successive points in time. In addition to that, the distribution of the temperature over the disk is shown.

CONCLUSIONS
Starting from a GENERIC-based formulation of large-strain thermo-viscoelasticity, we have developed alternative structure-preserving schemes based on the implicit mid-point rule. Depending on the choice of the thermodynamic variable, already the plain mid-point rule yields partially structure-preserving schemes (see Section 5.1). Despite this, these schemes cannot prevent from numerical instabilities as has been shown in the numerical examples. Specifically, choosing the internal energy density as thermodynamic variable leads to an EM scheme (called (EM) u ). For Hamiltonian mechanical systems, EM schemes typically provide enhanced numerical stability. However, this situation does no extend to thermomechanically coupled dissipative systems. Here, EME schemes such as those developed in Section 5.2 are required to prevent numerical instabilities. The three EME schemes newly proposed in the present work essentially rely on the following features: • The underlying GENERIC formulation in bracket form (Section 2) leads to characteristic expressions such as (14) and (17) for the first Piola-Kirchhoff stress tensor and the Mandel stress tensor, respectively.
• The transformation properties of the GENERIC description facilitate the use of alternative thermodynamic variables such as the absolute temperature, the internal energy density, and the entropy density used in the present work.
• The newly proposed material form of the inelastic dissipative bracket (Section 3) makes possible to include into the GENERIC formulation often used nonlinear evolution laws for the internal variables associated with inelastic deformations.
• The mixed finite element discretization in space (Section 4) has been shown to be GENERIC-consistent. This means that the evolution equations for the state variables of the semi-discrete system fit into the GENERIC framework for discrete systems. This is an essential prerequisite for the development of fully structure-preserving EME schemes.

ACKNOWLEDGEMENT
This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) 388118188 (BE 2285/13-1). This support is gratefully acknowledged. Open Access funding enabled and organized by Projekt DEAL.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.