Thermoviscoelasticity in Kelvin-Voigt rheology at large strains

The frame-indifferent thermodynamically-consistent model of thermoviscoelasticity at large strain is formulated in the reference configuration with using the concept of the second-grade nonsimple materials. We focus on physically correct viscous stresses that are frame indifferent under time-dependent rotations. Also elastic stresses are frame indifferent under rotations and respect positivity of the determinant of the deformation gradient. The heat transfer is governed by the Fourier law in the actual deformed configuration, which leads to a nontrivial description when pulled back into the reference configuration. Existence of weak solutions in the quasistatic setting, i.e.\ inertial forces are ignored, is shown by time discretization.


Introduction
For a long time, thermoviscoelasticity was considered as a quite difficult problem even at small strains, mainly because of the nonlinear coupling with the heat-transfer equation which has no obvious variational structure; hence special techniques had to be developed. It took about two decades after the pioneering work by C.M. Dafermos [Daf82] in one space dimension that first three-dimensional studies occurred (cf. e.g. [BlG00,BoB03,Rou09]). The basic new ingredient was the L 1 -theory for the nonlinear heat equation developed in [BD * 97,BoG89]. At large strains, in simple materials, the problem is still recognized to be very difficult even for the case of mere viscoelasticity without coupling with temperature, and only few results are available if the physically relevant frame-indifference is respected, as articulated by J.M. Ball [Bal77], see also [Bal02,Bal10]. In particular, localin-time existence [LeM13] or existence of measure-valued solutions [Dem00,DST01] are known for simple materials. Further examples in this direction are [Tve08] for a general three-dimensional theory, but not respecting frame indifference and the determinant constraints, or [MOS13] for a one-dimensional theory using the variation structure. While the static theory for large-strain elasticity developed rapidly after [Bal77], there are still only few result for time-dependent processes respecting frame indifference as well as the determinant constraint. The first cases were restricted to rate-independent processes, such as elastoplasticity (cf. [MaM09,MiR16]) or crack growth (cf. [DaL10], see [MiR15,Sec. 4.2] for a survey. Recently the case of viscoplasticity was treated in [MRS18]. The main features of the model discussed in this work can be summarized in brief as follows: the thermo-visco-elastic continuum is formulated at large strains in a reference configuration, i.e. the Lagrangian approach. The concepts of 2nd-grade nonsimple material is used, which gives higher regularity of the deformation. The heat transfer is modeled by the Fourier law in the actual deformed configuration, but transformed (pulled back) into the reference configuration for the analysis. Our model respects both static frame-indifference of the free energy and dynamic frame indifference for the dissipation potential. Moreover, the local non-selfpenetration is realized by imposing a blowup of the free energy if the determinant of the deformation gradient approaches 0 from above, however we do not enforce global non-selfpenetration. Also, we neglect inertial effect; cf. Remark 6.6 for more detailed discussion.
Let us highlight the important aspects of the presented model and their consequences: (α) The temperature-dependence of the free energy creates adiabatic effects involving the rate of the deformation gradient. To handle this, the Kelvin-Voigt-type viscosity is used to control the rate of the deformation gradient. In addition, we separate the purely mechanical part, cf. (2.15) below, which allows us to decouple the singularities of large-strain elasticity from the heat equation.
(β) The heat transfer itself (and also the viscosity from (α)) is clearly rate dependent and the technique of rate-independent processes supported by variationally efficient energetic-solution concept cannot be used (which also prevents us from excluding possible global selfpenetration).
(γ) The equations for the solid continuum need to be formulated and analyzed in the fixed reference configuration but transport processes (here only the heat transfer) happen rather in the actual configuration and the pull-back procedure needs the determinant of the deformation gradient to be well away from 0. To achieve this, we exploit the concept of 2nd-grade nonsimple materials together with the results of T.J. Healey and S. Krömer [HeK09], which allow us to show that the determinant for the deformation gradient is bounded away from 0, see Section 3.1.
(δ) The transport coefficients depend on the deformation gradient because of the reasons in point (γ). For this, measurability in time is needed and thus the concept of global quasistatic minimization of deformation (as in rate-independent systems [MiR15] or in viscoplasticity in [MRS18]) would not be satisfactory; therefore we rather control the time derivative of the deformation, which can be done either by inertia (which is neglected in our work) or by the Kelvin-Voigt-type viscosity from (α).
(ǫ) The viscosity from (α) must satisfy time-dependent frame indifference as explained in [Ant98], thus it is dependent on the rate of the right Cauchy-Green tensor rather than on the rate of the deformation gradient itself. However, the adiabatic heat sources/sinks involve terms where the rate of the deformation gradient occurs directly. To control the latter by the former, we exploit results of P. Neff [Nef02] in the extension by W. Pompe [Pom03] for generalized Korn's inequalities, see Section 3.2. Here, again the mentioned concept of 2nd-grade nonsimple materials is used to control determinant of the deformation gradient, see (γ).
As mentioned above, our model heavily relies on the strain-gradient theories to describe materials, referred as nonsimple, or also multipolar or complex. This concept has been introduced long time ago, cf. [Tou62] or also e.g. [FrG06,MiE68,Pod02,Šil85,TrA86,BaC11] and in the thermodynamical concept also [Bat76]. In the simplest scenario, which is also used here, the stored-energy density depends only on the strain F = ∇y and on the first gradient ∇F of the strain. This case is called 2nd-grade nonsimple material. Possible generalization using only certain parts of the 2nd in the spirit of [KPS19] still need to be explored.
The structure of the paper is as follows. In Section 2 we present the model in physical and mathematical terms. After the precise definition of our notion of solution, Theorem 2.2 presents the main existence result for global-in-time solutions for the large-strain thermoviscoelastic system, while Corollary 2.3 gives the corresponding existence result for viscoelasticity at large-strain and at constant temperature, which, to the knowledge of the authors, is also new. A related result for isothermal large-strain viscoelasticity is derived in [FrK18], but there the limit of small strains is treated.
In Section 4 we start the proof of the main result by introducing certain regularizations as well as a time-incremental approach that is particularly constructed in such a split (sometimes called staggered) way that the deformation is first updated at fixed temperature and then the temperature is updated, where in some terms the old and in others the new deformation is used. Another important step in the analysis is the usage of an energy-like variable w = w(∇y, θ) instead of temperature θ, which enables us to exploit the balance-law structure of the heat equation; cf. [Mie13,MiM18] for arguments for the preference of energy in favor of temperature. As an intermediate result Proposition 5.1 provides the existence of solutions (y ε , θ ε ) of the regularized problem.
In Section 6 we finally show that the limit ε k → 0 for (y ε k , θ ε k ) → (y, θ) can be controlled in such a way that (y, θ) are the desired solutions. We conclude with a few remarks concerning potential generalizations and further applications of the methods.

Modeling of thermoviscoelastic materials in the reference configuration
We will use the Lagrangian approach and formulate the model in the reference (fixed) domain Ω ⊂ R d being bounded with smooth boundary Γ . We assume d ≥ 2 although, of course, the rather trivial case d = 1 works too if p ≥ 2 is assumed additionally to p > d in (2.30) below. We will consider a fixed time horizon T > 0 and use the notation I := [0, T ], Q := I × Ω, and Σ := I × Γ . For readers' convenience, Table 1 summarizes the main nomenclature used throughout the paper. To introduce our model in a broader context, we may define the total free energy and y deformation, y(t, x) ∈ R d , θ absolute temperature, (·) . time derivative, ψ = φ + ϕ free energy, σ el = ∂ F ψ elastic stress, F ζ viscous stress, F = ∇y deformation gradient, G = ∇F = ∇ 2 y valued in R d×d×d , w heat part of internal energy, h el elastic hyperstress, c v = c v (F, θ) heat capacity, q heat flux, M = Φ el +H main mechanical energy, H hyperstress energy, Φ cpl coupling energy, Ψ = M + Φ cpl free energy, W thermal energy, E = M + W total energy, ζ potential of dissipative forces, ξ rate of dissipation (=heat production), K = K(θ) material heat conductivity, K = K(F, θ) pulled-back heat conductivity, C = F ⊤ F right Cauchy-Green tensor, κ heat-transfer coefficient on Γ , g : I×Ω → R d a time-dependent dead force, f : I×Γ N → R d a boundary traction, ℓ an external mechanical loading, Ω the reference domain, Γ the boundary of Ω, Γ = Γ D ∩ Γ N , Table 1. Summary of the basic notation used throughout the paper.
respectively. The mechanical evolution part can then be viewed as an abstract gradient flow D .
y R(y, . cf. also [Tve08,MOS13] for the isothermal case and [Mie11] for the general case. The sum of the conservative and the dissipative parts corresponds to the Kelvin-Voigt rheological model in the quasistatic variant (neglecting inertia). The notation " ∂ " is used for partial derivatives (here functional or later in Euclidean spaces), while (·) ′ will occasionally be used for functions of only one variable.
Writing (2.2) locally in the classical formulation, one arrives at the nonlinear parabolic 4th-order partial differential equation expressing quasistatic momentum equilibrium where the viscous stress is σ vi = σ vi (F, . F , θ) and the elastic stress is σ el = σ el (F, θ), while h el is a so-called hyperstress arising from the 2nd-grade nonsimple material concept, cf. e.g. [Pod02,Šil85,Tou62]. In view of the local potentials used in (2.2), we have where G ∈ R d×d×d is a placeholder for ∇F .
An important physical requirement is static and dynamic frame indifference. For the elastic stresses, static frame indifference means that σ el (RF, θ) = R σ el (F, θ) and h el (RG) = Rh el (G) (2.5a) for all R ∈ SO(d), F and G. For the viscous stresses ✿ , dynamic frame indifference means that for all smoothly time-varying R : t → R(t) ∈ SO(d), cf. [Ant98]. Note that R may depend on t but not on x ∈ Ω, since frame-indifference relates to superimposing time-dependent rigid-body motions.
In terms of the thermodynamic potentials ζ, ψ, and H , these frame indifferences read as for R, F and ∇F as above. These frame indifferences imply the existence of reduced potentialsψ,ζ, andĤ such that F . More specifically, denoting G = [G αij ] the placeholder for ∂ ∂x j F αi with F αi the placeholder for ∂ ∂x i y α , the exact meaning is The ansatz (2.7) also means that The simplest choice, which is adopted in this paper for avoiding unnecessary technicalities, is that the viscosity σ vi is linear in . C. This is the relevant modeling choice for non-activated dissipative processes with rather moderate rates (in contrast to activated processes like plasticity having nonsmooth potentials that are homogeneous of degree 1 in a small-rate approximation). This linear viscosity leads to a potential which is quadratic in .
While we will be able to handle general dependence on F , it will be a crucial restriction that Furthermore, the specific dissipation rate can be simply identified in terms ofζ as (2.10) For our choice (2.9), we simply have ξ(F, . C: . F , θ). In brief, the standard thermodynamical arguments start from the free energy density ψ and the definition of entropy via s = −∂ θ ψ (here H does play no role as it is chosen to be independent of θ) and the entropy equation with the dissipation rate ξ from (2.10) and the heat flux q. We further use the formula . F and the Fourier law formulated in the reference configuration which will be specified later in (2.24). Altogether, we arrive at the coupled system div σ vi (∇y, ∇ . θθ ψ(F, θ) and ξ from (2.10) (2.13b) on Q. We complete (2.13) by some boundary conditions. For simplicity, we only consider a mechanically fixed part Γ D time independent undeformed (i.e. identity) while the whole boundary is thermally exposed with a phenomenological heat-transfer coefficient κ ≥ 0: where n is the outward pointing normal vector, and θ ♭ is a given external temperature. Moreover, following [Bet86] the surface divergence "div S " in (2.14a) is defined as div S (·) = tr ∇ S (·) , where tr(·) denotes the trace and ∇ S denotes the surface gradient given by ∇ S v = (I − n⊗ n)∇v = ∇v − ∂v ∂ n n. See (2.29) for a short mathematical derivation of the boundary conditions (2.14a) and (2.14c), and [Ste15, for the mechanical interpretation in second-order materials.
y with w = w(F, θ) . (2.17) In particular, the purely mechanical stored energy ϕ does not occur in (2.16) and does not influence the heat production and transfer (2.17). The energetics of the system (2.13)-(2.14) can be best described by introducing additional energy functionals as follows: Using σ el = ∂ F ϕ + ∂ F φ and integrating in time leads to the relation (2.20) that will be very useful for obtaining a priori estimates in the following sections.
Next, we test the heat equation in its simplified form (2.17) together with the boundary conditions (2.14d) by the constant function 1 (i.e. we merely integrated over Ω) and add the result to (2.20). After major cancellations we obtain the total energy balance: (2.21) In particular, we see that the total energy is conserved up to the work induced by the external loadings or the flux of heat through the boundary. From the entropy equation (2.11), we can read the total entropy balance (the Clausius-Duhem inequality): (2.22) This articulates, in particular, the second law of thermodynamics that the total entropy in the isolated systems (i.e. here q = 0 on Γ ) is nondecreasing with time provided K = K(∇y, θ) is positive semidefinite and the dissipation rate is non-negative. It is certainly a very natural modeling choice that Fourier's law is formulated in the actual (also called the deformed) configuration in a simple form, namely the actual heat flux is given by with the heat-conductivity tensor K = K(x, θ) considered as a material parameter possibly dependent on x ∈ Ω. We transform (i.e. pulled-back) this Fourier law into the reference configuration via the heat flux q(x) = K(x)∇θ = K(∇y(x)) ⊤ ∇ z θ(y(x)) and q = (Cof F ⊤ ) q, because fluxes should be considered as (d−1)-forms. With (2.23) the usual transformation rule for 2nd-order contra-variant tensors yields the heat-conductivity tensor if det F > 0, whereas the case det F ≤ 0 is considered nonphysical, so K is then not defined. Here we used the standard shorthand notation F −⊤ = [F −1 ] ⊤ = [F ⊤ ] −1 and also the algebraic formula F −1 = (Cof F ⊤ )/ det F . In what follows, we omit explicit xdependence for notational simplicity. Let us emphasize that in our formulations ∇θ is not treated as a vector, but a contravariant 1-form. Starting from θ(x) = θ(y(x)) the chain-rule gives ∇(x) = ∇y(x) ⊤ ∇ Y θ(y(x)). It should be noted that (2.23) is rather formal argumentation, assuming injectivity of the deformation y and thus existence of y −1 , which is however not guaranteed in our model; anyhow, handling only local non-selfpenetration while ignoring possible global selfpenetration is our modeling approach often accepted in engineering, too. For the isotropic case K(θ) = κ(θ)I, relation (2.24) can also be written by using the right Cauchy-Green tensor C = F ⊤ F as K = det(F )κ(θ)C −1 , cf. e.g. [DSF10,Formula (67)] or [GoS93,Formula (3.19)] for the mass instead of the heat transport. In principle, K in (2.23) itself may also depend on C = F ⊤ F , which we omitted to emphasize that K in (2.24) will depend on F anyhow.
In what follows, we will use the (standard) notation for the Lebesgue L p -spaces and W k,p for Sobolev spaces whose k-th distributional derivatives are in L p -spaces and the abbreviation H k = W k,2 . The notation W 1,p D will indicate the closed subspace of W 1,p with zero traces on Γ D . Moreover, we will use the standard notation p ′ = p/(p−1). In the vectorial case, we will write L p (Ω; R n ) ∼ = L p (Ω) n and W 1,p (Ω; R n ) ∼ = W 1,p (Ω) n . Thus, for example, For the fixed time interval I = [0, T ], we denote by L p (I; X) the standard Bochner space of Bochner-measurable mappings I → X with X a Banach space. Also, W k,p (I; X) denotes the Banach space of mappings from L p (I; X) whose k-th distributional derivative in time is also in L p (I; X). The dual space to X will be denoted by X * . Moreover, C w (I; X) denotes the Banach space of weakly continuous functions I → X. The scalar product between vectors, matrices, or 3rd-order tensors will be denoted by " · ", " : ", or " . . . ", respectively. Finally, in what follows, K denotes a positive, possibly large constant.
We consider an initial-value problem, imposing the initial conditions y(0, ·) = y 0 and θ(0, ·) = θ 0 on Ω. (2.26) Having in mind the form (2.17) of the heat equation, we can now state the following definition for a weak solution: At first sight, it seems that (2.27a) is not suited to apply the test function z = .
y, which is the natural and necessary choice for deriving energy bounds. Obviously, we will not be able to obtain enough control on ∇ 2 . y. However, using the abstract chain rules provides in Section 3.3 this problem can be handled by extending H(y) = Ω H (∇ 2 y) dx to a lower semicontinuous and convex functional on H 1 (Ω; R d ) by setting it ∞ outside W 2,p (Ω; R d ), see the rigorous proof of (5.9) in Step 3 of the proof of Proposition 5.1.
It will be somewhat technical to see that the weak formulation (2.27a) is indeed selective enough, in the sense that for sufficiently smooth solutions one can indeed obtain the classical formulation (2.13) together with the boundary conditions (2.14), cf. also [Rou13, Sect. 2.4.4]. In particular, abbreviating σ = σ vi (∇y, ∇ . y, θ) + σ el (∇y, θ), integrating by part once, and using the boundary conditions (2.14a,c) yields . . . (∇z⊗ n) dS dt. (2.28) We now want to show how the strong form (2.13a) and the associated boundary conditions (2.14a,c) follow from (2.28). For this goal, we apply Green's formula in the opposite direction to remove ∇ in front of the test function z. Using also the orthogonal decomposition of ∇z = ∇ S z + ∂ ∂ n z ⊗ n involving the surface gradient ∇ S z and writing shortly h for h el (∇ 2 y) ∈ R d×d×d , relation (2.28) leads to the identity Using the surface divergence div S and the projection P S : A → A−A n⊗ n to the tangential part, we obtain the integration by parts formula (cf. [Bet86] or [Ste15, pp. 358-359]) where the surface Γ is now assumed to be sufficiently smooth. Using this with A = h n for the previous relation we find where we have used z = 0 on Σ D = Σ \ Σ N . Now, taking z's with a compact support in Q, we obtain the equilibrium (2.13a) in the bulk. Next taking taking z's with zero traces on Σ but general ∂z ∂ n , we obtain (2.14c). Note that the latter condition implies P S (h n) = h n − h : ( n⊗ n) ⊗ n = h n. Hence, taking finally general z's, we obtain (2.14a), as P S can be dropped because of (2.14c).
Moreover, also note that, from the integral identity (2.27b), one can read w(∇y(0), θ(0)) = w(∇y 0 , θ 0 ) from which θ(0) = θ 0 follows when taken the invertibility of w(F, ·) and y(0) = y 0 into account. Now we exploit the decomposition (2.15) of ψ into φ and ϕ, which allows us to impose coercivity assumptions for the purely elastic part φ that are independent of those for ϕ, namely where GL + (d) denotes the set of matrices in R d×d with positive determinant. The last assumption in (2.30c) means that c v together with c −1 v are bounded, which is a major restriction. However, it allows for a rather simple estimation in Lemma 6.3; for alternative, more general situations dealing with increasing c v (·) we refer to [KrR19, Sec. 8.3].
The function w = w(F, θ) defined in (2.16) satisfies w(F, 0) = 0 by (2.15). Moreover, we have ∂ θ w(F, θ) = −θ∂ 2 θ φ(F, θ). Hence assumption (2.30c) implies, for all F ∈ GL + (R d ), the two-sided estimateŝ (2.31) The assumptions (2.30b,c) make the thermomechanical coupling through φ rather weak in order to allow for a simple handling of the mechanical part independently of the temperature. These restrictive assumptions are needed for our specific and simple way of approximation method rather than with the problem itself. E.g. the assumption in (2.30b) is used to facilitate the estimate (4.12), which allows us to control the difference between Ω (∇y k , θ) dx and Ω (∇y k−1 , θ) dx in terms of M(y k ), M(y k−1 ), and ∇y k −∇y k−1 2 L 2 . Moreover, after having derived uniform bounds on |∇y k | it will be exploited to show that the thermo-coupling stress ∂ F φ is bounded. Finally, (2.30d,h) makes the stored energy finite at time t = 0.
We now state our main existence results, which will be proved in the following Sections 4 to 6. The method will be constructive, avoiding non-constructive Schauder fixedpoint arguments, however some non-constructive attributes such as selections of converging subsequences will remain. More specifically, the proof is obtained by first making the a priori estimate for time-discretized solutions in, see Proposition 4.2, and then deriving an existence result for time-continuous solutions of an ε-regularized problem, see Proposition 5.1. Finally, Proposition 6.4 provides convergence for ε → 0.
As mentioned in the introduction, a lot of publications are devoted to the simpler isothermal viscoelasticity at largestrain, yet, in the multi-dimensional case, they do not satisfy all the necessary physical requirements. It is therefore worthwhile to present a version of our existence result by restricting it to this simpler case, for which a lot of assumptions are irrelevant or simplify. In particular, (2.15) simplifies as ψ(F, θ) = ϕ(F ). Of course, our theory only works because we are using a non-degenerate secondgrade material, where H(y) := Ω H (∇ 2 y) dx generates enough regularity to handle the geometric and physical nonlinearities. To the best of the authors knowledge, even the following result for isothermal viscoelasticity is new.
A similar regularization approach to isothermal large-strain viscoelasticity was considered in [FrK18], where the H(y) is multiplied with a small parameter that vanishes slower than the loading. Hence, the authors are able to show that their solutions are sufficiently close to the identity which allows them to exploit a simpler Korn's inequality obtained by a perturbation argument. Hence, to the best of the author's knowledge the following result is the first that allows for truly largestrains.
Before going into the proof of our main result, we show that our conditions are general enough for a series of nontrivial applications: Example 2.4 (Classical thermomechanical coupling). The classical example of a free energy in thermomechanical coupling is given in the form i.e. φ(F, θ) involves a term in the product form −a(θ)ϕ 1 (F ). For the purely mechanical part we may take the polyconvex energy ϕ(F ) = c 1 |F | s + c 2 /(det F ) q for det F > 0 and ∞ otherwise. For the thermomechanical coupling we obtain c v (F, θ) = −θ∂ 2 θθ ψ(F, θ) = c + a ′′ (θ)ϕ 1 (F ), thus to have positivity of the heat capacity c v , we assume a ′′ (θ) ≥ 0 and ϕ 1 (F ) ≥ 0. Moreover, we have Thus, we see that all assumptions in (2.30) can easily be satisfied, e.g. by choosing a(θ) = (1+θ) −α with α > 0, which is smooth bounded and convex, and taking any φ 1 ∈ C 2 c (R d×d ).
Example 2.5 (Phase transformation in shape-memory alloys). An interesting example of a free energy ψ occurs in modeling of austenite-martensite transformation in so-called shape-memory alloys: cf. e.g. [Rou04] and references therein. Here a denotes the volume fraction of the austenite versus martensite which is supposed to depend only on temperature. Of course, this is only a rather simplified model. For, ψ 0 (θ) = cθ(1− log θ) it complies with the ansatz The heat capacity then reads as To ensure its positivity, ψ 0 is to be strictly concave in such a way that ψ ′′ 0 (θ) ≤ K/θ and then inf (F,θ) θa ′′ (θ)[ϕ A −ϕ M ](F ) + K > 0 is to (and can) be ensured by suitable modeling assumptions.
Example 2.6 (Thermal expansion). Multiplicative decomposition F = F el F th with the "thermal strain" F th = I/µ(θ) and the elastic strain F el which enters the elastic part of the stored energy ϕ. This leads to (2.33) Unfortunately, (2.33) is inconsistent with the ansatz (2.15) because the contribution ϕ which has been important for our analysis due to uniform coercivity, cannot be identified in (2.33).

A few auxiliary results
In this subsection we provide a series of auxiliary results that are crucial to tackle the difficulties arising from large-strain theory. First we show how the theory developed by Healey and Krömer [HeK09] which allows us to show that a bound for the elastic energy M(y, θ) provides lower bounds on the det ∇y. This can then be used to establish the validity of the Euler-Lagrange equations and useful λ-convexity result, which is needed for obtaining optimal energy estimates. Second we provide a version of Korn's inequality from Pompe [Pom03] that allows us to obtain dissipation estimates via D(y, . . Finally, in Section 3.3 we provide abstract chain rules as derived in [MRS13, Sec. 2.2] that allows us to derive energy balances like (2.20) from the corresponding weak equations.

Local invertibility and Euler-Lagrange equations
A crucial point in large-strain theory is the blow-up of the energy density ψ(F, θ) for det F ց 0. Thus, it is desirable to find a suitable positive lower bound for det ∇y(t, x). The following theorem is an adaptation of the result in [HeK09, Thm. 3.1].
Theorem 3.1 (Positivity of determinant). Assume that the functional M : W 2,p (Ω; R d ) → R ∞ satisfies the assumption (2.30a) and (2.30d). Then, for each C M > 0 there exists a C HK > 0 such that all y ∈ Y id with M(y) ≤ C M satisfy Proof. We give the full proof, since our mixed boundary conditions are not covered in [HeK09]. From M(y) ≤ C M and the coercivities of ϕ and H we obtain det ∇y ≥ 0 a.e.
in Ω and the a priori bounds M . Together with the Dirichlet boundary conditions in Y id we obtain an a priori bound for y in W 2,p (Ω; R d ) and hence also in C 1,λ (Ω; R d ), where λ = 1 − d/p > 0. This proves the first two assertions.
In particular, the function δ : x → det(∇y(x)) is Hölder continuous as well with Since Ω is a bounded Lipschitz domain, there exist a radius r * > 0 and a constant α * > 0 such that for all x ∈ Ω the sets B r * (x) ∩ Ω contains an interior cone Thus, using the Hölder continuity we can estimate as follows: where in the last estimate we crucially used the assumption q > pd/(p−d) which implies λq > d. Since in the last expression both exponents of δ(x) are positive, we obtain the explicit lower bound which gives the third assertion in (3.1). The last assertion follows via the implicit function theorem.

A generalized Korn's inequality
The following result will be crucial to show that the nonlinear viscosity depending on F = ∇y really controls the H 1 norm of of the rate .
y. It relies on Neff's generalization [Nef02] of the Korn inequality, in the essential improvement obtained by Pompe [Pom03].
Theorem 3.3 (Generalized Korn's inequality). For a fixed λ ∈ ]0, 1[ and positive constants K > 1 define the set Then, for all K > 1 there exists a constant c K > 0 such that for all F ∈ F K we have Proof. In [Pom03, Thm. 2.3] it is shown that (3.4) holds for any given F ∈ F K . Let us denote by c(F ) > 0 the supremum of all possible such constants for the given F . By a perturbation argument it is easy to see that the mapping F → c(F ) is continuous with respect to the L ∞ norm in C 0 (Ω; R d×d ). Since F K is a compact subset of C 0 (Ω; R d×d ) the infimum of c on F K is attained at some F * ∈ F K by Weierstraß' extremum principle. Because of c(F ) ≥ c(F * ), we conclude that (3.4) holds with c K = c(F * ).
We emphasize that estimate (3.4) is not valid if F is not continuous, see [Pom03,Thm. 4.2]. This shows that without the in W 2,p is crucial to control the rate of the strain ∇ . y, which is necessary to handle the thermomechanical coupling. The following corollary combines Theorems 3.1 and 3.3, by using the compact embedding W 2,p (Ω; R d ) ⊂ C 1,λ (Ω; R d ).
Corollary 3.4 (Uniform generalized Korn's inequality on sublevels). Given any C M > 0 there exists a c K > 0 such that for all y ∈ Y id with M(y) ≤ C M we have the generalized Korn inequality By ∂J we denote the Fréchet subdifferential which is defined by Proposition 3.6 (Chain rule for locally semiconvex functionals). Consider a separable reflexive Banach space, a q ∈ ]1, ∞[ with q ′ = q/(q−1), and J : X → R ∞ a lower semicontinuous and locally semiconvex functional. If the functions z ∈ W 1,q ([0, T ]; X) and Ξ ∈ L q ′ ([0, T ]; X * ) satisfy

Time discretization of a regularized problem
Before we construct solution by a suitable time-discretization, we introduce regularizations in two points. Firstly, we add a linear viscous damping which allows us to obtain simple a priori bounds for the strain rate ∇ .
y, because in the first steps of the construction we are not yet in the position to exploiting the generalized Korn inequality of Theorem 3.3. Secondly, we modify the creation of heat through the viscous damping, which in the physically correct form leads to an L 1 source term which can only be handled in the first steps of the construction either.
Similarly, the boundary value problem (4.3b) and (4.4c) for θ k ετ , where y k−1 ετ and y k ετ are given, has a variational structure. For this we define the functions φ C (F, θ) := θ 0 φ(F,θ) dθ and W (F, θ) = 2φ C (F, θ) − θφ(F, θ) to obtain the relation With ∂ 2 θ W (F, θ) = ∂ θ w(F, θ) = −θ∂ 2 θ φ(F, θ) ≥ǫ we see that W (F, ·) is uniformly convex by assumption (2.30c). Thus, we can obtain solutions θ k ετ of (4.3b) and (4.4c) via the minimization problem We emphasize that this staggered scheme is constructed in a very specific way by taking θ = θ k−1 ετ from the previous time step in the mechanics problem for y k ετ , see (4.5). For the construction of θ = θ k ετ from the heat equation we have to use sometimes the explicit (backward) approximations θ k−1 ετ and sometimes the implicit (forward) approximation θ k ετ . Clearly, the former is simpler and it is used in the heat conduction tensor K(∇y k−1 ετ , θ k−1 ετ ) and in the heat production ξ reg ε . It is tempting to use the explicit choice θ k−1 ετ also in the thermo-mechanical coupling term ∂ F φ(∇y k ετ , θ):∇δ τ y k ε (last term in (4.3b)) as it would simplify the energy balance, see Remark 6.1. However, as this term does not have a sign, we would not be able to guarantee positivity of θ k ετ . Thus we are forced to use the more involved implicit term θ → ∂ F φ C (∇y k ε , θ):∇δ τ y k ε in (4.7) instead of the simpler, linear choice θ → θ∂ F φ(∇y k ετ , θ k−1 ετ ):∇δ τ y k ε . This choice may introduce a nonconvexity, so that θ k ετ may not be unique.
The following result states that we can obtain solutions (y k ετ , θ k ετ ) of (4.3) and (4.4) by solving the minimization problems (4.5) and (4.7), alternatingly. For notational simplicity we have written the minimization problem (4.7) for θ with the constraint θ ≥ 0, however, for establishing the Euler-Lagrange (4.3b) and (4.4c) we need to show that non-negativity of θ comes even without imposing the constraint. This will be achieved by minimization over θ ∈ H 1 (Ω) after extending all functionals suitably for θ < 0.
In summary, we conclude that the extended functional in (4.7) is weakly lower semicontinuous and and coercive. Hence, a global minimizers θ * exist and moreover these minimizers solve the associated Euler-Lagrange equation as ∂ θ W (F, θ) = w(F, θ) and To show that all global minimizers are non-negative we test the Euler-Lagrange equation by the negative part θ − * := min{θ * , 0} of θ * , which is still an H 1 function: In the first estimate we have used w k−1 ετ = w(∇y k−1 ετ , θ k−1 ετ ) ≥ǫθ k−1 ετ ≥ 0, ξ reg ε ≥ 0, and θ k ♭,ε,τ ≥ 0 which gives the non-negativity of p 2 , p 4 , and p 7 , while the first and fifth term vanish identically since for θ * > 0 we have θ − * = 0 while for θ * < 0 we have w(F, θ * ) = 0 and ∂ F φ(F, θ * ) = 0 (here we crucially use the implicit structure). Thus, we conclude θ − * = 0 which is equivalent to θ * ≥ 0. Thus, choosing θ k ετ = θ * for any global minimizer of the extended functional we see that it is also a global minimizer of (4.7) and that the Euler-Lagrange equations hold.
The following result provides the basic energy estimates where we will crucially use the carefully chosen semi-implicit scheme defined through the staggered minimization problems (4.5) and (4.7). Here also we will essentially rely regularizing viscous term ε∆ .
y, as R cannot be used because of the missing a priori bound for y k ετ in W 2,p (Ω; R d ). Moreover, we will exploit the fact that we have global minimizers in (4.5) rather than arbitrary solutions of the Euler-Lagrange equations (4.3a). This latter argument works because we have neglected inertial terms in the momentum balance (2.27a) and hence in (4.3a). We refer to [KrR19] to cases where inertial effects are treated but in the isothermal case.
Proof. As y k ετ is a global minimizer, we can insert y = y k−1 ετ as testfunction in (4.5) to obtain the estimate (recall δ τ y k ε = 1 τ (y k ετ −y k−1 ετ )) The proof will be divided into three steps.
Step 2: Dissipation bound. We return to (4.13) and add all estimates from k = 1 to This provides the uniform bound for y ετ in H 1 (I; H 1 (Ω; R d )), and (4.9a) is established.
This finishes the proof of Proposition 4.2.
5 The limit τ → 0 in the regularized problem Using the above a priori estimates for the interpolants we will be able to extract convergent subsequences. First we will observe that the three different types of interpolants have to converge to the same limit. Next we want to pass to the limit in the discretized weak forms of the momentum balance and the heat equation. While most terms can be handled by compactness arguments or weak-convergence methods, there is one term that needs special attention namely the heat-source term ξ reg ε that is quadratic in ∇ . y ε . Thus, it will be a crucial step to show strong convergence of . y ετ in L 2 (I; H 1 (Ω)), which can be done by passing to the limit in a suitable discretized version of the mechanical energy balance (2.20). In this argument we will use the Λ-convexity derived in Proposition 3.2 to relate the mechanical energies M(y k−1 ετ ) and M(y k ετ ). With the definition (4.8) for the three types of interpolants, we see that the following discretized version (5.1) of the momentum balance and heat equations (4.1) and (4.2) holds for the discrete solutions constructed in Proposition 4.1: − div σ vi (∇y ετ , ∇ . y ετ , θ ετ ) + ε∇ .
Proof. The proof consists of five steps.
F and by Minty's trick for the monotone operator induced by h el = H ′ . For a reflexive Banach space X and a hemi-continuous, monotone operator H : X → X * Minty's trick means the implication We apply this for H defined by H(y), z = Q h el (∇ 2 y(t, x)) . . . ∇ 2 z(t, x) dx dt, where X = W 2,p (Q). Clearly, H is hemi-continuous and monotone. Choosing u τ = y ετ the weak equations (5.1a) and (5.2) are interpreted as H( y ετ +σ el (∇y ετ , θ ετ ) : ∇z dx dt + T 0 ℓ τ , z dt.
To use Minty's trick (5.8) we still need to check b τ , y ετ → b, y ε . However, as we have shown above b τ is bounded (and hence weakly converging to b) in L 2 I; H 1 D (Ω) * and y ετ → y ε in L 2 I; H 1 D (Ω) strongly (by (5.4a), the result follows immediately. Hence, we conclude H(y ε ) = b, which is nothing else than the regularized momentum balance (4.1a), (4.2a), and (4.2b).
Step 3: Balance of mechanical energy. For the limit passage in the heat equation we need strong L 2 -convergence of ∇ . y ετ due to the viscous dissipation ξ reg ε (F, . F , θ) that is nonlinear in . F . The strategy is to use the balance of mechanical energy as follows. Rewriting the regularized momentum balance (4.1a), (4.2a), and (4.2b) in the form D .
Step 4: Strong convergence of strain rates. The next step is now to derive a similar mechanical energy balance for the time-discretized solutions, which is better than the previously used estimate (4.11). Passing to the limit τ → 0 from the latter estimate we would arrive at an estimate like (5.9), but with 2R and ε replaced by R and ε/2, respectively.
To improve the discrete bounds used in Proposition 4.2 we can exploit the a priori estimates M(y k ετ ) ≤ K ε , which allow us to use the geodesic Λ-convexity result in Proposition 3.2. Instead of using the minimization property of y k ετ in (4.5) we test the Euler-Lagrange equation (4.3a) with boundary conditions (4.4a) and (4.4b) by y k ετ −y k−1 ετ to obtain τ 2R(y k−1 ετ , δ τ y k ε , θ k−1 ετ ) + τ ε ∇δ τ y k ε 2 where we have the correct factors 2R and ε. To recover the energy values M(y j ετ ) we now eliminate the term involving DM using the Λ-convexity estimate (3.3) with y (1) = y k ετ and y (2) = y k−1 ετ , which yields We now sum this inequality over k = 1, , . . . , N τ and using the interpolants we obtain the integral estimate M(y ετ (T )) + Using the the convergences (5.3) and (5.4) it is immediate to see that the all the terms on the right-hand side converge to the corresponding terms on the right-hand side in (5.9). Now denote the three terms on the left-hand side by I (j) ετ and set I ετ . Using lower semicontinuity arguments (use [FoL07,Thm. 7.5] once again for I (2) ετ ) we find Thus, passing to the liminf on the left-hand side and to the limit on the right-hand side in (5.10) and comparing with (5.9) we obtain Together with (5.11) we conclude that we must have equality in all three cases after "=⇒".
y ε in L 2 (Q; R d×d ) and imply the desired strong convergence ∇ . y ετ → ∇ .
This conclude the proof of Proposition 5.1.
6 Limit passage ε → 0 In this final step of the proof of Theorem 2.2 we have to pass to the limit with the regularization parameter ε → 0. As we are already in the time-continuous setting we are now able to make the formally derived total energy balance (2.21) for E rigorous for all ε > 0. From this we will be able to derive a priori bounds for (y ε , θ ε ) that are independent of ε.
Remark 6.1 (Missing discrete estimate for the total energy). The derivation of the total energy balance is achieved by testing the momentum balance by .
y and the heat equation by the constant function 1. The corresponding step on the time-discrete level would be the test (4.3a) by δ τ y k and (4.3b) by 1. We would be able to use the desirable cancellation of the dissipation, namely ξ reg ε − ξ ≤ 0; however for the coupling terms ∂ F φ(∇y k ετ , θ k−1 ετ ) : δ τ ∇y k ε and ∂ F φ(∇y k ετ , θ k ετ ) : δ τ ∇y k ε , which arise from (4.3a) and (4.3b) respectively, we do not have any way to estimate the first against the second. Recall that we were forced to use the explicit/forward value θ k ετ to maintain positivity of the temperature.
To exploit the balance of the total energy we have to strengthen the assumption on the leading ℓ(t), i.e. the functions g, and f , in (2.30g), namely g ∈ W 1,1 (I; L 2 (Ω; R d )), f ∈ W 1,1 (I; L 2 (Γ N ; R d )). (6.1) This implies that t → ℓ(t) lies in W 1,1 I; H 1 Γ D (Ω; R d ) * , which is what we will only need.
The new ε-independent estimates on ∇ .
y ε in L 2 (Q) will be obtain by exploiting the Pompe's generalized Korn's inequality (cf. [Pom03]) as prepared in Theorem 3.3 above.
with q from (2.30a), where again sym(·) denotes the symmetric part of a (d×d)-matrix.
Proof. We proceed in two steps that are close to estimates we have done in the timediscrete setting.
The importance is the cancellation of the term ∂ F φ : ∇ . y ε and that the difference of the dissipation integrals has a sign.
Step 2: Estimate for the strain rate ∇ . y ε . We return to the mechanical energy balance (5.9) on the interval I = [0, T ]. We recall that the dissipation function ξ(F, . F , θ) is assumed to control the symmetric part of F ⊤ .
F only, namely Using our a priori bounds on M(y ε (t)), we can apply the generalized Korn's inequality a prepared in Corollary 3.4 with y = y ε (t, ·) and v = .
For the deformation y ε we have all the estimates we need for passing to the limit. But we still need good a priori estimates for the temperature. Here the problem arises that the heating arising through the viscous dissipation ξ(∇y ε , ∇ . y ε , θ ε ) is only bounded in L 1 (Q). So, obtaining improved estimates we have to invoke special test functions developed by Boccardo and Gallouët [BoG89] for parabolic equations with measure-valued right-hand sides.
We are now in the position to pass to the limit ε → 0 in the regularized system (4.1)-(4.2), and thus provide the proof of our main existence result presented in Theorem 2.2. The approach is close to the convergence result presented in Proposition 5.1: first we extract converging subsequences and then pass to the limit in the mechanical momentum balance. This also provides the necessary strong convergence of the the strain rates that is needed to eventually pass to the limit in the heat heat equation.
Proof. The proof follows the lines of the proof of Proposition 5.1, so we do not repeat all details of the arguments.
Step 2: Convergence in the mechanical equation. The limit passage in the momentum balance (4.1a)-(4.2) works as before, again using the Minty trick (5.8). Of course, the additional regularizing viscosity term ε∇ .
We see that the first term converges by (6.10), while the second term converges by the weak convergence V ε ⇀ V and the strong convergence D(C ε , θ ε )V → D(C, θ)V (as D is bounded and the arguments converge pointwise). Similarly, δ(ε) → 0 by Lebesgue's dominated convergence theorem, and thus we conclude the strong convergence V ε −V L 2 (Q) → 0.

Because of
Step 4, we know V ε → V strongly in L 2 (Q; R d×d sym ). Hence, we have g ε → g := K|V | 2 in L 1 (Q) and may assume, after extracting another subsequence, V ε (t, x) → V (t, x) a.e. in Q. By the uniform/pointwise convergence of C ε and θ ε for any v ∈ C 0 (Q) we obtain y, θ)v a.e. in Q.
As the majorants g ε v L ∞ (Q) converge to g v L ∞ (Q) in L 1 (Q) the generalized dominated convergence theorem implies convergence of the second term in (6.11).
In the third term we have weak convergence of ∇ .
y ε and strong convergence of v∂ F φ(∇y ε , θ ε ). Similarly, the remaining four terms converge to the desired limits. Thus, we have shown that (y, θ) satisfy (2.27b), which finishes the proof of Proposition 6.4. Remark 6.5 (Strong convergence of y ετ and y ε ). Strengthening monotonicity of h el , cf.
y| 2 with a mass density ̺ = ̺(x) > 0 leads to an inertial force ̺ .. y in the momentum equation (2.13a), which would make the nonlinear problem hyperbolic. It is generally recognized as analytically very troublesome. Here, it would work for isothermal situation like in Corollary 2.3 if we would be able to work with weak convergence, i.e. H needs to be quadratic (p = 2). Staying with H depending on the second gradient ∇ 2 y we would be forced to give up the determinant constraint det ∇y > 0, which is indeed possible if heat conduction is not considered. Alternatively, one may take H quadratic but coercive in Hilbert space norms H s (Ω) with s > 1 + d/2, such that H s (Ω) still embeds into C 1,α for some α > 0, cf. also [KrR19,Ch. 9.3]. In the anisothermal situation, it seems difficult to ensure that the acceleration .. y ∈ L 2 (I; H 1+κ (Ω; R d×d ) stays in duality with the velocity .
y. The regularity seems difficult and the higher-order viscosity is inevitably very nonlinear to comply with frame-indifference while the corresponding generalization of Korn's inequality does not seem available.
y + ∇µ · M(x, ∇y, c)∇µ (6.12c) with σ vi as in (2.13a), c v (F, c, θ) = −θ∂ 2 θθ ψ(F, c, θ), and ξ from (2.10). In (6.12b), the variable µ is called a chemical potential. One can also augment the model by some inelastic (plastic or creep-type) strain like in [RoS18] where also the inertial forces have been involved and the viscosity ignored but the concept of small elastic strains imposed as a modeling assumption.