Microstructure-based finite deformation modelling of super-elastic NiTi shape memory alloy considering various inelastic mechanisms

In this study the mechanisms of transformation-induced plasticity (TRIP), detwinning-induced plasticity (DIP), and accumulation of residual martensite, are incorporated into a finite-deformation crystal plasticity model of NiTi SMA for the first time. The constitutive model is constructed at the single-crystal scale and also includes phase transformation and detwinning mechanisms. Using a proposed Helmholtz free energy, the driving forces for inelastic mechanisms are derived within the framework of thermodynamics. The constitutive model has been implemented in the Abaqus/Explicit finite-element program, using VUMAT subroutine to simulate a polycrystalline material. Considering various orientations for crystals, the effect of texture on tension-compression asymmetry is investigated. It is shown that different textures may cause stiffer, softer, or similar response in compression compared to tension. Due to the incorporation of the effect of residual martensite, the model provides accurate predictions of experimentally measured residual strain. The incorporation of the aforementioned inelastic deformation mechanisms is shown to accurately capture the key features related to cyclic loading. Finally, the effect of DIP in compressive cyclic loading of NiTi SMA is investigated. Considering the effect of DIP causes a less negative peak stress and a more negative residual strain to be observed in the strain-controlled cyclic compression-unloading tests.

The behavior of NiTi SMA is based on the temperature-dependent reversible phase transformation from the high symmetry austenite phase to the low symmetry martensite phase and vice versa (Lagoudas, 2008). At temperatures higher than the austenite finish-temperature, application of stress to austenite will cause strain generation, which is recoverable upon unloading (super-elasticity). Upon cooling, in the absence of stress, the high temperature austenite phase will transform to the low temperature martensite phase (austenite to twinnedmartensite transformation). Application of enough stress to twinned-martensite will cause strain generation during loading, which is not recovered upon unloading, and requires subsequent heating to be recovered (shape-memory effect).
Many constitutive models have been developed to predict the behavior of SMAs based on experimental observations. The models can be categorized into two main groups: macroscopic phenomenological models and microstructure-based models. Typical macroscopic constitutive models include Auricchio et al. (Auricchio et al., 2007), Lagoudas et al. (Lagoudas et al., 2012), Panico and Brinson (Panico and Brinson, 2007), Zaki and Moumni (Zaki and Moumni, 2007), Arghavani et al. (Arghavani et al., 2010), and Xu et al. (Xu et al., 2019). While macroscopic models do not consider the microscale mechanisms underlying SMA behavior, they are readily implementable in finite element codes for structural analysis and design.
Microstructure-based constitutive models have been proposed to incorporate the underlying microscale physical mechanisms underlying SMA behavior. The popularity of crystal plasticity-based models is due to considering various variants of martensite with different morphological features and their evolution during deformation. These models are developed at the single-crystal scale and then, adopting finite element method or scale-transition rules, transformed to the poly-crystal scale. Representative constitutive models include Thamburaja and Anand (Thamburaja and Anand, 2003, 2002, Anand and Gurtin (Anand and Gurtin, 2003), Thamburaja (Thamburaja, 2005), Thamburaja et al. , Wang et al. (Wang et al., 2008), Manchiraju and Anderson (Manchiraju and Anderson, 2010), Yu et al. (Yu et al., 2015(Yu et al., , 2013(Yu et al., , 2012, Xiao et al. (Xiao et al., 2018), and Dhala et al. (Dhala et al., 2019). It should be noted that some of the above-mentioned constitutive models do not consider the effect of large deformations in the modeling framework. In practical applications, the strains and rotations applied to SMA can be quite large such that the theories of small deformation cannot be used to reasonably predict the behavior of these materials. In recent years, the behavior of NiTi SMAs under finite deformations has been studied by some researchers. The effect of phase transformation on the mechanical and thermomechanical response of NiTi SMA has been investigated by Thamburaja and Anand (Thamburaja and Anand, 2003, 2002, and Anand and Gurtin (Anand and Gurtin, 2003). The effects of detwinning and reorientation have been incorporated into constitutive modeling of NiTi SMA by Thamburaja (Thamburaja, 2005). This model has been modified to consider austenite to martensite phase transformation, as well . Manchiraju and Anderson (Manchiraju and Anderson, 2010) considered the effects of phase transformation and austenite plasticity, which results from phase transformation at high temperatures. It should be noted that the residual strain observed in this modeling is due to austenite plasticity, occurs only at sufficiently high temperatures, and differs from the mechanism associated with the accumulation of residual martensite, which can be activated at any temperature. The effects of austenite-plasticity and martensite-plasticity on the mechanical behavior of super-elastic NiTi have been investigated by Dhala et al. (Dhala et al., 2019). During phase transformation, some martensite variants are locked, due to increase of defects, and cannot return to austenite phase. This process, called accumulation of residual martensite, is known to be one of the main mechanisms in the deformation of super-elastic NiTi SMA (Kang et al., 2009;Yu et al., 2013).
It should be noted that the effect of residual martensite has not been considered in the abovementioned finite deformation constitutive models, while it is evident in experimental observations. On the other hand, some kind of local stress field and plasticity can be activated in SMAs which is different from conventional plasticity (Lagoudas, 2008). The plastic deformation which is created to overcome the local stress field at the interface of asutenite and martensie phases is called transformation-induced plasticity (Yu et al., 2015), which has significant effect on predicting the behavior of NiTi SMA under cyclic loading (Lagoudas, 2008). In addition, the experimental results of Liu et al. (Liu et al., 1998) show a high density of dislocations at the interfaces of lattice correspondent variants during compressive loading, but no significant plastic deformation is detected in such areas under tension. Detwinning-induced plasticity is referred to that part of plastic deformation which is induced to relax the local stress field at the interfaces of lattice correspondent variants (Ebrahimi et al., 2020). As demonstrated by Ebrahimi et al. (Ebrahimi et al., 2020), DIP has considerable effect on compressive cyclic loading of NiTi SMAs. It should be noted that the cyclic behavior of NiTi SMA under large deformations, considering the crucial mechanisms of TRIP and DIP, has not been investigated in previous studies.
The goal of the current study is to predict the behavior of NiTi under monotonic and cyclic loadings considering large deformations. Using crystal plasticity formulation, diverse mechanisms of phase transformation, transformation-induced plasticity, detwinning, detwinning-induced plasticity, and accumulation of residual martensite have been incorporated into finite-deformation modeling of NiTi SMA. The single-crystal constitutive model has been implemented in the Abaqus/Explicit finite-element program, using VUMAT subroutine to obtain the polycrystalline response. The model is shown to capture various features observed during monotonic and cyclic loading of super-elastic NiTi.

Microstructure-based constitutive model
2.1. Definition of various deformation mechanisms in the framework of thermodynamics Using multiplicative decomposition, the deformation gradient tensor F can be decomposed into elastic, e F , and inelastic, p F , parts (Lubarda, 2004;Thamburaja, 2005): where p F is composed of several parts related to different mechanisms of phase transformation (tr), detwinning (det), transformation-induced plasticity (TRIP), and detwinning-induced plasticity (DIP).
The velocity gradient tensor, L , can also be divided into elastic, e L , and inelastic, p L , parts (Dhala et al., 2019;Erdle and Böhlke, 2017;Fohrmeister et al., 2019): The plastic velocity gradient tensor, p L , comprises various velocity gradient tensors related to phase transformation tr L , transformation-induced plasticity TRIP L , detwinning det L , and detwinning-induced plasticity DIP L .
Due to experimentally observed accumulation of residual martensite (Kang et al., 2009) where  , rev  , and ire  are the total, reversible, and irreversible martensite volume fractions, respectively.   is the martensite volume fraction for the -th  martensite variant.
The orientation tensor  P of the -th  martensite variant can be expressed by (Anand and Gurtin, 2003;Thamburaja and Anand, 2001): where b and m denote transformation direction and habit plane normal of the -th  martensite variant, respectively. The magnitudes of theses parameters, related to different martensite variants, are presented by (Wang et al., 2008).
In the absence of stress, by reducing the temperature to lower than the martensite finishtemperature, austenite will transform to twinned martensite. Consequently, some habit plane variants (HPVs) will be formed. Each HPV is composed of two lattice correspondent variants (LCVs) or sub-variants (Ebrahimi et al., 2020;Yu et al., 2014). The stress-induced movement of interfaces between sub-variants is called detwinning (Thamburaja, 2005 where a and w denote shear direction and unit vector normal to detwinning plane for the -th  martensite variant, respectively. The magnitudes of theses parameters, related to different martensite variants are presented by .
According to crystal plasticity (Dhala et al., 2019;Rengarajan and Rajagopal, 2001;Thamburaja et al., 2005), the inelastic velocity gradient tensor is proposed here such that: , and represents the rate change of sub-variants' volume fractions. Since TRIP is the plastic deformation occurred to relax the local stress field at the austenite-martensite interfaces, its orientation tensor is assumed to be the same as phase transformation mechanism (Yu et al., 2015(Yu et al., , 2013. The Helmholtz free energy of a representative volume element (RVE) is considered to be composed of different parts related to elastic energy, e , chemical energy, ch , the part related to the influence of internal stress, int , and the part reflecting the effects of hardening or softening during phase transformation, tr rev : The elastic part of free energy can be expressed by (Yu et al., 2015): where () is the fourth-order stiffness tensor, considered as a function of the volume fraction of martensite according to the rule of mixtures, AM . e E is the Green-Lagrange elastic strain tensor.
The chemical energy is used to account for the effect of temperature on the phase transformation (Thamburaja, 2005): T is the phase equilibrium temperature, at which the free energies of the austenite and martensite phases are equal. T h is the absorbed or released latent heat per unit reference volume at the temperature of T during phase transformation.
During cyclic loading, some reduction in the start stress of phase transformation will be observed by increasing number of cycles (Kang et al., 2009). It is concluded that local stress fields due to different inelastic mechanisms will assist the occurrence of phase transformation, and consequently lower the start stress of forward transformation in cyclic loading. Using this fact and considering the work of Yu et al. (Yu et al., 2015), the form of int is proposed here such that: B is an internal variable, known as internal stress.
tr rev is used to account for the hardening or softening effects during phase transformation. Hardening or softening will be observed in the response of shape memory alloy by any increase or decrease in the amount of tr rev , respectively (Yu et al., 2015). Since internal stress is related to local stress field caused by different variants of martensite, it can be divided into 24 components (Yu et al., 2015): where  B is the norm of  B (the internal stress caused by the -th  martensite variant). . rev is the accumulated reversible martensite volume fraction, which is shown in the rateform by (Yu et al., 2013): The experimentally observed transformation hardening in cyclic loading of SMA (Kang et al., 2009), is considered in the model by adopting X  as a function of the reversible martensite volume fraction (Yu et al., 2015). The Clausius inequality, in finite-deformation regimes, can be expressed by (Anand and Gurtin, 2003;Thamburaja, 2005): where S is the first Piola-Kirchhoff stress tensor. Considering quasi-static loading and neglecting the transformation latent heat and the heat from internal dissipation, the third term in Eq. (17) can be neglected. Therefore, using the equality ** : : (17) can be re-written as: (Thamburaja, 2005) is the second Piola-Kirchhoff stress tensor associated with the elastic part, and T is the Cauchy stress tensor.
Using Eqs. (10) and (11) (21) Since either positive or negative values can be assigned to the parameters e E and , according to the method proposed by Coleman and Noll (Noll et al., 1974), the necessary condition for Eq. (21) to be satisfied is that:  (Lagoudas and Entchev, 2004;Patoor et al., 2006), the evolution equations for the reversible martensite volume fraction are as follows: where Y is the resistance to phase transformation.
The magnitude of Y at the start of forward transformation is denoted by  Experimental observations reveal an asymmetry in the stress-strain response of NiTi SMA under compressive and tensile loadings. While a flat stress-plateau is observed during tension, a rapid strain hardening is reported during compression. This difference in material behaviour is addressed to the generation and migration of dislocations during compressive loading (Liu et al., 1998). On the other hand, it is known that residual martensite is referred to the part of martensite which is locked by defects and can not return to austenite phase (Yu et al., 2015).
Considering these points, it is concluded that the amount of residual martensite will be higher in compression than in tension. This conclusion is supported by the experimental results of Thamburaja and Anand (Thamburaja and Anand, 2001). To reflect the difference of residual where sat irev is the saturated value of irev , and 1 b controls the saturation value of irev in compressive loading.  Table 1 represents a summary of the formulation used in the proposed constitutive model.

Detwinning
The elastic Green-Lagrange strain tensor, in terms of the elastic deformation gradient tensor, is calculated such that (Anand and Gurtin, 2003): Using e E and * T (Eq. (22)), the Cauchy stress tensor can be calculated.

Material parameters identification
For austenite phase, with cubic crystal structure, three independent elastic coefficients are used to define the fourth-order stiffness tensor (Jamal et al., 2014): The corresponding independent elastic constants for austenite phase are denoted by 11 A , 12 A , and 44 A , and considered to be 130, 98, and 34 GPa, respectively (Brill et al., 1991;Thamburaja and Anand, 2001). The modulus of elasticity of austenite is usually two to three times greater than the modulus of elasticity of martensite, while the Poisson's ratios are almost the same. Therefore, the elastic constants of martensite ( 11 12 44 ,, MMM ) are considered equal to half of the corresponding constants of austenite (Thamburaja and Anand, 2001).
Sufficiently large values are selected for the parameters tr n and det n , to reproduce the rateindependent response of NiTi SMAs.
The evolution law of reversible martensite volume fraction during forward and reverse transformations is described by the parameter Y , while 5 b is responsible for difference in evolution of reversible martensite during loading and unloading. The width of hysteresis loop for the first and steady-state cycles is controlled by the parameters 0 Y and d Y , respectively.
Using the curves related to start-stress of detwinning, det K can be obtained.
As experimentally observed (Liu et al., 1998) and computationally predicted (Ebrahimi et al., 2020), DIP has no considerable effect in tension. Therefore, the parameters related to DIP cannot be obtained by using the data related to tensile loading of SMA. On the other hand, as commented by Kang et al. (Kang et al., 2009), TRIP affects both the accumulated peak strain and accumulated residual strain during uniaxial cyclic deformation, while the influence of residual martensite is mainly focused on the accumulated residual strain. Thus, the data of accumulated peak strain in tensile cyclic loading can be used to obtain the parameters   (Thamburaja and Anand, 2001).
It should be noted that the above-mentioned parameters are related to a single crystal and should be calibrated using the corresponding single crystal experimental data. In the lack of such experimental results for single crystal, a trial-and-error method is used to calibrate the parameters using the polycrystalline experimental data.
The adopted parameters for predicting the behavior of NiTi SMA by considering the effect of residual martensite, are listed in Table 2. , , The parameters used to predict the cyclic behavior of NiTi, considering the effects of DIP and TRIP are presented in Tables 3 and 4, respectively. , ,  , ,

Results and discussion
In polycrystalline shape memory alloys, manufacturing process can affect the crystallographic texture which is the main reason for tension-compression asymmetry observed in such materials (Thamburaja and Anand, 2001). As commented by Gall et al. (Gall et al., 2001), based on the method used to form an alloy, three different textures with different effects on the overall response can be formed: a   100 texture making the mechanical response stiffer in tension than in compression, a   111 texture causing softer response in tension than in compression, and a   110 texture which almost has no effect on the mechanical responses in tension or compression.
It should be noted that the parameters related to phase transformation direction, normal to the phase transformation plane, and detwinning are expressed in a crystal local coordinate system.
Therefore, a transformation of the deformation gradient tensor from poly-crystal to singlecrystal coordinate system is performed at the beginning of each increment. At the end of calculations, the calculated stress tensor should be transformed back to poly-crystal coordinate system, which is done by using the rotation matrix Q .
Appropriate ranges of Euler angles have been calculated in order to model different textures.
To this end, the single-crystal and poly-crystal coordinates are considered as local and global coordinate systems, respectively. According to the Bunge system, the orientation of each crystal can be identified using 3 Euler angles,  ,  , and  . In this system, the rotation matrix in 3D is expressed as (Depriester, 2018): cos cos sin sin cos sin cos cos sin cos sin sin ( , , ) cos sin sin cos cos sin sin cos cos cos cos sin sin sin cos sin In order to model crystals with [111] texture, the [111] direction of each crystal should be nearly parallel to the loading direction. A set of randomly generated Euler angles without any preferred orientation, is used to simulate a random texture/untextured polycrystalline material.
Using the proposed model, the behavior of a single-crystal for different textures under tensile and compressive strain-controlled loading-unloading is presented in Fig. 1. As shown in Fig. 1, for [111] texture the material behavior is harder in compression than in tension. In contrast, for the two other textures, and also for the randomly oriented crystals, the behavior in compression is softer in compression than in tension. Tension-compression differences are less pronounced for [110] and random textures. Based on the asymmetrical tension/compression response observed in the experimental test of Thamburaja and Anand (Thamburaja and Anand, 2001), the [111] texture is adopted to model the behavior of NiTi SMA in this work.
Since the number of crystals used to model the polycrystalline material can affect the predicted behavior, an appropriate number of crystals should be selected. The number of crystals should be as small as possible to reduce the computational cost and on the other hand should be large enough to prevent considerable change in material behavior by increasing number of crystals afterwards. Therefore, a sensitivity analysis is performed to determine the number of crystals in a representative volume element required for a converged stress-strain response. Results for 1, 8, 27, 64, and 125 crystals are shown in Fig. 2. Converged stress-strain behavior is obtained for 27, 64, and 125 crystals. Based on this sensitivity analysis 64 crystals will be used to model the polycrystalline behavior for the remainder of this study.  Fig. 4. It should be mentioned that the same sets of Euler angles have been used to capture [111] texture in both cases, but the Euler angles assigned to each crystal are randomly selected from that set in each case. In order to evaluate the proposed model, in Fig. 5 computed results for tensile and compressive loadings are directly compared with corresponding NiTi experimental data reported by Thamburaja and Anand (Thamburaja and Anand, 2001). The parameters used to model Nitinol's behavior are presented in Table 2. Several key features should be noted: (i) the experimental results show residual strain after unloading, in both compression and tension. This is accurately predicted by the proposed model due to the incorporation of the effect of residual martensite. Experimental results of Liu et al. (Liu et al., 1998) show a stressplateau in the stress-strain response during tensile loading, while the response is quickly strain hardened and no plateau is observed during compressive loading, which is related to the generation of dislocations (Liu et al., 1998). On the other hand, the increase of defects causes the residual martensite to be accumulated (Yu et al., 2015). Using these 2 facts and the experimental results of Thamburaja and Anand (Thamburaja and Anand, 2001), it is concluded that the amount of residual martensite should be higher in compression than in tension. In order to reflect the difference of residual martensite volume fraction in tension and compression, two separate values are assigned to the saturation value of this parameter for 2 loading cases in Eq. (34).
(ii) The proposed model accurately predicts the significant experimentally measured asymmetry between tension and compression, both in terms of loading and unloading behavior.
As shown in Fig. 1     The proposed model captures the main features observed in the experimental tests (Kang et al., 2009;Song et al., 2014). In summary, with increasing number of cycles, the start stress of austenite to martensite transformation decreases, the width of the hysteresis loop decreases, and the residual strain approaches a steady state, indicating that the accumulation of residual martensite approaches a saturation level.
In order to investigate the effect of DIP in large deformation, the stress-strain diagrams of the material under compressive loading to a minimum strain of -5% and unloading to zero stress are shown for the first, fifth and tenth cycles in Fig. 9. As shown in Fig. 9, the effect of detwinning-induced plasticity in the first cycle is negligible but is more pronounced with increasing number of cycles. Similar to the results presented by Ebrahimi et al. (Ebrahimi et al., 2020), considering the effect of DIP reduces the absolute value of the maximum stress and increases the absolute value of the residual strain in the straincontrolled cyclic loading of NiTi SMA.
The contour plots of shearing deformation caused by detwinning-induced plasticity for two variants of martensite are depicted in Fig. 10 Fig. 10. contour plots of shearing deformation, due to detwinning-induced plasticity, in two variants of martensite.

Conclusions
In this study a microstructure-based model for NiTi SMAs is developed by incorporating transformation-induced plasticity (TRIP), detwinning-induced plasticity (DIP), and accumulation of residual martensite into a finite-deformation crystal plasticity framework. The constitutive model is constructed at the single-crystal scale and also includes phase transformation and detwinning mechanisms. Using a proposed Helmholtz free energy, the driving forces for inelastic mechanisms are derived within the framework of thermodynamics.
The constitutive model has been implemented in the Abaqus/Explicit finite-element program, using VUMAT subroutine to simulate a polycrystalline material. Considering various orientations for crystals, the effect of texture on tension-compression asymmetry is investigated. It is shown that different textures may cause harder, softer, or similar response in compression compared to tension, which is in good agreement with the comments of Gall and Sehitoglu (Gall and Sehitoglu, 1999), Gao et al. (Gao et al., 2000), and Thamburaja and Anand (Thamburaja and Anand, 2001).
Due to the incorporation of the effect of residual martensite, the model provides accurate predictions of experimentally measured residual strain (Thamburaja and Anand, 2001). The incorporation of the inelastic DIP and TRIP deformation mechanisms and using appropriate form of free energy is shown to accurately capture the key features related to cyclic loading, including decrease in start stress of austenite to martensite transformation, decrease in width of the hysteresis loop, and saturation of residual strain (Song et al., 2014). Finally, the effect of detwinning-induced plasticity in compressive cyclic loading of NiTi SMA is investigated. DIP affects the results of strain-controlled cyclic compression-unloading simulations by causing the peak stress to be less negative and residual strain to be more negative in different cycles (Ebrahimi et al., 2020).  (Scheinert et al., 2005) and the analysis presented in the current study and the role of inelastic microstructural mechanisms in cyclic response of SMAs can potentially be used to guide the material and structural design of next-generation stents.