A three-dimensional constitutive model for the martensitic transformation in polycrystalline shape memory alloys under large deformation

TL;DR

This paper introduces a 3D constitutive model for polycrystalline shape memory alloys that accurately captures large deformations and rotations using logarithmic strain, improving upon existing models by reducing artificial stress errors.

Contribution

The paper develops a novel thermodynamics-based model utilizing logarithmic strain and rate, with detailed implementation for simulating large deformation behaviors in SMAs.

Findings

Model accurately predicts cyclic loading behavior of SMA components.

Implementation in finite element software shows reduced artificial stress errors.

Numerical experiments validate the model's effectiveness for large deformation analysis.

Abstract

This work presents a three-dimensional constitutive model for the martensitic transformation in polycrystalline Shape Memory Alloys (SMAs) under large deformation. By utilizing the logarithmic strain and rate, the model is able to account for large strains and rotations that SMA-based components may undertake, but also resolves the artificial stress errors caused by the non-integrable objective rates that are widely used in current commercial finite element software. The model is developed through classical thermodynamic laws combined with the standard Coleman-Noll procedure. The scalar martensitic volume fraction and the second-order transformation strain tensor are chosen as the internal state variables to capture the material response exhibited by polycrystalline SMAs. A detailed implementation procedure of the proposed model is described through a user-defined material subroutine. Numerical experiments considering SMA components including a bar, a beam, a torque tube and a solid flexible structure under stress/thermally-induced phase transformations are investigated via the proposed model, and the results under cyclic loading are compared against the predictions provided by the Abaqus nonlinear solver. The development framework of the proposed model and its implementation procedure can be extended to incorporate other nonlinear phenomena exhibited by SMAs, such as transformation-induced plasticity, viscoplasticity, and damage under large deformation.

Tables

Table 1. Table 2: A set of representative material parameters used for the parametric numerical study [ 1 , 29 ] .

Type	Parameter	Value	Parameter	Value
	$E_A$	60 [GPa]	$C_A$	8 [MPa/K]
	$E_M$	40 [GPa]	$C_M$	6 [MPa/K]
Key material parameters	${\nu }_A={\nu }_M$	0.3	$M_s$	333 [K]
12	${\alpha }_A={\alpha }_M$	1.0$\times$10-5 [K-1]	$M_f$	220 [K]
	$H^{\text{max}}$	3%, 5%, 8%	$A_s$	274 [K]
	$k_t$	0.02	$A_f$	370 [K]
Smooth hardening parameters	$n_1$	0.5	$n_3$	0.5
4	$n_2$	0.5	$n_4$	0.5

Table 2. Table 3: Elastic modulus and transformation strain calibrated based on engineering and true scale stress-strain curves.

	Engineering scale	Case 1	Case 2	Case 3	True scale	Case 1	Case 2	Case 3
	$E_A$ [GPa]	60.00	60.00	60.00	$E_A$ [GPa]	60.00	60.00	60.00
	$E_M$[GPa]	35.48	34.06	32.05	$E_M$[GPa]	40.00	40.00	40.00
	$H_{\text{eng}}^{max}$	$3.04\%$	$5.13\%$	$8.33\%$	$H^{max}$	$3.00\%$	$5.00\%$	$8.00\%$

Table 3. Table 4: Calibrated values of material parameters for equiatomic NiTi [ 29 ] .

Type	Parameter	Value	Parameter	Value
	$E_A$	60 [GPa]	$C_A$	7.8 [MPa/K]
	$E_M$	60 [GPa]	$C_M$	7.3 [MPa/K]
Key material parameters	${\nu }_A={\nu }_M$	0.3	$M_s$	333 [K]
12	${\alpha }_A={\alpha }_M$	1.0$\times$10-5 [K-1]	$M_f$	220 [K]
	$H^{\text{max}}$	4.7%	$A_s$	274 [K]
	$k_t$	0.021	$A_f$	370 [K]
Smooth hardening parameters	$n_1$	0.5	$n_3$	0.5
4	$n_2$	0.5	$n_4$	0.5

Table 4. Table 5: Calibrated values of material parameters for NiTi (50.8 at.% Ni) [ 29 ] .

Type	Parameter	Value	Parameter444The values of transformation temperatures $(M_s,M_f,A_s,A_f)$ are referenced from [81] in order to realize the self-expanding process within human body environment, the rest of values of the material parameters are taken from [29] as they are not provided from [81].	Value
	$E_A$	32.5 [GPa]	$C_A$	3.5 [MPa/K]
	$E_M$	23.0 [GPa]	$C_M$	3.5 [MPa/K]
Key material parameters	${\nu }_A={\nu }_M$	0.3	$M_s$	264 [K]
12	${\alpha }_A={\alpha }_M$	2.2$\times$10-5 [K-1]	$M_f$	160 [K]
	$H^{\text{max}}$	3.3%	$A_s$	217 [K]
	$k_t$	N/A	$A_f$	290 [K]
Smooth hardening parameters	$n_1$	0.17	$n_3$	0.25
4	$n_2$	0.27	$n_4$	0.35

Table 5. Table 6: Stress residuals obtained by hypoelastic equation ( 62 ) with different objective rates for the elastic square

	Normalized stress	Logarithmic rate	Jaumman rate	Green-Naghdi rate
	${\tau }_{11}/E$	-1.71e-5	0.0164	2.17 e-6
	${\tau }_{12}/E$	2.17e-6	0.0411	-0.023
	${\tau }_{22}/E$	1.71e-5	-0.0164	1.76e-5

Table 6. Table 7: Calibrated material parameters of a NiTi SMA used for the model validation.

Type	Parameter	Value	Parameter	Value
	$E_A$	41 [GPa]	$C_A$	5.5 [MPa/K]
	$E_M$	22 [GPa]	$C_M$	5.5 [MPa/K]
Key material parameters	${\nu }_A={\nu }_M$	0.33	$M_s$	237 [K]
12	${\alpha }_A={\alpha }_M$	1.0$\times$10-5 [K-1]	$M_f$	217.5 [K]
	$H^{\text{max}}$	3.35%	$A_s$	254 [K]
	$k_t$	N/A	$A_f$	282 [K]
Smooth hardening parameters	$n_1$	0.15	$n_3$	0.25
4	$n_2$	0.17	$n_4$	0.15

Equations

$$\mathbf{F}(\mathbf{x},t)=\frac{\partial \mathbf{x}}{\partial \mathbf{X}}$$

$$\mathbf{v}=\frac{d\mathbf{x}}{dt}=\dot{\mathbf{x}}$$

$$\mathbf{L}=\frac{\partial \mathbf{v}}{\partial \mathbf{x}}=\dot{\mathbf{F}}{\mathbf{F}}^{-1}$$

$$\mathbf{F}=\mathbf{RU}=\mathbf{VR}$$

$$\mathbf{C}={\mathbf{F}}^{\mathbf{T}}\mathbf{F}={\mathbf{U}}^2$$

$$\mathbf{B}=\mathbf{F}{\mathbf{F}}^{\mathbf{T}}={\mathbf{V}}^2$$

$$\mathbf{H}=\frac{1}{2}\ln (\mathbf{C})=\ln (\mathbf{U})$$

$$\mathbf{h}=\frac{1}{2}\ln (\mathbf{B})=\ln (\mathbf{V})$$

$$\mathbf{L}=\mathbf{D}+\mathbf{W},\mathbf{D}=\frac{1}{2}(\mathbf{L}+{\mathbf{L}}^{\mathbf{T}}),\text{vspace}5pt\mathbf{W}=\frac{1}{2}(\mathbf{L}-{\mathbf{L}}^{\mathbf{T}})$$

$${\mathring{\mathbf{h}}}^{log}=\dot{\mathbf{h}}+\mathbf{h}{\mathbf{\Omega }}^{log}-{\mathbf{\Omega }}^{log}\mathbf{h}=\mathbf{D}$$

$$\mathbf{\Omega}^{log}=\mathbf{W}+\sum_{i\neq j}^{n}\big{(}\frac{1+(\lambda_{i}/\lambda_{j})}{1-(\lambda_{i}/\lambda_{j})}+\frac{2}{\ln(\lambda_{i}/\lambda_{j})}\big{)}\mathbf{b}_{i}\mathbf{D}\mathbf{b}_{j}$$

$${\mathbf{\Omega }}^{log}={\dot{\mathbf{R}}}^{log}({\mathbf{R}}^{log}{)}^T$$

$$\mathbf{h}=\int_{\text{corot.}}\mathbf{D}~{}\text{d}t=(\mathbf{R}^{log})^{T}\bigg{(}\int_{0}^{t}\mathbf{R}^{log}\mathbf{D}^{e}(\mathbf{R}^{log})^{T}\text{d}t^{\prime}\bigg{)}\mathbf{R}^{log}$$

$$\mathbf{D}={\mathbf{D}}^e+{\mathbf{D}}^{tr}$$

$${\mathring{\mathbf{h}}}^{e\_log}={\mathbf{D}}^e;{\mathring{\mathbf{h}}}^{tr\_log}={\mathbf{D}}^{tr}$$

$${\mathring{\mathbf{h}}}^{log}={\mathring{\mathbf{h}}}^{e\_log}+{\mathring{\mathbf{h}}}^{tr\_log}$$

$${\mathbf{h}}^e$$

$${\mathbf{h}}^{tr}$$

$$\mathbf{h}={\mathbf{h}}^e+{\mathbf{h}}^{tr}$$

$$G(\text{uptau},\mathbf{h},T,\mathbf{{\rm Y}})=u-\frac{1}{{\rho }_0}\text{uptau}:\mathbf{h}-sT$$

$$\mathcal{D}=\text{uptau}:\mathbf{D}-{\rho }_0(\dot{u}-T\dot{s})⩾0$$

$${\mathring{G}}^{log}=\dot{u}-\frac{1}{{\rho }_0}{\mathring{\text{uptau}}}^{log}:\mathbf{h}-\frac{1}{{\rho }_0}\text{uptau}:{\mathring{\mathbf{h}}}^{log}-s\dot{T}-\dot{s}T$$

$$\dot{u}-\dot{s}T=\dot{G}+\frac{1}{{\rho }_0}\mathring{\text{uptau}}:\mathbf{h}+\frac{1}{{\rho }_0}\text{uptau}:\mathring{\mathbf{h}}+s\dot{T}$$

$$\mathcal{D}=-{\rho }_0\dot{G}-{\rho }_{0}s\dot{T}-\mathring{\text{uptau}}:\mathbf{h}⩾0$$

$$\mathring{G}=\frac{\partial G}{\partial \text{uptau}}:\mathring{\text{uptau}}+\frac{\partial G}{\partial T}\dot{T}+\frac{\partial G}{\partial \mathbf{{\rm Y}}}:\mathring{\mathbf{{\rm Y}}}$$

$$\mathcal{D}=-({\rho }_0\frac{\partial G}{\partial \text{uptau}}+\mathbf{h}):\mathring{\text{uptau}}-({\rho }_0\frac{\partial G}{\partial T}+s):\dot{T}-{\rho }_0\frac{\partial G}{\partial \mathbf{{\rm Y}}}:\mathring{\mathbf{{\rm Y}}}⩾0$$

$$\mathbf{h}=-{\rho }_0\frac{\partial G}{\partial \text{uptau}}$$

$$s=-{\rho }_0\frac{\partial G}{\partial T}$$

$$-{\rho }_0\frac{\partial G}{\partial \mathbf{{\rm Y}}}:\mathring{\mathbf{{\rm Y}}}⩾0$$

$$\displaystyle G=-\dfrac{1}{2\rho_{0}}\bm{\uptau}:\mathcal{S}\bm{\uptau}-\dfrac{1}{\rho_{0}}\bm{\uptau}:[~{}\bm{\alpha}(T-T_{0})+\mathbf{h}^{tr}]+c\Big{[}(T-T_{0})-T\ln(\dfrac{T}{T_{0}})\Big{]}-s_{0}(T-T_{0})+u_{0}+\dfrac{1}{\rho_{0}}f(\xi)$$

$$\mathcal{S}(\xi )={\mathcal{S}}^A+\xi ({\mathcal{S}}^M-{\mathcal{S}}^A)={\mathcal{S}}^A+\xi \Delta \mathcal{S}$$

Full text

1 INTRODUCTION

SMAs belong to a specialized subgroup of multifunctional materials known as active materials, and are capable of recovering their pre-defined geometry when subjected to a thermal stimulus above certain temperatures. This unique property of SMAs is achieved through a solid-to-solid state diffusionless phase transformation between the high-symmetry, high-temperature austenitic phase and the low-symmetry, low-temperature martensitic phase [1]. Since the discovery of shape memory effect, SMAs have been extensively investigated as sensors and actuators towards building smart systems integrated with adaptive and morphing features [2, 3]. Recently, aerospace researchers have considered to use the SMA-based actuators to reconfigure the shape of a supersonic aircraft to meet the noise and efficiency requirements in response to the real-time changing ambient environment, which has the great potential to realize a commercially viable overland civil supersonic flight in the near future [4, 5].

A substantial number of constitutive theories for SMAs have been proposed so far with the majority of them based on the infinitesimal strain theory for small deformation analysis. Thorough reviews can be found from [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. In general, constitutive models for SMAs can be approximately categorized into three different types: phase-field theory based models, crystal-plasticity theory based models, and classical $J_{2}$-flow theory based models. The phase-field models, in which order parameters are utilized to differentiate austenitic and martensitic phases, can track microstructure evolution, such as phases front movement, during the phase transformation process [11, 12, 20, 21, 22, 23, 24]. Therefore, phase-field theory based models are well suited to investigate the dynamic nucleation and phase morphology growth for martensitic phase transformation in SMAs. The enormous computational time, however, needed to solve the phase kinetic partial differential equations hinders its popularity to the extent that macroscopic structural response is concerned. In regard to crystal-plasticity theory based models [25, 26, 27], with the consideration of the effect of material microstructure (e.g., crystal orientation, texture, etc.), these models are able to capture the anisotropy in material response (e.g., tension-compression asymmetry) exhibited by textured polycrystalline SMAs. Similar to phase-field methods, the complex implementation procedure of these models to incorporate the microstructure information makes them computational costly. Following the legacy of $J_{2}$-flow theory, phenomenological SMA models have attracted attention in the engineering community where repetitive designs and optimization procedures on SMA components are needed to find target shapes. By introducing a set of internal state variables (such as volume fraction and transformation strain tensor), $J_{2}$ theory-based SMA models are able to simulate the macroscopic response of an SMA component in an efficient way. The simplicity of this model type and its well-established implementation procedure have allowed it to be widely used among real engineering applications [28, 29, 30, 31, 32, 33, 14, 34].

Constitutive models based on the infinitesimal strain theory are able to predict SMA response accurately under small deformation situations. However, it has been reported that SMAs can reversibly deform to a relatively large strain regime up to 8% [35, 36]. Also for specific boundary value problem such as fracture in SMAs, the strain levels close to the crack tip are well beyond 10% [37, 38] within the finite strain regime. In addition to such relatively large strain, SMA-based actuators (e.g., spring, beam, torque tube) may also undergo large rotation during their deployment. For example, an SMA beam component is used as a bending actuator in order to realize a morphing engine shape, where the SMA experienced large bending induced rotation [3]. Another example is the SMA tube component utilized as a torsional actuator to repeatedly rotate a deployable and retractable solar panel [39], where the SMA tube is subjected to cyclic large rotation. Combining the aforementioned two facts that SMAs may undergo large strain and rotation, it is necessary to develop a constitutive model based on a finite deformation framework to provide an accurate prediction for the response of SMAs.

Two kinematic assumptions are often employed in the finite deformation theory, i.e., the multiplicative decomposition of the deformation gradient and the additive decomposition of the strain rate tensor. In the multiplicative approach, the deformation gradient is usually decomposed into an elastic part multiplied with an inelastic part. Finite strain SMA models based on the multiplicative decomposition can be obtained from literature [40, 34, 41, 42, 43, 44, 45], among which some advanced capabilities are considered. For examples, Wang and coworkers [43] presented a finite strain SMA model with the fully thermomechanically coupled feature, the consideration of coexistence of different martensitic variants, and accounting for temperature effect on the hysteresis size. In the work of Stupkiewicz and Petryk [44], they proposed a finite strain SMA model to capture the tension-compression asymmetry phenomenon. Damanpack and coworkers [45] also developed an SMA model that considers anisotropic behaviors and reorientation in SMAs at finite deformation. However, it is known that finite strain model based on additive decomposition significantly reduces the complexities of model structure compared to multiplicative models, which in return facilitates the computational efficiency of the finite strain model as a 3-D design tool. Therefore, they are widely used in current available finite element softwares (e.g., Abaqus, ANSYS). However, to satisfy the principle of objectivity, additive models are required to use an objective rate in its rate form hypoelastic constitutive relation. A number of objective rates (such as Zaremba-Jaumann-Noll rate, Green-Naghdi-Dienes rate, and Truesdell rate) have been proposed to meet this goal. However, those objective rates are not essentially “objective” because of their failure to integrate the rate form hypoelastic relation to yield a free energy based hyperelastic stress-strain relation [46]. As a result, spurious phenomena (e.g., shear stress oscillation, artificial stress residuals, etc.) are often observed in the predicted response even for simple elastic materials through these objective rates. More details regarding this issue are addressed in A.

It was not until the logarithmic rate was proposed in the literature [47, 48, 49, 50, 51, 52, 53] that the previously mentioned self-inconsistency issue related to non-integrable objective rates was resolved. As was proved in the work [47], the logarithmic rate of Eulerian logarithmic strain $\mathbf{h}$ is exactly identical to the strain rate tensor $\mathbf{D}$, and the logarithmic strain is the only one among all other strain measures enjoying this important property. Therefore, the finite strain models using logarithmic strain and rate are not only able to capture large strain and large rotation but also are capable of resolving the aforementioned spurious phenomena. This new development in finite deformation theory not only provides solution to classical finite elastoplastic problems for conventional metallic materials [54, 55], but also sheds light on the finite strain model development for active materials such as SMAs. Few SMA models using additive approach can be found from [56, 57, 58, 59, 60, 61, 62], but some of the very important SMA phase transformation characteristics have not been addressed among them, such as the smooth transition during the phase transformation, the stress dependent transformation strain to account for the coexistence of oriented/self-accommodated martensitic variants, and a stress dependent critical driving force to consider the effect of applied stress levels on the size of hysteresis loop. To this end, this work presents a finite strain SMA model formulation based on the additive decomposition using the logarithmic strain and rate. As a continuous development from the infinitesimal SMA model [29], the proposed model has a simple model structure and considers three very important characteristics for SMA response as its infinitesimal counterpart does. These developments combined result in an improved computational efficiency and robustness for the proposed finite strain model to predict the SMA response at large deformation, without introducing additional intermediate state variables, such as Mandel stress, that are utilized in the multiplicative models. It is noted that the primary focus of this work is mainly on the formulation of a finite strain SMA model, rather than the development of a constitutive model that can capture the full complexities of the SMA thermomechanical deformation. Thus, tension-compression, latent-heat effects, reorientation between orientated and self-accommodated martensitic variants, cyclic evolution features (transformation-induced plasticity, two-way shape memory effect at stress-free conditions) are not included here for simplicity. Moreover, this work carefully examines the artificial stress errors caused by using other non-integrable objective rates in current commercial finite element packages. The capability of the proposed model to eliminate such stress errors shows significant importance for the analysis of SMA-based actuators, e.g., SMA beam and SMA torque tube subjected to cyclic large deformation.

This paper is organized as follows. Section 2 presents the kinematic preliminaries. Section 3 concentrates on the model development based on the logarithmic strain and logarithmic rate. The derivation of the consistent tangent stiffness matrix and the consistent thermal matrix are also provided. In section 4, the detailed implementation procedure for the proposed model is described by using a user-defined material subroutine (UMAT) through the finite element software Abaqus. Numerical examples are studied to demonstrate the capability of the proposed model in Section 5. Conclusions are presented in Section 6. A detailed calibration procedure for the material parameters used in this model is also provided in the B.

2 Preliminaries

2.1 Kinematics

Let material point $\mathcal{P}$ from body $\mathcal{B}$ be defined by a position vector $\mathbf{X}$ in the reference (undeformed) configuration at time $t_{0}$, and let vector $\mathbf{x}$ represent the position vector of that material point in the current (deformed) configuration at time $t$. Therefore ,the deformation process of point $\mathcal{P}$ between the reference configuration and the current configuration can be defined through the well-known deformation gradient tensor $\mathbf{F}(\mathbf{x},t)$:

[TABLE]

and the velocity field $\mathbf{v}$ can be defined as,

[TABLE]

based on the velocity field $\mathbf{v}$, the velocity gradient $\mathbf{L}$ can be derived as,

[TABLE]

the following polar decomposition equation for deformation gradient $\mathbf{F}$ is well known,

[TABLE]

where $\mathbf{R}$ is the rotation tensor, $\mathbf{U}$ and $\mathbf{V}$ are the right (or Lagrangian) and the left (or Eulerian) stretch tensors, respectively, by which the right Cauchy-Green tensor $\mathbf{C}$ and the left Cauchy-Green tensor $\mathbf{B}$ can be obtained, as follows,

[TABLE]

[TABLE]

where $\mathbf{I}$ is the second order identity tensor. The logarithmic strain (also called Hencky or true strain) of Lagrangian type $\mathbf{H}$ and Eulerian type $\mathbf{h}$ can thus be defined as,

[TABLE]

[TABLE]

It is also well known that the velocity gradient $\mathbf{L}$ can be additively decomposed into a symmetric part, the strain rate tensor $\mathbf{D}$, and an anti-symmetric part, the spin tensor $\mathbf{W}$,

[TABLE]

2.2 Logarithmic strain, logarithmic rate and logarithmic spin

As was mentioned in section 1, two widely accepted kinematic assumptions, i.e., the multiplicative decomposition of deformation gradient $\mathbf{F}$ and the additive decomposition of the strain rate tensor $\mathbf{D}$, are usually considered in finite deformation theory. The multiplicative models use a hyperelastic constitutive relation while a rate form hypoelastic constitutive equation is usually adopted for additive models. The rate form hypoelastic constitutive theory using objective rates has been criticized for its non-integrability because it can not well define an essential elastic material behavior [63], this includes many well known objective rates such as Zaremba-Jaumann rate, Green-Naghdi rate, Truesdell rate, etc.[46].

The aforementioned problems about objective rates were solved in the work by Xiao et al.[47, 48, 46], Bruhns et al.[49, 50, 51] and Meyers et al.[52, 53], where they proved that the logarithmic rate of the Eulerian logarithmic strain $\mathbf{h}$ is identical with the strain rate tensor $\mathbf{D}$, by which a hypoelastic model can be exactly integrated to a hyperelastic finite strain model [47]. This unique relationship between logarithmic strain $\mathbf{h}$ and the strain rate tensor $\mathbf{D}$ is expressed as follows,

[TABLE]

where $\mathbf{\Omega}^{log}$ is called the logarithmic spin introduced by [47] defined as,

[TABLE]

in which $\lambda_{i,j}(i,j=1,2,3)$ are the eigenvalues of left Cauchy-Green tensor $\mathbf{B}$ and $\mathbf{b}_{i},\mathbf{b}_{j}$ are the corresponding subordinate eigenprojections. As along as the logarithmic spin tensor $\mathbf{\Omega}^{log}$ is defined, the second order rotation tensor $\mathbf{R}^{log}$, associated with $\mathbf{\Omega}^{log}$, can be determined through the following differential equation (12). In general cases, the initial condition is assumed as $\mathbf{R}^{log}|_{t=0}=\mathbf{I}$.

[TABLE]

follow the corotational integration definition from [64], and assume the initial conditions $\mathbf{h}|_{t=0}=\mathbf{0}$, equation (10) yields the total logarithmic strain $\mathbf{h}$ after the logarithmic corotational integration,

[TABLE]

2.3 Additive decomposition of logarithmic strain

The kinematic assumption starts with the additive decomposition of the strain rate tensor $\mathbf{D}$ into an elastic part $\mathbf{D}^{e}$ plus a transformation part $\mathbf{D}^{tr}$,

[TABLE]

The elastic strain rate part $\mathbf{D}^{e}$ and the transformation strain rate part $\mathbf{D}^{tr}$ in equation (14) can be rewritten as $\mathring{\mathbf{h}}^{e\_log}$ and $\mathring{\mathbf{h}}^{tr\_log}$ by virtue of the relation in equation (10) respectively,

[TABLE]

By combining equations (10), (14) and (15), the following equation can be obtained,

[TABLE]

Similar to the results obtained from equation (13), equation (16) can yield the following relation after applying the logarithmic corotational integration,

[TABLE]

Based on the additive decomposition on the strain rate tensor, combing equations (14), (16) and (17), the following additive decomposition on the total logarithmic strain $\mathbf{h}$ can be achieved, i.e., the total logarithmic strain $\mathbf{h}$ can be additively split into an elastic strain like part $\mathbf{h}^{e}$ plus a transformation strain like part $\mathbf{h}^{tr}$.

[TABLE]

3 Model Formulation

3.1 Thermodynamic framework

The Gibbs free energy potential $G$ is defined to be a continuous function dependent on Kirchhoff stress tensor $\bm{\uptau}$ 111The relationship between Kirchhoff stress $\bm{\uptau}$ and Cauchy stress $\bm{\sigma}$ is $\bm{\uptau}=J\bm{\sigma}$, where $J$ is the determinant of the deformation gradient $\mathbf{F}$, i.e., $J=\text{det}|\mathbf{F}|$. Assuming phase transformation to be volume preserving, $J$ is approximately equivalent to 1, so $\bm{\uptau}\approx\bm{\sigma}$. Kirchhoff stress $\bm{\uptau}$ and Eulerian logarithmic strain $\mathbf{h}$, called an energetic conjugate pair [65], are usually paired up in the formation of free energy potentials., Eulerian logarithmic strain $\mathbf{h}$, temperature $T$ and a set of internal state variables $\Upsilon$.

[TABLE]

where $\rho_{0}$ is the density in the reference configuration, $s$ and $u$ are the specific entropy and internal energy respectively. From the $2^{nd}$ law of thermodynamics, the dissipation energy $\mathcal{D}$ can be written in the form of Clausius-Duhem inequality,

[TABLE]

The logarithmic rate of the Gibbs free energy is taken in equation (19). Note that a scalar subjected to an objective rate equals to its conventional time rate, the following equation is derived. An circle hat denotes the logarithmic rate in the following text for brevity.

[TABLE]

the following equation is obtained after rearrangement of equation (21),

[TABLE]

Substitute equation (22) into Clausius-Duhem inequality (20), the dissipation energy is rearranged as the following,

[TABLE]

Recall that the Gibbs free energy $G(\bm{\uptau},T,\bm{\Upsilon})$ is a continuous function, chain rule differentiation can be applied on the Gibbs free energy with respect to its independent state variables (i.e., Kirchhoff stress $\bm{\uptau}$, temperature $T$ and internal state variables $\bm{\Upsilon}$), which gives,

[TABLE]

Substitute equation (24) into equation (23), the following expression for the dissipation energy $\mathcal{D}$ is acquired,

[TABLE]

Following the standard Coleman-Noll procedure, all admissible values for $\mathring{\bm{\uptau}}$, $\dot{T}$ and $\mathring{\bm{\Upsilon}}$ have to comply with the dissipation inequality (25) regardless of thermodynamic paths, thereby the constitutive relationships between stress and strain, entropy and temperature can be inferred as,

[TABLE]

[TABLE]

Substitute equations (26) and (27) into equation (25), the following reduced form of the dissipation inequality is acquired,

[TABLE]

3.2 Constitutive modeling for SMAs

3.2.1 Thermodynamic potential

The formulation of the proposed model is based on the thermodynamic framework presented in section 3 and the early SMA model developed by Lagoudas and coworkers [7, 29] for small deformation analysis. The model is able to predict the pseudoelastic (isothermal) and actuation (isobaric) response under large deformation including both large strain and large rotation. A quadratic Gibbs free energy potential $G$ is introduced in equation (29), in which Kirchhoff stress tensor $\bm{\uptau}$ and temperature $T$ are chosen as the independent state variables. The martensitic volume fraction $\xi$ and the second order transformation strain tensor $\mathbf{h}^{tr}$ are chosen as internal state variables $\mathbf{\Upsilon}=\{\xi,\mathbf{h}^{tr}\}$ to capture the material response exhibited by polycrystalline SMAs. The Gibbs free energy potential $G$ is employed as follows,

[TABLE]

where $\mathcal{S}$ is the effective compliance tensor calculated by equation (30), $\mathcal{S}^{A}$ is the compliance tensor for austenitic phase while $\mathcal{S}^{M}$ is for martensitic phase, and $\Delta\mathcal{S}$ is the phase difference for the compliance tensor. The effective stiffness tensor $\mathcal{C}$ can be gained by taking the inverse of the effective compliance tensor, i.e., $\mathcal{C}=\mathcal{S}^{-1}$. $\bm{\alpha}$ is the second order thermoelastic expansion tensor, $c$ is the effective specific heat, $s_{0}$ and $u_{0}$ are the effective specific entropy and effective specific internal energy at the reference state. All the aforementioned effective variables are determined from equation (31) to (34). $T$ represents the temperature at current state while $T_{0}$ is the temperature at reference state.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A smooth hardening function $f(\xi)$ is proposed in equation (35) to account for the hardening effects in polycrystalline SMAs, such as the plastic strain accumulation after the training procedure, imperfections located at the grain boundary, and nano-precipitates hardening effects, etc.[1], where three additional intermediate material parameters $a_{1},a_{2},a_{3}$ and four curve fitting parameters $n_{1},n_{2},n_{3},n_{4}$ are introduced to better treat the smooth transition behaviors at the initiation and completion of phase transformation.

[TABLE]

Following the standard Coleman-Noll procedure described in section 3, the explicit form for constitutive relation (26) between stress and strain is derived as follows. Note that the nonlinearity in this constitutive relation is implied by the transformation strain $\mathbf{h}^{tr}$.

[TABLE]

the explicit form for constitutive relation (27) between entropy $s$ and temperature $T$ can also be derived as,

[TABLE]

the reduced form of the dissipation inequality (28) can be rewritten in terms of the chosen internal state variables $\bm{\Upsilon}=\{\xi,\mathbf{h}^{tr}\}$ as,

[TABLE]

3.2.2 Evolution equation of internal state variables

The evolution equation for the internal state variables $\bm{\Upsilon}=\{\xi,\mathbf{h}^{tr}\}$ is presented here. It is proposed that the logarithmic rate of the transformation strain $\mathbf{h}^{tr}$ is proportional to the rate change of the martensitic volume fraction $\xi$. This proportional evolution rule is adopted by following the principle of maximum dissipation such that among all the admissible thermodynamic paths, the one dissipating the most energy is chosen during the SMAs phase transformation process [66]. The idea of maximum dissipation for inelastic materials is not new, it was also widely employed for plastic deformed materials to derive the associated flow rule [67]. It is worth pointing out that the rate applied on the transformation strain is the logarithmic rate rather than the conventional time rate. The explicit evolution rule is as follows,

[TABLE]

where $\bm{\Lambda}^{\textit{fwd}}$ is called the forward transformation direction tensor and $\bm{\Lambda}^{\textit{rev}}$ is called the reverse transformation direction tensor. They are defined as,

[TABLE]

in which, $\bm{\uptau}^{{}^{\prime}}$ is the deviatoric part of Kirchhoff stress tensor calculated by $\bm{\uptau}^{{}^{\prime}}=\bm{\uptau}-{\small\frac{1}{3}}\textrm{tr}(\bm{\uptau})~{}\mathbf{1}$, $\mathbf{1}$ is the second order identity tensor. The effective Mises equivalent stress is given by $\bar{\bm{\uptau}}={\small\sqrt{{{\small\frac{3}{2}}\bm{\uptau}}^{{}^{\prime}}:\bm{\uptau}^{{}^{\prime}}}}$. $\mathbf{h}^{\textit{tr-r}}$ and ${\xi}^{r}$ represent the transformation strain value and martensitic volume fraction at the reverse transformation starting point. It is common among available SMA models that the magnitude of the inelastic recoverable transformation strain is the same for full transformation under any applied stress levels. This is true when the stress levels is high enough to generate maximum oriented martensitic variants. However, if the applied stress level is not sufficiently high, self-accommodated martensitic variants will be generated. This renders the value of transformation strain less than it is in the high stress level case (i.e., the stress dependency of the magnitude of the transformation strain). Therefore, an exponential function $H^{\textit{cur}}$ dependent on current stress levels is introduced to calculate the current transformation strain as shown in equation (41), where $H^{\textit{max}}$ is the maximum (or saturated) transformation strain and $\textit{k}_{t}$ is a curve fitting material parameter.

[TABLE]

3.2.3 Transformation function

The objective in this part is to define a proper transformation criterion to determine the occurrence of the phase transformation. Recall the reduced form for dissipation energy is given by inequility (38) and the relation between $\mathbf{h}^{tr}$ and $\xi$ is defined through evolution equation (39). Substituting evolution equation (39) into reduced form dissipation inequality (38), the following equation is obtained,

[TABLE]

the above equation implies that all the dissipation energy is directly a result of the change in the martensitic volume fraction. Based upon this, a scalar variable $\pi$, called the thermodynamic driving force conjugated to the martensitic volume fraction $\xi$, can thus be defined. Substitution of Gibbs free energy potential $G$ in equation (29) into equation (42) yields the explicit expression for $\pi$ as follows,

[TABLE]

where $\Delta\mathcal{S},\Delta\mathbf{\alpha},\Delta c,\Delta s_{0}$, and $\Delta u_{0}$ are the phase differences on compliance tensor, thermal expansion tensor, specific heat, reference entropy and reference internal energy, respectively. It can be observed that the thermodynamic driving force $\pi$ is a function of Kirchhoff stress $\bm{\uptau}$, temperature $T$ and martenstic volume fraction $\xi$. This indicates that the phase transformation process can be activated by two independent sources, namely either the stress or temperature, which correlates quite well with the experimentally observed stress-induced and thermally-induced phase transformations in SMAs. To proceed to the goal of defining a transformation criteria, it is assumed that whenever the thermodynamic driving force $\pi$ reaches a critical value $Y$ ($-Y$), the forward (reverse) phase transformation takes place. Therefore a transformation function $\Phi$ can be defined as the transformation criteria to determine the transformation occurrence as follows,

[TABLE]

In the infinitesimal strain theory based SMA model [29], a reference critical value $Y_{0}$ and an additional parameter $D$ were introduced into $Y$, through which the thermodynamical critical value $Y$ becomes a function dependent on applied stress levels, see in equation (45). Such treatment let the model consider the effect of applied stress levels on the size of hysteresis loop. This capability is provided through capturing the different slopes $C_{A},C_{M}$ in the effective stress-temperature phase diagram. The explicit derivation is provided from equation (70) to equation (74) at the model calibration part in B.

[TABLE]

As a consequence of the application of the principle of maximum dissipation [66], the so-called Kuhn-Tucker constraints are placed on the proposed model, which are stated as follows for the forward and reverse cases respectively,

[TABLE]

3.3 Consistent tangent stiffness and thermal matrix

In this section, a detailed derivation of the consistent tangent stiffness matrix and the thermal matrix is provided to complete the proposed model. For most typical displacement-based finite element softwares, such as Abaqus, the consistent tangent matrices are often required to be provided in the UMAT so that the finite element solver can achieve a fast convergence for the global equilibrium equations. Normally, consistent tangent matrices can be expressed in the rate form shown in equation (47), where $\mathcal{L}$ is called the consistent tangent stiffness matrix and $\Theta$ is the consistent thermal matrix.

[TABLE]

applying the logarithmic rate on constitutive equation (36) yields,

[TABLE]

taking chain rule differentiation on the transformation function equation (44) gives,

[TABLE]

substituting equation (48) back into equation (49) to eliminate $\mathring{\bm{\uptau}}$ and solving it for $\dot{\xi}$, the following expression for $\dot{\xi}$ can be obtained,

[TABLE]

substituting equation(50) back into the rate form constitutive equation(48) to eliminate $\dot{\xi}$, considering the phase difference of the thermal expansion coefficients can be ignored for martensite and austenite phase, the final explicit expression corresponding to equation (47) can be obtained as follows,

[TABLE]

in which the consistent tangent stiffness matrix ${\mathcal{L}}$ is,

[TABLE]

and the consistent thermal matrix $\Theta$ is,

[TABLE]

In order to fully determine the explicit values for $\mathcal{L}$ and ${\Theta}$ during the implementation section for the proposed model, the explicit expressions of the following terms $\partial_{\bm{\uptau}}\Phi$, $\partial_{\xi}\Phi$, $\partial_{T}\Phi$ used in above equations are derived in D.

4 Numerical implementation

This section focuses on the implementation of the proposed model within finite element (FE) solvers to solve boundary value problems (BVPs). The implementation flowchart is shown in figure 1. While typically stress and strain information are provided from FE solver, the initial input information used in this model are only the temperatures $T_{n},~{}\Delta T_{n}$ and deformation gradients at current step $\mathbf{F}_{n}$ and next step $\mathbf{F}_{n+1}$. The reason for using only these information is that other tensorial variables have been rotated by the finite element (FE) solver before they are used as inputs, in which the rotation tensor is calculated based on the other non-integrable objective rates. This consequently leads to the artificial stress errors described in section 1. During the implementation for the proposed model, a pre-calculation and a rotation procedure are employed before calling the main UMAT subroutine. In the pre-calculation procedure, the logarithmic strain at current step $\mathbf{h}_{n}$ and next step $\mathbf{h}_{n+1}$ are calculated based on $\mathbf{F}_{n}$ and $\mathbf{F}_{n+1}$. The incremental rotation tensor $\Delta\mathbf{R}_{n}^{log}$ based on the logarithmic rate can be calculated by using the exponential map scheme [63, 68, 54]. In the rotation procedure, tensorial variables including $\mathbf{h}_{n}$, $\mathbf{h}^{tr}_{n}$ and $\mathbf{\Lambda}_{n}$ are rotated from step $n$ configuration to the new configuration at step $n+1$ by using the obtained $\Delta\mathbf{R}_{n}^{log}$, thus, the so-called principle of objectivity is preserved.

The rest of implementation procedure consists of two steps, the first step is called the thermoelastic predictor and the second step is called the transformation corrector. During the initialization of thermoelastic step, the total strain $\mathbf{h}_{n+1}$ and current temperature $T_{n+1}=T_{n}+\Delta T_{n}$ are provided. The initial internal state variables $\mathbf{\Upsilon}_{n+1}^{(0)}$ are assumed to be the same as $\mathbf{\Upsilon}_{n}$ for the initial Kuhn-Tucker consistency checking, i.e., $\Phi^{(0)}_{n+1}\leqslant 0$. If the initial value of $\Phi^{(0)}_{n+1}$ satisfies the consistency checking, the n+1 step is detected as a thermoelastic response and the UMAT returns to global FE solver for next increment. In case of consistency condition violated, the second step called the transformation corrector is activated to find the updated internal state variables $\mathbf{\Upsilon}_{n+1}^{(k)}$ in order to regain the Kuhn-Tucker consistency. A detailed summary for the implementation procedure is listed in the table 1.

4.1 Thermoelastic Prediction

Take the $(n+1)^{th}$ step as an example to go through the thermoelastic prediction process. The total strain tensor $\mathbf{h}_{n+1}$ and the temperature $T_{n+1}$ are provided from Pre-Calculation procedure, and the initial internal state variables $\mathbf{\Upsilon}_{n+1}^{(0)}$ are assumed the same as $\mathbf{\Upsilon}_{n}$,

[TABLE]

Based on equation (54), the initial guess for Kirchhoff stress $\bm{\uptau}^{(0)}_{n+1}$ can be calculated through the constitutive equation (55). Here the integer in the upper parenthesis represents that how many iterations have been done during the transformation correction procedure, and integer zero means that this step is just an initial guess in the thermoelastic procedure. The initial calculation for stress $\bm{\uptau}^{(0)}_{n+1}$ can be obtained,

[TABLE]

After the calculation of $\bm{\uptau}^{(0)}_{n+1}$, the value of transformation function $\Phi^{(0)}_{n+1}$ can be evaluated based on equations (43) and (44) for the initial Kuhn-Tucker consistency checking,

[TABLE]

If the calculated value of transformation function $\Phi^{(0)}_{n+1}$ remains under the transformation surface (i.e., $\Phi^{(0)}_{n+1}\leqslant\text{'tol'}$, ’tol’ is usually set to be $10^{-6}$), step n+1 is detected as a thermoelastic response. Therefore the values of current state variables $\bm{\uptau}^{(0)}_{n+1}$ and $\mathbf{\Upsilon}^{(0)}_{n+1}$ are accepted as correct and the UMAT proceeds to the global FE solver for the next increment. In case the transformation surface is violated (i.e. $\Phi^{(0)}_{n+1}\geqslant\text{tol}$), the transformation corrector step is activated to find the updated state variables until the consistency equation (46) is preserved.

4.2 Transformation Correction

This part addresses the iterative procedures required for the transformation corrector to restore the Kuhn-Tucker consistency. In general, the transformation corrector is nothing but a set of Newton-Raphson iterations on equations (57) and (58) to find the updated internal state variables. Take the $k^{th}$ local iteration for example, the corrector is activated to find a set of $\mathbf{\Upsilon}^{(k)}_{n+1}$ which makes the residual terms $\mathbf{R}_{n+1}^{tr(k)}$ in equation (57) and transformation function $\Phi^{(k)}_{n+1}$ in equation (58) less than ’tol’.

[TABLE]

[TABLE]

This objective is equivalent to the following convergence conditions,

[TABLE]

Use the standard Newton-Raphson procedure 222The explicit expression for the Jacobian matrix during this Newton-Raphson iteration in equation (60) is quite complicated. The symbolic calculation tool in MATLAB is used here to find the Jacobian matrix, and the authors suggest interested readers to utilize this method to perform the tedious calculation. to solve equations (57) and (58),

[TABLE]

The following results on internal state variables at $(k+1)^{th}$ iteration can be obtained,

[TABLE]

Once the converged values of $\{\mathbf{h}^{tr(k+1)}_{n+1},\xi^{(k+1)}_{n+1}\}$ are found, the current transformation corrector step is labeled as finished and the UMAT proceeds to the next increment. Otherwise the Newton-Raphson procedure exits at this step after certain number of iterations and the current finite element increment step stops.

[TABLE]

5 Numerical Results

In this section, the proposed model is used to predict the stress/thermally-induced phase transformations in SMAs subjected to general three-dimensional thermo-mechanical loading. Several numerical examples are presented here to test the capabilities of this model to account for large strains and rotations, and also to resolve the artificial stress errors issue. First, a parametric analysis on a uniaxial SMA bar is studied to show that the proposed model is able to consider the geometry nonlinearity induced by large strains. Second, two BVPs, i.e., an SMA beam and an SMA torque tube subjected to stress-induced phase transformations, are tested as large rotation cases. To show the model is able to resolve the artificial stress errors issue, the cyclic response of the beam and the torque tube are obtained via the proposed model, and the results are compared against the predictions obtained by the Abaqus nonlinear solver333As the nonlinear solver is activated for implicit analysis (i.e., select NLGEOM on), Abaqus automatically use the logarithmic strain as its strain measure, and the Jaumman rate is the utilized objective rate to account for the large rotation[69].. Next, an isobaric BVP of an SMA torque tube subjected to varying thermal loading is investigated to predict the thermally-induced phase transformation. In the end, to show the model is able to capture the non-proportional local stress and strain evolutions, a 3-D solid flexible structure undergoing a self-expanding process is studied. The proposed model is anticipated to be further validated against experimental data of NiTi and NiTiHf SMAs under uniaxial and other non-uniform loading conditions. The ultimate objective is to validate the capability of the proposed model to predict the response of SMA-based actuators, such as SMA beams and torque tubes, which are intended to be integrated with the future supersonic transport aircrafts to realize the morphing capabilities to reduce the sonic boom noise.

5.1 SMA bar under isothermal loading

To test the capability of the proposed model to account for the effects of large strain, an SMA prismatic bar is studied under uniaxial isothermal loading condition. A parametric study is performed with the maximum transformation strain $H^{max}=3\%,5\%,8\%$ to represent three different loading cases. A group of representative material parameters (two material parameter groups combined) used in this example are listed in table 2 referenced from [1, 29]. The SMA prismatic bar has a length $L=100$ (mm) and an square cross section with an edge length $a=10$ (mm). It is subjected to a proportional force loading up to 120 (kN) then unloading to 0 (kN), the temperature is kept constant at $380$K throughout the process. Generally, the load-displacement curves provided from such uniaxial experiments are interpreted into the engineering scale stress-strain curves to facilitate the model calibration. However, when the materials experience a strain that is no longer considered small, the geometry nonlinearity due to such strain has to be taken into consideration. To demonstrate that the proposed model accounts for this, the calibrated values of elastic modulus ($E_{A},E_{M}$) and maximum transformation strain ($H^{max}$) based on the true stress-strain curve are compared against the values from its infinitesimal counterpart. Three sets of load-displacement curves are generated shown in figure 2(a). They are interpreted into stress-strain curves in two scales, i.e, the true stress (Cauchy stress) versus the true strain (logarithmic strain) curve and the engineering stress (nominal stress) versus engineering strain (infinitesimal strain) curve. By using the calibration procedure described in B, the calibrated values of $E_{A},E_{M}$ and $H^{max}$ summarized in table 3.

Table 3 shows that the values of $E_{A}$ are identical in both the two scales. However, the values of $E_{M}$ in engineering scale change from 35.48 GPa to 32.05 GPa, which indicates a material softening. Actually, such material softning is not real. Instead, it is the effect of disregarding the geometric nonlinearity induced by large strain as described previously. In this case, the geometric nonlinearity means that the bar needs to contract its cross section to compensate for its elongation to preserve the volume conservation. Disregarding the change of cross section results in an unreal decreasing on the values of $E_{M}$. By doing the calibration based on the true stress-strain curve instead of the engineering one, the proposed model is able to exclude the geometry nonlinearity induced by large strain, so that the calibrated values of $E_{M}$ remain the same in the three loading cases from true scale. Besides, The values of $H^{max}$ are also worth to be noted. Although $H^{max}$ shows different values in the two scales, a relationship exists between the true scale $H^{max}$ and the engineering scale $H^{max}_{\text{eng}}$, i.e., $H^{max}=\ln(1+H^{max}_{\text{eng}})$. Based on the results from this parametric study, it is shown that the infinitesimal strain assumption may no longer be considered as an accurate approximation when the strain regime is beyond 3%. In order to account for the effects caused by large strain, a finite strain model to consider the geometry nonlinearity is required even in a uniaxial case.

5.2 SMA beam under isothermal loading

The second BVP considered here is an SMA beam subjected to isothermal loading shown in figure 3. The SMA beam component has been investigated as bending actuators in [3] to realize a morphing variable-geometry chevron in order to change the outer engine shell shape to achieve specific aerodynamic characteristics. While only one loading cycle was considered in the previous study, this example examines the cyclic material and structural response. The studied beam has the same geometry as the SMA bar in section 5.1. Refer to figure 3, the beam is simply supported with one node being fixed to suppress the rigid body motion, and the upper face is subjected to a traction that ramps up to 24 (MPa) then decreases to 0 (MPa). Temperature is kept constant at $380$K throughout the whole numerical experiment. Material parameters used in this simulation are summarized in table 4. The cyclic material and structural response are obtained by the proposed model for a material point p (in figure 3) located at the middle bottom position, and are compared against the results obtained from the Abaqus nonlinear solver.

Figure 4 shows the obtained cyclic longitudinal stress-strain curve for material point p under tension. As shown in figure 4(a), the proposed model provided a stable material response, while the Abaqus nonlinear solver predicted a shifting, instead of stable, response shown in figure 4(b). The observation from these results indicate that the spurious material response is obtained due to the usage of non-integrable objective rates in Abaqus as discussed in the introduction. Although the initial several loading cycles are almost the same in the results provided by Abaqus nonlinear solver, the accumulation of artificially introduced stress errors, around -2 MPa for each cycle, gradually drifts the material response left downwards throughout the 100 loading cycles. In total, -200 MPa stress residuals together with -0.6% remnant strains are observed at the end. Such stress errors consist of almost 18% of the maximum stress levels experienced by material point p. As a comparison, figure 5 shows 100 stress-strain curves for another material point subjected to compression at the middle of beam upper surface. The result shows an opposite shifting trend in contrast to the results of point p. Again, a stable compressive stress-strain curve are predicted by the proposed model while the Abaqus nonlinear solver predicts a shifting one. In addition, figure 6 shows the obtained cyclic load-displacement curves for point p. It can be seen that the proposed model predicted a stable structural response while the Abaqus nonlinear solver predicted a shifting structural response. Based on these results, it is demonstrated that the Abaqus nonlinear solver can no longer produce reliable results for the SMA beam subjected to 100 bending cycles. Therefore, the proposed model with the capability to eliminate the stress errors is required for the SMA beam subjected to cyclic loading.

5.3 SMA tube under isothermal loading

In this subsection, the BVP of a three-dimensional SMA torque tube under torsion loading is studied. Refer to figure 7(a), the tube has an inner radius $r=3.0$(mm) and thickness $t/r=0.1$. In order to reduce the computational cost, a representative tube segment $L/r=2/3$ is analyzed here. Boundary conditions are depicted in figure 7(b), the tube left face is fixed and the right face is subjected to a torsion loading. The torque proportionally increases to 25 (N$\cdot$m) then unloads to 0 (N$\cdot$m), the temperature is kept constant at $380$ K. The torque tube undergoes a fully forward phase transformation from austenitic phase to martensitic phase followed by a reverse phase transformation from martensitic phase to austenitic phase. The material parameters used in this simulation are from table 4.

The cyclic shear stress-strain response of a material point from the tube outer surface at the right end is predicted by the proposed model and Abaqus nonlinear solver. As shown in figure 9(b), similar to the results observed in the SMA beam case, a shifting response is predicted by the Abaqus nonlinear solver due to the accumulation of shear stress errors. In contrast, a stable response is predicted by the proposed model shown in figure 9(a). More specifically, figure 8 shows the magnitude of stress residual accumulated after one loading cycle. The value of Mises stress residual predicted by the proposed model is almost zero compared to a value around 4 MPa predicted by the the Abaqus nonlinear solver. As a result of the accumulation of such shear stress errors, the shear stress levels required to start the forward phase transformation spuriously decreases in the case of Abaqus nonlinear solver shown in figure 9(b). Besides, the maximum shear stress levels at the end of forward transformation increases, and the shape of hysteresis loop also changes. The cyclic structural response of the torque tube is also provided in figure 10 by plotting the applied torque versus the twist angle $\theta_{z}$. It can be seen that a stable structural response is predicted by the proposed model shown in figure 10(a) compared to a shifting structural response predicted by the Abaqus nonlinear solver shown in figure 10(b). From the observation on these results, it is seen that the Abaqus nonlinear solver is not able to predict reliable results for the SMA torque tube subjected to 100 shearing cycles any more. Thus, the proposed model that can resolve the shear stress errors is required for the SMA torque tube subjected to cyclic torsion loading.

5.4 SMA tube under isobaric loading

In order to test the capability of the proposed model to predict the thermally-induced phase transformation in SMAs, a three-dimensional cylindrical SMA tube is studied under cyclic isobaric torsional loading, i.e., subjected to a constant torsion load with temperature variation cycles. The SMA torque tubes has been investigated as rotational actuators to realize a morphing wing during the plane take-off and cruise regime [70, 71, 72, 73]. The design and optimization of such SMA-based morphing structure requires a thorough understanding on the response of SMA torque tubes subjected to cyclic isobaric loading. To that end, the SMA tube component is analyzed under cyclic isobaric loading conditions. The model has the same geometry and material information as the tube simulation in section 5.3. The loading condition is as follows, a 3 N$\cdot$m torque load is applied to the tube right end and the temperature varies between $250$ K and $390$ K for 100 cycles. Cyclic shear strain-temperature and $\bm{\theta}_{z}$-temperature curves are obtained via the proposed model and the Abaqus nonlinear solver.

As it is shown in figure 11(a), a stable cyclic shear strain-temperature response for the tube is predicted by the proposed model. In contrast, a shifting cyclic strain-temperature response is predicted by using the Abaqus nonlinear solver shown in figure 11(b). More specifically, it can be seen that the isobaric response drifts downwards at $T=240$ K and is lifted up at $T=390$ K from cycle to cycle due to the stress error accumulation. Similar to the observation on the strain-temperature response, the cyclic $\bm{\theta}_{z}$-temperature response is stable in the case of proposed model, and it is a shifting response predicted by the Abaqus nonlinear solver. The comparison on $\bm{\theta}_{z}$-temperature response is plotted in figure 12. Based upon the analysis of SMA tube subjected to thermal loading cycles, it is shown that the accumulated stress errors from Abaqus nonlinear solver result in an shifting cyclic isobaric response, and such artificial stress errors can be eliminated by using the proposed model.

5.5 3-D analysis of a flexible SMA structure

In this section, a 3-D solid SMA flexible structure is studied to demonstrate that the proposed model can be used as a 3-D structural design tool. This SMA flexible structure can be used as a stent to provide a less invasive method for administering support to diseased arteries, veins or other vessels in the human body. They are crimped into a smaller shape outside the body then inserted into the diseased artery. After being delivered into the desired position, the flexible structure expands automatically by using the phase transformation of SMAs [74]. The loading path of this self-expanding structure can be described with three steps, i.e., 1) The flexure is firstly crimped outside the body by external constraints and attached to a constraint container device called catheter or cannula. 2) The flexure is inserted into the body while the temperature increases from the room temperature to the body temperature. 3) The structure recovers its original shape when the constraint is removed [75, 76, 77]. Although there are analysis for similar type of flexible structure by other researchers [78, 79], among which only a small unit cell of the structure is analyzed. Here a full scale 3-D SMA flexible structure designed for repairing aortic dissection [74, 80] is studied to provide the global response of the structure during the self-expanding process. The structure is 32 (mm) long, 25.4 (mm) in outer diameter with 4 struts and 0.5 (mm) thick in radial direction. Due to the curvature of the structure strut, the SMA flexure experiences large rotation and stress concentrations around the hinge part, which in return results in a complex local non-proportional stress and strain evolutions at the hinge. The loading path of this analysis is indicated by the red curve in figure.13(b). The material parameters used in this simulation are from table 5.

The expanded and crimped shapes of the SMA flexure during the analysis are shown in figure 13(a). As illustrated in figure 13(c), stress concentration due to the strut curved part is observed. The von Mises contour indicates that a local non-proportional stress field is evolved at the hinge location during the crimping process. As shown in figure 13(d), while the straight strut part is still in austenitic phase, the stress-induced martensitic phase transformation is activated by the stress concentration at the hinge location subjected to bending. The martensitic volume fraction contour for the flexure during the crimping process is shown in figure 13(d). The global structural response of the SMA flexure for a material point P (see figure 13(c)) is provided in terms of a 3-D stress-temperature-displacement curve in figure.13(b), in which the red curve indicates the actual response while the blue curve is the projection of blue curve on the stress-radial reduction ratio plane. This well captured non-proportional stress evolution and martensitic phase transformation in this example demonstrates that the proposed model can be used as an efficient tool for the 3-D analysis and design of complex SMA-based structures.

6 Concluding remarks

Based on the SMA model proposed by Lagoudas and coworkers for small deformation analysis, a three-dimensional constitutive model for martensitic transformation in polycrystalline SMAs accounting for large deformation has been proposed in this work. Three important characteristics in SMA reponse are considered, i.e., the smooth transition during the phase transformation, the stress dependent transformation strain to account for the coexistence of oriented/self-accommodated martensitic variants, and a stress dependent critical driving force to consider the effect of applied stress levels on the size of hysteresis loop. The proposed model is formulated based on the finite deformation framework that utilizes logarithmic strain and rate such that it not only accounts for the large strains and rotations that an SMA component may undertake, but also resolves the artificial stress errors that are caused by using other non-integrable objective rates. The proposed model is able to predict the stress-induced and thermally induced phase transformations in SMAs under general three-dimensional thermomechanical loading. In particular, it was shown in the example of an SMA bar that the proposed model accounts for the geometry nonlinearity induced by large strains, so that it corrects the spurious material softening in the results from its infinitesimal counterpart. In the numerical examples of an SMA beam and an SMA torque tube, it was demonstrated that the proposed model captures the large rotations that SMA-based components may undertake. By comparing the predicted cyclic response to the results obtained through the Abaqus nonlinear solver, the proposed model demonstrated that it effectively resolves the artificial stress errors. In the end, a 3-D solid flexible structure experiencing local, non-proportional stress and strain evolution was analyzed by the proposed model, which shows the proposed model can be used as an efficient tool for the 3-D analysis and design of complex SMA-based structures. The detailed formulation of the proposed model and its implementation procedures make it readily used by other researchers. The model can be further extended to incorporate additional nonlinear phenomena exhibited by SMAs, such as transformation-induced plasticity, viscoplasticity, and damage evolution.

Acknowledgments

The authors would like to acknowledge the financial support provided by the Qatar National Research Fund (QNRF) under grant number: NPRP 7-032-2-016, and the NASA University Leadership Initiative (ULI) project under the grant number: NNX17AJ96A. The conclusions in this work are solely made by the authors and do not necessarily represent the perspectives of QNRF and NASA.

Data availability

The user-defined material subroutine implemented in commercial finite element software Abaqus for the proposed finite strain SMA model is available to download from the GitHub: https://github.com/Aero-tomato/SMA-UMAT.

Appendix A Artificial stress residuals due to the non-integrable objective rates

In this appendix section, BVP is investigated to study the artificial stress residuals caused by using other non-integrable objective rates. As it is discussed in the introduction, the rate form hypoelastic constitutive theory has been criticized for its inconsistent choices on objective rates [63], this includes the well-known objective rates (e.g., Zaremba-Jaumann rate, Green-Naghdi rate, Truesdell rate, etc). In other words, the rate form hypoelastic constitutive equation fails to be integrated to deliver an algebraic hyperelastic constitutive equation via the so-called objective rates, because of which spurious phenomena (e.g., shear stress oscillation, artificial stress residuals, etc.) are often observed in a simple elastic deformation [46]. It was not until recently such self-inconsistent issues about hypoelastic constitutive models has been resolved by the logarithmic rate proposed by [47, 48, 46, 49, 50, 51, 52, 53]. As their work proved that the logarithmic rate of the logarithmic strain $\mathbf{h}$ of its Eulerian type is equivalent to the strain rate $\mathbf{D}$, by which a grade-zero hypoelastic model can be exactly integrated into an finite strain elastic model based on logarithmic strain.

A.1 The cyclic response of an elastic square

Refer to figure 14 for the BVP schematics, a two-dimensional elastic square with length $H$ is under a closed path cyclic loading. The upper line of the square is subjected to a displacement control circular deformation, the deformation over geometry ratio is as $r/H=0.2$ to induce large deformation strain. The stress components are examined by the hypoelastic equation (62) for 10 loading cycles. In equation (62), $\mathcal{\bm{C}}$ is the stiffness tensor and an circle over $\bm{\uptau}$ means the different objective rates adopted.

[TABLE]

Based on the results from [47], hypoelastic equation (62) can be self-consistently integrated to a hyperelastic constitutive equation (63) based on logarithmic strain through the logarithmic corotational integration [64].

[TABLE]

Kirchhoff stress components are obtained by equation (62) with three different objective rates, i.e. Jaumman rate, Green-Naghdi rate and Logarithmic rate. The predicted results are presented in figure 15, figure 16 and figure 17. The stress results are normalized by the material Young’s modulus $E$.

The stress residuals are examined at the end of the loading cycle and summarized in table 6. First, the stress components in all the three cases showed periodic oscillation. Since the material is confined to behave elastically, the deformation should be indissipative to anticipate that all the stress components should return to zero value in the end. However, the predicted stress components in the cases of Jaumman rate and Green-Naghdi rate showed artificial stress residuals are introduced. Refer to figure 15 for the case of Jaumman rate, stress residuals $\uptau_{11}$ is 0.0164, $\uptau_{22}$ is -0.0164 and $\uptau_{12}$ is 0.0411 after the 10 loading cycles. In the case of Green-Naghdi rate, although there are inconsiderable stress residuals for $\uptau_{11}$ and $\uptau_{22}$ components, the shear residuals $\uptau_{12}$ is -0.023. In contrast, all the stress residuals are almost negligible in the case of Logarithmic rate, which demonstrates that the hypoelastic constitutive equation utilizing logarithmic rate can be self-consistently integrated to deliver a hyperelastic equation based on the logarithmic strain. Interested readers are encouraged to further read [52, 82].

Appendix B Calibration of the material parameters

In this section, the material parameters utilized in the proposed model are identified from a set of one-dimensional experimental data. Note that the strain measure used here should be in the true (or logarithmic) scale rather than the engineering (or infinitesimal) scale. Material parameters used in the proposed model can be categorized into three groups, i.e. the key material parameters, smooth hardening parameters and intermediate parameters. First, the material constants such as elastic modulus $E_{A},E_{M}$, Poisson’s ratios $\nu_{A}$ and $\nu_{B}$, the thermal expansion tensors $\bm{\alpha}_{A}$ and $\bm{\alpha}_{M}$, stress influenced coefficients $C_{A}$ and $C_{M}$ from the phase diagram (or called clausius clapeyron coefficient), critical phase transformation temperatures $A_{s},A_{f},M_{s},M_{f}$ at stress free state are determined. Secondly, the hardening parameters describing the smooth transition feature are discussed. Finally, the intermediate parameters are derived based on the aforementioned two parameter groups. All the material parameters used in this model are summarized in table 2.

Because the data is provided in one-dimensional case, all tensorial variables of the proposed model have to be reduced into 1-D scalar value. For example, the stress tensor is reduced as $\bm{\uptau}\rightarrow\uptau_{11}=\uptau$; logarithmic strain tensor is reduced as $\mathbf{h}\rightarrow h_{11}=h$, etc. Constitutive equation (36) can be rewritten as one dimensional form as follows,

[TABLE]

where the effective elastic modulus $E$ is calculated by using the rule of mixture as follows,

[TABLE]

The evolution equation (39) is also reduced in one-dimensional form as,

[TABLE]

the thermodynamic driving force $\pi$ in one-dimensional case can thus be obtained,

[TABLE]

the transformation function in one-dimensional form can be calculated based on equation (67). Considering the phase difference for the thermal expansion $\Delta\alpha$ and specific heat $\Delta c$ are small enough to be ignored, the following transformation functions for the forward case and the reverse case can be obtained respectively,

[TABLE]

[TABLE]

As described in the first paragraph of this section, there are three sets of material parameters that need to be identified. First, let’s consider the material constants $(E_{A},E_{M},\nu_{A},\nu_{M},\alpha_{A},\alpha_{M})$. Elastic modulus $E_{A},E_{M}$ can be determined through a pseudoelastic stress and strain curve by calculating the slopes at martensitic phase and austenite phase. Poisson’s ratio is attained using a widely accepted value of $\nu_{A}=\nu_{M}=0.33$ found in [29]. The thermal expansion coefficient are usually considered as $\alpha_{A}=\alpha_{M}$, which can be calibrated through an isobaric actuation experiment. The maximum transformation strain $H^{max}$ can be determined from the pseudoelastic experimental and the value of parameter $k_{t}$ are chosen to best fit the $H^{cur}$ curve. The stress influence coefficients and the critical phase transformation temperatures $(C_{A},C_{M},M_{s},M_{f},A_{s},A_{f})$ can be calibrated through the phase diagram. Second, the material parameters related to the smooth hardening features are discussed. The Coefficients $n_{1},n_{2},n_{3},n_{4}$ without specific physical meanings are determined to best match the smoothness in corners of material response. Lastly, there are seven intermediate material parameters $(\rho_{0}\Delta s_{0},\rho_{0}\Delta u_{0},a_{1},a_{2},a_{3},Y_{0}~{}\text{and}~{}D)$ that need to be calculated to complete the model. Determination of such intermediate parameters requires a set of seven algebraic equations. The needed four equations come from transformation constraints as the Kuhn-Tucker condition Equation.46 (i.e. ${\Phi}_{rev}(\uptau,T,\xi)=0$). The fifth equation comes from the continuity of Gibbs free energy at the end of the forward transformation $(\xi=1)$. The needed five algebraic equations are summarized as follows,

1. (i)

Start of the forward transformation at zero stress $(\uptau=0;T=M_{s};\xi=0)$.

${\Phi}_{\textit{fwd}}(0,M_{s},0)=\rho_{0}\Delta s_{0}M_{s}-\rho_{0}\Delta u_{0}-a_{3}-Y_{0}=0$ 2. (ii)

Finish of the forward transformation at zero stress $(\uptau=0;T=M_{f};\xi=1)$.

${\Phi}_{\textit{fwd}}(0,M_{f},1)=\rho_{0}\Delta s_{0}M_{f}-\rho_{0}\Delta u_{0}-a_{1}-a_{3}-Y_{0}=0$ 3. (iii)

Start of reverse transformation under zero stress $(\uptau=0;T=A_{s};\xi=1)$.

${\Phi}_{\textit{rev}}(0,A_{s},1)=-\rho_{0}\Delta s_{0}M_{s}+\rho_{0}\Delta u_{0}+a_{2}-a_{3}-Y_{0}=0$ 4. (iv)

Finish of the reverse transformation under zero stress $(\uptau=0;T=A_{f};\xi=0)$.

${\Phi}_{\textit{rev}}(0,A_{f},0)=-\rho_{0}\Delta s_{0}A_{f}+\rho_{0}\Delta u_{0}-a_{3}-Y_{0}=0$ 5. (v)

The gibbs free energy continuity at the end of forward transformation $(\xi=1)$.

$f(\xi=1)|_{\dot{\xi}\geq 0}$ = $f(\xi=1)|_{\dot{\xi}\leq 0}$

The above five algebraic equations yield the following expression for the five out of seven intermediate model parameters,

[TABLE]

Another two equations are derived from the Kuhn-Tucker condition in order to complete the calculation. For a one-dimensional uniaxial experiment, the Kuhn-Tucker condition (46) requires equation (71) to hold true at any specific stress level $\uptau^{*}$,

[TABLE]

Evaluate $\text{d}\Phi$ at the start point of the forward phase transformation (i.e. $\xi=0$), and at the finish point of the forward phase transformation (i.e. $\xi=1$), the incremental part of martensitic volume fraction should be zero (i.e. $\text{d}\xi=0$) in both of the aforementioned cases. Therefore, the relationships between the stress temperature coefficients $C_{M}$, $C_{A}$ and the stress temperature slopes $\frac{\text{d}\uptau}{\text{d}T}$ can be obtained.

For the forward transformation case, $\dot{\xi}>0$,

[TABLE]

For the reverse transformation case, $\dot{\xi}<0$

[TABLE]

Using the equations (72) and (73), the rest two intermediate material parameters $\rho_{0}\Delta s_{0}$ and $D$ can thus be expressed as follows,

[TABLE]

[TABLE]

Appendix C Model validation

For the validation of proposed model against experimental results, a uniaxial pseudoelastic tensile test is performed on a NiTi SMA [1]. In the experiment, a NiTi SMA strip is loaded at a constant temperature of $320$ K, which is larger than the austenistic finish temperature $A_{f}$. The strip is subjected to a traction up to $600$ MPa then unloaded to [math] MPa. The material parameters listed in table 7 are calibrated based on the experimental result. Figure 18 shows the comparison between the experimental results and the simulations from the proposed model, it clearly demonstrates that the proposed model predicts the stress-strain response of the NiTi strip quite well, including the phase transformation starting and finishing points, and the size of the hysteresis loop. As the experimentally tested specimen is an untrained SMA, there is a amount of transformation-induced plastic strain remained at the end of the loading. This phenomenon can be captured by extending the current model with the consideration of the transformation-induced plasticity.

Appendix D Supplementary calculation for consistent tangent stiffness and thermal matrix

The consistent tangent stiffness and thermal matrix are derived in section 3.3. In order to determine the explicit values for $\mathcal{L}$ and ${\Theta}$ during the implementation of the proposed model, the explicit expressions of the following terms $\partial_{\bm{\uptau}}\Phi$, $\partial_{\xi}\Phi$, $\partial_{T}\Phi$ used in equation (51) are needed. First, the partial derivative of transformation function $\Phi$ with respect to stress $\bm{\uptau}$ can be obtained through differentiating equation (44) by $\bm{\uptau}$. Utilize the expression for $\pi$ in equation (43), it obtains,

[TABLE]

where the partial derivative of the thermodynamic driving force $\pi$ with respect to stress $\bm{\uptau}$ is,

[TABLE]

and the partial derivative of critical driving force value $Y$ with respect to stress $\bm{\uptau}$ is,

[TABLE]

based on the expression for the transformation direction tensor in equation (40), the partial derivative of $\bm{\Lambda}$ with respect to stress $\bm{\uptau}$ are provided for the forward and reverse transformation cases as follows,

[TABLE]

where the partial derivative of the term \Big{(}\dfrac{\bm{\uptau}^{{}^{\prime}}}{\bar{\bm{\uptau}}}\Big{)} with respect to stress $\bm{\uptau}$ is provided in the following equation, in which $\mathbf{I}$ is the forth order identity tensor and $\mathbf{1}$ is the second order identity tensor. It can be observed that $\partial_{\bm{\uptau}}\bm{\Lambda}$ only has value for the forward transformation case while it has value zero for the reverse transformation case.

[TABLE]

to calculate the partial derivative of the current maximum transformation strain $H^{cur}$ with respect to stress $\bm{\uptau}$, the following result can be obtained based on equation (41),

[TABLE]

Follow the similar procedure to obtain $\partial_{\bm{\uptau}}\Phi$, the partial derivative of the transformation function $\Phi$ with respect to martensitic volume fraction $\xi$, and the partial derivative of the transformation function $\Phi$ with respect to temperature $T$ can be calculated as follows,

[TABLE]

[TABLE]