A thermomechanically coupled finite deformation constitutive model for shape memory alloys based on Hencky strain

International Journal of Engineering Science 117 (2017) 51–77

Contents lists available at ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

A thermomechanically coupled finite deformation constitutive model for shape memory alloys based on Hencky strain Jun Wang a, Ziad Moumni a,b,∗, Weihong Zhang b,∗, Wael Zaki c a

IMSIA, UMR 8193 CNRS-EDF-CEA-ENSTA, Université Paris Saclay, 828 Boulevard des Maréchaux, 91762 Palaiseau Cedex, France Engineering Simulation and Aerospace Computing (ESAC), Northwestern Polytechnical University, 710072 Xi’an, Shaanxi, China c Khalifa University of Science, Technology, and Research P.O. Box 127788, Abu Dhabi, UAE b

a r t i c l e

i n f o

Article history: Received 9 November 2016 Revised 26 April 2017 Accepted 4 May 2017

Keywords: Shape memory alloy Constitutive model Finite deformation Thermomechanical coupling Smooth transition Numerical simulation

a b s t r a c t This paper presents a new thermomechanically coupled constitutive model for polycrystalline shape memory alloys (SMAs) undergoing finite deformation. Three important characteristics of SMA behavior are considered in the development of the model, namely the effect of coexistence between austenite and two martensite variants, the variation of hysteresis size with temperature, and the smooth material response at initiation and completion of phase transformation. The formulation of the model is based on a multi-tier decomposition of the deformation kinematics comprising, a multiplicative decomposition of the deformation gradient into thermal, elastic and transformation parts, an additive decomposition of the Hencky strain into spherical and deviatoric parts, and an additive decomposition of the transformation stretching tensor into phase transformation and martensite reorientation parts. A thermodynamically consistent framework is developed, and a Helmholtz free energy function consisting of elastic, thermal, interaction and constraint components is introduced. Constitutive and heat equations are then derived from this energy in compliance with thermodynamic principles. Time-discrete formulations of the constitutive equations and a Hencky-strain return-mapping integration algorithm are presented. The algorithm is then implemented in Abaqus/Explicit by means of a user-defined material subroutine (VUMAT). Numerical results are validated against experimental data obtained under various thermomechanical loading conditions. The robustness and efficiency of the proposed framework are illustrated by simulating a SMA helical spring actuator. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Shape memory alloys (SMAs) are a class of smart materials capable of inelastic deformation that can be recovered under the influence of appropriate thermomechanical loadings. The most salient features of SMA behavior are undoubtedly their pseudoelasticity (PE) and shape memory effect (SME), which take place at the crystal level by means of first-order diffusionless transformation between austenite and martensite solid phases (Otsuka & Wayman, 1999). These unique features are exploited in a variety of industrial applications, including orthodontic archwires, cardiovascular stents, robotic muscles, ∗

Corresponding authors. E-mail addresses: [email protected] (Z. Moumni), [email protected] (W. Zhang).

http://dx.doi.org/10.1016/j.ijengsci.2017.05.003 0020-7225/© 2017 Elsevier Ltd. All rights reserved.

52

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

shape morphing actuators, vibration dampers, etc. (Auricchio, Reali, & Tardugno, 2010; Gu, Zaki, Morin, Moumni, & Zhang, 2015; Saadat et al., 2002; Williams & Elahinia, 2008; Wu, Qi, Liu, Yang, & Wang, 2007; Yin, He, & Sun, 2014). A review of shape memory alloy research, applications and opportunities is reported by Mohd Jani, Leary, Subic, and Gibson (2014). Rao, Srinivasa, and Reddy (2015) presented an analysis and design methodology for SMA components such as wires, beams, and springs for different applications. Over the last decades, substantial research was dedicated to developing constitutive models that can accurately describe the thermomechanical behavior of SMAs (Auricchio, Bonetti, Scalet, & Ubertini, 2014; Kan & Kang, 2010; Kelly, Stebner, & Bhattacharya, 2015; Lagoudas, Hartl, Chemisky, Machado, & Popov, 2012; Shaw & Kyriakides, 1995; Souza, Mamiya, & Zouain, 1998; Zaki & Moumni, 2007b). A recent review of these models can be found in the work of Cisse, Zaki, and Zineb (2016a, b). The existing models can be classified, depending on the modeling scale, into microscopic, mesoscopic and macroscopic. Microscopic models are developed at the lattice or grain-crystal levels to describe microstructural features such as phase nucleation-completion, interface motion, martensite twin growth, etc. (Benafan et al., 2014; Guthikonda & Elliott, 2013; Kelly et al., 2015). Mesoscopic models are developed on the basis of micromechanics, macroscopic constitutive equations are then derived by the method of scale transition (Khater, Monnet, Terentyev, & Serra, 2014; Levitas & Ozsoy, 2009; Patoor, Lagoudas, Entchev, Brinson, & Gao, 2006). In contrast, macroscopic models are based on simplified description of the material behavior observed at the macroscopic scale and described by means of a finite set of state variables and associated constitutive equations. These models are widely considered the simplest and most convenient for structural analysis and engineering design applications (Arghavani, Auricchio, Naghdabadi, Reali, & Sohrabpour, 2010; León Baldelli, Maurini, & Pham, 2015; Lexcellent, Boubakar, Bouvet, & Calloch, 2006; Moumni, Zaki, & Nguyen, 2008; Peultier, Ben Zineb, & Patoor, 2006; Souza et al., 1998) A large number of phenomenological models have been proposed in the literature with the various degrees of sophistication. Some of these models were shown to accurately describe complex SMA material responses such as asymmetry (Mehrabi, Andani, Elahinia, & Kadkhodaei, 2014; Reedlunn, Churchill, Nelson, Shaw, & Daly, 2014; Rizzoni & Marfia, 2015; Zaki, 2010), multi-variant phase transformation (Auricchio et al., 2014; Idesman, Levitas, Preston, & Cho, 2005; Levitas & Preston, 2002), strain localization (Bechle & Kyriakides, 2016), variation of transformation hysteresis (Lagoudas et al., 2012; Sedlák, Frost, Benešová, Ben Zineb, & Šittner, 2012), cyclic deformation (Kan, Yu, Kang, Li, & Yan, 2016; Kimiecik, Jones, & Daly, 2016; Song, Kang, Kan, Yu, & Zhang, 2014; Yu, Kang, & Kan, 2014; 2015a; Zaki & Moumni, 2007a), crack growth (Baxevanis, Parrinello, & Lagoudas, 2015; Hazar, Zaki, Moumni, & Anlas, 2015), rate dependence (Andani & Elahinia, 2014; Hartl, Chatzigeorgiou, & Lagoudas, 2010; Yu, Kang, Kan, & Zhu, 2015b), etc. Moreover, Doraiswamy, Rao, and Srinivasa (2011) presented a model for simulating the pseudoelastic SMA response including internal hysteresis loops by using a Gibbs potential and a Preisach model. This model was then generalized to predict torsional response of SMA components such as wires, rods and springs at different twists and temperatures (Rao & Srinivasa, 2015), and further to address functional fatigue of these components (Rao & Srinivasa, 2013). However, many of these models disregard thermomechanical coupling (Moumni et al., 2008; Xiao, 2014) and are limited to infinitesimal deformations (León Baldelli et al., 2015; Peng, Pi, & Fan, 2008; Popov & Lagoudas, 2007). The vast majority of SMA models are developed under the assumption of infinitesimal deformations, which greatly simplifies the modeling framework and associated kinematics. However, these models lose accuracy with increasing strains and rotations (Christ & Reese, 2008). Several models using a finite strain formulation (FSF) have been proposed to account for potentially large SMA deformation (Arghavani, Auricchio, & Naghdabadi, 2011; Müller & Bruhns, 2006; Paranjape, Manchiraju, & Anderson, 2016; Reese & Christ, 2008; Teeriaho, 2013; Thamburaja, 2005; Xiao, 2014). These models are generally developed in analogy with finite-strain models for elastoplastic materials, the kinematics of which involve additive split of the stretching tensor D or multiplicative decomposition of the deformation gradient F into elastic and inelastic parts. The former approach is expressed directly in rate form and would be consistent with the latter approach if appropriate integrability conditions are satisfied. Müller and Bruhns (2006) and (Teeriaho, 2013) proposed SMA models for large isotropic deformation with integrable Eulerian rate formulation. Xiao (2014) proposed a finite J2 -flow elastoplastic model with nonlinear combined hardening, in which hysteresis loops may be generated via explicit procedures without reference to phase variables. Thamburaja and Anand (20 01); 20 02) developed a crystal-mechanics-based constitutive model for polycrystalline SMAs to predict the pseudoelastic response under proportional and nonproportional tension-torsion loadings. The model was then improved to account for martensitic reorientation and detwinning (Thamburaja, 2005). Reese and Christ (2008) developed a finite deformation constitutive model to describe the pseudoelasticity of SMAs and proposed a new integration scheme to preserve the incompressibility of material during phase transformation. Arghavani et al. (2011) presented a phenomenological finite-strain formulation featuring kinematic hardening and non-associative flow rules integrated using a solution algorithm involving a nucleation-completion step. Paranjape et al. (2016) developed a finite element framework in which the constitutive formulation accounts for phase transformations at the martensite correspondence variant (CV) scale and for rate-dependent crystal plasticity in austenite. However, these approaches do not always guarantee complete recovery of transformation strain when SMAs undergo multiaxial non-proportional loading, wherein phase transformation and martensite reorientation may take place simultaneously. In addition to the assumption of small strain, most SMA models neglect thermomechanical coupling, which makes them ineffective in solving thermomechanical boundary value problems involving heat generation and transfer within the SMA. Indeed, experimental evidence suggests that heat generation due to latent heat and intrinsic dissipation with a SMA significantly influences the process of phase transformation and therefore the mechanical behavior of the material

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

53

(Müller & Bruhns, 2006; Shaw & Kyriakides, 1995). Shaw and Kyriakides (1995) carried out a series of experiments on SMAs loaded at different strain rates and temperatures. The results show stronger strain rate dependence in air compared to water, indicating the dependence is due to thermomechanical coupling rather than viscosity. Peyroux, Chrysochoos, Licht, and Löbel (1998) used infrared thermography to investigate heat exchange during phase transformation in SMAs and concluded that intrinsic dissipation can be neglected before latent heat. Grabe and Bruhns (2008) conducted biaxial tension/torsion experiments under both isothermal and non-isothermal conditions and found that the stress-strain curves at various strain rates coincide under isothermal conditions, proving that the dependence is due to temperature effect rather than strain rate. Based on these experimental findings, a few models that take into account thermomechanical coupling were proposed (Christ & Reese, 2009; Müller & Bruhns, 2006). More recently, Morin, Moumni, and Zaki (2011) presented a generalized Zaki-Moumni model (Zaki & Moumni, 2007b) for SMAs accounting for thermomechanical coupling. Andani and Elahinia (2014) modified the model originally developed by Lagoudas et al. (2012) to capture the multi-axial pseudoelastic behavior of SMAs under quasi-static isothermal and dynamic loading conditions. However, these models are either limited to small deformations (Andani & Elahinia, 2014; Morin et al., 2011) or developed without considering thermal deformation (Christ & Reese, 2009; Müller & Bruhns, 2006). In the present work, a fully coupled thermomechanical constitutive model for large deformations of SMAs is developed, in which constitutive and heat equations are formulated in terms of Hencky strain. This strain measure, sometimes referred to as logarithmic or nature strain, is used in constitutive modeling of solids due to its remarkable properties in large deformations (Arghavani et al., 2011). Among the usual finite strain measures, the Hencky strain has two important features: (i) only Hencky strain maps the volumetric and isochoric parts of the deformation gradient onto the pure spherical and deviatoric strain measures, which allows straightforward additive split of the total strain (Xiao, Bruhns, & Meyers, 2004) (ii) the logarithmic corotational rate of the Eulerian Hencky strain is the stretching tensor D (Reinhardt & Dubey, 1996), which makes it conjugate to the Cauchy stress in finite strain formulations. The model formulation is based on the following new multi-tier decomposition of deformation kinematics: (i) multiplicative decomposition of deformation gradient F into thermal Fθ , elastic Fe and transformation Ft parts; (ii) additive decomposition of the Hencky strain h into spherical δ and deviatoric h¯ parts; (iii) additive decomposition of the transformation stretching tensor Dt into phase change χ˙ S and martensite reorientation DN parts. In addition, the proposed model incorporates the following three important characteristics of SMA behavior that have not been concurrently addressed in previous work: (i) effect of coexistence between austenite and two martensite variants; (ii) variation of the hysteresis size with temperature; (iii) smooth material response at the initiation and completion of phase transformation. In most phenomenological SMA model, the choice of internal state variables includes a scalar variable representing the volume fraction of martensite and a tensorial variable containing information about the orientation and reorientation of martensite variants (or just one of them) (Lagoudas et al., 2012; Paranjape et al., 2016; Zaki & Moumni, 2007b), wherein the difference between multi-variant and single-variant martensite is neglected. As a result, these models do not account properly for thermomechanical SMA behavior in situations involving coexistence of the three phases. At high stress level, phase transformation takes place only between austenite and single-variant martensite as temperature varies, whereas for low to moderate stress level austenite transforms into both single-variant and multi-variant martensites and vice versa (Wu, Sun, & Wu, 2003). The coexistence of these three phases leads to the unsaturation of transformation strain, which could impact the utilization of SMA, e.g., as actuators. Efforts to model this effect consisted in considering separate volume fractions for self-accommodated and reoriented martensite or relating the magnitude of the transformation strain to the applied stress level (Auricchio et al., 2014; Lagoudas et al., 2012; Zaki & Moumni, 2007b). In the present model, the modeling approach consists in introducing two scalar variables to represent the volume fractions of multi-variant and single-variant martensites. The associated yield functions and evolution rules are derived that govern changes in these variables due to thermomechanical loading. In addition to disregarding the possible coexistence between austenite and two martensite variants, the variation with temperature of the hysteresis size is also rarely addressed by current models despite experimental evidence (Shaw & Kyriakides, 1995). The hysteresis area can be physically interpreted as the amount of energy dissipated due to irreversible processes during phase transformation. In some models, the variation of this area with temperature is accounted for by introducing different temperature influence coefficients for forward (austenite to martensite) and reverse (martensite to austenite) transformation (Lagoudas et al., 2012; Saleeb, Padula, & Kumar, 2011; Sedlák et al., 2012). However, the underlying thermodynamics of such variation have been largely left unaddressed. In this work, the hysteresis size is controlled by the entropy difference ηt between different phases, which is related to latent heat and transformation temperatures. The forward entropy difference ηf and the reverse entropy difference ηr are determined, respectively, in terms of the forward transformation temperature θ AM and the reverse transformation temperature θ MA . Finally, most SMA models predict abrupt transition between elastic and phase transformation regimes. Such transition is experimentally observed in single-crystalline and untrained polycrystalline SMAs (Patoor et al., 2006; Shaw & Kyriakides, 1995), but not in trained polycrystalline SMAs in which phase transformation initiates and completes smoothly and gradually (Grabe & Bruhns, 2008; Sittner, Hara, & Tokuda, 1995). To account for the smooth transition, several transformation hardening functions are proposed in the literature, including exponential (Tanaka, Kobayashi, & Sato, 1986), trigonometric (Liang & Rogers, 1990), power law functions (Lagoudas et al., 2012), etc. In many cases, however, the introduction of these functions result in significant increase in the number of material parameters without achieving the desired accuracy. In the

54

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

present work, a unique tangential transformation hardening function ft featuring only 5 parameters is introduced to describe smooth transformation in polycrystalline SMAs to a good degree of accuracy. The proposed model, incorporating the features described above, is presented in the following sections. In particular, the constitutive equations are derived in Section 2, which includes a description of the proposed Helmholtz free energy density and the thermodynamic framework used for the derivation. Time-discrete formulations of the constitutive equations as well as a Hencky-strain return-mapping integration algorithm are proposed in Section 3. Numerical examples are then carried out and validated against experimental data in Section 4, in which the model is also utilized to simulate the behavior of a SMA helical spring actuator subjected to complex thermomechanical loading and boundary conditions. A conclusion is finally provided in Section 5. 2. Thermomechanically coupled constitutive model based on Hencky strain In this section, on the basis of a multiplicative decomposition of deformation gradient F into elastic Fe , transformation Ft and thermal Fθ parts, a thermodynamically consistent framework and a Helmholtz free energy potential are introduced. Then, a thermomechanically coupled constitutive model in terms of Hencky strains is derived to describe SMA behavior in presence of complex thermomechanical loading and boundary conditions. 2.1. Kinematics In most SMA models (Arghavani et al., 2011; Auricchio et al., 2014; Chatziathanasiou, Chemisky, Chatzigeorgiou, & Meraghni, 2016; Reese & Christ, 2008; Xiao, 2014), the deformation gradient F is simply decomposed into elastic and inelastic parts. However, extensive experimental investigations have shown that phase transformation in SMA is a thermomechanically coupled process, which involves heat production due to latent heat and intrinsic dissipation (Grabe & Bruhns, 2008; Shaw & Kyriakides, 1995). Therefore, an extended multiplicative decomposition of deformation gradient F into elastic Fe , transformation Ft and thermal Fθ parts is assumed here, such that

F = Fe Ft Fθ ,

J = Je Jt Jθ ,

F¯ = F¯e F¯t F¯θ ,

(1)

where Fe is defined with respect to a local unstressed intermediate configuration, Ft is defined with respect to a thermally expanded intermediate configuration and Fθ is defined with respect to the reference configuration. Je = det Fe , Jt = det Ft and Jθ = det Fθ denote the volume-changing deformations, F¯e , F¯t and F¯θ denote the volume-preserving deformations and det F¯e = det F¯t = det F¯θ = 1. Substituting Eq. (1) into Eq. (A.6) gives the following additive split of the velocity gradient L:

L = L¯ +

   1 ˙ 1 δ 1 = L¯ e + F¯e L¯ t + F¯t L¯ θ F¯t−1 F¯e−1 + δ˙ e + δ˙ t + δ˙ θ 1, 3 3

(2)

where

L¯ e = F¯˙ e F¯e−1 ,

δ˙ e =

Je˙ , Je

L¯ t = F¯˙ t F¯t−1 , J˙ Jt

δ˙ t = t , δ˙ θ =

L¯ θ = F¯˙ θ F¯θ−1 , Jθ˙ . Jθ

(3)

Using Eq. (A.7), the elastic, transformation, thermal stretching and spin tensors are defined as

1 ˙ δe 1 and We = skew(Le ), 3 1 Dt = sym(Lt ) = D¯ t + δ˙ t 1 and Wt = skew(Lt ), 3 1 ¯ Dθ = sym(Lθ ) = Dθ + δ˙ θ 1 and Wθ = skew(Lθ ). 3 De = sym(Le ) = D¯ e +

(4)

With regard to the thermal and transformation deformations, the following three kinematic assumptions are made: (i) First, the thermal deformation is assumed to be an isotropic thermal expansion, so that 1

Fθ = Jθ3 1.

(5)

Hence, using Eq. (4)3 ,

Dθ =

1 ˙ δ 1 3 θ

and Wθ = 0.

(6)

(ii) Second, the transformation deformation is assumed incompressible, that is

Jt = det Ft = 1

and Ft = F¯t .

(7)

Hence, using Eq. (3)5 and (4)2 ,

δ˙ t = 0 and Dt = D¯ t .

(8)

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

55

(iii) Last, for isotropic materials, the assumption of the irrotationality of inelastic flow is widely used in modern elastoviscoplastic constitutive theories (Anand, Ames, Srivastava, & Chester, 2009; Dafalias, 1984; Gurtin & Anand, 2005). Therefore, the transformation flow is assumed irrotational (zero transformation spin)

Wt = 0.

(9)

Hence, the transformation stretching tensor

Dt = Lt = F˙t Ft−1

and F˙t = Dt Ft .

(10)

Using the above three kinematic assumptions, the stretching and spin tensors, D and W, are written as





D = D¯ e + sym F¯e Dt F¯e−1 +





 1 ˙ δe + δ˙ θ 1, 3

(11)

W = We + skew F¯e Dt F¯e−1 . Using Eq. (A.4), the elastic, transformation and thermal right and left Cauchy-Green deformation tensors are written as 2

Ce = FeT Fe = Je3 C¯e , Ct = FtT Ft ,

2

be = Fe FeT = Je3 b¯ e ,

bt = Ft FtT ,

(12)

2 3

Cθ = bθ = Jθ 1, where det C¯e = det b¯ e = 1 and det Ct = det bt = 1. Moreover, using Eq. (A.5), the elastic, transformation and thermal Lagrangian and Eulerian Hencky strain tensors are written as

1 ln Ce = H¯ e + 2 1 Ht = ln Ct , ht = 2 1 Hθ = hθ = δθ 1, 3 He =

1 1 1 δe 1, he = ln be = h¯ e + δe 1, 3 2 3 1 ln bt , 2

(13)

where

H¯ e =

1 ln C¯e , 2

1 h¯ e = ln b¯ e , 2

δe = ln Je , δθ = ln Jθ .

(14)

By making use of Eq. (A.22), the logarithmic tensor function maps the unimodular tensors C¯e , b¯ e , Ct and bt onto the traceless tensors H¯ e , h¯ e , Ht and ht , i.e., tr(H¯ e ) = tr(h¯ e ) = tr(Ht ) = tr(ht ) = 0. 2.2. Thermodynamic framework The model is derived based on the principle of virtual power (PVP) in compliance with the first (energy balance) and second (entropy inequality) laws of thermodynamics. A reduced form of the energy balance between internal energy E, internal stress power Pi and heat supply Q is written as



  ˙ ev dV − Jσ : D dV +

B   B E˙

 ∂ B

 Pi

 q · n dS −



B

hv dV = 0,

(15)



Q˙

where ev and hv are, respectively, the internal energy and the heat source per unit volume in spatial configuration, q is the heat flux vector per unit area and n is the outward unit vector normal to the boundary in the spatial configuration. Using Eqs. (A.17) and (11)1 , the internal stress power can be written as

Jσ : D = pδ˙ e + pδ˙ θ + s : D¯ e + F¯eT sF¯e−T : Dt .

(16)

Substituting Eq. (16) into Eq. (15) and localizing the obtained integral equation gives the following local form of the energy balance:





e˙ v − pδ˙ e + pδ˙ θ + s : D¯ e + F¯eT sF¯e−T : Dt − hv + ∇ · q = 0.

(17)

The local form of the Clausius–Duhem entropy inequality is expressed as

1

η˙ v θ − hv + ∇ · q − q · ∇θ ≥ 0, θ

(18)

 = ev − ηv θ and ˙ = e˙ v − η˙ v θ − ηv θ˙ .

(19)

where ηv is the entropy per unit volume in the spatial configuration. The Helmholtz free energy  per unit volume and its time derivative ˙ in the spatial configuration are given by

56

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

Combining the local energy balance in Eq. (17), the local entropy inequality in Eq. (18) and the free energy in Eq. (19) gives the following inequality:

 1    p δ˙ e + δ˙ θ + s : D¯ e + F¯eT sF¯e−T : Dt − ˙ + ηv θ˙ − q · ∇θ ≥ 0.

(20)

θ

The Helmholtz free energy per unit volume is assumed here to depend on the elastic deformation gradient Fe through the elastic Lagrangian Hencky strain He , on the transformation deformation gradient Ft through the transformation Lagrangian Hencky strain Ht , on the thermal deformation gradient Fθ through the thermal Lagrangian Hencky strain Hθ , on scalar internal variables χ M , χ S , and on temperature θ as

   = ψ He , Ht , Hθ , χ M , χ S , θ , χM

(21)

χS

where and denote, respectively, volume fractions of multi-variant and single-variant martensites. The free energy  is assumed to be an isotropic function of He , Ht and Hθ , which allows  to be expressed as the function of their argument invariants. As the Lagrangian and Eulerian Hencky strains share the same invariants, the Helmholtz free energy  can be written in terms of Eulerian Hencky strains as

    ψ He , Ht , Hθ , χ M , χ S , θ = ψ δe , h¯ e , ht , δθ , χ M , χ S , θ .

(22)

Hence, the time derivatives of  gives

˙ =

∂ψ ˙ ∂ψ ¯˙ ∂ψ ˙ ∂ψ ˙ ∂ψ M ∂ψ S ∂ψ ˙ δ + :h + :h + δ + χ˙ + χ˙ + θ. ∂δe e ∂ h¯ e e ∂ ht t ∂δθ θ ∂χ M ∂θ ∂χ S

(23)

Because of isotropy, the derivatives of the isotropic functions ∂ψ and ∂∂ψ are symmetric and coaxial with h¯ e and ht . ht ∂ h¯ e Using the definition of the logarithmic corotational rate of the Eulerian Hencky strain in Eqs. (A.12) and (A.16),

∂ψ ¯˙ ∂ψ ¯˚ L ∂ψ ¯ : he = : he = : De , ∂ h¯ e ∂ h¯ e ∂ h¯ e ∂ψ ˙ ∂ψ ˚ L ∂ψ :h = :h = : Dt , ∂ ht t ∂ ht t ∂ ht

(24)

where h¯˚ Le and h˚ tL are, respectively, the logarithmic corotational rates of the elastic and transformation Eulerian Hencky strains. Substituting Eqs. (23) and (24) into inequality (20) gives









∂ψ ˙ ∂ψ ˙ ∂ψ ∂ψ T ¯ −T ¯ ¯ p− δ + p− δ + s− : De + Fe sFe − : Dt − ∂δe e ∂δθ θ ∂ ht ∂ h¯ e

1 ∂ψ M ∂ψ S ∂ψ χ˙ − χ˙ − + ηv θ˙ − q · ∇θ ≥ 0. ∂θ θ ∂χ M ∂χ S

(25)

For arbitrary thermodynamic processes, inequality (25) is guaranteed by the following choice of the constitutive equations with respect to p, s and ηv :

p=

∂ψ ∂ψ ∂ψ ∂ψ = , s= , ηv = − , ¯ ∂δe ∂δθ ∂θ ∂ he

(26)

and by the following inequalities, which ensure nonnegative intrinsic dissipation in arbitrary evolutions of the internal state variables Dt , χ M and χ S , so that



F¯eT sF¯e−T −

∂ψ ∂ ht



: Dt −

∂ψ M ∂ψ S χ˙ − χ˙ ≥ 0, ∂χ M ∂χ S

(27)

as well as nonnegative dissipation due to heat conduction, so that

1 − q · ∇θ ≥ 0.

(28)

θ

Using the definition of mandel stress M in Eq. (A.19), the following deviatoric mandel stress tensor is obtained with respect to the intermediate unstressed configuration:

M¯ e = F¯eT sF¯e−T = C¯e S¯e ,

(29) F¯e−1 sF¯e−T

where C¯e = F¯eT F¯e is the isochoric elastic right Cauchy-Green tensor and S¯e = is a new second Piola-Kirchhoff stress tensor defined with respect to the intermediate unstressed configuration. As s and V¯ e commute, M¯ e can be written as

M¯ e = ReT V¯ eT sV¯ e−1 Re−T = ReT sRe ,

(30)

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

57

which is a symmetric and traceless tensor by rotating the deviatoric Kirchhoff stress back into the intermediate unstressed configuration. Substituting M¯ e into inequality (27) gives



M¯ e −

∂ψ ∂ ht



: Dt −

∂ψ M ∂ψ S χ˙ − χ˙ ≥ 0. ∂χ M ∂χ S

(31)

Substituting Eq. (26) into Eq. (23) gives the following first Gibbs relation of the free energy  and the temperature θ as

  ∂ψ ˙ ∂ψ M ∂ψ S ˙ + ηv θ˙ = p δ˙ e + δ˙ θ + s : h¯˙ e + :h + χ˙ + χ˙ , ∂ ht t ∂χ M ∂χ S

(32)

  e˙ v − θ η˙ v = p δ˙ e + δ˙ θ + s : h¯˙ e +

(33)

with Eq. (19)2 , which yields the following second Gibbs relation of internal energy e and the entropy ηv :

∂ψ ˙ ∂ψ M ∂ψ S :h + χ˙ + χ˙ . ∂ ht t ∂χ M ∂χ S

Combining Eqs. (33) and (17) gives the following entropy relation:



∂ψ ∂ψ M ∂ψ S θ η˙ v = hv − ∇ · q + M¯ e − : Dt − χ˙ − χ˙ . ∂ ht ∂χ M ∂χ S

(34)

The internal energy and the entropy densities are assumed to depend on the same state variables as the Helmholtz free energy:

 δe , h¯ e , ht , δθ , χ M , χ S , θ ,   ηv = η δe , h¯ e , ht , δθ , χ M , χ S , θ . ev = e



(35)

The specific heat capacity is defined as

  ∂ e δe , h¯ e , ht , δθ , χ M , χ S , θ cv = , ∂θ def

(36)

with Eqs. (19) and (26)3 , which can be written as

cv = θ

    ∂η δe , h¯ e , ht , δθ , χ M , χ S , θ ∂ 2 ψ δe , h¯ e , ht , δθ , χ M , χ S , θ = −θ . ∂θ ∂θ 2

(37)

Therefore, from Eqs. (35), (26)3 and (37),



∂η ¯˙ ∂η ˙ ∂η ˙ ∂η M ∂η S ∂η ˙ θ η˙ v = θ δ + :h + :h + δ + χ˙ + χ˙ + cv θ˙ ∂δe e ∂ h¯ e e ∂ ht t ∂δθ θ ∂χ M ∂χ S 2

∂ 2 ψ ¯˙ ∂ 2ψ ˙ ∂ 2ψ ˙ ∂ 2ψ M ∂ 2ψ S ∂ ψ ˙ = −θ δ + :h + :h + δ + χ˙ + χ˙ + cv θ˙ . ∂ θ ∂ δe e ∂ θ ∂ h¯ e e ∂ θ ∂ ht t ∂ θ ∂ δθ θ ∂ θ ∂ χ M ∂θ∂χS

(38)

Substituting Eq. (38) into Eq. (34) gives the following temperature evolution equation:



cv θ˙ = hv − ∇ · q +

M¯ e −



∂ψ ∂ ht



∂ψ M ∂ψ S χ˙ − χ˙ + ∂χ M ∂χ S  

: Dt −

D=intrinsic dissipation

∂ ψ ∂ 2ψ ˙ ∂ ψ ˙ θ δe + : D¯ e + δθ + ∂ θ ∂ δe ∂ θ ∂ δθ ∂ θ ∂ h¯ e   

2

2

X =thermoelastic coupling

∂ 2ψ M ∂ 2ψ S ∂ 2ψ θ : Dt + χ ˙ + χ ˙ . ∂ θ ∂ ht ∂θ∂χM ∂θ∂χS   

(39)

Z=thermo-transformation coupling

At this stage of the development of the theory and the experimental investigation, the thermoelastic coupling terms in Eq. (39) that give rise to a temperature variation due to variations in the state variables δ e , δ θ and h¯ e are not well characterized (Ames, Srivastava, Chester, & Anand, 2009; Bouvard et al., 2013). Hence, as approximations, the thermoelastic coupling is assumed to be neglected and instead a factor ω = (1 + X /D ) is introduced to represent the ratio of intrinsic dissipation converted into heat. Consequently, Eq. (39) reduces to

cv θ˙ = hv − ∇ · q + ωD + Z,

(40)

where the thermo-transformation coupling term Z is associated with the latent heat of phase transformation. The heat flux q is taken to be governed by Fourier’s law

q = −k∇θ , where k is the nonnegative thermal conductivity.

(41)

58

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

2.3. Helmholtz free energy The Helmholtz free energy function  in Eq. (21) is additively split as follows:

 =  e +  θ +  int +  cst ,

(42)

where  e is the elastic energy,  θ is the thermal contribution to free energy (Lagoudas et al., 2012; Morin et al., 2011),  int is the interaction energy between single-variant martensite, multi-variant martensite and austenite phases (Moumni et al., 2008),  cst is the potential energy due to the internal constraints. By applying the rule of mixtures (Auricchio et al., 2014; Lagoudas et al., 2006; Moumni et al., 2008), the elastic energy  e and the thermal energy  θ in Eq. (42) can be expressed in the following form:

 e    = 1 − χ M − χ S  e,A + χ M  e,M + χ S  e,S ,  θ = 1 − χ M − χ S  θ ,A + χ M  θ ,M + χ S  θ ,S ,

(43)

where  ∗, A ,  ∗, M and  ∗, S are the free energies of austenite, multi-variant and single-variant martensites, respectively. Considering a Reuss rheological model (Auricchio et al., 2014; Moumni et al., 2008; Reuss, 1929), the elastic free energy of the mixture in Eq. (43)1 is given by

1 2

 e = K δ 2 + μh¯ e 2 − 3α K δ (θ − θ0 ),

(44)

where δ = δe + δθ is the volumetric part of the Hencky strain, θ 0 is the reference temperature, α is the thermal expansion coefficient. The bulk modulus K is assumed to be equal for all phases and shear modulus μ are determined from the Reuss model (Auricchio et al., 2014; Wagner & Windl, 2008) using the relations:



K = KA = K M = K S , 1 S M 1 M 1 S 1 μ = 1−χ −χ μA + χ μM + χ μS .

(45)

The thermal energy  θ in Eq. (42), sometimes referred to as chemical energy (Auricchio et al., 2014; Sedlák et al., 2012), is related to changes in the internal energy and entropy. The thermal energies in Eq. (43)2 for the three different phases,  θ , ∗ , are defined as



 θ ,∗ = e∗0 − η0∗ θ + cv∗



θ (θ − θ0 ) − θ ln θ0



,

(46)

where ∗ denotes either of the phases A, M and S, e0 and η0 are reference internal energy and entropy at reference temperature θ 0 , cv is the specific heat capacity. Substituting Eq. (46) into Eq. (43)2 gives the following thermal energy  θ :



  θ  θ =eA0 − η0A θ + cvA (θ − θ0 ) − θ ln − χ M eAM + χ S eAS + θ0

 M AM   M AM  θ S AS S AS , χ η + χ η θ − χ cv + χ cv (θ − θ0 ) − θ ln θ0

(47)

where eA∗ = eA0 − e∗0 , ηA∗ = η0A − η0∗ and cvA∗ = cvA − cv∗ are the differences of internal energy, entropy and heat capacity between austenite and multi-/single-variant martensites. Regarding the three different material phases, the following assumptions are introduced to simplify the formulation of the model : (i) since many DSC experiments showed that the heat capacities of austenite and martensite are almost equal (Qidwai & Lagoudas, 20 0 0), cvA = cvM = cvS = cv and cvAM = cvAS = 0; (ii) equal internal energy and entropy between multi- and single-variant martensites, i.e. eAM = eAS = Lh and ηAM = ηAS = ηt = Lh /θt , where Lh is the latent heat of phase transformation between austenite and martensite, the transformation temperature θ t and the entropy difference ηt are given by



θt =

θAM θMA

A → M, S, and ηt = M, S → A,



η f = Lh /θAM , ηr = Lh /θMA ,

(48)

where θ AM and θ MA are forward and reverse phase transformation temperatures, ηf and ηr are forward and reverse entropy differences. With these assumptions, the thermal energy in Eq. (47) reduces to



  θ  θ = eA0 − η0A θ + cv (θ − θ0 ) − θ ln + ηt χ M + χ S (θ − θ0 ). θ0

(49)

The interaction energy  int in Eq. (42) represents the crystal boundary energy between austenite, multi- and singlemartensite phases, which is often derived from micro-mechanical or crystallographic considerations and finally expressed as

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

59

a function of internal state variables (Auricchio et al., 2014; Zaki & Moumni, 2007b). In the present work, the interaction energy  int is assumed to depend on the multi-variant χ M , single-variant χ S martensite volume fractions and transformation strain ht as

  1  int =  int χ M , χ S , ht = gt (χ ) + μt ht 2 ,

(50)

2

where χ = χ M + χ S is the volume fraction of total martensite. The first term gt on the right-hand side of Eq. (50) represents the energy of interaction between austenite and martensite, while the second term 12 μt ht 2 represents increase in interaction due to orientation/reorientation of martensite plates (multi-variant martensite to single-variant martensite). The difference between the two martensite variants is neglected with respect to interaction with austenite. The internal state variables χ M , χ S and ht used for modeling obey the following physical constraints: • The volume fractions of multi-variant χ M , single-variant χ S and total χ M + χ S martensites are restricted to the [0, 1] interval,

0 ≤ χ M , χ S , χ M + χ S ≤ 1,

(51)

which can be rewritten as the following three active constraints:

  χ M ≥ 0, χ S ≥ 0, 1 − χ M + χ S ≥ 0.

(52)

• The magnitude of transformation strain due to orientation of martensite variants must remain lower than a materialspecific limit h0 , i.e.

h0 − ht  ≥ 0.

(53)

The above physical constraints on the internal state variables are guaranteed by the following Lagrangian potential as the constraints contribution to free energy  cst in Eq. (42):

   cst = −ζ M χ M − ζ S χ S − ζ MS 1 − χ M − χ S − ζ t (h0 − ht  ), ζ M,

ζ S,

ζ MS

(54)

ζt

where and are nonnegative Lagrangian multipliers. The constraint energy (54) is complemented by the following classical Kuhn–Tucker conditions:

⎧ M χ ≥ 0, ⎪ ⎨ S χ ≥ 0,  M S ⎪ ⎩1 − χ + χ ≥ 0, h0 − ht  ≥ 0,

ζ M ≥ 0, ζ S ≥ 0, ζ MS ≥ 0, ζ t ≥ 0,

ζ M χ M = 0; ζ S χS = 0;  ζ MS 1 − χ M − χ S = 0; ζ t (h0 − ht  ) = 0.

In conclusion, the Helmholtz free energy  is constructed as



1  = K δ 2 + μh¯ e 2 − 3α K δ (θ − θ0 ) + eA0 − η0A θ + cv (θ − θ0 ) − θ ln 2

  1 ηt χ M + χ S (θ − θ0 ) + gt + μt ht 2 − 2  M M    ζ χ + ζ S χ S + ζ MS 1 − χ M − χ S + ζ t (h0 − ht  ) .

(55)



θ θ0

+ (56)

2.4. Constitutive equations Substituting the Helmholtz free energy function  (56) into the thermodynamic framework in Section 2.2, the hydrostatic stress p, deviatoric stress s and entropy ηv are derived as

⎧ ⎨p =

∂  = K δ − 3α ( θ − θ ) , [ 0 ] ∂δ ∂ ¯ e, s=  = 2 μ h ∂ h¯ e     ⎩ ηv = − ∂∂θ = 3α K δ + η0A + cv ln θθ0 − ηt χ M + χ S .

(57)

Using Eqs. (57)2 and (30), the mandel stress M¯ e in the intermediate unstressed configuration is written as

M¯ e = 2μH¯ e . Moreover, the partial derivatives of  with respect to the internal state variables ht ,

⎧ ∂ h ⎨ ∂ ht = μt ht + ζ t htt  , ∂  = η (θ − θ ) + f t − ζ M + ζ MS , t 0 M ⎩ ∂χ ∂  = η (θ − θ ) + f t − ζ S + ζ MS , t 0 S ∂χ

(58)

χM

and

χS

are given by

(59)

where f t = ∂ gt /∂χ is the phase transformation hardening function. In order to capture the smooth phase transformation response, as experimentally demonstrated (Hartl et al., 2010; Lagoudas et al., 2006) and theoretically formulated

60

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

(Auricchio et al., 2014; Lagoudas et al., 2012) in the literature for polycrystalline SMAs, a tangential hardening function ft is considered as

f t = κ tanm

1 2

 (ξa χ + ξb )π + r,

(60)

where κ , m, ξ a , ξ b and r are material parameters, which obey the following constraints:

κ > 0 , 0 < m < 1 , ξ a > 0 , ξb > 0 , ξa + ξb < 1 .

(61)

χ S,

The single-variant martensite volume fraction according to the literature (Lagoudas et al., 2012; Moumni et al., 2008; Reese & Christ, 2008), can be expressed in terms of the transformation strain as

χS =

ht  h0

.

(62)

In addition, the transformation strain ht can be equivalently expressed as

ht = ht N t ,

(63)

where N t is the direction of transformation strain. Using Eqs. (62) and (63), the time derivative of ht is given by

h˙ t = h0 N t χ˙ S + h0 χ S N˙ t ,

(64)

where the first term represents the phase transformation between austenite and single-variant martensite, the second term represents the reorientation of single-variant martensite. Substituting Eq. (64) into Eq. (A.12) gives the following split of logarithmic corotational rate of the transformation Hencky strain:





h˚ tL = h0 N t χ˙ S + h0 χ S N˙ t − tL N t + N t tL .

(65)

Considering the coaxiality of ht and N t , the following corotational rate of N t is introduced as L

DN = N˚ t = N˙ t − tL N t + N t tL .

(66)

Therefore, the transformation stretching tensor Dt is rewritten as

Dt = h0 χ˙ S N t + h0 χ S DN .

(67)

Substituting Eq. (67) into inequality (31) gives



h0 χ



S

M¯ e −



∂ψ ∂ ht





∂ψ M ∂ψ ∂ψ ¯ : DN − χ˙ + h0 Me − : Nt − χ˙ S ≥ 0. ∂ ht ∂χ M ∂χ S        AM

AN

(68)

AS

where the thermodynamic forces AN associated to DN , AM associated to χ˙ M and AS associated to χ˙ S are given using Eqs. (58) and (59) as

⎧   ⎨AN = h0 χ S 2μH¯ e − μt ht − ζ t N t , A = −ηt (θ − θ ) − f t + ζ M − ζ MS , ⎩AM = 2μh H¯ : N 0− μ h2 χ S − η (θ − θ ) − f t − h ζ t + ζ S − ζ MS . t t 0 t S 0 e 0 0

(69)

The nonnegative dissipation of multi-variant and single-variant martensite phase transformations is guaranteed by the following choice of the evolution equations:

χ˙ M = γ˙ M

AM

and

|AM |

χ˙ S = γ˙ S

AS

|AS |

,

(70)

where γ M and γ S are nonnegative multipliers. The associated yield functions for phase transformations are given by

FM = |AM | − YM

and FS = |AS | − YS ,

where YM and YS control the thresholds of the phase transformations. Kuhn–Tucker consistency conditions:

 γ˙ M ≥ 0, γ˙ S ≥ 0,

FM ≤ 0, FS ≤ 0,

(71)

χ˙ M ,

FM ,

χ˙ S ,

and FS are subject to the following

γ˙ M FM = 0; γ˙ S FS = 0.

(72)

The nonnegative dissipation of martensite reorientation is guaranteed by the following associative evolution equation:

DN = γ˙ t

AN

 AN 

,

(73)

where γ t is a nonnegative multiplier. The associated yield function for martensite reorientation is given by

Fr = AN  − (χ S )2Yr ,

(74)

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

61

where Yr defines the martensite reorientation threshold. γ t and Fr are subject to the Kuhn–Tucker conditions:

γ˙ t ≥ 0, Fr ≤ 0, γ˙ t Fr = 0.

(75)

Moreover, from the time derivative of N t  = 1 and the Eq. (A.16), one gets

N t : N˙ t = N t : DN = 0.

(76)

Substituting Eq. (73) into Eq. (76) allows writing the Lagrangian multiplier ζ t in Eq. (69) as

ζ t = M¯ e : N t − μt h0 χ S .

(77)

Then, the thermodynamic force AN is written as









AN = h0 χ S M¯ e − M¯ e : N t N t = h0 χ S (I − N t  N t )M¯ e ,

(78)

where I is the fourth-order identity tensor. (I − N t  N t )M¯ e represents the component normal to N t of the mandel stress M¯ e . Substituting Eqs. (70)2 and (73) into Eq. (67) gives the following evolution equation for the transformation strain:

Dt = h0



γS



AS AN N + γ tχS . |AS | t  AN 

(79)

Finally, using Eqs. (39), (40), (70) and (73), the temperature evolution equation is written as

cv θ˙ = hv − ∇ · q + ω

   t  A A γ AN  + γ M |AM | + γ S |AS | + θ ηt γ M M + γ S S , |AM | |AS |

(80)

3. Numerical implementation In this section, the numerical implementation of the constitutive model in an finite element analysis (FEA) framework is carried out. Time-discrete formulation of the proposed model is presented by using a exponential integrator for the transformation evolution equation. Then, a Hencky-strain return-mapping algorithm is introduced, which greatly increases the numerical simplicity and efficiency. 3.1. A time-discrete framework The deformation-driven initial boundary value problem of the proposed thermomechanically coupled finite-strain model is described as: within the time interval [tn , tn+1 ], given the deformation gradient Fn , the internal state variables N tn , χnM , χnS , the temperature θ n at time tn and the deformation gradient F at time tn+1 , determine the state variables and temperature change at time tn+1 satisfying the time-discrete constitutive equations. The evolution Eq. (70)1 of the multi-variant martensite phase fraction χ M is written using a backward-Euler integration scheme as

χ M = χnM + γ M S (AM ),

(81)

γ M

where is the multiplier at time tn+1 , the sign function S (· ) is used to extract the sign of AM . The numerical approximation of the evolution Eq. (79) of the transformation strain is based on an exponential integration scheme as

F t = exp (T )Fnt , with



T = h 0 γ S

(82)



AS AN N + γ t χ S , |AS | t  AN 

(83)

where T is the flow vector of the transformation strain, γ S and γ t are the multipliers at time tn+1 . Regarding the temperature change in Eq. (80), the following time-discrete temperature evolution equation is given as

1 (hv − ∇ · q + ωD + Z ), cv

(84)

M S D = γ t  A  N  M+ γ |AM | + Sγ |AS |, Z = θn ηt γ S (AM ) + γ S (AS ) .

(85)

θ = θn + with



62

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

The time-discrete evolution Eqs. (81), (82) and (84), together with the yield functions (71) and (74), construct the following system of time-discrete constitutive equations:

⎧ M R = χ M − χnM − γ M S (AM ) = 0, ⎪ ⎪ t t t ⎪ ⎪ ⎨Rθ = F − exp (1T )Fn = 0, R = θ − θn − c v ( h v − ∇ · q + ω D + Z ) = 0 , FM = |AM | − YM = 0, ⎪ ⎪ ⎪ ⎪ ⎩FS = S|AS | − YS = 0S, Fr /χ = AN /χ − χ SYr = 0,

(86)

which is complemented by the following time-discrete Kuhn–Tucker conditions:

 M γ ≥ 0 , γ S ≥ 0 , γ t ≥ 0 ,

γ M F M = 0 ; γ S F S = 0 ; γ t F r = 0 .

FM ≤ 0, FS ≤ 0, Fr ≤ 0,

(87)

3.2. A Hencky-strain return-mapping algorithm By using the multiplicative decomposition Ft = F¯e−1 F¯ and Eq. (A.23), the Eq. (82) is equivalently expressed as

F¯ e = F¯ e

trial

exp (−T ),

(88)

where F¯ e trial = F¯ F¯ne is the elastic trial deformation gradient, F¯ is the incremental deformation gradient. Then, the postmultiplication of each side of F¯ e exp(T ) = F¯ e trial by its transpose, together with the use of Eq. (A.24), gives







T V¯ e exp 2Re T Re V¯ e = V¯ e



trial 2

.

(89)

In terms of the Eulerian Hencky strains, h¯ e and h¯ e

 e

exp h¯



e

exp 2R T R

 eT

 e

exp h¯



trial ,

¯ e trial

= exp 2h



Eq. (89) can be written as

.

 e

From the series representation of the tensor exponential (Eq. (A.20)), the term exp h¯ as

 

(90) in above equation can be expressed

 

exp h¯ e = 1 + h¯ e + o h¯ e ,

 e

where o h¯



(91)

is a second-order term in the elastic Hencky strain. Therefore, Eq. (90) is written as

exp 2h¯ e

trial





= exp 2Re T Re

T











 

T T + h¯ e exp 2Re T Re + exp 2Re T Re h¯ e + o h¯ e ,

(92) ht

with the following definition of Green-Lagrange-type tensor of incremental transformation strain  :

ht =

   1 T exp 2Re T Re − 1 , 2

(93)

which can be rewritten equivalently as



exp 2h¯ e

trial







 

= 1 + 2 ht + h¯ e + h¯ e ht + ht h¯ e + o h¯ e .

(94)

By making use of Eq. (A.21), the logarithm of both sides of Eq. (94) gives

2h¯ e

trial

=

∞    k+1 (−1 )k   t 2 h + h¯ e + h¯ e ht + ht h¯ e + o h¯ e . k+1

(95)

k=0

According to the multinomial theorem, Eq. (95) can be approximately expressed as

2h¯ e

trial

=

∞ ∞      (−1 )k  t k+1  (−1 )k  ¯ e k+1 2h + 2h + o ht + o h¯ e , k+1 k+1 k=0

(96)

k=0

where o(ht ) is the second-order term in the incremental transformation strain. In view of Eqs. (91) and (93), Eq. (96) then is written as

h¯ e

trial





 

T = h¯ e + Re T Re + o ht + o h¯ e .

(97)

Since the exponential map integrator (89) is first-order accurate and the elastic strain is assumed to be infinitesimal, the second-order terms can be neglected in Eq. (97) leading to the following approximate formula:

Rt = h¯ e

trial

T

− h¯ e − Re T Re = 0,

(98)

which is used as a substitute for Eq. (86)2 to improve numerical efficiency. Finally, the system (86) is solved using a NewtonRaphson method.

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

63

Table 1 Material parameters used in isothermal pseudoelastic uniaxial tensile simulations. Parameter K

μMS ηf κ ξa r

θ0

YS

Value 23,400 9038 0.147 5.6 0.998 −6.65 240 2.1

Unit

Parameter

Value

Unit

MPa MPa MPa · K−1 MPa – MPa K MPa

μ μt ηr

12,500 500 0.133 0.2 0.001 0.035 2.1 110

MPa MPa MPa · K−1 – – – MPa MPa

A

m

ξb h0 YM Yr

Table 2 Material parameters used in unstressed thermally-induced phase transformation simulation. Parameter K

μMS ηf κ ξa r

θ0 YS

Value 41,667 11,538 0.588 9.6 0.998 10.2 313 8.1

Unit

Parameter

Value

Unit

MPa MPa MPa · K−1 MPa – MPa K MPa

μ μt ηr

26,923 300 0.408 0.2 0.001 0.06 8.1 110

MPa MPa MPa · K−1 – – – MPa MPa

A

m

ξb h0 YM Yr

4. Numerical simulations and experimental validations The proposed model and the corresponding numerical algorithm presented in the previous sections are implemented into the finite element analysis software Abaqus/Explicit by means of a user-defined material subroutine VUMAT. To demonstrate capabilities of the model, numerical simulations are carried out and compared to experimental data under various thermomechanical loading conditions. Finite element simulation of a SMA helical spring actuator undergoing large deformations and temperature variation is also performed. The model parameters used in the simulations are calibrated from the experimental data. A detail calibration procedure is shown in Appendix B, wherein the three-dimensional constitutive equations are simplified to a one-dimensional case, and then a MATLAB routine is used to determine the parameter values. 4.1. Isothermal pseudoelasticity The first set of simulations is dedicated to pseudoelasticity under isothermal uniaxial tensile loading at temperatures of 298 K, 303 K and 313 K. Experimental data for a polycrystalline NiTi wire (50.8 at.% Ni, provided by Memry Corporation) reported by Lagoudas et al. (2012) is utilized for validation. The material parameters used for this set of simulations are calibrated from the experimental data and listed in Table 1. In the simulations, the applied stress is increased from 0 to a maximum value of 700 MPa then removed, with the temperature maintained at a constant value. Fig. 1 shows comparisons between the model predictions and the experimental data reported by Lagoudas et al. (2012). The pseudoelastic stress-strain responses at constant temperatures of 298 K, 303 K and 313 K are, respectively, plotted in Fig. 1(a)–(c). In addition, the evolution of the single-variant martensite volume fraction χ S is presented in Fig. 1(d). From the figures, it is seen that the proposed model captures the pseudoelastic behavior of the considered polycrystalline SMA with good accuracy. In particular, complete shape recovery is achieved upon unloading, a smooth response is observed at the initiation and completion of phase transformation, and the elastic modulus is found to depend on phase composition. 4.2. Thermally-induced phase transformation This section is devoted to the investigation of SMA response to temperature variation under zero-stress condition. Extensive experimental investigations have shown that phase transformation of SMA is a thermomechanically coupled process. The latent heat and intrinsic dissipation contribute to heat generation and absorption during phase transformation, which can be measured by Differential Scanning Calorimetry (DSC) thermoanalytical technique. The difference in the amount of heat required to increase the temperature of a SMA sample compared to a reference specimen is measured as a function of temperature in DSC analysis. The thermally-induced phase transformation and the corresponding heat generation/absorption are predicted by the proposed constitutive model. Experimental date reported by Popov and Lagoudas (2007) is utilized for model validation. The material parameters used in the simulation are determined from the uniaxial stress-strain response and the DSC data, listed in Table 2. In the simulation, an unconstrained SMA component is first uniformly heated from 273 K to 393 K and then uniformly cooled back to 273 K.

64

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

Fig. 1. Model prediction vs. eperimental data under isothermal uniaxial tensile loading: (a) θ = 298 K, (b) θ = 303 K, (c) θ = 313 K, (d) evolution of the single-variant martensite volume fraction.

Fig. 2. Unstressed thermally-induced phase transformation: (a) heat generation and absorption (latent heat), (b) evolution of multi-variant martensite volume fraction.

Fig. 2(a) shows a comparison between model predictions and the experimental DSC data reported by Popov and Lagoudas (2007), wherein the relative heat flow is plotted as a function of temperature. Heat absorption from martensite to austenite phase transformation results in a ridge along the heating part of the curve, while heat generation from austenite to martensite phase transformation results in a valley along the cooling part. The bottom and peak points on the cooling and heating parts, corresponding to temperatures θMA = 303 K and θAM = 339 K, denote average temperatures for forward and reverse phase transformation respectively. It is worth noting that the ridge is short and wide, while the valley is deep and narrow. This is because the area under the curve in each case, representing phase transformation latent heat, is assumed to be the same. In this particular case, the temperature range of reverse transformation is larger than that of forward transformation, resulting in a low and wide ridge versus a deep and narrow valley. Overall, the heat absorption/generation during thermally-induced phase transformation is well described by the proposed model despite a small deviation at completion of

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

65

Table 3 Material parameters used in the simulations of isobaric thermal cycling tests. Parameter K

μMS ηf κ ξa r

θ0

YS

Value 50,0 0 0 23,077 0.2585 7.52 0.998 -2.58 297 5.64

Unit

Parameter

Value

Unit

MPa MPa MPa · K−1 MPa – MPa K MPa

μ μt ηr

23,077 300 0.282 0.08 0.001 0.047 5.64 110

MPa MPa MPa · K−1 – – – MPa MPa

A

m

ξb

h0 YM Yr

Fig. 3. Comparisons between model predictions and the experimental data reported by Wu et al. (2003): transformation strain vs. temperature at constant stress levels 160 MPa and 360 MPa.

the reverse phase transformation. Thanks to the tangential hardening function ft defined in Eq. (60), the smooth experimental DSC curves are simulated with good accuracy. In addition to simulating the heat exchange during phase transformation, the evolution of multi-variant martensite volume fraction is simulated and shown in Fig. 2(b). According to the simulation, the reverse phase transformation of martensite to austenite is found to initiate at temperature 305 K and complete at temperature 353 K, while the forward phase transformation of austenite to martensite initiates at temperature 313 K and completes at temperature 300 K. 4.3. Isobaric thermal cycling tests In isobaric thermal cycling tests, SMA samples are preloaded uniaxially to different constant stress levels, which are maintained while the material undergoes thermally-induced transformation cycles (actuation cycles). The evolution of transformations train with temperature and its dependence on the level of applied stress are predicted by the proposed model. The experimental data reported by Wu et al. (2003) for equiatomic NiTi (50.0 at.% Ni) is utilized to validate the model. The material parameters used in the simulations are calibrated from the experimental data and listed in Table 3. Two constant stress levels of 160 MPa and 360 MPa are considered in the simulations, wherein temperature is decreased uniformly from 393 K to 253 K and then increased uniformly back to 393 K. Fig. 3 shows comparisons between model predictions and the experimental data reported by Wu et al. (2003) for stress levels of 160 MPa and 360 MPa. Since the elastic SMA response was not addressed in the source work, the comparisons shown in Fig. 3 is based on the variation of the transformation strain with temperature. At 160 MPa, the yield functions associated with both single-variant martensite transformation and multi-variant martensite transformation are activated. Thus, austenite transforms simultaneously into these two martensite variants, leading to unsaturated transformation strain of 0.021. In contrast, the stress at 360 MPa is too high to activate the yield function associated with multi-variant martensite transformation. In this case, austenite transforms entirely into single-variant martensite and the maximum transformation strain reaches the saturated value h0 = 0.047, shown as Table 3. It is worth noting that the material response shows abrupt transitions at the completion of forward transformation and at the initiation of reverse transformation at 160 MPa. This phenomenon is due to an assumption in the proposed model that the activation of the yield function associated with multivariant martensite transformation will shut down the transformation to single-variant martensite. This deficiency can be addressed by considering the coupling effect between these two types of martensite transformation. Overall, the evolution of transformation strain with temperature and the variation of the maximum transformation strain with the applied stress are successfully predicted by the proposed model.

66

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77 Table 4 Material parameters used in butterfly-shaped non-proportional simulation. Parameter K

μMS ηf κ ξa r

θ0

YS

Value 52,083 24,038 0.198 6 0.98 −4.8 288 2.4

Unit

Parameter

Value

Unit

MPa MPa MPa · K−1 MPa – MPa K MPa

μ μt ηr

24,038 100 0.198 0.2 0.01 0.03 2.4 110

MPa MPa MPa · K−1 – – – MPa MPa

A

m

ξb h0 YM Yr

Fig. 4. Comparisons between simulations using the present model and experimental data reported by Grabe and Bruhns (2009): (a) butterfly-shaped axial-shear loading path, (b) shear stress vs. axial stress, (c) axial stress-strain response, (d) shear stress-strain response.

4.4. Multiaxial non-proportional loading tests In this section, a biaxial strain-controlled butterfly-shaped loading test is investigated to demonstrate the reliability of the model in presence of multiaxial non-proportional loading conditions. Experimental data reported by Grabe and Bruhns (2009) is utilized as references. The material parameters used in √ √ the simulation are calibrated from the experimental data and listed in Table 4. In the source work, γ = γ / 3 and τ = 3τ are used as shear strain and stress measures for simplicity of a von Mises-type equivalence. Fig. 4 shows the butterfly-shaped strain input and the stress output in both axial and shear directions. The maximum √ axial and shear strain magnitudes reached in each case are ε = γ / 3 = 0.015, as shown in Fig. 4(a). The results of numerical simulation are compared to the reference experimental data in Fig. 4(b–d). Specifically, Fig. 4(b) shows the √ √ scaled shear stress 3τ versus axial stress σ , Fig. 4 (c) shows axial stress σ versus axial strain ε and Fig. 4 (d) shows 3τ versus the √ scaled shear strain γ / 3. The overall agreement with the experimental data is satisfactory despite small deviation in presence of dominant shear and compression loadings. This is because the yield criteria for phase transformation and martensite reorientation of the von Mises type, which does not account for the experimentally observed asymmetry in NiTi behavior in tension, compression and shear. A better simultaneous fit to tension, compression and shear data is achievable by means of more sophisticated criteria.

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

67

Table 5 Material parameters used in non-isothermal simulations at different loading rates. Parameter K

μMS ηf κ ξa r

θ0

YS k

ω

Value 52,083 24,038 0.33 3.92 0.98 1.96 288 7.84 18 0.86

Unit

Parameter

Value

Unit

MPa MPa MPa · K−1 MPa – MPa K MPa W · m−1 · K−1 –

μ μt ηr

24,038 100 0.33 0.2 0.01 0.049 7.84 110 440 50

MPa MPa MPa · K−1 – – – MPa MPa J · kg−1 · K−1 W · m−2 · K−1

A

m

ξb

h0 YM Yr cv h

Fig. 5. Comparisons between model predictions and the experimental data reported by Shaw and Kyriakides (1995) at different loading rates: (a) stressstrain response, (b) temperature variation.

4.5. Pseudoelasticity under non-isothermal conditions To demonstrate the influence of thermomechanical coupling on SMA behavior, simulations of pseudoelastic uniaxial tensile tests under non-isothermal boundary conditions are carried out. Experimental data reported by Shaw and Kyriakides (1995) is used to validate the model. In the referenced work, pseudoelastic uniaxial tensile tests were performed on polycrystalline NiTi specimens (50.1 at.% Ni) at three different strain rates of 0.0 0 04 s−1 , 0.0 04 s−1 and 0.04 s−1 . In each case, heat exchange with the surrounding takes place by means of natural air convection with the external temperature maintained at 343 K. The material parameters used in the simulations are calibrated from the experimental data and listed in Table 5. In the simulations, the applied stress is increased from zero to the maximum of 10 0 0 MPa then decreased to zero for each strain rate. The finite element model is defined with a uniform initial temperature field of 343 K and natural convection at the boundary with a convection heat transfer coefficient h = 50 W m−2 K−1 . Fig. 5 compares the model predictions with the experimental data reported by Shaw and Kyriakides (1995) at three different strain rates of 0.0 0 04 s−1 , 0.0 04 s−1 and 0.04 s−1 . In particular, Fig. 5(a) illustrates the pseudoelastic stress-strain response and Fig. 5(b) the variation of temperature with the strain. The pseudoelastic stress-strain response significantly depends on temperature and therefore is conjointly determined by the energy released or absorbed due to latent heat and intrinsic dissipation as well as the heat transfer by convection into surrounding medium. At low strain rate (0.0 0 04 s−1 ), the rate of heat transfer by convection into the surrounding medium is sufficient to evacuate dissipative and latent phase transformation heat resulting in minor variation in temperature. This variation increases at moderate strain rate (0.004 s−1 ) for which the released heat is not convected into the surrounding medium timely, resulting in increased hysteresis. In this particular case, the finish temperature for forward phase transformation increased by 7 K and for reverse phase transformation decreased by 12 K with respect to the initial temperature of 343 K. Finally, at high strain rate (0.04 s−1 ), the material behavior approaches approximately adiabatic conditions, with temperature increasing to 368 K upon loading then decreasing to 340 K after unloading. The rate of the released energy becomes dominant in this case, resulting in increased temperature and therefore higher thermodynamic stability of the austenite phase. As a result, it becomes increasingly difficult for forward transformation to proceed whereas the reverse transformation becomes easier. A manifestation of this effect is the increased slope of the stress-strain plateaus corresponding to forward and reverse transformations.

68

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

Fig. 6. Simulation of a SMA helical spring actuator: (a) geometry and mesh, (b) force and temperature loading paths.

Fig. 7. SMA helical spring actuator at maximum loading: (a) von Mises stress contour plot, (b) martensite volume fraction contour plot.

4.6. Simulation of a SMA helical spring actuator To demonstrate the usefulness of the proposed model in analyzing complex structural SMA components subjected to non-trivial thermomechanical loading conditions, a SMA helical spring actuator is simulated in this section. The simulation is carried out considering separate loading cases to illustrate stress-induced pseudoelasticity and temperature-induced actuation. Previous work on the simulation of SMA actuators can be found in literature (Auricchio et al., 2014; Hartl et al., 2010; Lagoudas et al., 2012; Saleeb, Dhakal, Hosseini, & Padula II, 2013; Toi, Lee, & Taya, 2004). The SMA helical spring actuator adopted in the simulation has a coil diameter of 20 mm, a wire diameter of 2 mm and a spring pitch of 10 mm. Fig. 6(a) shows the initial geometry, meshed using 3312 coupled temperature-displacement reduced integration hexahedral elements (Abaqus/Explicit C3D8RT). The element highlighted in red is chosen for studying the local material response, i.e. stress-strain behavior and phase transformation. The mesh of the cross-section adjacent to the red element is shown in the figure. Loading is applied at the nodes on end surfaces of the spring. The thermomechanical loading cases considered in the simulations are shown in Fig. 6(b). In the first simulation, the applied forces are increased from zero to the maximum value of 15 N at the time of 5 s and then removed at the time of 10 s, with a prescribed temperature 380 K (black lines in Fig. 6(b)). It is worth noting that the force is applied in a smooth manner to eliminate dynamic effect for quasi-static analysis in Abaqus/Explicit. Fig. 7 shows the contour plots of the von Mises stress and the martensite volume fraction χ in the spring at the maximum force magnitude of 15 N. The maximum von Mises stress is achieved in central part of the spring, shown as Fig. 7(a). Regarding its distribution through the cross-section, the von Mises stress is greater towards the outside compared to the core. In Fig. 7(b), SDV2 is the second state variable used in VUMAT subroutine to represent the volume fraction of martensite. The volume fraction is highest on the inner surface of the central segment of the spring. The local stress response and the phase transformation in the chosen element (red in Fig. 6(a)) are shown in Fig. 8(a). With increasing applied force, the von Mises stress on the element increases as well. Phase transformation from austenite to martensite then takes place when the von Mises stress exceeds 670 MPa. When the load is removed, martensite transforms back to austenite. The global

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

69

Fig. 8. Stress-induced pseudoelasticity of the SMA helical spring actuator: (a) martensite volume fraction vs. von Mises stress curve for the selected element, (b) force vs. displacement curve at the point of application of the load on the end surface.

Fig. 9. SMA helical spring actuator at temperature 280 K: (a) von Mises stress contour plot, (b) martensite volume fraction contour plot.

force-displacement curve at the point of application of the load on the end surface is plotted in Fig. 8(b). The figure shows complete recovery of inelastic strain upon removal of the load. In the second simulation, the spring is first preloaded with a force of magnitude 4 N in a smooth manner over a time period of 5 s at a constant temperature of 380 K. It is then cooled uniformly to 280 K before being uniformly heated back to 380 K over another time period of 5 s (blue lines in Fig. 6(b)). Fig. 9 shows contour plots of the von Mises stress and the martensite volume fraction χ in the spring at 280 K. The von Mises stress is maximum near the inner surface of the central portion of the spring as shown in Fig. 9(a). In the same portion, the martensite volume fraction reaches its maximum value of magnitude 1 along toward the outside of the section, indicating complete transformation of austenite to martensite. The local stress response and phase transformation in the highlighted element (red in Fig. 6(a)) are shown in Fig. 10(a). In the first half of the loading sequence, the von Mises stress increases with increased applied force to nearly 480 MPa. However, as the temperature is decreased from 380 K to 280 K at constant applied force, the von Mises stress decreases to 230 MPa and austenite completely transforms to martensite. When temperature is increased again to 380 K, the martensite finally transforms back to austenite. The displacement of the point of application of the load on the end surface is plotted versus temperature in Fig. 10(b). At 380 K, the applied force of 4 N generates approximately 4 mm of elastic displacement. when temperature is decreased to 320 K, the formation of single-variant martensite results in significant transformation strain. The spring retrieves its original shape as the temperature is increased back to 380 K. The complete actuation cycle of the SMA helical spring actuator is shown in the figure. 5. Conclusion A thermomechanically coupled, finite-strain constitutive model for shape memory alloys has been proposed in this work. The formulation of the model is based on a multi-tier decomposition of deformation kinematics comprising: (i) a multiplica-

70

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

Fig. 10. Temperature-induced actuation of the SMA helical spring actuator: (a) martensite volume fraction vs. von Mises stress curve for the selected element, (b) displacement vs. temperature curve at the point of application of the load on the end surface.

tive decomposition of the deformation gradient into thermal, elastic and transformation parts, (ii) an additive decomposition of the Hencky strain into spherical and deviatoric parts, (iii) an additive decomposition of the transformation stretching tensor into phase transformation and martensite reorientation parts. A thermodynamically consistent framework and a Helmholtz free energy function including elastic, thermal, interaction and constraint components were constructed, from which constitutive and heat equations were derived. Three important characteristics of SMA behavior are considered in formulating the model. First, the effect of phase coexistence between austenite and two martensite variants is accounted for by deriving separate evolution laws of volume fractions for single-variant χ S and multi-variant χ M martensites, which can be active simultaneously. Then, the variation of the hysteresis size with temperature is characterized by means of the forward entropy difference ηf and the reverse difference ηr . Finally, the smooth material behavior at initiation and completion of phase transformation is accounted for using a unique tangential transformation hardening function ft . Time-discrete formulation of the constitutive equations were derived and integrated using a Hencky-strain return-mapping algorithm. Whilst the proposed three-dimensional model was reduced to a one-dimensional version to perform the calibration of the model parameters. The model was then implemented into Abaqus/Explicit using a user-defined material subroutine VUMAT. Numerical simulations are carried out and validated against experimental data for various loading cases, including proportional, non-proportional, isothermal, non-isothermal conditions. Finally, a simulation example is proposed in which a SMA helical spring is first subjected to mechanical loading resulting in a pseudoelastic response then to thermal loading resulting in a temperature-induced actuation. The simulation results illustrate the capability of the proposed model to describe the behavior of SMA devices subjected to complex thermomechanical loadings. Acknowledgments Prof Ziad Moumni would like to thank the 10 0 0 Talent Plan and Northwestern Polytechnical University for their financial support. Prof Weihong Zhang would like to thank the National Natural Science Foundation of China (11432011,11620101002,51521061). Jun Wang would like to thank the financial support from CSC (China Scholarship Council). Appendix A. Mathmatical and mechanical preliminaries A1. Kinematics Let X be an arbitrary material point of a homogeneous body B. The motion of B is described with a smooth one-to-one mapping x = χ(X, t ). The deformation gradient F, velocity v, and velocity gradient L of the motion χ are given by

F = ∇χ,

v = χ˙ , L = ∇x v = F˙ F −1 ,

(A.1)

where the notation ∇ denotes the gradient of a tensor, the notation ∇ x denotes the gradient of a tensor with respect to x, the dot is used to indicate derivative with respect to time. def

By making use of the definition J = det F , F is multiplicatively decomposed into volumetric and isochoric components: 1

F = J 3 F¯ ,

(A.2)

1 3

where J and F¯ denote, respectively, the volume-changing and volume-preserving deformations. Using polar decomposition, F is written as

F = RU = V R,

(A.3)

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

71

where the proper orthogonal tensor R is the local rotation tensor and the symmetric positive definite tensors U and V are right and left stretch tensors, respectively. The right and left Cauchy-Green deformation tensors, C and b, are defined by 2

2

2

C = F T F = U 2 = J 3 U¯ 2 = J 3 F¯ T F¯ = J 3 C¯ , 2 2 2 b = F F T = V 2 = J 3 V¯ 2 = J 3 F¯ F¯ T = J 3 b¯ ,

(A.4)

where the isochoric components U¯ , C¯ , V¯ and b¯ represent volume-preserving deformations, i.e., deformations for which det C¯ = det U¯ = det b¯ = det V¯ = 1. The Lagrangian and Eulerian Hencky strain tensors are respectively defined as

1 1 ln C = H¯ + δ 1, 2 3 1 1 h = ln V = ln b = h¯ + δ 1, 2 3 H = ln U =

(A.5)

where H¯ = ln U¯ and h¯ = ln V¯ are the deviatoric components and δ = ln J the volumetric component of the Hencky strain. Substituting Eq. (A.2) into Eq. (A.1)3 gives the following isochoric and volumetric split of L:

L = L¯ +

1 ˙ δ 1, 3

(A.6)

˙ where L¯ = F¯˙ F¯ −1 and δ˙ = JJ = tr(L ) denote, respectively, the isochoric and volumetric components of L. The stretching and spin tensors, D and W, are defined as the symmetric and skew parts of L. They are given by

1 1 (L + LT ) = D¯ + δ˙ 1, 2 3 1 T ¯ W = skew(L ) = (L − L ) = W , 2 D = sym(L ) =

(A.7)

¯ = skew(L¯ ) are the isochoric components of the stretching and spin tensors. where D¯ = sym(L¯ ) and W Time derivative of (A.4)1 yields the following important relations:

C˙ = 2F T DF

and C¯˙ = 2F¯ T D¯ F¯ .

(A.8)

A2. Corotational rates The corotational rate of a symmetric Eulerian tensor A is defined by





A˚ = Q Q T ˙AQ Q T = A˙ − A + A,

(A.9)

where Q is a time dependent orthogonal tensor denoting the rotation of the frame, A˙ is the material time derivative and  = Q˙ Q T is an antisymmetric tensor denoting the frame-spin. Xiao, Bruhns, and Meyers (1998) showed that the spin tensor  in Eq. (A.9) can be written in the following general form:



n 

=W +

h

i, j=1,i= j

λi λi I

,

I

Pi DPj ,

(A.10)

where λi and Pi are the distinct eigenvalues and eigenprojections and I = trb the first invariant of the left Cauchy-Green deformation tensor b, I = trb, h is a continuous antisymmetric function denoting the particular corotational objective rate. Moreover, for the Eulerian Hencky strain measure h, a logarithmic spin tensor is considered that has the following expression:

L = W +

n  i, j=1,i= j



λi + λ j 2 + P DP , λi − λ j ln λi − ln λ j i j

(A.11)

which yields the following relation between the logarithmic corotational rate of Eulerian Hencky strain h˚ L and the stretching tensor D:

h˚ L = h˙ − L h + hL = D.

(A.12)

Khan and Huang (1995) proposed that the corotational integration for any spatial tensor A can, in general, be written as

 C

Adt = Q





m

Q T AQ dt Q T ,

(A.13)

72

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

where ∫C and ∫m denote, respectively, the corotational integration and the material integration. When the spatial tensor A itself is a corotational rate, i.e., the logarithmic corotational rate h˚ L , substituting Eq. (A.12) into Eq. (A.13) gives



h=

LC

h˚ L dt =



LC

Ddt = Q L

 

QL

T

DQ L dt



QL

T

m

,

(A.14)

where ∫LC denotes the logarithmic corotational integration. The integration scheme (A.14) indicates that the logarithmic corotational integration of the stretching tensor D is the Eulerian Hencky strain h. By using the property of the skew symmetric tensor  = −T , the inner product of A˚ and an arbitrary second-order symmetric tensor B is written as

A˚ : B = A˙ : B + 2 : AB.

(A.15)

If A and B are coaxial, AB is also symmetric, indicating that the inner product of the symmetric tensor AB and the skew tensor  becomes null. Therefore, the following identity is obtained:

A˚ : B = A˙ : B,

(A.16)

when the symmetric tensors A and B are coaxial for any arbitrary choice of the spin tensor . A3. Stress measures The Kirchhoff stress tensor τ and its additive decomposition are given by

τ = Jσ = s + p1,

(A.17)

where σ is the Cauchy stress tensor, s and p = tr(τ )/3 are the deviatoric and hydrostatic parts of the Kirchhoff stress. Hence, the second Piola–Kirchhoff stress tensor is given by





S = F −1 τ F −T = J − 3 pC¯ −1 + F¯ −1 sF¯ −T . 2

(A.18)

The mandel stress tensor is given by

M = F T τ F −T = CS = C¯ F¯ −1 sF¯ −T + p1.

(A.19)

A4. Tensor exponential and logarithm For a generic tensor X, the tensor exponential function can be explicitly expressed by means of its series expansion as

exp (X ) =

∞  1 n X . n!

(A.20)

n=0

Moreover, the series representation of the tensor logarithm function is written as

ln (1 + X ) =

∞  (−1 )n n+1 X . n+1

(A.21)

n=0

The tensor exponential function possesses the following important properties: (i) The determinant of the tensor exponential satisfies

det [exp(X )] = exp [tr(X )].

(A.22)

(ii) The inverse of the tensor exponential satisfies

[exp(X )]−1 = exp(−X ).

(A.23)

(iii) For any invertible tensor B, the tensor exponential satisfies





exp BX B−1 = B exp (X )B−1 .

(A.24)

Appendix B. Calibration of the model parameters This appendix is dedicated to the calibration procedure for estimating values of the model parameters. To this end, the proposed three-dimensional (3D) model is first reduced to a one-dimensional (1D) version. Then, the 1D model is implemented into MATLAB using an explicit integration scheme. Finally, the model parameters are calibrated from the experimental data. The 1D Helmholtz free energy takes the following simplified form:



=

1 σ2 − 3α Kh(θ − θ0 ) + eA0 − η0A θ + cv (θ − θ0 ) − θ ln 2 E



θ θ0



+

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

73

Fig. B1. Calibration of the model parameters from the experimental stress-strain curves: (a) identification of EA , EM , ηf , ηr and YS from the pseudoelastic stress-strain curve, (b) identification of Yr , μt and h0 from the stress-strain curve of martensite orientation.

  1 ηt χ M + χ S (θ − θ0 ) + gt + μt ht2 − 2  M M    ζ χ + ζ S χ S + ζ MS 1 − χ M − χ S + ζ t (h0 − ht ) ,

(B.1)

where σ is the uniaxial stress and E is the elastic modulus defined in terms of the elastic moduli of austenite EA , multivariant martensite EM and single-variant martensite ES as follows:

1 1  1 1 = 1 − χS − χM A + χM M + χS S . E E E E

(B.2)

The thermodynamic forces associated to the multi-variant χ M and single-variant χ S martensites are given by



AM = −ηt (θ − θ0 ) − f t + ζ M − ζ MS , AS = h0 σ − h20 χ S − ηt (θ − θ0 ) − f t − h0 ζ t + ζ S − ζ MS ,

(B.3)

74

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

Fig. B2. Effects of the model parameters κ , m, ξ a , ξ b and r in the tangential hardening function on the stress-strain behavior, arrows indicate the value increase of the respective parameters.

where the entropy difference ηt takes different quantities ηf for the forward and ηr the reverse phase transformation, the tangential hardening function ft is given by

f t = κ tanm

1 2

 (ξa χ + ξb )π + r.

(B.4)

Then, the 1D evolution laws for multi-variant and single-variant martensite phase transformation reduce to

χ˙ M = γ˙ M S (AM ) and χ˙ S = γ˙ S S (AS ),

(B.5)

where the sign function S (· ) is used to extract the signs of AM and AS . The associated yield functions are expressed as



FM = | − ηt (θ − θ0 ) − f t + ζ M − ζ MS | − YM = 0, FS = |h0 σ − h20 χ S − ηt (θ − θ0 ) − f t − h0 ζ t + ζ S − ζ MS | − YS = 0.

(B.6)

The thermodynamic force AN associated with martensite reorientation in the 3D model is given by

AN = h0 χ S (I − N t  N t )M¯ e ,

(B.7)

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

75

where (I − N t  N t )M¯ e represents the component normal to N t of the Mandel stress M¯ e . In the uniaxial stress state, N t and M¯ e are coaxial, indicating null value of this component, and therefore zero of the thermodynamic force:

AN = 0,

(B.8)

in 1D case, which indicates that no martensite reorientation takes place in the uniaxial stress state. The 1D model is then implemented into MATLAB to calibrate the model parameters from the uniaxial experimental data. In the following, the model parameters are determined from the uniaxial experimental data with illustration, whilst physical interpretation of the model parameters are provided. • The elastic moduli EA and EM are determined from the experimental stress-strain curve of pseudoelasticity, as shown in Fig. B.11(a); • The forward ηf and reverse ηr entropy differences govern the variation of the hysteresis size with temperature. This influence is illustrated in Fig. B.11(a) with two pseudoelastic stress-strain curves obtained at temperature θ (black) and θ + θ (blue). The temperature increase of θ results in stress increases on forward σ f and reverse σ r transformation plateaus. Using Eq. (B.6), these two model parameters are calculated as follows:

η f = h 0

σ f θ

and

η r = h 0

σ r . θ

(B.9)

• YM and YS controls the thresholds of the multi-variant and single-variant martensite transformations. For the sake of calculus simplicity, it is assumed YM = YS . These two yield thresholds determine the height of the hysteresis loop σ h on the pseudoelastic stress-strain curve at reference temperature, as shown in Fig. B.11(a). They are calculated by

YM = YS =

1 h 0 σh . 2

(B.10)

• The martensite reorientation threshold Yr , the kinematic hardening modulus μt and the magnitude limit of transformation strain h0 are determined from the stress-strain curve of martensite orientation experiment, wherein the initial multi-variant martensite is completely oriented to the single-variant martensite under a uniaxial tensile loading. Yr takes the initial yield stress, μt the slope of orientation plateau and h0 the magnitude of the residual strain, as shown in Fig. B.11(b). • The model parameters κ , m, ξ a , ξ b and r in the tangential hardening function ft control the smooth material behavior at initiation and completion of phase transformation. Fig. B.12 shows how the parameters κ , m, ξ a , ξ b and r affect the stress-strain behavior of polycrystalline SMAs. The parameters for the three-dimensional model that were determined by following the calibration procedure described in this appendix are listed in Tables 1–5. References Ames, N. M., Srivastava, V., Chester, S. A., & Anand, L. (2009). A thermo-mechanically coupled theory for large deformations of amorphous polymers. part II: Applications. International Journal of Plasticity, 25, 1495–1539. doi:10.1016/j.ijplas.20 08.11.0 05. Anand, L., Ames, N. M., Srivastava, V., & Chester, S. A. (2009). A thermo-mechanically coupled theory for large deformations of amorphous polymers. part i: Formulation. International Journal of Plasticity, 25, 1474–1494. doi:10.1016/j.ijplas.20 08.11.0 04. Andani, M. T., & Elahinia, M. (2014). A rate dependent tension-torsion constitutive model for superelastic nitinol under non-proportional loading; a departure from von mises equivalency. Smart Materials and Structures, 23, 15012. doi:10.1088/0964-1726/23/1/015012. Arghavani, J., Auricchio, F., & Naghdabadi, R. (2011). A finite strain kinematic hardening constitutive model based on hencky strain: General framework, solution algorithm and application to shape memory alloys. International Journal of Plasticity, 27, 940–961. doi:10.1016/j.ijplas.2010.10.006. Arghavani, J., Auricchio, F., Naghdabadi, R., Reali, a., & Sohrabpour, S. (2010). A 3-d phenomenological constitutive model for shape memory alloys under multiaxial loadings. International Journal of Plasticity, 26, 976–991. doi:10.1016/j.ijplas.20 09.12.0 03. Auricchio, F., Bonetti, E., Scalet, G., & Ubertini, F. (2014). Theoretical and numerical modeling of shape memory alloys accounting for multiple phase transformations and martensite reorientation. International Journal of Plasticity, 59, 30–54. doi:10.1016/j.ijplas.2014.03.008. Auricchio, F., Reali, A., & Tardugno, A. (2010). Shape-memory alloys: Effective 3d modelling, computational aspects and design of devices. International Journal of Computational Materials Science and Surface Engineering, 3, 199–223. doi:10.1504/IJCMSSE.2010.033154. Baxevanis, T., Parrinello, A., & Lagoudas, D. (2015). On the driving force for crack growth during thermal actuation of shape memory alloys. Journal of the Mechanics and Physics of Solids, 89, 255–271. doi:10.1016/j.jmps.2015.12.011. Bechle, N. J., & Kyriakides, S. (2016). Evolution of localization in pseudoelastic NiTi tubes under biaxial stress states. International Journal of Plasticity, 82, 1–31. doi:10.1016/j.ijplas.2016.01.017. Benafan, O., Noebe, R. D., Padula II, S. A., Brown, D. W., Vogel, S., & Vaidyanathan, R. (2014). Thermomechanical cycling of a NiTi shape memory alloymacroscopic response and microstructural evolution. International Journal of Plasticity, 56, 99–118. doi:10.1016/j.ijplas.2014.01.006. Bouvard, J. L., Francis, D. K., Tschopp, M. A., Marin, E. B., Bammann, D. J., & Horstemeyer, M. F. (2013). An internal state variable material model for predicting the time, thermomechanical, and stress state dependence of amorphous glassy polymers under large deformation. International Journal of Plasticity, 42, 168–193. doi:10.1016/j.ijplas.2012.10.005. Chatziathanasiou, D., Chemisky, Y., Chatzigeorgiou, G., & Meraghni, F. (2016). Modeling of coupled phase transformation and reorientation in shape memory alloys under non-proportional thermomechanical loading. International Journal of Plasticity, 82, 192–224. doi:10.1016/j.ijplas.2016.03.005. Christ, D., & Reese, S. (2008). Finite-element modelling of shape memory alloys-a comparison between small-strain and large-strain formulations. Materials Science and Engineering A, 481–482, 343–346. doi:10.1016/j.msea.2006.11.174. Christ, D., & Reese, S. (2009). A finite element model for shape memory alloys considering thermomechanical couplings at large strains. International Journal of Solids and Structures, 46, 3694–3709. doi:10.1016/j.ijsolstr.2009.06.017. Cisse, C., Zaki, W., & Zineb, T. B. (2016a). A review of constitutive models and modeling techniques for shape memory alloys. International Journal of Plasticity, 76, 244–284. doi:10.1016/j.ijplas.2015.08.006. Cisse, C., Zaki, W., & Zineb, T. B. (2016b). A review of modeling techniques for advanced effects in shape memory alloy behavior. Smart Materials and Structures, 25, 103001. doi:10.1088/0964-1726/25/10/103001.

76

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

Dafalias, Y. F. (1984). The plastic spin concept and a simple illustration of its role in finite plastic transformations. Mechanics of Materials, 3, 361. doi:10. 1016/0167- 6636(84)90036- X. Doraiswamy, S., Rao, A., & Srinivasa, A. R. (2011). Combining thermodynamic principles with Preisach models for superelastic shape memory alloy wires. Smart Materials and Structures, 20, 85032. doi:10.1088/0964-1726/20/8/085032. Grabe, C., & Bruhns, O. T. (2008). On the viscous and strain rate dependent behavior of polycrystalline NiTi. International Journal of Solids and Structures, 45, 1876–1895. doi:10.1016/j.ijsolstr.2007.10.029. Grabe, C., & Bruhns, O. T. (2009). Path dependence and multiaxial behavior of a polycrystalline NiTi alloy within the pseudoelastic and pseudoplastic temperature regimes. International Journal of Plasticity, 25, 513–545. doi:10.1016/j.ijplas.20 08.03.0 02. Gu, X., Zaki, W., Morin, C., Moumni, Z., & Zhang, W. (2015). Time integration and assessment of a model for shape memory alloys considering multiaxial nonproportional loading cases. International Journal of Solids and Structures, 54, 82–99. doi:10.1016/j.ijsolstr.2014.11.005. Gurtin, M. E., & Anand, L. (2005). The decomposition f= f e f p, material symmetry, and plastic irrotationality for solids that are isotropic-viscoplastic or amorphous. International Journal of Plasticity, 21, 1686–1719. doi:10.1016/j.ijplas.20 04.11.0 07. Guthikonda, V. S., & Elliott, R. S. (2013). Modeling martensitic phase transformations in shape memory alloys with the self-consistent lattice dynamics approach. Journal of the Mechanics and Physics of Solids, 61, 1010–1026. doi:10.1016/j.jmps.2012.12.003. Hartl, D. J., Chatzigeorgiou, G., & Lagoudas, D. C. (2010). Three-dimensional modeling and numerical analysis of rate-dependent irrecoverable deformation in shape memory alloys. International Journal of Plasticity, 26, 1485–1507. doi:10.1016/j.ijplas.2010.01.002. Hazar, S., Zaki, W., Moumni, Z., & Anlas, G. (2015). Modeling of steady-state crack growth in shape memory alloys using a stationary method. International Journal of Plasticity, 67, 26–38. doi:10.1016/j.ijplas.2014.08.018. Idesman, A. V., Levitas, V. I., Preston, D. L., & Cho, J. Y. (2005). Finite element simulations of martensitic phase transitions and microstructures based on a strain softening model. Journal of the Mechanics and Physics of Solids, 53, 495–523. doi:10.1016/j.jmps.2004.10.001. Kan, Q., & Kang, G. (2010). Constitutive model for uniaxial transformation ratchetting of super-elastic NiTi shape memory alloy at room temperature. International Journal of Plasticity, 26, 441–465. doi:10.1016/j.ijplas.20 09.08.0 05. Kan, Q., Yu, C., Kang, G., Li, J., & Yan, W. (2016). Experimental observations on rate-dependent cyclic deformation of super-elastic NiTi shape memory alloy. Mechanics of Materials, 97, 48–58. doi:10.1016/j.mechmat.2016.02.011. Kelly, A., Stebner, A. P., & Bhattacharya, K. (2015). A micromechanics-inspired constitutive model for shape-memory alloys that accounts for initiation and saturation of phase transformation. Journal of the Mechanics and Physics of Solids, 1–28. doi:10.1016/j.jmps.2016.02.007. Khan, A. S., & Huang, S. (1995). Continuum theory of plasticity. John Wiley & Sons. Khater, H. A., Monnet, G., Terentyev, D., & Serra, A. (2014). Dislocation glide in fe-carbon solid solution: From atomistic to continuum level description. International Journal of Plasticity, 62, 34–49. doi:10.1016/j.ijplas.2014.06.006. Kimiecik, M., Jones, J. W., & Daly, S. (2016). The effect of microstructure on stress-induced martensitic transformation under cyclic loading in the SMA nickel-titanium. Journal of the Mechanics and Physics of Solids, 89, 16–30. doi:10.1016/j.jmps.2016.01.007. Lagoudas, D., Hartl, D., Chemisky, Y., Machado, L., & Popov, P. (2012). Constitutive model for the numerical analysis of phase transformation in polycrystalline shape memory alloys. International Journal of Plasticity, 32–33, 155–183. doi:10.1016/j.ijplas.2011.10.009. Lagoudas, D. C., Entchev, P. B., Popov, P., Patoor, E., Brinson, L. C., & Gao, X. (2006). Shape memory alloys, part II: Modeling of polycrystals. Mechanics of Materials, 38, 430–462. doi:10.1016/j.mechmat.20 05.08.0 03. León Baldelli, A., Maurini, C., & Pham, K. (2015). A gradient approach for the macroscopic modeling of superelasticity in softening shape memory alloys. International Journal of Solids and Structures, 52, 45–55. doi:10.1016/j.ijsolstr.2014.09.009. Levitas, V. I., & Ozsoy, I. B. (2009). Micromechanical modeling of stress-induced phase transformations. part 1. thermodynamics and kinetics of coupled interface propagation and reorientation. International Journal of Plasticity, 25, 239–280. doi:10.1016/j.ijplas.20 08.02.0 04. Levitas, V. I., & Preston, D. L. (2002). Three-dimensional Landau theory for multivariant stress-induced martensitic phase transformations. i. austenitemartensite. Physical review B, 66, 134206. doi:10.1103/PhysRevB.66.134206. Lexcellent, C., Boubakar, M., Bouvet, C., & Calloch, S. (2006). About modelling the shape memory alloy behaviour based on the phase transformation surface identification under proportional loading and anisothermal conditions. International Journal of Solids and Structures, 43, 613–626. doi:10.1016/j.ijsolstr. 20 05.07.0 04. Liang, C., & Rogers, C. A. (1990). One-dimensional thermomechanical constitutive relations for shape memory materials. Journal of Intelligent Material Systems and Structures, 1, 207–234. doi:10.1177/1045389X970 080 0402. Mehrabi, R., Andani, M. T., Elahinia, M., & Kadkhodaei, M. (2014). Anisotropic behavior of superelastic NiTi shape memory alloys; an experimental investigation and constitutive modeling. Mechanics of Materials, 77, 110–124. doi:10.1016/j.mechmat.2014.07.006. Mohd Jani, J., Leary, M., Subic, A., & Gibson, M. A. (2014). A review of shape memory alloy research, applications and opportunities. Materials and Design, 56, 1078–1113. doi:10.1016/j.matdes.2013.11.084. Morin, C., Moumni, Z., & Zaki, W. (2011). A constitutive model for shape memory alloys accounting for thermomechanical coupling. International Journal of Plasticity, 27, 748–767. doi:10.1016/j.ijplas.2010.09.005. Moumni, Z., Zaki, W., & Nguyen, Q. S. (2008). Theoretical and numerical modeling of solidâ;;solid phase change: Application to the description of the thermomechanical behavior of shape memory alloys. International Journal of Plasticity, 24, 614–645. doi:10.1016/j.ijplas.2007.07.007. Müller, C., & Bruhns, O. T. (2006). A thermodynamic finite-strain model for pseudoelastic shape memory alloys. International Journal of Plasticity, 22, 1658– 1682. doi:10.1016/j.ijplas.2006.02.010. Otsuka, K., & Wayman, C. M. (1999). Shape memory materials. Cambridge university press. Paranjape, H. M., Manchiraju, S., & Anderson, P. M. (2016). A phase field - finite element approach to model the interaction between phase transformations and plasticity in shape memory alloys. International Journal of Plasticity, 80, 1–18. doi:10.1016/j.ijplas.2015.12.007. Patoor, E., Lagoudas, D. C., Entchev, P. B., Brinson, L. C., & Gao, X. (2006). Shape memory alloys, part i: General properties and modeling of single crystals. Mechanics of Materials, 38, 391–429. doi:10.1016/j.mechmat.2005.05.027. Peng, X., Pi, W., & Fan, J. (2008). A microstructure-based constitutive model for the pseudoelastic behavior of NiTi SMAs. International Journal of Plasticity, 24, 966–990. doi:10.1016/j.ijplas.20 07.08.0 03. Peultier, B., Ben Zineb, T., & Patoor, E. (2006). Macroscopic constitutive law of shape memory alloy thermomechanical behaviour. application to structure computation by FEM. Mechanics of Materials, 38, 510–524. doi:10.1016/j.mechmat.2005.05.026. Peyroux, R., Chrysochoos, A., Licht, C., & Löbel, M. (1998). Thermomechanical couplings and pseudoelasticity of shape memory alloys. International Journal of Engineering Science, 36, 489–509. doi:10.1016/S0 020-7225(97)0 0 052-9. Popov, P., & Lagoudas, D. C. (2007). A 3-d constitutive model for shape memory alloys incorporating pseudoelasticity and detwinning of self-accommodated martensite. International Journal of Plasticity, 23, 1679–1720. doi:10.1016/j.ijplas.2007.03.011. Qidwai, M. A., & Lagoudas, D. C. (20 0 0). On thermomechanics and transformation surfaces of polycrystalline NiTi shape memory alloy material. International journal of plasticity, 16, 1309–1343. doi:10.1016/S0749-6419(0 0)0 0 012-7. Rao, A., & Srinivasa, A. R. (2013). Experiments on functional fatigue of thermally activated shape memory alloy springs and correlations with driving force intensity. SPIE smart structures and materials + nondestructive evaluation and health monitoring. International Society for Optics and Photonics. 86890T–86890T Rao, A., & Srinivasa, A. R. (2015). A three-species model for simulating torsional response of shape memory alloy components using thermodynamic principles and discrete preisach models. Mathematics and Mechanics of Solids, 20, 345–372. doi:10.1177/1081286514545917. Rao, A., Srinivasa, A. R., & Reddy, J. N. (2015). Design of shape memory alloy (SMA) actuators. Springer. Reedlunn, B., Churchill, C. B., Nelson, E. E., Shaw, J. A., & Daly, S. H. (2014). Tension, compression, and bending of superelastic shape memory alloy tubes. Journal of the Mechanics and Physics of Solids, 63, 506–537. doi:10.1016/j.jmps.2012.12.012.

J. Wang et al. / International Journal of Engineering Science 117 (2017) 51–77

77

Reese, S., & Christ, D. (2008). Finite deformation pseudo-elasticity of shape memory alloys constitutive modelling and finite element implementation. International Journal of Plasticity, 24, 455–482. doi:10.1016/j.ijplas.20 07.05.0 05. Reinhardt, W. D., & Dubey, R. N. (1996). Application of objective rates in mechanical modeling of solids. Journal of applied mechanics, 63, 692–698. doi:10. 1115/1.2823351. Reuss, A. (1929). Berechnung der fließgrenze von mischkristallen auf grund der plastizitätsbedingung für einkristalle.. ZAMM Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 9, 49–58. doi:10.10 02/zamm.19290 090104. Rizzoni, R., & Marfia, S. (2015). A thermodynamical formulation for the constitutive modeling of a shape memory alloy with two martensite phases. Meccanica, 50, 1121–1145. doi:10.1007/s11012-014-0078-8. Saadat, S., Salichs, J., Noori, M., Hou, Z., Davoodi, H., Bar-on, I., Suzuki, Y., & Masuda, a. (2002). An overview of vibration and seismic applications of NiTi shape memory alloy. Smart Materials and Structures, 11, 218–229. doi:10.1088/0964-1726/11/2/305. Saleeb, A. F., Dhakal, B., Hosseini, M. S., & Padula II, S. A. (2013). Large scale simulation of NiTi helical spring actuators under repeated thermomechanical cycles. Smart Materials and Structures, 22, 94006. doi:10.1088/0964-1726/22/9/094006. Saleeb, A. F., Padula, S. A., & Kumar, A. (2011). A multi-axial, multimechanism based constitutive model for the comprehensive representation of the evolutionary response of SMAs under general thermomechanical loading conditions. International Journal of Plasticity, 27, 655–687. doi:10.1016/j.ijplas. 2010.08.012. Sedlák, P., Frost, M., Benešová, B., Ben Zineb, T., & Šittner, P. (2012). Thermomechanical model for NiTi-based shape memory alloys including r-phase and material anisotropy under multi-axial loadings. International Journal of Plasticity, 39, 132–151. doi:10.1016/j.ijplas.2012.06.008. Shaw, J. A., & Kyriakides, S. (1995). Thermomechanical aspects of NiTi. Journal of the Mechanics and Physics of Solids, 43, 1243–1281. doi:10.1016/ 0 022-5096(95)0 0 024-D. Sittner, P., Hara, Y., & Tokuda, M. (1995). Experimental study on the thermoelastic martensitic transformation in shape memory alloy polycrystal induced by combined external forces. Metallurgical and Materials Transactions A, 26, 2923–2935. doi:10.1007/BF02669649. Song, D., Kang, G., Kan, Q., Yu, C., & Zhang, C. (2014). Non-proportional multiaxial transformation ratchetting of super-elastic NiTi shape memory alloy: Experimental observations. Mechanics of Materials, 70, 94–105. doi:10.1016/j.mechmat.2013.12.003. Souza, A. C., Mamiya, E. N., & Zouain, N. (1998). Three-dimensional model for solids undergoing stress-induced phase transformations. European Journal of Mechanics - A/Solids, 17, 789–806. doi:10.1016/S0997-7538(98)80 0 05-3. Tanaka, K., Kobayashi, S., & Sato, Y. (1986). Thermomechanics of transformation pseudoelasticity and shape memory effect in alloys. International Journal of Plasticity, 2, 59–72. doi:10.1016/0749-6419(86)90016-1. Teeriaho, J. P. (2013). An extension of a shape memory alloy model for large deformations based on an exactly integrable eulerian rate formulation with changing elastic properties. International Journal of Plasticity, 43, 153–176. doi:10.1016/j.ijplas.2012.11.009. Thamburaja, P. (2005). Constitutive equations for martensitic reorientation and detwinning in shape-memory alloys. Journal of the Mechanics and Physics of Solids, 53, 825–856. doi:10.1016/j.jmps.20 04.11.0 04. Thamburaja, P., & Anand, L. (2001). Polycrystalline shape-memory materials: Effect of crystallographic texture. Journal of the Mechanics and Physics of Solids, 49, 709–737. doi:10.1016/S0 022-5096(0 0)0 0 061-2. Thamburaja, P., & Anand, L. (2002). Superelastic behavior in tension-torsion of an initially-textured ti-ni shape-memory alloy. International Journal of Plasticity, 18, 1607–1617. doi:10.1016/S0749-6419(02)0 0 031-1. Toi, Y., Lee, J.-B., & Taya, M. (2004). Finite element analysis of superelastic, large deformation behavior of shape memory alloy helical springs. Computers & Structures, 82, 1685–1693. doi:10.1016/j.compstruc.2004.03.025. Wagner, M.-X., & Windl, W. (2008). Lattice stability, elastic constants and macroscopic moduli of NiTi martensites from first principles. Acta Materialia, 56, 6232–6245. doi:10.1016/j.actamat.2008.08.043. Williams, E., & Elahinia, M. H. (2008). An automotive SMA mirror actuator: Modeling, design, and experimental evaluation. Journal of Intelligent Material Systems and Structures, 19, 1425–1434. doi:10.1177/1045389X07087328. Wu, W., Qi, M., Liu, X.-P., Yang, D.-Z., & Wang, W.-Q. (2007). Delivery and release of nitinol stent in carotid artery and their interactions: A finite element analysis. Journal of Biomechanics, 40, 3034–3040. doi:10.1016/j.jbiomech.2007.02.024. Wu, X. D., Sun, G. J., & Wu, J. S. (2003). The nonlinear relationship between transformation strain and applied stress for nitinol. Materials Letters, 57, 1334–1338. doi:10.1016/S0167- 577X(02)00983- 7. Xiao, H. (2014). An explicit, straightforward approach to modeling SMA pseudoelastic hysteresis. International Journal of Plasticity, 53, 228–240. doi:10.1016/ j.ijplas.2013.08.010. Xiao, H., Bruhns, O., & Meyers, a. (2004). Explicit dual stress-strain and strain-stress relations of incompressible isotropic hyperelastic solids via deviatoric hencky strain and cauchy stress. Acta Mechanica, 168, 21–33. doi:10.10 07/s0 0707-0 04-0 074-5. Xiao, H., Bruhns, O. T., & Meyers, A. (1998). On objective corotational rates and their defining spin tensors. International journal of solids and structures, 35, 4001–4014. doi:10.1016/S0020- 7683(97)00267- 9. Yin, H., He, Y., & Sun, Q. (2014). Effect of deformation frequency on temperature and stress oscillations in cyclic phase transition of NiTi shape memory alloy. Journal of the Mechanics and Physics of Solids, 67, 100–128. doi:10.1016/j.jmps.2014.01.013. Yu, C., Kang, G., & Kan, Q. (2014). Study on the rate-dependent cyclic deformation of super-elastic NiTi shape memory alloy based on a new crystal plasticity constitutive model. International Journal of Solids and Structures, 51, 4386–4405. doi:10.1016/j.ijsolstr.2014.09.006. Yu, C., Kang, G., & Kan, Q. (2015a). A micromechanical constitutive model for anisotropic cyclic deformation of super-elastic NiTi shape memory alloy single crystals. Journal of the Mechanics and Physics of Solids, 82, 97–136. doi:10.1016/j.jmps.2015.05.012. Yu, C., Kang, G., Kan, Q., & Zhu, Y. (2015b). Rate-dependent cyclic deformation of super-elastic NiTi shape memory alloy: Thermo-mechanical coupled and physical mechanism-based constitutive model. International Journal of Plasticity, 72, 60–90. doi:10.1016/j.ijplas.2015.05.011. Zaki, W. (2010). An approach to modeling tensile-compressive asymmetry for martensitic shape memory alloys. Smart Materials and Structures, 19, 25009. doi:10.1088/0964-1726/19/2/025009. Zaki, W., & Moumni, Z. (2007a). A 3d model of the cyclic thermomechanical behavior of shape memory alloys. Journal of the Mechanics and Physics of Solids, 55, 2427–2454. doi:10.1016/j.jmps.2007.03.011. Zaki, W., & Moumni, Z. (2007b). A three-dimensional model of the thermomechanical behavior of shape memory alloys. Journal of the Mechanics and Physics of Solids, 55, 2455–2490. doi:10.1016/j.jmps.2007.03.012.