An experimentally-fitted thermodynamical constitutive model for polycrystalline shape memory alloys

A phenomenological model for polycrystalline NiTi shape-memory alloys with a refined dissipation function is here enhanced by a thermomechanical coupling and rigorously analyzed as far as existence of weak solutions and numerical stability and convergence of the numerical approximation performed by a staggered time discretization. Moreover, the model is verified on one-dimensional computational simulations compared with real laboratory experiments on a NiTi wire.

1. Introduction. Shape-memory alloys (SMA) are intermetallic materials with an extra-ordinary response to mechanical and thermal loadings. Namely, they are able to recover large strains (commonly 5-7 %) in the high-temperature regime and seem to react plastically at low temperatures. But, interestingly, upon heating they are able to recover their original shape; which is referred to as the shape-memory effect. At the heart of these remarkable properties is a diffusionless phase transition accompanied with a reversible rearrangement of the crystal lattice. During this phase transition, that can be induced by thermal and/or mechanical loads, the material transits between a higher symmetric (typically cubic) high-temperature phase austenite and a low temperature phase martensite. The martensite phase arises in a number of variants which allows for formation of microstructure; then, by rearrangement, or reorientation, of this microstructure the material may compensate the stress due to applied loads, which ultimately leads to the remarkable response described above.
The unique features of SMA as well as their multiscale character make them appealing from the point of view of mathematical modeling so that models at different scales and with different modeling focuses have been proposed, see e.g. [57,17,18] for overviews. Nevertheless, from the point of view of applications, in particular in medical devices [30], automotive industry [37] or in space exploration [19], models for polycrystalline materials that allow for numerically stable and convergent discretization schemes, efficient computer implementation, and for usage of fast algorithms, as well as for easy adaptation of material parameters are particularly desired. Hence, a large number of such models has been proposed to date, e.g. [2,16,25,26,72,66,61,45,15]. In [62,27] a phenomenological model for NiTi (commercially called Nitinol) SMA was proposed that features an elaborated dissipation mechanism which seems to be crucial for an accurate prediction of the response of a SMA specimen to non-proportional loads. The dissipation mechanism can be derived either by comparison with experiment [62] or by microscopically motivated arguments [27]. In [27] also the existence of energetic solutions (see also e.g. [50]) has been shown. Yet, the model analyzed in [27] contained the temperature only as a "control" field; full thermomechanical coupling of the evolution of the mechanic and internal variables to the heat equation was not included (a similar approach has been followed also in e.g. [49,48]). However, the combined effect of latent heat released during the phase transformation and the strong coupling between stress and temperature on its initiation and propagation substantially influence the ensuing behaviors of the SMA material. This is particularly manifested in experiments involving higher loading rates and/or impeded heat transfer from (or to) the transforming material [69,54]. Moreover, even at quasi-static loading rates at stable temperature, common NiTi specimens (e.g. wires, ribbons, tubes) tend to transform in a highly inhomogeneous, localized manner [64,8,65,63]. Combination of all aforementioned effects gives rise to a complex pattern of transformation bands visualizable, e.g., by digital image correlation technique or thermal camera [73,54,56]. Naturally, when modeling such situations, the temperature cannot be prescribed but rather has to be calculated in each material point using the heat equation that takes heat production by latent heat as well as dissipation into account, cf. [35,59,5,40,3,42,71]. This is also the main goal of this contribution: that is, we couple the model formulated in [62,27,63] to the heat equation, provide a sound mathematical analysis and assess the performance of the model in simple numerical tests.
To this goal, we will augment the rate-independent dissipation potential from [62] by a (presumably) small viscosity in the internal variable that directly couples with temperature in the free energy (i.e. the volume fraction of martensite ξ in our case). This makes the model and its analysis simpler, although the thermal coupling brings own specific demands.
The devised model is consistent with the basic requirements from thermodynamics and maintains the important properties of the dissipation mechanism proposed in [62,27]. We moreover prove that it is mathematically well-posed in the sense that it possesses weak solutions. We prove the existence of weak solutions by time discretization (Rothe method) which, at the same, provides a guide for numerical implementation. This is also performed in 1D setting for demonstration purposes. Simple numerical simulations are conducted and supplemented by an experimental Ω a domain in R d (a specimen) K bulk modulus Γ the boundary of Ω G = G(ξ) (a part of) shear modulus u displacement G = G(ξ, e) shear strain energy ε(u) total (small) strain = 1 2 (∇u) ⊤ + 1 2 ∇u cv heat capacity, e, π elastic and inelastic strains K heat-transfer tensor ξ volume fraction of martensite θ D prescribed outer temperature θ temperature u D prescribed boundary displacement s entropy sam entropy of A/M-transformation w heat part of the internal energy θtr the transformation temperature η a regularization parameter for approximation τ a time step for time discretisation Table 1. Main nomenclature for variables and for the data of the model.
visualization. This article is structured as follows: First, in Section 2, we formulate the mathematical model and present its thermodynamics. In Section 3 we then perform the analysis as far as the existence is concerned. Here we devise an efficient staggered time discretization with a regularization and show its numerical stability (a-priori estimates) and convergence. Eventually, in Section 4 we illustrate the model on one-dimensional computational simulations of a NiTi polycrystalline wire compared with real laboratory experiments.
2. The thermomechanical model. In this section, we will introduce the studied thermomechanical model. Even though all quantities will be explained at the particular spot where they appear first, we summarize them also in Table 1 for a better orientation.

Basic ingredients.
We formulate the present model in the framework of socalled generalized standard models (GSM, also called generalized standard materials or generalized standard solids) developed in [34]. Within this framework, we have to prescribe the free energy ψ, the dissipation function r (or sometimes called dissipation potential) and a Fourier law for the heat flux to build-up the constitutive model. Let us note that the GSM framework is one of many approaches to non-equilibrium thermodynamics (see also [53], [24]) that is very convenient in continuum mechanics of solids [41].
The free energy depends on the state variables while the dissipation function may depend on the state variables as well as on their rates. Here, we have to respect the constraint that the dissipation function has to be non-negative, convex in the rates of the state variables and vanishing if the rates are zero; given this constraint we can assure that the devised model satisfies the second law of thermodynamics in the case of the absence of heat conduction, cf. [51]. If heat conduction is present, also the Fourier law needs to be chosen in a thermodynamically consistent way.
We denote by Ω ⊂ R d the domain occupied by the specimen and, as the primary state variables, we choose the displacement u : Ω → R d , together with the corresponding small strain tensor ε(u) = 1 2 (∇u) ⊤ + 1 2 ∇u, and the temperature θ : Ω → R. Moreover, to capture the multiscale behavior of SMA, we introduce two internal variables: the volume fraction of martensite ξ : Ω → [0, 1] and the inelastic strain tensor π : Ω → R d×d . The approach of introducing internal variables is common within the GSM framework [47].
The inelastic strain enters the (small) strain decomposition as follows where e represents the common elastic part of the strain tensor. Thus, in our model, the inelastic strain stores all the information about the microscopic internal structure of martensite. This can be seen even clearer, if we rewrite here e tr represents a "homogenized" mean transformation strain of martensite and was used as a primary variable [72,16,59,42]. However, working directly with the inelastic strain will allow us to work with a convex free energy in all the variables except the ξ-variable. Notice that it follows from (1) that other contributions (than those coming from the reversible transformation) to the inelastic strain (like plasticity) are not included in the model at this point. As far as the inelastic strain is concerned, we assume it to be traceless and restricted as follows: where tr(x) denotes the trace of the tensor x and · : R d×d → R + is a positively 1-homogeneous convex function that is to be specified for the individual material. Indeed, the relative volume change associated with martensitic transformation is very small [36] which justifies the vanishing trace of the inelastic strain. In addition, crystallographic considerations show [52] that there exists a maximum value of strain that is attainable due to phase transformation; so, the transformations strain is considered to lie in a bounded convex set; in our modeling this set is given by e tr ≤ 1. In view of the relation between the inelastic and transformation strain, this gives the constraint in (3). We now introduce the free energy and dissipation functions. Let us start with the free energy, which is assumed to be of the form ψ(ε(u), π, ξ, θ) = 1 2 K(tr e) 2 + G(ξ, e) + ζ(π, ξ) + ϕ(ξ, θ) where K is the bulk modulus so that 1 2 K(tr e) 2 the bulk elastic strain energy while G(ξ, e) is the shear elastic strain energy being influenced by the volume fraction ξ. Later, we will need it to have at most linear growth in e so that, instead of a "first-choice" quadratic ansatz G(ξ, e) = G(ξ)|dev e| 2 with G = G(ξ) being the shear modulus, we will have in mind rather with a small regularizing parameter α > 0 and with G M > 0, cf. also [43,Example 7.5.6]. The subsequent term ζ(π, ξ) is a non-convex term that we assume convex in π and of double-well type in the volume fraction ξ. Including such terms is motivated by the so-called localization of martensitic transformation, a common phenomenon in NiTi SMA, in which the austenite-martensite transformation proceeds not in a homogenized but rather in a localized manner; cf. the recent works [4,1,38,55]. The term ϕ(ξ, θ) is the part of the chemical energy driving the thermomechanical coupling, and it is chosen in such a way that ϕ(ξ, 0) = 0; the other terms (e.g. internal stored energies of various kinds, cf. [16,72,42]) are formally included in ζ(π, ξ). Further, δ S (·) is the indicator function of a convex set S, i.e. δ S (π, ξ) = 0 if (π, ξ) ∈ S while δ S (π, ξ) = ∞ otherwise. The indicator function will assure the constraint (3) as well as the fact that the volume fraction of martensite takes only values between 0 and 1. Thus, the mentioned convex set is given by Finally, the terms ν 2 |∇π| 2 and ν 2 |∇ξ| 2 are needed for mathematical reasons but, since they incorporate information on internal length scales of the material (via ν, cf. [33]), they also open the space for convenient handling with localization processes [7,39].
As far as the dissipation function is concerned, we shall use a form that consists of a rate-independent contribution r RI and a viscous one, denoted r VI . For the dominant, rate-independent part of the dissipation function, we choose the form used in [62] that was shown to produce predictions that correspond well to experiments. Indeed, as demonstrated in [31], in quasi-static loading under strictly isothermal conditions, NiTi SMA behave in a rate-independent manner. However, to conform also to high-driving-rates situations, the model shall be extended by a rate-dependent contribution; this is justified by the fact that phase boundaries require an unlimitedly increasing driving force to propagate at speeds reaching towards some sound speed [10]. For the rate-dependent part, we choose the simplest possible option and take (uncoupled) terms quadratic in the rates. Finally, we get: where a AM (·) > 0, a MA (·) > 0, and a reo (·) > 0 are given coefficients depending continuously on their arguments; such a form has already been used in, e.g., [9]. Note also that r TOT (π, ξ, θ; ·, ·, ·) is convex and non-negative, and r RI (π, ξ, θ; ·, ·) is positively homogeneous (of degree 1); this is why r RI is called the rate-independent (part of the) dissipation potential. The exponent "2" in the rate-dependent contribution (7c) was considered as a "material parameter" p > 2 in [60], see also Fig. 2 therein.
Here, however, from analytical purposes below (namely the L 2 -estimate of ∆π and ∆ξ needed for the energy conservation), we will use p = 2. Let us also remark that, for a small viscosity parameter µ > 0, this model is not substantially different from the purely rate-independent one (at least for lower strain rates) which would, however, exhibit analytical difficulties with measurability-in-time of the displacement field. The free energy specified in (4) as well as the dissipation function in (7) are given in each material point; therefore, they can be understood as densities. However, in order to derive the balance equations, we need to work with their global forms. This reads for the total free energy as where we used that tr ε(u) = div u. Moreover, we introduced in (8) an outer mechanical loading by a prescribed displacement u D (t). The modeling idea here is that the body is attached to a hard device, described by u D (t), by (infinitesimal) elastic "spring" with the elastic moduli κ 1 .
The total dissipation function is given by The functional derivative (i.e. Gâteaux differential) of the convex functional (u, π, ξ) → G (u, π, ξ, θ) then drives the mechanical part of the system, which then can be written abstractly in the form of a doubly nonlinear generalized gradient flow: with R = R(π, ξ, θ; . π, . ξ) given by (9) and actually independent of . u so that ∂ ( . u) R = 0. Another view on (10) is to understand it as a force balance.

Remark 1 (Energy mixing).
By introducing an individual internal variable for the net (averaged) transformation strain (see (2)), we cannot trace the individual martensitic variants; i.e. we work with a single tensorial variable capturing the strain state of martensite and one scalar variable capturing its volume fraction as already used in [59,46,67,32,15,42]. Thus, we consider the martensitic structure in an averaged sense and abandon the concept of a mixture of austenite and several martensitic variants; the latter can be historically traced back, although usually not explicitly mentioned in literature with only rather rare exceptions [6,72], to the phenomenological model devised by M. Frémond [25,26], largely scrutinized in mathematical community, cf. e.g. [12,20,21,22]; see also [57,Rem.4.4].

Remark 2.
Usually, the rule of mixtures is used for derivation of the free energy of the austenite-martensite mixture [44]. E.g., for the chemical part of the free energy, one gets so that the heat capacities in pure martensite and in austenite are −θϕ ch M ′′ (θ) and −θϕ ch A ′′ (θ), respectively. A simple rearrangement allows the general form Hence, ϕ 1 (θ) controls thermal induction of martensite (under stress-free condition) and the entropy reads

Partial differential equations and inclusions.
Plugging in the particular forms (8), (9) and (7) into (10), we obtain the balance of (linear) momentum as well as the so-called flow rule governing the evolution of the internal variables. In order to close the thermomechanical system, we still need to add a balance equation for the entropy s defined by the Gibbs relation In the absence of radiative heat sources, the change of heat in a control volume is only due to heat flux and the dissipation. Thus, we have that with the heat flux q given phenomenologically by the Fourier law q = −K(ξ, θ)∇θ with K(ξ, θ) the heat-conductivity tensor. The right-hand side of (14) is the dissipation rate.
To simplify the form of the heat-transfer equation but without restricting substantially possible applications, we assume ϕ(ξ, 0) = 0. Substituting the Fourier law and the Gibbs relation (13) into (14), we obtain the heat-transfer equation with In the classical formulation, considering also the heat-transfer equation (15), we thus are to solve the following system of equations/inclusions where N S denotes the normal cone to the set S from (6) and r is from (14). We accompany the system (17) by boundary conditions ∇π·⃗ n = 0, ∇ξ·⃗ n = 0, (18b) where ⃗ n denotes the outward unit normal to the boundary Γ , θ D is the outer temperature, and κ 2 ≥ 0 is a phenomenological heat transfer coefficient at the boundary. Moreover, we prescribe the initial conditions Note that we did not prescribe any condition on the deformation at the initial time. Instead, we assume that even at t = 0 the system is in equilibrium so that the initial displacement is calculated from (17a).

Remark 3 (The meaning of w).
At times, we shall refer to w or C ϕ (ξ, θ) as the thermal part of the internal energy. This is justified by the following observation: Realizing that u = ψ + θs is the density of the internal energy, we can calculate the difference from w = C ϕ (ξ, θ) as so that C ϕ differs from ψ + θs only by function dependent exclusively on the mechanical variables.
Remark 4 (Coupling of π and ξ). Let us note that the evolution of the internal variables π and ξ is mutually coupled both through the dissipation term ∂ r RI and through the constraint N S (·). Therefore, we cannot deduce separate flow-rules for π and ξ but instead have to keep it coupled as in (17b).
Remark 5 (Coupling of θ and e). Let us note that, above, the balance of linear momentum (17a) does not directly depend on temperature. This is because we assume isotropic temperature-independent elasticity of the material. Let us note that according to a recent investigation on polycrystalline specimens [68], elastic properties of both phases exhibit anisotropy (pronounced in case of stress-oriented martensite) and temperature dependence. However, a consistent and robust way to capture the complex evolution of total (anisotropic) elasticity tensor with both internal variables and temperature has not been established yet; first attempts in this direction have been made in [29]. Moreover, as in [13,60,72,62,42], we neglect thermal expansion; thermal strains for temperature variation within the range of 100 K do not exceed 0.1% [70], which is much less than the typical maximum transformation strains relevant for NiTi SMA applications (5-7 %). However, it can be easily adopted in the model [2,45,29] while the analysis becomes more involved due to the adiabatic effects in the heat-transfer equation, cf. e.g. [ . ξ), respectively, and integrating over Ω. Thus, we obtain from (17a) together with boundary condition (18a) that Furthermore, multiplying (17b) with ( . ξ), integrating over Ω with using the Green formula, and exploiting also the calculus ∂ ∂t ϕ(ξ, θ) = ∂ ξ ϕ(ξ, θ) .
. ξ, and the boundary condition (18a) gives that Summing (20) and (21) up, and integrating over a time interval [0, t] gives the mechanical energy balance: where the "stored" energy E is the mechanical part of the free energy G from (8), namely Finally, integrating (17c) over Ω and using boundary condition (18c) leads to and adding this to (22) and using the Gibbs relation leads to the internal energy balance: work of external loading and heating (24) where U is the total internal energy As already mentioned, in order to satisfy the second law of thermodynamics, we have to choose the dissipation function non-negative, zero at zero rates and convex in the rate variables; this assures that the total dissipated energy is non-negative. Moreover, we assume that the heat-conductivity tensor K is positive semi-definite, which is a natural requirement. Then, we see from (14)  . π, 3. Analysis of the system (17)- (19). In this section, we show that the model devised in Section 2 is mathematically well-posed in the sense that it possesses weak solutions and simultaneously we devise a stable and convergent approximation.When scrutinizing a suitable notion of weak solutions, we have to deal with doubly-nonlinear structure of the balance equations for the mechanic variables; cf. (10). In isothermal situations, analysis of such structures has been pioneered in [23] but here the situation is more general because we allow for a state-dependence of the dissipation function and thermal coupling. We now consider a bounded domain Ω ⊂ R d with a Lipschitz boundary Γ = ∂Ω. Beside the standard notation for the Lebesgue L p -spaces and W k,p for Sobolev spaces whose k-th distributional derivatives are in L p -spaces, we will use the abbreviation H k = W k,2 . Moreover, we use the standard notation p ′ = p/(p−1), and p * for the Sobolev exponent We consider a fixed time interval I = [0, T ] and we denote by L p (I; X) the standard Bochner space of Bochner-measurable mappings I → X with X a Banach space. Also, W k,p (I; X) denotes the Banach space of mappings from L p (I; X) whose k-th distributional derivative in time is also in L p (I; X). Using the shorthand notation Q = I × Ω, we will often use the isomorphism L p (I; L p (Ω)) ∼ = L p (Q).
In order to cast a suitable definition of weak solutions to (17)- (19), we rewrite (10) in the form of a system for some thermodynamical driving force F = (0, Π, Ξ) and the external timedependent loading f = f (t) arising from the boundary conditions, namely ⟨f (t), v⟩ := ∫ Γ u D ·v dS, and then use the definition of the convex subdifferentials. Thus, together with the standard weak formulation of the heat-transfer problem, we arrive at the following definition: (17)- (19)). The seven-tuple (u, π, ξ, θ, w, Π, Ξ) will be called a weak solution to the initial-boundary-value problem (17)
We will assume the following data qualification which are fitted to application both towards the above specified attributes and towards a relatively simple analysis of the problems, without claiming that they are the weakest possible. More specifically, we assume C ϕ , ∂ ξ ϕ, r RI , G, and K continuous, satisfying K : [0, 1] × R + → R d×d bounded and uniformly positive definite, (28b) r RI (π, ξ, θ; ·, ·) convex, 0 ≤ r RI (π, ξ, θ; . π, with R d×d dev = {e ∈ R d×d ; tr e = 0}. Let us note that the coercivity of C ϕ in (28a) allows for the heat capacity ∂ θ C ϕ bounded (not growing with increasing temperature), which is a physically relevant assumption. Furthermore, we will qualify the initial conditions and the boundary loading as We will prove existence of weak solutions in the sense of Definition 1 to the system (17)- (19) in two steps. First, using a small parameter η > 0, we regularize the right-hand side of the heat equation to make it uniformly bounded. This will allow us to deduce a-priori estimates even if we decouple the mechanical inclusions from the heat equation.
Once this regularization is performed, we will discretize the system in time, which will lead (after another discretization in space) to a conceptually implementable, numerically stable, and convergent scheme. To this goal, we consider a time step τ > 0 which does not vary within particular time levels and such that T /τ is integer, leading to an equidistant partition of the considered time interval. Let us emphasize, however, that a varying time-step and non-equidistant partitions can be easily implemented because we will always consider only first-order time differences and one-step formulas. The mentioned regularized system is: r TOT π η , ξ η , θ η ; . .

Lemma 1 (Existence of approximate solutions). For every
Proof. Since the approximate problem is decoupled, we solve first the problem for (u k ητ , π k ητ , ξ k ητ ) and then for (θ k ητ , w k ητ ). Unfortunately, the boundary-value problem (32a,b)-(33a,b) does not seem to admit for a potential and existence of its weak solution is ensured only nonconstructively through the Browder-type fixed-point theorem exploiting at least the potentiality of the nonsmooth part arising from the convex dissipation potential and from the convex constraints due to the set S, cf. [58,Sec.5.3].
The boundary-value problem (32c,d)-(33c) is even simpler, having a potential and admitting thus direct-method arguments.
Having established the existence of regularized and discretized weak solutions, our goal will be to establish suitable a-priori estimates in order to be able to pass to the limit τ → 0. Here we will exploit that, in the heat equation, the dissipation as well as the term related to latent heat are uniformly bounded (albeit depending on η that we will keep fixed for the first limit passage). Later, we will pass with the regularization η → 0.
and every such (u η , π η , ξ η , θ η , w η ) is a weak solution to the regularized system (30) with the corresponding boundary and initial conditions.
Proof. The selection of a subsequence weakly* convergent in the topologies indicated in (37) is by the Banach selection principle. The strong convergence (40a) can be seen from the estimate where ϵ is from (28f). This estimate can be obtained by testing (35a) by u ητ − u η ; on the right-hand side we employ the weak convergence of u ητ towards u η as well as π ητ towards π η and the strong convergence of ξ ητ due to the Aubin-Lions theorem in combination with (28e). The strong convergence (40b,c) is quite delicate. It is important to have the estimates (37f) which will be inherited for the limit, i.e. we know that ∆π η ∈ L 2 (Q; R d×d ) and ∆ξ η ∈ L 2 (Q). We can then use (38) and, by the convexity of r(π, ξ, θ; ·, ·) and of the stored energy, we can estimate Q r π η , ξ η , θ η ; .
The strong convergence of w ητ is due to the compact embedding, based on the estimates (37d) by the Aubin-Lions theorem and an interpolation by the Gagliardo-Nirenberg inequality. In fact, the latter estimate in (37d) is to be modified for w ητ by replacing L 1 (I) by measures on I, and the Aubin-Lions theorem is to be generalized appropriately, cf. [58,Cor. 7.9].
The convergence (in terms of the above selected subsequences) of the weak solutions of (35)-(36) towards the weak (30) with the corresponding boundary conditions is then easy. (30)). Let (28)- (29) hold. Then we have the following estimates
By the boundedness and coercivity of the stored energy, this test yields the L ∞estimates in (43a-c).
Then we perform this test once again but with (30c,d) tested by some small positive number less than 1 in order to see also the dissipation on the left-hand side.
The estimate of ∇θ η in L r (Q; R d ) in (43d) needs a nonlinear test of the heat equation and an interpolation by Gagliardo-Nirenberg inequality, see [11] or also e.g. [43,Prop. 8.2.1] or [58]. Here we need the growth (28c) of ∂ ξ ϕ which, when involved in the mentioned interpolation game, ensures ∂ ξ ϕ(ξ η , θ η ) bounded in L 2 (Q) so that the adiabatic-heat (also called latent-heat) term ∂ ξ ϕ(ξ η , θ η ) . ξ η is bounded in L 1 (Q). We refer to [43,Prop. 8.3.1] for details about the interpolation and the estimation. Furthermore, the estimate of ∇w η in L r (Q; R d ) in (37e) follows from the calculus like (39), i.e. (44) and the already obtained estimates on ∇ξ η and ∇θ η . Then the latter estimate in (43d) follows by comparison when using (30c). Eventually, for the estimates (43f) we still test the particular inclusions in (35b) by ∆π η and ∆ξ η and use analogous arguments as in the proof of Lemma 2.
The proof of the following assertion then copies the arguments from the proof of Proposition 1.

Numerical implementation and 1D illustrative simulations.
In this section, we illustrate the coupling of phase transformation with heat transfer on a one-dimensional reduction of the model. Our computation domain will be a line segment representing a piece of an SMA wire of the length l and of the radius ρ subjected to tensile loading. Thus, the dimension d = 1 and Ω will be the interval (0, l).

4.1.
Implementation of the discrete model. Here, we specify the material functions and parameters introduced in a general form in Section 2. We complement the staggered implicit time discretization suggested in (32) by a suitable spatial discretization and approximate the functional spaces utilizing the finite-element method. Particular form of the constitutive model and material functions: As for G(ξ, e), we use the specific ansatz (5) with the prefactor G(ξ) as follows so that the shear energy term proposed in [62] (corresponding to the Reuss model) is recovered for α = 0; let us note that (28f) is then satisfied (by choosing, e.g., ϵ = G M ). Starting from assumption (11) and following the procedure described in [62], we arrive at a common form of the chemical energy as, e.g., in [62,66,40,45]: where θ tr and s am are positive physical constants, namely the phase-equilibrium temperature and the entropy variation related to the phase transition (at θ tr ), respectively. ϑ is a constant, see [62] for details on physical background of applied approximations. Now, (47) is compatible with ϕ(ξ, 0) = 0, by setting ϕ(ξ, θ) = ϕ ch (ξ, θ) + s am θ tr ξ, adding the term −s am θ tr ξ to ζ(π, ξ) and setting ϑ − c v θ tr = 0. Let us note that this is possible due to the specific form of (47); then it also holds ∂ ξ ϕ(ξ, 0) = 0 as assumed above, cf. (28c). Moreover, for simulations of simple stretching of a wire in our one-dimensional situation, the constraint on π in (6) determined by function · from [62] degenerates into the relation 0 ≤ π ≤ f tens ξ with the material parameter f tens > 0 (maximum transformation strain in tension).  Functions a reo (θ), a AM (ξ), a MA (ξ) in (7) are linear in their arguments, and take the form [62] a reo (θ) = σ reo tr are material parameters and Σ reo is a material constant. We define ζ(π, ξ) so that the non-convex part (accounting for localization effects) is adopted from [63]: where E int , E nl , β are non-negative material parameters.
In the first set of numerical simulations, we put E nl = 0 and omit the gradient terms from (4). In the second set of simulations, we choose E nl > 0. For the implementation of the π-and ξ-gradient terms, we made (with a slight shortcut) a nonlocal regularization, which anyhow did not substantially influence the results when taking the coefficient ν > 0 small.
Computational domain, initial and boundary conditions: We use a equidistant partition of a line segment of length l to n parts, each of length l n , via n+1 uniformly distributed nodes. We consider Dirichlet boundary conditions for the mechanical loading, fix displacement at one end of the 1D specimen (a straight wire), and prescribe time-dependent displacement at the other one. For the heat-transfer equation, we consider Neumann boundary conditions with prescribed heat fluxes at both ends. We prescribe a spatially constant initial temperature θ 0 > A f so that the wire is completely in austenite at t = 0 and, consistently, ξ 0 = 0, π 0 = 0 everywhere.
Discretization: Piecewise linear (P1) elements are used for spatial discretization of displacement u and temperature θ, whereas the internal variables ξ, π are discretized via piecewise constant (P0) elements. With respect to the fractional-step strategy, we decouple the system at each time step k = 1, . . . , T /τ so that we solve first for (u k ητ , ξ k ητ , π k ητ ) from the mechanical part (temperature fixed) via searching for a critical point (cf. [27]) and then temperature θ k ητ is recovered from the heat equation (fixed displacement and internal variables); then we move to the next time-level. The equidistant time partition is used.

Illustrative 1D simulations and comparison to experiments.
The discrete staggered system derived in the spirit of (32) was implemented in a MATLAB computational environment. The material parameters are summarized in Table 2; they were set to fit the experimental data (see below) and correspond to typical values for NiTi alloys used in the literature [35,14,62].
In the first set of numerical simulations, we compare the isothermal and the adiabatic response of the wire to a simple axial stretching first up to 7%, then down to 0%. Since E nl = 0, any physical quantity is distributed homogeneously along the wire. 1 Figure 1 shows the dependency of the stress on the total strain for the isothermal response (with temperature prescribed and constant) at θ iso = 20 • C and for the adiabatic response (completely thermally insulated domain, κ 2 = 0) with initial temperature θ 0 = 20 • C. Three different rates of change of total strain are presented in the adiabatic case, namely v 1 = 10 −3 s −1 , v 2 = 5·10 −2 s −1 , and v 3 = 5·10 −1 s −1 . For the highest rate, results are shown for three different time increments -namely 10 −2 s, 5·10 −3 s, and 10 −3 s -to demonstrate that with a time-step refinement the results are closer to each other. An analogous tendency was observed also for other strain rates, hence, the results from one time-discretization are presented.
In all situations, the phase transformation from austenite to martensite is induced during stretching and the reverse one is initiated during unloading. Comparing isothermal and adiabatic simulations, substantial growth of "strain-hardening" (increase of stress with increasing strain in the "plateau") is observed for the latter ones with respect to the former. This is predominantly due to the latent heat released/absorbed within the transforming material and the coupled influence of stress and temperature on coexistence of a phase mixture (analogy to the Clausius-Clapeyron relation): as phase transformation is underway, released/absorbed heat heats/cools the wire which in turn leads to variation of phase-equilibrium stress. This is documented on Figure 2, which shows the evolution of temperature and volume fraction of martensite in an arbitrary element of the domain with respect to the relative progress of stretching (horizontal axis: 0 denotes start of stretching, 0.5 maximum tension, 1 complete unstraining). Finally, Figure 3 shows the cumulative sum of the respective contributions to the total heat, i.e. latent heat, energy dissipated in rate-independent processes and in rate-dependent (viscous) processes, in the course of stretching (cf. integral contributions in (22)). Since the latent heat is released during forward transformation and absorbed during the reverse one, the net contribution is zero at the end of simulation for rates v 1 , v 2 . However, the temperature at the end of the simulation is slightly higher than the initial one due to the contribution of dissipated energy. Thanks to a relatively small contribution of viscous dissipation in comparison with the rate-independent dissipation, results for v 1 and v 2 are close to each other. In case of the highest rate, v 3 , the contribution of viscous dissipation is comparable to the rate-independent one, which results in several effects: i) the thermodynamic force needed for initiation of transformation is higher, hence transformation occurs at higher stress, ii) at 7% strain, slightly less martensite is induced in the material, iii) the maximum attained as well as final temperatures are higher, iv) reverse transformation is not finished at 0% strain (cf. Figs. 1 and 2).
In the second set of simulations, we include the gradient regularization of internal variables and study the localization phenomenon in an SMA wire. The computational domain now represents a half of the wire, the symmetric boundary conditions are imposed by fixing displacement and prescribing zero heat flux at one end of the domain. At the other end, strain is uniformly increased. The Neumann (convective heat transfer) boundary condition with a high heat transfer coefficient (κ 2 = 4 kW/(m 2 • C), cf. [14]) mimics fixation of this end of the wire to conductive grips of the loading device. The ambient temperature, θ amb = 20 • C, is fixed within lat. heat, v 3 Figure 3. Evolution of cumulative sum of three respective contributions to the total heat -latent heat, energy dissipated in rateindependent (RI) processes and in rate-dependent (RD) processes, respectively -integrated over the whole wire at three different strain rates. Relative progress of stretching on the horizontal axis (0 denotes start of stretching, 0.5 maximum tension, 1 complete unstretching). The RI contribution for v 2 basically coincides with the RI contribution of v 1 .  the simulations. To control the spot where the phase transformation initiates, a slightly decreased value of the parameter E nl is prescribed in the boundary element at the fixed end of the wire. Let us note that various types of such an artificial initial nucleation, which mimics fixation to the grip, are common in localization models [63,35]. Distance along the wire from one end is on the vertical axis, relative time from the start of the experiment/simulation on the horizontal one. The color at any point within the plot represents the local strain. Ambient temperature of air as well as initial temperature of the wire in the experiment was set to 20 • C, the same were set in simulations. Let us remind that symmetry of solutions was assumed in simulations, hence mirror-symmetric plots are presented.
For v s (Fig. 4), two localized transformation bands are initiated at the grips and move towards the center of the wire, both in the experiment and in the simulation. In contrast, during faster stretching at v f (Fig. 5), additional transformation bands occur along the wire in a clearly separated manner. This is because of the local increase of temperature in and close to nucleation sites, which impedes spreading of transformation bands along the wire. Specifically in the simulation, the transformation band nucleated at the end of the wire (approx. in time t 1 = 0.2 s) starts to propagate towards the middle of the wire, increasing the temperature in its close neighborhood. The rise in temperature impedes further transformation (recall the Clausius-Clapeyron relation), hence, the band halts and a new transformation band appears in the the middle of the wire (approx. in time t 2 = 0.5 s), where the temperature growth is slowest. This again leads to a local heating of the material nearby the new band and the situation repeats (approx. t 3 = 0.55 s) until the temperature field over the wire again favors propagation of the bands over nucleation of additional ones (approx. since t 4 = 0.7 s). The same process occurs in real experiments; however, locations of initial nucleation sites is (in addition to the temperature field) also influenced by local irregularities of the wire, e.g. lower cross-section area, local inhomogeneity of material properties or internal residual stresses. Combination of such material, geometrical and temperature effects gives rise to seemingly "haphazard" locations of nucleation. In contrast, perfect material and geometry is considered the simulation, hence, only the heat transfer influences band spacing. One may conclude, that the model is able to reproduce qualitatively different kinetics of localized transformation at different strain rates thanks to heat transfer effects.

Remark 6.
As shown in the work [63], localized transformation fronts in NiTi SMA wires exhibits quite a fine morphology with complex, inhomogeneous spatial distribution of stress and phases within the cross-section. Hence, the aim of presented one-dimensional simulations is mainly to demonstrate qualitatively the influence of thermal coupling on occurrence of phase fronts in the wire. However, mathematical treatment in Section 3 allows for a direct extension to three-dimensional situations in future. To this end, the research on further refinement of free energy contributions related to strain-softening, which would cover experimental peculiarities summarized in [65], is underway [28].