A polynomial constitutive model of shape memory alloys based on kinematic hardening

This paper describes the derivation of a phenomenological model for shape memory alloys under the framework of classical plasticity theory. The proposed model combines the Souza constitutive approach with kinematic hardening; the model requires solving only one nonlinear equation rather than several nonlinear ones, therefore increasing the computational efficiency and convergence. Moreover, the original Souza model is improved by adding an odd polynomial function to describe the phase transformation of the shape memory alloys, making it possible to use a lower number of parameters for the inverse identification of the constitutive properties of SMAs from simple tensile tests. A tangent stiffness formulation is also derived to simulate the variation of the elastic modulus during the phase transformation. The tangent stiffness formulation proposed here extends the one used in classical plasticity and improves the convergence of the proposed model. The reliability and fidelity of the model described in this work are benchmarked against experimental data and other models. The numerical results show that the proposed phenomenological approach can describe well the pseudoelasticity and shape memory effect of shape memory alloys. The formulation described in this paper can be readily generalized to finite strains and other formulations based on existing formulations related to classical plasticity theory.

(Some figures may appear in colour only in the online journal) * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Shape memory alloys (SMAs) have been widely used in aerospace [1,2], instrumentation [3,4] and biomedical applications [5,6] because of their pseudoelasticity and shape memory effect. In recent decades, scholar and engineers have extensively developed constitutive models and prototypes of devices involving the mechanical behavior of SMAs [7,8].
Studies have also shown that the pseudoelasticity and shape memory effect of SMAs are caused by the phase changes of twinned and detwinned martensite and austenite under different temperatures and stresses [7,[9][10][11].
SMA constitutive models can be approximately divided into two categories: meso-crystal ones and macroscopic thermodynamic phenomenological approaches [12]. We introduce first the theoretical fundamentals of the meso-crystal constitutive model and its advantages and disadvantages. The meso-crystal constitutive model [13][14][15] is based on the lattice deformation of martensite variants. The single-crystal constitutive relationship is derived by applying the WLR crystallography theory, which is used to predict the crystallographic characteristic parameters, with the multi-crystal constitutive relationship being then derived by combining a selfconsistent method. Meso-crystal models aim to describe the microstructural behaviors of SMA grains (i.e. nucleation and interface movements), which can help to understand the mechanisms underpinning the mechanics of SMAs. Anand, Wang, and Xiao et al [16][17][18][19] considered the volume fraction of a group of martensite variants as the internal variables and introduced an interaction energy matrix to describe the interaction between martensite variants. The Gibbs free energy function can then be obtained by using the meso-mechanical theory, resulting in the derivation of an SMA constitutive model. This model can describe well the reorientation [20,21] of martensite, detwinning-induced plasticity, and cyclic degeneration [22][23][24]. However, the material parameters required for the meso-crystal approach need to be extracted from specialistic experiment tests. Moreover, the meso-crystal model involves the use of a significant number of internal variables, leading to very time-consuming calculations that are not suitable for implementations within finite elements frameworks. The phenomenological constitutive model is more suitable to be incorporated with finite element models because of its simplicity in predicting the thermomechanical response of SMA actuators [25]. The phenomenological model only uses macroscopic state variables to describe the overall mechanical behavior of the SMA materials, with no need of using microscale parameters. This approach is therefore well suited to practical engineering simulations. The material parameters associated to phenomenological models can are usually determined using classical experimental measurements, and the constitutive equations can be also embedded within finite elements codes.
Phenomenological models are broadly categorized in two classes, the Lagoudas and the Souza models. These models involve the derivation of the free energy function by selecting appropriate internal variables [26]. The fundamentals about the Lagoudas model are described here first. The model developed by Lagoudas and Boyd et al [26][27][28] is based on the classical Tanaka model [29]. In this model, the Gibbs free energy function of the SMA is derived with the martensite volume fraction and the phase transformation strain as the internal variables. The governing equations of the control variables and the thermodynamic driving force of martensite are derived using the Clausius-Duhem inequality. The phase transformation yield condition used in the Lagoudas model states that when the thermodynamic driving force reaches the critical value, the SMA will undergo phase transformation. The thermodynamic driving force is obtained by taking the partial derivative of the Gibbs free energy with respect to martensite volume fraction, which makes the physical meaning of thermodynamic driving force and critical driving force not very intuitive and easy to understand. Lagoudas and Brinson et al [30][31][32] also proposed several phase change hardening functions in exponential and trigonometric forms to better characterize the mechanical behaviors of the SMA phase transformation. Based on the temperature-stress phase diagram, Popov and Lagoudas [33] introduced the volume fraction of twinned martensite and detwinned and austenite to describe the transformation between different metal phases and complex mechanical behaviors in SMAs under various thermal and static loading paths. Zaki and Moumni [34,35] developed the Lagoudas model by introducing the martensite orientation strain, which enables to describe the reorientation and plastic deformation [36] of microstructures in SMAs. On the basis of Lagoudas model, Chemisky, Lagoudas, Hartl and Zhang et al have performed several analyses on the tensioncompression asymmetry [37,38], finite strain [39,40], large deformation [41], non-proportional loading [39,42], load rate influence [43][44][45], thermomechanical coupling [46][47][48][49], cyclic loading [32,50,51] and plastic fatigue [52][53][54][55][56][57] of SMAs, which have extended the assessment of the feasibility of this model.
The second class of SMA phenomenological models is the Souza theoretical approach. The basic theory underpinning this model and its developments are described as follows. Based on the use of temperature-stress phase diagram, Souza et al [58] constructed Helmholtz free energy using the transformation strain as the internal variable and derived a 3D macroscopic constitutive model, in which the phase transformation driving force is obtained by taking the partial derivative of the Helmholtz free energy with respect to the transformation strain. The phase transformation driving force has a similar meaning to stress, which makes its physical significance intuitive and easy to understand. On this basis, Souza introduced a transformation yield function, which was further improved in Souza's follow-up work [59]. This yield function describes the occurrence of the phase transformation within SMAs. This phase transformation yield function also indicates that SMAs will undergo phase transformation when the driving force of the phase transformation goes beyond the elastic region [60]. Souza et al [59] also derived an evolution criteria based on the Drucker-Prager inequality in terms of the deviatoric transformation stress and the mean transformation stress. The model can describe the shape memory effect and pseudoelasticity of SMAs. Arghavani et al [61] expressed the transformation strain as the product of a scalar and an orientation tensor, which also enabled the Souza model to describe the reorientation phenomenon. Auricchio and Gu et al [62,63] have also optimized the numerical integration algorithm adopted in the original Souza model to provide a more stable convergence. Altas et al [64] introduced exponential functions to describe the phase transformation hardening characteristics of SMAs during the forward and reverse phase transformations. Based on the Souza model, Chen et al [65] considered the plastic deformation after the completion of phase transition and the coupling between the plastic strain and the transformation strain. Kang and Kan et al [66][67][68] first observed ratchetting deformations of super-elastic NiTi alloy under uniaxial stress-controlled cyclic tests and then evaluated the dependence of the ratchetting versus the applied stress, stress rates, temperatures and types of loading. The Authors indicated that the phenomenon of the transformation ratchetting is caused by transformation-induced plasticity (TRIP) and plastic yield. They therefore introduced energy dissipation into the equations describing the damage evolution and established a predictive multi-axial fatigue life model for NiTi-SMAs [69]. Kang and Petrini et al [70,71] introduced the cumulated induced-martensite volume fraction variable into the phenomenological model to describe the cyclic ratchetting behavior. Their model can comprehensively consider the evolution of TRIP, tension-compression asymmetry, the dependence of transformation ratchetting on the stress and phase transformation hardening. Similarly, Auricchio, Zaki, Moumni, and Petrini et al have proposed several refinements to the Souza model in terms of martensite reorientation [72][73][74][75], tension-compression asymmetry [76], cyclic behavior [77], finite strain [72,[78][79][80], and plastic coupling [81].
The Souza model has a clear physical implication, and its constitutive equations are similar to those used in classical plastic mechanics theory. The numerical solution of the model is however difficult, because a set of nonlinear equations needs to be solved to obtain the transformation strain. To deal with this problem, this paper proposes a macroscopic phenomenological constitutive model for SMAs within the framework of classical plasticity theory. This approach is based on the use of the Souza model and on adopting kinematic hardening principles. The second section of this paper describes the internal relation between temperature-stress phase diagrams and stress-strain curves and derives the transformation relation related to the different internal variables. In the third section, a 1D SMA macroscopic constitutive mode like a kinematic hardening one is firstly derived by using Helmholtz free energy function. Then, an odd polynomial function is introduced to describe the mechanical response during the phase transformation and an analytical expression related to a consistent tangent stiffness is given. In the fourth section of this paper, the midpoint integration algorithm is used to establish the solution for the numerical integral of the polynomial macroscopic constitutive model. When phase transformation occurs, the proposed algorithm requires to solve one nonlinear equation only to complete the stress correction; this can help to overcome the intrinsic difficulty of the Souza model, which requires the solution of a set of nonlinear equations instead. In the fifth section, numerical solutions are provided to verify the reliability of the proposed model to accurately describe the complex mechanical behaviors of SMAs, such as pseudoelasticity and shape memory effect.

Transform stress-strain curve into kinematic hardening curve
The model described in this paper is established based on the temperature-stress (σ − T) phase diagram, and the material parameters can be obtained by standard uniaxial tensile and compression tests. Therefore, the parameters of the phenomenological model can be measured by simple mechanical tests. The (σ − T) phase diagram contains the key information at the beginning and the completion of the SMA phase transformation. This key information has a strong inherent mechanical relation with the stress-strain curve, representing the constitutive properties of the SMA. Therefore, the inherent relation between the σ − T phase diagram and the stress-strain (σ − ε) curve will be discussed in this Section. How to transform the stress-strain curve of the SMA material into a quasikinematic hardening curve is also presented.
The significance of the typical stress-temperature phase diagram of a SMA (figure 1) is described briefly. The metal phases of SMAs can be divided into twinned martensite (M t ), detwinned martensite (M d ) and austenite (A), which are related to temperature and stress. The σ − T phase diagram mainly includes three types of parameters: the start and the finish of the complete critical transformation temperature for martensite and austenite (T M s , T M f , T A s , T A f ), the forward and the reverse critical transformation stresses for martensite and austenite (σ M s , σ M f , σ A s , σ A f ), and the slope of the critical phase transformation stresses (C M and C A ). R is the radius of the elastic region in the σ − T phase diagram. The term τ m is the so-called Maxwell stress.
SMAs display different mechanical behaviors at low and high temperatures. SMAs behaves in a pseudoelastic manner (PE) at temperatures higher than T A f (see figure 2). When the temperature is lower than T M f , the SMA behaves under shape memory effect (SME). There is a strong internal relation between the σ − T phase diagram and the σ − ε curve. The σ − ε curves of SMAs at different temperatures can be determined based on the σ − T phase diagram, the elastic modulus of the material at martensite and austenite states, and the maximum phase transformation strain. The internal mechanical relations between the phase diagram and stress-strain curves are discussed in appendix. The SMA phenomenological model described in this paper is derived by referring to the classical plastic mechanics theory. Therefore, the stress-strain curve of the SMA material is transformed into a quasi-kinematic hardening curve. Taking the pseudoelastic behavior of a shape memory alloy as an example (figure 3), the coordinate point (ε τm el , τ m ) of the original stress-strain curve is shifted to the origin to obtain a stress-strain curve which is similar to the one of a kinematic hardening model. ε τm el is the elastic strain corresponding to the Maxwell stress τ m under elastic loading The term ε L tr represents   the maximum phase transformation strain. We are calling this model as the quasi-kinematic hardening one. The red dashed line in the figure is the kinematic hardening track (typical back stress). The amplitudes of the back stress in the forward and the reverse transformation process can be expressed as: During the shape memory state, the stress τ m is equal to zero. In this state, the amplitude of the back stress in the forward transformation process can also be expressed as in equation (1). The back stress amplitudes of the forward and the reverse transformations are denoted as α amp . To ensure that the forward and the reverse back stress are consistent, we assume that: Hence, the maximum back stress amplitude is . Based on the above transformation process following a kinematic hardening model, the function between the stress τ m and the temperature is where the constant β is defined as: represents the positive function, which is defined as (x) + = (x + |x|)/ 2. The Maxwell stress is 0 when T * = T. At that particular instant of time, the stress-strain curve of SMA is the one shown in figure 2(a). The phase transformation of the SMA is related to the historical loading path. Moreover, the proposed model can record the volume fractions of twinned martensite, detwinned martensite and austenite during the loading history. Hence, only twinned martensite and austenite phases undergo forward phase transformation. No reverse phase transformation occurs at T * = T when the load is applied to SMA. The relationship between the radius of the elastic region R and temperature of the quasi-kinematic hardening model (figure 4) is

One-dimensional odd polynomial phenomenological model
In this Section, the driving force of the SMA phase transformation, composed of the stress and the back stress, is derived according to the Helmholtz free energy. An odd polynomial function is constructed to describe the phase transformation process and obtain its tangential modulus. Based on the tangential modulus and the phase transformation yield function, it can be further derived the evolution law of the transformation strain and back stress. Finally, the tangent and thermal stiffnesses are deduced by the consistency condition of the yield function and help to improving the convergence of the numerical solution.

Constitutive relations
The local state method is used here to derive the constitutive equation of SMAs based on the Clausius-Duhem inequality. This use of the method first requires the value of the specific free energy (Helmholtz free energy or Gibbs free energy) at any point of the state of the material, which can be represented by the state variable. Here, the phenomenological model of an SMA is derived in terms of Helmholtz free energy. The selected state variables are the total strain ε, the test temperature T, the phase transformation strain ε tr and the phase transformation hardening parameter E h . Among them, ε and T are observable variables that can be directly measured from mechanical tests. The parameters ε tr and E h are internal variables. The term ε el is the elastic strain. The constitutive model in this paper does not consider the effect of the plastic strain during phase transformation. In addition, the thermal strain caused by the change of temperature is relatively small compared with the phase transformation strain. For the sake of simplicity, the influence of the thermal strain is ignored in this model. Therefore, the total strain of a SMA can be decomposed into an elastic and a phase change part: From a physical point of view, the Helmholtz free energy of a representative unit is composed of the elastic potential energy Ψ el , the thermal energy Ψ T and the phase transformation energy Ψ tr . Hence, the Helmholtz free energy of a SMA can be written as where E is the elastic modulus of the NiTi SMA. ρ is the density of NiTi SMA. The temperature T 0 is the equilibrium temperature at which the Helmholtz free energy of the austenite and martensite phases are the same. The term c v is the specific heat at constant volume. The other terms of the Helmholtz free energy are: The equations include the terms E (ε tr ) = E A + εtr εL ∆E and ∆E = E M − E A . The term f (ε tr ) is the phase transformation stress during the phase transform. The complete differential equation of the Helmholtz free energy function is According to the Clausius-Duhem inequality expressed in terms of Helmholtz free energy (σ : dε − ρ where q is the heat flux and ∇T is the temperature gradient. The above inequality holds for any values of ε, ε tr and T. Furthermore, it is assumed that the SMA material is strongly dissipative, so the stress and entropy are respectively: Therefore, the phase transformation and the heat dissipation potential energy are respectively: Also, according to equation (12), the driving force of the SMA phase transformation is: The right end of equation (14) is involves the use of the stress minus the term related to the temperature and phase transformation strain. By comparing the kinematic hardening model with figure 3, α = τ m + f (ε tr ) can be regarded as the back stress of the SMAs. Thus, the constitutive relations of the SMA material are established as:

Odd polynomial function during phase transformation
The stress-strain curves during phase transformation are like odd polynomial functions (for example the segment BC of the curve in figure 10). Therefore, we introduce an odd polynomial function to describe the phase transformation process, and its form is as follows: where m = 3, 5, 7, · · · . The strain ε L is the maximum strain generated during the phase transformation stage, which is equal to the maximum phase transformation strain plus the elastic strain during the phase transformation stage: The coefficients A, B and C in equation (16) need to be determined. The coefficient B can be identified from the slope of the midpoint of the stress-strain curve at the phase transformation stage. The coefficient C can be set as half of the maximum back stress, C = α amp /2. In addition, the coefficient A can be obtained combining the starting point Therefore, the tangential modulus of the phase transformation stage can be obtained by taking the derivative of equation (16) versus the strain ε: Let us take the shape memory effect SME as an example. The elastic modulus of the SMA remains unchanged during the phase transformation stage. The stress-strain curves and its tangent modulus change characteristics during the phase transformation (figure 5). When the exponential of an odd polynomial function is large, the middle segment of the curve is flatter. Moreover, the minimum tangential modulus of the odd polynomial function is equal to the coefficient B. The layout of the stress-strain curve described by equation (16) is intuitive and simple. At first, only the start and the finish points of the phase transformation stage and the slope of the intermediate point need to be determined. Then different indices are selected according to the flatness of the stress-strain curve of the SMA phase transformation. All these parameters can be obtained by standard uniaxial tension and compression tests. The resulting odd polynomial function provides a good fit to the stress-strain curve during transformation.

Evolution of the transformation strain and back stress
The evolution law of the SMA constitutive model is derived by referring to the kinematic hardening model of plastic mechanics theory. The back stress and plastic strain need to be identified when the kinematic hardening model is applied. To obtain the evolution law of the phase transformation strain, the phase transformation yield function of the SMA should be established first. The SMA model used in this paper is related to one-dimensional structures and the model only involves uniaxial tension and compression, which can be represented by positive and negative signs. Therefore, the yield function of the phase transformation can be assumed to follow the form: The physical meaning of equation (19) is that when the stress exceeds the interval with the radius R of the center point as τ m + α, the SMA will undergo phase transformation to pull the stress back to the yield surface. For the case of associated plasticity, the evolution law of the transformation strain is expressed as: The evolution law of the back stress requires the hardening modulus of the phase transformation process. However, the explicit expression of the hardening modulus E h cannot be obtained directly by using the stress-strain relation of the phase transformation mentioned above. Hence, the function describing the hardening modulus needs to be deduced first according to the mechanics principles. The total differential form of equation (10) is: where dE = ( dε tr /ε L tr ) ∆E. In addition, the stress increment can also be obtained by the tangential modulus and the strain increment during the transformation stage: By combining equations (21) and (22) and the chain rule, the following relationship of the hardening modulus is obtained The first term on the right end of equation (23) is the classical expression of the hardening modulus used in plastic mechanics theory when the elastic modulus is constant. The second term on the right is the correction caused by the difference of the elastic modulus between martensite and austenite during phase transformation. Obviously, ∆E is equal to zero when the SMA is converted from twinned martensite to detwinned martensite. In this case, the second term at the right end is zero, and the above equation degenerates into the classical expression of the hardening modulus. Thus, when the phase transformation occurs, the increment of back stress is:

Consistent tangent stiffness
The tangent and thermal stiffnesses of a SMA can be deduced by the total differential form of equation (10) and the consistency condition of the yield function (19). The tangential stiffness here is the consistent tangent stiffness matrix adopted in standard implicit commercial finite element codes like ANSYS or ABAQUS. In the process of phase transformation, the state variables always remain on the yield surface, and the yield function is always equal to zero. In other words, the value of the yield function does not change with the state variables, i.e. the yield function of phase transformation meets the consistency condition: where Equation (21) is arranged and expressed as: where Q = E − ∆Eε el ε L tr . Substituting equation (27) into equation (25), the expression of transformation strain increment is: By combining equations (25) and (28), the tangent and thermal stiffnesses of a SMA can be obtained as: (29) and (30) can be simplified as:

Numerical integration algorithm framework
In this paper, the semi-implicit midpoint integration algorithm is used to solve the proposed SMA phenomenological model. The constitutive model here is based on the principles of the kinematic hardening model, so its numerical integration algorithm takes inspiration from the schemes adopted in plastic kinematic hardening. At the beginning of the current increment step, the known state variables are σ n , α n , ε n el , ε n tr , ∆ε, ∆T. Firstly, it is necessary to solve the plastic multiplier ∆λ and then solve the increments of each variable ∆ε tr , ∆α, ∆σ and obtain each state variable under the current increment step. The numerical integration framework of the 1D SMA phenomenological model ( figure 6) is performed as follows: (a) Calculate the trail stress: (b) Substitute equation (32) into equation (25) to obtain the value of the yield function: 1. If F ( σ trail,n+1 , ε n tr , T n ) ⩽ 0, the SMA is in an elastic state. Therefore, the trail stress is the actual stress, and each state variable is 2. F ( σ trail,n+1 , ε n tr , T n ) > 0, the SMA undergoes phase transformation. It is therefore necessary to solve the transformation strain increment ∆ε tr and then carry out the stress correction. (c) If the value of the yield function is positive, the stress correction will be carried out. It is not necessary to consider the directions of the stress or of the transformation strain in the 1D model. Hence, the actual stress is: Set The increment of back stress can be expressed as: where: Therefore, the back stress can be obtained as: In addition Substituting equations (35)-(42) into equation (19), the nonlinear equation about ∆ε tr can be obtained as follows: The Newton-Rapson method is used to solve the nonlinear equation (43) iteratively to obtain the transformation strain increment ∆ε tr . The key calculation formulas of the iterative process are described as: (d) Substitute ∆ε tr back to equations (34), (35), (37), and (41) to obtain ε n+1 tr , σ n+1 , ∆α, α n+1 . (e) Finally, obtain the tangent stiffness and thermal stiffness of SMA as: where

Numerical results and discussion
In this section, the accuracy of the polynomial phenomenological model is verified for the pseudoelasticity and the shape memory effect in shape memory alloys. The research object here is the SMA wire. The one end of the SMA wire is fixed and the tension is applied to the other end of it. The SMA wire is heated or cooled through a temperature control box to control its temperature. Firstly, the stress-strain curves of the proposed model at different temperatures are simulated and compared with the predictions provided by the Zaki and Moumni model and the experimental data of the Zaki and Moumni's test [33], in order to verify the pseudoelastic manner of the polynomial phenomenological model. The martensite volume fraction and transformation strain of the SMA are also analyzed during loading and unloading process at different temperatures. Secondly, the SMA is heated up and then cooled down to simulate the shape memory effect of the proposed model. The mechanical characteristics of the SMA under complex loading paths are predicted by the proposed model at last.

Pseudoelasticity (PE)
The stress-strain curves obtained by Zaki and Moumni's tensile tests [33] (loading and unloading) of an SMA wire at different temperatures are used to verify the validity of the proposed model. The mechanical material parameters are shown in table 1. The integral solver of the polynomial phenomenological model has been compiled and implemented on a Matlab platform. The program is only solved for a single integral point. The mechanical material parameters in table 1 are inputed into the integral solver. The loading conditions are as follows: at different temperatures, the integral point of the SMA is first loaded until the forward phase transformation is completed, and then the force load is unloaded to zero. The resulting stress-strain curves are shown in figure 7. Note that the exponent of the polynomial function of the phenomenological model is set to 7.
As shown in figure 7, the phenomenological model proposed in this paper can reflect more accurately the actual mechanical properties of the SMA, compared with the predictions provided by the Zaki and Moumni model. However, there are some differences between the calculated results of the proposed constitutive model and the experimental results. The reason behind these discrepancies is that the polynomial phenomenological model is based on the σ − T phase diagram, which simplifies the actual situation to some extent. The σ − T phase diagram assumes a linear relationship between the critical stress and temperature. In addition, the critical stress of a twinned martensite transforming into detwinned martensite is constant, but the critical stress and temperature are not actually linear.
The numerical prediction of the critical transformation stress is lower than the experimental data at 30 • C. The reason is that the transformation forms R phase in SMA wire used in Zaki and Moumni's test [33], a phenomenon affect the transformation but is not accounted in the proposed model. When the SMA wire is stretched at this temperature, the R phase of the SMA is generated first, and then it becomes twinned martensite. The critical stress of the SMA wire is increased due to this physical process. Therefore, the experimental data are larger than the predictions provided by the numerical model. When the temperature is higher than 30.0 • C, the numerical calculation of the critical stress has a small offset around the experimental results. Because the proposed model assumes that the critical stress is proportional to the temperature, i.e. the stress influence coefficients C M and C A of the martensite and austenite transformations are constant. However, these coefficients will have uncertain and small changes with the change of temperature, which will lead to the discrepancy of the numerical prediction around the experimental results. On the whole, the stress-strain curves obtained by the polynomial phenomenological model are however in good agreement with the experimental data obtained by Zaki and Moumni above 40 • C; this is a confirmation of the overall validity of the pseudoelastic modeling approach for the SMA.
The stress-strain curves shown in figure 8 describe the behavior of the model with the temperature. When the ambient temperature is lower than T A s , only the forward phase transformation occurs during the loading and unloading process, but no inverse phase transformation happens. Twinned martensite or the mixture of twinned martensite and austenite is transformed into detwinned martensite. The transformation strain generated by the loading cannot be recovered during unloading. When the temperature is higher than T A f , the forward phase transformation occurs during the loading process, and the austenite changes to detwinned martensite. During unloading, the inverse phase transformation occurs. The detwinned martensite changes back to austenite and the transformation strain disappears completely. However, when the ambient temperature is between T A s and T A f , such as at 37.0 • C, a forward phase transformation occurs during the loading process and only part of the transformation strain undergoes an inverse phase transformation during the unloading process. Therefore, residual transformation strain is present after unloading. Figure 7 also shows the critical stresses at different temperatures during the forward (inverse) phase transformation. Each critical stress is linearly dependent versus the temperature. Figure 8 shows that the proposed model allows to predict the mechanical properties of SMA actuators within a broad range of temperatures. In addition, due to the close relationship existing between the proposed model and the phase diagram, one can develop from this simulation approach a  [33].   higher fidelity version for practical engineering applications by simply determining the key parameters of the phase diagram. Hence, the proposed model could be more convenient to use for the general design and evaluation of SMA actuators. The proposed phenomenological model can describe the physical characteristics of a SMA in different phase transformations. Figure 9 shows the variation of volume fraction of twinned and detwinned martensite during loading and unloading at three typical temperatures (10 • C: twinned martensite, 20 • C: twinned martensite and austenite mixture, and 60 • C: austenite). Figure 8 also shows the change of the transformation strain and of the total strain at the corresponding temperature.

Constitutive parameter
As shown in figure 9(a), when the ambient temperature is 10 • C (which is lower than T M f −291.0 K), only twinned martensitic phase exists in the SMA. The twinned martensite phase changes to detwinned martensite during loading. The unloading process is elastic unloading of the twinned martensite without any phase transformation. At the beginning of the phase transformation, the volume fraction of the twinned martensite is ξ t = 1.0, and the volume fraction of the detwinned martensite is ξ dt = 0. The twinned martensite begins to transform into detwinned martensite with increasing load. The volume fraction of the twinned martensite decreases, while the volume fraction of the detwinned martensite increases. Correspondingly, the transformation strain shown in figure 9(d) begins to increase. Before the loading step reaches 0.5 (i.e. the stress corresponding to the load almost reaches 275.0 MPa.), the phase transformation is completed and at this point, ξ t = 0.0, ξ dt = 1.0, while the phase change strain reaches 4%. During the phase transformation, the increment of the volume fraction of the detwinned martensite is equal to the volume fraction decrease of the twinned martensite, dξ dt = −dξ t , i.e. the dissociated twinned martensite is transformed into twinned martensite. After that point, the load continues to increase, the volume fraction of twinned and detwinned martensite do not change and the phase change strain remains at 0.04. So does the subsequent unloading process, which is only related to the elastic unloading process of the twinned martensite.
It must be noticed that the change of volume fraction of twinned and detwinned martensite is not linear with the change of load during the phase transformation. At the start and finish of the phase transformation, the volume fraction of twinned and detwinned martensite changes relatively slowly. The volume fraction changes rapidly in the middle of the phase transformation. The trend of the variation of the martensite volume fraction shows the characteristics of slowly increasing at first and then rapidly increasing, followed by a soft increase. This is because the stresses and strains show an odd polynomial function relationship in the phase transformation. At the beginning of phase transformation, the stress increases quickly and the total strain increases slowly, making the elastic strain the predominant component of the strain. However, the increase of the transformation strain and the martensite volume fraction is relatively small. The same is true when the phase transformation is finished. In the middle of the phase transformation, the stress increases very little while the total strain increases rapidly. The result is a small proportion of elastic strain and a large proportion of transformation strain. Therefore, the volume fraction of the twinned and detwinned martensite increases rapidly. The change characteristics of the twinned and detwinned martensite volume fractions reasonably reflect the physical characteristics during the phase transformation process.
When the load is zero and the ambient temperature is 20 • C (i.e. between T M s and T M f ), the SMA wires are made of a mixture between twinned martensite and austenite. Figure 9(b) shows that the volume fraction of the twinned martensite is about 0.43 without loading; the volume fraction of the twinned martensite is 0, the austenite volume fraction is about 0.57 and the martensite volume fraction is close to 0.43. The mixture begins to transform into detwinned martensite when the load exceeds a critical stress. Before the load step reaches 0.5, the phase transformation is completed. At this time, ξ t = 0, ξ dt = 1.0, and the twinned martensite and austenite completely become detwinned martensite. The phenomenological model proposed in this work assumes that twinned martensite and austenite transform into detwinned martensite in equal proportions during the phase transformation. When the phase transformation is finished, the transformation strain is equal to 4.0%. Thereafter, the transformation strain remains unchanged. Figure 9(d) also shows that the critical strain at the phase transformation start at 10 • C is larger than the analogous one at 20 • C. The reason is that the SMA wire is pure twinned martensite phase at 10 • C, but the mixture between twinned martensite and austenite phase is present at 20 • C. The elastic modulus of the former is less than the one of the latter. Therefore, the critical strain of the former is larger than the SMA with twinned martensite and austenite when the critical stresses are equal.
When the ambient temperature is higher than 60 • C (i.e. larger than T A f ), only austenite exists within the SMA wire (figure 9(c)). At the initial loading time, ξ t = ξ dt = 0 and the austenite volume fraction is 1.0. When the load exceeds a certain value, the austenite starts to transform into twinned martensite. Before the load step reaches 1.0, the phase transformation is finished and at this point, ξ dt = 1.0, ξ t = 0, the austenite volume fraction is 0 and the transformation strain is 0.04. The transformation process is only a forward transformation from austenite to twinned martensite. From the completion of the phase transformation to the peak loading, it is the elastic loading process of the twinned martensite. Load step 2 represents the unloading process. Elastic unloading of the twinned martensite occurs first. Then, the twinned martensite transforms into austenite when the stress is lower than the initial critical one. This process is a reverse phase transformation. Before the loading step reaches 1.75, the reverse phase transformation of the SMA wire is completed. The austenite volume fraction is 1.0, ξ dt = 0, ξ t = 0 and the phase transformation strain is 0. The following process is represented by the elastic unloading during austenite. When the stress drops to 0, the total strain of the SMA wire is 0; during loading and unloading, only the forward and reverse phase transformations between austenite and twinned martensite occur in the SMA wires.
As shown in figure 9, the proposed model can track the change of the volume fraction of each material phase (twinned and detwinned martensite, and austenite) during the loading history. In other words, the proposed model can represent the dependence of the phase states of the SMA during its loading history. This can help understand the internal changes of the shape memory alloy and provide guidance for the design and optimization of actuators made with this smart material.

Shape memory effect (SME)
To assess how the proposed model represents the shape memory effect of the constitutive model, a load corresponding to 450.0 MPa was applied at 283.15 K, which is a lower temperature than T M f . Then, keeping the load constant, the temperature was increased to 370.0 K and finally lowered to 283.15 K.
As shown in the stress-strain curves of figure 10, during the initial loading stage, the A-B segment is the elastic loading of the twinned martensite. When the stress exceeds point B, the SMA begins to undergo forward phase transformation until the finish of the phase transformation (point C). The C-D segment is the elastic loading of the twinned martensite. The slopes of the A-B and C-D points are equal because the elastic modulus of the twinned and detwinned martensite are equal. After point D, by keeping the load unchanged and the temperature rising to 370.0 K, the SMA undergoes a reverse phase transformation starting from point E and the total strain decreases accordingly-i.e. the point E moves towards point F.  This is because the twinned martensite is transformed into austenite during the reverse phase transformation, and the transformation strain decreases to 0. Then, the forward phase transformation occurs when the SMA wire is cooled to 283.15 K. The austenite is transformed into twinned martensite again, and the point F goes to point H through point G. During this process, transformation strain and total strain increase. Figures 11 and 12 clearly show the shape memory effect of the SMA wire when the temperature is modified. The total strain corresponding to point D is 0.061, the transformation strain is 0.04, and the elastic strain is 0.021. The temperature rises from 283.15 K. The strain remains constant, and no phase transformation occurs until the temperature rises to 347.1 K. When the ambient temperature is higher than 347.1 K, the reverse phase transformation starts. The total strain and transformation strain decrease slowly at first, then rapidly, and finally again slowly. When the temperature reaches 351.1 K, the total and elastic strains are 1.1%, and the transformation strain is 0. Then the temperature is increased to 370.0 K and the strain remains unchanged. During cooling, the strain does not change until the temperature drops to 332.1 K. After the temperature is lower than 332.1 K, the forward phase transformation starts and the total and phase transformation strains increase again until 325.6 K, when the forward phase transformation is completed. Now we analyze the mechanical characteristics of the SMA under complex loading paths as predicted by the model developed in this work. The loading paths 1 and 2 in figure 13(a) show that the SMA is loaded to 450.0 MPa at the initial temperature of 283.15 K (twinned martensite) and 323.15 K (austenite), respectively; the load is then kept constant and the temperature is increased up to 370.0 K; the SMA is finally unloaded to 0 MPa stress at the same temperature.  forward, and the critical transformation stress of the twinned martensite is lower than that of the austenite. Both loading paths lead to a complete transformation into twinned martensite after loading to 450.0 MPa. The reverse phase transformations of loading paths 1 and 2 in the heating process are also coincident, i.e. they reproduce a shape memory effect. In this process, the twinned martensite is reversed into austenite. The slope of the G(g)-H(h) segment is equal to that of the a-b segment (austenite elastic modulus) and larger than that of the A-B segment (martensite elastic modulus). After the SMA wire subjected to loading paths 1 and 2 is completely unloaded, the total strain also becomes 0.
The loading paths 3 and 4 in figure 13(b) are slightly different from the loading paths 1 and 2. After the temperature rises to 370.0 K, the load is not unloaded directly, but the SMA is cooled to the initial temperature (283.15 K or 323.15 K), and then the load is reduced to 0. During the cooling process, loading paths 3 and 4 also experience a forward phase transformation process. The forward phase transformation process of the two loading paths coincides, which is H(h)-I(i) segment. In the loading path 3, the temperature is reduced to 283.15 K, and then SMA is unloaded. The unloading section is J-K, which is the elastic unloading of twinned martensite. After complete unloading, the residual transformation strain of SMA is 0.04. In the loading path 4, the temperature is lowered to 323.15 K and then SMA is unloaded. The unloading stage is j-k-l-m, which includes elastic unloading of twinned martensite, reverse phase transformation and austenite elastic unloading. After complete unloading, the SMA returns to its original length and the total strain is 0. Note that loading paths 3 and 4 exhibit two typical mechanical behaviors of SMA: pseudoelasticity and shape memory effect.
According to figure 13 and the analysis above, the proposed model appears suitable for describing a variety of complex thermal-mechanical loading conditions. The proposed model can effectively predict the mechanical properties of SMAs during their loading history and characterize the pseudoelasticity and shape memory effect of SMAs themselves. This model is therefore expected to identify the restoring force and study the mechanical properties of SMA actuators under complex loading conditions.

Conclusion
This paper has first introduced an internal transformation relationship between temperature-stress phase diagram and stress-strain curve of shape memory alloy materials. Combined with temperature-stress phase diagram, the stressstrain curve of a SMA is described by using odd polynomial functions, which make the proposed model more intuitive and simpler to mimic the mechanical behavior of SMA. In this way, the mechanical parameters of the proposed model can be measured directly by simple uniaxial tensile tests. A 1D thermodynamic phenomenological model of the SMA is then derived and established by referring to the kinematic hardening model from classical plasticity theory. Due to the difference between elastic moduli at martensite and austenite, the elastic modulus of a SMA will change during the phase transformation. Hence, the SMA constitutive models established in the past are approximated by the martensite volume fraction to represent the consistent tangent stiffness matrix or the Jacobi matrix of a set of nonlinear equations. In this paper, the explicit expression of the equivalent hardening modulus has been derived, and its expression consists of the classical hardening modulus term plus a correction term due to the change of the elastic modulus. Furthermore, the exact analytical formula of the consistent tangent stiffness matrix of the SMA phase transformation has been derived in this paper and helps to improving the convergence of the numerical solution. This makes the tangent stiffness of classical plasticity theory suitable for more general cases.
A semi-implicit midpoint integration algorithm is used to establish the numerical integration framework of the proposed model. By comparing the numerical results with the experimental data and the model of Zaki and Moumni, it is proved that the proposed polynomial phenomenological model can accurately describe the pseudoelasticity of shape memory alloy wires. The numerical results also prove that the proposed model can reliably describe shape memory effects in SMAs. Moreover, the phenomenological model proposed in this work has a strong scalability. For example, if the plastic strain is considered, the phenomenological model can be directly extended to the mixed hardening model of classical plastic theory.
There are still currently some limitations in the proposed phenomenological model. Further work needs to be done to make the model reproduce the 3D pseudo-elasticity and shape memory effect. An improved model should be also able to mimic the cyclic degradations present in pseudo-elasticity and shape memory effects under the cyclic thermo-mechanical loading conditions, such as the so-called transformation ratchetting and re-orientation ratchetting identified by Kang et al [66][67][68][69].

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).