An extended thermo-mechanically coupled algorithm for simulation of superelasticity and shape memory effect in shape memory alloys

Abstract

Thermo-mechanical coupling in shape memory alloys is a very complicated phenomenon. The heat generation/absorption during forward/reverse transformation can lead to temperature-dependent variation of its mechanical behavior in the forms of superelasticity and shape memory effect. However, unlike the usual assumption, slow loading rate cannot guarantee an isothermal process. A two-dimensional thermo-mechanically coupled algorithm is proposed based on the original model of Lagoudas to efficiently model both superelasticity and shape memory effects and the influence of various strain rates, aspect ratios and boundary conditions. To implement the coupled model into a finite element code, a numerical staggered algorithm is employed. A number of simulations are performed to verify the proposed approach with available experimental and numerical data and to assess its efficiency in solving complex SMA problems.

Keywords

shape memory alloy / thermo-mechanical coupling / superplasticity / shape memory effect

Introduction

Shape memory alloys are categorized among the smart materials, which exhibit a direct interaction between thermal and mechanical fields. Unique behaviors of shape memory alloys (SMAs), including superelasticity and shape memory effect (SME), have made them suitable for various high-tech applications in different areas such as biotechnology and aerospace industry. These behaviors are based on the reversible phase transformation phenomenon between the austenite (A) and martensite (M) phases; associated with the high and low temperature phases, respectively. The phase transformation of SMA is accompanied with heat generation/absorption. As a result, it is trivial that SMAs exhibit a temperature-sensitive behavior; and therefore, the thermo-mechanical coupling effects should be taken into account.

It is accepted that the forward phase transformation (A to M) and the reverse transformation (M to A) in SMAs are exothermic and endothermic processes, respectively [ 1]. In other words, in the case of high loading rates, the material does not have enough time to disperse the heat flow, and its temperature will consequently increase/decrease during the forward/reverse transformation. Therefore, the necessity of considering thermo-mechanical coupling is evident. The influence of high loading rates, particularly strain rates, on thermo-mechanical coupling has been discussed in several studies [ 1, 2]. For instance, Morin and Moumni [ 3] studied the isothermal superelastic behavior of SMAs and the effect of loading rates using the Helmholtz- based ZM constitutive model. However, recent experimental and analytical studies showed that other parameters such as boundary conditions and the size of the SMA device are also important [ 4, 5].

In this study, to capture both the shape memory effect and the superelasticity behaviors and to study the effects of various parameters, the original Lagoudas model [ 6] for thermo-mechanically coupled behavior of SMAs is extended and implemented in a staggered time-stepping algorithm, which uncouples the original problem into two smaller problems. In each iteration, after solving the isothermal mechanical problem, the heat equation is solved to determine the temperature variations due to phase transformation. This process continues until the problem convergences.

This paper consists of three main sections. First, a brief theory of adopted constitutive model is presented. Then, an implementation of the constitutive model in a time-stepping staggered algorithm is presented within the finite element framework. The last part is dedicated to simulation and discussion of several numerical examples to assess the performance of the developed methodology.

Constitutive model

First, a brief discussion of the constitutive model for SMAs is presented. The explicit form of the Gibbs free energy is given by [ 7]:

(1) G(\sigma ,\mathit{T},\xi ,{\epsilon }^t)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\rho }}\sigma :\mathit{S}:\sigma -{\displaystyle \frac{\mathrm{1}}{\rho }}\sigma :[\mathit{a}(\mathit{T}-{\mathit{T}}_{\mathrm{0}})+{\epsilon }^t]+\mathit{c}\left[(\mathit{T}-{\mathit{T}}_{\mathrm{0}})-\mathit{T}\ln \left({\displaystyle \frac{T}{T_{\mathrm{0}}}}\right)\right]-s_{\mathrm{0}}T+u_{\mathrm{0}}+{\displaystyle \frac{\mathrm{1}}{\rho }}f(\xi ),

where the state variables \sigma ,\xi ,{\epsilon }^t,\mathit{T},{\mathit{T}}_{\mathrm{0}} are the Cauchy stress, the martensitic volume fraction, the transformation strain, the temperature and the reference temperature, respectively. The function f(\xi ) is a representative form for the transformation hardening due to phase interaction. The material parameters \rho ,S,\mathit{a},\mathit{c},{\mathit{s}}_{\mathrm{0}} and u_{\mathrm{0}} are the density, the fourth-order effective compliance tensor, the second-order effective thermal expansion tensor, the effective specific heat, the effective specific entropy at the reference state, and the effective specific internal energy at the reference state, respectively. Assuming isotropic and linear elastic behavior for each phase, the rule of mixtures allows for determining these material parameters,

(2) \mathit{S}(\xi )={\mathit{S}}^A+\xi ({\mathit{S}}^M-{\mathit{S}}^A)={\mathit{S}}^A+\xi \mathrm{\Delta }S,

(3) \mathit{a}(\xi )={\mathit{a}}^A+\xi ({\mathit{a}}^M-{\mathit{a}}^A)={\mathit{a}}^A+\xi \mathrm{\Delta }a,

(4) \mathit{c}(\xi )={\mathit{c}}^A+\xi ({\mathit{c}}^M-{\mathit{c}}^A)={\mathit{c}}^A+\xi \mathrm{\Delta }c,

(5) {\mathit{s}}_{\mathrm{0}}(\xi )={\mathit{s}}_{\mathit{A}}^{}+\xi ({\mathit{s}}_{\mathit{M}}^{}-{\mathit{s}}_{\mathit{A}}^{})={\mathit{s}}_{\mathit{A}}^{}+\xi \mathrm{\Delta }{\mathit{s}}_{\mathrm{0}},

(6) {\mathit{u}}_{\mathrm{0}}(\xi )={\mathit{u}}_{\mathit{A}}^{}+\xi ({\mathit{u}}_{\mathit{M}}^{}-{\mathit{u}}_{\mathit{A}}^{})={\mathit{u}}_{\mathit{A}}^{}+\xi \mathrm{\Delta }{\mathit{u}}_{\mathrm{0}},

where superscripts A and M indicate the austenitic and martensitic phases, respectively and Δ represents the difference of a particular parameter between the mentioned phases. The constitutive relation for the total strain can be derived from substituting Eq. (1) into the first and second laws of thermodynamics:

(7) \epsilon =-\rho {\displaystyle \frac{\partial G}{\partial \sigma }}=S:\sigma +a(\mathit{T}-{\mathit{T}}_{\mathrm{0}})+{\epsilon }^t.

Equation (7) is not solely capable of fully describing the material behavior so the necessity of describing a phase transformation rule is evident. The transformation rule is defined as the relation between the evolution of transformation strain tensor and the evolution of martensitic volume fraction:

(8) {\dot{\epsilon }}^t=\mathrm{\Lambda }\dot{\xi },

where \mathrm{\Lambda } is the transformation tensor:

(9)

where the material parameter H is associated with the maximum uniaxial transformation strain, {\sigma }^{\prime } is the deviatoric stress tensor, and {\epsilon }^{tr-R} is the transformation strain at the reversal point. To solve the constitutive relation numerically, an implicit scheme based on the return mapping and the Newton-Raphson method is exploited. For further information, refer to [ 7, 8].

To discuss the discrete form of the heat equation, first the fully coupled form of the heat equation of SMAs is presented [ 9]:

(10) Ta:\dot{\sigma }+\rho c\dot{T}+\left(-\mathrm{\pi }+\mathit{{\rm T}}\mathrm{\Delta }\mathit{\alpha }:\sigma -\rho \mathrm{\Delta }c\ln \left({\displaystyle \frac{T}{T_0}}\right)+\rho \mathrm{\Delta }{\mathit{s}}_{0}T\right)\dot{\xi }=-\nabla .q,

where π is the thermodynamic force

(11) \mathrm{\pi }=\sigma :\mathrm{\Lambda }-\rho {\displaystyle \frac{\partial G}{\partial \xi }}.

The first and third terms of the left-hand side of Eq. (10) express how the temperature varies due to stress variations, and due to a change in martensitic volume fraction, respectively. The second term of the left-hand side of Eq. (10) represents the heat capacity. The first and second terms on the right-hand side of Eq. (10) are related to the heat transfer process by the heat flux q.

The thermal conductivity, the effective specific heat, and the thermal expansion coefficient of the two phases are identical for most SMAs [ 8]. By this assumption and using the Fourier’s law of heat conduction q=-\mathit{k}\nabla T , Eq. (10) yield to:

(12) Ta:\dot{\sigma }+\rho c\dot{T}+\left(-\mathrm{\pi }\rho \mathrm{\Delta }{s}_{0}T\right)\dot{\xi }=\nabla .\left(\mathit{k}\nabla T\right).

Now, the general type of convection boundary condition is imposed:

(13) q.n=-(\mathit{k}\nabla T).\mathit{n}=h(T-T_{\mathrm{e}\mathrm{x}\mathrm{t}}),

where n is the outward normal vector, h is the convection coefficient and T_ext is the external temperature. Temperature distribution at t = 0 is given as the initial condition T= T_ext.

By considering the Dirichlet and Neumann temperature boundary conditions and applying Stokes’ theorem along with the conventional finite element discretization, Eq. (12) yields to

(14) (\mathit{K}+{\mathit{K}}_h+{\mathit{K}}_{\sigma }+{\mathit{K}}_{\xi }).\mathit{T}\mathit{+}\mathit{M}\mathit{.}\dot{\mathit{T}}\mathit{=}\mathit{F},

or in an incremental form:

(15) (\mathit{K}+{\mathit{K}}_h+{\mathit{K}}_{\sigma }^{n+\mathrm{1}}+{\mathit{K}}_{\xi }^{n+\mathrm{1}}).{\mathit{T}}^{n+\mathrm{1}}\mathit{+}\mathit{M}\mathit{.}{\displaystyle \frac{T^{n+\mathrm{1}}-T^n}{\mathrm{\Delta }t}}\mathit{=}{\mathit{F}}^{n+\mathrm{1}},

where the superscript n represents the increment, \mathrm{\Delta }t is the time of each increment, K is the heat conduction matrix, {\mathit{K}}_h is the heat convection matrix, {\mathit{K}}_{\sigma } is the matrix of latent heat due to stress changes, {\mathit{K}}_{\xi } is the matrix of latent heat due to martensite volume fraction changes, M is the heat capacity matrix and F is a load vector:

(16) {(\mathit{K})}_{ij}\mathit{=}{\displaystyle \underset{\mathrm{\Omega }}{\int }k\nabla N_i}\nabla N_j\mathrm{d}\mathrm{\Omega },

(17) {({\mathit{K}}_h)}_{ij}\mathit{=}{\displaystyle \underset{S}{\int }h{N}_i}{N}_j\mathrm{d}S,

(18) {({\mathit{K}}_{\sigma })}_{ij}\mathit{=}{\displaystyle \underset{\mathrm{\Omega }}{\int }a:\stackrel{\cdot }{\sigma }{N}_i}{N}_j\mathrm{d}\mathrm{\Omega },

(19) {({\mathit{K}}_{\xi })}_{ij}\mathit{=}{\displaystyle \underset{\mathrm{\Omega }}{\int }\rho \mathrm{\Delta }{s}_{\mathrm{0}}\dot{\xi }{N}_i}{N}_j\mathrm{d}\mathrm{\Omega },

(20) {(M)}_{ij}\mathit{=}{\displaystyle \underset{\mathrm{\Omega }}{\int }\rho c{N}_i}{N}_j\mathrm{d}\mathrm{\Omega },

(21) {\left(F\right)}_i\mathit{=}{\displaystyle \underset{S}{\int }h{T}_{\mathrm{e}\mathrm{x}\mathrm{t}}{N}_i}\mathrm{d}\mathit{S}{\displaystyle \underset{\mathrm{\Omega }}{\int }\mathrm{\pi }\dot{\xi }{N}_i}\mathrm{d}\mathrm{\Omega },

where N_i and N_j are the finite element shape functions. Ultimately, the temperature at the increment n + 1 is obtained from:

(22) T^{n+\mathrm{1}}={({\mathit{K}}_{im}^{n+\mathrm{1}})}^{-\mathrm{1}}({\mathit{F}}_{im}^{n+\mathrm{1}}),

where {\mathit{K}}_{im}^{n+\mathrm{1}} and {\mathit{F}}_{im}^{n+\mathrm{1}} are:

(23) {\mathit{K}}_{im}^{n+\mathrm{1}}=\left(\mathit{K}+{\mathit{K}}_h+{\mathit{K}}_{\sigma }^{n+\mathrm{1}}+{\mathit{K}}_{\xi }^{n+\mathrm{1}}+{\displaystyle \frac{M}{\mathrm{\Delta }t}}\right),

(24) {\mathit{F}}_{im}^{n+\mathrm{1}}=\left({\mathit{F}}^{n+\mathrm{1}}+{\displaystyle \frac{M}{\mathrm{\Delta }t}}{T}^n\right).

Numerical algorithm

To implement the thermo-mechanical coupled model into a finite element code, a staggered time-step numerical algorithm is employed. According to this algorithm, in each iteration, first the isothermal mechanical problem is solved. To numerically solve the constitutive equation, the return mapping algorithm based on an elastic prediction and phase transformation corrector is employed. In the elastic prediction part, the phase transformation strain is assumed constant, and other variables are calculated. Then, in the phase transformation step, the strain and other variables are corrected. After solving the mechanical part, the thermal effects are calculated and considered for the next mechanical step. This step-by-step algorithm is given bellow. Note that in this algorithm n is the step number and k represents the iteration number in each step.

1) Solve the balance equation and calculate the displacement vector

2) Calculate the strain tensor from the consistency equation

3) Return mapping algorithm

(a) Set initial values

(25) k=\mathrm{0},\hspace{1em}{\xi }_{n+\mathrm{1}}^{(\mathrm{0})}={\xi }_n,\hspace{1em}{\epsilon }_{n+\mathrm{1}}^{tr(\mathrm{0})}={\epsilon }_n^{tr},\hspace{1em}{S}_{n+\mathrm{1}}^{(\mathrm{0})}=S_n,\hspace{1em}{a}_{n+\mathrm{1}}^{(\mathrm{0})}=a_n.

(b) Elastic prediction and checking the phase transformation condition

(26) {\sigma }_{n+\mathrm{1}}^{(k)}=C_{n+\mathrm{1}}^{(k)}:({\epsilon }_{n+\mathrm{1}}-a_{n+\mathrm{1}}^{(k)}({\mathit{T}}_{n+\mathrm{1}}-{\mathit{T}}_{\mathrm{0}})-{\epsilon }_{n+\mathrm{1}}^{tr(k)}),

(27) {\phi }_{n+\mathrm{1}}^{(k)}=\phi ({\sigma }_{n+\mathrm{1}}^{(k)},{\mathit{T}}_{n+\mathrm{1}},{\xi }_{n+\mathrm{1}}^{(k)}).

(c) If \left|{\phi }_{n+\mathrm{1}}^{(k)}\right|\le to{l}_{\mathrm{1}} , return to step 4

(d) Calculate variations of martensitic volume fraction \xi and phase transformation strain {\epsilon }^{tr}

(28) \mathrm{\Delta }{\xi }_{n+\mathrm{1}}^{(k)}={\displaystyle \frac{{\mathit{\Phi }}_{n+\mathrm{1}}^{(k)}}{\pm {\partial }_{\sigma }{\mathit{\Phi }}_{n+\mathrm{1}}^{(k)}:C_{n+\mathrm{1}}^{(k)}:{\partial }_{\sigma }{\mathit{\Phi }}_{n+\mathrm{1}}^{(k)}-{\partial }_{\epsilon }{\mathit{\Phi }}_{n+\mathrm{1}}^{(k)}}},

(29) \mathrm{\Delta }{\epsilon }_{n+\mathrm{1}}^{tr(\mathit{k})}=\mathrm{\Delta }{\xi }_{n+\mathrm{1}}^{(k)}{\mathrm{\Lambda }}_{n+\mathrm{1}}^{(\mathit{k})},

where+ and − are for the forward and reverse transformations, respectively.

(e) Updating state variables; phase transformation strain and martensitic volume fraction

(30) {\xi }_{n+\mathrm{1}}^{(k+\mathrm{1})}={\xi }_{n+\mathrm{1}}^{(k)}+\mathrm{\Delta }{\xi }_{n+\mathrm{1}}^{(k)},

(31) {\epsilon }_{n+1}^{tr\left(\mathit{k}\mathit{+}\right)}={\epsilon }_{n+1}^{tr\left(\mathit{k}\right)}+\mathrm{\Delta }{\epsilon }_{n+1}^{tr\left(\mathit{k}\right)}.

(f) Updating material parameters; elastic and thermal moduli

(32) C(\xi )=C^A+\xi (C^M-C^A)=C^A+\xi \mathrm{\Delta }C,

(33) a(\xi )=a^A+\xi (a^M-a^A)=a^A+\xi \mathrm{\Delta }a,

where C^A and C^M are the stress influence coefficients of austenite and martensite phases, respectively (see Eq. (36)-(39)).

(g) Set k = k + 1 and return to (a).

4) Calculate the temperature T_{n+\mathrm{1}}^{k+\mathrm{1}}

5) Calculate the residual internal force (R)

(34) R=f_{\mathrm{e}\mathrm{x}\mathrm{t}}-f_{\mathrm{int}},

where f_{\mathrm{e}\mathrm{x}\mathrm{t}} is the equivalent applied load and f_{\mathrm{int}} is

(35) f_{\mathrm{int}}={\displaystyle \underset{\mathrm{\Omega }}{\int }{B}^T\sigma \mathrm{d}\mathrm{\Omega }},

where B is the matrix of shape function derivatives.

6) If \left|R\right|\le to{l}_{\mathrm{2}} or \left|T_{n+\mathrm{1}}^{k+\mathrm{1}}-T_{n+\mathrm{1}}^k\right|\le to{l}_{\mathrm{3}} , return to step 1

Numerical simulations

Tensile test

The first simulation is performed to assess the accuracy of the proposed model in three different loading rates in comparison with the loading-unloading in a pseudo-elasticity experimental problem, previously reported by [ 10]. The simulations are carried out on a mesh of 450-quadrilateral elements to model a thin strip subjected to the displacement control condition, as depicted in Fig. 1. For this simulation L, w, and t, are chosen 30, 2.6, and 0.5, respectively. The applied displacement is set to u = 2.4 mm. The specimens are fixed at both ends and are elongated up to the strain of 0.08. The original tests were conducted in the air, where the ambient and initial temperatures were equal to 23°C. Due to the fact that some of the parameters of the model were not reported in [ 10], a model calibration is required. The transformation surfaces can be defined as [ 9]:

(36) {\sigma }_{MS}=C^M(\mathit{T}-{\mathit{T}}_{MS}),

(37) {\sigma }_{Mf}=C^M(\mathit{T}-{\mathit{T}}_{Mf}),

(38) {\sigma }_{AS}=C^A(\mathit{T}-{\mathit{T}}_{AS}),

(39) {\sigma }_{Af}=C^A(\mathit{T}-{\mathit{T}}_{Af}),

where T is the material temperature. {\sigma }_{MS} , {\sigma }_{Mf} , {\sigma }_{AS} and {\sigma }_{Af} are the start and finish stresses of forward and reverse transformations, respectively, which are captured from the stress versus strain diagram provided in [ 10]. Consequently, the unknown material parameters T_{Ms} , T_{Mf} , T_{As} , T_{Af} can be calculated from Eqs. (36) − (39). Table 1 lists the adopted material properties. The average convection coefficient in air is taken from Ref. [ 11] and the conduction coefficient is adopted for a typical Ni Ti [ 8]. The convection coefficient in water is assumed 300 W/(m²K), which is required to study the effects of convection boundary conditions.

The stress- strain curves at different strain rates of 1.1 × 10⁻¹, 1.1 × 10⁻² and 1.1 × 10⁻⁴ are presented in Fig. 2 and are in good agreement with the experimental data. The slight difference at the onset of forward transformation is due to nucleation phenomenon which is a consequence of unstable phase transformation. According to the fact that in practice, SMAs are often trained to shake out this slight instability, particularly for engineering purposes, it is neglected in this study.

Note that there is a high heat exchange between the specimen and the grips which can lead to faster decrease in the temperature of the material during the unloading, causing higher slopes of the reverse phase transformation plateau. Hence, by increasing the strain rate, the numerical and experimental curves tend to deviate from each other during the reverse transformation. However, at slow strain rates (Fig. 2(c)), which is near the adiabatic condition, the specimen would not have enough time to dissipate the latent heat to the environment; consequently, slopes of the reverse transformation plateau in the experiment and numerical simulation become nearly identical.

Convective boundary condition

To study the effect of convective boundary condition, simulations are now performed in different environments; air and water. The results are demonstrated in Fig. 3. Comparing the stress versus strain diagrams of air and water in different strain rates, shows that the water boundary generates much less variations than the air. Due to the higher value of convection coefficient of water, heat exchange with the environment is higher. The increase in the convection coefficient leads to a decrease in the slope of the stress-strain curve. Thereby, a near-isothermal material response of SMAs in highly convective environments such as water is trivially expected (Fig. 3(e)).

Effect of the strain rate

Now all the simulations are performed in air. As shown in Fig. 4, the increase in strain rate will result in higher stress levels during the forward transformation. The reason can be attributed to the fact that during the forward transformation, the temperature of material increases monotonically with the increase of the strain rate, resulting in higher stress levels. However, the increase in strain rate does not necessarily result in lower stress levels during the reverse transformation.

According to Fig. 5, with the increase of the strain rate, the stress level of reverse transformation is decreased. However, further increase in the strain rate results adversely. The reason is that at these strain rates, there is not enough time for the released latent heat to completely transfer to the environment. Hence, the temperature of the material becomes much higher than the ambient temperature at the start of the reverse transformation, which leads to higher stress levels. Consequently, depending on the value of strain rate, the temperature of the material may become higher or lower than or even equal to the ambient temperature after the unloading.

Note that the behavior of SMAs at different strain rates can be explained by a comparison between the dominance of heat exchange or internal heat generation/absorption during the forward/reverse phase change [ 9]. At intermediate strain rates, the amount of latent heat generation/ absorption due to the forward/reverse transformation is comparable to the amount of heat exchange with the environment due to convection. Hence, as it can be seen in Fig. 2, a slightly nonlinear behavior is observed in the transformation region of stress-strain curves. However, the increase/decrease in the temperature will make the heat exchange or heat generation/absorption dominant, and the transformation region of the stress-strain would consequently become nearly linear.

SME versus superelasticity for a tensile SMA

The same specimen with dimensions of 18 mm×10 mm×1 mm is now simulated at two temperatures of-76ºC and-10ºC, and is subjected to an elongation equivalent to 9.5% strain. The material properties are changed to Table 2 [ 8] to allow for better capturing the differences between the shape memory effect (SME) and superelasticity.

According to the results of Fig. 6, the loading paths in superelasticity and SME are respectively shown with A-B-C-D and A-F-H-D, while the unloading paths follow D-E-G-A and D-E-I, respectively. In superelasticity, as expected, the generated strain during the loading is fully recovered upon the unloading. In SME, the unloading process cannot recover the entire strain and a strain equal to 6.7% is remained. Further heating can initiate the reverse transformation and the strain will become fully recovered under zero stress condition (I-A).

The effect of aspect ratio

To study the effect of aspect ratio, three thin strips with different length/width ratios of 5, 15, and 100 are simulated with 248- regular quadrilateral element specimens under the pure tension condition. The material properties are chosen as Table 2. The specimens are elongated up to a given strain of 0.085 at the constant strain rate of 1.1 × 10⁻² under the displacement control condition. Note that the temperature is fixed at the right and left sides of the strips and is equal to 55°C.

The stress- strain responses for different aspect ratios of 5, 15, and 100 are presented in Fig. 7. During the forward transformation, the material response changes from nearly isothermal to nearly adiabatic condition by the change of length/width ratio. To explain this phenomenon, one should note that different heat transfers in various regions can cause a non-uniform distribution of temperature (higher temperature at inner regions), which due to thermo-mechanical coupling results in non-uniform distribution of stress (higher stress at inner regions). However, the lower stress plateau in Fig. 7 indicates that the non-monotonic stress level changes during the reverse transformation do not necessarily lead to lower stress levels.

Bending test

The second simulation is a three point bending test to assess the necessity of considering thermo-mechanical coupling. The geometry, loading, and boundary conditions are adopted from [ 12] and are shown in Fig. 8. The original test was conducted in the air and on a wire with the length and diameter of 20 and 1.49 mm, respectively. The constant thickness of finite element of the present simulation is chosen in such a way that the same moment of inertia is achieved. Material parameters are adopted from [ 8] and are given in Table 3. The plane-stress element is elongated up to 5.2 mm with a constant strain rate of 1.1 × 10⁻¹ under the displacement control condition.

Variations of the force versus displacement of midpoint of the specimen in comparison with numerical and experimental results of reference [ 12] are shown in Fig. 9. Auricchio et al. [ 12] used an uncoupled scheme to solve the 322-element model with the exponential hardening function; while neglecting the existing strong thermo-mechanical coupling may lead to inaccuracy. Here, the simulation is performed on a model of 915 quadrilateral elements and the thermo-mechanical coupling is taken into account. As expected, the results of this study shows better agreement with the experimental data.

To demonstrate the change of phase transformation patterns during the uniform 40-step increasing bending loading, the phase transformation profiles at 8-step intervals are shown in Fig. 10. At first, the material is fully austenite, and the martensitic volume fraction is zero. By increasing the load, the material starts to behave elastically. Further increase in loading initiates the phase transformation and consequently the martensitic volume fraction begins to increase. The temperature boundary condition is assumed fixed at the right and left sides of the structure which causes lower temperature for points closer to these sides compared with the inner regions. Considering this fact and occurrence of stress localization under the point load, it can be concluded that the phase transformation should start in the middle of the beam, right under the point load. Eventually, the value of the martensitic volume fraction reaches 1 and the phase transformation completes.

Strip with hole

The third simulation is performed on a thin strip with a central circular hole. The geometry, loading, and boundary condition are shown in Fig. 11(a). Material parameters are adopted from [ 8] and are given in Table 2. The strip is elongated up to 5% strain under the displacement control condition. Due to the symmetry, only a quarter of specimen is simulated, as depicted in Fig. 11(b).

The change of temperature and stress contours are demonstrated in Figs. 12 and 13, respectively. The contours are at five equal loading steps; in other words, each contour represents a 20% progression in the forward phase transformation. Clearly, the temperature evolution and the stress generation are similar. They initiate from the point A on the edge of the circular hole, where the stress concentration is maximum. By increasing the load, the phase transformation propagates diagonally through the specimen. This phenomenon is accompanied with a rise in temperature in the transformed region, also associated with the stress increase. As expected, the temperature and stress variation are minimum at the point B on the edge of the circular hole due to the minimum stress concentration in that part of the specimen.

Conclusion

In this contribution, a staggered algorithm for solving thermo-mechanical coupling in SMAs is presented. The presented framework is capable of simulating the thermo-mechanical coupling for both superelasticity and SME phenomena. Several numerical simulations are performed to investigate the sensitivity of effective parameters and to examine accuracy and efficiency of the model. It is demonstrated that an isothermal process cannot be guaranteed with the slow strain rate. To analyze the sensitivity of material response, other variables such the convective boundary conditions should accurately be taken into consideration.

References

[1] Shaw J, Kyriakides S. Thermomechanical aspects of NiTi. Journal of the Mechanics and Physics of Solids, 1995, 43(8): 1243–1281

[2] Chang B. Shaw, Iadicola M. Thermodynamics of shape memory alloy wire: Modeling, experiments, and application. Continuum Mechanics and Thermodynamics, 2006, 18(1−2): 83–118

[3] Morin C, Moumni Z. Thermomechanical coupling in shape memory alloys under cyclic loadings: Experimental analysis and constitutive modeling. International Journal of Plasticity, 2011, 27: 1959–1980

[4] Desroches R, McCormick J, Delemont M. Cyclic properties of superelastic shape memory alloy wires and bars. Journal of Structural Engineering, 2004, 130(1): 38–46

[5] Mirzaeifar R, Desroches R, Yavari A. Analysis of the rate-dependent coupled thermo-mechanical response of shape memory alloy bars and wires in tension. Continuum Mechanics and Thermodynamics, 2011, 23: 363–385

[6] Lagoudas D C, Bo Z, Qidwai M A. A unified thermodynamic constitutive model for SMA and finite element analysis of active metal matrix composites. Mechanics of Composite Materials and Structures, 1996, 3: 153–179

[7] Qidwai M A, Lagoudas D C. Numerical implementation of shape memory alloy thermomechanical constitutive model using return mapping algorithm. International Journal for Numerical Methods in Engineering, 2000, 47: 1123–1168

[8] Lagoudas D C. Shape Memory Alloys, Modeling and Engineering Applications. Springer, 2008

[9] He Y J, Sun Q P. On non-monotonic rate dependence of stress hysteresis of superelastic shape memory alloy bars. International Journal of Solids and Structures, 2011, 48: 1688–1695

[10] Zhange X H, Feng P, He Y J, Yu T X, Sun Q P. Experimental study on rate dependence of macroscopic domain and stress hysteresis in NiTi shape memory alloy strips. International Journal of Mechanical Sciences, 2010, 52: 1660–1670

[11] Morin C, Moumni Z, Zaki W. A constitutive model for shape memory alloys accounting for thermomechanical coupling. International Journal of Plasticity, 2011, 27: 748–767

[12] Auricchio F, Taylor R L, Lubliner J. Shape-memory alloys: macromodelling and numerical simulations of the superelastic behavior. Computer Methods in Applied Mechanics and Engineering, 1997, 146(3−4): 281–312

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg