Existence results for a coupled viscoplastic-damage model in thermoviscoelasticity

In this paper we address a model coupling viscoplasticity with damage in thermoviscoelasticity. The associated PDE system consists of the momentum balance with viscosity and inertia for the displacement variable, at small strains, of the plastic and damage flow rules, and of the heat equation. It has a strongly nonlinear character and in particular features quadratic terms on the right-hand side of the heat equation and of the damage flow rule, which have to be handled carefully. We propose two weak solution concepts for the related initial-boundary value problem, namely `entropic' and `weak energy' solutions. Accordingly, we prove two existence results by passing to the limit in a carefully devised time discretization scheme. Finally, in the case of a prescribed temperature profile, and under a strongly simplifying condition, we provide a continuous dependence result, yielding uniqueness of weak energy solutions.


1.
Introduction. This paper focuses on two phenomena related to inelastic behavior in materials, namely damage and plasticity. Damage can be interpreted as a degradation of the elastic properties of a material due to the failure of its microscopic structure. Such macroscopic mechanical effects take their origin from the formation of cracks and cavities at the microscopic scale. They may be described in terms of an internal variable, the damage parameter, on which the elastic modulus depends, in such a way that stiffness decreases with ongoing damage. Plasticity produces residual deformations that remain after complete unloading.
Recently, models combining plasticity with damage have been proposed in the context of geophysical modeling [36,38] and, more in general, within the study of the thermomechanics of damageable materials under diffusion [37]. Perfect plasticity is featured in [36,38], where the evolution of the damage variable is governed by viscosity, i.e. it is rate-dependent. Conversely, in [37] damage evolves rate-independently, while the evolution of plasticity is rate-dependent. In a different spirit, a fully rate-independent system for the evolution of the damage parameter, coupled with a tensorial variable which stands for the transformation strain arising during damage evolution, is analyzed in [6]. Finally, let us mention the coupled elastoplastic damage model from [1,2], analyzed in [9,10,11]. While the first two papers deal with the fully rate-independent case in [9], in [11], the model is regularized by adding viscosity to the damage flow rule, while keeping the evolution of the plastic tensor rate-independent. A vanishing-viscosity analysis is then carried out, leading to an alternative solution concept ('rescaled quasistatic viscosity evolution') for the rateindependent elastoplastic damage system. A common feature in [36,38,9,11] is that the plastic yield surface, and thus the plastic dissipation potential, depends on the damage variable.
In this paper we aim to bring temperature into the picture. Thermoplasticity models, both in the case of rate-independent evolution of the plastic variable and of rate-dependent one, have been the object of several studies, cf. e.g. [24,25,26,3,4,35,20,32]. In recent years, there has also been a growing literature on the analysis of (rate-independent or rate-dependent) damage models with thermal effects: We quote among others [5,33,30,19,31,27]. As for models coupling plasticity and damage with temperature, one of the examples illustrating the general theory developed in [37] concerns geophysical models of lithospheres in short time scales, which couple the (small-strain) momentum balance, damage, rate-dependendent plasticity, the heat equation, as well as the porosity and the water concentration variables. Here we shall neglect the latter two variables and tackle the weak solvability and the existence of solutions, for a (fully rate-dependent) viscoplastic (gradient) damage model, with viscosity and inertia in the momentum balance (the first according to Kelvin-Voigt rheology), and with thermal effects encompassed through the heat equation, whereas in [37] the enthalpy equation was analyzed, after a transformation of variables. We plan to address the vanishing-viscosity and inertia analysis for our model, and discuss the weak solution concept thus obtained, in a future contribution.
In what follows, we shortly comment on the model. Then, we illustrate the mathematical challenges posed by its analysis, motivate a suitable regularization, and introduce the two solution concepts, for the original system and its regularized version, for which we will prove two existence results.
which provides a decomposition of the linearized strain tensor E(u) = 1 2 (∇u+∇u ) into the sum of the elastic and plastic strains e and p. In fact, e ∈ M d×d sym (the space of symmetric (d×d)-matrices), while p ∈ M d×d D (the space of symmetric (d×d)-matrices with null trace).
-The momentum balance with the stress given by σ = D(z)ė + C(z)e − C(z)Eϑ in Ω × (0, T ) according to Kelvin-Voigt rheology for materials subject to thermal expansion (with E the matrix of the thermal expansion coefficients). Here, F is a given body force. Observe that both the elasticity and viscosity tensors C and D depend on the damage parameter z, but we restrict to incomplete damage. Namely, the tensors C and D are definite positive uniformly w.r.t. z, meaning that the system retains its elastic properties even when damage is maximal. -The damage flow rule for z ∂R(ż) +ż + A s (z) + W (z) − 1 2 C (z)e : e + ϑ in Ω × (0, T ), where the dissipation potential (density) R : R → [0, +∞] defined by R(η) := |η| if η ≤ 0, +∞ otherwise, encompasses the unidirectionality in the evolution of damage. We denote by ∂R : R ⇒ R its subdifferential in the sense of convex analysis. As in several other damage models, we confine ourselves to a gradient theory. However, along the lines of [23] and [11], we adopt a special choice of the gradient regularization, i.e. through the s-Laplacian operator A s , with s > d 2 . This technical condition ensures the compact embedding H s (Ω) C 0 (Ω) for the associated Sobolev-Slobodeckij space H s (Ω) := W s,2 (Ω). Furthermore, a key role in the analysis, especially for the regularized system with the flow rule (3c) ahead, will be played by the linearity of the operator A s . As in [11], the term W (z) will have a singularity at z = 0, which will enable us to prove the positivity of the damage variable. Combining this with the unidirectionality constraintż ≤ 0 a.e. in Ω×(0, T ), we will ultimately infer that all z emanating from an initial datum z 0 ≤ 1 take values in the physically admissible interval [0, 1].
-The flow rule for the plastic tensor reads ∂ṗH(z, ϑ;ṗ) +ṗ σ D in Ω × (0, T ), with σ D the deviatoric part of the stress tensor σ. Here, the plastic dissipation potential (density) H depends on the plastic strain rateṗ but also (on the space variable x), on the temperature, and on the damage variable; the symbol ∂ṗH denotes its convex subdifferential w.r.t. the plastic rate. Typically, one may assume that the plastic yield surface decreases as damage increases, although this monotonicity property will not be needed for the analysis developed in this paper. Observe that (1e) is in fact a viscous regularization of the flow rule ∂ṗH(z, ϑ;ṗ) σ D in Ω × (0, T ) of perfect plasticity. -The heat equatioṅ ϑ − div(κ(ϑ)∇ϑ) = G + Dė :ė − ϑC(z)E :ė + R(ż) + |ż| 2 − ϑż + H(z, ϑ;ṗ) +ṗ :ṗ (1f) in Ω×(0, T ), with κ ∈ C 0 (R + ) the heat conductivity coefficient and G a given, positive, heat source. We will supplement system (1) with the boundary conditions where n is the outward unit normal to ∂Ω, e Γ Dir , Γ Neu the Dirichlet/Neumann parts of the boundary, respectively, and with initial conditions. In fact, system (1) can be seen as the extension of the model considered in [11], featuring the static momentum balance and a mixed rate-dependent/rateindependent character in the evolution laws for the damage/plastic variables, respectively, to the case where viscosity is included in the plastic flow rule and in the momentum balance, the latter also with inertia, and the evolution of temperature is also encompassed.
System (1,2) can be rigorously derived following, e.g., the thermomechanical modeling approach by M. Frémond [15,Chap. 12]. In this way one can also verify its compliance with the first and second principle of Thermodynamics, hence its thermodynamical consistency.
1.2. Analytical challenges and weak solution concepts. Despite the fact that the 'viscous' plastic flow rule (1e) does not bring all the technical difficulties attached to perfect plasticity (cf. [12], see also [35] for the coupling with temperature), the analysis of system (1, 2) still poses some mathematical difficulties. Namely, (1): The overall nonlinear character of (1) and, in particular, the quadratic terms on the right-hand side of the damage flow rule (1d), and on the right-hand side of the heat equation (1f). Both sides are, thus, only estimated in L 1 (Ω×(0, T )) as soon asė,ż, andṗ are estimated in L 2 (Ω×(0, T ); M d×d sym ), L 2 (Ω), and L 2 (Ω×(0, T ); M d×d D ), respectively, as guaranteed by the dissipative estimates associated with (1).
Observe that the particular character of the momentum balance, where the elasticity and the viscosity contributions only involve the elastic part of the strain e and its rateė, instead of the full strain E(u) and strain rate E(u), does not allow for elliptic regularity arguments which could at least enhance the spatial regularity/summability of the right-hand side of the damage flow rule.
(2): Another obstacle is given by the presence of the unbounded maximal monotone operator ∂R in the damage flow rule. Because of this, no comparison estimates can be performed. In particular, a pointwise formulation of (1d) would require a separate estimate of the terms A s (z) and of (a selection in) ∂R(ż). This cannot be obtained by standard monotonicity arguments due to the nonlocal character of the operator A s . All of these issues shall be reflected in the weak solution concept for system (1,2) proposed in the forthcoming Definition 2.3 and referred to as 'entropic solution'. This solvability notion consists of the so-called entropic formulation of the heat equation, and of a weak formulation of the damage flow rule, in the spirit of the Karush-Kuhn-Tucker conditions. The entropic formulation originates from the work by E. Feireisl in fluid mechanics [13] and has been first adapted to the context of phase transition systems in [14], and later extended to damage models in [31]. It is given by an entropy inequality, formally obtained by dividing the heat equation by ϑ and testing the resulting relation by a sufficiently regular, positive test function (cf. the calculations at the beginning of Section 2.3), combined with a total energy inequality. The weak formulation of the damage flow rule has been first proposed in the context of damage modeling in [17,18]: the subdifferential inclusion for damage is replaced by a one-sided variational inequality, with test functions reflecting the sign constraint imposed by the dissipation potential R, joint with a (mechanical) energy-dissipation inequality also incorporating contributions from the momentum balance and the plastic flow rule.
Clearly, one of the analytical advantages of the entropy inequality for the heat equation, and of the one-sided inequality for the damage flow rule, is that the troublesome quadratic terms on the right-hand sides of (1d) and of (1f) feature as multiplied by a negative test function, cf. (39a) and (42) ahead, respectively. This allows for upper semicontinuity arguments in the limit passage in suitable approximations of such inequalities. Instead, the total and mechanical energy inequalities can be obtained by lower semicontinuity techniques.
We will also consider a regularized version of system (1,2), where the damage flow rule features the additional term A s (ż), modulated by a positive constant ν, and, accordingly, the term νa s (ż,ż) (with a s the bilinear form associated with A s ) occurs on the right-hand side of the heat equation. This leads to the following regularized thermoviscoelastoplastic damage system, posed in Ω × (0, T ), ϑ − div(κ(ϑ)∇ϑ) =G + Dė :ė − ϑC(z)E :ė + R(ż) + |ż| 2 +νa s (ż,ż) − ϑż + H(z, ϑ;ṗ) +ṗ :ṗ (3e) withν = ν/|Ω|, supplemented with the boundary conditions For this regularized system we will be able to show the existence of an enhanced type of solution. It features a conventional weak formulation of the heat equation with suitable test functions, and the formulation of the damage flow rule as a subdifferential inclusion in H s (Ω) * . To obtain the latter, a key role is played by the regularizing term A s (ż), ensuring thatż ∈ H s (Ω) a.e. in (0, T ). Thanks to this feature, it is admissible to test the subdifferential inclusion rendering (3c) byż itself. This is at the core of the validity of the total energy balance. That is why, also in accordance with the nomenclature from [31,32], we shall refer to these enhanced solutions as weak energy solutions.
1.3. Our results. We will prove the existence of entropic solutions, see Theorem 2.4, and of weak energy solutions, see Theorem 2.5, to (the Cauchy problems for) systems (1,2) and (3,4), respectively, by passing to the limit in a time-discretization scheme carefully devised in such a way as to ensure the validity of discrete versions of the entropy and energy inequalities along suitable interpolants of the discrete solutions. One of our standing assumptions will be a suitable growth of the heat conductivity coefficient κ, namely κ(ϑ) ∼ ϑ µ for some µ > 1.
Mimicking the calculations from [14,31], we will exploit (5) in a key way to derive an estimate for ϑ in L 2 (0, T ; H 1 (Ω)) by testing the (discrete) heat equation by a suitable negative power of ϑ.
Under a more restrictive condition on µ, in fact depending on the space dimension d, (i.e., µ ∈ (1, 2) if d = 2, µ ∈ (1, 5 3 ) if d = 3), we will also be able to obtain a BVin-time estimate for ϑ, with values in a suitable dual space, which will be at the core of the proof of the enhanced formulation of the heat equation for the weak energy solutions to system (3,4). Concerning the physical interpretation of our growth conditions, we refer to [22] for a discussion of experimental findings suggesting that a class of polymers exhibit a subquadratic growth for κ.
Finally, with Proposition 1 we will provide a continuous dependence estimate, yielding uniqueness, for the weak energy solutions to (3,4) in the case of a prescribed temperature profile, and with a plastic dissipation potential independent of the state variables z and ϑ. Plan of the paper. In Section 2 we fix all our assumptions, motivate and state our two weak solvability notions for systems (1,2) & (3,4), and finally give the existence Theorems 2.4 & 2.5, and the continuous dependence result Proposition 1. In Section 3 we set up a common time discretization scheme for systems (1,2) and (3,4), and prove the existence of discrete solutions, while Section 4 is devoted to the derivation of all the a priori estimates on the approximate solutions, obtained by interpolation of the discrete ones. In Section 5 we conclude the proofs of Thms. 2.4 & 2.5 by passing to the time-continuous limit, while in Section 6 we perform the proof of Prop. 1.
We conclude by fixing some notation that shall be used in the paper.
Notation (General notation). Throughout the paper, R + shall stand for (0, +∞). For a given z ∈ R, we will use the notation (z) + for its positive part max{z, 0}. We will denote by M d×d (M d×d×d×d ) the space of (d×d) ((d×d×d×d), respectively) matrices. We will consider M d×d endowed with the Frobenius inner product A : B := ij a ij b ij for two matrices A = (a ij ) and B = (b ij ), which induces the matrix norm | · |. Therefore, we will often write |A| 2 in place of A : A. The symbol M d×d sym stands for the subspace of symmetric matrices, and M d×d D for the subspace of symmetric matrices with null trace. We recall that M d×d sym = M d×d D ⊕ RI (I denoting the identity matrix), since every η ∈ M d×d sym can be written as η = η D + tr(η) d I with η D the orthogonal projection of η into M d×d D . We will refer to η D as the deviatoric part of η.
For a given Banach space X, the symbol ·, · X will stand for the duality pairing between X * and X; if X is a Hilbert space, (·, ·) X will denote its inner product. For simpler notation, we shall often write · X both for the norm on X, and on the product space X × . . . × X. With the symbol B 1,X (0) we will denote the closed unitary ball in X. We shall use the symbols Finally, throughout the paper we will denote various positive constants depending only on known quantities by the symbols c, c , C, C , whose meaning may vary even within the same line. Furthermore, the symbols I i , i = 0, 1, ..., will be used as place-holders for several integral terms (or sums of integral terms) occurring in the various estimates: we warn the reader that we will not be self-consistent with the numbering, so that, for instance, the symbol I 1 will have different meanings.
2. Setup and main results for the thermoviscoelastoplastic damage system. After fixing the setup for our analysis in Section 2.1, in Sec. 2.2 we motivate the notion of 'weak energy' solution to system (3,4)  We will denote by |Ω| the Lebesgue measure of Ω. On the Dirichlet part Γ Dir , assumed with H d−1 (Γ Dir ) > 0, we shall prescribe the displacement, while on Γ Neu we will impose a Neumann condition on the displacement. The trace of a function v on Γ Dir or Γ Neu shall be still denoted by the symbol v. Sobolev spaces, the s-Laplacian, Korn's inequality. In what follows, we will use the notation H 1 Dir (Ω; R d ) := {u ∈ H 1 (Ω; R d ) : u| Dir = 0}. The symbol W 1,p Dir (Ω; R d ), p > 1, shall denote the analogous W 1,p -space. Further, we will use the notation and analogously for W 1,p − (Ω). Throughout the paper, we shall extensively resort to Korn's inequality (cf. [16]): for every 1 < p < ∞ there exists a constant C K = C K (Ω, p) > 0 such that there holds We will denote by H s (Ω) the Sobolev-Slobodeckij space W s,2 (Ω), with s ∈ d 2 , 2 . We will also use the notation H s + (Ω) = {z ∈ H s (Ω) : z ≥ 0 in Ω} and the analogously defined notation H s − (Ω). We recall that H s (Ω) is a Hilbert space, with inner product (z 1 , z 2 ) H s (Ω) := (z 1 , z 2 ) L 2 (Ω) + a s (z 1 , z 2 ), where the bilinear form a s (·, ·) is defined by We denote by A s : H s (Ω) → H s (Ω) * the associated operator Kinematic admissibility and stress. Given a function w ∈ H 1 (Ω; R d ), we say that a triple (u, e, p) is kinematically admissible with boundary datum w, and write (u, e, p) ∈ A(w), if As for the elasticity and viscosity tensors, we will suppose that where Lin(M d×d sym ) denotes the space of linear operators from M d×d sym to M d×d sym . Furthermore, we will suppose that for every x ∈ Ω the map z → C(x, z) is continuously differentiable on R and fulfills Finally, for technical reasons (cf. Remark 2 later on) it will be convenient to require that the map z → C(x, z) is convex, i.e. for every A ∈ M d×d sym there holds It follows from the convexity (2.(C, D)3) that whence (C (x, z 1 )−C (x, z 2 ))(z 1 −z 2 )A : A ≥ 0. In particular, due to the first of (2.(C, D)3), we find that Finally, we also suppose that the thermal expansion tensor fulfills Observe that with (2.(C, D, E)) and (2.E) we encompass in our analysis the case of an anisotropic and inhomogeneous material. External heat sources. For the volume and boundary heat sources G and g we require Indeed, the positivity of G and g is crucial for obtaining the strict positivity of the temperature ϑ.
Body force and traction. Our basic conditions on the volume force F and the assigned traction f are For technical reasons, in order to allow for a non-zero traction f , we will need to additionally require a uniform safe load type condition. Observe that this kind of assumption usually occurs in the analysis of perfectly plastic systems. In the present context, it will play a pivotal role in the derivation of the First a priori estimate for the approximate solutions constructed by time discretization, cf. the proof of Proposition 3 later on as well as [32,Rmk. 4.4] for more detailed comments. Namely, we impose that there exists a function : [0, T ] → L 2 (Ω; M d×d sym ), with ∈ W 1,1 (0, T ; L 2 (Ω; M d×d sym )) and D ∈ L 1 (0, T ; L ∞ (Ω; M d×d D )), solving for almost all t ∈ (0, T ) the following elliptic problem When not explicitly using (2.L 2 ), to shorten notation we will incorporate the volume force F and the traction f into the induced total load, namely the function for all u ∈ H 1 Dir (Ω; R d ), which fulfills L ∈ L 2 (0, T ; H 1 Dir (Ω; R d ) * ) in view of (2.L 1 ). Dirichlet loading. We will suppose that the hard device w to which the body is subject on Γ Dir is the trace on Γ Dir of a function, denoted by the same symbol, fulfilling w ∈ L 1 (0, T ; W 1,∞ (Ω; R d )) ∩ W 2,1 (0, T ; H 1 (Ω; R d )) ∩ H 2 (0, T ; L 2 (Ω; R d )) . (2.w) Also condition (2.w) will be used in the proof of Prop. 3 ahead; we again refer to [32,Rmk. 4.4] for more comments. The weak formulation of the momentum balance. The variational formulation of (1b), supplemented with the boundary conditions (2a), reads for almost all t ∈ (0, T ) We will often use the short-hand notation −div Dir for the elliptic operator defined by The plastic dissipation potential. Our assumptions on the multifunction K : involve the notions of measurability, lower semicontinuity, and upper semicontinuity for general multifunctions. For such concepts and the related results, we refer, e.g., to [8]. Hence, we suppose that is continuous for almost all x ∈ Ω.
(2.K 1 ) D for all (z, ϑ) ∈ R × R + and for almost all x ∈ Ω, Therefore, the support function associated with the multifunction K, i.e.
Finally, we also introduce the plastic dissipation potential H : L 1 (Ω)×L 1 (Ω; R + )× L 1 (Ω; M d×d D ) given by From now on, throughout the paper we will most often omit the x-dependence of the tensors C, D, E, and of the dissipation density H. Nonlinearities in the damage flow rule. Along the footsteps of [11], we will suppose that W ∈ C 2 (R + ) is bounded from below and fulfills The latter coercivity condition will play a key role in the proof that the damage variable z takes values in the feasible interval [0, 1]. In this way, we will not have to include the indicator term I [0,1] in the potential energy. This will greatly simplify the analysis of the damage flow rule. Furthermore, we shall require that Observe that (2.W 2 ) is equivalent to imposing that the function z → W (z) + λ W 2 |z| 2 =: β(z) is convex. Therefore, we have the convex/concave decomposition W (z) = β(z) − λ W 2 |z| 2 with β ∈ C 2 (R + ), convex, and fulfilling z 2d β(z) → +∞ as z ↓ 0.
Let us point out that (19) will be expedient in devising the time discretization scheme for the (regularized) thermoviscoelastoplastic damage system, in such a way that its solutions comply with the discrete version of the total energy inequality. We refer to Remark 2 ahead for more comments. Cauchy data. We will supplement the thermoviscoelastoplastic damage system with initial data ϑ 0 ∈ L 1 (Ω), fulfilling the strict positivity ∃ ϑ * > 0 : inf and such that log(ϑ 0 ) ∈ L 1 (Ω). (20d) In the remainder of this section, we shall suppose that the functions C, . . . , W , the data G, . . . , w, and the initial data (u 0 ,u 0 , e 0 , p 0 , z 0 , ϑ 0 ) fulfill the conditions stated in Section 2.1. We will first address the weak solvability of the regularized system in Sec. 2.2, and then turn to examining the non-regularized one in Sec. 2.3. We will then state our existence results for both in Sec. 2.4.

2.2.
Energetics and weak solvability for the (regularized) thermoviscoelastoplastic damage system. Prior to stating the precise notion of weak solution for the regularized thermoviscoelastoplastic damage system in Definition 2.1 ahead, we formally derive the mechanical & total energy balances associated with systems (1, 2) and (3,4) (in the ensuing discussion, we shall take the parameter ν ≥ 0). The mechanical and total energy balances. The free energy of the system is given by The total energy balance can be (formally) obtained by testing the momentum balance (3b) by (u−ẇ), the damage flow rule (3c) byż, the plastic flow rule (3d) byṗ, the heat equation (3e) by 1, adding the resulting relations and integrating in space and over a generic interval (s, t) ⊂ (0, T ). Indeed, the tests of (3b), (3c), and (3d) yield Now, taking into account that E(u) =ė +ṗ by the kinematical admissibility condition, rearranging some terms one has that We substitute this in (22) and note that t s Ω (D(z)ė + C(z)e − ϑC(z)E) :ṗ dx dr = t s Ω σ D :ṗ dx dr sinceṗ ∈ M d×d D , so that the last term on the right-hand side of (22) cancels out. Furthermore, by the chain rule we have that Collecting all of the above calculations, we obtain the mechanical energy balance, featuring the kinetic, dissipated, and mechanical energies (23) Let us highlight that (23) will also have a significant role for our analysis.
Hence, we sum (23) with the heat equation (3e) tested by 1 and integrated in time and space. We observe the cancelation of some terms, in particular noting thatν t s Ω a s (ż,ż) · 1 dx dr = ν t s a s (ż,ż) dr .
All in all, we conclude the total energy balance Weak energy solutions for the regularized system. With the following definition (where the conditions from Sec. 2.1 are tacitly assumed), we fix the properties of the weak solution concept for the regularized thermoviscoplastic damage system. Let us mention in advance that, in addition to the conventional weak formulations of the momentum balance and of the heat equation (in the latter case, with test functions with suitable regularity and summability properties), we will require the validity of the plastic flow rule pointwise (almost everywhere) in Ω × (0, T ), and that of the damage flow rule as a subdifferential inclusion in H s (Ω) * . It will be in fact possible to obtain the latter by exploiting the additional, strongly regularizing term A s (ż) in the damage flow rule (3c). In this connection, we now introduce the dissipation potential, defined on H s (Ω) and induced by R, namely and its 'viscous' regularization, where R is augmented by the (squared) L 2 (Ω)-norm We will denote by ∂R : H s (Ω) ⇒ H s (Ω) * and ∂R 2 : H s (Ω) ⇒ H s (Ω) * the the convex analysis subdifferentials of R and R 2 , respectively.
Definition 2.1 (Weak energy solutions to the regularized thermoviscoelastoplastic damage system). Given initial data (u 0 ,u 0 , e 0 , z 0 , p 0 , ϑ 0 ) fulfilling (20), we call a quintuple (u, e, z, p, ϑ) a weak energy solution to the Cauchy problem for system (3,4), supplemented with the boundary conditions (2a, 4), if (u, e, z, p, ϑ) comply with the initial conditions and with -the kinematic admissibility condition -the weak formulation (13) of the momentum balance (3b); -the feasibility and unidirectionality constraints and the subdifferential inclusion for damage evolution -the plastic flow rule -the strict positivity of ϑ: and the weak formulation of the heat equation (3e) for every test function ϕ ∈ W 1,∞ (Ω): Since the 'viscous' contributionż Nevertheless, as the spaces (H s (Ω), L 2 (Ω), H s (Ω) * ) form a Hilbert triple, in what follows we will omit the symbol J. Therefore, (31) rewrites as a.e. in (0, T ). We refer to the previously developed calculations, leading to the mechanical and total energy balances (23) and (24), for the proof of the following result. Lemma 2.2. Let (u, e, z, p, ϑ) be a weak energy solution to (the Cauchy problem for) system (3,4). Then, for every 0 ≤ s ≤ t ≤ T the functions (u, e, z, p, ϑ) comply with the mechanical energy balance (23) and with for almost all t ∈ (0, T ) and for t = 0, and in that case (36) coincides with the total energy balance (24).

Entropic solutions for the thermoviscoelastoplastic damage system.
For the (Cauchy problem associated with the) non-regularized system (1, 2), we will be able to prove an existence result only for a solution concept containing much less information than that from Def. 2.1. In particular, we will notably weaken the formulations of the heat equation (1f), given in terms of an entropy inequality joint with the total energy inequality, and of the damage and plastic flow rules. In order to motivate Definition 2.3 ahead, we develop some preliminary considerations on the weak formulation of the damage and plastic flow rules, and on the entropy inequality. The latter will be formally obtained from the heat equation (1f) assuming the strict positivity of the temperature ϑ, which shall be rigorously proved (cf. Prop. 2 ahead). The entropy inequality. It can be formally obtained by multiplying the heat equation (1f) by ϕ/ϑ, with ϕ a smooth and positive test function. Integrating in space and over a generic interval (s, t) ⊂ (0, T ) leads to the identity The entropy inequality (42) ahead is indeed the ≥-estimate in (37), with positive test functions, and where the first integral on the left-hand side is integrated in time.
Weak solvability of the damage and plastic flow rules. Setting ξ := ϑ − 1 2 C (z)e : e −ż − A s (z) − W (z), the subdifferential inclusion (1d) for damage evolution reformulates as ξ ∈ ∂R(ż) in Ω × (0, T ). By the 1-homogeneity of R, the latter is in turn equivalent to the system of inequalities Therefore, the damage flow rule (1d) could be formulated in terms of the constraints (30), of the integrated (in space) version of (38a) with test functions ζ ∈ H s − (Ω), and of the integrated (in space and time) version of (38b), which can be interpreted as an energy-dissipation inequality. Namely, To our knowledge this formulation, inspired by the Karush-Kuhn-Tucker conditions, was first proposed in the context of damage in [17]. Analogously, the plastic flow rule (1e) can be weakly formulated in terms of the system of inequalities Entropic solutions for the non-regularized system. We are now in the position to give our (extremely) weak solution concept for the initial-boundary value problem associated with system (1), where, in addition to the entropic formulation of the heat equation, the damage and plastic flow rules shall be formulated only through the Kuhn-Tucker type variational inequalities (39a) and (40a), combined with the mechanical energy inequality (41) ahead. Observe that the latter corresponds to the sum of the energy-dissipation inequalities for damage and plastic evolution (39b) and (40b), with the weak formulation of the momentum balance tested by (u−ẇ) and integrated in time.

Remark 1.
A few comments on the above definition are in order: 1. While for weak energy solutions it is possible to a posteriori deduce the validity of mechanical and total energy balances via suitable tests, here the upper energy inequalities (41) and (43) have to be both claimed, as neither of them follows from the other items of the definition.
2. Observe that, subtracting the weak momentum balance (13), (legally) tested by (u−ẇ) and integrated in time, from the mechanical energy inequality (41), it would be possible to deduce the joint energy-dissipation inequality for damage and plastic evolution only under the validity of the chain rule However, (45) can be only formally written: indeed, observe that the summability properties z ∈ H 1 (0, T ; L 2 (Ω)) and e ∈ L ∞ (0, T ; L 2 (Ω; M d×d sym )) do not ensure that C (z)że:e ∈ L 1 (Q).
That is why, in Definition 2.3 we only claim the validity of the full inequality (41), which shall be obtained via lower semicontinuity arguments, by passing to the limit in a discrete version of it.

Combining the information that
. The regularity and summability ϕ ∈ L ∞ (0, T ; W 1,∞ (Ω)) ∩ H 1 (0, T ; L 6/5 (Ω)) required for every admissible test function for the entropy inequality (42) guarantee the integrals log(ϑ)φ dx dr and Ω log(ϑ)ϕ dx are well defined since, in particular, L 6/5 (Ω)) is the dual of L 6 (Ω), which is the smallest Lebesgue space into which H 1 (Ω) embeds in d = 3. Furthermore, with (42) we are also tacitly claiming the summability properties 4. We refer to [31,Rmk. 2.6] for some discussion on the consistency between the entropic (consisting of the entropy and total energy inequalities) and the classical formulations of the heat equation.
2.4. Existence results. We start by stating the existence of entropic solutions to the (non-regularized) system (1, 2) under a mild growth condition on the thermal conductivity κ. Observe that, with (47) below we will exhibit a precise lower bound for the temperature in terms of quantities related to the material tensors C, D, E.
Then, for every (u 0 ,u 0 , e 0 , z 0 , p 0 , ϑ 0 ) satisfying (20) there exists an entropic solution (u, e, z, p, ϑ) to the Cauchy problem for system (1,2) such that, in addition, 1. there exists ζ * > 0 such that 2. ϑ complies with the positivity property where ϑ * > 0 is from (20d) and C * : Under a more restrictive growth condition on κ, we are able to establish the existence of weak energy solutions for the regularized thermoviscoelastoplastic damage system. Let us also point out that we will be able to enhance the temporal regularity of the temperature and obtain a variational formulation of the heat equation with a wider class of test functions. We will also show the validity of the entropy inequality (42): this is a result on its own, as (42) cannot be inferred from the weak formulation (34) of the heat equation, not even in the enhanced form established with Thm. 2.5.
The proof of Theorems 2.4 & 2.5 shall be developed throughout Secs. 3-5 by passing to the limit in a carefully devised time discretization scheme, along the footsteps of the analysis previously developed in [31,32].
As we will illustrate in Remark 3, it would also be possible to prove the existence of entropic solutions to the (Cauchy problem for the) thermoviscoelastoplastic system (1, 2) by an alternative method. Namely, we could pass to the limit as the regularization parameter ν ↓ 0 in the weak energy formulation of the regularized system, featuring a family (κ ν ) ν of thermal conductivities fulfilling (2.κ 2 ) and suitably converging as ν ↓ 0 to a function κ that only complies with (2.κ 1 ). However, to avoid overburdening the paper, we have chosen not to develop this asymptotic analysis.
2.5. Continuous dependence on the external and initial data in the case of a prescribed temperature profile. Let us now confine the discussion to the regularized system to the case the temperature profile is prescribed. Namely we consider the the PDE system consisting of the momentum balance (3b), of the regularized damage flow rule (3c), and of the plastic flow rule (3d), with a given temperature In this context, weak energy solutions fulfill the weak momentum balance (13), the subdifferential inclusion for damage evolution (31) (i.e., (35)), and the pointwise formulation (32) of the plastic flow rule. We aim at providing a continuous dependence estimate for weak energy solutions in terms of the initial and external data, in particular obtaining their uniqueness. To this end we shall have to introduce a further, quite strong simplification. Namely, we shall assume that the plastic dissipation potential H neither depends on the temperature, nor on the damage variable, and we thus restrict to a functional 1-positively homog., and fulfilling (17a). (51) Indeed, while the dependence of H on the fixed temperature profile could be kept, proving continuous dependence/uniqueness results in the case of a state-dependent dissipation potential is definitely more arduous (cf. e.g. [7], [29] for some results in the context of abstract hysteresis/rate-independent systems), and outside the scope of the present contribution. Finally, for technical reasons we will have to strengthen the regularity of C and D and require that In this context we have the following result, where we now write (2.(C, D, E)) as a place-holder for ( be two sets of external and initial data for the regularized viscoplastic damage system (3b, 3c, 3d), with boundary conditions (4) and with given temperature profiles (20). Let (u i , e i , z i , p i ), i = 1, 2, be corresponding weak energy solutions to the initial-boundary value problem for system (3b, 3c, 3d).
Then, there exists a positive constant C P depending on P such that In particular, the initial boundary value problem for the regularized viscoplastic damage system with prescribed temperature admits a unique solution.

Time discretization of the thermoviscoelastoplastic damage system(s).
In all of the results of this section, we will tacitly assume all of the conditions listed in Section 2.1.
3.1. The time discrete scheme. We will consider a unified discretization scheme for both the regularized thermoviscoelastoplastic damage system (3, 4) and for system (1,2). Therefore, within this section, the parameter ν modulating the viscous regularizing contribution to (3c) shall be considered as ν ≥ 0. Given a partition of [0, T ] with constant time-step τ > 0 and nodes t k τ := kτ , k = 0, . . . , K τ , we approximate the data F , f , G, and g by local means: for all k = 1, . . . , K τ . From the terms F k τ and f k τ one then defines the elements L k τ , which are the local-mean approximations of L. Hereafter, given elements (v k τ ) k=1,...,Kτ , we will use the notation We construct discrete solutions to the (regularized) thermoviscoelastoplastic system by recursively solving an elliptic system, cf. the forthcoming Problem 3.1, where the weak formulation of the discrete heat equation features the function space and, for k ∈ {1, . . . , K τ }, the elliptic operator Furthermore, for technical reasons (cf. Remark 2 ahead), we shall add the regularizing term −τ div(|e k τ | γ−2 e k τ ) to the discrete momentum equation, as well as τ |p k τ | γ−2 p k τ to the discrete plastic flow rule, respectively. We will take γ > 4. That is why, we will seek for discrete solutions with e k τ ∈ L γ (Ω; M d×d sym ) and p k τ ∈ L γ (Ω; M d×d D ), giving E(u k τ ) ∈ L γ (Ω; M d×d sym ) by the kinematic admissibility condition and thus, via Korn's inequality (7), u k τ ∈ W 1,γ Dir (Ω; R d ).
Because of these regularizations, it will be necessary to supplement the discrete system with approximate initial data By consistency with the kinematic admissibility condition at time t = 0, we will also approximate the initial datum u 0 with a family (u 0 The data (u 0 τ ) τ may be constructed by a perturbation technique. In connection with the regularization of the discrete momentum balance, we will have to approximate the Dirichlet loading w by a family (w τ ) τ ⊂ W ∩ W 1,1 (0, T ; W 1,γ (Ω; R d )), where we have used the place-holder W : . We will require that w τ → w in W as τ ↓ 0, as well as We will then consider the discrete data For technical reasons related to the proof of Prop. 2 (cf. (63) ahead), it will be expedient to replace the argument of the elasticity tensor C with its positive part. We will proceed in this way in the thermal expansion terms contributing to the momentum balance and to the heat equation. Since we will ultimately prove that the discrete damage solutions are confined to the admissible interval [0, 1], cf. (86) in Prop. 3 ahead, the restriction to the positive part in the argument of C will "disappear" in the end.
Finally, in the discrete version of the damage flow rule (where we will stay with the notation (35)), we will resort to the convex-concave decomposition W (z) = β(z) − λ W 2 |z| 2 from (19), with λ W > 0 and β ∈ C 2 (R + ) convex. For shorter notation, in what follows we will use the place-holders (Ω) to indicate the state spaces for the solutions to system (60) below.
-the kinematic admissibility (u k τ , e k τ , p k τ ) ∈ A(w k τ ) (in the sense of (10)); -the discrete momentum balance a.e. in Ω; (60c) -ϑ k τ ∈ X θ and the discrete heat equation Remark 2. The discrete system (60) has been designed in such a way as to ensure the validity of a discrete total energy inequality, cf. (84) ahead. The latter will be proved by exploiting suitable cancellations of the various terms contributing to (60), as well as the convex-concave decomposition (19) of W in the discrete damage flow rule (60a), where the contribution β (z k τ ) from the convex part has been kept implicit, while the term −λ W z k−1 τ related to the concave part is explicit. The convexity of z → C(z) will also be a key ingredient in the proof of (84), cf. the calculations in the proof of Lemma 3.2.
Several terms in (60) have been kept implicit, not only towards the validity of (84), but also in view of the strict positivity property (61) ahead for the discrete temperature. Our proof of (61) requires that ϑ k τ is implicit in the thermal expansion coupling term on the r.h.s. of the discrete heat equation, cf. the calculations leading to (62), which also rely on the truncation of the elasticity tensor C. Therefore it has to be implicit in the corresponding terms in the discrete momentum balance and in the discrete plastic flow rule, which cannot be thus decoupled one from another. Instead, still compatibly with the proof of (61), the discrete damage flow rule is decoupled from the other equations. This will greatly simplify the proof of existence of solutions to (60).
Because of this implicit character of the thermoviscoplastic subsystem, in order to prove the existence of solutions (to an approximate version of it), we will have to resort to a (nonconstructive) existence result, of fixed point type, for elliptic systems involving pseudo-monotone coercive operators. The regularizing terms −τ div(|e k τ | γ−2 e k τ ) and τ |p k τ | γ−2 p k τ , guaranteeing enhanced integrability properties, have to ensure the coercivity of the pseudo-monotone operator underlying the (approximate versions of the) discrete momentum balance, plastic flow rule, and heat equation. These terms will vanish in the limit as τ ↓ 0. Let us point out that, thanks to them, the right-hand side of the discrete heat equation is indeed an element in H 1 (Ω) * . Thus all the calculations performed in the proof of Proposition 3 and involving suitable tests of the discrete heat equation will be rigorous.
In proving the existence of solutions to Problem 3.1, we will perform these steps: Step 1. While the discrete damage flow rule will be solved by a variational argument, we will approximate the discrete thermoviscoplastic subsystem by truncating the heat conductivity coefficient in the elliptic operator. On the one hand, this will allow us to apply the existence result from [34], based on the theory of pseudomonotone operators, for proving the existence of solutions. On the other hand, due to this truncation we will no longer be able to exploit the growth of κ in order to handle the thermal expansion term on the r.h.s. of (60d). Therefore, in this term we will replace ϑ k τ by its truncation T M (ϑ k τ ). We shall do the same for the corresponding coupling terms in the discrete momentum balance and plastic flow rule.
Step 2. Existence of solutions to the approximate discrete thermoviscoplastic subsystem.
Step 3. A priori estimates on the discrete solutions, uniform with respect to the truncation parameter M Step 4. Limit passage as M → ∞.
In completing Steps 2-4 we will most often have to adapt analogous arguments developed in [31,32], to which we will refer for all details.
Step 1. The approximate discrete system will feature the truncation operator where we suppose that M ∈ N \ {0}. We thus introduce the truncated heat conductivity and, accordingly, the approximate elliptic operator We are now in the position to introduce the approximate discrete system (68). For notational simplicity, we will omit to indicate the dependence of the solution quintuple on the index M . ) ∈ X and z k τ ∈ H s (Ω), find (u k τ , e k τ , p k τ , ϑ k τ ) ∈ B fulfilling -the kinematic admissibility (u k τ , e k τ , p k τ ) ∈ A(w k τ ); -the approximate discrete momentum balance a.e. in Ω; (68b) -the approximate discrete heat equation (68c) Step 2. Existence of solutions to system (68): we have the following result.
Proof. The positivity property (61) (note that the constant providing the lower bound for ϑ k τ is independent of the truncation parameter M ), follows from the analogue of estimate (62), with the same comparison argument. Let us now address the existence of solutions.
First of all, we find a solution z k τ to (60a) via the minimum problem With the direct method in the calculus of variations, it is easy to check that this minimum problem, whose Euler-Lagrange equation is (60a), has a solution z k τ . Let us now briefly address the solvability of the approximate discrete thermoviscoplastic system (68), for fixed k ∈ {1, . . . , K τ } with z k τ given. To this end, we reformulate system (68) in the form with the dissipation potential Ψ k : V → [0, +∞) defined by Ψ k (ũ, p, ϑ) = Ψ k (p) := H(z k τ , ϑ k−1 τ ; p−p k−1 τ ), the space V := W 1,γ Dir (Ω; R d ) × L γ (Ω; M d×d D ) × H 1 (Ω), and the operator A k : V → V * given component-wise by where −div Dir is defined by (14), while the vector B k ∈ B * on the right-hand side of (69) has components The arguments therefore reduce to proving the existence of a solution to the abstract subdifferential inclusion (69 The calculations for [32,Lemma 3.4] show the key role of the regularizing terms −τ div(|e k τ | γ−2 e k τ ) and τ |p k τ | γ−2 p k τ , added to the discrete momentum equation and plastic flow rule, for proving the above estimate.
Step 3. A priori estimates on the solutions of system (68): in order to pass to the limit as M → +∞ in Problem 3.2, for fixed k ∈ {1, . . . , K τ } and z k τ solving the discrete damage flow rule (60a), we need to establish suitable a priori estimates on a family (u k M,τ , e k M,τ , p k M,τ , ϑ k M,τ ) M of solutions to system (68). Along the footsteps of [32], we shall derive them from a discrete version of the total energy inequality (43), cf. (73) below, featuring the discrete free energy (recall that the energy functional E was defined in (21)) E τ (ϑ, e, z, p) := E(ϑ, e, z) + τ γ Ω (|e| γ + |p| γ ) dx .

Proof. Inequality (73) follows by testing (60a) by
), (68b) by p k M,τ − p k−1 M,τ , and by multiplying (68c) by τ and integrating it in space. We add the resulting relations. We develop the following estimates for the terms arising from the test of the momentum balance (68a) where have exploited the kinematic admissibility condition in (75b) and (75c), as well as elementary convexity inequalities to establish estimates (75a), (75c), and (75d). As for the terms arising from the test of the discrete damage flow rule (68b), we observe that again by convexity arguments, also relying on (2.(C, D)3). In particular, in (76d) we have exploited the convex-concave decomposition (19) of W . We now observe several cancellations. Indeed, the two terms on the right-hand side of (75b) respectively cancel with the second term on the r.h.s. of the heat equation (68c), multiplied by τ , and with the analogous term deriving from (68b), tested by p k M,τ − p k−1 τ . As for (75c), the second term on its r.h.s. cancels out with the first summand on the r.h.s. of (76e); the third term on its r.h.s. cancels with the one deriving from (68b), and so does the third term on the r.h.s. of (75d). Also the terms on the r.h.s. of (76a) and (76b) cancel out with the right-hand side of (68c) multiplied by τ : in particular, observe that In fact, with the exception of τ G k τ , all the terms on the r.h.s. of (68c) cancel out. Thus, straightforward calculations lead to (73).
We refer to the proofs of [32, Lemma 3.5] and [31,Lemma 4.4] for all the detailed calculations leading to estimates (74): let us only mention that one has to first test the discrete heat equation (60d) by T M (ϑ k M,τ ) and then by ϑ k M,τ , and exploit the growth properties of κ.
Step 4. Limit passage as M → +∞: We again refer to [32] (cf. Lemma 3.6 therein) for the proof of the following result on the limit passage in the approximate discrete thermoviscoplastic subsystem (68). Let us only mention here that the strong convergences (77a)-(77c) arise from standard lim sup-arguments, developed by testing the discrete momentum balance by u k M,τ − w k τ and the discrete plastic flow rule by p k M,τ .
With this result, we conclude the proof of Proposition 2.

4.
A priori estimates. Again, throughout this section we will tacitly assume all of the conditions listed in Section 2.1. We start by fixing some notation for the approximate solutions.
The main result of this section collects all the a priori estimates on the approximate solutions.
Step 4. Third a priori estimate. We consider the mechanical energy inequality (85) written for s = 0. Since most of the terms on its right-hand side can be handled by the very same calculations developed for the right-hand side terms of (84), we refer to the proof of [32,Prop. 4.3] for all details and only mention how to estimate the third and the fourth integral terms. We use that 0 Ω |ė τ | 2 dx dr + C δ ϑ τ 2 L 2 (0,T ;L 2 (Ω)) via Young's inequality, the previously obtained bound (93), and the fact that |E(x)| ≤ C, and that with the constant δ > 0 chosen in such a way as to absorb the terms |ė τ | 2 and |ż τ | 2 into the left-hand side of (85). Since ϑ τ L 2 (Q) , ϑ τ L 2 (Q) ≤ C thanks to (87o), we conclude a uniform bound for all the terms on the right-hand side of (85). Therefore, estimates (87e), (87h), (87i), and (87l) ensue. We then obtain (87b)(1) and (87c)(1) via kinematic admissibility. Furthermore, (87l) clearly implies estimate (87k), and then (87a) again by kinematic admissibility.
Step 5. Fourth a priori estimate. Let us now shortly sketch the argument for (87r). The very same calculations as in the proof of [31,Prop. 4.10] lead us to deduce, from the discrete entropy inequality (83) written on the generic sub-interval [s i−1 , s i ] of a partition 0 = s 0 < s 1 < . . . < s J = T of [0, T ], the following estimate for the total variation of log(ϑ τ ): for all ϕ ∈ W 1,d+ (Ω), with > 0 arbitrary. Here we have used the place-holder The second, third, and fourth terms on the r.h.s. of (96) can be estimated by adapting the computations from the proof of [31,Prop. 4.10], taking into account the previously obtained bounds. All in all, from (96) we obtain that for every ϕ ∈ W 1,d+ (Ω) with ϕ W 1,d+ (Ω) ≤ 1, where for the last estimate we have used the bound for log(ϑ τ ) in L ∞ (0, T ; L 1 (Ω)) from (87p). Thus, (87r) ensues.

5.
Proofs of Theorems 2.4 and 2.5. We start by fixing the compactness properties of a family (u τ k , u τ k , u τ k , e τ k , e τ k , e τ k , z τ k , z τ k , p τ k , p τ k , ϑ τ k , ϑ τ k , ϑ τ k , ω τ k , ζ τ k ) k , of approximate solutions in the following result, where we again tacitly assume the validity of the conditions from Sec. 2.1. We will only distinguish the case where we only require (2.κ 1 ), from that where (2.κ 2 ) is also imposed and we are thus in the position to enhance the convergence properties of the temperature variables by virtue of the additional estimate (87s).
We refer to [31] and [32] (for a slight refinement of) the proof of the following compactness result.
and suppose in addition that where, for given ∈ B([0, T ]; Y * ) and ϕ ∈ Y we set Then, there exist a (not relabeled) subsequence ( k ) k and a function ∈ L p (0, T ; Furthermore, for almost all t ∈ (0, T ) and any sequence (t k ) k ⊂ [0, T ] with t k → t there holds For expository reasons, in developing the proofs of our existence results we will reverse the order with which we have presented them. More precisely, we will start with the proof of Theorem 2.5 and develop the existence of weak energy solutions and, in addition, establish the validity of the entropy inequality. Let us consider a null sequence (τ k ) k and, correspondingly, a sequence (u τ k , u τ k , u τ k , e τ k , e τ k , e τ k , z τ k , z τ k , p τ k , p τ k , ϑ τ k , ϑ τ k , ϑ τ k , ω τ k , ζ τ k ) k , of solutions to the approximate thermoviscoelastoplastic damage system (79) along which convergences (99) to a seventuple (u, e, z, p, ϑ, ω, ζ) hold. Exploiting them we will pass to the limit in the time-discrete version of the momentum balance. In order to take the limit of the damage and plastic flow rules, and of the temperature equation, we will need to enhance the convergence properties of the approximate solutions. We will do so by establishing the validity of the mechanical energy balance. Therefrom we will strengthen some weak convergence properties, derived by compactness, via standard "limsup"-arguments. We will thus show that the quintuple (u, e, z, p, ϑ) is a weak energy solution of the (initial-boundary value problem for the) regularized thermoviscoelastoplastic damage system. We will also deduce the validity of the entropy inequality.
Step 0. ad the initial conditions (28) and the kinematic admissibility (29). We use the uniform convergences (99c), (99i), (99n), and (99r), as well as the pointwise convergences (99d), (101b), to pass to the limit in the discrete initial conditions (59), also taking into account convergences (57) for the approximate initial data (e 0 τ k , p 0 τ k ) k . Also exploiting (78d) for (w τ k ) k , we pass to the limit in the discrete version of the kinematic admissibility condition. We thus conclude (29).
We pass to the limit of the first integral term on the left-hand side of (83) relying on convergence (99u) for log(ϑ τ k ). For the second integral term, we prove that withδ = α µ and α ∈ [0∨(2−µ), 1), by repeating the very same arguments from the proofs of [31,Thm. 1] and [32,Thm. 1], to which we refer the reader for all details.
To conclude the plastic flow rule (32), it thus remains to show that ζ ∈ ∂ṗH(x, z, ϑ;ṗ) a.e. in Q, H(z(t), ϑ(t); η(t)) dt for all η ∈ L 2 (Q; M d×d D ), (115a) featuring the plastic dissipation potential H from (18). With this aim, we will pass to the limit in the inequalities satisfied at level k using the discrete subdifferential previously established lim inf-inequality (118). We use that to deduce from (99i) that In the same way, we prove that Therefore, again via the Ioffe theorem we have that lim inf k→∞ Q(e τ k (t), z τ k (t)) ≥ Q(e(t), z(t)); by the lower semicontinuity of the functional G we also infer that lim inf k→∞ G(z τ k (t)) ≥ G(z(t)) (where the functionals Q and G are from (21)).

(124)
Let us mention that the very same convergence arguments as in the above lines also give the upper total energy inequality on (0, t).
Step 5. ad the one-sided variational inequality (39a). Using the 1-homogeneity of R (cf. again Sec. 2.2), we reformulate the discrete damage flow rule as the system We now pass to the limit in (125a), integrated along the interval (0, t), using convergences (99k)-(99n), as well as (99w). Let us only comment on the fact that, since β ∈ C 2 (R + ) and the functions (z τ k ) k , taking values in the interval [ζ * , 1] by (86), converge to z uniformly in Q, there holds Moreover, by the Ioffe Theorem (recall that ζ ≤ 0 in Ω and the positivity property (11b) of C ), we have We have thus established for all ζ ∈ H s − (Ω). A localization argument yields the one-sided variational inequality (39a).
We now consider (115a), with the test function η = χ (0,t)ṗ (and χ (0,t) the characteristic function of (0, t)), and deduce that Finally, we test the momentum balance byu −ẇ, integrate on (0, t), and add the resulting relation with (127) and (128). We observe that some terms cancel out and repeat the very same calculations as those leading to (23). We thus conclude that the converse of the mechanical energy inequality (41) holds. This establishes the mechanical energy balance (23) on the interval (0, t). Therefore, all inequalities hold as equalities, the lower and upper limits coincide. Moreover, taking into account the lim inf-inequalities previously observed in Step 3, with a standard argument we conclude that each of the terms on the left-hand side of (85) does converge to its analogue on the left-hand side of (23). Thus, we have in particular lim k→∞ tτ (t) 0ν a s (ż τ k (r),ż τ k (r)) dr = t 0ν a s (ż(r),ż(r)) dr, (129e) lim k→∞ tτ (t) 0 H(z τ k (r), ϑ τ k (r);ṗ τ k (r)) dr = t 0 H(z(r), ϑ(r);ṗ(r)) dr , In particular, from (129b), repeating the very same arguments from the proof of [27,Thm. 5.3], which substantially rely on the uniform positive definiteness of the tensor D, we infer thatė Then, e τ k → e in C 0 ([0, T ]; L 2 (Ω; M d×d sym )). In view of (120) we then infer that e τ k , e τ k → e in L ∞ (0, T ; L 2 (Ω; M d×d sym )) .
Let us also record that (129d) & (129e) yielḋ which gives, by dominated convergence, Finally, (129g) givesṗ from which we deduce, taking into account convergences (99n), (99w), the continuity of H, and repeating the same arguments leading to (117), that Furthermore, in view of the strong convergence (134), the lim sup-inequality (119) immediately follows. This establishes the plastic flow rule.
Step 7. ad the damage flow rule. It follows from (2.(C, D)2) and (99n) that Combining this with (131) we find that Combining this with convergences (99k)-(99o), (99w), and (126), we are thus able to pass to the limit in an integrated version of the discrete damage flow rule (79b), with test functions in L ∞ (0, T ; H s (Ω)) (which embeds continuously into L ∞ (Q)).
Proof of Theorem 2.4 (Entropic solutions). It is immediate to check that Steps 0-5 of the proof of Thm. 2.5 do not rely on the more stringent growth condition (2.κ 2 ) on κ. Therefore, they carry over in the setting of Theorem 2.4. We thus immediately establish the validity of the kinematic admissibility, of the weak momentum balance, of the one-sided damage variational inequality (39a), of the entropy inequality, and of the upper mechanical energy estimate (41).
From (115a) and (113) we also infer the plastic variational inequality (40a). The upper total energy estimate (43) follows from passing to the limit in (84) with the very same arguments used for obtaining the mechanical energy inequality.
This concludes the proof.

Remark 3.
A few more comments on the proof of Theorem 2.4 are in order. 1.
Step 6 in the proof of Thm. 2.5 (and, consequently, Steps 7 & 8) does not carry over to the setting of Thm. 2.4. Indeed, testing the one-sided variational inequality for damage (39a) byż is no longer possible, as, setting ν = 0 we have lost the information thatż ∈ H s (Ω). Hence, all the lim sup-arguments from Step 7, which strengthened the convergences of the approximate solutions, cannot be repeated. 2. As already hinted at the end of Sec. 2.4, it would be possible to establish the existence of entropic solutions to the (non-regularized) thermoviscoelastoplastic damage system by passing to the limit as ν ↓ 0 in its regularized version, featuring a family (κ ν ) ν of heat conductivity functions converging to some κ that only complies with (2.κ 1 ). In fact, the a priori estimates (87a)-(87r), performed on the time-discrete version of the regularized system, are inherited by the weak energy solutions: in particular, it is easy to check that the bounded variation estimate (87r) is preserved by lower semicontinuity arguments. Then, a close perusal of Steps 0-5 reveals that all the arguments performed for the time-discrete to continuous limit carry over to the limit passage ν ↓ 0.
6. Proof of Proposition 1. Throughout this section we will work under the strongly simplifying condition that H neither depends on the damage variable nor on the temperature. Let (u i , e i , z i , p i ), i = 1, 2, be two (weak energy) solutions to the regularized viscolelastoplastic damage system with a given temperature profile Θ, supplemented with initial data (u 0 i , e 0 i , z 0 i , p 0 i ) and external data (L i , w i ), i = 1, 2 (where L i subsume the volume forces F i and applied tractions f i ). We will use the place-holders u := u 1 − u 2 ,ẽ := e 1 − e 2 ,z := z 1 − z 2 ,p := p 1 − p 2 ,σ := σ 1 − σ 2 , where σ i := D(z i )ė 1 + C(z i )e i − Θ i C(z i )E for i = 1, 2. Throughout the proof, we will also often use the short hand · L p for · L p (Ω) , · L p (Ω;R d ) and so forth.
In order to prove (52), we start by subtracting the weak momentum balance (13) for u 2 from that for u 1 , test the resulting relation by ∂ tũ −∂ t (w 1 −w 2 ), and integrate on an arbitrary time interval (0, t). Elementary calculations (cf. (22)) lead to : E(∂ t (w 1 −w 2 )) dx dr + ρ We estimate where the bound for I 3 follows from (140), and the constant η > 0 therein, on which C η > 0 depends, will be chosen suitably. As for the second term on the left-hand side of (141), it follows from the kinematic admissibility condition that : ∂ tp dx dr . = I 5 + I 6 .
We mention in advance that I 6 will cancel out with a term arising from the test of the plastic flow rules. As for I 5 , we have that (cf. (140)) where the constant M 2 depends on z 1 L ∞ , z 2 L ∞ , and e 2 L 2 . The positive constant (note that C( , M 2 ) depends on and M 2 ), will be specified later. Next, we subtract the subdifferential inclusion (35) for z 2 from that for z 1 , test the resulting relation by ∂ tz , and integrate on (0, t). We thus obtain (W (z 2 )−W (z 1 )) ∂ tz dx dr + t 0 Ω 1 2 C (z 2 )|e 2 | 2 − 1 2 C (z 1 )|e 1 | 2 ∂ tz dx dr . = I 7 + I 8 + I 9 + I 10 , with ω i (t) ∈ ∂R(ż i (t)) for almost all t ∈ (0, T ). Observe that I 7 ≤ z 0 1 −z 0 The constant M 3 depends on z 1 L ∞ , z 2 L ∞ , and on the constants ζ 1 * , ζ 2 * > 0 such that for i = 1, 2 there holds ζ i * ≤ z i (x, t) ≤ 1 for every (x, t) ∈ Q, observing that the restriction of W to [min{ζ 1 * , ζ 2 * }, 1] is of class C 2 . The constant M 4 depends on z 1 L ∞ , z 2 L ∞ , e 1 L 2 , and e 2 L 2 . The positive constant λ will be specified later; C(λ, M 4 ) depends on λ and M 4 .
Finally, we subtract the pointwise plastic flow rule (32) for p 2 from that for p 1 , test the resulting relation by ∂ tp , and integrate in time. This leads to with ζ i ∈ ∂H(ṗ i ) a.e. in Q fulfilling the plastic flow rules for i = 1, 2. By monotonicity (it is crucial that H does not depend on the other variables), we clearly have I 10 ≤ 0, while the second integral coincides with I 6 .
In the end, we sum up (141), (142), and (143). Taking into account all of the previous calculations, the cancellation of the last term on the right-hand side of (143) with that arising from the test of the momentum balance, and choosing the constants η, , and λ in such a way as to absorb the term Q |∂ tẽ | 2 into its analogue on the left-hand side, and to absorb t 0 ∂ tz 2 L ∞ into t 0 a s (∂ tz , ∂ tz ) dr, we obtain where we have also used the continuous embedding H s (Ω) ⊂ L ∞ (Ω). Hence, applying the Gronwall Lemma we obtain a continuous dependence estimate in terms of the norms ∂ tũ L ∞ (0,t;L 2 (Ω;R d )) , ∂ tẽ L 2 (0,t;L 2 (Ω;M d×d sym )) , ∂ tz L 2 (0,t;H s (Ω)) , z L ∞ (0,t;H s (Ω)) , ∂ tp L 2 (0,t;L 2 (Ω;M d×d sym )) . Then, estimate (52) follows by easy calculations. This concludes the proof.