Plasticity-Induced Heating: Revisiting the Energy-Based Variational Model

Temperature evolution during plastic deformation is of great importance for the design of manufacturing processes, as well as for the analysis and prediction of tool wear. However, the results from experimental- and numerical-type research are still often contradictory. In this paper, we analyze methods for estimating plasticity-induced heating directly from displacement fields that can be recorded during experiments or extracted from simulation results. In terms of computational methodology, the thermodynamically motivated energy-based variational formulation of the coupled thermo-mechanical boundary-value problem is adapted to the problem at hand. Since an analysis of this variational formulation exhibits challenges and distinct inconsistencies with respect to the problem at hand, an alternative approach is proposed. This alternative approach is essentially a purely thermal finite element simulation, and it is conducted using a heat source term that is empirically based on the fraction of irreversible deformation work converted to heat. Our approach estimates plasticity-induced heating based on the strain and strain rate data derived from displacement fields. We therefore incorporate thermo-visco-plastic constitutive behavior (Johnson–Cook) with a thermodynamically motivated model that specifies the fraction of plastic work converted to heat (the Taylor–Quinney coefficient).


Introduction
When a metallic material is subjected to an irreversible deformation, only a part of the deformation work is dissipated as heat [1,2].The rest is stored as potential energy in the crystal lattice structure of the material.Since the storage of potential energy depends on various state variables, of which temperature is an obvious and significant influence, it is common to distinguish between cold and hot work.In the case of hot work, the temperature has to be above a certain limit, which is specific to each material, so that very little or no potential energy is stored.For the process of cold work, on the other hand, this the storage of energy is characteristic [3].Thus, it is quite common to refer to the stored potential energy as the stored energy of cold work [4].However, in the derivation of the variational model in Section 4, no explicit distinction is made between hot and cold work.Instead, a single temperature-dependent term for the stored potential energy is introduced.As a side note, thermoelastic phenomena also affect the temperature evolution during deformation.Due to the Joule-Thomson effect, elastic deformation can result in either a heating or cooling of the workpiece [5].In the present work, the aforementioned phenomenon or degradation is not considered.Experimental studies have described the influence of elastic deformation on the temperature evolution as negligible, and this assumption is supported in the results of [3,4,6].
On a microscopic scale, the storage of potential energy is a result of the creation and rearrangement of crystal imperfections, namely dislocations, point defects, stacking faults and twins [4].These crystal imperfections cause a local distortion of the lattice resulting in a microstressed state and thus store energy within the lattice structure of the material [3].The dissipation during the process of irreversible deformation is mainly caused by the generation, movement and annihilation of point defects and dislocations [7,8].Which of them, the point defects or the dislocations, dominates the heat generation depends on the material, especially on the structure of the crystal lattice, the temperature and the deformation state [3].While the microscopic processes are revealing and allow a better understanding of the overall process, they will not be elaborated on as this work advocates a macroscopic approach.

State-of-the-Art Works
As introduced in Section 1, only a fraction of the work on irreversible deformation, which is more commonly referred to as plastic work, is converted into heat.What fraction of the plastic work has actually dissipated is a crucial assumption for calculating the temperature evolution, for example.Commonly referred to as β and called the Taylor-Quinney coefficient, the calculations required to determine the fraction of plastic work that is converted to heat can be traced back to Taylor and Quinney [2].In their pioneering work, they found that β is typically between 0.8 and 1.0.
It has been common to assume that the coefficient β is a constant, e.g., as in [9][10][11][12][13].However, most of this information was obtained from quasi-static experiments, which also showed significant discrepancies in the data for nominally identical materials [4].Experimental investigations, e.g., as in [6,[14][15][16][17], also clearly question the limited accuracy of the assumption.Thus, the assumption that the Taylor-Quinney coefficient (TQC) is a constant is an "approximation of dubious validity" [4].
Leaving this assumption aside, there are two general types of approaches to modeling β.These are either approaches developed from a macroscopic point of view, e.g., as in [4,[18][19][20] or a microscopic point of view based on the above-mentioned generation and evolution of crystal defects, e.g., as in [21][22][23].However, even the more prominent ones, such as the model derived by Rosakis et al. [4], have some significant flaws.The model created by Rosakis et al. [4] is discussed in more detail in Section 4.5.
The energy-based variational model introduced by [24], which will be the guiding framework throughout this work, refrains from explicitly modeling the TQC.Within this framework, the dissipation fraction follows directly from the definition of the state potential and a dissipation pseudo-potential [24].This model is discussed in detail in Section 4. Strictly following mathematical conventions, the term 'coefficient' refers only to fixed scalar parameters, i.e., the cases where β is assumed to be constant.For simplicity, all approaches to modeling the fraction of plastic work converted to heat as a function will also be referred to as TQC throughout this paper.
For the sake of brevity, we should also distinguish between two factors, β int and β diff , as introduced by Rittel [25].The first factor β int is defined as follows where ρ is the mass density, c p is the specific heat capacity, T is the temperature, and A p is the irreversible work of deformation per unit of undeformed volume.Note that the irreversible work of deformation is more commonly referred to as W p , which has already been used to describe the stored plastic energy in this work.The factor β int is exactly the factor β introduced by Taylor and Quinney [2].The second factor β diff is defined as and it is used for inclusion in the coupled heat equation, which is discussed in more detail in Section 4.5.In the illustrative plot of thermal work as a function of mechanical work, β int is the secant slope and β diff is the tangent slope [25].
Both versions have been widely used in the literature, e.g., β diff by Mason et al. [16] and β int by Kapoor and Nemat-Nasser [26], to name two prominent examples.However, the specific implications are rarely explicitly explained or emphasized.
In the remainder of this paper, the distinction between β int and β diff will be maintained.What exactly may be called the Taylor-Quinney coefficient will be clarified in a more detailed examination in Section 4.5.Furthermore, a convention will be followed that no precise differentiation is required when the term Taylor-Quinney coefficient is used without further specification.In this context, the term TQC simply refers to a factor that is used to calculate the temperature evolution by multiplying it by the corresponding work or power.

Objective
There is no general agreement on what proportion of the irreversible deformation work is dissipated in heat.This fraction, commonly referred to as the Taylor-Quinney coefficient, has been the subject of research since its introduction by Taylor and Quinney [2].
This lack of agreement is what makes the energy-based variational formulation of the coupled thermo-mechanical boundary value problem proposed by Yang et al. [24] so attractive.They claim that this variational formulation is able to model the constitutive behavior in a thermodynamically consistent way without an explicit assumption for the modeling of the Taylor-Quinney coefficient [27].Thus, this model is the starting point for all further considerations within this work (see Section 4 for more detail).
The aim of this work is to investigate the energy-based variational model in the context of in situ measured displacement field data, which are readily available for a wide range of experiments and processes [28,29].This becomes even more important with the increasing accuracy and precision of full-field optical measurement techniques [30][31][32][33].The first goal of this work was to adapt the energy-based variational model to the special case of given displacement fields.This was followed by a corresponding verification of the validity of this model.Based on these results (see Section 5 for more detail), an alternative approach is presented in Section 6. Section 7 contains a summary and links the main findings of the previous sections.An outlook in Section 8 concludes the work.

The Energy-Based Variational Model
In general, calculus is a branch of mathematical analysis that deals with the real functions of functions, called functionals.Functionals are often expressed as integrals of functions and their derivatives.The goal is to find stationary functions that, apart from singularities, lead to the maxima, minima, or saddle points of the functional.Some classical problems in mathematics and physics can be solved elegantly with the help of functionals.
In this way, Yang et al. [24] derived a thermodynamically consistent variational formulation of the coupled thermo-mechanical boundary value problem for general dissipative solids, which is simply referred to in this work as the energy-based variational model (EBVM).In this section, the mathematical basis for the EBVM is established, starting from Yang et al. [24] and then accompanied by additional insights from Stainier and Ortiz [27], Stainier [34], and Su and Stainier [35].Further, our own explanatory intermediate steps, assumptions, and modifications are also included.Separating them into a separate chapter would lead to unreasonable complication and would severely impair the overall comprehensibility.This section is divided into five subsections, of which Section 4.1 illustrates the approach for general dissipative solids and Section 4.2 treats the adaptation to the case of finite thermo-visco-plasticity. Subsequently, the Johnson-Cook plasticity model is incorporated into the EBVM framework (see Section 4.3 for more detail).In Section 4.4 the EBVM is further adapted to handle the given displacement fields.The introduction and adaptation is concluded with a more detailed examination of the TQC in Section 4.5.

Derivation for General Dissipative Solids
Fundamental to this framework is the assumption of a Helmholtz free energy where W is the Helmholtz free energy per unit of undeformed volume, U is the internal energy per unit of undeformed volume, T is the absolute temperature, and N is the entropy per unit of undeformed volume [24,36] (pp.112-114).
All motions and thermodynamic processes were performed by a continuum body with the reference configuration Ω ⊂ R dim , dim ∈ {2, 3}.The motions were described by a timedependent deformation mapping φ : [a, b] is the corresponding time interval.The deformation gradient is introduced as F = ∇φ.Furthermore, a collection Z of general internal variables was assumed, the specific form of which depended on the material class and could not be defined universally for all solids.This led to the reformulation of Equation (3) (as in [24]): Applying the Legendre-Fenchel transformation to the internal energy U leads to [34] U(F, N, Z) = sup in terms of which the equilibrium relations can be stated as Therein, P e are the equilibrium stresses and Y is the thermodynamic driving forces that conjugate with the internal variables Z [24].Furthermore, the first Piola-Kirchhoff stress tensor P can be split as where P v is the viscous stress.To ensure thermodynamic consistency, the first principle of thermodynamics in its local form must hold.Therein, R is the mass density per unit of undeformed volume, Q the distributed heat source per unit mass, and H the outward heat flux.With the theorem of Coleman and Noll [37], T = ∂ N U, and Equations ( 7)-( 9), this can be restated as [24] T Thus, it becomes clear that the equilibrium relations must be supplemented with an appropriate kinetic relation that allows the determination of P v , Ż, and H in order to obtain a closed set of governing equations that define a well-posed initial boundary-value problem [24].A general kinetic potential is a function ∆( Ḟ, Ż, G, F, N, Z) for that [24,27] For example, a material with uncoupled viscosity, rate sensitivity and heat conduction has a general kinetic potential of the form where ϕ * , τ * and χ are the viscous potential, the kinetic potential and a Fourier potential, respectively [24].The superscript (•) * refers to the use of the dual potential, the implications of which will be explained in Section 4.2.Regarding the heat conduction, the Fourier potential χ(G), commonly known as Biot's dissipation function, was assumed such that [38] with With the above equations, the constitutive equations of continuum mechanics can be restated in a thermodynamically consistent, variational framework, as was explained Ortiz and Stainier [39] or Yang et al. [24].This eventually leads to a power density function where ϑ(F, N, Z) = ∂ N U is the equilibrium temperature in contrast to the external temperature field T. This distinction follows from the variational formulation.The optimality condition with respect to Ṅ yields T = ∂ N U = ϑ [24,27].It is important to note that this power density function does not include inertial effects or distributed heat sources.
Considering the problem at hand, namely the estimation of the temperature field evolution from displacement field data, the neglect of inertial effects seems acceptable at first.As for heat sources, of which there were none anyway, the heating of the material was due to the dissipation of plastic work.With respect to Equation (4), resulting in U = Ẇ + ṄT + N Ṫ, Equation ( 18) is restated as When the above is expressed in a time-incremental framework and with integration over the whole domain, it yields the functional [24,35] with Again, it is important to note that this equation does not include the boundary integrals.Therefore, it only allows an implicit inclusion of the Dirichlet boundaries or an adiabatic Neumann boundary.Although the extension to general Neumann or mixed boundaries is rather straightforward [35], there will be no need for their incorporation in this work due to the examination of the adiabatic case in Section 5.
With this functional, the incremental coupled thermo-mechanical boundary value problem can be stated as a variational principle [24,35] [φ n+1 , T n+1 , Z n+1 ] = arg inf It should be noted that, unlike [24], the variation used in this paper does not involve minimization with respect to N n+1 , which is a consequence of reformulating the Equation (18) in the manner of Equation (19).In addition, Φ n does not depend on the gradient of Z n+1 , which allows for a local enforcement of the minimization [24].Restating the functional Φ n as eventually yields This variational formulation has two major advantages over the conventional formulations of thermo-mechanical problems, and both are due to its variational nature.First, is the symmetry associated with all variational formulations; second, is the applicability of tools from optimization theory and mathematical programming [35].Furthermore, this approach is characterized by its extensibility.The incorporation of material behavior such as thermal softening or strain hardening-as well as the physical implications of the process or experiment such as friction, damage, and fracture-can be accomplished by adding or changing the respective potentials.

Adjustment to Thermo-Visco-Plasticity
For metals, it is generally accepted that Helmholtz free energy is composed of elastic, plastic, and thermal energies, as expressed in [35]: When considering the large plastic deformation, it is quite common to neglect the elastic deformation in the numerical simulations, e.g., as in Li and Peng [40] or Zhao et al. [41].As already mentioned and justified in Section 1, this assumption will also be made.This implies not only that W e = 0, but also that F = F p , which is the plastic deformation.The stored heat energy W th is defined as where c p is the specific heat capacity, T r is the reference temperature, and ρ 0 is the mass density.Note that, in Section 4.1, R was used for the mass density per unit of undeformed volume.This replacement emphasizes that from this point on, the common assumption of a constant mass density is made.Equation (26) corresponds to the expression used by Su [42].
It is similar to the one derived by Miehe [43] but reduced by a constant summand, which becomes obsolete anyway when using the incremental framework.
The term for the plastic energy W p is specific to the material model used, and it is defined in Section 4.3.Therefore, the expression for N = ρ 0 η = −∂ T W, with η as the specific entropy, must also be derived specifically for each material model.
The set of internal variables for thermo-visco-plastic materials is defined by Stainier and Ortiz [27] as Z = {F p , εp }, where εp is a scalar internal variable measuring the cumulative or equivalent plastic strain.Similar to the work of Stainier [34], the plastic part of the deformation gradient F p in this work is replaced by another strain measure as an internal variable.However, in contrast to Stainier [34], the true plastic strain was chosen instead of the engineering strain.To the best of the authors' knowledge, this did not affect the validity of the above variational principles as long as the Piola-Kirchhoff stress tensor P was also replaced by the Cauchy stress tensor σ, which is commonly referred to as the true stress.
The equivalent true plastic strain εp and the plastic part of the deformation gradient F p were then connected by ε = 1  2 ln(F T • F) and thus with neglecting the elastic strains ε p = 1 2 ln(F pT • F p ) and the rate relation via the classical von Mises (J 2 ) plasticity, as shown in [35], However, under the assumption that the elastic strains are negligible, ε p is not needed as an internal variable to calculate the flow direction and thus the decomposition of F into F p and F e and ε into ε p and ε e , respectively.As such, the reduction to εp as the only internal variable was acceptable.
Still using the assumption of uncoupled potentials as in Equation ( 15), the general kinetic potential of thermo-visco-plastic materials is defined as ∆ = ψ * + χ, where ψ * is a dual dissipation pseudo-potential [27].This is in agreement with [44], who state that the constitutive behavior of irreversible processes is governed by two potentials: a state potential and a dissipation pseudo-potential.Although it is common to denote the heat conduction potential as part of the dissipation pseudo-potential [45] (pp.61-65), χ will be treated separately for the sake of clarity when considering Furthermore, this work is restricted to isotropic materials with thermal conductivity λ, thus leading to an exemplary Fourier potential in the style of [35].
For the dual dissipation pseudopotential ψ * , the starting point is the traditional approach of deriving the evolution of the internal variable εp from a dissipation pseudopotential ψ(Y) as εp = ∂ Y ψ.Note that Y is still the force thermodynamically conjugate to the internal variable, but with the assumption Z = εp , now a scalar.By means of a Legendre-Fenchel transformation, the dual dissipation pseudo-potential can be defined as or-alternatively and assuming that εp is determined only by the local thermodynamic state, as Note that this is consistent with Equation ( 13).Also note that ψ must be defined convexly in Y.The convexity of ψ(Y) then implies the convexity of ψ * ( εp ), which ensures the positivity of the intrinsic dissipation [27].This positivity is needed to ensure the verification of the Clausius-Duhem inequality, which is commonly referred to as the second principle of thermodynamics [24].
The specific form of ψ * again depends on the material model used and is therefore introduced in Section 4.3.
Regardless of its specific form, there are several ways to integrate ψ * or χ as part of (23).A fairly straightforward approach to integrating χ would be to use Su and Stainier [35] as which is similar but not identical to the use of an implicitly computed right Riemann sum.
For the integration of ψ * , however, a more sophisticated and accurate approach was used, as was suggested by Stainier [46], Here the strain and temperature are calculated as where α ∈ [0, 1] is an algorithmic parameter [46].For other ways of integration, which may not exactly satisfy the second principle of thermodynamics, the interested reader is referred to [24].

Adjustment to the Johnson-Cook Plasticity Model
Because of its simplicity and because it is one of the most commonly used ratedependent models, the Johnson-Cook (JC) model was used throughout this thesis to describe visco-plastic material behavior [47,48].This model was deliberately introduced in a dedicated chapter to highlight its implications and to allow for a fairly straightforward adaptation of the overall approach to other constitutive models.In fact, an additional constitutive model will be treated in Section 5.2.Since this model will be restricted to this subsection only, an introduction in this section is deliberately omitted in order not to introduce unnecessary complexity.The JC model gives a constitutive equation that is suitable for metals subjected to large strains, high strain rates, and high temperatures, such as in [35] σy with the nondimensional temperature θ being defined as Here, σy is the equivalent stress corresponding to the yield stress σ y .The parameters A, B, C, n, and q are material constants; ε0 is a reference strain rate; T m is the melting temperature; and T r is the reference temperature (which is usually room temperature).The term in the first set of parentheses in Equation (33) gives a stress as a function of strain, where B is the strain hardening coefficient and n is the strain hardening exponent.The terms in the second and third sets of parentheses represent the effects of strain rate and temperature, respectively [48].
Based on Equations ( 7), ( 9) and ( 12)-as well as through considering the neglect of an elastic energy W e and the fact that the stored thermal energy W th and the Fourier heat conduction potential χ do not depend on strains or strain rates-it can be stated that To the best of the authors' knowledge, this relation has never been explicitly stated in the EBVM literature.When using the assumed neglect of the elastic stored energy W e , it even contradicted what was stated by Stainier and Ortiz [27].However, it was consistent with the terms for W p and ψ * used by Su and Stainier [35] (see also Su [42] (p.78)).Furthermore, Section 5.2 will show that the stress data that were calculated using Equation (35) were in good overall agreement with the stress data reported by Stainier and Ortiz [27].Regarding its physical interpretation, this division of the stress was considered to be in accordance with the general approach presented by Ziegler and Wehrli [49], who divided the stress into a quasi-conservative and a dissipative part.To further substantiate the similarity, the quasi-conservative stress therein resulted from a partial derivative of the free energy with respect to the strains, and the dissipative part was related to a partial derivative of a dissipation function with respect to the strain rate [49,50].
Assuming that, due to the neglect of elastic strains, the equivalent stress σ coincides with the equivalent yield curve σy , Equation ( 33) leads to a possible description of the plastic stored energy potential as [35] and to the dissipation pseudopotential in the form of [42] Here, as well as according to Su and Stainier [35], the parameters A s , B s , A d , and B d correspond to the stored and dissipative parts of the parameters A and B. Thus, they can be, respectively, expressed as A s + A d = A and B s + B d = B [35].It should be noted that, however, the strain rate-dependent part of the constitutive equation was modeled exclusively in the dissipation pseudopotential, which includes A d and B d as part of A and B, respectively.Thus, the physical implications were not as distinct as the naming suggests.Nevertheless, the relative values of these parameters have a strong influence on the distribution of the plastic work when they are converted into heat.However, their ratio is not identical to the Taylor-Quinney coefficient, as will be explained in Section 4.5.
With respect to the specific entropy introduced in Section 4.2, the above equations lead to

Implications of Given Displacement Fields
The given time-resolved displacement fields also contained information about strain and strain rate fields.Equation ( 27) already stated how to calculate the equivalent plastic strain rate εp .By neglecting the path dependence of the deformation, the equivalent strains could be calculated in a similar manner as However, this neglect only gives a reasonable approximation for time steps ∆t → 0. Furthermore, a strain measure that is independent of the path can lead to negative increments ∆ε p = εp n+1 − εp n , e.g., when the loading direction changes.Considering a rather simple explicit approximation of the temperature increment in the style of Equation (1) as ∆T ≈ β int σ∆ε ρc p ∆t , the effect of such negative increments is particularly evident: in the illustrative example of cyclic loading, this would mean that the specimen first heats up and then cools down when the loading direction is reversed, thus tending toward the initial temperature as it approaches its initial geometric state.Mathematically, in the presented variational model, a negative equivalent strain increment can lead to imaginary components in the potentials or pseudo-potentials, which inhibit the finding of the stationary temperature field.Therefore, and despite being aware of the introduced inaccuracy, the authors recommend replacing ∆ε p ∆t with εp in Equation (31) in cases where there are no path-dependent equivalent strains available.To improve the accuracy of the results, the equivalent strain fields used throughout this paper were calculated with respect to path dependence.For details on the methodology, the interested reader is referred to Hartmann et al. [51].
Since the displacement fields are given and εp is the only internal variable as stated in Section 4.2, the optimization with respect to it in Equations ( 23) and ( 24) becomes obsolete.Therefore, the optimization problem is reduced to with the function Φ n represented as To bypass the optimization and allow for a direct solution of a system of nonlinear equations, the functional Φ n could be derived analytically with respect to T n+1 in future research, thus eventually leading to

Calculation of the Taylor-Quinney Coefficient
Based on the EBVM introduced in the previous subsections, the TQC β and factor β diff can be calculated a posteriori [27].As mentioned in Section 2, β diff was the factor used in the heat equation.It is common to refer to β diff as the differential form of the TQC.Strictly speaking, however, these two factors, β and β diff , must be distinguished in the case of a temperature-dependent W p = W p (ε p , T).As introduced in Section 2, the TQC β int = β is commonly used to describe the fraction of plastic work converted into heat and, for simplicity by maintaining a limited set of variables, is also used to describe the fraction of plastic work rate Ȧp converted into heat rate Qp as with Ȧp = σ εp [4].If W p decreases solely due to a change in temperature, the resulting difference ∆W p also contributes to an increase in temperature.As a consequence, in the case of a complete dissipation of the currently introduced deformation work rate Ȧp with a simultaneous decrease in W p due to temperature, β diff can take values of β diff > 1.
However, the fraction of the plastic work rate Ȧp in the heating rate Qp cannot, by definition, be greater than one, i.e., β ̸ > 1.
It should be noted that neither Rosakis et al. [4] nor Stainier and Ortiz [27], on whose work this section is largely based, explicitly distinguished between β, β int , and β diff , but rather they made a tacit assumption.For the sake of clarity and consistency, however, this distinction in nomenclature will be maintained throughout the rest of this paper.
To derive the corresponding relation for β for the EBVM under consideration, this work follows Stainier and Ortiz [27], who provide equations for the dissipation as D int = Y εp and D int = β diff S • D p , where Y is the scalar thermodynamic driving force, S is the Mandel stress tensor, and D p is the plastic deformation rate.Using the scalar value corresponding to the Cauchy stress tensor σ and consistently neglecting elastic strains, the latter equation can be rewritten as D int = β diff σ εp , or due to the work conjugate σ εp = σ εp as D int = β diff σ εp .Combining the two equations gives β σ εp = Y εp .Rearranging and using Equations ( 13) and (33) in the form Y = ∂ψ * ∂ εp finally yields By including Equation (33) and the derivative of Equation ( 37) with respect to εp , we obtain Here it becomes obvious that β does not depend on the temperature T. Thus the a posteriori calculation of the TQC becomes an a priori calculation in the case of given displacement fields.As a side note, this cancellation of the temperature, when assuming a similar decomposition into W p and ψ * as for the JC model presented, holds for any approach that models the flow stress with a multiplicative incorporation of the temperature or any term exclusively depending on it as a variable, e.g., the model proposed by Klopp et al. [52] or the KHL constitutive models [53,54].However, it is not generally applicable.The equation for β diff is derived as by Stainier and Ortiz [27] where σW p = ∂W p ∂ε p and σψ * = ∂ψ * ∂ εp are the respective parts of the equivalent stress σ.The superscript (•) SO refers to the authors of the respective paper and is used for distinction.It will be shown in Section 5 that the inclusion of this factor β SO diff in the adiabatic heat equation results in the same temperature evolution as the EBVM, at least for when the JC parameter q = 1.It should be mentioned that, from an energy conservation point of view, the case of a non-linear dependence of W p on T with the simultaneous use of a Helmholtz thermal free energy W th leads to an additional term in the heat equation.The derivation can be found in Appendix A.1.For the presented JC version of the EBVM, this additional term vanishes for q = 1, which may explain the discrepancies between the temperature resulting from the EBVM and from Equation (48) for q ̸ = 1.With a rearrangement of Equation ( 44) as β = 1 − ∂W p ∂ε p / σ, the proximity to the approach proposed by Rosakis et al. [4] becomes evident.For the common assumption that the only internal variable on which the stored plastic energy depends is the plastic strain, as well as for treating the one-dimensional case without using equivalent strains and strain rates, they proposed where E is the stored cold work; E ′ is its derivative with respect to ε p ; and σ is a stress function that, for the most general case, can depend on ε p , εp , and T.
In cases where the plastic free energy does not depend on the temperature W p ̸ = W p (T), the term T ∂ σW p ∂T in Equation ( 46) vanishes, thereby resulting in Thus, the most important difference is the exclusive dependence of E on the cumulative plastic strain, which is only an assumption of Rosakis et al. [4], and this is in contrast to the dependence of W p on the cumulative plastic strain, as well as on the temperature, in the EBVM.The second difference is the assumption of the stress function σ, which is also used to model E. This function σ does not coincide with the constitutive equation for the yield stress σ y [4].
In contrast, W p is a result of the partition of σ y (according to Equation ( 33)).The specific form of this split is itself an assumption, since, to the best of the authors' knowledge, the only additional requirements are the convexity of ψ * ( εp ), as was required in Section 4.2, and the fact that-for the sake of thermodynamic consistency-the stored potential energy W p must not depend on the strain rate [4].Thus, the formulation of β diff derived from the EBVM in this paper appears to be more general while requiring fewer assumptions.Despite the improved consistency of the approach, it should still be carefully noted that the assumption of the explicit modeling of β diff is ultimately just replaced by another assumption of the expressions for W p and ψ * .
Due to the lack of an illustrative physical interpretation of the definition of β diff in Equations ( 46) and ( 47), a descriptive approximation β approx diff was introduced as This approximation, despite its reduced accuracy, allows a clear interpretation, i.e., the full amount of plastic deformation work in each incremental step ∆A p that is not stored as plastic potential energy W p is dissipated as heat and thus contributes to an increase in temperature.This approximation also takes into account the exclusively temperaturedependent decrease in W p .

Analysis of the Energy-Based Variational Model
In the following subsections, the EBVM is evaluated by comparing the resulting temperature with the reference temperatures derived from the adiabatic heat equation, i.e., the heat equation that neglects conductive heat transfer, The specific formulation of the β diff factors was introduced as appropriate.The implementation and all corresponding calculations were performed with The Mathworks, Inc. [55].Furthermore, for the sake of clarity, only a single material point equipped with a synthetic data set for strains and strain rates will be considered.The corresponding integral over the domain Ω in Equation ( 41) therefore vanishes.It should be noted that, in general, there is no need to calculate this integral if heat conduction is neglected.When taking the displacement fields as given, the temperatures at the spatial points were coupled only by the heat conduction term χ.Without this conduction, Equation (41) could be rewritten as with the auxiliary function If the integral over Ω is now approximated as a sum over m infinitesimally small elements, the functional Φ n can be written as where V i is the corresponding infinitesimal volume of element i. Due to additivity, the temperature T n+1,i for each element i is optimized to find the stationary point of W i with respect to T n+1,i .In other words, the stationary value of the integrated potential Φ n is the sum of the stationary values of the subpotentials W i multiplied by the corresponding infinitesimal volume V i .Since the values of V i have no influence, and as the total potential Φ n serves only as a vehicle while its value is irrelevant, the temperatures can be obtained by local optimization as long as the conductive coupling is neglected and the strain and strain rate fields are given.The synthetic data set mentioned above consists of equivalent plastic strains εp = [0, 0.4] and a constant equivalent plastic strain rate εp = 1000 s −1 .The calculation uses num steps = 200 time steps and thus, assuming path independence, a corresponding time step size ∆t = max(ε p ) num steps • εp = 1 × 10 −6 s.Furthermore, it is irrelevant whether the equivalent strain and strain rate values, εp and εp , or the true strains and strain rates, ε p and εp , are used in the case of uniaxial loading.This equivalence is derived in Appendix A.2, and it is crucial for the calibration of the parameters.Although irrelevant for the evaluation, it should be noted that the assumption of adiabaticity, i.e., λ = 0 or χ = 0 in the case of a single material point, is indeed a valid approximation for strain rates as high as those used [6].

Examination with the Johnson-Cook Constitutive Model
Although the resulting temperatures will not be discussed in terms of their absolute values, this section will nevertheless use actual JC material parameters for the titanium alloy Ti-6Al-4V, as obtained in a prominent paper by Lee and Lin [48].In addition to the parameters A = 724.7 MPa, B = 683.1 MPa, C = 0.035, n = 0.47, q = 1.0, and ε0 = 1 × 10 −5 s −1 , a mass density ρ 0 = 4430 kg m −3 was used and the specific heat capacity c p was converted from 0.13 cal(g • C) −1 to 518 Jkg −1 K −1 .Since Lee and Lin [48] did not explicitly specify the melting temperature T m -i.e., the actual melting temperature of Ti-6Al-4V-T m = 1604 • C was used [56].
As indicated in Section 4.3, splitting the parameters A and B into A s , A d and B s , B d had a significant effect on the fraction of plastic work converted into heat, although it was found to be neither identical to the TQC β nor to the factor β diff .To illustrate the justification of the EBVM, a split factor was introduced as The actual values of this split and the fact that the split was identical for A and B were arbitrary in nature and only enhanced the comprehensibility of the overall evaluation.
For this evaluation, four different temperature histories were calculated: (1) T (a) , which results from the EBVM by solving Equation ( 40); (2) T (b) , which results from using the factor β SO diff , as in Equation (46), in the adiabatic heat Equation ( 50); (3) T (c) , which results from a similar use of β, as defined in Equations ( 44) and ( 45); (4) T (d) , which is obtained from β approx diff , as in Equation (49).The superscripts (•) (•) are retained for clarity.The exact calculation for the temperature evolution T (b) to T (d) will be explained in the following paragraphs.
The evolution of T (b) was determined based on β SO diff = β SO diff (T (a) , εp , εp ) using T (a) as an input to accurately describe the β diff factor resulting from the EBVM.For each incremental step, the temperature T (b) n+1 was then calculated as The use of σ(b) (T (b) , εp , εp ) instead of σ(a) (T (a) , εp , εp ) may be controversial.Aside from not making a significant difference, this choice was made to correctly describe the current plastic work rate Ȧp as a function of the current stress σ, which-in turn-depends on the current temperature T and could potentially show a discrepancy between T (a) and T (b) .
The choice of which indices, n or n + 1, to use to calculate T n+1 is a matter of the time integration method used.This topic is treated in more detail in Section 6, but suffice it to say here that the calculation as shown above, but using σ(b) (T (b) n+1 , εp n+1 , εp n+1 ) instead, would be the usual implicit or backward Euler method, thus introducing a nonlinearity.Equation ( 54) is therefore only an approximation, with the results being between those for the explicit and implicit Euler schemes.However, due to the small synthetic time steps used, the choice of time integration method was insignificant but the authors considered this choice to be the most illustrative.Recalling Equations ( 1) and ( 2), it becomes apparent that Equation (54) uses β diff instead of β int in combination with ∆A p = ∆t σ εp .This is due to the fact that ∆A p /∆t refers to a time step that is deliberately small enough to approximate the rate ∆A p /∆t ≈ Ȧp .
To determine T (c) , β is calculated a priori from the strain and strain rate data (as in Equation ( 45)), which is followed by a stepwise approximate solution of the adiabatic heat Equation (50).For this approximate solution, a similar time integration method was used as for the determination of T (b) , which finally led to as defined in Equation (49).Like β SO diff , β approx diff depends on the state variables of the EBVM and is therefore calculated a posteriori as The temperature n+1 is then calculated as As for the calculation of T (b) , it is again controversial to replace σ(d) with σ(a) .However, the above argument in favor of σ(d) is similarly valid for this case.Plotting the evolution of the temperature rise ∆T = T − T 0 for the calculated temperatures T (a) to T (d)  over the equivalent strains εp , yields the Figure 1 for β split = 1.It can be clearly observed that all four temperatures were nearly congruent.The scenario where no plastic Helmholtz free energy is stored, i.e., W p = 0, was also investigated by Yang et al. [24], Su [42] and Su and Stainier [35].
Choosing, e.g., β split = 0.5 gives the Figure 2. It can be seen that the temperature rises for T (a) (which resulted from the EBVM) and T (b) (which resulted from a calculation with β SO diff ) were nearly identical, as was already predicted in Section 4.5.However, the temperature evolution of T (c) and T (d) were significantly different.Another interesting feature was the close proximity of these two curves.The choice of β split = 0, as shown in Figure 3, further emphasizes these qualitative observations.This has two important implications.First, the congruence of T (a) and T (b) suggests the correctness of the adjustments made to the original EBVM throughout this work and the implementation, or at least its consistency, since the EBVM yielded the results predicted by Stainier and Ortiz [27] with their definition of β SO diff .Second, the adapted EBVM exhibited incorrectness, as T (a) and T (b) differ significantly from T (d) .The latter was, indeed, only an approximation, but with an undeniable claim to approximate correctness based on the physically motivated derivation of the underlying β approx diff in Equation ( 49).These problems of the EBVM are, by the way, in their general form and are a order of magnitude independent of a myriad of small modifications.A brief overview of the variations examined without significant impact can be found in Appendix A.3.In addition, Appendix A.4 contains the corresponding plots of the β factors, and Appendix A.5 contains a plot for the JC parameter q ̸ = 1, which-in this case-shows the predicted discrepancy between T (a) and T (b) .
To put the term "incorrectness" into perspective, it should be noted that the temperature evolution obtained by the EBVM is not necessarily wrong if the model is calibrated accordingly.However, this would still imply a lack of physical interpretability if the postulated a posteriori determination of the TQC and the factor β diff were not consistent with the calculated stored plastic energy W p .
The above observations together with the theoretical considerations in Section 4.5 suggest that the discrepancies and deficiencies of the EBVM are related to a temperaturedependent W p , i.e., for the EBVM with the JC model and the current split to W p and ψ * if β split ̸ = 1.Despite the congruence of T (a) and T (b) , the errors in the adaptation regarding the case of given displacement fields and neglected elastic strains (or in the implementation) cannot be excluded.Therefore, in the following section, the EBVM will be further evaluated with another constitutive model.

Examination with the Stainier-Ortiz Constitutive Model
This constitutive model was deliberately not introduced in Section 4, since it will be used exclusively in this section and a double assignment of the variables would have introduced an avoidable additional complexity into the matter.The model is taken from a paper by Stainier and Ortiz [27] and will therefore be referred to as the Stainier-Ortiz (SO) model.Since this model is entirely based on their work, it has been deliberately omitted to provide the same reference for each individual equation, see [27].In contrast to the original model, the SO model in this work again omits the elastic strains ε e and thus also the elastic Helmholtz free energy W e .The plastic stored energy W p in this model is composed of a power law term and an exponential saturation term as The dual dissipation pseudopotential is composed of a rate-independent term, i.e., it is homogeneous of order 1 in εp , and a rate-dependent power-law term as The rate-independent term itself is, like W p , composed of a power-law term and an exponential saturation term as Using Equation (35) gives an equation for the yield stress as Entropy can be derived as and the TQC β has the form . ( 63) The critical stresses in the above equations are defined as To ensure good comparability, the nomenclature utilized was identical to that used by Stainier and Ortiz [27] (except for the use of σy instead of σ y to avoid overlapping with the yield stress).The assignment of other variables, e.g., n (which was also used in the JC model), should remain sufficiently clear since the use of the SO model is limited to this section.It was also pointed out that all σ variables should be written as σ when calculating and calibrating their values with equivalent strains and strain rates in the case of a fully consistent approach regarding the nomenclature.Since the true strain σ is interchangeable with the equivalent true strain σ in the case of uniaxial loading (see Appendix A.2), this correction has been omitted.
Table 1 lists the SO parameters that were used by Stainier and Ortiz [27] for the aluminum alloy 2023-T3 with an initial and reference temperature T 0 = 293 K, a specific heat capacity c p = 875.0Jkg −1 K −1 , and a mass density ρ 0 = 2780.0kg m −3 .The reduced number of variables in the Table 1 was due to the rate-independent modeling of the material behavior of the aluminum alloy.
Figure 4 shows the resulting temperature increases ∆T for the temperatures T (a) to T (d) , which was calculated in the same way as in Section 5.1.Again, these curves were found to be nearly congruent, which is consistent with the above assumptions since W p is independent of the temperature for the material parameters that are presented in Table 1.In addition, Figure 4 also shows the temperature increase that were calculated by Stainier and Ortiz [27] as a dotted line.Furthermore, Figure 5 shows the evolution of the β-factors and Figure 6 the evolution of the stress, which were calculated using Equation (61).These two figures also include the evolution of the state variables that was reported by Stainier and Ortiz [27].The discrepancies between the calculated and reported values of the variables may have several causes.As can be seen in particular in Figure 6, the elastic strain range was not excluded but only ignored due to its limited extent and the fact that these graphs were only for illustrative purposes.Furthermore, the exact calculation method was explicitly and deliberately not specified by Stainier and Ortiz [27].This did not only include the time step ∆t or the specific approach of evaluating the integral t n+1 t n ψ * ′ dt ′ , as was the approach in Equation (31), which was introduced by Stainier [46] only one year later.Their results cannot be obtained at all by the incremental approach of the EBVM since the calculation of W p with the given SO material parameters for tantalum, which is also part of their research, would imply an impermissible division by zero.
Nevertheless, the calculated results show a relatively good agreement with the evolution of the state variables extracted from the work of Stainier and Ortiz [27] in all Figures 4-6.This reinforces the assumption of a correct adjustment of the EBVM in Section 4, as well as a correct corresponding implementation.In particular, the good agreement of the calculated stress values seems to validate Equation (35).
However, if the stored plastic energy W p becomes temperature-dependent-e.g., by arbitrarily choosing ω 0 = ω0 = 0.001, as shown in Figures 7 and 8 with the corresponding temperature and β factor evolution)-again a discrepancy between the temperatures T (a) (which results from the EBVM) and T (d) (which results from β approx diff ) becomes evident.Also, a similar proximity of T (c) and T (d) , as already noted in Section 5.1, can be observed.To further substantiate the assumption of a correct adaptation and implementation, Appendix A.6 contains similar plots for temperature, β factors, and stress for α titanium-a that was also studied by Stainier and Ortiz [27].A comparison for tantalum is not possible, as mentioned above, because the incremental version of the EBVM does not work for the listed set of parameters.
In conclusion, the considerations of this section suggest a correct adaptation of the EBVM to given displacement fields, as well as to neglect elastic strains, the corresponding elastic Helmholtz free energy, and a correct implementation.However, these potential sources of error cannot be completely ruled out.More likely, however, is the thermodynamically inconsistent modeling of material behavior by the EBVM in the case of a temperature-dependent plastic stored energy W p .A closed theoretical derivation of the source of the discrepancy was not possible within the scope of this work due to the limited time available.However, these potential sources of error cannot be completely ruled out.More likely, however, is the thermodynamically inconsistent modeling of material behavior by the EBVM in the case of a temperature-dependent plastic stored energy W p .A closed theoretical derivation of the source of the discrepancy was not possible within the scope of this work due to the limited time available.However, this is a crucial step for further research.Possible starting points might be a reconsideration of the fundamental thermodynamic admissibility of a temperature-dependent plastic stored energy W p or the definition of stress as σ = ∂W p ∂ε p + ∂ψ * ∂ εp , as shown in Equation (35).Regarding the stress definition, Ziegler and Wehrli [49], who followed a similar division of the stress into a quasi-conservative part and a dissipative part, used an additional scalar factor multiplied by the derivative of the dissipation function with respect to the strain rate to ensure an orthogonality condition.However, its implications and the transfer to the EBVM under consideration were beyond the scope of this work [49,50].It has also been pointed out that there is a similar, but different, splitting of the stress into viscous and plastic stresses, as proposed by [57].As they also introduced this split for the JC model, this could be a promising starting point for further research.
When considering how to determine the temperature evolution from displacement field data, two fundamentally different approaches can be followed.First, a temperature independent W p = W p (ε p ) can be defined.This implies that, by the way (due to the modified splitting of W p and ψ * , which is still coupled by Equation ( 35)), the term in Equation ( 46) vanishes.However, the temperature-dependent term in β in Equation ( 45) no longer vanishes.This prevents an a priori determination of β diff for the case of given displacement fields.This choice of W p = W p (ε p ) raises concerns about the calibration process.Ideally, such a calibration would be based on stress-strain curves with corresponding temperature-strain curves for the same experiments.If a temperature-strain curve was available, the direct optimization of the material parameters would be possible if the optimization with respect to temperature is replaced by an explicit analytical derivative, as in Equation (42).However, such simultaneously measured temperature-strain curves are rare.If only stress-strain curves are given, the stress and temperature over the entire deformation range must be determined iteratively with the incremental EBVM, which is then followed by a subsequent calculation of the residuals to the given stress-strain curves.This calculation process would have to be performed for all parameter combinations within an expected range, e.g., starting within [0, 1500] MPa for A s .Then, the range should be strategically narrowed, e.g., with a classical bisectional approach.These concerns, which were formulated for illustrative purposes for the special case of W p ̸ = W p (T), were also applicable to the unrestricted general case of the EBVM.In order to be able to use classical optimization approaches for the constitutive parameters, i.e., eventually fitting a curve to the stress-strain curves, one has to make the assumption that, for W p (ε p ), the results of the EBVM are congruent with the results of the heat equation when using β diff = 1 − ∂W p ∂ε p / σ with sufficient accuracy.With this assumption, which is supported by the above considerations, the temperature evolution can be predicted explicitly and is a priori-based on stress-strain curves, thus allowing for an efficient determination of the constitutive parameters.However, in recalling Section 4.5, this special case is identical in its result to the approach proposed by [4,27].This implies the need to justify the practical applicability of the EBVM in its current form with respect to the computational effort involved.This is a particularly crucial consideration, since its attractive advantage of modularity, i.e., the simple extension by arbitrary potentials to incorporate further physical effects (e.g., an extension by a damage term), is not used in this work again due to a limited prescribed processing time.
The second basic approach would be to assume that β diff = β = ∂ψ * ∂ εp / σ and obtain the temperature evolution from solving the heat equation.This assumption, β diff = β, is not theoretically motivated by thermodynamics but rather empirically motivated by the already good agreement with the temperatures obtained by using β approx diff .It can be seen as an initially arbitrary choice of a modeling approach for β diff = β diff (ε p , εp ) without any inherent claim to physical interpretability.With an appropriate calibration, however, it is considered to sufficiently describe the constitutional behavior.Despite a temperatureindependent β diff , the heat equation with a consideration of heat transfer by conduction [4] is still a nonlinear partial differential equation.This is because, for the sake of consistency, the determination of T n+1 , σn+1 = σ(ε p n+1 , εp n+1 , T n+1 ) with its nonlinear dependence on T n+1 must also be used.Nevertheless, the author assumes that the implementation of this nonlinear partial differential equation, or even a corresponding linearization, could be associated with a significant reduction in computational cost for the problem under consideration when compared to the EBVM.Considering this estimated computational advantage, the authors also strongly recommend to evaluate the overall justification of the practical applicability of the energy-based variational model in terms of the trade-off between the increased accuracy of the results by simultaneous solution of the displacement and temperature fields and the computational effort involved.Even for the general case of non-given displacement fields and a temperature-dependent β diff , the results could be approximated by successive evaluation: first, by solving a more classical purely mechanical problem for the displacement field, and then by a purely thermal problem for the temperature evolution.When evaluating the proof of applicability, it should be recapitulated that even the EBVM presented above is a mere approximation due to its incremental character.
Considering the improved suitability for calibration without an additional, though probably valid, assumption (and especially a computational advantage), the modeling of ∂ εp / σ with a corresponding solution of heat Equation (69) will serve as the basic idea of the alternative approach for the determination of the temperature evolution shown in Section 6.

Alternative Approach
As motivated in Section 5, the alternative approach presented here is based on the assumption that β diff = β = ∂ψ * ∂ εp / σ.With this assumption, the partial differential heat equation is solved using the finite element method (FEM).Furthermore, it is simplified that the stress σ always depends on the temperature of the previous time step T n , no matter for which step, either n or n + 1, σ is evaluated.This approximation eliminates the nonlinearity introduced by σ = σ(ε p , εp , T) in the general case of q ̸ = 1 for the JC model used.In future research, this assumption should obviously be omitted, perhaps accompanied by an appropriate linearization to ensure computational efficiency.It should also be noted that both the thermal conductivity λ and the specific heat capacity c p were modeled as a function of temperature T in the actual calculations.Despite this temperature dependence, both were treated as constants in previous sections when calculating the derivatives with respect to temperature T. This simplification should have only a small impact on the resulting temperature fields, but it effectively avoids an unreasonably high level of complexity.To further minimize the impact and reduce the computational cost, both variables λ and c p , as well as σ, are were inserted into the respective equations for the calculation of the temperature increment with their values at the previous time step, i.e., T n+1 = f (λ n , c p,n ).
For the sake of transparency and reproducibility, this section contains a brief derivation of the implemented finite element approach.Since FEM is one of the most widely used methods for solving engineering problems, its basic approach and mathematical foundations are assumed to be well known.Therefore, the derivation is deliberately short and not every single step is provided with a corresponding reference to the literature.In case of uncertainties, the interested reader is referred to the standard literature, such as the books by Belytschko et al. [58] or Zienkiewicz et al. [59].In general, the basic approach of FEM is to divide the whole system into smaller, simpler parts: the so-called finite elements.The necessary discretization in space and, eventually, in time is already given with the displacement fields as input data.Incorporating the shape functions into the governing equations leads to a system of algebraic equations for the elements, which in turn leads to a system of equations for the entire problem.
To allow a discretization of heat Equation (70), it is necessary to integrate this equation over the respective domain, which is achieved by multiplying it by arbitrary weighting functions v ∈ V = {v(x, t) | v ∈ C 1 (x)} to ensure its validity for every single point in the domain Ω as This derivation implicitly assumes that the material is isotropic with respect to the thermal conductivity λ.Using the Gaussian divergence theorem and defining the source term ζ as as well as still using the definition of the nondimensional temperature θ from Equation (34), yields The entire Γ boundary is assumed to be an adiabatic Neumann boundary with ∇T • n = 0. Thus, the corresponding term vanishes in Equation (73).A more detailed treatment of the boundary conditions will follow at the end of this section.A corresponding finite element discretization in index notation yields with shape functions N, thus resulting in T = N • T. T, which are the corresponding discrete values of the temperature field T. The same is true for v.A rearrangement, and the fact that Equation (74) must hold for all arbitrary values of the weighting functions v, leads to Applying a Gaussian quadrature, including a transformation from the physical space x ∈ Ω to the reference space ξ ∈ Ω (re f ) , and summing over all elements results in The indices e and k indicate that the functions and parameters are evaluated and taken, respectively, for the corresponding element or Gaussian integration point.Similarly, w k is the Gaussian weight for the integration point k.The functions c p (T), λ(T), and ζ(ε p , εp , T) are defined as follows The nomenclature εp and εp were used to indicate the discrete nature of strains and strain rates, respectively.The measurement data from experiments eventually yielded a discrete field of displacements, strains, and strain rates at a set of points in the region of interest.This set of points served as the nodes in the FEM approach.The subdivision of this discrete field into quadrilateral elements, as used throughout this work, was performed with a basic algorithm that subsequently numbers and assigns the nodes in both spatial directions.It should be noted that the condition of the overall resulting system of equations can be improved in future research by using a more sophisticated approach to this discretization.
As a further side note, it may be advisable to replace the Jacobian determinant in Equation ( 76) by its absolute value | det J(ξ k )| for the purpose of the actual implementation, since this prevents the calculation from being affected if the transformation into parameter space unintentionally reverses the orientation of the axes.It should also be noted that, since quadrilateral elements and bilinear shape functions are used throughout this work, two Gaussian integration points per spatial direction were sufficient to ensure the accuracy of the integration.
The use of Told in Equations ( 77)-(79) should only emphasize the above introduced approximation that the functions do not depend on the temperature Tn+1 to be determined.Which temperatures are actually used depends on the time integration method used.The one used in this paper to solve an ordinary differential equation (ODE) of the form ∂T ∂t = f (T(t), t) was a classical one-step theta (OST) method, which resolves to [60]  (pp.312-315) where Θ is a numeric parameter to be specified later.Where Equation ( 76) is an ODE of the form M ∂ T ∂t = B T + C, the OST finally leads to a linear system of equations like where the coefficients are calculated as The use of Tn in Equation (82) was again due to the approximate assumption made to avoid the nonlinearity.For the choices Θ = 0 and Θ = 1 of the numerical parameter, the OST method was found to be essentially a forward (explicit) and a backward (implicit) Euler method, respectively.Choosing Θ = 0.5 in the OST method resulted in the Crank-Nicolson method [61].This method is particularly popular for solving diffusion equations, such as the heat equation of [62].Its popularity is due to its stability and improved accuracy compared to the forward or backward Euler methods [63] (pp.[29][30][31]237).Again, in considering the approximate nature of Equation (82), Θ = 0.5 does not exactly yield a Crank-Nicolson approach.Nevertheless, it will remain the value of choice within this work.
This section concludes with some remarks on the physical boundary conditions and their theoretical incorporation.Since the deformation was given, no mechanical boundary conditions were imposed.For the temperature distribution, however, thermal boundary conditions (BCs) were found to be crucial to build a realistic model.For example, heat conduction within the specimen but outside the ROI can be approximated using a Robin boundary condition.Its introduction, which is also accompanied by some further considerations on its incorporation into the FEM approach, can be found in Appendix A.7.This extension of the mesh was again assumed to be adiabatic with respect to its other boundaries.Thus, the neglect of the boundary integral in Equation ( 73) is still valid.

Discussion and Conclusions
The energy-based variational model, which was introduced in Section 4 and evaluated in Section 5, is a promising approach through which to model the constitutive behavior in a thermodynamically consistent way.In particular, its modularity and its extension by arbitrary potentials to include other physical phenomena, such as damage, are attractive features.However, its fully implicit and symmetric formulation, which ensures unconditional stability and excellent convergence properties, as claimed by Stainier and Ortiz [27], are diminished by its adaptation to given displacement fields.These advantages are also offset by the significant computational disadvantages of the implementation in its current form, i.e., as a nonlinear optimization, when compared to the presented alternative approach.Furthermore, the apparent advantage of the EBVM in not using an explicit modeling of the TQC or β diff factor is of limited significance as it is only replaced by the choice of the formulation of the plastic stored energy W p and the dual dissipation pseudopotential ψ * , i.e., the respective splitting of the stress σ in Equation (35).However, the main disadvantage of the EBVM is clearly its apparent inconsistency and limited validity for the temperature-dependent modeling of the stored plastic energy W p (T), as discussed in Section 5.
The alternative approach, which was introduced in Section 6, is a fast, efficient, and comprehensible way through which to determine the temperature evolution for given displacement fields.However, the ALT approach itself is not fully consistent because it approximates the nonlinearity of heat Equation (70) by using temperatures from previous time steps for the source term in Equation (79).It also lacks a full theoretical derivation based on thermodynamic fundamentals and must therefore be considered an empirically motivated approach.
In conclusion, this work shows that the energy-based variational model seems to correctly and consistently describe material behavior only in the special case where the stored energy W p in the plastic is independent of temperature.An exhaustive theoretical investigation of this question is still pending.Furthermore, an efficient and tractable alternative approach for the determination of the temperature evolution for given displacement fields is proposed.

Outlook
Regardless of whether the EBVM or ALT approach is pursued, several issues need to be addressed to improve the accuracy of the calculated temperature evolution.First, the influence of macroscopic friction should be incorporated into the respective models for a more general analysis of different experiments or processes.One approach would be to include an accompanying mechanical simulation that provides, for example, any contact forces, shear, and slip velocity in the experiment from which the displacement field data are taken.From these values, a thermal boundary condition could be derived for a purely thermal ALT approach.
Further work could also consider elasticity.This area of research is less about the direct influence of the Joule-Thomson effect on the temperature evolution but more about the high values of hydrostatic pressure encountered.A separation of this hydrostatic pressure could allow for a distinction between elastic and plastic strains.This could lead to a more accurate evaluation of the stress state associated with the dissipative part of the irreversible deformation work, thus resulting in a temperature evolution.
Furthermore, the range of materials studied could be extended, and the EBVM and ALT approaches could also be applied to other constitutive models in a relatively straightforward manner.
Regarding the energy-based variational model, the obvious next step is the theoretical evaluation of the root cause of the problems arising for temperature-dependent plasticstored energy terms W p (T).If it is not possible to eliminate the root cause, it might be recommendable to evaluate alternative incremental variational formulations, starting from the continuous formulation proposed by Yang et al. [24], e.g., the one published in Canadija and Mosler [64] (albeit they used internal energy instead of free energy).According to Stainier [46], this type of approach leads to a significantly higher complexity in the construction of analytical expressions.If the root cause can be eliminated, a next step might be a full derivation with respect to the temperature, as shown in Equation (42), and this would be accompanied by a linearization and a corresponding implementation with an FEM or Rayleigh-Ritz approach to reduce the computational effort.As a side note, a nonlinear FEM approach for the EBVM was developed and implemented during the preparation of this work.Its results were close to those obtained with the ALT approach for a stored plastic energy term that was independent of temperature.However, the computational efficiency was significantly inferior to the ALT approach due to the incorporated nonlinearity.
With a valid EBVM, another possible next step would be to extend it to other constitutive phenomena, e.g., to include a damage potential.This damage potential could be in the style of Lemaitre, who introduced a description of damage mechanics based on a scalar damage variable D ranging between D = 0 for the undamaged state and D = 1 for the rupture of the respective volume elements.A corresponding translation into potentials has been published by Lemaitre et al [45] (pp.350, 396-409).In the present version of the EBVM and the ALT approach, the softening effect of material damage was implicitly and undifferentiatedly included by calibration with actual experimental results.
Furthermore, the EBVM in its present form for the use of displacement field data could be extended toward dynamics.A basis for such an adaptation is provided by Stainier [34].
Since the theoretical root cause analysis for the shortcomings of the EBVM is still missing, an efficient and comprehensible alternative approach is presented in Section 6.The most important next step is the establishment of a solid theoretical foundation since the ALT approach in its current form is, rather, an empirically motivated approach.However, this may be a by-product of the in-depth theoretical evaluation of the EBVM.Furthermore, for the sake of consistency, the source term in Equation (79) should be evaluated with a dependence on the correct temperature, i.e., it should not be approximated by using the temperature of a previous time step.In its current form, this would introduce a nonlinearity for the general case of the JC parameter q ̸ = 1.In order to avoid an unreasonable increase in the computational effort, a corresponding linearization is recommendable.Although a holistic calibration-i.e., including the thermal material behavior, which is based exclusively on stress-strain curves and may be a central advantage of the presented approach-its mathematical implications, especially with respect to the suspected local minima, should be subjected to a more profound theoretical analysis.Here, the experimental work on the time-resolved identification of the temperature-dependent plastic-stored energy terms W p (T) is especially promising for the purpose of further validating the approaches.
With regard to practical applications, the presented approaches take up the increasing capabilities of optical measurement systems for displacement and strain field identification.Using the presented methods enables one to directly derive parameters for coupled thermomechanical problems, e.g., in manufacturing processes like metal forming or machining.
In the first case, it was assumed that there was no separate Helmholtz free energy, i.e., W = W p .When utilizing the usual equation ρc p = T ∂N ∂T [65]  The term multiplied by εp is exactly the factor β SO diff .In the second case, there was thermal Helmholtz free energy (as was the case in the EBVM), i.e., W = W p + W th .The adiabatic heat equation can then be written as Without further expansion of the term, it was stated that T ∂ 2 W p ∂T 2 Ṫ = 0 ∀ q = 1.Thus, this additional term could be the source of the observed discrepancies for the general case of q ̸ = 1. of brevity, it should be emphasized that the above considerations are by no means an exhaustive mathematical proof; rather, they should be seen as considerations from an engineering point of view.