Localization of the stress-induced martensitic phase transformation plays an important role in the fatigue behavior of shape memory alloys (SMAs). The phenomenon of return-point memory that is observed during the subloop deformation of a partially-transformed SMA is a clear manifestation of the interaction between localized phase transformation and degradation of the functional properties. The present study aims to demonstrate this structure–material interaction in the modeling of return- point memory. It seems that this crucial aspect has been overlooked in previous modeling studies. For this purpose, we developed a gradient-enhanced model of pseudoelasticity that incorporates the degradation of functional properties in its constitutive description. The model is employed to reproduce the hierarchical return-point memory in a pseudoelastic NiTi wire under isothermal uniaxial tension with nested subloops. Additionally, a detailed analysis is carried out for a NiTi strip with more complex transformation pattern. Our study highlights the subtle morphological changes of phase transformation under different loading scenarios and the resulting implications for return-point memory. Keywords: Shape memory alloys; Phase transformation; Functional degradation; Propagating instabilities; Subloop deformation; Modeling 1. Introduction The practical interest in shape memory alloys (SMAs), especially NiTi, stems from their ability 1 to withstand and recover large strains. This ability is exhibited through mechanical loading and 2 unloading at sufficiently high temperatures (pseudoelasticity) or through mechanical loading and 3 unloading followed by heating (shape memory effect). The underlying mechanism is the crystallo- 4 graphically reversible martensitic phase transformation that occurs between the austenitic parent 5 phase (stable at higher temperatures, possessing higher crystal symmetry) and the martensitic prod- 6 uct phase (stable at lower temperatures, possessing lower crystal symmetry) [1]. By leveraging the 7 unique characteristics of SMAs, they have found a broad range of applications across various fields, 8 from micro-scale biomedical devices to macro-scale industrial components [2, 3]. The operational 9 ∗Corresponding author. Email addresses: mrezaee@ippt.pan.pl (Mohsen Rezaee-Hajidehi), maciej.rys@ncbj.gov.pl (Maciej Ry´ s) †Published in Int. J. Mech. Sci. , 2024, doi: 10.1016/j.ijmecsci.2024.109569 Preprint submitted to Elsevier July 21, 2024 lifespan of SMAs in most of the applications involves enduring cyclic mechanical/thermal loadings, 10 which highlights the great importance of identifying their fatigue behavior. It is well-recognized 11 that, due to the martensitic phase transformation, fatigue in SMAs is more complex than in common 12 engineering metals and is mainly classified into two aspects: degradation of functional properties 13 (such as recoverable strain, transformation stress, and hysteresis loop area), known as functional 14 fatigue, and the evolution of damage in the material, known as structural fatigue [4]. This complex 15 nature demands special attention and, as a result, has prompted a tremendous number of stud- 16 ies that focus on the fatigue characterization of SMAs from a variety of perspectives and on the 17 underlying micromechanical processes [3–13]. 18 Stress-induced martensitic transformation in pseudoelastic NiTi appears (typically, in tension- 19 dominated loadings) as localized instabilities in the form of martensite bands, and subsequently 20 progresses via propagation of the instabilities in the form of patterned interfaces (macroscopic trans- 21 formation fronts) that separate the domains of low-strained austenite and high-strained martensite, 22 e.g., [14–17]. Due to the high strain incompatibilities that exist within the transformation front 23 and the ensuing large local stresses, it can be reasonably inferred that propagating instabilities can 24 vitally influence both the functional fatigue and structural fatigue of the material. Despite the 25 longstanding recognition of this crucial aspect [4, 18–21], its direct validation was provided only a 26 few years ago in the experiments conducted by Zheng et al. [8, 22, 23]. It was demonstrated that 27 in view of the repetitive nucleation and propagation of the localized transformation in NiTi strips 28 under cyclic uniaxial tension, a rapid degradation of pseudoelasticity occurs that accelerates the 29 fatigue crack initiation and fatigue failure. 30 An interesting manifestation of the interaction between propagating instabilities and functional 31 degradation is found in the subloop deformation of a partially-transformed SMA specimen under 32 displacement-control loading. The subloop deformation has been extensively investigated experi- 33 mentally, notably for NiTi [18, 19, 23–26] but also for other SMAs [27, 28]. Fig. 1(a), reproduced 34 from Tobushi et al. [24], depicts the mechanical response of a NiTi wire subjected to subloop paths. 35 For a more intuitive description of the phenomenon, hypothetical schematics of the corresponding 36 transformation front evolution are provided in Fig. 1(b). As the front propagates along the wire, 37 it leaves behind transformation-induced microstructural defects, such as dislocations and stabilized 38 (locked-in) martensite. During the subloop unloading (for instance, the first subloop, which starts 39 at point A), the front travels backward over an already swept zone (from A to B), hence intensify- 40 ing the generated defects. Accordingly, during the subloop reloading, the propagation of the front 41 over the twice-swept zone occurs with a lower stress level compared to the original transformation 42 plateau. Upon entering the pristine zone (at point A), which is virtually free of transformation- 43 induced defects, the front experiences the transformation-onset stress characteristic to the initial 44 material state. This causes the stress to catch up with the original plateau by passing through the 45 subloop unloading point. This trait is known as the return-point memory. The process repeats in 46 the subsequent subloops and culminates in an intriguing hierarchical return-point memory. 47 Motivated by the experimental results, numerous attempts have been made to develop SMA 48 models capable of capturing the phenomenon of return-point memory during the subloop deforma- 49 tion, either through incorporating the permanent strain contribution and degradation of functional 50 properties [29–32] or by merely refining the constitutive equations of (non-cyclic) model of pseudoe- 51 2 (a) (b) 0 3 6 90175350525700 Axial strain [%]Axial stress [MPa ]A BC DE FG Figure 1: Return-point memory in NiTi wire subjected to uniaxial tension with three nested subloop paths: (a) the structural stress–strain response, and (b) hypothetical schematics of the corresponding transformation front evolution. The stress–strain response in panel (a) is reproduced from Tobushi et al. [24] (courtesy of R. Matsui). The red arrows in panel (b) indicate the trajectory of the front propagation, and the color scales quantify the recurrence of the front’s traversal over the wire’s segments. lasticity [25, 33, 34]. In fact, a physically-relevant approach for modeling the return-point memory 52 should hinge on the interaction between the propagating instabilities (structural inhomogeneities) 53 and the functional degradation of the material. Nevertheless, most of the existing models (including 54 those referenced above) postulate a homogeneous martensitic phase transformation, while address- 55 ing a problem with a transformation of localized nature. Albeit this simplifies the computations 56 significantly, it is not a plausible assumption in the present context. To the best of our knowledge, 57 the only related modeling study that has accounted for this structure–material interaction is the 1D 58 model of Bartel et al. [32]. In their model, however, instabilities do not originate from a softening- 59 type intrinsic material response but are rather treated as weak displacement discontinuities that 60 separate the transformed and untransformed material points (indeed, experiments, e.g., [35, 36], 61 have confirmed that the true intrinsic response of NiTi is characterized by a significant soften- 62 ing branch). It should be remarked that recently Xiao and Jiang [37, 38] have acknowledged this 63 structure–material interaction in their simulations, however, their applications did not specifically 64 pertain to the subloop deformation and return-point memory. 65 In light of the above premise, this work aims to provide a detailed analysis of the phenomenon 66 of return-point memory by accounting for the interaction between propagating instabilities and the 67 degradation of the functional properties of the material. To achieve this, a gradient-enhanced model 68 of pseudoelasticity with functional degradation is developed in this work. The model is formulated 69 within the small-strain theory. The basic structure of the model follows the non-gradient model of 70 pseudoelasticity developed by Stupkiewicz and Petryk [39] and is based on the energy minimization 71 principle. The gradient-enhancement, micromorphic regularization, and thermomechanical coupling 72 are adopted from our previously-developed gradient-enhanced model [40, 41]. This previous model 73 has demonstrated its capability in reproducing the complex patterns of phase transformation in 74 3 pseudoelastic NiTi specimens under uniaxial tension [41], including the effect of loading rate and 75 latent heat of transformation on martensite domain formation, and in pseudoelastic NiTi tubes 76 under combined tension–torsion [42, 43]. The main advancement of the model in the present work 77 compared to the previous version lies in the incorporation of permanent inelastic strain and the 78 enrichment of the constitutive equations with functional degradation effects. Consequently, given 79 its ability to treat localization effects via gradient-enhancement and micromorphic regularization, 80 the model can be considered a suitable tool for addressing problems where both cyclic loading and 81 transformation localization are at play. 82 In what follows, we first introduce the model in Section 2. The model is employed to analyze the 83 problem of subloop deformation in NiTi wire and strip under uniaxial tension. The corresponding 84 results are presented and discussed in Section 3. In addition, a simplified version of the model is 85 provided in Appendix A. 86 2. A small-strain model of pseudoelasticity with functional degradation The present model falls in the category of phenomenological models. Accordingly, the consti- 87 tutive relations are tailored, in a simple phenomenological manner, to mimic the pseudoelasticity 88 degradation effects. Since the focus of this study is on the analysis of the return-point memory, 89 which is relevant at the macroscopic scale, a phenomenological description seems to adequately 90 fulfill the intended purpose. In Section 2.1, we introduce the constitutive model in an isothermal 91 format. Subsequently, in Section 2.2, micromorphic regularization, thermomechanical coupling, and 92 finite-element implementation are briefly discussed. 93 2.1. Constitutive model We begin the model description by noting that functional fatigue in SMAs is typically attributed 94 to a number of mechanisms. Among them, generation of dislocation slip [4, 44], formation of 95 stabilized martensite [21, 45] and non-transforming austenite [45, 46] are the most likely involved 96 mechanisms. In the present model, a subdivision into the possible mechanisms and their mutual 97 interaction is not attempted, instead, they are unitedly represented by phenomenological evolution 98 equations, and are directly linked to the martensitic phase transformation through the accumulated 99 martensite volume fraction ηacc. In line with this notion, the inelastic mechanism responsible for 100 functional degradation is herein denoted as transformation-induced plasticity (TRIP). 101 The material state at each point is characterized by two quantities, namely the total strain 102 ε=1 2 ∇u+ (∇u)T , with uas the displacement vector, and the martensite volume fraction η. 103 The total strain is additively decomposed into 104 ε=εe+εt+εp, (1) where εedenotes the elastic contribution, εtdenotes the martensitic transformation contribution 105 andεpis the permanent strain associated with TRIP. At the same time, it is assumed that during 106 the martensitic transformation a fraction of martensite stabilizes and does not transform back 107 to austenite. Hence, the martensite volume fraction ηis split into the reversible part ηrevand 108 4 irreversible part ηir, viz., 109 η=ηrev+ηir, (2) and the following inequality constraints hold, 110 0≤ηir≤η≤1 = ⇒ 0≤ηrev≤1−ηir. (3) The material is in the fully austenitic state when η=ηrev= 0 and is in the fully martensitic state 111 when η= 1. Nevertheless, once the material starts transforming to martensite from a pristine 112 austenitic state, ηirbecomes immediately nonzero, as indicated by Eqs. (4)–(6) below, and thereby, 113 a fully austenitic state will not be recoverable. 114 It has been repeatedly observed in the experiments that the degradation of pseudoelasticity in 115 conventional polycrystalline NiTi are mostly pronounced during the first tens of cycles, gradually 116 diminishing and eventually reaching saturation as the material passes the so-called shakedown 117 stage, e.g., [22, 47, 48]. In view of this general consensus, we adopt the assumption that both the 118 irreversible volume fraction ηirand the permanent strain εpfollow exponential-type evolution laws. 119 Note that this assumption is not unique to the present model and has been exploited in various 120 SMA models that account for functional degradation, e.g., [37, 49–51]. With this assumption in 121 place, we first introduce the accumulated volume fraction ηaccas 122 ˙ηacc=|˙ηrev|=⇒ ηacc=Zt 0|˙ηrev|dτ, (4) where the overdot denotes the rate of change of the variable and tdenotes the time. The evolution 123 equation for the irreversible volume fraction ηiris then explicitly postulated as 124 ηir=hsat ir(1−exp(−Cpηacc)), (5) which results from the time-integration of the following rate equation (with ηacc t=0= 0 and 125 ηir t=0= 0, as for the initial conditions), 126 ˙ηir=hsat irCpexp(−Cpηacc) ˙ηacc. (6) Analogously, the evolution equation for the permanent strain εpis postulated as 127 ˙εp=ϵsat pCpexp(−Cpηacc) ˙ηaccNp. (7) In Eqs. (5)–(7), hsat irandϵsat prepresent the respective saturation values for irreversible volume 128 fraction and permanent strain, Cpis the degradation rate, and Npis the direction tensor which is 129 defined such that the rate of the permanent strain ˙εpis aligned with the martensitic transformation 130 strain εt, i.e., 131 Np=εt ∥εt∥,∥εt∥=p tr(εt)2. (8) Note that, in view of the definition of the accumulated volume fraction ηacc, the variables ηirand 132 εpevolve continuously during both the forward and backward transformations. 133 5 Martensitic transformation in SMAs usually exhibits negligible volumetric change [1]. The 134 transformation strain εtis therefore assumed to be deviatoric (i.e., tr εt= 0). Moreover, since 135 the stress-induced transformation renders the martensite variants to be oriented in the direction 136 of the applied stress, martensite is here considered to appear in a fully-oriented state so that the 137 transformation strain εtis defined as a function of the reversible volume fraction ηrevand the 138 transformation strain of fully-oriented martensite ¯εt, 139 εt=ηrev¯εt,¯εt∈¯P={¯εt:g(¯εt) = 0}. (9) The set ¯Pdefines the admissible limit transformation strain tensors characterized by the surface 140 g(¯εt) = 0 which is expressed in the following form [52], 141 g(¯εt) =h (−I2)3/2−bI3−cI3 4i1/3 −a. (10) In Eq. (10), I2andI3denote the principal invariants of the limit transformation strain tensor ¯εt142 while I4denotes a mixed invariant, defined as 143 I2=−1 2tr(¯εt)2, I 3= det ¯εt, I 4=m·¯εtm, (11) where mis the axis of the transverse isotropy. The parameters a,bandccharacterize the shape 144 and size of the surface g(¯εt) = 0 and are specified as 145 a=ϵTh3√ 3 4(1 + α3)i1/3 , b =√ 3 69α3β3−7α3+ 7β3−9 (1 +α3)(1 + β3), c =2√ 3 3α3−β3 (1 +α3)(1 + β3),(12) with ϵTas the maximum transformation strain in tension, αas the tension–compression asymmetry 146 ratio in the direction along the axis of transverse isotropy (i.e., parallel to m), and βas the tension– 147 compression asymmetry ratio in the direction perpendicular to the axis of transverse isotropy (i.e., 148 perpendicular to m). 149 It is noteworthy that the deviatoric nature of the transformation strain εtdictates, in accordance 150 with the definition of the direction tensor Np, see Eq. (8), that the permanent strain εpis also 151 deviatoric. Models within the present context often postulate that the permanent inelastic strain 152 evolves in the direction of stress deviator, e.g., [37, 49, 50]. In the present formulation, it can be 153 easily shown that the stress deviator is perpendicular to the surface g(¯εt) = 0, see [39], and thereby, 154 the transformation strain εtdepends on the direction of stress deviator. This, however, does not 155 imply that the transformation strain εt, and accordingly the permanent strain rate ˙εp, are colinear 156 with the stress deviator. 157 Another important aspect to highlight is that the accumulation of the irreversible volume frac- 158 tionηirand its impact on the reversible volume fraction ηrevcause the magnitude of the transfor- 159 mation strain εt, which serves as the actual transformation strain measure in the present model, 160 to decrease. However, the surface g(¯εt) = 0 and so the limit transformation strain ¯εtremain in- 161 tact throughout the cyclic transformation. This represents an underlying modeling assumption in 162 the present framework regarding the interaction between phase transformation and cyclic degra- 163 dation. It reflects the notion that the inherent characteristics of the transformation strain are not 164 6 affected during the cyclic degradation. Instead, it is the accumulation of TRIP and the decrease in 165 the amount of transformable (reversible) martensite that lead to the contraction of the maximum 166 attainable transformation strain. 167 We now elaborate on the Helmholtz free energy function and the dissipation potential, both customized to incorporate the degradation effects. Assuming an isothermal process, the Helmholtz free energy ϕis composed of the following contributions: the chemical energy ϕchem, the elastic strain energy ϕel, the austenite–martensite interaction energy ϕint, the energy of the diffuse interface ϕgrad, and the energy contribution ϕdegrelated to the pseudoelasticity degradation, i.e., ϕ(ε,¯εt,εp, ηrev,∇ηrev, ηir) = ϕchem(ηrev, ηir) +ϕel(ε,¯εt,εp, ηrev) +ϕint(ηrev) +ϕgrad(∇ηrev) +ϕdeg(ηrev, ηir).(13) Among the contributions to the Helmholtz free energy ϕ, only ϕdegis specific to the present model. 168 The remaining contributions are rather standard and adhere to the non-cyclic model of pseudoe- 169 lasticity [39–41] and are formulated as 170 ϕchem(ηrev, ηir) = (1 −η)ϕa 0+ηϕm 0=ϕa 0+ ∆ϕ0η, (14) 171 ϕel(ε,¯εt,εp, ηrev) =µtr(εe dev)2+1 2κ(trεe)2,εe=ε−ηrev¯εt−εp, (15) 172 ϕint(ηrev) =1 2Hint(ηrev)2, (16) 173 ϕgrad(∇ηrev) =1 2G∇ηrev· ∇ηrev. (17) Here, ∆ ϕ0=ϕm 0−ϕa 0is the phase transformation chemical energy, µis the elastic shear modulus 174 and is calculated via applying the Reuss averaging scheme based on the total volume fraction ηto 175 the shear moduli of austenite µaand martensite µm(i.e., 1 /µ= (1−η)/µa+η/µm),κis the elastic 176 bulk modulus (assumed constant), Hintis the parameter that characterizes the material response 177 (softening- or hardening-type) within the transformation regime, and G >0 is the gradient energy 178 coefficient. Note that the parameter Hintcan be adapted such that it reflects a loading-dependent 179 material response (typically, a softening-type response in tension and hardening-type response in 180 compression), e.g., [42]. However, for simplicity, Hintis here considered as a constant parameter. 181 Given that the simulations in this study involve predominantly tensile loading, see Section 3, this 182 simplification does not pose a serious limitation. Note also that the interaction energy ϕintis 183 a quadratic function of the volume fraction ηrev, resulting in a tri-linear intrinsic stress–strain 184 response, as illustrated in Fig. 2. This choice is also made for simplicity and can be readily adapted 185 to more complex functions to achieve a more realistic response [43]. 186 On the other hand, the degradation contribution ϕdegtakes the following form 187 ϕdeg(ηrev, ηir) =Adegηirηrev+1 2Hdegηir(ηrev)2, (18) where AdegandHdegrepresent the degradation parameters. The contribution ϕdegis specifically 188 tailored to address two primary effects of pseudoelasticity degradation: it accounts for the reduction 189 7 of the transformation-onset stress (described by the term Adegηirηrev) and the conversion of the 190 mechanical response towards a hardening-type response (described by the term1 2Hdegηir(ηrev)2). 191 In line with the evolution of ηir, Eq. (5), both effects progress exponentially. Note that the approach 192 of incorporating the cyclic degradation effects into the free energy function has been also used in 193 other SMA models in the literature, e.g., [53, 54]. 194 Finally, a rate-independent dissipation potential is adopted in the following form 195 D( ˙ηrev, ηacc) =fc(ηacc)|˙ηrev|, (19) where fc(ηacc), which is called the critical thermodynamic driving force, controls the width of the 196 hysteresis loop in the stress–strain response. To capture the decrease in the hysteresis loop area 197 (i.e., the dissipated energy) during the cyclic transformation, the parameter fcis defined in relation 198 to the accumulated volume fraction ηacc. Similar to the permanent strain εpand the irreversible 199 volume fraction ηir, Eqs. (5)–(7), fcevolves exponentially as follows 200 fc(ηacc) =ffin c+ (fini c−ffin c) exp(−Cfηacc), (20) where fini candffin crepresent, respectively, the initial and final values of fc, and Cfdenotes the 201 corresponding evolution rate. 202 To formulate the incremental energy minimization problem, we derive the time-discrete version 203 of the constitutive equations by employing the backward Euler scheme. Having known the variables 204 related to the previous time step tn, the variables related to the current time step tn+1=tn+∆tare 205 sought. We begin by approximating the incremental evolution equation for the irreversible volume 206 fraction ηirand the permanent strain εp, 207 ∆t˙ηir≈∆ηir=hsat irCpexp(−Cpηacc)∆ηacc,∆t˙εp≈∆εp=ϵsat pCpexp(−Cpηacc)∆ηaccNp,(21) where 208 ηacc=Ztn+1 0∆ηaccdτ,∆ηacc=|∆ηrev|,∆ηrev=ηrev−ηrev n, (22) with ηrev nas the value of the reversible volume fraction from the previous time step tn. At the same 209 time, the incremental form of the rate-independent dissipation potential is obtained as 210 ∆D(∆ηrev, ηacc) =fc(ηacc)|∆ηrev|. (23) The solution of the problem is determined via the incremental energy minimization principle 211 [39, 41, 55]. A global incremental potential Π is defined by summing up the increment of the 212 total Helmholtz free energy ∆Φ (where Φ =R BϕdV), the global dissipation potential ∆ D(where 213 ∆D=R B∆DdV) and the potential of the external loads ∆Ω, and is subsequently minimized with 214 respect to the unknowns u,¯εtandηrev, i.e., 215 Π = ∆Φ + ∆Ω + ∆ D → min u,¯εt,ηrev(24) which is subject to the inequality constraints on the reversible volume fraction ηrev, Eq. (3), and 216 8 (a) (b) referencecycle 1cycle 2cycle 3 0 3 6 90175350525700 Axial strain, ε[%]Axial stress, σ[MPa ] cycle 5cycle 10cycle 20cycle 50 0 3 690175350525700 Axial strain, ε[%]Axial stress, σ[MPa ]Figure 2: The intrinsic stress–strain response of the model under full-transformation cycles of uniaxial-tension: (a) the first three cycles, and (b) the first 50 cycles. The dashed curve in panel (a), denoted as ‘reference’, represents the pseudoelastic intrinsic response with no degradation effects. The model parameters adopted to produce the intrinsic responses are the same as those in the main simulations, see Section 3. to the constraint related to the limit transformation strain surface, Eq. (9). At the same time, 217 ηirandεp, which contribute directly to the minimization problem, are explicitly evaluated from 218 Eq. (21). To provide a clearer idea of the structure of the minimization problem and the underlying 219 constitutive behavior of the model, a simplified 1D version of the model is elaborated in Appendix 220 A. 221 Fig. 2 showcases the intrinsic stress–strain response predicted by the model under cyclic tensile 222 loading. Two cases are highlighted: the pseudoelasticity degradation effects observed within the 223 first three cycles, relevant to the problem of subloop deformation investigated in this study, and 224 the degradation effects observed within 50 cycles, which provides a more holistic view of the model 225 behavior. Note that the material parameters adopted to generate the intrinsic response in Fig. 2 226 are the same as those adopted in the main simulations in Section 3. 227 2.2. Further extensions and finite-element implementation The model presented in Section 2.1 is now enriched with micromorphic regularization and is 228 made thermomechanically coupled. Both extensions have been thoroughly discussed in our previous 229 works [40, 41]. Hence, we only briefly discuss them here. 230 The purpose of adopting the micromorphic regularization is to facilitate the finite-element imple- 231 mentation of the gradient-enhanced model by restructuring the minimization problem in a way that 232 the constitutive complexities are transferred to the local level (for instance, at the Gauss points) 233 where they can be treated in a more efficient way. To do so, a new degree of freedom ˘ ηis introduced 234 and is coupled with the volume fraction ηrevthrough the following penalization term ϕpenwhich is 235 added into the Helmholtz free energy function ϕ, see Eq. (13), 236 ϕpen(ηrev,˘η) =1 2χ(ηrev−˘η)2, (25) 9 with χas the penalty parameter. The gradient energy ϕgrad, see Eq. (17), is then redefined in terms 237 of the gradient of the new variable ˘ η, i.e., 238 ϕgrad(∇˘η) =1 2G∇˘η· ∇˘η. (26) Following this modification, the volume fraction ηrevcan be considered as a local quantity and the 239 respective evolution equation can be solved (together with that of ¯εt) at the local level. For further 240 details regarding the micromorphic regularization, interested readers are referred to [56, 57]. 241 To arrive at a thermomechanically-coupled model, two most important couplings are taken 242 into consideration. First, the chemical energy ϕchem, Eq. (14), is extended to reflect the effect of 243 temperature on the mechanical response (the Clausius–Clapeyron relation), i.e., 244 ϕ0(ηrev, ηir, T) =ϕa 0(T) + ∆ ϕ0(T)η,∆ϕ0(T) = ∆ s∗(T−Tt), (27) where ∆ s∗represents the transformation entropy change, Tis the temperature, and Ttis the 245 transformation equilibrium temperature. Next, the internal heat source ˙Ris defined to encompass 246 the latent heat of transformation and the heat release by mechanical dissipation, viz., 247 ˙R= ∆s∗T˙ηrev+fc(ηacc)|˙ηrev|. (28) Eq. (28) is then introduced into the (isotropic) heat conduction equation 248 ϱ0c˙T+∇ ·Q=˙R,Q=−K∇T, (29) where Qis the heat flux, ϱ0cis the specific heat, and the scalar Kis the heat conduction coefficient. 249 It follows from Eq. (28) that the internal heat generation is influenced during the cyclic phase 250 transformation. This influence is manifested by both the latent heat of transformation and the 251 mechanical dissipation and operates through the reversible volume fraction ηrevand the hysteresis 252 parameter fc, cf. Eqs. (3) and (20). 253 The full thermomechanically-coupled model comprises three global unknown fields: the displace- 254 mentu, the micromorphic variable ˘ ηand the temperature T; and two local unknown variables: the 255 reversible volume fraction ηrevand the limit transformation strain ¯εt. Recall that the irreversible 256 volume fraction ηirand the permanent strain εpare explicitly integrated by using Eq. (21). The 257 finite-element discretization of the displacement field uis performed by using 20-noded quadratic 258 hexahedral (Serendipity) elements with reduced Gauss integration rule (2 ×2×2). On the other 259 hand, 8-noded linear hexahedral elements with standard Gauss integration rule (2 ×2×2) are 260 used for ˘ ηandT. For the 2D axisymmetric wire problem discussed in Section 3.2, the respective 261 discretizations have been done by 8-noded quadratic elements and 4-noded linear elements. The 262 resulting global–local problem is structured as a nested iterative-subiterative scheme and is solved 263 at both the global and local levels by using the Newton method. Notably, a fully-coupled monolithic 264 scheme is adopted so that the problem is solved simultaneously with respect to all unknowns. 265 It is worth noting that the local minimization problem of ηrevis non-smooth, in view of the 266 rate-independent dissipation potential, see Eqs. (19) and (23). To address this issue, the aug- 267 10 mented Lagrangian method is utilized, which handles adeptly both the non-smoothness of the 268 rate-independent dissipation potential and the inequality constraints on the reversible volume frac- 269 tionηrev, i.e., 0 ≤ηrev≤1−ηir, see Eq. (3). The local problem has an additional constraint to be 270 satisfied, namely the equality constraint of the limit transformation strain surface, g(¯εt) = 0, see 271 Eq. (9). The latter is addressed by using a standard Lagrange multiplier method. For brevity, the 272 related technical details are not discussed here, see [39]. 273 The model is transformed into a finite-element code using the automatic differentiation tool 274 AceGen [58, 59], thanks to which the residual vector and the tangent matrix are derived automat- 275 ically, and thereby, the quadratic convergence of the Newton method is ensured. The simulations 276 are carried out in the finite-element environment AceFEM. 277 3. Simulations This section is devoted to the analysis of the simulation results. Section 3.1 begins with a 278 presentation of the simulation setup and calibration of the model parameters, and concludes with 279 a brief discussion on the results for NiTi wire under full loading–unloading cycles. Our main 280 modeling study concerns a NiTi specimen subjected to uniaxial tension with subloop loading paths. 281 Two scenarios are explored. First, in line with the experimental study of Tobushi et al. [24], the 282 subloop deformation behavior of a NiTi wire is analyzed, see Section 3.2. The loading program in 283 this scenario encompasses three nested subloops with increasing strain amplitudes, as depicted in 284 Fig. 3(a). As shown later, this setup enables us to reproduce neatly the hierarchical return-point 285 memory. Section 3.2 concludes with a supplementary analysis of the TRIP evolution under a large 286 number of subloops. Next, in Section 3.3, we extend our analysis to a NiTi strip, where we elucidate 287 how the subloop behavior is influenced by the complexity of the pattern of propagating instabilities. 288 This scenario is then examined under two additional loading programs, see Fig. 3(b,c). 289 3.1. Preliminaries In all simulations, the loading is exerted in a displacement-control mode at a (constant) low 290 strain rate of 1 .67×10−4s−1. The NiTi wire has a diameter of 0.75 mm and a total length of 291 L0= 20 mm. To facilitate this analysis, the wire is justifiably reduced to a 2D axisymmetric 292 geometry. The corresponding 2D problem is then discretized by a uniform finite-element mesh 293 consisting of equiaxed elements with an edge size of 0.01 mm. This resulted in 76 000 elements 294 and approximately 620 000 degrees of freedom. Meanwhile, the NiTi strip is treated as a full 3D 295 problem. The strip has a thickness of 0.4 mm, a width of 10 mm and a total length of L0= 100 mm. 296 It is discretized by a uniform mesh consisting of elements with an in-plane edge size of 0.2 mm and 297 a through-thickness size of 0.4 mm (i.e., only one element is used through the thickness). This 298 mesh led to 25 000 elements and nearly 640 000 degrees of freedom. In both problems, the following 299 boundary conditions are imposed. The displacements at the bottom edge of the specimen are fully 300 constrained. At the top edge, the axial displacement δis prescribed and the lateral displacements 301 are constrained. At the same time, the temperature at both top and bottom edges is set equal to the 302 ambient temperature, i.e., T=T0= 353 K, which is the actual ambient temperature maintained 303 during the experiment [24]. Finally, the heat convection effect is neglected. 304 11 (a)Loading program 1 (b)Loading program 2 (c)Loading program 3 subloop ε_ max[%]ε_ min[%] 1 5.0 3.5 2 5.75 2.75 3 6.5 2.0 0 500 1000 1500 2000036912 Time [s]Average elongation, ε_=δ/L0[%]subloop ε_ max[%]ε_ min[%] 1 6.5 2.0 2 5.75 2.75 3 5.0 3.5 0 500 1000 1500 2000036912 Time [s]Average elongation, ε_=δ/L0[%]subloop ε_ max[%]ε_ min[%] 1 3.5 1.5 2 5.5 3.5 3 7.5 5.5 0 400 800 1200 1600036912 Time [s]Average elongation, ε_=δ/L0[%]Figure 3: The loading programs used in the simulations. All loading programs represent displacement-control uniaxial tension (with a low strain rate of 1 .67×10−4s−1) and incorporate three subloops. Loading program 1 consists of three nested subloops with increasing strain amplitudes. Loading program 2 employs subloops in a reverse order compared to loading program 1. Loading program 3 consists of three equally-spaced subloops with a constant strain amplitude. The model parameters adopted in the simulations are summarized in Tab. 1. Except for the 305 gradient energy parameter G, all the model parameters are identical in the wire and strip problems. 306 The parameter Gsets the length-scale associated with the phase transformation front and can be 307 linked to the geometry and micromechanical characteristics [60, 61]. Thus, Gtakes different values 308 in each problem. To calibrate G, first, an assumption ought to be made regarding the theoretical 309 thickness of the macroscopic interface, λ. Subsequently, Gis determined through the analytical 310 relation G=−Hintλ2/π2, which is derived from the solution of the 1D small-strain model of 311 pseudoelasticity [40]. The identification procedure for the remaining model parameters which are 312 unrelated to TRIP has been thoroughly discussed in our recent study [43], see Section 2.3 and 313 Appendix E therein, and is not repeated here. 314 Identification of some TRIP-related parameters is guided by the indications obtained from the 315 structural stress–strain response from the experiment, see Fig. 1(a). These include the significant 316 decrease in the level of the upper stress plateau during the hierarchical subloop deformation and 317 the value of the residual strain at the end of the experiment. Accordingly, the parameters Adeg= 318 −45 MPa, ϵsat r= 0.4ϵT= 2.4% (recall that ϵTdenotes the maximum transformation strain, Eq. (12)) 319 andCp= 0.05 have been calibrated to produce similar effects. We, however, acknowledge that there 320 exists a degree of uncertainty in the identification of the remaining parameters, for which we lack 321 definitive experimental evidence. With this in mind, the parameters Hdeg= 40 MPa, ffin c= 2 MPa, 322 Cf=Cp= 0.05 and hsat ir= 0.4 are selected such that the changes in the stress–strain response 323 under a large number of loading cycles (in particular, as it concerns the transition to a hardening- 324 type response, decrease in the hysteresis loop area and decrease in the extent of the transformation 325 strain) align with the trends observed in the experiments, e.g., [47, 62–64], see also the discussion 326 below. The intrinsic response of the model resulting from the adopted parameters is illustrated in 327 Fig. 2. 328 12 Table 1: Model parameters adopted in the simulations. Category Parameter Value Elasticity κ Bulk modulus 130 GPa µa Shear modulus for austenite 21 GPa µm Shear modulus for martensite 9 GPa Martensitic phase transformation ∆ s∗Specific entropy difference 0.24 MPa/K Tt Transformation equilibrium temperature 222 K fini c Hysteresis loop parameter (initial) 10 MPa Hint Austenite–martensite interaction parameter -10.5 MPa ϵT Maximum tensile transformation strain 6% α Tension–compression asymmetry ratio 1.4 β Transverse isotropy parameter 1.0 Macroscopic transformation front G Gradient energy parameter (wire problem) 0.04 MPa mm2 G Gradient energy parameter (strip problem) 0.4 MPa mm2 χ Micromorphic regularization parameter 100 MPa Heat transfer ϱ0c Specific heat 2.86 MJ/(m3K) K Heat conductivity 18 W/(m K) TRIP Adeg Pseudoelasticity degradation parameter -45 MPa Hdeg Pseudoelasticity degradation parameter 40 MPa ϵsat p Saturated permanent strain 0.4 ϵT= 2.4% hsat ir Saturated irreversible volume fraction 0.4 ffin c Hysteresis loop parameter (final) 2 MPa Cp Degradation rate 0.05 Cf Hysteresis loop degradation rate 0.05 13 (a) cycle 1cycle 5cycle 10cycle 20cycle 50 0 3 6 90175350525700 Average elongation, ε_=δ/L0[%]Average axial stress, σ_=P/A0[MPa ](b) 0 2 4 6 80150300450600 Axial strain [%]Axial stress [MPa ] cycle 1 cycle 30(c) 024680200400600800 Axial strain [%]Axial stress [MPa ] cycle 1 cycle 50 (d) 02468100150300450600 Axial strain [%]Axial stress [MPa ] cycle 1 cycle 20(e) 02468100300600900 Axial strain [%]Axial stress [MPa ]cycle 1 cycle 100Figure 4: (a) NiTi wire subjected to 50 loading–unloading cycles of uniaxial tension: structural stress–elongation (¯σ–¯ε) response. The average axial stress ¯ σand the average elongation ¯ εare calculated, respectively, as the reaction force Pdivided by the initial cross-section area A0, and the axial displacement δdivided by the initial length L0. The intrinsic responses associated with these structural responses are illustrated in Fig. 2. (b)–(e) Typical cyclic responses of NiTi specimens observed in the experiments, taken from (b) Wang et al. [47], (c) Morin et al. [62], (d) Kan et al. [63], and (e) ˇSittner et al. [64]. It is worth noting that in all the simulations, as a way to trigger the phase transformation 329 instability, a geometric imperfection in the form of a slight indent is applied to the specimen. The 330 indent is located at a distance equal to the diameter/width of the wire/strip from its lower end. 331 Before entering into the main analysis of subloop deformation, a simulation is performed for 332 the NiTi wire subjected to 50 loading–unloading cycles of uniaxial tension. Fig. 4 illustrates the 333 structural response of the wire. Here, as well as in the figures in the following subsections, the 334 structural response is represented in terms of the average axial stress ¯ σ=P/A 0versus average 335 elongation (engineering strain) ¯ ε=δ/L 0, where Pdenotes the reaction force and A0denotes the 336 initial cross-section area. Recall that δandL0are the axial displacement and the initial length, 337 respectively. As it is evident in the structural response, the wire undergoes a complete phase 338 transformation within each cycle. Initially, the wire exhibits a localized phase transformation, 339 characterized by a stress drop at the transformation onset and a subsequent stress plateau. The 340 localized transformation persists for about 15 cycles. Thereafter, the transformation proceeds in a 341 more homogeneous manner, and the structural response displays a mild hardening. As the number 342 of cycles increases, the slope of the hardening branch also increases. The cyclic behavior captured 343 in the simulation is in a qualitative agreement with the typical cyclic behavior of NiTi specimens 344 observed in experiments [47, 62–64], which underscores the reliability of the simulation results. 345 14 3.2. NiTi wire subjected to subloop deformation The results pertaining to the subloop behavior of the NiTi wire are presented in Figs. 5 and 6. 346 The phase transformation evolution in Fig. 5(a) and TRIP evolution in Fig. 5(b) are displayed via, 347 respectively, the distribution of the reversible volume fraction ηrevand irreversible volume fraction 348 ηir. Note that, for a more natural visualization, the results of the 2D axisymmetric wire are post- 349 processed and presented in a 3D configuration. As anticipated, the transformation initiates at the 350 position of the geometric imperfection. Throughout the entire loading stage of the global cycle 351 (hereinafter, to avoid confusion with the subloops, we use the term ‘global’), the transformation 352 maintains a single propagating front. Interestingly, while the front appears to be a flat (and visibly 353 diffuse) interface in the 3D-wire configuration, e.g., [65], it takes on a spherical-shaped appearance 354 (or ‘cone-shaped’ as described in [60, 66]), as can be conceived from the corresponding pattern in the 355 axisymmetric planes (not shown here). During the global unloading, the backward transformation 356 commences from the wire’s central part. As shown in Fig. 5(b) and discussed below, the highest 357 amount of irreversible volume fraction ηir, thus the highest TRIP, is accumulated within the central 358 part, making it a favorable site for the nucleation of the austenitic band. At the same time, due to 359 a slight asymmetry in the distribution of ηirwith respect to the wire’s midpoint, the two evolved 360 fronts do not propagate concurrently. More specifically, the top front reaches the boundary first 361 and annihilates, which manifests as an abrupt stress rise in the structural stress-elongation (¯ σ–¯ε) 362 response, occurring at an average elongation of about ¯ ε= 4% (see Fig. 6). Subsequently, the bottom 363 front follows suit. 364 Within each subloop path, the front retreats downward during unloading and advances upward 365 during reloading. This cyclic movement prompts the material points inside the front’s sweeping 366 zone to undergo backward-then-forward transformation, and thereby, gives rise to the accumulation 367 of TRIP within the sweeping zone, while the material points beyond it remain unaffected. Note that 368 the loading program adheres to a fixed nominal mean strain (set at ¯ ε= 4.5%, which corresponds 369 to the front’s proximity to the wire’s midpoint) but an increasing strain amplitude, see Fig. 3(a). 370 Thus, the sweeping zone expands successively from subloop 1 to subloop 3, and at the same time, 371 the sweeping zone of each subloop encompasses that of the previous one. This therefore results in 372 the highest concentration of TRIP within the central part of the wire and its step-wise decreasing 373 trend as it moves away from it, as can be clearly seen in Fig. 5(b). 374 The hierarchical return-point memory, which is an outcome of the cyclic traversal of the front 375 across the boundaries of the sweeping zones, is correctly reproduced in the structural stress– 376 elongation response in Fig. 6. The reproduced feature is in a reasonable agreement with the 377 experimental result of Tobushi et al. [24], see Fig. 1 and the accompanying discussion. In view 378 of the exponential nature of the pseudoelasticity degradation effects, the reduction in the level of 379 the upper stress plateau is at the highest within the first level of hierarchy (of about ∆¯ σ= 21 380 MPa) and diminishes to its lowest within the last level of hierarchy (of about ∆¯ σ= 17 MPa). It is 381 worth noting that in this scenario, where the strain rate corresponds to nearly isothermal conditions 382 (i.e., the temperature variation lies within the range of −2 K to 2 K), the stress, upon reaching 383 the return-point, appears to catch up closely with the corresponding stress plateau before applying 384 the subloop. As shown in Section 3.3 and also observed in the experiments [25, 26], such a close 385 catching up does not occur when thermal effects are at play. 386 15 subloop 1 subloop 2 subloop 3 global unloading global loading subloop 1 subloop 2 subloop 3 global unloading global loadingFigure 5: NiTi wire subjected to loading program 1: (a) snapshots of reversible volume fraction ηrevillustrating the phase transformation evolution, and (b) snapshots of irreversible volume fraction ηirillustrating the TRIP evolution. For a natural visualization, the axisymmetric wire is presented in full 3D configuration. 16 (a) (b) (c) subloop 1 0 3 6 90175350525700 Average elongation, ε_=δ/L0[%]Average axial stress, σ_=P/A0[MPa ] subloop 2 0 3 690175350525700 Average elongation, ε_=δ/L0[%]Average axial stress, σ_=P/A0[MPa ] subloop 3 0 3 690175350525700 Average elongation, ε_=δ/L0[%]Average axial stress, σ_=P/A0[MPa ]Figure 6: NiTi wire subjected to loading program 1: structural stress–elongation (¯ σ–¯ε) response. A notable observation from the experimental curve in Fig. 1(a) is the absence of the return-point 387 memory during the unloading stages of the subloops. Instead, the lower stress plateau seems to 388 shift slightly downward from one subloop to the next. Anyway, no attempt was made to adjust the 389 material parameters to replicate the observed behavior, which is, however, present in the results of 390 NiTi strip reported in Section 3.3. 391 In concluding the discussion in this section, we present the results of a supplementary analysis 392 on the NiTi wire subjected to 12 subloops. The aim of this analysis is to illustrate the evolution 393 of TRIP and subloop deformation behavior over a large number of subloops. The results, as 394 depicted in Fig. 7, follow an expected trend. However, two specific observations deserve further 395 comment. Firstly, the stress plateau in a number of subloops exhibits irregularities, specifically a 396 second stress drop appears ahead of the return-point memory. This is because in these subloops 397 the transformation during subloop reloading does not proceed by the propagation of the existing 398 front. Instead, a second front emerges at the opposite end of the sweeping zone and eventually 399 merges with the original front, thereby, leading to the observed effects. Secondly, as shown in 400 Fig. 7(b), the sweeping zone continuously expands from one subloop to the next. This is explained 401 by the accumulation of TRIP within the sweeping zone, which reduces the amount of transformable 402 martensite. Consequently, since the applied strain amplitude of the subloops is held fixed, the front 403 gradually moves towards the untransformed segments of the wire to compensate for the reduced 404 transformation. 405 3.3. NiTi strip subjected to subloop deformation We begin this section by analyzing the NiTi strip under loading program 1. The primary aim is 406 to examine the subloop behavior in a notably more involved scenario than the NiTi wire discussed 407 earlier, arising mainly from a more complex transformation pattern and heightened thermal effects. 408 The simulation results are presented in Figs. 8 and 9. A quick look at Fig. 9 immediately indicates 409 that the return-point memory is only observable in the trajectories that lead to the global stress 410 plateau, while the hierarchical return-point memory is lost. This is undoubtedly an outcome of 411 17 (a) (b) 0 3 6 90175350525700 Average elongation, ε_=δ/L0[%]Average axial stress, σ_=P/A0[MPa ] Timeε_1 2 12 Figure 7: NiTi wire subjected to uniaxial tension with 12 subloop paths: (a) structural stress–elongation (¯ σ–¯ε) response, and (b) snapshots of irreversible volume fraction ηirillustrating the TRIP evolution at the end of each subloop. The inset in panel (a) represents the loading program used for this simulation. The applied subloops have the same mean strain and the same strain amplitude, i.e., ¯ εmax= 5.75% and ¯ εmin= 2.75%. the nontrivial pattern of phase transformation and resulting TRIP distribution within the strip. 412 In contrast to the NiTi wire, where a single phase transformation front remained active during 413 all subloops, the strip features multiple transformation fronts, each presenting a less predictable 414 pattern of activation. Below, we provide a more detailed account of the unfolding events. 415 The phase transformation initiates with the nucleation of a single martensite band at the location 416 of the geometric imperfection. The band is oriented at approximately 54◦with respect to the 417 longitudinal axis, in agreement with the experimental observations [14] and theoretical analysis 418 [67], and somewhat changes as loading progresses. At an average elongation of about ¯ ε= 3%, 419 another martensite band emerges at the opposite end, and henceforth, the two transformation 420 fronts propagate towards each other. This non-synchronous double nucleation has been commonly 421 observed in the experiment of NiTi specimens at relatively low strain rates, e.g., [14, 68, 69]. 422 Within subloop 1, the two fronts exhibit a short back-and-forth movement, manifesting a clear 423 return-point memory in the structural stress–elongation response in Fig. 9(a). During subloop 424 2, not only the hitherto active fronts but also the fronts near the boundaries become engaged in 425 the transformation evolution. As a consequence, TRIP is induced via all fronts. This behavior is 426 reflected in the structural response which takes on an irregular appearance characterized by few 427 sudden stress changes, thus spoiling the return-point memory in the inner part (Fig. 9(b)). A 428 similar process recurs within subloop 3, albeit with a more complex phase transformation evolution 429 during the reloading stage and also more distinct stress events in the structural response. During the 430 global unloading, the backward transformation proceeds predominantly in a criss-cross mode, which 431 persists until an average elongation of about ¯ ε= 3%. Subsequently, the fronts reconfigure into sharp 432 inclined interfaces that move towards each other until the complete annihilation of the (reversible) 433 18 subloop 1 subloop 2 subloop 3 global unloading global loading subloop 1 subloop 2 subloop 3 global unloading global loadingFigure 8: NiTi strip subjected to loading program 1: (a) snapshots of reversible volume fraction ηrevillustrating the phase transformation evolution, and (b) snapshots of irreversible volume fraction ηirillustrating the TRIP evolution. 19 (a) (b) (c) subloop 1 0 3 6 90175350525700 Average elongation, ε_=δ/L0[%]Average axial stress, σ_=P/A0[MPa ] subloop 2 0 3 6 90175350525700 Average elongation, ε_=δ/L0[%]Average axial stress, σ_=P/A0[MPa ] subloop 3 0 3 690175350525700 Average elongation, ε_=δ/L0[%]Average axial stress, σ_=P/A0[MPa ]Figure 9: NiTi strip subjected to loading program 1: structural stress–elongation (¯ σ–¯ε) response. martensite domain. The reconfiguration of the fronts between criss-cross mode and sharp interfaces 434 has been also observed in experimental and previous modeling studies, e.g., [36, 70, 71]. It is worth 435 remarking that the distribution of the irreversible volume fraction ηirwithin the entire strip at the 436 end of the global unloading remains consistent with that at the end of the global loading, while its 437 magnitude increases uniformly. 438 Upon inspecting the return-point memory in Fig. 9, a slight difference can be noticed concerning 439 the level of the global stress plateau before and after a subloop path. This difference stems from the 440 thermal effects. Specifically, compared to the NiTi wire, a more pronounced temperature variation 441 is produced across the specimen during the forward transformation (for instance, of about 10 K 442 immediately before subloop 1), resulting in a more visible thermal hardening that sustains a higher 443 stress for the propagation of the front. Within the subloop path, the transformation latent heat is 444 initially absorbed during the backward transformation (self-cooling) and is subsequently released 445 when the forward transformation resumes (self-heating). Accordingly, as the front reaches the pris- 446 tine material, the temperature variation across the specimen is reduced compared to the state before 447 the subloop (for instance, of about 5 K immediately after subloop 1). Thereby, thermal hardening 448 diminishes, necessitating a lower stress for interface propagation. Note also that as a result of the 449 cyclic transformation of the material points, and thus the accumulation of irreversible martensite, 450 a smaller martensite volume fraction is transformed during the subloop reloading compared to the 451 state before the subloop, and this contributes to the reduction of the latent heat generation [20], 452 see Eq. (28). 453 We now proceed with the analysis of the NiTi strip under two additional loading programs, one 454 consisting of nested subloops with decreasing strain amplitudes, i.e., subloops are applied in a reverse 455 order compared to loading program 1, and the other consisting of three equally-spaced distinct 456 subloops with a constant strain amplitude, see Fig. 3(b,c). The corresponding results are shown in 457 Figs. 10, 11 and 12. The comparison of the snapshots of the reversible volume fraction ηrevin the two 458 additional cases to those of loading program 1 reveals noticeable morphological differences, which are 459 20 subloop 1 subloop 2 subloop 3 global unloading global loading subloop 1 subloop 2 subloop 3 global unloading global loadingFigure 10: NiTi strip subjected to loading program 2: (a) snapshots of reversible volume fraction ηrevillustrating the phase transformation evolution, and (b) snapshots of irreversible volume fraction ηirillustrating the TRIP evolution. 21 subloop 1 subloop 2 subloop 3 global unloading global loading subloop 1 subloop 2 subloop 3 global unloading global loadingFigure 11: NiTi strip subjected to loading program 3: (a) snapshots of reversible volume fraction ηrevillustrating the phase transformation evolution, and (b) snapshots of irreversible volume fraction ηirillustrating the TRIP evolution. 22 (a) (b) 0 3 6 90175350525700 Average elongation, ε_=δ/L0[%]Average axial stress, σ_=P/A0[MPa ] 0 3 6 90175350525700 Average elongation, ε_=δ/L0[%]Average axial stress, σ_=P/A0[MPa ]Figure 12: Structural stress–elongation (¯ σ–¯ε) response of NiTi strip subjected to (a) loading program 2 and (b) loading program 3. beyond the differences arising solely from the loading-dependent transformation evolution pathways. 460 The differences mainly concern the varying number of martensite domains formed during the global 461 loading stage and the activation pattern of the fronts within the subloops. More specifically, unlike 462 loading program 1, loading programs 2 and 3 exhibit only two martensite domains during the global 463 loading. In loading program 2, all four fronts remain consistently active within all subloops, resulting 464 in a clear demonstration of the hierarchical return-point memory in the stress–elongation response, 465 as shown in Fig. 12(a). In loading program 3, however, while the involvement of the fronts near the 466 boundaries is eye-catching within subloop 1, overall, the interior fronts are prominently active. In 467 this case, the front sweeping zones in the subloops do not interact with each other (as can be also 468 recognized from the snapshots of ηirin Fig. 11), and the resulting subloops are independent, see 469 Fig. 12(b). During the global unloading, all cases show a similar transformation evolution pattern 470 characterized by two active fronts retracting in a criss-cross manner.. 471 We conclude this discussion by addressing TRIP accumulation within the strip in relation to 472 the loading program. Similar to the martensitic transformation, TRIP exhibits characteristics that 473 are specific to the applied loading program. Given that loading programs 1 and 2 have a reverse 474 arrangement of the subloops but are otherwise identical, one would intuitively expect that the 475 resulting TRIP accumulations, in terms of both the pattern and the intensity, would be the same 476 after applying all the three subloops. A comparison of the snapshots of the irreversible volume 477 fraction ηir(Figs. 8 and 10) indeed confirms that TRIP hotspots in these two cases are located in 478 nearly the same regions, with two hotspots near the boundaries and two within the interior of the 479 strip, corresponding to the regions with the highest activity of the fronts. Yet, minor discrepancies 480 can be observed, particularly concerning the intensity of TRIP within the hotspot regions. On the 481 other hand, loading program 3 demonstrates a rather distinct TRIP accumulation characterized by 482 several regions with mild intensity within the interior and localized hotspots near the boundaries. 483 As previously noted, this particular TRIP distribution results from the lack of interaction among 484 the fronts sweeping zones of the independent subloops. As a summary of this discussion, Fig 13 485 23 Loading program 1 0 25 50 75 10000.050.10.15 Position [mm]ηir[-]Loading program 1Loading program 2Loading program 3Figure 13: Distribution of the irreversible volume fraction ηiralong the entire length of the strip (taken in the reference configuration) at the end of the global loading stage. The graphs correspond to the midsection of the strip, as indicated by the white dashed curve overlaid on the snapshot. compares the distribution of ηiralong the entire length of the strip for various loading programs. 486 4. Conclusions The phenomenon of return-point memory that appears during the subloop deformation of pseu- 487 doelastic SMA is an outcome of the interaction between the structural instabilities of phase trans- 488 formation and the degradation of functional properties. It seems that this crucial aspect has been 489 generally overlooked in existing modeling approaches. The goal of our study is to demonstrate this 490 structure–material interaction by modeling the phenomenon of return-point memory. To achieve 491 this, we have developed a gradient-enhanced model of pseudoelasticity. The developed model rep- 492 resents an advancement over previous versions [39–41], extending the constitutive description to 493 incorporate pseudoelasticity degradation. The capabilities of the model in reproducing the essen- 494 tial aspects of pseudoelasticity degradation have been shown for NiTi under cyclic uniaxial tension. 495 We examine an illustrative example of a NiTi wire subjected to nearly isothermal uniaxial tension 496 with nested subloops. The obtained results clearly correlate with the experimental observations of 497 Tobushi et al. [24], especially regarding the hierarchical return-point memory. The accumulation 498 of TRIP and its distribution during subloop deformation underline the intertwined evolution of 499 inhomogeneous phase transformation and cyclic degradation. 500 The study is then extended to a more involved scenario of a NiTi strip, where a detailed analysis 501 is performed by examining three different loading programs. The impact of the loading program on 502 the evolution of phase transformation and TRIP has been highlighted through the activation pat- 503 tern of phase transformation fronts within the subloops, and its implications on the phenomenon 504 of return-point memory have been pointed out. In addition, the results hint at the visible con- 505 tribution of the thermomechanical coupling effects within the subloops, stemming from the self- 506 cooling/heating process of the transforming material. 507 Acknowledgement. The authors wish to thank Prof. Stanis law Stupkiewicz for helpful discussions 508 on the model formulation and reviewing the manuscript. This work has been supported by the 509 24 National Science Centre (NCN) in Poland through the Grant No. 2021/43/D/ST8/02555. Maciej 510 Ry´ s acknowledges the support from the European Union Horizon 2020 research and innovation 511 program under Grant Agreement No. 857470 and from the European Regional development Fund 512 under the program of the Foundation for Polish Science International Research Agenda PLUS, 513 grant No. MAB PLUS/2018/8, and the initiative of the Ministry of Science and Higher Education 514 ’Support for the activities of Centers of Excellence established in Poland under the Horizon 2020 515 program’ under agreement No. MEiN/2023/DIR/3795. For the purpose of Open Access, the authors 516 have applied a CC-BY public copyright license to any Author Accepted Manuscript (AAM) version 517 arising from this submission. 518 Appendix A. A simplified 1D demonstration of the model In this appendix, we derive the governing equation of the transformation stress for a simplified 519 isothermal 1D model. In 1D setting, the model features four variables, namely the total strain 520 ε=∇u, the reversible volume fraction ηrev, the irreversible volume fraction ηirand the permanent 521 strain εp. The Helmholtz free energy is thus expressed as follows 522 ϕ(ε, εp, ηrev, ηir) =ϕchem(ηrev, ηir) +ϕel(ε, εp, ηrev) +ϕint(ηrev) +ϕdeg(ηrev, ηir) +I(ηrev),(A.1) where the indicator function Ipertains to the inequality constraints on the reversible volume fraction 523 ηrev(I= 0 if 0 ≤ηrev≤1−ηirandI=∞otherwise). Note that the gradient energy associated 524 with the austenite–martensite diffuse interface, ϕgrad, is disregarded here. 525 The elastic strain energy ϕelis formulated as 526 ϕel(ε, εp, ηrev) =1 2E(ε−εt−εp)2, εt=ηrevϵT, (A.2) where Eis the Young’s modulus (for simplicity, Eis assumed constant and independent of η) and 527 the constant ϵTis the maximum transformation strain. The remaining components of the free 528 energy, as well as the dissipation potential, are identical to those of the general 3D model, see 529 Eqs. (14), (16), (18) and Eq. (23). Moreover, the evolution equations for the permanent strain εp530 and the irreversible volume fraction ηirare postulated as (cf. Eqs. (5)–(7)) 531 ηir=hsat ir(1−exp(−Cpηacc)), εp=ϵsat p(1−exp(−Cpηacc)). (A.3) For a given total strain ε, the volume fraction ηrevcan be determined by minimizing the local 532 potential π= ∆ϕ+∆D, cf. Eq. (24). It is immediate to see that the local potential πis non-smooth, 533 due to the presence of the rate-independent dissipation ∆ Dand the indicator function I. In line 534 with [40], the minimization of πwith respect to ηrevis written as a differential inclusion, given by 535 fηrev∈∂ηrev¯D(ηrev, ηacc) (A.4) where ¯D= ∆D+Iencompasses the non-smooth components of πandfηrevis the thermodynamic 536 25 driving force associated with ηrevand is expressed as 537 fηrev=−∂ϕ ∂ηrev+∂ϕ ∂εp∂εp ∂ηrev+∂ϕ ∂ηir∂ηir ∂ηrev . (A.5) During the forward/backward transformation, i.e., when the bound constraints are inactive, the 538 inclusion (A.4) yields 539 fηrev=±fc, (A.6) and gives the following equation for the transformation stress σt ±(σt +for the forward transformation 540 andσt −for the backward transformation), 541 σt ±=∆ϕ0k1±fc+Hintηrev+Adegk2+Hdegηrevk3 k4, (A.7) where fcis defined in Eq. (20) and kiare expressed as 542 k1= 1 +∂ηir ∂ηrev, k 2=ηir+ηrev∂ηir ∂ηrev, k 3=ηir+1 2ηrev∂ηir ∂ηrev, k 4=ϵT+∂εp ∂ηrev. (A.8) It is important to highlight that the necessary condition for the minimum of πwith respect 543 toηrev, which leads to the transformation criteria (A.7), is not computed in a standard manner. 544 This stems from the state-dependence of the dissipation potential D, i.e., the dependence of fcon 545 the accumulated volume fraction ηacc, see Eq. (20). Having the minimization problem formulated 546 in rates (not shown here), it becomes apparent that fcis treated as a constant when evaluating 547 the necessary condition for the rate ˙ ηrev. In the incremental setting, to maintain consistency with 548 the rate-problem, the increment of the martensite volume fraction, ∆ ηrev, present in the current 549 unknown ηrev= ∆ηrev+ηrev nis distinguished from the increment upon which the evolution equation 550 forfcrely. Despite the two increments coincide, the latter is considered as constant when evaluating 551 the necessary condition. Accordingly, the minimization problem does posses the structure of a quasi- 552 optimization problem and not a genuine optimization problem. To avoid the complexity in the model 553 presentation, this issue is not elaborated here. It should be remarked that upon assuming the same 554 increment ∆ ηrevfor the current unknown ηrevand the evolving parameter fc, extra differentiation 555 terms arise in the transformation criteria (A.7). However, our auxiliary simulations showed that 556 these extra terms only marginally contribute to the results. 557