Multiscale analysis of the effect of interfacial energy on non-monotonic stress-strain response in shape memory alloys

The effect of formation and evolution of stress-induced martensitic microstructures on macroscopic mechanical properties of shape memory alloys in the pseudoelastic regime is investigated with account for size-dependent energy of interfaces. A quantitative relationship is established between the changes in free energy and dissipation on the interfaces at three microstructural scales and the overall mechanical characteristic of the material under tensile loading. The multiscale analysis carried out for a polycrystalline NiTi shape memory alloy has revealed that the interfacial energy storage and dissipation can strongly affect the shape and width of the stress–strain hysteresis loop. The predicted non-monotonic stress–strain response for the material of a selected grain size shows a remarkable similarity to the experimental one extracted from a tensile test of a laminate by Hallai and Kyriakides (2013). By the classical Maxwell construction, the non-monotonic response for a material element results in a commonly observed stress plateau for a tensile specimen, which is associated with the propagation of phase transformation fronts. This behaviour is confirmed with striking accuracy by 3D finite-element computations performed for a macroscopic tensile specimen, in which propagating instability bands are treated explicitly. © 2020 Institute of Fundamental Technological Research PAS. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
The unique behaviour of shape memory alloys (SMAs), notably the shape memory effect and pseudoelasticity, results from the martensitic phase transformation ( Otsuka and Wayman, 1998 ). The transformation is reversible in the sense that the overall inelastic strain can be induced and fully recovered upon adequate thermomechanical loading. This is accompanied by formation and evolution of martensitic microstructures at several spatial scales so that the martensitic transformation in SMA polycrystals is a truly multi-scale phenomenon. The present paper deals with quasistatic, stress-induced and nearly isothermal transformation in the pseudoelastic regime, with the emphasis put on rarely studied effects of the interfacial energy and dissipation at lower-scale interfaces.
The overall behaviour of SMAs is described by numerous phenomenological models (e.g., Raniecki et al., 1992;Auricchio and Petrini, 2002;Lagoudas et al., 2012;Sedlák et al., 2012;Stupkiewicz and Petryk, 2013;Qiao and Radovitzky, 2016;Jiang et al., 2016 ), and many others, where basic properties of SMAs are assumed at the macroscopic level. However, several important fea- JID: SAS [m5G;April 28, 2020;20:1 ] macroscopic volume element depend on the characteristic dimensions of microstructural elements, and thus the overall energy balance includes size-dependent contributions. As a result, the interfacial energy can govern the overall behaviour of SMAs in a size-dependent manner, see  for a general theory and  for its application to a simplified case of a rank-three laminate. The conclusion that the interfacial energy can significantly influence the overall response of SMAs can be traced back to Müller and Xu (1991) , but the respective macroscopic free energy factor has remained unknown. The aim of this work is to develop a hierarchical multi-scale framework for predictive modelling of the interfacial energy effects in SMA polycrystals. The specific focus of this work is on pseudoelastic behaviour of SMA polycrystals under quasi-static tension, with the restriction to the stress-induced transformation and isothermal conditions. In this regime, the material is initially in the austenitic state, and application of a mechanical load, after an initial stage of elastic response, induces the martensitic transformation that is accompanied by an overall inelastic strain. Upon unloading, the inelastic strain vanishes, and the initial undeformed state is recovered. The pseudoelastic response exhibits then two characteristic features, namely the hysteresis and frequently the stress plateau on the nominal stress-elongation diagram. The origin of the plateau and hysteresis has been widely discussed in the literature, but apparently not yet fully clarified. The present paper is intended to make a progress in this direction, although with the limitation of the analysis to the tensile loading/unloading. Nowadays, the aforementioned stress-elongation plateau is, or should be, no longer understood as a physical property of the material, but rather as the effect of propagating instability due to a negative slope on a part of the stress-strain diagram for the material . This has been clearly shown in the experiment by Hallai and Kyriakides (2013) , in which the intrinsic material response of NiTi under tension has been extracted from a laminate of a NiTi strip and steel face-strips that practically enforce a uniform deformation of the NiTi strip. Such a uniform deformation cannot be observed for a free-standing strip because the negative slope of the uniaxial stress-strain diagram leads to path instability in the form of strain localization. The transformation proceeds then through propagation of macroscopic transformation fronts that separate the transformed and untransformed domains. The fronts propagate at an approximately constant load under end-displacement control, hence a plateau is observed on the load-elongation diagram which should not be interpreted as a stress-strain diagram for the material. Fig. 1 reproduced from ( Hallai and Kyriakides, 2013 ) shows both the intrinsic material response and the specimen response.

ARTICLE IN PRESS
The former exhibits a negative slope on a significant part of the stress-strain diagram, the latter shows a load plateau.
The scenario discussed above (from propagating instability to load plateau) is in fact commonly observed in experiment. Strain localization and propagation of Lüders-like bands have been observed in NiTi strips and wires under tension ( Shaw and Kyriakides, 1997;Pieczyska et al., 20 06;Churchill et al., 20 09;Zhang et al., 2010;Sedmák et al., 2016 ), and instabilities leading to complex transformation patterns have been observed in NiTi tubes under tension and combined tension-torsion ( Sun and Li, 2002;Reedlunn et al., 2014;Bechle and Kyriakides, 2016 ) and under bending ( Bechle and Kyriakides, 2014;Reedlunn et al., 2014 ). The stress plateau has also been observed in numerous experiments without direct characterization of strain localization.
The stress plateau during the forward transformation under tension is often accompanied by a lower stress plateau during the reverse transformation upon unloading, with the related rateindependent hysteresis (e.g., Shaw and Kyriakides, 1997;Churchill et al., 2009;Bechle and Kyriakides, 2014;Reedlunn et al., 2014 ). This suggests that the respective material stress-strain response is also non-monotonic, however, it has never been determined experimentally.
From a mathematical point of view, the effect of a nonmonotonic local response on macroscopic behaviour in 1D setting has been widely studied in the literature, following Ericksen (1975) . Several macroscopic models of phenomenological type have also been developed for SMAs and implemented numerically in order to simulate the related phenomena in 1D setting ( Chang et al., 2006;Alessi and Bernardini, 2015;León Baldelli et al., 2015;Rezaee-Hajidehi and Stupkiewicz, 2018 ) and in 2D/3D setting ( Duval et al., 2011;Jiang et al., 2017a;2017b;Badnava et al., 2018;Rezaee-Hajidehi et al., 2020 ). A simplified form of the non-monotonic local response is then usually adopted, often a piecewise-linear one, although the actual response may exhibit a visible non-linearity, cf. Fig. 1 . Note that a non-monotonic local response results also from phase-field approaches applicable to SMA single crystals (e.g. Esfahani et al., 2018;T ůma et al., 2018 ). The models mentioned above involve some kind of regularization (gradient enhancement, non-local formulation, or viscous effects) that introduces a characteristic length-scale into the model and improves the robustness of the corresponding computational schemes. As an exception, the only regularization in the model of Jiang et al. (2017a,b) is that due to the natural 3D effects at the macroscopic transformation fronts.
In this paper, a multiscale framework for a quantitative analysis of the effect of interfacial energy on non-monotonic stressstrain response of pseudoelastic SMA polycrystals is developed, apparently for the first time. The framework involves three essentially independent models, two of them are micromechanical models themselves, and the multi-scale coupling between the models is performed in a hierarchical manner, i.e. the results of the model at a lower scale are transferred to the model at a higher scale in the form of constitutive functions or parameters that are identified using the lower-scale data. The multiscale framework is schematically shown in Fig. 2 , and the structure of the paper is presented below. The whole approach is built on the previous extensive work by the authors, which is referred to in the sequel. Each specific model used in this work relies on the incremental energy minimization, and the related framework that includes the interfacial energy effects is briefly recalled in Section 2 .
At the lowest scale, we apply the model of a spherical subgrain with a rank-two laminate of twinned martensite plates within the austenite matrix to determine the size-dependent interfacial energy contributions coming from the interfaces at three scales: at the twin boundaries, at the austenite-martensite interfaces, and at the boundary of the sub-grain. The model is adopted from ( Stupkiewicz and Petryk, 2010b ) and is briefly presented in Section 3.1 ; the basic equations are provided in Appendix A . In this model, the results of yet another micromechanical model (which is not indicated in Fig. 2 nor discussed in detail in the paper) are utilized, namely the estimate of the elastic micro-strain energy at the austenite-twinned martensite interface in NiTi is taken from ( Stupkiewicz et al., 2012 ). The outcome of the sub-grain model, i.e. closed-form expressions for the size-dependent interfacial energy contributions, constitutes then the input to the model of a SMA polycrystal to be considered at a higher scale. At the intermediate scale, the overall response of a SMA polycrystal is obtained by applying a respective multiscale model that is enhanced with the interfacial energy effects. This is presented in Section 4 . First, in Section 4.1 , the interfacial energy effects are introduced to the model of a single grain, which is done in a quite general setting without specifying the explicit form of the bulk free energy and dissipation contributions. The interfacial energy contributions are evaluated for a representative spherical sub-grain, and the respective results are taken from Section 3.1 . This general framework is then combined with a multiscale model of a SMA polycrystal. Here, the bi-crystal aggregate model ( Stupkiewicz and Petryk, 2010a ) is employed, which is briefly introduced in Section 4.2 and the equations are provided in Appendix B . Application of the complete scheme to polycrystalline NiTi under uniaxial tension is then reported in Section 4.3 . The predicted stress-strain response is non-monotonic during both loading and unloading and exhibits a significant dependence on the grain size, which is due to the interfacial energy effects.
Considering that the non-monotonic stress-strain response is expected to cause the path instability and strain localization, we finally carry out 3D finite-element simulations of the tension test in order to study the related effects quantitatively. To this end, in Section 5 , the macroscopic finite-strain model of pseudoelasticity ( Stupkiewicz and Petryk, 2013 ) is extended in order to accurately represent the intrinsic material response predicted by the micromechanical model in Section 4 . Specifically, by following the hierarchical multi-scale approach, the dissipation and free energy functions are fitted to the respective data resulting from the lowerscale analysis. To efficiently treat the instabilities and nonuniform transformation patterns, the finite-element implementation of the model employs the gradient regularization that has been recently developed by Rezaee-Hajidehi et al. (2020) . The finite-element simulations deliver the expected load-elongation response with two plateaus, as well as the transformation pattern in the form of a growing martensite band with inclined macroscopic transformation fronts.

Incremental energy minimization framework
It is well known that the stress-induced martensitic transformation in SMAs in the pseudoelastic regime proceeds on the microscopic scale by formation and evolution of laminate microstructures composed of crystallographic variants of martensite. In view of a large number of possible combinations of different variants of martensite, the micromechanical modelling of SMAs requires a method of selection of their physically favourable combination for a single grain. Under the assumption of isothermal conditions and an additional symmetry restriction imposed on dissipation, this is offered by the incremental energy minimization method ( Petryk, 2003 ) extended to the interfacial contributions to free energy and dissipation by  . On the time scale adopted in the rate-independent modeling, a negative increment of the interfacial energy, as in Eq. (3) below, is locally related to the annihilation of interfaces. This is associated with a sudden release of the interfacial energy that can hardly be reverted into the bulk free energy. Therefore, in the model it contributes predominantly to the rate-independent dissipation .
In outline, the incremental energy minimization framework with interfacial energy effects on different scales indexed by s can be summarized as follows : (i) Split of the total free energy into bulk and interfacial contributions,  (ii) Split of the dissipation increment into bulk and interfacial contributions, (iii) Negative increments of interfacial energy contribute to dissipation, (iv) Microstructure evolution is determined by succesive minimization of the incremental energy supply, The aforementioned scheme can be applied at any level of a hierarchical model, either separately or jointly for all levels of a polycrystalline aggregate. The theory is described by Petryk (2003) and  . Further details of the micromechanical model are discussed in the next sections.

Compound interfacial energy in a spherical sub-grain
3.1. Interfacial energy contributions at three scales (Stupkiewicz and Petryk, 2010b) The model developed by Stupkiewicz and Petryk (2010b) is employed here in order to estimate the interfacial energy contributions in a spherical domain. With reference to the micromechanical model of a polycrystalline aggregate discussed in Section 4 , consider a spherical domain of diameter d that is identified with a sub-grain , i.e. a part of a grain occupied by a single family of internally twinned martensite plates (habit-plane variants) within the austenite matrix, see Fig. 3 . The microstructure is assumed sufficiently fine so that, in view of the assumed periodicity, it is fully characterized by specifying the volume fraction of martensite η, the twin spacing h tw and the plate thickness M (or equivalently the plate spacing H = M/η).
Following Stupkiewicz and Petryk (2010b) , by exploiting the assumptions of separation of scales and periodicity of the microstructure, the total free energy density (per unit volume of the sub-grain) is decomposed into bulk and interfacial contributions, where ε is the average strain in the sub-grain in the smalldeformation setting. We note that the entire interfacial contribution φ i does not depend on the strain ε . Assumption (3) , which we adopt here, means that interfacial energy release is not fully converted into the bulk elastic energy but is at least partially dissipated. The interfacial contribution to dissipation is expressed through the negative increments of the interfacial contributions to the free energy, and thus it does not depend on the strain either. Accordingly, by applying the incremental energy minimization (4) , the evolution of the microstructural length parameters (here, M and h tw ) along a transformation path parameterized by η can be determined by minimizing the interfacial energy contributions independently of the bulk part. The corresponding model developed for the spherical sub-grain under consideration is briefly described below, and the set of respective equations is provided in Appendix A . For more details the reader is referred to Stupkiewicz and Petryk (2010b) . The analyzed microstructure, cf. Fig. 3 , involves three types of interfaces, namely the twin boundaries, the austenite-martensite interfaces and the sub-grain boundary, each associated with the corresponding spacial scale indicated by the subscript s = tw, am, gb . The interfacial energy contributions from the three scales are specified by the respective densities φ i s per unit volume of the sub-grain, cf. Eqs. (A .2) -(A .4) .
Two sources of interfacial energy are accounted for, namely the atomic-scale energy of phase boundaries and the energy of elastic micro-strains. The atomic-scale interfacial energy is considered at the twin boundaries and at the austenite-martensite interfaces, and the corresponding energy densities per unit area are denoted by γ a tw and γ a am , respectively. The transformation-induced changes in the atomic-scale energy of grain (sub-grain) boundaries are assumed negligible, and hence the respective atomic-scale interfacial energy is not considered.
The energy of elastic micro-strains results from local incompatibility of transformation strains. At the austenite-twinned martensite interfaces, the individual martensite variants are usually not compatible with the austenite, while compatibility in the average sense can be achieved due to twinning ( Bhattacharya, 2003 ). The energy of the elastic strains that accommodate the local incompatibility, when integrated over the volume and referred to the nominal interface area, is interpreted as an interfacial energy denoted by γ e am . Note that this energy is not a material parameter. It can be shown that γ e am is proportional to the twin spacing h tw , thus γ e am = e am h tw , where e am is a size-independent energy factor that can be determined by energy minimization for each representative microstructure along with the corresponding corrugated interface shape. This approach has been developed in a series of papers cited in ( Stupkiewicz et al., 2012 ), where specific numerical results for representative microstructures for a NiTi alloy can be found. Elastic micro-strain energy is also associated with the martensite plates that terminate at the sub-grain boundary, and the related energy At the initial stage of transformation ( η ≈ 0), nucleation of the first martensitic plates is associated with the increase of the interfacial energy, and hence the interfacial contribution to dissipation is then equal to zero, cf. Eq. (3) . Energetically optimal microstructural parameters h 0 tw and M 0 , cf. Eq. (A.5) , corresponding to the initial stage of transformation are then obtained by minimizing the free energy alone, see Eq. (A.6) for the respective formula for the compound interfacial energy at η ≈ 0. Note the characteristic square-root scaling rule for the parameters h 0 tw and M 0 in Eq. (A.5) with the characteristic lengths am and gb defined in terms of the interfacial energy parameters, which is typical for rank-one laminates ( Khachaturyan, 1983 ). The subsequent evolution of the plate thickness M = ηH is obtained by minimizing φ i for given η, while the twin spacing is kept constant, h tw = h 0 tw , cf. Eqs. (A .7) -(A .8) .
The predictive capabilities of the approach described above have been confirmed, both qualitatively and quantitatively, by the phase-field computations carried out by T ůma et al. (2016) with no a-priori assumptions concerning the microstructure. Specifically, the problem of a cylindrical grain with an austenite-twinned martensite microstructure has been analyzed, and it has been shown that the phase-field computations agree reasonably well with the square-root scaling of M , as predicted by the present approach, although with a somewhat different coefficient.

Application to NiTi alloy
The procedure developed by Stupkiewicz and Petryk (2010b) and summarized in Section 3.1 and Appendix A is now applied to a spherical sub-grain of a NiTi alloy. The computations are here performed for three sub-grain sizes d = 13 , 38 , 113 μm that correspond to the three grain diameters d gr = 20 , 60 , 180 μm considered in Section 4.3 .
To estimate the elastic micro-strain energy factor e am , we use the results of the computations carried out by Stupkiewicz et al. (2012) for the eight crystallographically distinct microstructures at the austenite-twinned martensite interface in NiTi, and we adopt the smallest value of this factor predicted for microstructure C-I-2, namely e am = 9 . 7 MJ/m 3 , which is close to the earlier simple estimate in Stupkiewicz and Petryk (2010b) . The estimate of the elastic micro-strain energy at the sub-grain boundary is given by Eq. (A.4) and results from a simplified micromechanical model which, unlike the model used to estimate e am , assumes elastic isotropy and small-strain kinematics. The corresponding material parameters are the elastic shear modulus μ = 26 . 3 GPa, Poisson's ratio ν = 0 . 33 and shape-strain vector modulus b = 0 . 1343 , as in ( Stupkiewicz and Petryk, 2010b ). Likewise, following Stupkiewicz and Petryk (2010b) , the atomic-scale interfacial energies are assumed equal to γ a tw = 0 . 014 J/m 2 and γ a am = 0 . 3 J/m 2 . Parameters κ s are assumed to be close to unity, specifically, κ tw = 1 , κ am = 0 . 95 , and κ gb = 0 . 9 .
In Fig. 4 , the results of computations are presented in terms of the evolution of plate thickness M during a complete forward and reverse transformation cycle and in terms of the number of plates per sub-grain, treated here as a continuous variable equal to d / H . The predicted microstructural parameters at the initial stage of transformation, cf. Eq. (A.5) , are provided in Table 1 .
The closed-form formulae for the coefficients φ max tw and φ max am+gb have been derived by Stupkiewicz and Petryk (2010b) , and the values of those parameters corresponding to the three sub-grain diameters considered are provided in Table 1 .
The total interfacial contribution to the thermodynamic driving force for phase transformation is It follows that for the austenite-to-martensite transformation during loading, the transformation driving force f i increases linearly with η, which corresponds to the transformation stress decreasing linearly with η. The latter property is frequently introduced as a constitutive assumption into phenomenological models of pseudoelastic transformation, with the proportionality factor taken arbitrarily. In contrast, the proportionality factor here takes the value expressed by the closed-form formula (7) derived from the interfacial energy considerations alone. This remarkable conclusion will have its consequences below when studying a more complex behaviour of a polycrystalline aggregate.

Interfacial energy contributions for a single grain
A typical scheme of micromechanical modelling of a SMA polycrystal involves a certain model of overall behaviour of a single grain that is combined with a suitable micro-macro transition from the scale of a single grain to that of a polycrystal. Representative examples include ( Thamburaja and Anand, 2001;Hackl and Heinen, 2008;Sengupta and Papadopoulos, 2009;Stupkiewicz and Petryk, 2010a;Xiao et al., 2019 ). The variables involved in the constitutive description of a single grain, in addition to the overall stress and strain, include the volume fractions of martensite variants, usually represented by the so-called habit-plane variants (HPVs) that account for the compatibility of average transformation strains at austenite-martensite interfaces ( Bhattacharya, 2003 ). Such a model does not involve any characteristic dimension of the microstructure and thus is incapable of describing size effects. Our aim here is to develop a model that, in a simplified manner, accounts for the size-dependent interfacial energy effects by exploiting the results of the model described in Section 3 .
Assume that the size-independent bulk contribution to the overall Helmholtz free energy density of a single grain in a polycrystalline aggregate depends on the average strain ε in the grain and on the volume fractions η k , k = 1 , . . . , N, of martensite variants (HPVs) in the grain, thus φv = φv ( ε , η k ) . As a rule, only several η k will be nonzero in a given grain. We further assume that the microstructure that develops within the grain during the stressinduced martensitic transformation divides the grain into several sub-grains such that each sub-grain is occupied by one family of parallel, internally twinned martensite plates corresponding to one HPV. Since individual volumes of the sub-grains can hardly be determined, a simplifying assumption is introduced that the interfacial energy contributions are only estimated for a representative spherical sub-grain , of diameter d being a specified fraction of a diameter d gr of the actual grain. Consequently, the volume fraction of the parallel martensite plates within the representative subgrain is assumed equal to the average volume fraction of martensite within the whole grain, η = k η k . As a result, the total free energy density per unit volume of the grain can be written in the following form where the interfacial energy contributions φ i s are determined using the model described in Section 3 , with the sub-grain diameter d as the input parameter. Recall that the model delivers the evolution of microstructural parameters and the related interfacial energies along a transformation path parameterized by the volume fraction of martensite, and thus each individual interfacial energy contribution φ i s can be expressed as a function of η. The complete evolution problem for a grain is now again formulated as the minimization problem for the incremental energy supply E , cf. Section 2 . In a strain-controlled process, the overall strain ε is prescribed and minimization is performed with respect to the volume fractions η k . The incremental minimization problem then reads By a similar argument, during the reverse transformation ( η < 0 ), we have  . 6. Bi-crystal aggregate model: a polycrystal is treated as an aggregate of bicrystals, and the corresponding additional level is introduced into the typical sequential averaging scheme ( Stupkiewicz and Petryk, 2010a ).
Since minimization (10) does not involve evolution of dimensional quantities, the approximate explicit formulae (6) for φ i tw and φ i am + φ i gb can be used with η in place of η, thus leading to a particularly simple formulation. The model of a grain, as presented above, can now be combined with a suitable micromechanical model of a SMA polycrystal.

4.2.
Bi-crystal aggregate model (Stupkiewicz and Petryk, 2010a) The transition from the level of a crystallographic lattice up to the scale of a SMA polycrystal undergoing stress-induced martensitic transformation is in this work performed using the bi-crystal aggregate model ( Stupkiewicz and Petryk, 2010a ). This multiscale model is illustrated graphically in Fig. 6 , and the underlying basic equations are provided in Appendix B . The model is presented here only in outline; more details can be found in the reference, along with the discussion of the assumptions and results, and comparison with experiment.
A distinctive feature of the model is that the scale transition from the scale of individual grains to that of a polycrystal is performed by introducing an intermediate scale of bi-crystals, cf. Fig. 6 . The polycrystal is thus treated as an aggregate of randomly oriented bi-crystals, each bi-crystal being formed by two neighbouring grains separated by a planar interface (grain boundary). Accordingly, grain interaction is accounted for by enforcing the compatibility conditions (B.5) on the average strains and stresses of each grain of the bi-crystal. Introduction of the intermediate scale of bi-crystals and consideration of the respective grain interaction mechanism significantly improves, as shown in Stupkiewicz and Petryk (2010a) , the performance of the simple Taylor (Voigt) and Sachs (Reuss) grain-to-polycrystal scale transition schemes when applied as the transition scheme between the bi-crystals and the polycrystal. As a result, those simple transition schemes can be used with a higher confidence, and specifically the Taylor scheme is employed in this work.
A characteristic feature of the model is the quantitative description of the transformation at all scales using only a small number of parameters that remain to be assumed arbitrarily, which will allow us to predict the changes in the free energy and dissipation during transformation. This is done by step-by-step minimization of the incremental energy supply with respect to the variables at all levels of the description, cf. Appendix B .

The effect of interfacial energy on non-monotonic stress-strain characteristic
A multiscale model of a SMA polycrystal with interfacial energy effects is now obtained in a straightforward manner by includ-ing in the bi-crystal model the scheme developed in Section 4.1 . Specifically, an enhanced model corresponding to the level of a single grain is obtained by combining the bulk contribution φv to the free energy, cf. Eq. (B.2) , and the bulk contribution D v to the dissipation function, cf. Eq. (B.4) , with the respective interfacial energy contributions according to the general scheme outlined in Eqs. (9) -(12) . The resulting single-grain model enhanced with the interfacial energy effects can then be seamlessly employed in the overall scheme of the bi-crystal model. Importantly, the quadraticprogramming structure of the actual computational scheme is preserved when the closed-form formulae (6) are used for φ i tw and φ i am + φ i gb . The results reported below have been obtained for a NiTi polycrystal for three values of an average grain diameter, d gr = 20 , 60 , 180 μm. Following Stupkiewicz and Petryk (2010b) , it has been assumed that the diameter of a representative spherical sub-grain is such that its volume is equal to one fourth of the volume of the grain, thus d = 4 −1 / 3 d gr . This choice is justified by the analysis of the results of the bi-crystal aggregate model ( Stupkiewicz and Petryk, 2010a ) showing that typically three to five martensite variants (HPV's) appear within one grain during proportional loading. Accordingly, the estimates of the interfacial energy contributions have been computed in Section 3.2 for three sub-grain diameters d = 13 , 38 , 113 μm. The respective material parameters used as the input to those computations are provided in Section 3.2 , and the resulting interfacial energy contributions, which serve as the input to the present computations, are provided in Table 1 .
The remaining material parameters pertinent to the bicrystal model are the following. The elastic properties are as in Section 3.2 , and the chemical energy, which is correlated with the temperature and influences the level of the transformation stress, is assumed equal to am φ 0 = 19 MPa. The bulk contribution to dissipation is specified by the critical driving force for transformation f c = 3 MPa and by the critical driving force for martensite reorientation f r = f c / 2 . The transformation strains ε t k of martensite variants (HPV's) of NiTi are determined using the crystallographic theory of martensite ( Bhattacharya, 2003 ), and N = 48 variants are used in the present (small-strain) computations of the uniaxial tension test, see Stupkiewicz and Petryk (2010a) for details. A drawing texture is assumed with the 111 poles of the cubic austenite preferably aligned with the drawing direction (along which tension is applied) such that the angle between the 111 poles and the drawing axis is randomly chosen between 0 and 30 • . Following the procedure developed by Stupkiewicz and Petryk (2010a) , the total of 256 distinct grain orientations has been generated according to the assumed drawing texture. The grains have been arranged in pairs in a random manner, and the orientation of the corresponding grain boundary has also been generated in a random manner. The relative volume fraction of grains in each bi-crystal has been assumed equal to 0.5. Overall, the bi-crystal averaging has been performed for an aggregate of 128 randomly generated bi-crystals. It has been checked that this number is sufficient, and further increase of the number of bi-crystals does not influence the results visibly. Fig. 7 a shows the predicted pseudoelastic response of a representative element of a NiTi polycrystal under uniaxial tension. As expected, due to including the interfacial energy effects, the response exhibits a significant dependence on the grain size. The dashed line in Fig. 7 a indicates the underlying response predicted by accounting only for the bulk contributions to the free energy and dissipation, i.e. in the limit of vanishing interfacial energy contributions for d gr → ∞ .
Several effects can be observed in Fig. 7 a. With decreasing grain size, the initiation of transformation is shifted to a higher stress. The stress must increase to compensate the extra thermodynamic driving force associated with the interfacial contribution to the free energy, f i = −∂ φ i / ∂ η. At the initial stage of transformation, the interfacial contribution to the free energy increases, cf. Fig. 5 and formula (8) , hence f i is then negative. The initially negative driving force f i increases during transformation, which leads to a decrease of the stress. Averaging over the polycrystalline aggregate causes the stress to vary nonlinearly with the overall tensile strain. At the later stage, the further decrease of the stress is suspended because the negative increments in interfacial free energy are dissipated, and ultimately the stress-strain slope becomes strongly positive in the final stage of transformation. The final slope tends to the elastic modulus since martensite detwinning has been excluded form considerations. This provides an explanation, coming from the interfacialenergy considerations, why the stress-strain relationship during the forward transformation is non-monotonic. For the grain size small enough, a significant branch of negative slope is observed, in agreement with the experimental relationship quoted after Hallai and Kyriakides (2013) in Fig. 1 . The maximum stress and the subsequent stress drop increase with decreasing grain size. Interestingly, the predicted shape of the upper curve for grain size d gr = 20 μm almost coincides with the experimental curve for the material shown in Fig. 1 after Hallai and Kyriakides (2013) ; unfortunately, the underlying average grain size of the material tested was not provided in the reference. The response during the reverse transformation upon unloading is similar to that during the forward transformation, except that the stress is lower so that a hysteresis loop is observed. The area of the hysteresis loop, which is equal to the energy dissipated in the complete loading-unloading cycle (in isothermal conditions), also increases with decreasing grain size. This illustrates the important effect of the interfacial contribution to dissipation.
Note that the 'softening' response, as predicted by the present model, cannot be observed in a uniaxial tension experiment unless special means are taken to avoid the nonuniform deformation, as in the experiment of Hallai and Kyriakides (2013) . Indeed, it is commonly observed that the stress-induced transformation, particularly in NiTi in tension, proceeds through nucleation of localized transformation bands followed by propagation of macroscopic transformation fronts (e.g., Shaw and Kyriakides, 1997;Sun and Li, 2002;Zhang et al., 2010;Sedmák et al., 2016 ). This is typically accompanied by a stress plateau that is associated with propagation of the fronts and does not represent the material behaviour at a material point. Fig. 7 b shows hypothetical specimen responses corresponding to the intrinsic material responses shown in Fig. 7 a. The uniaxial stress-strain diagram with a negative slope leads to instability of a uniform deformation path and hence cannot be reproduced by the specimen ( Petryk, 20 0 0 ). Here, the classical Maxwell construction has been used. The plateau segment connects the monotonically increasing segments of the original curve such that the areas above and beneath the plateau and limited by the original stress-strain curve are equal. The plateau stress can be shown to be the lowest thermodynamically admissible one for propagation of the transformation zone in which the original stress-strain curve is traversed at every material point. Additionally, it has been assumed here that the hypothetical specimen response is formed by the loading branch of the original stress-strain curve up to the maximum stress, followed by an instanteneous stress drop, a stress plateau corresponding to the Maxwell stress, and the final hardening part towards the end of transformation. The unloading part of the response is constructed analogously. The maximum stress and the subsequent instanteneous stress drop are supposed to correspond to nucleation of the first martensite band. Similar features are commonly observed in experiments, although the actual stress overshoot is highly sensitive to imperfections, material inhomogeneity, etc. The minimum stress followed by a stress increase at the beginning of the reverse transformation would correspond to nucleation of the austenite band in a completely transformed specimen. If the reverse transformation does not require nucleation of the austenite band, but proceeds by a reverse propagation of an existing macroscopic transformation front, then the corresponding feature is not expected since no energy barrier must then be overcome and the plateau can be followed from the beginning.
It is not evident in advance that the hypothetical behaviour described above will agree with 3D simulations of the transformation process. To investigate this, a detailed analysis of strain localization and nonuniform transformation is carried out in the next section using a phenomenological finite-strain model of pseudoelasticity.  JID: SAS [m5G;April 28, 2020;20:1 ]

Macroscopic modelling of nonuniform phase transformation
The micromechanical model developed in the previous sections predicts the stress-strain response that exhibits a negative slope during stress-induced transformation. The related scenario of strain localization in the macroscopic specimen has then been analyzed qualitatively by employing the Maxwell construction, cf. Fig. 7 b. In this section, nonuniform phase transformation is analyzed quantitatively using a macroscopic phenomenological model of pseudoelasticity implemented in a finite-element code. The model is briefly described in Section 5.1 , its finite-element implementation is commented in Section 5.2 , and the results of finite-element computations are reported in Section 5.3 .

Phenomenological finite-strain model of pseudoelastic SMAs
The model used in this work is an extension of the model developed recently by Rezaee-Hajidehi et al. (2020) by introducing a gradient regularization to the finite-strain 3D model of pseudoelasticity of Stupkiewicz and Petryk (2013) . The gradient regularization has been introduced in order to adequately treat strain localization and propagation of macroscopic transformation fronts. The model is formulated by specifying the free energy function and the dissipation function, and the complete evolution problem is then formulated within the incremental energy minimization framework. The reader is referred to Stupkiewicz and Petryk (2013) for a detailed presentation of the constitutive model and to Rezaee-Hajidehi et al. (2020) for the aspects related to the gradient regularization. Below we describe very briefly the extensions introduced to the model in order to accurately represent the intrinsic material response predicted by the micromechanical model.
The modifications introduced in the present paper concern both the free energy function and the dissipation function. In general, the two functions cannot be uniquely fitted using the stress-strain curve alone. To avoid ambiguity, the contributions coming from the free energy and dissipation functions have been separated and fitted individually by adequate postprocessing of the results of the micromechanical model.
The modification introduced to the free energy function concerns the interaction energy term φ int ( η), which in the original macroscopic model ( Stupkiewicz and Petryk, 2013 ) was adopted as a quadratic function of η, the volume fraction of martensite, with a constant coefficient playing the role of a hardening or softening modulus. Actually, the interaction energy influences the stressstrain response only through its derivative φ int (η) . It is thus convenient to perform the fitting directly for the derivative φ int (η) , rather than for the function φ int ( η) itself, and this approach has been adopted here. The rate-independent dissipation function in the macroscopic model is assumed as ˆ where f c > 0 is the critical thermodynamic driving force for phase transformation. The modification introduced here amounts to specifying the critical driving force separately for the forward and reverse transformation, and basically consists in introducing the dependence of the critical driving force on η, thus with the meaning of angular brackets as in Eq. (3) .
Each of the constitutive functions φ int (η) , f + c (η) and f − c (η) has been fitted using a Bernstein polynomial of degree 12. The impact of the individual contributions on the pseudoelastic stress-strain response is illustrated in Fig. 8 a for the case of the grain diameter d gr = 20 μm. The adopted procedure yields a good fitting of the response predicted by the micromechanical model, as shown in Fig. 8 b for the three grain diameters examined in Section 4.3 .
As a part of the identification procedure, the transformation strain in tension has been determined as equal to ε T = 0 . 055 .

Finite-element model
The finite-element formulation of the gradient-enhanced model relies on the micromorphic-type regularization, as proposed by Rezaee-Hajidehi and Stupkiewicz (2018) , see also Mazière and Forest (2015) . The computer implementation follows exactly that developed by Rezaee-Hajidehi et al. (2020) , where all the details can be found. The resulting computational model involves, as the global unknown fields, the displacement field, the micromorphic counterpart of the volume fraction of martensite, and temperature. The thermomechanical coupling is introduced to the model because it provides an additional, physically-based regularization of the non-monotonic response. Indeed, the thermomechanical coupling improves the robustness of the computational model, which is particularly beneficial when the softening is high. The thermomechanical coupling renders the response rate-dependent so that the experimentally observed loading-rate effects can also be studied (e.g., Rezaee-Hajidehi et al., 2020 ), which, however, is not pursued here.
In the present computations, isoparametric hexahedral elements are used with triquadratic shape functions for the displacement and with trilinear shape functions for the remaining unknowns. The computer implementation is carried out using AceGen , a symbolic code generation system, and the finite-element computations are performed using AceFEM , a finite-element code interfaced with AceGen ( Korelc and Wriggers, 2016 ).
Finite-element computations have been carried out for a dogbone specimen loaded in tension. The specimen of the total length of 60 mm contains the gauge segment of the length L = 24 mm and uniform cross-section of 5.2 × 0.15 mm 2 . The specimen width gradually increases in the transition segments between the gauge segment and the gripping segments, the latter of the width of 20 mm. The axial displacement is prescribed at the specimen ends in a way that no bending moment is transmitted to the specimen by allowing the gripping segments to rotate. The load is applied with a low constant nominal loading rate of 10 −5 s −1 in order to obtain a nearly isothermal response. The elongation δ/ L is determined as the averaged relative axial displacement of the two crosssections at the ends of the gauge segment normalized by the gauge length L . A small geometrical imperfection is introduced in the upper-right part of the gauge segment to trigger a non-symmetric mode of strain localization.
The finite-element mesh within the gauge length is adopted such that the in-plane element size is 0.13 mm, and through-thethickness element size is 0.15 mm, so that the gauge segment is discretized into 185 × 40 × 1 elements. The in-plane element size gradually increases towards the gripping segments. The total number of degrees of freedom exceeds 550 0 0 0. For completeness, the model parameters used in the finite-element computations are provided in Appendix C . Fig. 9 shows the normalized force-elongation diagram computed for the grain diameter d gr = 20 μm. Also shown are the snapshots illustrating the evolution of the transformation pattern at selected instants, as indicated by the markers and labels on the force-elongation diagram. Each snapshot represents the specimen in the deformed configuration with the colour indicating the volume fraction of martensite.

Finite-element simulations of uniaxial tension test
Initially, the transformation is nearly uniform within the gauge segment, and this stage corresponds to the initial 'hardening'  branch of the intrinsic stress-strain relationship for the material.
At the elongation of about 1% a sudden load drop is observed which is associated with the nucleation of a thin inclined band of martensite (snapshot 1). The band nucleates at the imperfection. The transformation proceeds then through propagation of macroscopic transformation fronts at an approximately constant load. The obtained nominal stress-elongation curve is visibly smoother than the experimental one shown in Fig. 1 . The wiggles in the experimental curve may be related to abrupt events, like nucleation of new interfaces, changes in orientation of interfaces, and formation of criss-cross patterns, or to material inhomogeneity, e.g., due to the grain microstructure. In the present computations, the abrupt events are not observed during loading, but they are observed during unloading, thus leading to small wiggles at the initial stage of unloading, as commented below. On the other hand, material inhomogeneity is not included in the model, except for that resulting from the finite-element discretization. However, the finite-element mesh is here sufficiently fine, in particular, with respect to the thickness of the diffuse macroscopic transforma-tion fronts, so that the related nominal stress fluctuations are very small and thus not visible. The macroscopic transformation fronts are diffuse interfaces that separate the domains of low and high volume fraction of martensite equal to 0.037 and 0.94, respectively. Those values are close to the limit values of 0 and 1, and their deviation from the limit values is thus hardly visible. The end of the load plateau corresponds to the instant at which the martensite domain extends over the entire gauge length (snapshot 3). With further elongation, the load increases which is accompanied by a uniform transformation within the gauge segment until the transformation is completed.
Upon unloading, the reverse transformation proceeds through a reverse motion of the macroscopic transformation fronts. At the initial stage, small oscillations of the load are observed which are associated with the development and evolution of a crisscross pattern at the macroscopic transformation fronts (snapshot 5). Subsequently, two inclined interfaces form and propagate in a stable manner until the martensite band annihilates. The event of annihilation of the macroscopic transformation fronts is accompanied by a small increase of the load. It has been checked that the plateau stress, both during the forward and reverse transformation, is in excellent agreement with the Maxwell stress. This is illustrated for d gr = 20 μm in Fig. 10 in which the horizontal dashed line denotes the Maxwell stress, determined in the standard way by equating the shaded areas above and below the Maxwell line. It can be seen from Fig. 10 that the actual specimen response can be remarkably well approximated by the hypothetical specimen response obtained by employing the Maxwell construction, cf. Fig. 7 . A minor and expected difference between the two is that the stress minimum followed by a sudden stress increase at the beginning of the reverse transformation is not observed in the actual specimen response because the reverse transformation proceeds through a reverse motion of already existing macroscopic transformation fronts. Another small difference is that the nucleation of the martensite band starts slightly before the maximum stress point, and the nucleation stress is influenced by the imperfection. Moreover, the subsequent stress drop is characterized by a finite slope, which has been found to depend mostly on the specimen length, while an infinite slope is assumed in the hypothetical response.
The specimen responses computed for the three grain sizes are shown in Fig. 11 . The qualitative features are here the same regardless of the grain size. Also, it has been checked that in all cases the plateau stresses agree very well with the respective Maxwell stresses. As a result, the width of the pseudoelastic hysteresis loop shows a significant grain-size dependence which is inherited from the respective intrinsic material responses.

Conclusion
A multiscale analysis of pseudoelastic behavior of SMAs has been carried out, starting from crystallographic lattice rearrangements up to the scale of a polycrystalline specimen. The standard analysis of a rank-two laminate domain has been enhanced by including a hierarchy of interfacial energies on three scales: of planar twin interfaces, of corrugated interfaces between austenite and twinned martensite, and of domain boundaries. To determine the actual microstructure evolution among plenty of possibilities, the incremental energy minimization technique ( Petryk, 2003 ) has been applied. Upon subsequent passage to the scale of a polycrystalline aggregate, it has been shown, apparently for the first time, that the interfacial energy accumulation and release can be predominant in predicting a non-monotonic uniaxial stress-strain response and the accompanying hysteresis. The effects are dependent on the grain size, which has been examined quantitatively for a NiTi SMA with the conclusion that the smaller the grain size (in the relevant grain-size range) the stronger the calculated effect. Finally, the finite-element study of a tensile specimen in the finite-deformation setting has confirmed that the simple Maxwell construction for the non-monotonic material response in tension accurately predicts the plateau stress in loading and another one in unloading.
The simulated stress-elongation hysteresis loop with two plateaus for a NiTi specimen subjected to uniaxial tensile loading and unloading is both qualitatively and quantitatively close to that commonly observed in the quasi-static tension tests performed at a constant temperature corresponding to the pseudoelastic regime. Of course, the mean stress level depends significantly on the temperature, and can be straightforwardly fitted in the calculations by adjusting the value of chemical free energy. In turn, the calculated width of the hysteresis loop is found to depend on the grain size through the interfacial energy effects and also directly on the assumed value of the critical driving force f c that specifies the bulk contribution to dissipation. Validation of the present model in this respect is more difficult since experimental studies of grain-size dependence of the hysteresis loop are scarce. An increase of hysteresis with decreasing grain size has been observed experimentally in CuAlBe polycrystal ( Montecinos et al., 2008 ) and in CuAlNi microwires ( Chen and Schuh, 2011 ), and this qualitative effect is correctly reproduced by our model. On the other hand, the experimental data for the grain sizes within a nanometer range ( Ahadi and Sun, 2015;Sun and He, 2008 ) show an opposite effect. Note, however, that the present model deals with domains of micrometer size to allow formation of rank-two austenite-martensite laminates, while different transformation mechanisms may operate in materials with nanometer-sized grains ( Waitz et al., 2007 JID: SAS [m5G;April 28, 2020;20:1 ] same time, the shape of the non-monotonic stress-strain curve in loading predicted for the NiTi material of grain size d gr = 20 μm shows a remarkable similarity to the experimental one quoted in Fig. 1 after Hallai and Kyriakides (2013) . As mentioned above, that shape is obtained here solely from the interfacial energy considerations. All effort s have been made to leave the number of parameters of unverifiable value in the analysis as small as possible. In result, practically only one such parameter, f c , has remained. This has been done at the cost of a number of simplifying assumptions which have been explained in the text and in the cited references. These assumptions impose limitations on the applicability of the present micromechanical approach, for instance, to predominantly proportional loading paths. Further work in needed to clarify and possibly overcome such limitations in the future.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement
This work has been partially supported by the National Science Center ( NCN ) in Poland through Grant No. 2018/29/B/ST8/00729 .

Appendix A. Size-dependent microstructure evolution within a spherical sub-grain
In this appendix we provide the complete set of the governing equations of the model of microstructure evolution in a spherical sub-grain, cf. Fig. 3 . The model is briefly described in Section 3.1 , while the detailed derivation can be found in ( Stupkiewicz and Petryk, 2010b ).

A1. Interfacial energy at three scales
The total free energy density (per unit volume of the sub-grain) comprises the interfacial energy contributions from three scales, where, as a specification of Eq. (3) , the interfacial contribution to dissipation is equal to Since the interfacial energy release can hardly be expected to be converted to the elastic bulk energy and is thus likely to be dissipated, it is reasonable to take values of κ s closer to 1 than to 0.

Appendix B. Bi-crystal aggregate model
In this appendix we provide the complete set of the governing equations of the bi-crystal aggregate model. The model, formulated in the small-strain setting, is briefly commented in Section 4.2 , and the details can be found in ( Stupkiewicz and Petryk, 2010a ).

B1. Rank-two laminated sub-grain
Martensite is assumed to appear in the form of internallytwinned plates (i.e. habit-plane variants, HPVs) that form a ranktwo austenite-martensite laminate within the corresponding subgrain, cf. Figs. 3 and 6 . The overall transformation strains ε t k of the martensite plates are compatible with unstressed austenite, thus . . . , N, (B.1) where the habit-plane normal m k and shape-strain vector b k can be obtained from the classical crystallographic theory of martensite ( Bhattacharya, 2003 ).

B2. Single grain
The Helmholtz free energy of a grain composed of laminated sub-grains of uniform elastic properties is assumed in the following form, where η k denotes the volume fraction of martensite variant k , η is the total volume fraction of martensite in the grain, am φ 0 is the chemical energy of transformation, and L is the elastic stiffness tensor. The interaction energy φv int is here neglected, φv int = 0 , which corresponds to a constant-stress averaging scheme within a single grain. As usual, we have σ = ∂ φ / ∂ ε = L ( ε −ε t ) .
The rate-independent dissipation function for phase transformation ( f c > 0) and martensite reorientation ( f r > 0) is adopted in the following form, (B.4) The original model of Stupkiewicz and Petryk (2010a) does not account for the interfacial energy effects, thus the free energy φ and the dissipation function D comprise only the bulk contributions φv and D v , respectively.

B3. Bi-crystal composed of two adjacent grains
The polycrystalline material is considered as an aggregate of bicrystals , where each bi-crystal is composed of a pair of neighbouring grains separated by a planar grain boundary of orientation n . The following compatibility conditions at the grain boundary are imposed on the (average) strains and stresses in the grains, ε (2) −ε (1) = 1 2 (c n + n c ) , ( σ (2) −σ (1) ) n = 0 , where c is an unknown vector. The average strain and stress within the bi-crystal are given by ε b = ζε (1) + (1 − ζ ) ε (2) , σ b = ζσ (1) + (1 − ζ ) σ (2) , (B.6) and the average free energy and dissipation densities read where ζ and 1 − ζ denote the volume fractions of the grains in the bi-crystal, and ζ = 1 2 is assumed in the simulations. An analytical expression for φ b has been derived in ( Stupkiewicz and Petryk, 2010a ).

B4. Polycrystal as an aggregate of differently oriented bi-crystals
The macroscopic response of a SMA polycrystal is obtained by averaging the responses of bi-crystals of random grain orientations and random grain-boundary orientations, possibly characterized by a crystallographic texture. The respective averaging operation is denoted by { ·} so that the macroscopic strain and stress, and the macroscopic free energy and dissipation densities are given by To close the model, the Taylor averaging scheme is adopted so that the average strain in each bi-crystal is constrained to be equal to the macroscopic strain,

B5. Incremental energy minimization
For a strain-controlled process, the macroscopic response is obtained by minimizing the incremental energy supply, E = + D → min for prescribed E . (B.10) In view of the Taylor constraint (B.9) , the minimization problem can be solved separately for each bi-crystal, (B.11) The minimization problem (B.11) is a non-smooth minimization problem with 2 N unknown increments of the volume fractions of martensite variants in the bi-crystal. The problem can be transformed to a quadratic programming problem with 4(N + 1) unknowns, which can be efficiently solved using the interior-point method.