Finite strain constitutive modeling for shape memory alloys considering transformation-induced plasticity and two-way shape memory effect
Lei Xu, Alexandros Solomou, Theocharis Baxevanis, Dimitris Lagoudas
FFinite strain constitutive modeling for shape memory alloysconsidering transformation-induced plasticity and two-wayshape memory effect
Lei Xu a , Alexandros Solomou a , Theocharis Baxevanis c , Dimitris Lagoudas a,b a Department of Aerospace Engineering, Texas A&M University, College Station, TX77843, USA b Department of Materials Science & Engineering, Texas A&M University, College Station, TX77843, USA c Department of Mechanical Engineering, University of Houston, Houston, TX77204, USA
Abstract
This work presents a three-dimensional constitutive model for shape memory alloys consider-ing the TRansformation-Induced Plasticity (TRIP) as well as the Two-Way Shape MemoryEffect (TWSME) through a large deformation framework. The presented logarithmic strainbased model is able to capture the large strains and rotations exhibited by SMAs undergeneral thermomechanical cycling. By using the martensitic volume fraction, transforma-tion strain, internal stress, and TRIP strain tensors as internal state variables, the modelis capable to capture the stress-dependent TRIP generation when SMAs are subjected to amultiaxial stress state, as well as the TWSME for thermomechanically trained SMAs underload-free conditions. A detailed implementation procedure of the proposed model is presentedthrough a user-defined material subroutine within a finite element framework allowing forsolving different Boundary Value Problems (BVPs). Comprehensive instruction on calibrat-ing the model parameters as well as the derivation of continuum tangent stiffness matrixare also provided. In the end, the simulated cyclic pseudoelastic and actuation responsesby the presented model for a wide range of SMA material systems under both uniaxial andmultiaxial stress states are compared against experimental results to validate the proposedmodeling capabilities.
Keywords:
Shape memory alloys, Large deformation, Constitutive modeling,Transformation-induced plasticity, Two-way shape memory effect ∗ Corresponding author
Email address: [email protected] (Lei Xu)
Preprint submitted to ArXiv January 15, 2020 a r X i v : . [ c ond - m a t . m t r l - s c i ] J a n . INTRODUCTION Shape Memory Alloys (SMAs) represent an active/smart material with the ability torecover their pre-defined shape via a diffusionless phase transformation between its high-symmetry, high-temperature austenitic phase and low-symmetry, low-temperature marten-sitic phase. Due to the high output energy density of SMAs, up to 500 MPa actuation stressand 8% recoverable strain (Otsuka & Wayman 1999), compared to other active materialssuch as shape memory polymers and piezoelectrics, their current and potential applicationsin the biomedical, aerospace, automobile and civil engineering fields are expanding rapidly(Hartl & Lagoudas 2007, Song et al. 2010, Peraza-Hernandez et al. 2013, Jani et al. 2014,Karakalas et al. 2019). SMAs have been extensively researched as solid-state actuators toenable adaptive and morphing structures. For examples, an SMA-based beam componenthas been used as a bending actuator to morph the engine outer shell geometry so that a de-sired aerodynamic conditions can be achieved during the airplane take-off and cruise regime(Hartl & Lagoudas 2007). SMA-based torque tubes have been used as rotation actuators todeploy and retract solar panels for small satellites (Wheeler et al. 2015), and also used asrotational actuators to rotate the trailing edge wing flap during an airplane on-fly test (Mabeet al. 2014, Calkins & Mabe 2016).The majority of engineering applications require SMAs experiencing a large number ofloading cycles involving repeated phase transformations, which brings the increasing necessityto understand SMA material response under thermomechanical cycling. Many experimentalresults (Strnadel et al. 1995 a , b , Lagoudas & Bo 1999, Lagoudas & Entchev 2004, Wheeleret al. 2014 a ) indicate that SMAs exhibit an evolving rather than stable material response un-der cyclic loading. More specifically, transformation characteristics of SMAs, e.g., the shapeof stress/thermal hysteresis, transformation temperatures, transformation strain magnitude,usually shift from one cycle to another, and large irrecoverable strains are usually accumu-lated. Such irrecoverable strains, often called TRIP, are caused by the microstructure changesas a result of the repeated phase transformation. These microstructure changes, including ac-cumulation of dislocation bands, initiation and growth of micro-voids and micro-cracks, anddamage accumulation, then effectively results in an observable macroscopic plastic strain,which occurs at an effective stress level much lower than the conventional plastic yielding2oint (Lagoudas & Entchev 2004). In addition, TRIP strain evolves with different ratesthroughout the entire material fatigue life state. It has been shown in the SMA fatiguetest results (Wheeler et al. 2014 b ) that the TRIP strain grows drastically during the veryfirst hundreds of loading cycles then tends to increase in a stabilized trend until the ma-terial reaches the failing point at the end. As the most of engineering applications requireactuators functioning in stable material behavior, SMAs are usually subjected to a trainingprocess, i.e., repeated thermal/stress cycling, to stabilize their behavior before being used asactuation components. As of now, many modeling efforts have been devoted to predictingthe evolving characteristics of SMAs under cyclic thermomechanical loading.A large number of legacy models have been proposed to predict the stable SMA mate-rial response. A thorough review of these works can be found from literature (Leclercq &Lexcellent 1996 a , Boyd & Lagoudas 1996, Patoor et al. 1996, Thamburaja & Anand 2001,Patoor et al. 2006, Zaki & Moumni 2007, Lagoudas 2008, Arghavani et al. 2011, Lagoudaset al. 2012, Kelly et al. 2016, Wang et al. 2017, Xu et al. 2017 b , 2018). In general, SMAmodels considering irrecoverable strains can be categorized into two types. One type of SMAmodels describe conventional plasticity due to the activation of slip systems at sufficientlyhigh-stress levels in either pure austenite or martensite phase. Modeling efforts fall into thistype can be obtained from publication (Wang et al. 2008, Hartl & Lagoudas 2009, Jiang &Landis 2016, Scalet et al. 2019). Another type of models concerns irrecoverable strains asTRIP caused by repeated phase transformations wherein the stress levels are much lowerthan the material yielding point. It is noted that the focus of this work falls into the secondtype. Many commonly cited SMA models have been proposed to capture such evolving re-sponse feature, and a subset of them are briefly reviewed here. The earliest models dealingwith TRIP were presented by Lim & McDowell (1994) and Tanaka et al. (1995) to capturethe cyclic loading effect on the SMA phase transformation characteristics. Later on, basedon the micromechanics averaging method on an SMA representative volume element, Bo &Lagoudas (1999 a , b ) proposed a model accounting for the one-dimensional TRIP strain ac-cumulation and the generation of TWSME under actuation cycling. Lexcellent et al. (2000)further extended their early SMA model describing stable material response (Leclercq &Lexcellent 1996 b ) to consider the irrecoverable strains by introducing two additional internal3tate variables, i.e., the volume fraction of self-accommodated and oriented martensite. Af-terward Lagoudas & Entchev (2004) proposed a model accounting for TRIP as well as theshape and size of the hysteresis during pseudoelastic cycling. Also, Zaki & Moumni (2007)proposed a model considering the TRIP in the case of cyclic pseudoelastic loading by usingadditional internal state variables, such as internal stress, TRIP strain, and accumulatedmartensite volume fraction. Other similar and recent modeling efforts can also be foundfrom the works (Auricchio et al. 2007, Saint-Sulpice et al. 2009, Yu et al. 2013, Barrera et al.2014, Yu et al. 2015, Xu et al. 2017 a , 2019 b ).Although many of the proposed models have enabled researchers to study the evolvingmaterial behaviors of SMAs, the majority of them are insufficient in their capacity to considerthe following critical features. (i) The first feature of the currently available models in need ofimprovement is their small deformation assumption based on infinitesimal strain theory. Thisassumption may be acceptable for SMA material systems, such as Ni-rich or NiTiHf SMAs,where the summation of all strains, including elastic, transformation, and TRIP strains, isbelow 3%. However, it has been reported that nearly 30% or even higher TRIP strains areobserved during the lifetime of near equiatomic NiTi SMA-based actuators (Wheeler et al.2014 a ). Also with the presence of cracks in SMAs, the strain regime in front of the cracktip can easily go beyond 10% (Haghgouyan et al. 2016, 2019). In the presence of such largestrain, a finite strain model is needed to account for the exhibited large strains to providean accurate structural response of SMA-based multifunctional components. (ii) The secondimportant aspect of many current models in need of improvement is the TWSME charac-teristic exhibited by trained SMAs at load-free conditions. Because of the required trainingprocess to stabilize the response of as-received SMAs before used as actuators, permanentchanges such as dislocation bands, accumulation of defects/damage and retained marten-site variants are usually introduced into the material’s microstructure, which then results inthe generation of an internal stress field oriented in the same direction as the applied load.Thereafter the generated internal stress field is capable to induce the oriented phase transfor-mation under pure thermal cycling without applying any external mechanical loads, i.e., theTWSME. Such unique TWSME property of SMAs have tremendous engineering potentials.For instance, it allows for mounting and dismounting of SMA-based connectors and couplers4 train Figure 1: Evolution of the summation of transformation and TRIP strain in the loading direction for anotched NiTi plate subjected to thermal cycling under constant load. Referenced from Wheeler (2017). in an easy procedure by just heating and cooling without pre-stressing (Niccoli et al. 2017,Tabesh et al. 2018). (iii) The last major contribution of this work is the consideration ofstress-dependent TRIP evolution under multiaxial stress state. The majority of applicationsrequire the functionality of actuators under multiaxial stress state originated from geometrycurvatures or installment required discontinuities, such as notches and holes. TRIP strainunder such multiaxial stress state evolves quite differently compared to that in the uniaxialloading case. Refer to Fig. 1 for the Digital Image Correlation (DIC) strain results of anotched NiTi plate under cyclic actuation loadings, it revealed that the TRIP strain evolvedin a much faster rate at the stress concentration region than the less stressed part. Despitethe importance of the above mentioned facts, they have rarely been addressed among existingmodels from available literature through a unified modeling framework.In order to address the three aforementioned critical features, a three-dimensional fi-nite strain constitutive model for SMAs considering TRIP and TWSME is proposed in thiswork. The presented modeling effort is developed based on the legacy SMA model (Boyd &Lagoudas 1996, Lagoudas et al. 2012, Xu et al. 2019 a ) and largely inspired by its continuousdevelopment considering TRIP (Bo & Lagoudas 1999 b , Lagoudas & Entchev 2004, Xu et al.2019 b ). The proposed logarithmic strain based model is able to capture the large strainsand rotations exhibited by SMAs under general thermomechanical cycling. By using the5artensitic volume fraction, transformation strain, internal stress, and TRIP strain tensorsas internal state variables, the model is also able to capture the stress-dependent TRIP evolu-tion when SMAs are under a multiaxial stress state, and the TWSME for thermomechanicallytrained SMAs due to the generation of internal stresses under load-free conditions.In summary, the paper is organized as follows. In Section 2, kinematic preliminariesused in the model formulation are presented. Section 3 focuses on the model developmentthat incorporates the three important features as mentioned earlier. In section 4, the de-tailed implementation procedure for the proposed model is described by using a user-definedmaterial subroutine (UMAT) through the finite element software Abaqus. Thereafter, nu-merical examples are analyzed to demonstrate the proposed modeling capabilities in Section5. Conclusions are summarized in Section 6. In the Appendix A, the explicit derivationof the continuum tangent stiffness matrix is provided, and Appendix B presents a detailedinstruction on how to calibrate all the model parameters based on available experimentaldata.
2. PRELIMINARY
Let the position of a material point P from body B defined by a vector X in the reference(undeformed) configuration at time t , and vector x represents the position of that materialpoint in the current (deformed) configuration at time t . It is well-known that the deformationprocess of this material point P from reference to current configurations can be described byusing the deformation gradient tensor F ( x , t ), F ( x , t ) = ∂ x ∂ X (1)also, the velocity field v of material point P can be defined as, v = d x dt = ˙ x (2)the velocity gradient L can be derived based on the velocity field v as, L = ∂ v ∂ x = ˙FF − (3)6he polar decomposition for the deformation gradient F is also known as follows, F = UR = VR (4)in which R is the orthogonal rotation tensor, U and V are called the right (or Lagrangian)and the left (or Eulerian) stretch tensors respectively. The right Cauchy-Green tensor C andthe left Cauchy-Green tensor B can be obtained as follows, C = F T F = U (5) B = FF T = V (6)the so called logarithmic (or Hencky/true) strain can be expressed in its Lagrangian form H and its Eulerian form h as, H = 12 ln( C ) = ln( U ) (7) h = 12 ln( B ) = ln( V ) (8)It is also known that the velocity gradient L can be additively decomposed into a sym-metric part, i.e., the rate of deformation tensor D , plus an anti-symmetric part, i.e., the spintensor W , L = D + W , D = 12 ( L + L T ) , W = 12 ( L − L T ) (9) Two kinematic assumptions, i.e., the multiplicative decomposition of deformation gra-dient F and the additive decomposition of the rate of deformation tensor D , are usuallyused in finite deformation theory. Hyperelastic constitutive relation is often used in multi-plicative models while hypoelastic constitutive equation is utilized for additive models. Fora long time, the rate form hypoelastic constitutive theory has been criticized for its fail-ure to be fully integrable to describe a simple recoverable elastic behavior. Many spurious7henomena, such as shear stress oscillation and residual stress errors are observed in simpleelastic deformation using objective rates, this includes many well-known objective rates suchas Zaremba-Jaumann rate, Green-Naghdi rate, and Truesdell rate, etc.(Xiao et al. 2006).However, such aforementioned issues regarding objective rates are resolved by the logarith-mic rate proposed by Xiao et al. (1997 a , b , 2006), Bruhns et al. (1999, 2001 a , b ), Meyers et al.(2003, 2006). As proved in their work, the logarithmic rate of Eulerian logarithmic strain h isidentical with the rate of deformation tensor D , by which a hypoelastic model can be exactlyintegrated into an finite strain elastic model (Xiao et al. (1997 a )). This unique relationshipbetween the logarithmic strain h and the rate of deformation tensor D can be expressed as,˚ h log = ˙ h + hΩ log − Ω log h = D (10)where Ω log is the logarithmic spin introduced by Xiao et al. (1997 a ), Ω log = W + n (cid:88) i (cid:54) = j (cid:0) λ i /λ j )1 − ( λ i /λ j ) + 2ln( λ i /λ j ) (cid:1) b i Db j (11)in which λ i,j ( i, j = 1 , ,
3) are the eigenvalues of left Cauchy-Green tensor B , b i , b j are thecorresponding subordinate eigenprojections. Additionally, the second-order rotation tensor R log can be obtained by using the logarithmic spin tensor Ω log after solving the followingtensorial differential equation. The initial condition R log | t =0 = I is often assumed. Ω log = ˙ R log ( R log ) T (12)Furthermore, equation (10) can be integrated to the following equation via the corota-tional integration scheme with the initial condition h | t =0 = , (Khan & Huang (1995)), h = (cid:90) corot. D d t = ( R log ) T (cid:18) (cid:90) t R log D e ( R log ) T d t (cid:48) (cid:19) R log (13) This part addresses the kinematic assumption of additive decomposition of logarithmicstrain. First, the rate of deformation tensor D is additively decomposed into three parts, The thermal strain part is small compared to the other major strains, hene it is omitted here for thesake of simplicity. D can be interpreted as the stress power provided from outside is splitinto an elastic part stored inside the material, a dissipative part associated with phase trans-formation process, and another dissipative part related to transformation-induced plasticdeformation. D = D e + D tr + D tp (14)Based on equation (10), the elastic part D e , transformation part D tr and plastic part D tp can be rewritten as ˚ h e log , ˚ h tr log and ˚ h tp log respectively,˚ h e log = D e ; ˚ h tr log = D tr ; ˚ h tp log = D tp (15)the following equation can be obtained by combining equations (10), (14) and (15),˚ h log = ˚ h e log + ˚ h tr log + ˚ h tp log (16)applying corotational integration on equation (16) gives, h e = (cid:90) corot. D e d t = ( R log ) T (cid:18) (cid:90) τ R log D e ( R log ) T d τ (cid:19) R log (17a) h tr = (cid:90) corot. D tr d t = ( R log ) T (cid:18) (cid:90) τ R log D tr ( R log ) T d τ (cid:19) R log (17b) h tp = (cid:90) corot. D tp d t = ( R log ) T (cid:18) (cid:90) τ R log D tp ( R log ) T d τ (cid:19) R log (17c)Combing equations (13), (14), (16) and (17), the following kinematic equation on totallogarithmic strain can be received. Namely, the total logarithmic strain is additively splitinto an elastic part, a transformation part, as well as a TRIP part h = h e + h tr + h tp (18)
3. MODEL FORMULATION
Based on the classical thermodynamic framework for dissipative materials, the devel-opment of the proposed model starts with the formulation of an explicit thermodynamic9otential. To that end, a quadratic form Gibbs free energy is proposed as a continuousfunction based on Kirchhoff stress tensor τ , temperature T , and a set of internal state vari-ables Υ = { ξ, h tr , h tp , β } , in which they are martensitic volume fraction ξ , transformationstrain tensor h tr , TRIP strain tensor h tp , and internal stress tensor β . The martensiticvolume fraction ξ ranging 0 (cid:54) ξ (cid:54) ξ = 0 represents pure austenitic phase while ξ = 1 indicates puremartensitic phase. The h tr accounts for the inelastic yet recoverable strain associated withthe phase transformation, h tp is used to represent the irrecoverable transformation-inducedplastic strain, and β is used to consider the internal stress field generated inside the materialas a result of the training process. The following explicit Gibbs free energy expression for G is given as, G = − ρ τ : S τ − ρ τ : [ α ( T − T ) + h tr + h tp ] − ρ (cid:90) ξ ( β : ∂ h tr ∂τ ) dτ + c (cid:104) ( T − T ) − T ln( TT ) (cid:105) − s T + u + 1 ρ f ( ξ ) (19)In which, S is the fourth-order effective compliance tensor that can be calculated by usingthe rule of mixture as per equation (20), S A is the compliance tensor for austenitic phasewhile S M is for martensitic phase, ∆ S represents the phase difference of the compliancetensor. Additionally, the effective stiffness tensor C can be obtained by taking the inverseof the above effective compliance tensor, i.e., C = S − . α is the second-order thermoelasticexpansion tensor, c is the effective specific heat, s and u are the effective specific entropyand effective specific internal energy at the reference state. All the aforementioned effectivevariables are determined by the rule of mixture from equation (21) to (24). T represents thetemperature at current state, while T is the temperature at reference state. S ( ξ ) = S A + ξ ( S M − S A ) = S A + ξ ∆ S (20) α ( ξ ) = α A + ξ ( α M − α A ) = α A + ξ ∆ α (21) c ( ξ ) = c A + ξ ( c M − c A ) = c A + ξ ∆ c (22)10 ( ξ ) = s A + ξ ( s M − s A ) = s A + ξ ∆ s (23) u ( ξ ) = u A + ξ ( u M − u A ) = u A + ξ ∆ u (24)A smooth hardening function f ( ξ ) is included in the Gibbs free energy to account for thepolycrystalline hardening effects, such as interactions between different phase variants, im-perfections located at the grain boundaries, and the presence of nano-precipitates (Lagoudas2008). Three intermediate material parameters ( a , a , a ) are introduced in this hardeningfunction, they can be determined as derived in Appendix B by using the known materialparameters . Besides, the other four smoothing parameters ( n , n , n , n ) are introduced toeffectively treat the smooth transition characteristics at the initiation and completion duringphase transformation. The complete form of hardening function is as follows, f ( ξ ) = a (cid:16) ξ + ξ n n +1 + (1 − ξ ) n n +1 (cid:17) + a ξ , ˙ ξ > , a (cid:16) ξ + ξ n n +1 + (1 − ξ ) n n +1 (cid:17) − a ξ , ˙ ξ < h = − ρ ∂G∂ τ = S τ + α ( T − T ) + h tr + h tp (26) s = − ρ ∂G∂T = 1 ρ τ : α + c ln( TT ) − s (27)The following reduced form of dissipation inequality (28) can be derived by substitutionof the above constitutive relationships into the Clausius-Planck inequality, − ρ ∂G∂ h tr : ˚ h tr − ρ ∂G∂ h tp : ˚ h tp − ρ ∂G∂ξ ˙ ξ (cid:62) .2. Evolution law for transformation strain This part focuses on proposing the evolution law for transformation strain. Following themaximum dissipation principle, the evolution of transformation strain that tends to dissipatethe most energy among all the admissible thermodynamic paths is chosen during the phasetransformation process (Boyd & Lagoudas 1996, Qidwai & Lagoudas 2000 b ). Therefore, it isassumed that the rate change of transformation strain is proportional to the rate change ofthe martensitic volume fraction ξ , and the direction of which is along the deviatoric part ofthe effective stress tensor. Be noted that in the following evolution law the rate applied ontop of the transformation strain is the logarithmic rate in order to account for the principleof objectivity, the ’log’ symbol is neglected for the sake of simplicity. Finally, the explicitevolution law for transformation strain is described as follows,˚ h tr = Λ ˙ ξ, Λ = Λ fwd , ˙ ξ > , Λ rev , ˙ ξ < , (29)where, Λ fwd is the forward transformation direction tensor, and Λ rev is the reverse transfor-mation direction tensor. They are defined explicitly as follows, Λ fwd = 32 H cur τ eff (cid:48) ¯ τ eff ; Λ rev = h tr-r ξ r . (30)in the above equations, τ eff is the effective stress tensor defined as equation (31) using thesummation of Cauchy stress and internal stress. It is worthy to point out that the majordifference here compared to available SMA model concerning stable material response is theintroduction of this effective stress tensor. It can be seen later on that the TWSME ofthermomechanically trained SMAs under load-free conditions can be achieved by using thiseffective stress term. The evolution law for internal stress tensor β is introduced in the latercontext shortly. τ eff = τ + β (31)in addition, the deviatoric part of effective stress tensor is defined as τ eff (cid:48) = τ eff − tr( τ eff ) ,where is the second-order identity tensor. The von Mises equivalent effective stress ¯ τ eff iscalculated as follows,¯ τ eff = (cid:114) τ eff (cid:48) : τ eff (cid:48) (32)12t is common in many of the available models that the magnitude of the recoverable trans-formation strain is assumed to be the same under different stress levels. Such considerationis valid when the applied stress levels are high enough to generate fully oriented marten-sitic variants. However, self-accommodated martensitic variants are also generated when thestress levels are not sufficiently high, which renders the value of transformation strain to bea stress-dependent variable. Therefore, the following exponential H cur function is introducedto calculate the value of current transformation strain given an effective stress state, where H max is the maximum (or saturated) value of transformation strain, H min corresponds toan observable TWSME strain for pre-trained SMAs or some SMAs experiencing particularproduction process such as extrusion and aging under stress. Besides, τ crit denotes a criticalstress value below which H cur = H min , and k t is a curve-fitting material parameter. Theexplicit form of H cur can be found as follows, H cur (¯ τ eff ) = H min + ( H max − H min )(1 − e − k t (¯ τ eff − τ crit ) ); ¯ τ eff > τ crit ,H min ; ¯ τ eff < τ crit , (33)It is also seen from experimental results (Atli et al. 2015) that a degradation sometimesexists for the value of maximum transformation strain as a result of the accumulation ofretained martensite. In order to extend the model capability to capture this phenomenon,a degradation law is proposed for the maximum transformation strain in equation (34),where H maxi and H maxf represent the value of H max before and after the cyclic loading. Inaddition, λ is a material parameter governing the degradation trend of H max , and ζ d calledthe accumulation of orientated martensitic volume fraction is introduced shortly in the latersubsection. H max = H maxf + ( H maxi − H maxf ) e − λ ζ d (34) An evolution law for the TRIP strain is proposed in this subsection in order for themodel to capture the irrecoverable strain exhibited by SMAs under cyclic thermomechanicalloading. Before the detailed formulation is presented, a major assumption is postulated, i.e.,among the total martensitic phase transformation, only the oriented transformation portion13ontributes to the generation of TRIP. This assumption is built upon the observation onthe experimental results (Lagoudas 2008) that no macroscopic irrecoverable deformationsare perceived for untrained SMAs under load-free thermal cycling. An early one-dimensionalform of TRIP evolution law was suggested in the work of Bo & Lagoudas (1999 b ), and athree-dimensional form was proposed by Lagoudas & Entchev (2004). As it was discussed inthe introduction, the generation of TRIP strain is highly stress-dependent under multiaxialstress state, and is totally different compared to that in the uniaxial case. However, none ofthe above mentioned TRIP evolution laws have well addressed this critical feature, nor theother available models to the authors’ best knowledge. Therefore, the following evolutionlaw for TRIP strain considering the effect of stress multiaxiality is suggested,˚ h tp = Λ tp ˙ ξ, Λ tp = Λ tpfwd , ˙ ξ > , Λ tprev , ˙ ξ < , (35)In the above equation, Λ pfwd is the forward TRIP direction tensor, and Λ prev is the reverseone. They have explicit definitions in equation (36). Note that the rate applied on top of h tp is again the logarithmic rate. Because TRIP is strongly driven by the phase transformations,it is proposed that the TRIP strain also evolves along the deviatoric part of the effectivestress as the transformation strain described as the following equation, Λ tpfwd = 32 ( H cur H max ) τ (cid:48) eff ¯ τ eff C p C p C p ζ d ; Λ tprev = − ( H cur H max ) h tr − r ξ r C p C p C p ζ d (36)In the preceding evolution equation, the material parameters C p and C p play the majorrole in dictating the magnitude and evolving trend for TRIP during cyclic thermo-mechanicalloading. The usage of the accumulation of oriented martensitic volume fraction ζ d is to beconsistent with the previously proposed assumption that only the oriented phase transfor-mation contributes to the TRIP strain generation. For a better illustration on the TRIPevolving trend, the rate form equation (36) can be integrated as the following algebraic formusing a logarithmic function. h tp = 32 H cur H max C p τ (cid:48) eff ¯ τ eff ln (cid:0) C p ζ d (cid:1) (37)As it can be seen from the above equation, the stress-dependent feature for TRIP straingeneration is incorporated by multiplying model parameter C p with the ratio ( H cur / H max ).14s it is shown in the latter results section, such consideration together with the usage of ζ d inthe evolution law enable the model’s capability to capture the highly stress-dependent TRIPgeneration under multiaxial stress condition. It is also worth to point out that the proposedlogarithmic function based evolution law indicates there is no saturation point on TRIP, whichis closely aligned with the experimental results from actuation cycling (Wheeler 2017). Whilein the previous work (Lagoudas & Entchev 2004), an exponential function based evolution lawwas adopted to be aligned with pseudoelastic cycling experimental results, which indicatesTRIP is expected to reach a saturation value, and eventually a stable pseudoelastic responseis anticipated to be attained after a certain number of mechanically cycling. Lastly, theaccumulation of oriented martensitic volume fraction ζ d used in this equation is defined asfollows, ζ d = (cid:90) t | ˙ ξ d ( t ) | dt (38)in which the oriented martensitic volume fraction ξ d is calculated as, ξ d = H cur H max ξ (39)the relation between the accumulation of oriented martensitic volume fraction ζ d and theaccumulation of total martensitic volume fraction ζ is, ζ d = H cur H max ζ (40) In order to meet the proposed goal of capturing the TWSME at load-free conditions forthermomechanically trained SMAs, a second-order internal stress tensor is introduced in thismodel. As discussed to some extents in the introduction section, during cyclic thermome-chanical loading SMAs experience microstructure changes, such as pileup of dislocation bandsdriven by phase transformations, initiation and growth of micro-voids and micro-cracks,and damage accumulation, at a stress level below the yielding point. These microstructurechanges gradually introduce a local stress field inside the SMA materials, which is describedby the proposed internal stress tensor via an effective manner. It is reasonably assumed theinternal stress may never go beyond the material yielding point, thus its magnitude is satu-rated at a maximum point after a certain number of loading cycles. Besides, the evolution15irection of the internal stress is determined as the same direction of the applied stress fielddue to the fact that the induced microstructure changes are caused by the external mechan-ical loads. On the basis of the above discussions, the following exponential evolution law isproposed for the internal stress as, β = σ b H cur H max τ eff ¯ τ eff (1 − e − λ ζ d ) (41)in which, σ b is a model parameter representing the maximum (or saturated) magnitude ofthe internal stress. Similar as the TRIP evolution law, the ratio ( H cur / H max ) and the term ζ d are included here to consider the stress-dependent effect on the internal stress evolution,and λ is a material parameter controlling the evolution trend. In this subsection, a transformation function and an associated transformation criterionare defined, upon which the initiation and completion of phase transformation can be deter-mined. By substitution of the evolution law for transformation strain (30) and TRIP strain(35) into the reduced form of dissipation inequality (28), the following equation is obtained, (cid:0) τ : Λ + β : Λ + τ : Λ tp − ρ ∂G∂ξ (cid:1) ˙ ξ = π ˙ ξ (cid:62) π is called the general thermodynamic driving force conjugated to ξ in the proposed model. The product by substitution of Gibbs free energy (19) into theabove equation (42) yields the explicit expression for π in the following equation, where ∆ S , ∆ α , ∆ c, ∆ s , ∆ u each represents the phase difference between martenstie and austenitefor that variable. π = ( τ + β ) : Λ + τ : Λ tp + 12 τ : ∆ S τ + τ : ∆ α ( T − T ) + ρ ∆ s T − ρ ∆ c (cid:2) T − T − T ln( TT ) (cid:3) − ρ ∆ u − ∂f∂ξ (43)Following the development in the early models (Boyd & Lagoudas 1996, Lagoudas et al.2012), it is assumed that the forward (reverse) phase transformation initiates whenever thethermodynamic driving force π ( − π ) reaches a critical value Y ( − Y ), and π ( − π ) is al-ways below such critical value as long as the forward (reverse) phase transformation is not16ompleted, upon which the following transformation function Φ is defined,Φ = π − Y, ˙ ξ > , − π − Y, ˙ ξ < , (44)A further improvement in the critical value Y is suggested from Lagoudas et al. (2012),see equation (45), wherein Y is constructed to be stress-dependent quantity with a referencecritical value Y and an additional parameter D . Such consideration enables the model tocapture the different size of hysteresis loops when SMA materials are subjected to various ap-plied stress levels. This capability is provided through capturing the different slopes C A , C M in the effective stress-temperature phase diagram. Y ( σ ) = Y + D σ : Λ fwd , ˙ ξ > ,Y + D σ : Λ rev , ˙ ξ < , (45)In order for the defined transformation functions and the corresponding transformationcriteria to satisfy the principle of maximum dissipation, the following Kuhn-Tucker constraintconditions have to be met,˙ ξ (cid:62)
0; Φ( τ , T, ξ ) = π − Y (cid:54)
0; Φ ˙ ξ = 0;˙ ξ (cid:54)
0; Φ( τ , T, ξ ) = − π − Y (cid:54)
0; Φ ˙ ξ = 0; (46)
4. IMPLEMENTATION
This section discusses the implementation of the above proposed model into a numericalenvironment through a user-defined material subroutine (UMAT), the adopted finite elementtool Abaqus equipped with UMAT is then used for solving different BVPs. As shown in Fig.2, all the UMAT variables used in the this model are presented via a flowchart. This imple-mentation procedure follows the knowledge of Return Mapping Algorithm (RMA) presentedin the available publications (Qidwai & Lagoudas 2000 a , Lagoudas 2008). The general goalof the UMAT is, given an increment of strain and temperature from the finite element solver,to calculate the output stress values by using the constitutive relationships (26) and (27)with all the internal state variables conforming to the Kuhn-Tucker consistency constraints1746). The implementation in general consists of two major procedures, one is called theThermoelastic-predictor, and the other is called the Transformation-corrector. It is worthyto point out that the input variables used for this implementation collected from the globalfinite element solver are the temperatures ( T n , T n +1 ), and deformation gradients at the cur-rent and next step ( F n , F n +1 ). As it was discussed in the work of Xu et al. (2019 a ), suchnon-trivial considerations on the kinematics allow for the implementation of the model toget rid of the accumulated stress errors as a result of using other non-integrable objectiverates. The effects of such accumulated stress errors on the cyclic thermomechanical responseof SMAs are analyzed in detail from Xu et al. (2019 a ).18 quilibrium Solver Kinematics Global FE Solver
Pre-Calculation
Thermoelastic predictor
Rotation Procedure
Transformation correctorUpdate variables Continuum tangent stiffness matrix
Figure 2: UMAT flowchart for all the variables used in the model during the implementation within the finiteelement framework. able 1: The summary of the UMAT implementation. Initialization • Conduct pre-calculation and rotation procedures. • k = 0; ξ (0) n +1 = ξ n ; h tr (0) n +1 = h trn ; h tp (0) n +1 = h tpn ; β (0) n +1 = β n Thermoelastic Predictor • τ (0) n +1 = C (0) n +1 [ h n +1 − h tr (0) n +1 − h tp (0) n +1 − α ( T n +1 − T )] • Calculate Φ ( k ) n +1 . • IF Φ (0) n +1 (cid:54) tol , GOTO 4 (thermoelastic response). • IF Φ (0) n +1 > tol , GOTO 3 (transformation happens).3. Transformation Corrector • Calculate residual matrix R tr ( k ) n +1 = − h tr ( k ) n +1 + h trn + Λ ( k ) n +1 ( ξ ( k ) n +1 − ξ n ) R tp ( k ) n +1 = − h tp ( k ) n +1 + h tpn + Λ tp ( k ) n +1 ( ξ ( k ) n +1 − ξ n )Φ ( k ) n +1 = Φ( τ ( k ) n +1 , T n +1 , ξ ( k ) n +1 ) • Perform Newton-Raphson iterations in equation (54) to obtain∆ ξ ( k ) n +1 , ∆ h tr ( k ) n +1 , ∆ h tp ( k ) n +1 . • Update variables ξ ( k +1) n +1 , h tr ( k +1) n +1 , h tp ( k +1) n +1 , S ( k +1) n +1 ξ ( k +1) n +1 = ξ ( k ) n +1 + ∆ ξ ( k ) n +1 h tr ( k +1) n +1 = h tr ( k ) n +1 + ∆ h tr ( k ) n +1 h tp ( k +1) n +1 = h tp ( k ) n +1 + ∆ h tp ( k ) n +1 S ( k +1) n +1 = S A + ξ ( k +1) n +1 ∆ S• IF Φ ( k +1) n +1 (cid:62) tol , GOTO step 3, next local iteration k = k + 1.ELSE GOTO step 4, EXIT4. Calculate continuum tangent stiffness L . • L = C + [ C (∆ S τ + Λ + Λ tp )] ⊗ [ C ∂ τ Φ] ∂ ξ Φ − ∂ τ Φ : C (∆ S τ + Λ + Λ tp )5. Update σ n +1 , ζ dn +1 and β n +1 • σ n +1 = J τ n +1 , wherein τ n +1 is calculated based on equation (26) • ζ dn +1 = ζ dn + H curn +1 H max | ξ n +1 − ξ n |• β n +1 = σ b τ eff n +1 ¯ τ eff n +1 (1 − e − λ ζ dn +1 )6.Exit UMAT and proceed to the global FE solver for the next increment20 pre-calculation and a rotation procedure are employed before calling the thermoe-lastic prediction and transformation correction steps. In the pre-calculation procedure, thetotal strain at current and next step ( h n , h n +1 ) are calculated based on ( F n , F n +1 ) usingequation (8). The incremental rotation tensor ∆ R logn based on the logarithmic rate can becalculated by using the exponential map scheme described in Simo & Hughes (2006). Inthe rotation procedure, the tensorial variables stored as solution-dependent quantities in-cluding h n , h trn , h tpn , β n , Λ n , and Λ tpn are rotated from the previous n th configuration tothe current ( n + 1) th configuration using ∆ R logn , thus, to preserve the so-called principle ofobjectivity. To proceed with the thermoelastic-predictor step, the internal state variables Υ (0) n +1 = { h tr (0) n +1 , h tp (0) n +1 , β (0) n +1 , ξ (0) n +1 } at ( n + 1) th step are assumed to be Υ n for the initial ther-moelastic evaluation. In the case that the Kuhn-Tucker consistency condition is violated,i.e., Φ (0) n +1 (cid:62)
0, the transformation correction procedure is initiated to attain updated inter-nal state variables to regain consistency. Otherwise, the current ( n + 1) th step is detected as athermoelastic response with Υ n +1 = Υ (0) n +1 , and the UMAT skips the transformation correc-tion step and continues the rest procedures. A more detailed description of implementationsteps is presented in the following context. This part is a detailed description of the thermoelastic prediction step. Take the ( n +1) th step for example, the total strain h n +1 and temperature T n +1 are obtained from the Pre-calculation procedure, and the initial internal state variables Υ (0) n +1 are assumed the same as Υ n for the initial consistency evaluation, i.e., h tr (0) n +1 = h trn ; h tp (0) n +1 = h tpn ; β (0) n +1 = β n ; ξ (0) n +1 = ξ n (47)on the basis of the above information, the guessed stress value τ (0) n +1 is calculated through theconstitutive equation (26), τ (0) n +1 = C (0) n +1 (cid:104) h n +1 − h tr (0) n +1 − h tp (0) n +1 − α (0) n +1 ( T n +1 − T ) (cid:105) (48) ( · ) ( k ) represent the local value of that variable at the k th iteration in the transformation correctionprocedure, here k = 0 means that this is an initial guess value in thermoelastic prediction procedure. (0) n +1 in thermoelastic procedure can be evaluatedbased on equations (43) and (44) as follow,Φ (0) n +1 = Φ( τ (0) n +1 , T n +1 , Υ (0) n +1 ) (49)If the transformation consistency constraints are preserved, i.e., the initial value of trans-formation function Φ (0) n +1 (cid:54) , no phase transformation is initiated inside SMAs under thecurrent material state, and the current step is detected as a thermoelastic step. If the trans-formation criterion is violated, i.e. Φ (0) n +1 (cid:62)
0, the transformation correction procedure isactivated in order to restore the consistency constraints (46) via retaining updated internalstate variables Υ ( k ) n +1 . This part focuses on the iterative transformation correction procedure to find a set ofupdated internal state variables to regain the transformation consistency conditions, i.e.,Φ ( k ) n +1 (cid:54)
0. Take the k th iteration as an example, the objective is to solve a system ofnonlinear equations summarized in equations (50), (51) and (52), where the residuals in therate form evolution equations (29) and (35) for transformation strain and TRIP strain canbe reformulated into the discretized form as equations (51) and (52),Φ ( k ) n +1 = Φ( τ ( k ) n +1 , T n +1 , ξ ( k ) n +1 ) (50) R tr ( k ) n +1 = − h tr ( k ) n +1 + h trn + Λ ( k ) n +1 ( ξ ( k ) n +1 − ξ n ) (51) R tp ( k ) n +1 = − h tp ( k ) n +1 + h tpn + Λ tp ( k ) n +1 ( ξ ( k ) n +1 − ξ n ) (52)The goal to regain the consistency condition then becomes to satisfy the following con-vergence inequalities, in which ’tol’ means a small convergence value usually can be set as Usually a small value ’tol’ is used for Φ (0) n +1 (cid:54) ’tol’ evaluation, ’tol’ is acceptable to be 10 − or a evensmaller value. − . | Φ ( k ) n +1 | (cid:54) tol ; | R tr ( k ) n +1 | (cid:54) tol ; | R tp ( k ) n +1 | (cid:54) tol (53)A standard Newton-Raphson iteration procedure can be utilized to solve the abovenonlinear system of equations. The following equation is usually the explicit form to beiterated wherein the first matrix term in the right-hand side is the inverse of the so-calledJacobian matrix for the nonlinear system of equations. ∆ ξ ( k ) n +1 ∆ h tr ( k ) n +1 ∆ h tp ( k ) n +1 = − ∂ Φ ( k ) n +1 ∂ξ ∂ Φ ( k ) n +1 ∂ h tr ∂ Φ ( k ) n +1 ∂ h tp ∂ R tr ( k ) n +1 ∂ξ ∂ R tr ( k ) n +1 ∂ h tr ∂ R tr ( k ) n +1 ∂ h tp ∂ R tp ( k ) n +1 ∂ξ ∂ R tp ( k ) n +1 ∂ h tr ∂ R tp ( k ) n +1 ∂ h tp − Φ ( k ) n +1 R tr ( k ) n +1 R tp ( k ) n +1 (54)During each k th iteration of the Newton-Raphson precedure, the following updated valuesfor the next k + 1 th iteration are obtained for the internal state variables, ξ ( k +1) n +1 h tr ( k +1) n +1 h tp ( k +1) n +1 = ξ ( k ) n +1 h tr ( k ) n +1 h tp ( k ) n +1 + ∆ ξ ( k ) n +1 ∆ h tr ( k ) n +1 ∆ h tp ( k ) n +1 (55)The Newton-Raphson procedure iterates a maximum number of iterations un-til the convergence criterion equation (53) is satisfied. Once the converged values { h tr ( k +1) n +1 , h tp ( k +1) n +1 , ξ ( k +1) n +1 } are accepted as the final value for the current material state atthe k th iteration, the current transformation correction step is completed.
5. RESULTS AND DISCUSSIONS
In this section, four BVPs are analyzed to test the proposed modeling capabilities of thismodel considering TRIP and TWSME. These BVPs include both pseudoelastic and actuation23ycling for a wide range of SMA material systems, incorporating Ni Ti Cu , Ni . Ti . ,Ni Ti , and Ni . Ti . Hf with their compositions all in atomic percentage (at.%). Be-cause the phase transformation characteristics, such as the transformaiton temperatures andstrains, shape and size of hysteresis, TRIP strain accumulations, are very various for differentmaterial compositions, the fidelity of the proposed model is therefore tested over a variety ofSMAs. Three out of the four BVPs are in uniaxial loading conditions while the last one is amultiaxial one.The first BVP is a Ni Ti Cu (at.%) SMA wire under uniaxial cyclic pseudoelasticloading, in which the accumulation of TRIP strain and the stress levels required to initiatethe phase transformation after loading cycles are analyzed. Secondly, two BVPs investigate aNi . Ti . (at.%) and a Ni . Ti . Hf (at%) SMA dogbone specimens under uniaxial cyclicactuation loading. Apart from the accumulated TRIP strain, the load-free TWSME responsesafter the training procedure are also analyzed. Finally, a Ni Ti (at.%) plate actuator witha centric hole is chosen to test the model’s capability to capture the stress-dependent TRIPevolution under multiaxial stress state. It should be noted that the logarithmic (or ture)strain is adopted in both the simulation and experimental results, and the calibration forthe model parameters are based on the true stress-strain curves instead of the engineeringone. Under large strain conditions, using an engineering scale stress-strain curve to calibratethe model is expected to cause some discrepancies on material parameters such as Young’smodulus and transformation strain as discussed in Xu et al. (2019 a ). Lastly, the simulation ofthese BVPs are obtained through the commercial finite element software Abaqus, into whichthe constitutive response of SMAs described above is implemented as a UMAT discussed inthe previous section. The first BVP describes a Ni Ti Cu (at.%) SMA wire under uniaxial pseudoelasticcyclic loading. A bias mechanical load is applied in the longitudinal wire direction from 0 to550 MPa and then unloaded, the temperature is kept constant at 360 K throughout the entireprocedure. Such loading and unloading process is repeated for 50 cycles. The NiTiCu SMAwire experienced a stress-induced forward phase transformation during the loading regimefollowed by a reverse phase transformation at the unloading. The material parameters used24 able 2: Model parameters used for NiTiCu SMA under uniaxial cyclic pseudoelastic loading Type Parameter Value Parameter Value E A
70 [GPa] C A E M
50 [GPa] C M ν A = ν M M s
264 [K]12 α A = α M × − [K − ] M f
160 [K] H max i A s
217 [K] H max f A f
290 [K] H min k t τ crit Smooth hardening parameters n n n n σ b
100 [MPa] λ C p × − C p st cycle to the 50 th cycle. In addition, it can be seen from Fig. 4 that the TRIP strain growsdrastically within the initial 30 loading cycles then stabilizes in a stationary increasing trendthereafter. The accumulated TRIP strain in the 1 st cycle is about 0.6% and grows to around6% after 30 cycles. The predicted TRIP accumulation result is shown in good agreementwith the experimental data. Apart from capturing the irrecoverable TRIP strain, the modelis also able to capture the experimentally observed decreasing stress levels required to initiatethe forward phase transformation. This modeling capability is achieved by the introductionof internal stress term into the model formulation. As the internal stress accumulates in anexponential manner described as evolution equation (41) when the SMA material is subjectedto mechanically cycling, it decreases the stress required from outside to kick off the phasetransformation. In summary, it is shown in this BVP that the proposed model is able topredict the TRIP strain accumulation during cyclic pseudoelastic loading, and also is capableto capture the decreasing stress level needed to initiate the phase transformations after the25echanically cycling. S t r e ss [ M P a ] Strain
Expriment Simulation
Figure 3: The 1 st and 50 th response of a NiTiCu SMA wire subjected to 50 unaxial tensile stress cycling.The experimental data used for comparison are referenced from Strnadel & Miyazaki (2011). In this section, two BVPs for SMAs with different material compositions, i.e., Ni . Ti . (at.%) and Ni . Ti . Hf (at.%), subjected to uniaxial actuation cycling are analyzed. TheNiTi material is a classical type of SMA having critical phase transforming temperaturesclose to room temperature, while the NiTiHf belongs to the so-called high-temperature SMAcategory that can function under very extreme environments with their phase transforma-tion temperatures around 100 ◦ C. These two BVPs are chosen to check the model’s fidelityover multiple SMA material systems, also to check the proposed capability of capturing theload-free TWSME for thermomechanically trained SMAs. The load-free thermal cycling re-sponse of SMAs are examined and compared to available experimental data to validate suchcapability.In the NiTi case, the SMA dogbone specimen was subjected to a 100 thermal cyclingbetween 310 and 440 K under a constant stress 200 MPa. After the 100 training cycles, thebias load was removed, and the final thermal cycling was performed to check the TWSMEresponse under load-free condition. The material parameters used here to simulate this ex-26
10 20 30 40 50 600.000.020.040.060.08
Experiment Simulation T R I P S t r a i n Number of cycles
Figure 4: Accumulation of TRIP strain for the NiTiCu SMA wire subjected to 50 unaxial tensile stresscycling. The experimental data used for comparison are referenced from Strnadel & Miyazaki (2011). periment are listed in table 3. Similar to the experimental procedure in the case of NiTi, theNiTiHf dogbone specimen was also loaded under 200 MPa and experienced a 100 thermalcycling between 310 and 580 K. The load-free TWSME was also checked after the thermo-mechanical training cycles. The material parameters used to simulate the NiTiHf actuationresponse are listed in table 4. The experimental data used for comparison were initiallyreported in Atli et al. (2015).The simulated cyclic actuation response by the proposed model are summarized in Fig.5 for NiTi SMA and in Fig. 6 for NiTiHf SMA. In the NiTi case, the results for selectedtraining cycle 1, cycle 10, cycle 20, cycle 40, cycle 70, and cycle 100 are plotted and comparedto the experiment results. It is shown that the simulated transformation characteristics inthese selected training cycles, including transformation temperatures, transformation trainmagnitude, and TRIP strain accumulation, are in good agreement with the experiment data.More specifically, the NiTi SMA accumulated about 11% TRIP strain after the 100 thermaltraining cycle. Additionally, the simulated load-free TWSME curve for NiTi SMA alsocorrelates well with the experimental one. In the case of NiTIHf, it can be seen that thehigh-temperature SMA material system shows a quite different actuation response comparedto that of NiTi. Specifically, the phase transformation temperatures are much higher, and27uch less TRIP strain is accumulated in the end. Although those characteristics are quitedifferent compared to NiTi, the simulated results also agree well with the experimental data.The comparison between simulation results and experimental data demonstrates that theproposed model is enabled to capture the load-free TWSME response for thermomechanicallytrained SMAs, and also the fidelity of the proposed modeling capabilities over multiple SMAmaterial systems.
Table 3: Model parameters used for NiTi SMA under uniaxial cyclic actuation loading.
Type Parameter Value Parameter Value E A C A
15 [MPa/K] E M C M ν A = ν M M s
330 [K]12 α A = α M × − [K − ] M f
300 [K] H max i A s
351 [K] H max f A f
375 [K] H min k t τ crit Smooth hardening parameters n n n n σ b
80 [MPa] λ C p × − C p able 4: Model parameters used for NiTiHf SMA under uniaxial cyclic actuation loading. Type Parameter Value Parameter Value E A
70 [GPa] C A E M
70 [GPa] C M ν A = ν M M s
441 [K]12 α A = α M × − [K − ] M f
430 [K] H max i A s
460 [K] H max f A f
466 [K] H min k t τ crit n n n n σ b
35 [MPa] λ C p × − C p
00 350 400 450 S t r a i n [ % ] Temp [K]
Sim_TWSMEExp_TWSME S t r a i n [ % ] Temp [K]
Sim_Cycle 100Exp_Cycle 100 S t r a i n [ % ] Temp [K]
Sim_Cycle 70Exp_Cycle 70 S t r a i n [ % ] Temp [K]
Sim_Cycle 40Exp_Cycle 40 S t r a i n [ % ] Temp [K]
Sim_Cycle 20Exp_Cycle 20 S t r a i n [ % ] Temp [K]
Sim_TWSMEExp_TWSME S t r a i n [ % ] Temp [K]
Sim_Cycle 10Exp_Cycle 10
300 350 400 450 S t r a i n [ % ] Temp [K]
Sim_Cycle 1Exp_Cycle 1 (b)(c)(e) (d)(f)(g) (h)(a)
Figure 5: Selected thermal cycling and load-free TWSME response for the NiTi SMA under a 200 MPaconstant load. (a) Combined, (b) Cycle 1, (c) Cycle 10, (f) Cycle 20, (e) Cycle 40, (f) Cycle 70, (g) Cycle100, (h) TWSME cycle. The experimental data used for comparison are referenced from Atli et al. (2015).
00 400 500 600 S t r a i n [ % ] Temp [K]
Sim_Cycle 1Exp_Cycle 1
300 400 500 600 S t r a i n [ % ] Temp [K]
Sim_Cycle 10Exp_Cycle 10
300 400 500 600 S t r a i n [ % ] Temp [K]
Sim_Cycle 40Exp_Cycle 40
300 400 500 600 S t r a i n [ % ] Temp [K]
Sim_Cycle 70Exp_Cycle 70
300 400 500 600 S t r a i n [ % ] Temp [K]
Sim_Cycle 100Exp_Cycle 100
300 400 500 600 S t r a i n [ % ] Temp [K]
Sim_TWSMEExp_TWSME
300 400 500 600 S t r a i n [ % ] Temp [K]
Sim_Cycle 20Exp_Cycle 20
300 400 500 600 S t r a i n [ % ] Temp [K]
Exp_TWSMESim_TWSME (b)(c)(e) (d)(f)(g) (h)(a)
Figure 6: Selected thermal cycling and load-free TWSME responses for the NiTiHf SMA under a 200 MPaconstant load. (a) Combined, (b) Cycle 1, (c) Cycle 10, (f) Cycle 20, (e) Cycle 40, (f) Cycle 70, (g) Cycle100, (h) TWSME cycle. The experimental data used for comparison are referenced from Atli et al. (2015). = 0.5 mmd = 2 mm a = 100 mmb = 10 mm Figure 7: The geometry of the Ni Ti (at.%) plate actuator with a centric hole.Table 5: Model parameters used for the Ni Ti (at.%) plate actuator. Type Parameter Value Parameter Value E A
70 [GPa] C A
22 [MPa/K] E M
70 [GPa] C M
22 [MPa/K]Key material parameters ν A = ν M M s
318 [K]12 α A = α M × − [K − ] M f
298 [K] H max i A s
332 [K] H max f A f
352 [K] H min k t τ crit n n n n σ b λ N/ATRIP parameters C p × − C p After the analysis of SMAs under uniaxial loading condition, this part analyzes theBVP when SMAs are subjected to multiaxial stress conditions. As it was mentioned in theintroduction, it is imperative to understand the cyclic response of SMAs under the multiaxialloading conditions, as the majority of applications require the functionality of SMA-basedcomponents under multiaxial stress state that are caused by geometry discontinuities, suchas notches and holes coming from installment requirements. It has been reported that TRIPstrain evolves differently in the multiaxial stress state compared to that in the uniaxial one.In order to study the effect of the stress multiaxiality, a plate actuator with a centric holemade from Ni Ti (at.%) subjected to cyclic actuation loading is chosen as the BVP to beinvestigated. The geometry of the plate actuator is shown in Fig. 7, which has 100 mm in the32ength, 10 mm in the width, and 0.5 mm in the thickness. The centric hole has a diameter of 2mm. In the experiment, a nominal load of 136 MPa was applied in the longitudinal directionof the plate, thereafter it was subjected to a thermal cycling between 280 and 400 K for 100cycles while the bias load was maintained constant. The TRIP strain in the longitudinaldirection was recorded at the end of each thermal cycling by using the Digital CorrectionImage (DIC) technique. A simulation was performed by using the presented model whereinthe loading conditions are the modeled the same as the experiment. The material parametersused in this simulation are listed in table 5.As it can be seen from the principal stress contour in Fig. 8(a), the non-uniform multi-axial stress field is caused due to the presence of the hole, a much larger stress field in thedesignated red region while a much smaller stress field in the rest part. As a result of suchstress multiaxiality, the phase transformation of the plate actuator is largely redistributed.As it is shown Fig. 8 (b), the stress concentrated part starts the phase transformation earlierwhile in the less stressed part it is later initiated. Furthermore, refer to Fig. 9 for the TRIPstrain accumulation results, the experimental DIC result is shown in the bottom row whilethe prediction is shown in the upper row of the figure. The experimental results show thatthe multiaxial stress state has a significant effect on the TRIP strain evolution from cycleto cycle. Specifically, the TRIP strain accumulates much faster in the stress concentratedpart compared to the less stress concentrated part. In cycle 10, There is around 0.17% TRIPstrain accumulated in the stress concentration part while the rest of the plate actuator ismuch less than 0.1%. In cycle 100, there is 0.32% TRIP strain generated in the red areawhile the less stressed part accumulates around 0.15%, and the rest part of the plate underuniform stress field accumulates around 0.25%. The observation on the experimental resultsindicated that the generation of TRIP strain is highly stress-dependent. Although this fea-ture seems very intrinsic to model, the proposed model, however, did a quite good job onpredicting that. The simulation results not only well captured the overall distribution shapefor the accumulated TRIP strain, the magnitude in the stress concentrated and less stressedregions are also in good agreement with the experiments. In this BVP, the proposed modelis demonstrated to have the capability to predict the highly stress-dependent TRIP straingeneration when SMAs are subjected to a multiaxial stress state.33 Max. Principal Stress (MPa)
Mart. Vol. Fract. (a) (b)
Figure 8: The simulated stress distribution and martensitic volume fraction contour for the Ni Ti (at.%)plate actuator. (a) Maximum principal stress, (b) martensitic volume fraction. Cycle 100
Cycle 70Cycle 40Cycle 20Cycle 10
TRIP strain h (%) Figure 9: TRIP strain evolution contour for the Ni Ti (at.%) plate actuator subjected to cyclic thermo-mehcanical loading. The first row in this figure is simulation while the second row is experimental DIC result.The experimental data is referenced from Wheeler (2017). . CONCLUSIONS In this work, a three-dimensional finite strain constitutive model for SMAs consideringthe stress-dependent TRIP evolution under multiaixal stress state, as well as the TWSMEat load-free condition is proposed. This model is developed based upon the early work byLagoudas et al. (2012), and largely inspired by its consideration of TRIP (Bo & Lagoudas1999 b , Lagoudas & Entchev 2004), as well as its recent extension into large deformationframework (Xu et al. 2019 a ). By consideration of the martensitic volume fraction, transfor-mation strain, internal stress, and TRIP strain tensors as internal state variables, the modelis enabled to capture the following important characteristics for polycrystalline SMAs. (1)The model formulated based on a finite deformation theory is inherently able to capture thelarge strains and rotations exhibited by SMAs under cyclic thermomechanical loading. (2)Through the introduction of an accumulating internal stress tensor, the model is capable ofpredicting the TWSME for thermomechanically trained SMAs under load-free conditions,and also capturing the decreasing stress level to initiate the phase transformation duringpseudoelastic cycling. (3) By proposing a new TRIP strain evolution law, the model is ableto predict the stress-dependent TRIP strain generation when SMAs are subjected to a mul-tiaxial stress state. A detailed implementation procedure of the proposed model is presentedthrough a user-defined material subroutine within a finite element framework allowing forsolving complex BVPs, and a comprehensive instruction on calibrating the model parametersas well as the derivation of continuum tangent stiffness matrix are also provided. The pro-posed modeling capabilities have been validated through both unaxial and multiaxial BVPsfor a wide range of SMA material systems under different thermomehcnaical loading paths,which demonstrates the presented model’s fidelity as an efficient design and analysis tool forthe future applications of SMA-based multifunctional components. ACKNOWLEDGMENT
The authors would like to acknowledge the financial support provided by the QatarNational Research Fund (QNRF) under the grant number: NPRP 7-032-2-016, and theNational Aeronautics and Space Administration (NASA) through the University LeadershipInitiative (ULI) program under the grant number: NNX17AJ96A. The authors would also35ike to thank Drs. Ibrahim Karaman, Robert Wheeler, Kadri Can Atli, and Mr. GlenBigelow for providing the invaluable experimental data for model validation. At last, theauthors would like to thank Dr. Anargyros Karakalas to proof reading the manuscript andprovide his invaluable comments.
DATA AVAILABILITY
The user-defined material subroutine (in the form of Abaqus UMAT) described in thispaper for the proposed SMA constitutive model considering TRIP and TWSME is availableto access from GitHub: https://github.com/Aero-tomato/SMA-UMAT.
References
Arghavani, J., Auricchio, F. & Naghdabadi, R. (2011), ‘A finite strain kinematic hardeningconstitutive model based on hencky strain: general framework, solution algorithm andapplication to shape memory alloys’,
International Journal of Plasticity (6), 940–961.Atli, K., Karaman, I., Noebe, R., Bigelow, G. & Gaydosh, D. (2015), ‘Work productionusing the two-way shape memory effect in niti and a ni-rich nitihf high-temperature shapememory alloy’, Smart Materials and Structures (12), 125023.Auricchio, F., Reali, A. & Stefanelli, U. (2007), ‘A three-dimensional model describing stress-induced solid phase transformation with permanent inelasticity’, International Journal ofplasticity (2), 207–226.Barrera, N., Biscari, P. & Urbano, M. F. (2014), ‘Macroscopic modeling of functional fatiguein shape memory alloys’, European Journal of Mechanics-A/Solids , 101–109.Bo, Z. & Lagoudas, D. C. (1999 a ), ‘Thermomechanical modeling of polycrystalline smasunder cyclic loading, part i: theoretical derivations’, International Journal of EngineeringScience (9), 1089–1140.Bo, Z. & Lagoudas, D. C. (1999 b ), ‘Thermomechanical modeling of polycrystalline smasunder cyclic loading, part iii: Evolution of plastic strains and two-way shape memoryeffect’, International Journal of Engineering Science (9), 1175–1203.36oyd, J. G. & Lagoudas, D. C. (1996), ‘A thermodynamical constitutive model for shapememory materials. part i. the monolithic shape memory alloy’, International Journal ofPlasticity (6), 805–842.Bruhns, O., Xiao, H. & Meyers, A. (1999), ‘Self-consistent eulerian rate type elasto-plasticity models based upon the logarithmic stress rate’, International Journal of Plastic-ity (5), 479–520.Bruhns, O., Xiao, H. & Meyers, A. (2001 a ), ‘Large simple shear and torsion problems inkinematic hardening elasto-plasticity with logarithmic rate’, International journal of solidsand structures (48), 8701–8722.Bruhns, O., Xiao, H. & Meyers, A. (2001 b ), ‘A self-consistent eulerian rate type model forfinite deformation elastoplasticity with isotropic damage’, International Journal of Solidsand Structures (4), 657–683.Calkins, F. & Mabe, J. (2016), Flight test of a shape memory alloy actuated adaptivetrailing edge flap, in ‘ASME 2016 Conference on Smart Materials, Adaptive Structuresand Intelligent Systems’, American Society of Mechanical Engineers, pp. V001T04A007–V001T04A007.Haghgouyan, B., Hayrettin, C., Baxevanis, T., Karaman, I. & Lagoudas, D. C. (2019), ‘Frac-ture toughness of niti–towards establishing standard test methods for phase transformingmaterials’, Acta Materialia , 226–238.Haghgouyan, B., Shafaghi, N., Aydıner, C. C. & Anlas, G. (2016), ‘Experimental and com-putational investigation of the effect of phase transformation on fracture parameters of ansma’,
Smart Materials and Structures (7), 075010.Hartl, D. J. & Lagoudas, D. C. (2007), ‘Aerospace applications of shape memory alloys’, Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace En-gineering (4), 535–552. 37artl, D. & Lagoudas, D. (2009), ‘Constitutive modeling and structural analysis consider-ing simultaneous phase transformation and plastic yield in shape memory alloys’,
SmartMaterials and Structures (10), 104017.Jani, J. M., Leary, M., Subic, A. & Gibson, M. A. (2014), ‘A review of shape memory alloyresearch, applications and opportunities’, Materials & Design , 1078–1113.Jiang, D. & Landis, C. M. (2016), ‘A constitutive model for isothermal pseudoelasticitycoupled with plasticity’, Shape Memory and Superelasticity (4), 360–370.Karakalas, A. A., Manolas, D. I., Machairas, T. T., Riziotis, V. A. & Saravanos, D. A. (2019),‘Active load alleviation potential of adaptive wind turbine blades using shape memory alloyactuators’, Wind Energy (5), 620–637.Kelly, A., Stebner, A. P. & Bhattacharya, K. (2016), ‘A micromechanics-inspired constitu-tive model for shape-memory alloys that accounts for initiation and saturation of phasetransformation’, Journal of the Mechanics and Physics of Solids , 197–224.Khan, A. S. & Huang, S. (1995), Continuum theory of plasticity , John Wiley & Sons.Lagoudas, D. C. (2008),
Shape memory alloys: modeling and engineering applications ,Springer Science & Business Media.Lagoudas, D. C. & Bo, Z. (1999), ‘Thermomechanical modeling of polycrystalline smas undercyclic loading, part ii: Material characterization and experimental results for a stabletransformation cycle’,
International Journal of Engineering Science (9), 1141–1173.Lagoudas, D. C. & Entchev, P. B. (2004), ‘Modeling of transformation-induced plasticity andits effect on the behavior of porous shape memory alloys. part i: constitutive model forfully dense smas’, Mechanics of Materials (9), 865–892.Lagoudas, D., Hartl, D., Chemisky, Y., Machado, L. & Popov, P. (2012), ‘Constitutive modelfor the numerical analysis of phase transformation in polycrystalline shape memory alloys’, International Journal of Plasticity , 155–183.38eclercq, S. & Lexcellent, C. (1996 a ), ‘A general macroscopic description of the thermome-chanical behavior of shape memory alloys’, Journal of the Mechanics and Physics of Solids (6), 953–980.Leclercq, S. & Lexcellent, C. (1996 b ), ‘A general macroscopic description of the thermome-chanical behavior of shape memory alloys’, Journal of the Mechanics and Physics of Solids (6), 953–980.Lexcellent, C., Leclercq, S., Gabry, B. & Bourbon, G. (2000), ‘The two way shape memoryeffect of shape memory alloys: an experimental study and a phenomenological model’, International Journal of Plasticity (10-11), 1155–1168.Lim, J. T. & McDowell, D. L. (1994), Degradation of an ni-ti alloy during cyclic loading, in ‘Smart Structures and Materials 1994: Smart Materials’, Vol. 2189, International Societyfor Optics and Photonics, pp. 326–341.Mabe, J., Brown, J. & Calkins, F. (2014), Flight test of a shape memory alloy actuated adap-tive trailing edge flap, part 1, in ‘Proceedings of SMST 2014 the International Conferenceon Shape Memory and Superelastic Technologies’.Meyers, A., Xiao, H. & Bruhns, O. (2003), ‘Elastic stress ratcheting and corotational stressrates’, Tech. Mech , 92–102.Meyers, A., Xiao, H. & Bruhns, O. (2006), ‘Choice of objective rate in single parameterhypoelastic deformation cycles’, Computers & structures (17), 1134–1140.Niccoli, F., Garion, C., Maletta, C., Sgambitterra, E., Furgiuele, F. & Chiggiato, P. (2017),‘Beam-pipe coupling in particle accelerators by shape memory alloy rings’, Materials &Design , 603–611.Otsuka, K. & Wayman, C. M. (1999),
Shape memory materials , Cambridge university press.Patoor, E., Eberhardt, A. & Berveiller, M. (1996), ‘Micromechanical modelling of superelas-ticity in shape memory alloys’,
Le Journal de Physique IV (C1), C1–277.39atoor, E., Lagoudas, D. C., Entchev, P. B., Brinson, L. C. & Gao, X. (2006), ‘Shape memoryalloys, part i: General properties and modeling of single crystals’, Mechanics of materials (5), 391–429.Peraza-Hernandez, E., Hartl, D., Galvan, E. & Malak, R. (2013), ‘Design and optimiza-tion of a shape memory alloy-based self-folding sheet’, Journal of Mechanical Design (11), 111007.Qidwai, M. & Lagoudas, D. (2000 a ), ‘Numerical implementation of a shape memory al-loy thermomechanical constitutive model using return mapping algorithms’, InternationalJournal for Numerical Methods in Engineering (6), 1123–1168.Qidwai, M. & Lagoudas, D. (2000 b ), ‘On thermomechanics and transformation surfacesof polycrystalline niti shape memory alloy material’, International Journal of Plasticity (10), 1309–1343.Saint-Sulpice, L., Chirani, S. A. & Calloch, S. (2009), ‘A 3d super-elastic model for shapememory alloys taking into account progressive strain under cyclic loadings’, Mechanics ofmaterials (1), 12–26.Scalet, G., Niccoli, F., Garion, C., Chiggiato, P., Maletta, C. & Auricchio, F. (2019), ‘Athree-dimensional phenomenological model for shape memory alloys including two-wayshape memory effect and plasticity’, Mechanics of Materials p. 103085.Simo, J. C. & Hughes, T. J. (2006),
Computational inelasticity , Vol. 7, Springer Science &Business Media.Song, G., Patil, D., Kocurek, C. & Bartos, J. (2010), Applications of shape memory alloysin offshore oil and gas industry: a review, in ‘Earth and Space 2010: Engineering, Science,Construction, and Operations in Challenging Environments’, pp. 1551–1567.Strnadel, B. & Miyazaki, S. (2011), ‘Modelling residual strains during cycling of ti–ni andti–ni–cu shape memory alloys in a pseudoelastic range of behaviour conditions’, Strain , e457–e466. 40trnadel, B., Ohashi, S., Ohtsuka, H., Ishihara, T. & Miyazaki, S. (1995 a ), ‘Cyclic stress-strain characteristics of ti-ni and ti-ni-cu shape memory alloys’, Materials Science andEngineering: A (1-2), 148–156.Strnadel, B., Ohashi, S., Ohtsuka, H., Ishihara, T. & Miyazaki, S. (1995 b ), ‘Effect of mechan-ical cycling on the pseudoelasticity characteristics of ti-ni and ti-ni-cu alloys’, MaterialsScience and Engineering: A (1-2), 187–196.Tabesh, M., Boyd, J., Atli, K. C., Karaman, I. & Lagoudas, D. (2018), ‘Design, fabrication,and testing of a multiple-actuation shape memory alloy pipe coupler’,
Journal of IntelligentMaterial Systems and Structures (6), 1165–1182.Tanaka, K., Nishimura, F., Hayashi, T., Tobushi, H. & Lexcellent, C. (1995), ‘Phenomeno-logical analysis on subloops and cyclic behavior in shape memory alloys under mechanicaland/or thermal loads’, Mechanics of Materials (4), 281–292.Thamburaja, P. & Anand, L. (2001), ‘Polycrystalline shape-memory materials: effect ofcrystallographic texture’, Journal of the Mechanics and Physics of Solids (4), 709–737.Wang, J., Moumni, Z., Zhang, W. & Zaki, W. (2017), ‘A thermomechanically coupled finitedeformation constitutive model for shape memory alloys based on hencky strain’, Interna-tional Journal of Engineering Science , 51–77.Wang, X., Xu, B. & Yue, Z. (2008), ‘Micromechanical modelling of the effect of plastic defor-mation on the mechanical behaviour in pseudoelastic shape memory alloys’,
InternationalJournal of Plasticity (8), 1307–1332.Wheeler, R., Saunders, R., Pickett, K., Eckert, C., Stroud, H., Fink, T., Gakhar, K., Boyd, J.& Lagoudas, D. (2015), Design of a reconfigurable sma-based solar array deployment mech-anism, in ‘ASME 2015 Conference on Smart Materials, Adaptive Structures and IntelligentSystems’, American Society of Mechanical Engineers, pp. V001T02A010–V001T02A010.Wheeler, R. W. (2017), Actuation Fatigue Characterization Methods and Lifetime Predic-tions of Shape Memory Alloy Actuators, PhD thesis.41heeler, R. W., Hartl, D. J., Chemisky, Y. & Lagoudas, D. C. (2014 a ), Characterizationand modeling of thermo-mechanical fatigue in equiatomic niti actuators, in ‘ASME 2014Conference on Smart Materials, Adaptive Structures and Intelligent Systems’, AmericanSociety of Mechanical Engineers, pp. V002T02A009–V002T02A009.Wheeler, R. W., Hartl, D. J., Chemisky, Y. & Lagoudas, D. C. (2014 b ), Characterizationand modeling of thermo-mechanical fatigue in equiatomic niti actuators, in ‘ASME 2014Conference on Smart Materials, Adaptive Structures and Intelligent Systems’, AmericanSociety of Mechanical Engineers, pp. V002T02A009–V002T02A009.Xiao, H., Bruhns, I. O. & Meyers, I. A. (1997 a ), ‘Logarithmic strain, logarithmic spin andlogarithmic rate’, Acta Mechanica (1-4), 89–105.Xiao, H., Bruhns, O. & Meyers, A. (1997 b ), ‘Hypo-elasticity model based upon the logarith-mic stress rate’, Journal of Elasticity (1), 51–68.Xiao, H., Bruhns, O. & Meyers, A. (2006), ‘Elastoplasticity beyond small deformations’, ActaMechanica (1-2), 31–111.Xu, L., Baxevanis, T. & Lagoudas, D. (2017 a ), A finite strain constitutive model consideringtransformation induced plasticity for shape memory alloys under cyclic loading, in ‘8thECCOMAS Thematic Conference on Smart Structures and Materials’, pp. 1645–1477.Xu, L., Baxevanis, T. & Lagoudas, D. (2018), A three-dimensional constitutive model forpolycrystalline shape memory alloys under large strains combined with large rotations, in ‘ASME 2018 Conference on Smart Materials, Adaptive Structures and Intelligent Systems’,American Society of Mechanical Engineers Digital Collection.Xu, L., Baxevanis, T. & Lagoudas, D. C. (2017 b ), A finite strain constitutive model formartensitic transformation in shape memory alloys based on logarithmic strain, in ‘25thAIAA/AHS Adaptive Structures Conference’, p. 0731.Xu, L., Baxevanis, T. & Lagoudas, D. C. (2019 a ), ‘A three-dimensional constitutive modelfor the martensitic transformation in polycrystalline shape memory alloys under largedeformation’, Smart Materials and Structures (7), 074004.42u, L., Baxevanis, T. & Lagoudas, D. C. (2019 b ), A three-dimensional constitutive modelingfor shape memory alloys considering two-way shape memory effect and transformation-induced plasticity, in ‘AIAA Scitech 2019 Forum’, p. 1195.Yu, C., Kang, G., Kan, Q. & Song, D. (2013), ‘A micromechanical constitutive model basedon crystal plasticity for thermo-mechanical cyclic deformation of niti shape memory alloys’, International Journal of Plasticity , 161–191.Yu, C., Kang, G., Kan, Q. & Zhu, Y. (2015), ‘Rate-dependent cyclic deformation of super-elastic niti shape memory alloy: thermo-mechanical coupled and physical mechanism-basedconstitutive model’, International Journal of Plasticity , 60–90.Zaki, W. & Moumni, Z. (2007), ‘A three-dimensional model of the thermomechanical behaviorof shape memory alloys’, Journal of the Mechanics and Physics of Solids (11), 2455–2490. 43 ppendix A. Continuum tangent stiffness and thermal moduli In order for the displacement-based Finite element solver to attain a converged solutionfor the global equilibrium equations, the so-called tangent stiffness matrix must be defined inthe UMAT. In general, the aforementioned matrix can be derived using either the continuumor the time-discretized equations of stress and transformation function. The tangent stiffnessmatrix derived using the time-discretized equations is commonly referred as consistent stiff-ness matrix, which is the consistent one used in the UMAT framework, and therefore leadsto a quadratic convergence rate in the Newton-Raphson process within the finite elementsolver. The tangent stiffness matrix, derived using the continuum equations, is called thecontinuum stiffness matrix and it may lead to a slower convergence rate in comparison to thetime discretized form, especially when large analysis time steps are used. However, it shouldbe noted that both approaches will lead to the same solution once the Newton-Raphson pro-cess converges. An extensive discussion about this matter can be found in Lagoudas (2008).In the current work, the continuum tangent matrix is the opted one due to the simplicity ofits derivations. Although we acknowledge the aforementioned potential disadvantages, thatcould be addressed in a future work to maximize the algorithm performance.A detailed derivation of the continuum tangent stiffness matrix is provided in the fol-lowing part. For the sake of completeness, the presented derivation has also included thederivation of the tangent thermal moduli. Although it should be noted that the tangentthermal moduli is required to be defined in the UMAT, only if the coupling between thethermal and mechanical equilibrium equations is considered, which has not been consideredin the current work. In general, these continuum tangent tensors can be formulated as shownin equation (A.1), where L is called the tangent stiffness matrix and Θ is the tangent thermalmoduli. Be noted that the Vogit notation of tensors is used in the actual implementation.˚ τ = L ˚ h + Θ ˙ T (A.1)Applying the logarithmic rate on constitutive equation (26) yields,˚ τ = C [˚ h − α ˙ T − (∆ S τ + ∆ α ∆ T + Λ + Λ tp ) ˙ ξ ] (A.2)Taking the chain rule differentiation on the transformation function equation (44) gives,˙Φ = ∂ τ Φ : ˚ τ + ∂ T Φ ˙ T + ∂ ξ Φ ˙ ξ = 0 (A.3)44ubstitute equation (A.2) into equation (A.3) to cancel ˚ τ and solve it for ˙ ξ , the followingexpression for ˙ ξ can be obtained,˙ ξ = − ∂ τ Φ : C ˚ h + ( ∂ T Φ − ∂ τ Φ : C α ) ˙ T∂ ξ Φ − ∂ τ Φ : C (∆ S τ + Λ + Λ tp + ∆ α ∆ T ) (A.4)As the phase difference of thermal expansion coefficients ∆ α is quite small usually, it isneglected in the derivation for the sake of brevity. The final explicit expression in the formas equation (A.1) can be obtained as follows by substitution of equation (A.4) into the rateform constitutive equation (A.2),˚ τ = (cid:104) C + [ C (∆ S τ + Λ + Λ tp )] ⊗ [ C ∂ τ Φ] ∂ ξ Φ − ∂ τ Φ : C (∆ S τ + Λ + Λ tp ) (cid:105) ˚ h + (cid:104) − C α + C (∆ S τ + Λ + Λ tp )( ∂ T Φ − ∂ τ Φ : C α ) ∂ ξ Φ − ∂ τ Φ : C (∆ S τ + Λ + Λ tp ) (cid:105) ˙ T (A.5)in which the continuum tangent stiffness matrix L is, L = C + [ C (∆ S τ + Λ + Λ tp )] ⊗ [ C ∂ τ Φ] ∂ ξ Φ − ∂ τ Φ : C (∆ S τ + Λ + Λ tp ) (A.6)and the continuum thermal moduli Θ is,Θ = C (∆ S τ + Λ + Λ tp )( ∂ T Φ − ∂ τ Φ : C α ) ∂ ξ Φ − ∂ τ Φ : C (∆ S τ + Λ + Λ tp ) − C α (A.7)In order to fully determine the explicit expressions for L and Θ during the implementa-tion, the following terms ∂ τ Φ, ∂ ξ Φ, ∂ T Φ used in above equations can be calculated by usingthe symbolic calculation toolbox provided in MATLAB.
Appendix B. Model calibration
This section presents a detailed instruction on how to identify all the model parametersfrom a set of calibration tests. In general, those parameters can be categorized into threemajor groups, i.e., the key material parameters, smooth hardening parameters, TRIP andinternal stresses parameters. Note that the strain measure used here is logarithmic (or true)strain rather than engineering (or infinitesimal) strain.First, material constants such as the elastic modulus ( E A , E M ) of austenite and marten-site can be obtained by calculating the slopes at martensitic phase and austenite phase from45niaxial mechanical loading-unloading experiment, i.e., the pseudoelastic response as shownin Fig. 10(a). Poisson’s ratios ν A and ν B are usually assumed as 0 .
33 for metallic materials.In order to construct the phase diagram, thermal cycling of the material under constantuniaxial tensile stress, i.e., actuation response as shown in Fig. 10(b), is performed, by whichcritical transformation temperatures ( M τ s , M τ f , A τ s , A τ f ) and transformation strain H cur ( τ ) fora specific applied stress level τ are measured. Such experiments are performed at threedifferent stress levels ( τ , τ , τ ). By using these collected temperature and transformationstrain information, the phase diagram and the H cur curve can be constructed, see Fig. B.11.Therefore, phase diagram related stress influence coefficients ( C A , C M ), phase transformationtemperatures ( M s , M f , A s , A f ) at zero stresses, are determined. In addition, model param-eters related to H cur curve ( H min , H max , k t , τ crit ) are also obtained. Secondly, the smoothhardening related coefficients ( n , n , n , n ) are determined by best matching the smooth-ness of response curve corner at the initiation and completion during phase transformation.The thermal expansion coefficients, usually assumed there is no phase differences α A = α M ,can also be determined through an actuation response.In order to calibrate the TRIP related model parameters ( C p , C p ), the thermal cyclingof SMAs, cyclic actuation response shown as Fig. 12(a), subjected to a constant stress level isneeded, in which the stress magnitude should exceed a saturation point beyond which no moretransformation strain is generated even the stress further increases. From such experiment,the collected TRIP strain with respect to the number of loading cycle are used to calibratethe TRIP parameters. Integrating the rate form TRIP evolution equation (36) can obtainthe algebraic form equation (B.1). Here the ratio term H cur /H max is to incorporate stress-dependency effect on C p when SMAs are subject to multiaxial stress state. h tp = H cur H max C p ln(1 + C p ζ d ) (B.1)Under the condition that the material is fully transformed into oriented martensite phaseat a saturated stress value, H cur /H max = 1 and the accumulation of orientated martensiticvolume fraction is reduced to the total one, i.e., ζ d = ζ . Therefore, equation (B.1) can befurther reduced into equation (B.2). TRIP related model parameters ( C p , C p ) can thus becalibrated by best fitting the TRIP strain evolution curve using the logarithmic function in46quation (B.2) as shown in Fig. 12(b). h tp = C p ln(1 + C p ζ ) (B.2) Strain E A E M S t r e ss (a) S t r a i n ( % ) Temperature (°C)
TRIP Strain T r a n s f o r m a ti on S t r a i n S t r a i n A 𝑠 τ A 𝑓 τ M 𝑠 τ M 𝑓 τ stress: τ Temperature (b)Figure B.10: Experiments utilized for the calibration of model parameters. (a) Uniaxial pseudoelastic ex-periment used for the calibration of elastic modulus of austenite and martensit. (b) Uniaxial actuationexperiments under three different stress levels used for the calibration of phase diagram related model pa-rameters. t r e ss Temperature M f M s A s A f C M C A τ τ τ (a) T r a n s f o r m a t i o n s t r a i n Stress τ τ τ H max H cur ( σ b ) σ b (b)Figure B.11: Calibration of model parameters. (a) Phase diagram. (b) H cur curve. Next, internal stress related model parameters ( σ b , λ ) can be calibrated from the samethermal cycling experiment as shown in Fig.12(a) with additional information on TWSMEresponse. As the internal stress is responsible for the generation of TWSME at load-freecondition, it can be inversely determined by using the TWSME curve. Specifically, themagnitude of transformation strain H cur ( σ b ) for the TWSME curve can be measured asshown in Fig.12(a). Based on the H cur curve from Fig. 11(b), in order to produce thismuch transformation strain, the effective stress value of τ eff can be calculated. Also, it isstraightforward to get τ eff = σ b under load-free condition. Thus the value of σ b can beobtained by using the following derivation, H cur ( σ b ) = H min + ( H max − H min )(1 − e − k t ( σ b − τ crit ) ) (B.3)The explicit calculation of σ b is as follow, σ b = τ crit + 1 k t ln (cid:0) H max − H min H max − H cur ( σ b ) (cid:1) (B.4)The additional model parameter λ is a curve fitting parameter which controls how fastthe quantities, such as the internal stress or H max , evolve into their saturated values. It canbe determined by best fitting the internal stress evolution law (41) or H max degradation law48 rainingTWSME T R I P s t r a i n H c u r ( σ b ) i n T W S M E Temperature S t r a i n (a) Exp
Fit