Numerical Investigation of Preload Process of Bolted Joint with Superelastic Shape Memory Alloy

A phenomenological constitutive model is developed to describe the uniaxial transformation ratcheting behaviors of the superelastic shape memory alloy (SMA) by employing a cosine–type phase transformation equation with the initial martensite evolution coefficient that can capture the feature of the predictive residual martensite accumulation evolution and the nonlinear hysteresis loop on a finite element (FE) analysis framework. The effect of the applied loading level on transformation ratcheting is considered in the proposed model. The evolutions of transformation ratcheting and transformation stresses are constructed as the function of the accumulated residual martensite volume fraction. The FE implementation of the proposed model is carried out for the numerical analysis of transformation ratcheting of the SMA bar element. The integration algorithm and the expression of consistent tangent modulus are deduced in a new form for the forward and reverse transformation. The numerical results are compared with those of existing models; experimental results show the validity of the proposed model and its FE implementation in transformation ratcheting. Finally, a FE modeling is established for a repeated preload analysis of SMA bolted joint.


Introduction
SMA has been widely applied in Micro-Electro-Mechanical System (MEMS), actuators, biomedical devices, organ transplantation, transportation, aerospace, and civil engineering due to its well-known superelasticity, shape memory effect, excellent biocompatibility, and wear resistance [1][2][3][4][5][6].Over the last two decades, with a deeper understanding of thermo-mechanical coupling behavior of SMA, more and more scholars have been paying attention to the response of thermo-mechanical coupling behavior of the alloy undergoing cyclic loading.The cyclic deformation behavior of SMA, as described in earlier research, was followed with interested by the research group led by Miyazaki [7,8].The group's research concluded that the residual irrecoverable deformation of SMA would increase gradually with the increase of cycle times under mechanical loading cycles.The phase transformation stress would also decrease with the increase of cycle numbers.Strnadel et al. [9,10] studied the effect of alloy elements and components on mechanical cyclic deformation behavior of superelastic NiTi alloy by experiments, and revealed the inhibition effect of a high nickel content on cyclic residual deformation.Nemat-Nasser also observed the deformation features of superelastic NiTi alloy under cyclic loadings through experiments [11].It was revealed that the increment of residual martensite strain and the dissipation energy decreases with the increase of cycle numbers, and that both tend to be stable following several cycles.
The superelasticity of NiTi SMA gradually deteriorates during the process of cyclic deformation.This phenomenon could mainly be expressed as four aspects: the residual martensite strain accumulation, the gradual decrease of transformation stresses, the gradual change of phase transformation modulus, and the gradual decrease of dissipation energy with the cyclic loadings [8,[12][13][14][15][16][17][18][19].The research done by [20][21][22][23][24][25][26] shows that the mechanical behavior of the alloy is strongly dependent on the loading rate.A non-monotonic relationship can be found between the dissipation energy and the loading rate.The mechanical dissipation and the latent heat of phase transformation cause heat generation in the alloy.Kan et al. [27] also studied experimentally the cyclic deformation behavior of the untrained superelastic NiTi SMA at different loading rates.The results show that all of the phase transformation modulus, the starting stress of martensitic transformation, the hysteresis loop, and the residual strain of NiTi SMA are strongly dependent on the loading rate.
Based on existing experimental studies, in recent years, a series of phenomenological constitutive models have been built by different scholars [28][29][30][31] to describe the thermo-mechanical coupling deformation characteristics of SMA.As a further work, the FE implementation of the proposed model was carried out to calculate the accurate stress-strain responses for SMA-based devices.In the generalized plasticity frame, the return-mapping algorithm was utilized for analysis of the superelastic NiTi alloy [32,33].However, it was not yet possible to describe the deterioration of the superelasticity and shape memory effect of untrained materials during cyclic deformation.In view of the shortcomings of the above work, some scholars have further developed the above constitutive model.Based on the existing experimental observations with cyclic deformation, the corresponding cyclic constitutive model was established by Lagoudas et al. [34] to describe the accumulation of interfacial defects and residual martensite in cyclic deformation.Though Lagoudas's model shows good capability in the prediction of the nonlinear hysteresis loop shapes with a lot of material parameters, that model can only give a reasonable description for the cyclic deformation characteristic of a specific loading conditions, instead of predicting the dependence on the applied loading level.Kan and Kang et al. [35,36] made a further expansion of the cyclic constitutive model based on the experimental observation by Kang et al. [37] that is able to reasonably describe the applied loading level related to the superelastic deterioration phenomenon.Furthermore, the Kan-Kang's model accounts for the asymmetric effect of tension and compression of the material on the phase transformation ratcheting behavior and the effect of applied loading level with just a small number of material parameters.Nevertheless, the linear phase transformation hardening rule employed by Kan-Kang's model cannot predict the nonlinear hysteresis loop shape well.
This paper thus develops a phenomenological constitutive model to describe the uniaxial transformation ratcheting behaviors of superelastic SMA.A cosine-type phase transformation equation with the initial martensite evolution coefficient that can capture the feature of the predictive residual martensite accumulation evolution and the nonlinear hysteresis loop is employed in the proposed model on a FE analysis framework.The return-mapping algorithm and the consistent tangent modulus are deduced in a new form for the phase transformation.The validity of the proposed model and its FE implementation in transformation ratcheting is finally examined by a comparison between the proposed model and the existing model and experimental results.Finally, a numerical example is given to analyze the repeated preload process of the SMA bolted joint.

Constitutive Modeling for NiTi SMA under Cyclic Loading
In the last decade, many researchers focused on constitutive modeling to describe the cyclic deformation of NiTi SMA.The phenomenological constitutive model represents a good candidate to be integrated into the structure computational methods, such as the FE method, to predict the cyclic deformation of the SMA structure.However, these models with FE implementation seem to be Metals 2018, 8, 730 3 of 24 ill-suited to efficiently analyzing the complex structure due to relatively low computation efficiency and nonlinearity.
Therefore, in this work, a one-dimensional constitutive model for phase transformation ratcheting of superelastic SMA is proposed based on the Brinson-model, and then implemented into the FE model of a one-dimensional bar element to describe the cyclic deformation of some mechanical structures, such as SMA bolt joints, SMA washers, and so on.

Constitutive Equation and Internal Variables
In the proposed model, the total strain ε can be decomposed into an elastic strain tensor ε e and an inelastic strain tensor ε in with infinitesimal strain assumption.
The residual deformation of SMA under cyclic loadings is considered to be attributed to the residual martensite deformation due to the phase transformation ratcheting.In order to characterize this deformation mechanisms, two internal variables are introduced here in Helmoltz free energy, which is assumed to be additively decomposed into elastic and inelastic parts, as follows: The internal variable ξ as the martensite volume fraction depicts the stress-induced martensite phase transformation, and is constrained by 0 ≤ ξ ≤ 1. δ characterizes the accumulative martensite transformation, including the accumulated martensite volume fraction δ c .
From the principle of thermodynamics, Clausius-Duhem inequality can be expressed as: . .
The incomplete phase transformation between martensite and austenite could be observed during the cyclic phase transformation [35][36][37].Furthermore, the amount of residual martensite would increase with the increasing number of loading cycles.Therefore, the total induced-martensite volume fraction ξ that represents the progressive increase of residual martensite strain is set as an internal variable and divided into two parts, i.e., reversible martensite volume fraction ξ r and an irreversible residual one ξ ir .
According to the experimental observations, the phase transformation deformation evolves with the increasing number of loading cycles, and reaches a stable value after a certain cycle [35][36][37].In order to account for this evolution process, the internal variable δ c is picked as the accumulated martensite volume fraction ξ c that represents the evolution process of some variables and material parameters with the increasing number of cycles, and is written by: where, t is a kinematic time.Based on the frame of generalized plasticity, the total inelastic strain ε in in Equation ( 1) equals to the transformation strain ε tr from stress-induced martensite phase and its reverse.
The SMA constitutive model used here is based on a model originally formulated by [28,29], which is a phenomenological macro-scale one-dimensional constitutive model, and can be written as where E s is the Young's modulus, Ω is the transformation coefficient, Θ is the thermal elastic coefficient; T is the temperature.The subscript '0' indicates the initial values.Young's modulus E s and transformation coefficient Ω are the function of the martensite volume fraction ξ, which are given as: It is assumed that the equivalent transformation strain and the recoverable martensite volume fraction are in a proportional relationship (irreversible martensitic transformation does not participate in the transformation process).Then, the following relationship can be obtained: Equation ( 5) illustrates that the total induced-martensite volume fraction ξ consists of reversible martensite volume fraction ξ r and an irreversible residual one ξ ir .Since the irreversible martensitic transformation does not take part in the phase transformation process, the stress is only the function of reversible martensite volume fraction, which is irrelevant to the residual martensite volume fraction in the phase transformation.The expression forms of ξ r are as follows: (i) Transformation to martensite phase: (ii) Transformation to austenite phase: where, σ cr s and σ cr f are the starting and finishing stresses of martensite transformation, C A and C M are the slope for the relation between critical transformation stress and temperature, A s and A f are the starting and finishing temperatures of austenite transformation, M s and M f are starting and finishing temperature of martensite transformation, ξ r f A→M is the volume fraction of martensite transformation at the end of transformation to martensite phase, and ξ ir f M→A is the residual volume fraction of martensite Metals 2018, 8, 730 5 of 24 transformation at the end of transformation to austenite phase.The parameters a A and a M are expressed as: Now, Equation ( 8) can be rewritten for the no initial values case by considering the irrecoverable feature as follows:

Evolution Law of Parameters Governed by Accumulated Martensite Volume Fraction
As described in [34][35][36], the evolution of transformation ratcheting and transformation stresses with a number of loading cycles are related to the accumulated martensite volume fraction δ c .It was found by Lagoudas et al. [34] that the evolution of peak strain and valley strain was due to transformation ratcheting, and provided an evolution equation taking the accumulated martensite volume fraction as a governing variable.However, further research by Kan et al. [35] found that the parameters of transformation ratcheting and transformation stress were associated with the loading stress level.Detailed evolution equations can be referred to from previous work [35].The evolution equations are outlined as follows, to satisfy the integrity of the content in this work.
(i) Evolution equation for residual martensite volume fraction: .
where ξ ir max = ε ir max /ε L is the maximum irreversible residual martensite volume fraction produced in the cyclic tension-unloading test with a peak stress equal to the finish stress of phase transformation from austenite to martensite, i.e., at the end of the stress plateau.ε ir max is the maximum irreversible residual martensite strain corresponding to the maximum load of stable cycle phase transformation.ε L is the uniaxial maximum phase transformation strain.The material parameter b is to govern the saturation rate of residual martensite volume fraction ξ ir .The revised function c AM (σ) is introduced in the evolution law associated with the loading stress level to determine the correlation between that level and the residual martensite fraction, and is written as: where AM is the value of c AM at the endpoint of the forward transformation.The material parameter n is to describe the nonlinear relationship between the residual martensite volume fraction and the loading stress level.It can be found according to Equation (16a) that, in the forward phase, the value of c AM increases with the increase of the loading stress level in the interval of Q AM s ≤ σ ≤ Q AM f , and reaches its maximum at the endpoint of forward transformation.
(ii) Evolution law of transformation stress Due to incomplete phase transformation during loading cycles, the superelastic NiTi SMA shows the mixture state of the austenite phase and residual martensite phase.The transformation stresses decrease with the increasing number of loading cycles.Thus, based on the experimental observations [37], the evolution equations with an exponential formulation were proposed by [35] to describe the progressive evolution of the transformation stresses with an increasing number of loading cycles from their initial values to stable ones, and are introduced here as: Metals 2018, 8, 730 6 of 24

Numerical Integration Algorithm
Based on the infinitesimal strain assumption, the increment of substep n + 1 of total strain ∆ε n+1 can be defined as the summation of the increments of elastic strain, and transformation strain, According to Equation (10), the increment formulation of transformation strain can be written as: Based on Equation ( 14), the stress in the incremental calculation can be obtained as: in which, ε n+1 = ε n + ∆ε n+1 is the strain in the current increment step, and ε n is the strain in the substep n.
In this research, the total martensite volume fraction is divided into two parts, i.e., reversible martensite volume fraction ξ r and irreversible residual one ξ ir .Hence, the increment of the reversible martensite volume fraction can be written as: According to Equations ( 7), (11), and ( 12), the phase transformation conditions can be introduced as: Substituting Equation (20) into Equation ( 22) gives: in which, σ s (ξ) = σ AM s (ξ) for forward phase transformation, and σ s (ξ) = σ MA s (ξ) for reverse phase transformation.
Solving Equation ( 23) by the Newton-Raphson method yields: in which, H tr = H f or for forward transformation, H tr = H rev for reverse transformation.The increment of total martensite volume fraction can be updated in substep k + 1 of Newton iteration by Equation ( 20) until the convergence conditions is reached, i.e., c tr | k+1 /∆ξ (n+1) | k+1 < Toler.
Derived from Equations ( 5), ( 9) and ( 15), the incremental formulation of some variables in Equation ( 24) can be obtained as: After ∆ξ (n+1) and ∆ξ r (n+1) are obtained, the increment of transformation strain and the stress and can be updated by Equations ( 19) and (20).

Solving for Incremental Stiffness Matrix
In order to construct the FE model of a one-dimensional bar element of superelastic NiTi SMA ratcheting behavior, the constitutive equation as Equation ( 14) should be expressed as the incremental formulation.The increment stiffness matrix of bar element should be further derived.Therefore, Equation ( 14) needs first to be expressed as the relationship of the change in nodal force and the elemental deformation in length.Then it is carried through a resolution of nonlinear equation by the variational method to obtain the stiffness matrix.The relationship between of force and deformation is expressed as: where, F s = σ s A s = σ s A s x L s is the force vector of bar element, σ s is the axial stress of SMA bolt bar, A s is the cross sectional area of bar element, L s 0 is the unstressed initial length of bar element, L s is the real-time length of bar element after deformation; x is the nodal coordinates of bar element, and x/L s is the nodal position vectors of bar element, and σ s = σ s A s x/L s is the vector formula of axial stress.
Since the phase transformation strain is a relatively large value, the total strain is calculated in the logarithmic formulation ε = ln L s /L s 0 to approximate the real one.The variation formula of F s can be expressed as: According to the variation formula of Equation ( 28), the following expression can be obtained as: and From Equation ( 9), it can be derived: Now, the stiffness matrix K sma U cannot be obtained yet unless the increment of the phase transition dσ/dξ r is given.The derivation process of the increment of the phase transition will be discussed as follow.
(i) Transition to martensite phase According to the model discussed in [28], the forward transition stresses can be expressed as: Equation ( 11) can be rewritten as: The differential formulation of Equation ( 34) can be further expressed as: in which,

(ii) Transition to austenite phase
For reverse transition, Equation ( 12) can be rewritten as: With the description in [28], the reverse transition stresses can be expressed as: Another expression can be written as: Therefore, Equation ( 13) can be expressed as: Therefore, the differential formulation of Equation ( 36) can be derived as: in which, Combined with Equations ( 29)-( 32), (35), and (40), the stiffness matrix K sma U of single bar element of SMA can be obtained.The return mapping algorithm is described in Algorithm 1, and calculated by MATLAB (Matlab 2015, Produced by MathWorks Company, Natick, MA, USA).

Begin
Step 1 Initialize variables, assumed ξ = 0 at the beginning Step 2 Given K sma U and ∆x, calculate

Numerical Simulation and Model Verification
In this section, the cyclic responses for ratcheting associated with the one dimensional cyclic behavior of SMA by the current FE model are verified by some typical experimental results [37] which are outlined in this section, and additional experimental observations made in the current work to describe the transformation ratcheting during the stress-controlled cyclic loading at room temperature.The typical tensile-unloading stress-strain curve of superelastic SMA is shown in Figure 1.An apparent superelastic feature of SMA is shown.However, its curve presents a slight difference from the description in the referred literature for the NiTi SMA manufactured by other companies (e.g., SMA, San Jose, CA, USA).It exhibits an apparent hardening behavior during the stress-induced martensite transformation.After unloading, a relative high residual strain (ε r = 4.5%) could be seen to remain.The high residual strain implies that there is an incomplete phase transformation from the stress-induced martensite to original austenite after unloading, which leads to some remaining amount of martensite.It was further found that the amount of remaining martensite increases progressively with the cyclic loadings.It is different from the ratcheting of ordinary metals without phase transformation; this phenomenon of accumulation deformation has been named as "phase transformation ratcheting", and is discussed in more detail in [34][35][36][37].It illustrates that the pure austenite replaced by the mixture of austenite and remained martensite in the metal after the cycle loadings.Certainly, the stresses are no longer the phase transformation stresses of the pure austenite and martensite.to remain.The high residual strain implies that there is an incomplete phase transformation from the stress-induced martensite to original austenite after unloading, which leads to some remaining amount of martensite.It was further found that the amount of remaining martensite increases progressively with the cyclic loadings.It is different from the ratcheting of ordinary metals without phase transformation; this phenomenon of accumulation deformation has been named as "phase transformation ratcheting", and is discussed in more detail in [34][35][36][37].It illustrates that the pure austenite replaced by the mixture of austenite and remained martensite in the metal after the cycle loadings.Certainly, the stresses are no longer the phase transformation stresses of the pure austenite and martensite.=   .This parameter reflects the damping feature of NiTi SMA, which is a unique property of the metals, and has been extensively applied in the engineering applications.
The capacity of the proposed model to describe the uniaxial transformation ratcheting of NiTi SMA at the temperature with pure austenite phase is firstly verified by comparing the predicted results with the corresponding experimental ones by Kang et al. [37] and the simulated ones by Kan et al. [35,36].The material parameters of the used NiTi SMA cited from Kan et al. [35].
The results obtained with various peak stresses, e.g., 450, 500, and 550 MPa are shown in Figures 2-4.It is seen from the figures that the proposed model provides reasonable predictions to the uniaxial transformation ratcheting of superelastic NiTi SMA compared with the experiments observed by Kang et al. [37], as shown in Figures 2a, 3a and 4a.The simulated results obtained from simulated results with various peak stresses are shown in Figures 2b, 3b and 4b.The figures show the effectiveness of the current FE model at predicting the uniaxial transformation ratcheting of This parameter reflects the damping feature of NiTi SMA, which is a unique property of the metals, and has been extensively applied in the engineering applications.
The capacity of the proposed model to describe the uniaxial transformation ratcheting of NiTi SMA at the temperature with pure austenite phase is firstly verified by comparing the predicted results with the corresponding experimental ones by Kang et al. [37] and the simulated ones by Kan et al. [35,36].The material parameters of the used NiTi SMA cited from Kan et al. [35].
The results obtained with various peak stresses, e.g., 450, 500, and 550 MPa are shown in Figures 2-4.It is seen from the figures that the proposed model provides reasonable predictions to the uniaxial transformation ratcheting of superelastic NiTi SMA compared with the experiments observed by Kang et al. [37], as shown in Figures 2a, 3a and 4a.The simulated results obtained from simulated results with various peak stresses are shown in Figures 2b, 3b and 4b.The figures show the effectiveness of the current FE model at predicting the uniaxial transformation ratcheting of superelastic SMA, including the hysteresis loop curve, and the predicted peak and valley strain.The peak and residual strains of superelastic SMA increase progressively during the initial several loading cycles, and then reach a stable value.The near-totally closed hysteresis loop means that the strain increment happens in the tension loading is completely recoverable under the following unloading.The closed hysteresis loop that still becomes smaller and smaller during the further cyclic loading means that the dissipation energy decreases with the increase number of loading cycles.Compared with the Kan-Kang's model and with linear hardening law, as described in [35], the introduction of the cosine-type nonlinear function in the constitutive model enables the predicted hysteresis loop with an apparent nonlinear feature that is closer to the experimental results.Compared with the Kan-Kang's model and with linear hardening law, as described in [35], the introduction of the cosine-type nonlinear function in the constitutive model enables the predicted hysteresis loop with an apparent nonlinear feature that is closer to the experimental results.
(a) (b)  Compared with the Kan-Kang's model and with linear hardening law, as described in [35], the introduction of the cosine-type nonlinear function in the constitutive model enables the predicted hysteresis loop with an apparent nonlinear feature that is closer to the experimental results.
(a) (b)   Due to the introduced power function c f AM (σ) that is associated with the applied loading level in the phase transformation ratcheting simulation, the model can reasonably predict the variation of peak and valley strain with different applied loading levels.It should be noted that the introduction of the power equation with the coefficient of n, as described by Equation (16a), can reasonably improve the prediction of peak strain and valley residual strain changing with different applied loading level, though differences exist between the simulated and experimental results, especially for larger loading levels.This indicates that the power equation introduced in the proposed model dose still does not fully reflect the highly nonlinear relationship between the transformation ratcheting and the applied loading level.To achieve better prediction results, a more complex function formulation with added more material parameters should be introduced.

Due to the introduced power function ( )
f AM c σ that is associated with the applied loading level in the phase transformation ratcheting simulation, the model can reasonably predict the variation of To evaluate the current FE model, some experiments under the uniaxial cyclic loading conditions are carried out, and then used to verify the validity of the proposed model in the following discussions.The materials used in the new experiments are the superelastic NiTi SMA bars (Ni, 55.89%, Ti, 44.11% at mass, Xi'an Saite Metal Materials Development Co., Ltd., Xi'an, China).The heat treatment process is different from the SMA used in Kang et al. [37], and then the start temperature of martensite transformation A f is 280 K (7 • C for celsius temperature scale), tested by Differential Scanning Calorimetry (DSC, Mettler-Toledo Ltd., Melbourne, Australia) with 5 • C/min for heating and cooling as show in Figure 5, lower than the test temperature 298 K, and the original phase of the alloy is the pure austenite phase.Three different kinds of the uniaxial nominal are controlled by axial load under cyclic tension-unloading with positive mean loads, including cyclic tension-unloading tests with various applied peak stresses, e.g., 325, 365, and 405 at room temperature.The number of cycles is prescribed as 50, and the stress rate is prescribed as 20 MPa/s.The material parameters used in the current FE model are determined by a trial and error method from the test data, as described in Section 2.2; the parameters are listed in Table 1 for the current FE simulation.peak and valley strain with different applied loading levels.It should be noted that the introduction of the power equation with the coefficient of n, as described by Equation (16a), can reasonably improve the prediction of peak strain and valley residual strain changing with different applied loading level, though differences exist between the simulated and experimental results, especially for larger loading levels.This indicates that the power equation introduced in the proposed model dose still does not fully reflect the highly nonlinear relationship between the transformation ratcheting and the applied loading level.To achieve better prediction results, a more complex function formulation with added more material parameters should be introduced.
To evaluate the current FE model, some experiments under the uniaxial cyclic loading conditions are carried out, and then used to verify the validity of the proposed model in the following discussions.The materials used in the new experiments are the superelastic NiTi SMA bars (Ni, 55.89%, Ti, 44.11% at mass, Xi'an Saite Metal Materials Development Co., Ltd., Xi'an, China).The heat treatment process is different from the SMA used in Kang et al. [37], and then the start temperature of martensite transformation Af is 280 K (7 °C for celsius temperature scale), tested by Differential Scanning Calorimetry (DSC, Mettler-Toledo Ltd., Melbourne, Australia) with 5 °C/min for heating and cooling as show in Figure 5, lower than the test temperature 298 K, and the original phase of the alloy is the pure austenite phase.Three different kinds of the uniaxial nominal are controlled by axial load under cyclic tension-unloading with positive mean loads, including cyclic tension-unloading tests with various applied peak stresses, e.g., 325, 365, and 405 at room temperature.The number of cycles is prescribed as 50, and the stress rate is prescribed as 20 MPa/s.The material parameters used in the current FE model are determined by a trial and error method from the test data, as described in Section 2.2; the parameters are listed in Table 1 for the current FE simulation.Figures 6a, 7a and 8a show the experimental results of the uniaxial tension-unloading, which is used to understand the basic performance of the NiTi SMA.It can be found from Figure 7 that the start and finish stresses of forward transformation show a decrease from 285 MPa to 225 MPa with the increase number of loading cycles, whereas the start and finish stresses of reverse transformation are from 345 MPa to 320 MPa, respectively.After the unloading, a small residual strain is observed; this is mainly caused by the martensite residual during phase transformation.It is seen from the Figures 6b, 7b and 8b, that the proposed model provides reasonable predictions to the uniaxial transformation ratcheting of superelastic NiTi SMA.The predicted peak and residual strains, and their evolution characteristics, are in fairly good agreement with the experimental results obtained in the cyclic tension-unloading tests by current tests.Also, the dependence of transformation ratcheting on the applied peak stress was reasonably well predicted by the model due to the employment of stress-dependent power function.

Material Parameters
From Figures 6c,d, 7c,d and 8c,d, it was found that the peak and valley accumulation strains increase with an exponential function law during the loading cycles, whereas the dissipation energy W d decreases.The proposed model also provides a reasonable prediction of the evolution of transformation stresses during the cyclic tension-unloading, as shown in Figure 7e.The relatively high degree of uniformity illustrates that the proposed model could reasonably simulate the transformation ratcheting of SMA.After a certain number of cycles, the accumulation strains, transformation stresses, and dissipation energy show apparent changes, and quickly reach a stable value.These findings also agree with the conclusions simulated by Kan-Kang's model [35].(e)  It is seen from the Figures 6b, 7b and 8b, that the proposed model provides reasonable predictions to the uniaxial transformation ratcheting of superelastic NiTi SMA.The predicted peak and residual strains, and their evolution characteristics, are in fairly good agreement with the experimental results obtained in the cyclic tension-unloading tests by current tests.Also, the dependence of transformation ratcheting on the applied peak stress was reasonably well predicted by the model due to the employment of stress-dependent power function.To further demonstrate the reasonability of the proposed constitutive model to predict the transformation ratcheting of superelastic NiTi alloy, two representative constitutive models developed by [34,35] are discussed in terms of their capability of predicting the transformation ratcheting in Figure 9, respectively.The experimental data is obtained from [10].The proposed model shows the capability of predicting the nonlinear feature of a hysteresis loop by comparing it with Kan-Kang's model with linear hardening features and an identical number of material parameters.Compared with the Lagoudas's model, the proposed model can also predict the nonlinear feature of hysteresis loop with the smaller number of material parameters.In addition, the proposed model can further give a reasonable description for the cyclic deformation characteristic on a specific applied loading level.
It should be noted that the exponential function introduced in the proposed model correlated the phase transformation ratcheting to the loading stress level still cannot fully reflect the high nonlinear correlation between this ratcheting behavior and the loading stress level (e.g., Figures 3 and 8).In order to achieve better prediction results, it is necessary to consider the introduction of more complex functional forms, and increase more material parameters.The evolution of the phase transformation ratcheting strain is the focus of this paper.Therefore, this paper adopts a simpler exponential function which can describe its general rule.
For some cases with unsatisfactory simulation results, the following are explained: (1) In the phase transformation process, the mismatch of the internal strain of austenite and martensite produces large levels of local stress between the interface of austenite and martensite.The local stress promotes the dislocation slip of austenite to reach the starting critical stress, that is, the induced plasticity due to phase transformation is produced.However, because the grain size of the selected material for experiment is relatively large, only the austenite near the interface between austenite and martensite can induce plasticity due to phase transformation at certain stress levels.
When the outer stress level reaches the starting critical stress of austenite dislocation slip, that slip of the whole grain can be motivated.Therefore, by the comparison between the experimental results of Figures 2a, 3a and 4a, it is known that the valley value strain of the first circle shows significant nonlinear growth; the same phenomenon can be observed from the experimental results of Figures 6-8.
value.These findings also agree with the conclusions simulated by Kan-Kang's model [35].
To further demonstrate the reasonability of the proposed constitutive model to predict the transformation ratcheting of superelastic NiTi alloy, two representative constitutive models developed by [34,35] are discussed in terms of their capability of predicting the transformation ratcheting in Figure 9, respectively.The experimental data is obtained from [10].The proposed model shows the capability of predicting the nonlinear feature of a hysteresis loop by comparing it with Kan-Kang's model with linear hardening features and an identical number of material parameters.Compared with the Lagoudas's model, the proposed model can also predict the nonlinear feature of hysteresis loop with the smaller number of material parameters.In addition, the proposed model can further give a reasonable description for the cyclic deformation characteristic on a specific applied loading level.It should be noted that the exponential function introduced in the proposed model correlated the phase transformation ratcheting to the loading stress level still cannot fully reflect the high nonlinear correlation between this ratcheting behavior and the loading stress level (e.g., Figures 3  and 8).In order to achieve better prediction results, it is necessary to consider the introduction of more complex functional forms, and increase more material parameters.The evolution of the phase (2) phenomenological constitutive model proposed in this paper is reasonable, as the amount of the induced plasticity is small.However, it does not consider the critical motivation conditions of the starting stress of the dislocation slip of austenite.Therefore, when the stress level of the external load reaches a certain value-for example, more than 500 MPa-a certain difference between the proposed model and the experimental value in the prediction of the hysteresis loop exists.
(3) Because of the above reasons, the prediction of ratcheting strain has a certain deviation in the initial cycles when the external load is larger, and a better prediction result is obtained.Although the proposed model has a certain numerical difference from the experimental results in the prediction of individual conditions, the main characteristics of the described phase transformation ratcheting behavior are in accordance with the experimental results.For this kind of cyclic deformation with such strong nonlinearity, the accurate description of the hysteresis loop should increase a lot of material parameters, which is not conducive to the application of the engineering for the constitutive model.

Numerical Examples for the Analyses of Preload Force of SMA Bolted Joint
The preload force of SMA bolted joint is an important factor for the service properties of bolted joint, which also plays a crucial role in protecting the static and dynamic characteristics and the locking and sealing properties of the whole SMA joint.However, in some conditions, it is inevitable that the preload force of the SMA bolt must be reloaded.This repeated process might decay the excellent superelasticity of SMA bolt due to its phase transformation ratcheting of materials, especially for large applied preload force.In this paper, based on the one-dimensional FE model of the superelastic SMA, a simplified FE model to analyze the repeated preloading process of SMA bolt is derived and used to calculate the repeated preloading process of SMA bolted joint as a numerical example.The member stiffness and the repeated preloading numbers are the influence factors considered in this research.
A preset interference method is adopted in this research to simulate the bolt preloading, which has successfully been used to simulate the preloading process bolt in [38][39][40].The preset interference amount of ∆u is built by shortening the distance of the bolt head and the nut in model, as shown in Figure 10.That means the preload force will be produced between the bearing surfaces of the bolt head or nut and the member that will be integrated into one surface after bolt preloading.

FE Modeling for the Preload Process of SMA Bolted Joint
For simplicity, the FE model is only established by a half of the whole bolted joint structure in this research due to its symmetrical characteristic on the whole structure size and load-bearing feature of bolted joint along the contact surface of the contacted members.The nodes (e.g., node numbers of 2 and 4) close to the side of bearing surface of one member are fixed at the restrained end.For the specific preloading process, the initial parameters for simulation should be inputted into model first, i.e., the elastic modulus of material, cross sectional area, and effective length of bolt bar, to calculate the tensile stiffness of bolt, and the compression stiffness of members.The FE model with bar element is setup based on the preset interference value 0 u Δ , which is the initial value of u Δ .
The target of the iterative calculation is to reach the force equilibrium conditions, i.e., for bolt and member, respectively.
For the absolute coordinate values of nodes 2 and 4 for each iterative sub-step calculation, the following deformation equation is used: in which, n is number of iterative sub-step.
The amount of interference u Δ decreases gradually with iterative calculation in progress.For the case of 0 u Δ = , the amount of penetration between the bearing surfaces of the bolt or nut and the member is zero, and the iteration is terminated.It should be noted that the target preloading force

FE Modeling for the Preload Process of SMA Bolted Joint
For simplicity, the FE model is only established by a half of the whole bolted joint structure in this research due to its symmetrical characteristic on the whole structure size and load-bearing feature of bolted joint along the contact surface of the contacted members.The nodes (e.g., node numbers of 2 and 4) close to the side of bearing surface of one member are fixed at the restrained end.For the specific preloading process, the initial parameters for simulation should be inputted into model first, i.e., the elastic modulus of material, cross sectional area, and effective length of bolt bar, to calculate the tensile stiffness of bolt, and the compression stiffness of members.The FE model with bar element is setup based on the preset interference value ∆u 0 , which is the initial value of ∆u.The target of the iterative calculation is to reach the force equilibrium conditions, i.e., F b = F m = F p0 from initial force conditions F b = F m = 0.For each sub-step of iterative calculation, a sub tension force is applied to the bolt, and the sub compression force is used on the member.The amount of deformation for each iterative sub-step is set as ∆u b and ∆u m for bolt and member, respectively.For the absolute coordinate values of nodes 2 and 4 for each iterative sub-step calculation, the following deformation equation is used: in which, n is number of iterative sub-step.The amount of interference ∆u decreases gradually with iterative calculation in progress.For the case of ∆u = 0, the amount of penetration between the bearing surfaces of the bolt or nut and the member is zero, and the iteration is terminated.It should be noted that the target preloading force F p0 should continuously adjust the amount of initial interference ∆u 0 by trial and error.The specific implementation process is shown in the flow chart of Figure 11.It should be noted that the SMA bolt element and the member element are independent of each other in FE modeling.This means that the force increments of iteration for The equilibrium equation of SMA bolt or member can be expressed using the standard FE assembly operation: where ΔF indicates the increment of nodal force vectors of bar elements of the FE model, Δu is the nodal displacement vectors, and U K is the stiffness matrices.
The FE model of a member is shown in Figure 10.It can be seen that the member elements consist of the nodes 1 and 2. The node 2 is the restrained end, and the node 1 is the force applied end of the member element.Thus, Equation ( 42) can be further decomposed into: where, 1 u , 2 u and 1 F , 2 F are the displacements and the forces of the nodes 1 and 2, respectively.Since the member is equivalent to a bar element in the FE model, the relationship between force vector and displacement vector of the node i can be written as: where, m A is the equivalent cross section area of the member element.0 m L is the initial unstressed length of member element, m E is the elastic modulus of member element, m k is the compression It should be noted that the SMA bolt element and the member element are independent of each other in FE modeling.This means that the force increments of iteration for ∆F m n and ∆F e n are applied on the SMA bolt element and the member element, respectively, in spite of ∆F m n = ∆F e n .Then, the deformations for each element are calculated independently until the condition of convergence, ∆u = 0, is satisfied.
The equilibrium equation of SMA bolt or member can be expressed using the standard FE assembly operation: where ∆F indicates the increment of nodal force vectors of bar elements of the FE model, ∆u is the nodal displacement vectors, and K U is the stiffness matrices.
The FE model of a member is shown in Figure 10.It can be seen that the member elements consist of the nodes 1 and 2. The node 2 is the restrained end, and the node 1 is the force applied end of the member element.Thus, Equation ( 42) can be further decomposed into: where, u 1 , u 2 and F 1 , F 2 are the displacements and the forces of the nodes 1 and 2, respectively.F 1 and F 2 are a pair of equilibrium forces with equal magnitude in opposite directions.K 11  u , K 12 u , K 13 u , K 21 u and K 22 u are the partitioned matrices of matrix K m U .Since the member is equivalent to a bar element in the FE model, the relationship between force vector and displacement vector of the node i can be written as: where, A m is the equivalent cross section area of the member element.L m0 is the initial unstressed length of member element, E m is the elastic modulus of member element, k m is the compression stiffness of member element.F m is the force vector of member element, x = x i − x j depicts the relative nodal position denoted by x i and x j in the global frame.i, j = 1 or 2 are the node numbers of member element.L m = x i − x j T x i − x j is the stressed length of the member element after deformed.With the first order Taylor expansion of Equation (45), it can be derived as: Since ∆F j = −∆F i , the incremental relationship between nodal force and element length can be expressed as:

Conclusions
In this work, a phenomenological constitutive modeling and its FE implementation of superelastic SMAs is carried out to simulate the transformation ratcheting behaviors of the superelastic SMA undergoing cyclic loading.The conclusions are as follows: (1) The assumed cosine-type phase transformation function considering the initial martensite evolution is verified to be a candidate for predicting the contributions of cyclic accumulation of residual martensite caused by incomplete reverse transformation, and the nonlinear behavior of hysteresis loop undergoing cyclic loading.The evolutions of transformation ratcheting and transformation stresses can be controlled by an accumulated martensite volume fraction as a internal variable.The correlation between the applied loading level and the transformation ratcheting is also able to be captured by the proposed model.
(2) A FE model derived by a return-mapping method is then used to analyze the attenuation law of the clamping force of SMA bolt experiences as numerical example.It was found that a great attenuation of the clamping force of the SMA bolt occurs during the initial several loading cycles, and tends to be stable with a very small value thereafter.The conclusions that the larger initial preload force of the bolt shows a higher residual clamping force after experiencing cycle preloading, and that the higher member stiffness would improve the reduction of preload force under repeated preloading cycles, will help the engineer to improve his design.

Conclusions
In this work, a phenomenological constitutive modeling and its FE implementation of superelastic SMAs is carried out to simulate the transformation ratcheting behaviors of the superelastic SMA undergoing cyclic loading.The conclusions are as follows: (1) The assumed cosine-type phase transformation function considering the initial martensite evolution is verified to be a candidate for predicting the contributions of cyclic accumulation of residual martensite caused by incomplete reverse transformation, and the nonlinear behavior of hysteresis loop undergoing cyclic loading.The evolutions of transformation ratcheting and transformation stresses can be controlled by an accumulated martensite volume fraction as a internal variable.The correlation between the applied loading level and the transformation ratcheting is also able to be captured by the proposed model.
(2) A FE model derived by a return-mapping method is then used to analyze the attenuation law of the clamping force of SMA bolt experiences as numerical example.It was found that a great attenuation of the clamping force of the SMA bolt occurs during the initial several loading cycles, and tends to be stable with a very small value thereafter.The conclusions that the larger initial preload force of the bolt shows a higher residual clamping force after experiencing cycle preloading, and that the higher member stiffness would improve the reduction of preload force under repeated preloading cycles, will help the engineer to improve his design.
stresses of stable phase transformation.c AM s , c AM f , c MA s and c MA f are the parameters to govern the saturated rates of the transformation stresses.

Figure 1 .
Figure 1.Typical cyclic tension-unloading stress-strain curve of materials.Some parameters should be defined prior to discussing the ratcheting deformation of superelastic SMA, as shown in Figure 1.Those parameters include the elastic modulus of austenite E A and martensite E M , the start stresses of austenite to martensite σ AM s and martensite to austenite σ MA s , and the finish stress austenite to martensite σ AM f and martensite to austenite σ AM f.The subscript '0' here indicates the first phase transformation cycle of cyclic tension-unloading.It should be noted that the parameters of elastic modulus and transformation stress are nominal variables due to the existence of residual martensite and its change during cyclic loadings.Moreover, the dissipation energy W d is defined as the area around by the stress-strain curve in each loading-unloading cycle, W d = σdε.This parameter reflects the damping feature of NiTi SMA, which is a unique property of the metals, and has been extensively applied in the engineering applications.The capacity of the proposed model to describe the uniaxial transformation ratcheting of NiTi SMA at the temperature with pure austenite phase is firstly verified by comparing the predicted results with the corresponding experimental ones by Kang et al.[37] and the simulated ones by Kan et al.[35,36].The material parameters of the used NiTi SMA cited from Kan et al.[35].The results obtained with various peak stresses, e.g., 450, 500, and 550 MPa are shown in Figures2-4.It is seen from the figures that the proposed model provides reasonable predictions

Figure 2 .Figure 3 .
Figure 2. Experiments and simulations for cyclic tension-tension with loading stress of 450 MPa: (a) experimental result; (b) simulated result by proposed model.

Figure 2 .
Figure 2. Experiments and simulations for cyclic tension-tension with loading stress of 450 MPa: (a) experimental result; (b) simulated result by proposed model.

Figure 2 .
Figure 2. Experiments and simulations for cyclic tension-tension with loading stress of 450 MPa: (a) experimental result; (b) simulated result by proposed model.

Figure 3 .
Figure 3. Experiments and simulations for cyclic tension-tension with loading stress of 500 MPa: (a) experimental result; (b) simulated result by proposed model.

Figure 3 .
Figure 3. Experiments and simulations for cyclic tension-tension with loading stress of 500 MPa: (a) experimental result; (b) simulated result by proposed model.

Figure 3 .
Figure 3. Experiments and simulations for cyclic tension-tension with loading stress of 500 MPa: (a) experimental result; (b) simulated result by proposed model.

Figure 4 .
Figure 4. Experiments and simulations for cyclic tension-tension with loading stress of 550 MPa: (a) experimental result; (b) simulated result by proposed model.

Figure 4 .
Figure 4. Experiments and simulations for cyclic tension-tension with loading stress of 550 MPa: (a) experimental result; (b) simulated result by proposed model.

Figure 6 .
Figure 6.Experiments and simulations for cyclic tension-unloading transformation ratchetting with constant loading stress of 325 MPa: (a) experimental results; (b) simulated results; (c) curves of peak and valley strains vs. cycle numbers; (d) curves of dissipation energy vs. cycle numbers.

Figure 6 .
Figure 6.Experiments and simulations for cyclic tension-unloading transformation ratchetting with constant loading stress of 325 MPa: (a) experimental results; (b) simulated results; (c) curves of peak and valley strains vs. cycle numbers; (d) curves of dissipation energy vs. cycle numbers.

Figure 6 .
Figure 6.Experiments and simulations for cyclic tension-unloading transformation ratchetting with constant loading stress of 325 MPa: (a) experimental results; (b) simulated results; (c) curves of peak and valley strains vs. cycle numbers; (d) curves of dissipation energy vs. cycle numbers.

Figure 7 .
Figure 7. Experiments and simulations for cyclic tension-unloading transformation ratchetting with constant loading stress of 365 MPa: (a) experimental results; (b) simulated results; (c) curves of peak and valley strains vs. cycle numbers; (d) curves of dissipation energy vs. cycle number; (e) transformation stresses vs. cycle numbers.

Figure 7 .
Figure 7. Experiments and simulations for cyclic tension-unloading transformation ratchetting with constant loading stress of 365 MPa: (a) experimental results; (b) simulated results; (c) curves of peak and valley strains vs. cycle numbers; (d) curves of dissipation energy vs. cycle number; (e) transformation stresses vs. cycle numbers.

Figure 7 .
Figure 7. Experiments and simulations for cyclic tension-unloading transformation ratchetting with constant loading stress of 365 MPa: (a) experimental results; (b) simulated results; (c) curves of peak and valley strains vs. cycle numbers; (d) curves of dissipation energy vs. cycle number; (e) transformation stresses vs. cycle numbers.

Figure 8 .
Figure 8. Experiments and simulations for cyclic tension-unloading transformation ratchetting with constant loading stress of 405 MPa: (a) experimental results; (b) simulated results; (c) curves of peak and valley strains vs. cycle numbers; (d) curves of dissipation energy vs. cycle number.

Figure 8 .
Figure 8. Experiments and simulations for cyclic tension-unloading transformation ratchetting with constant loading stress of 405 MPa: (a) experimental results; (b) simulated results; (c) curves of peak and valley strains vs. cycle numbers; (d) curves of dissipation energy vs. cycle number.

Figure 9 .
Figure 9. Stress-strain response of SMA to cyclic loading up to a constant value of stress: curves for the first and 50th cycles: (a) experimental results; (b) simulated results by Lagoudas's model; (c) simulated results by Kan-Kang's model; (d) simulated results by proposed model.

Figure 9 .
Figure 9. Stress-strain response of SMA to cyclic loading up to a constant value of stress: curves for the first and 50th cycles: (a) experimental results; (b) simulated results by Lagoudas's model; (c) simulated results by Kan-Kang's model; (d) simulated results by proposed model.

Figure 10 .
Figure 10.Loading process and FE modeling of bolted joint: (a) before preloaded; (b) after preloaded.
each sub-step of iterative calculation, a sub tension force is applied to the bolt, and the sub compression force is used on the member.The amount of deformation for each iterative sub-step is set as

F
should continuously adjust the amount of initial interference 0 u Δ by trial and error.The specific implementation process is shown in the flow chart of Figure 11.

Figure 10 .
Figure 10.Loading process and FE modeling of bolted joint: (a) before preloaded; (b) after preloaded.

Figure 11 .
Figure 11.Flow chart of analysis of the repeated cyclic preloading.
SMA bolt element and the member element, respectively, in spite of each element are calculated independently until the condition of convergence, 0 u Δ = , is satisfied.

1 F and 2 F
are a pair of equilibrium forces with equal magnitude in opposite directions.

Figure 11 .
Figure 11.Flow chart of analysis of the repeated cyclic preloading.

Figure 13 .
Figure 13.Clamping force reduction of bolt with increasing loading cycles under different load cases.

Figure 13 .
Figure 13.Clamping force reduction of bolt with increasing loading cycles under different load cases.

Step 3 else goto Step 5 endif Step 3
Calculate ∆ξ r Calculate σ, E s and ξ if ξ = 0, goto Step8 elseif T > A S and C A (T − A f ) < σ < C A (T − A s ) if c tr | k+1 /∆ξ (n+1) | k+1 < toler obtaining ∆ξ, and calculating dσ/dξ r , ξ r M→A , ξ c endif Step 7 Calculate stiffness matrix K sma U and force increment ∆F s , goto Step 8 Step 8 Updating status variables End The stiffness matrix of one node can be obtained by the differential calculation as: