Analytical study on phase transition of shape memory alloy wire under uniaxial tension

TL;DR

This paper analytically investigates the stress-induced phase transition and instability in shape memory alloy wires under uniaxial tension, revealing detailed transition mechanisms and their dependence on material and geometric parameters.

Contribution

It introduces a simplified 1D analytical model for phase transitions in SMA wires, providing explicit solutions and insights into the transition process and associated instabilities.

Findings

Phase transition occurs over a small region, not at a single point.

The width of the transition region depends on material and geometric parameters.

Stress-strain curves show sharp drops and bifurcations consistent with experiments.

Abstract

This paper considers the stress-induced phase transitions of shape memory alloy slender cylinder, and analytically studies the phase transition process and the associated instability. A three-dimensional (3D) phenomenological model with an internal variable is adopted, which is simplified to a 1D system of two strains in three regions (austenite, martensite and phase transition regions). Suitable boundary conditions and interface conditions are proposed. Theoretically, it is a free boundary problem, as the position and shape of phase interfaces are unknown. We then focus on planar interfaces (which are energetically favored), and a symmetric case when phase transition occurs in the middle. For given applied stress, two-region solutions, three-region solutions and the connecting solutions between them are obtained analytically or semi-analytically, including many period-k solutions. Two-region solutions show that phase transition does not take place at one point, but simultaneously in a small region. The width of phase transition region is found analytically, revealing the roles of the material and geometrical parameters. Three-region solutions represent localized inhomogeneous deformations in experiments, and capture that the stress stays almost at the Maxwell stress during propagation of transformation front. For displacement-controlled process, the transition process is demonstrated by the stress-strain curve, which captures the stress drop/rise and the transition from homogeneous deformations to period-1 localized inhomogeneous deformations. When the radius is smaller than a critical value (given by material constants), the stress drop is very sharp due to transition of solutions in a snap-back bifurcation. These features show good agreement with experimental observations and shed light on the difficulties of numerical simulations.

Equations

$$\mathbf{\Sigma }=\frac{\partial \Phi (\mathbf{F},\alpha ,T)}{\partial \mathbf{F}},\xi =-\frac{\partial \Phi (\mathbf{F},\alpha ,T)}{\partial \alpha }\dot{\alpha },$$

$${\mathbf{G}}_{\alpha }:=(1-\alpha )\mathbf{I}+\alpha \mathbf{G},0\le \alpha \le 1,$$

$$\mathbf{G}=\mathrm{diag}(1-s_1,1-s_1,1+s_2),$$

$$\Phi (\mathbf{F},\alpha ,T)=$$

$$+B\alpha (1-\alpha )+[(1-\alpha ){\varphi }_1(T)+\alpha {\varphi }_2(T)],$$

$$\xi =\{\begin{array}{ll}{A}^{+}(\alpha )|\dot{\alpha }|=(k^{+}\alpha +Y^{+})\dot{\alpha }& \mathtt{if}\dot{\alpha }\ge 0,\\ A^{-}(\alpha )|\dot{\alpha }|=-(k^{-}(1-\alpha )+Y^{-})\dot{\alpha }& \mathtt{if}\dot{\alpha }\le 0,\end{array}$$

$$-A^{-}(\alpha )<-\frac{\partial \Phi }{\partial \alpha }<A^{+}(\alpha )\Rightarrow \dot{\alpha }=0,$$

$$\dot{\alpha }\ne 0\Rightarrow -\frac{\partial \Phi }{\partial \alpha }=\{\begin{array}{ll}{A}^{+}(\alpha )& \mathtt{if}\dot{\alpha }>0,\\ -A^{-}(\alpha )& \mathtt{if}\dot{\alpha }<0.\end{array}$$

$$r=R+V(R,Z),\theta =\Theta ,z=Z+W(R,Z),$$

$$\frac{\partial {\Sigma }_{Rr}}{\partial R}+\frac{\partial {\Sigma }_{Zr}}{\partial Z}+\frac{{\Sigma }_{Rr}-{\Sigma }_{\Theta \theta }}{R}=0,\frac{\partial {\Sigma }_{Rz}}{\partial R}+\frac{\partial {\Sigma }_{Zz}}{\partial Z}+\frac{{\Sigma }_{Rz}}{R}=0,$$

$${\Sigma }_{Rr}{|}_{R=a}=0,{\Sigma }_{Rz}{|}_{R=a}=0.$$

$${\Sigma }_{Zz}{|}_{Z=0,1}=\gamma ,{\Sigma }_{Zr}{|}_{Z=0,1}=0,$$

$$\tilde{V}(R,Z)=\sum _{k=0}^{\infty }\frac{{\tilde{V}}_k(Z)R^{2k}}{2(1+k)},W_Z(R,Z)=\sum _{k=0}^{\infty }{W}_{kZ}(Z)R^{2k},$$

$${\mathbf{U}}_0(Z)=[W_{0Z}(Z),{\tilde{V}}_0(Z){]}^T,$$

$$a^2{\mathbf{U}}_{0ZZ}+\mathbf{H}{\mathbf{U}}_0={\mathbf{f}}_0,\mathrm{in}\mathrm{AR},$$

$$a^2{\mathbf{U}}_{0ZZ}+\bar{\mathbf{H}}{\mathbf{U}}_0=\bar{\mathbf{f}},\mathrm{in}\mathrm{PTR},$$

$$a^2{\mathbf{U}}_{0ZZ}+\mathbf{H}{\mathbf{U}}_0={\mathbf{f}}_1,\mathrm{in}\mathrm{MR},$$

$$-{\mathbf{b}}^{(2)}{\mathbf{U}}_0{|}_{Z=0,1}=\gamma ,{\mathbf{b}}^{(2)}{\mathbf{U}}_{0Z}{|}_{Z=0,1}=0,\mathrm{in}\mathrm{AR},$$

$$-{\mathbf{b}}^{(2)}{\mathbf{U}}_0=\gamma +s_2{\xi }_1-2s_1{\xi }_2,{\mathbf{b}}^{(2)}{\mathbf{U}}_{0Z}=0,\mathrm{in}\mathrm{MR},$$

$$\mathrm{or}$$

$${\gamma }^{\ast }=\frac{\gamma }{E},D_k^{\pm }=\frac{2B-k^{\pm }}{E},\Delta \varphi ={\varphi }_2-{\varphi }_1,$$

$$D_{\varphi }^{+}=\frac{B+Y^{+}+\Delta \varphi }{E},D_{\varphi }^{-}=\frac{B-Y^{-}-k^{-}+\Delta \varphi }{E},$$

$${\tilde{V}}_{0Z}^A(Z_0)=$$

$$W_{0Z}^{A,P}(Z_0)=$$

$${\tilde{V}}_0^{A,P}(Z_0)=$$

$${\tilde{V}}_{0Z}^P(Z_2)=$$

$$W_{0Z}^{P,M}(Z_2)=$$

$${\tilde{V}}_0^{P,M}(Z_2)=$$

$${\mathbf{U}}_{0Z}(0.5)=0,\iff W_{0ZZ}(0.5)=0,{\tilde{V}}_{0Z}(0.5)=0,$$

$${\mathbf{U}}_0^A(Z)=$$

$${\mathbf{U}}_0^M(Z)=$$

$$\displaystyle\mathbf{U}_{0c}(Z;B_{1},B_{2},B_{3},B_{4})=B_{1}e^{-d_{1}\tilde{Z}}\left(\begin{array}[]{c}\cos(d_{2}\tilde{Z})\\ q_{1}\cos(d_{2}\tilde{Z})+q_{2}\sin(d_{2}\tilde{Z})\\ \end{array}\right)$$

$$\displaystyle+B_{2}e^{-d_{1}\tilde{Z}}\left(\begin{array}[]{c}-\sin(d_{2}\tilde{Z})\\ -q_{1}\sin(d_{2}\tilde{Z})+q_{2}\cos(d_{2}\tilde{Z})\\ \end{array}\right)$$

$$\displaystyle+B_{3}e^{d_{1}\tilde{Z}}\left(\begin{array}[]{c}\cos(d_{2}\tilde{Z})\\ q_{1}\cos(d_{2}\tilde{Z})-q_{2}\sin(d_{2}\tilde{Z})\\ \end{array}\right)$$

$$\displaystyle+B_{4}e^{d_{1}\tilde{Z}}\left(\begin{array}[]{c}\sin(d_{2}\tilde{Z})\\ q_{1}\sin(d_{2}\tilde{Z})+q_{2}\cos(d_{2}\tilde{Z})\\ \end{array}\right),$$

$${\mathbf{U}}_0^P(Z)={\mathbf{U}}_{0p}^P+C_5{\bar{\mathbf{q}}}_1{e}^{d_3\tilde{Z}}+C_6{\bar{\mathbf{q}}}_1{e}^{-d_3\tilde{Z}}+C_7{\bar{\mathbf{q}}}_2\cos (d_4\tilde{Z})+C_8{\bar{\mathbf{q}}}_2\sin (d_4\tilde{Z}),$$

$${\mathbf{U}}_{0p}^P={\left(\frac{D_{\varphi }^{\pm }{s}_2}{D_k^{\pm }}+\frac{{\gamma }^{\ast }\left(D_k^{\pm }-s_2^2\right)}{D_k^{\pm }},-\frac{2D_{\varphi }^{\pm }{s}_1}{D_k^{\pm }}-\frac{2{\gamma }^{\ast }\left(\nu D_k^{\pm }-s_1{s}_2\right)}{D_k^{\pm }}\right)}^T,$$

$${\bar{\mathbf{q}}}_1=[1,{\bar{q}}_1{]}^T,{\bar{\mathbf{q}}}_2=[1,{\bar{q}}_2{]}^T.$$

Full text

Two-region solutions show that phase transition does not take place at one point, but simultaneously in a small region. The width of phase transition region is found analytically, revealing the roles of the material and geometrical parameters. Three-region solutions represent localized inhomogeneous deformations in experiments, and capture that the stress stays almost at the Maxwell stress during propagation of transformation front. For displacement-controlled process, the transition process is demonstrated by the stress-strain curve, which captures the stress drop/rise and the transition from homogeneous deformations to period-1 localized inhomogeneous deformations. When the radius is smaller than a critical value (given by material constants), the stress drop is very sharp due to transition of solutions in a snap-back bifurcation. These features show good agreement with experimental observations and shed light on the difficulties of numerical simulations.

1 Introduction

As a smart material, shape memory alloys (SMAs) such as Ni-Ti have been utilized for various real-world applications in aerospace, automotive and biomedical industries, and oil exploration [1]. One impressive example is the Ni-Ti self-expanding stent used in the minimally invasive surgery of heart disease. Many applications are motivated by two remarkable properties of SMAs: shape memory effect and pseudo-elasticity. This paper focuses on pseudo-elasticity, which means that SMAs can undergo large recoverable strain without permanent damage during the loading-unloading cycle at relatively high temperatures [2]. The mechanism of this property is closely related to the stress-induced martensitic phase transitions [3].

SMAs have attracted great research interests in the past few decades, including both experiments and theoretically modelling. The macroscopic behaviour of SMAs, including pseudo-elasticity, are characterized by stress-strain curves, and a large number of experiments [3, 4, 5, 6, 7] have been conducted on SMA wires at different temperatures under tension. For a class of SMAs (such as Ni-Ti), many experiments show strain softening behaviors in stress-strain (stress-elongation) curves. As a result, phase transitions are realized by nucleation of martensite band and subsequent propagation of transformation fronts along the wire [8, 4, 9], also called Lüders-like behavior. The distinction and roles of nucleation and propagation stresses have been observed in experiments and analytically studied in [10]. At nucleation of a new phase, initial homogeneous deformations will be replaced by localized inhomogeneous deformations, where the strain varies rapidly across the transformation front or the phase transition region (PTR). Therefore, capturing the transition process and how the geometrical and material parameters affect the process are of fundamental importance both theoretically and practically. However, systematic analytical investigations on the transition process and the associated instability under the three-dimensional (3D) framework are few, due to the complexity of the governing equations and the presence of free boundaries.

Much research effort has been devoted to the constitutive models for the phase transition of SMAs. At the microscopic level, the phase field (Ginzburg-Landau) models [11, 12, 13, 14] concentrate on the formation and evolution of martensitic microstructure, which are described by a set of kinetic equations for order parameters. They have shown success in dealing with many variants of martensite phase simultaneously and capturing the sharp interfaces. Nevertheless, the specimen treated in experiments are relatively large (e.g. 10-100 mm) for pseudo-elasticity, and so the size of PTR is much larger than that simulated in those models. It seems more suitable to investigate the above Lüders-like behavior at the macroscopic level.

At the macroscopic level, the models can be put into two categories: sharp interface models and diffuse interface models. Some thermoelastic models [15, 16, 17, 18] treat the phase interface as a sharp singular surface, where further conditions such as nucleation criterion and kinetic relation are proposed. For simulations with moving phase interface, see also [19]. In diffuse interface models, the transformation front is often identified with a small PTR, which smoothly connects the low-strain austenite phase and high-strain martensite phase with large gradients. Some models [20, 21, 22, 23, 24] include extra gradient terms (e.g., gradient of strain), whose coefficient sets a characteristic length scale for the width of PTR. One disadvantage might be that it is not an easy task to determine the artificial coefficient beforehand. The phenomenological models [25, 26, 27] introduce the martensite volume fraction as an internal variable into the Helmholtz free energy together with a dissipation function, then both the strain and martensite volume fraction are continuous across PTR.

On the continuum level, the transition from homogeneous deformations to localized inhomogeneous deformations is due to the material instability. This is reflected by the existence of multiple solutions at certain stress or elongation, leading to the transition to a more stable solution by the energy criterion. The current work adopts phenomenological model by [25]. One advantage of such a model is that all the parameters have clear physical meaning (no artificial coefficient) and the influence of each parameter on the transition process can then be easily identified.

Many analytical results in the literature are based on 1D models without strain softening behaviour, and investigated stress-strain-temperature response of various structures (e.g., wires [28, 29], layer structures [30] and actuators [31, 32, 33]). However, the high-dimensional effect and strain softening behaviour are crucial to capturing the stress drop and determining the width of PTR during phase transition, therefore it is desirable to take these into account. Under 3D frameworks, some analytical results [34] are obtained for thick-walled cylinder subjected to external pressure, but with Tresca transformation criterion and linear hardening. For an infinitely-long cylinder, some analytical works [35, 36] treat the martensite band as inclusion in a 3D inclusion-matrix system with planar sharp interfaces, capturing some important features such as peak and propagation stresses.

Another class of analytical study was first attempted in [37, 38] by adopting a non-convex hyperelastic 3D energy and reducing it to a 1D rod model. It has been generalized to study the geometric size effect and hysteresis loop in [39, 40] for SMA wires. By combining with the model in [25], it has been extended to study SMA layer [41, 42] and geometrically graded SMA strip [43, 44]. Our previous work [10] has studied piecewise homogeneous deformations of an SMA wire, and has obtained analytical results for peak and propagation stresses [45]. Subsequently for a 3D infinitely long cylinder [46], we have studied localized inhomogeneous deformations by keeping high-dimensional effects, capturing the propagation stress plateau. The present work aims to investigate the phase transition of an SMA wire with finite length, which is a natural continuation of infinite case [46]. This is not an easy task, since one has to deal with free boundaries in response to applied stress or total elongation. In this work, we aim to systematically study the transition process from the homogeneous deformation to the localized inhomogeneous deformation, and to reveal the roles of geometric and material parameters in this process.

We start from a 3D formulation [25], and the mechanical system is simplified to a 1D system for three different regions, i.e. austenite region (AR), martensite region (MR), and PTR, by utilizing coupled series-asymptotic technique [47, 48]. Since the shape and positions of interfaces between different regions are not known, mathematically this is a tricky free boundary problem. We concentrate on the planar interfaces, which are energetically favored (see [46]) and make analytical results available. Then we consider the symmetric case, where phase transition takes place in the middle of the cylinder. First, two-region solutions with AR and PTR are obtained analytically, with many period-k solutions. Second, for three-region solutions with AR, PTR and MR, determining the solutions is reduced to a system of two nonlinear algebraic equations. The case of infinitely long cylinder in [46] can be recovered as a special case of such solutions. Third, connecting solutions are constructed for the gap between two-region and three-region solutions, which is due to the interaction between phase interface and middle (boundary) surface. Finally, the stress-elongation curves are obtained based on these solutions, and the transition process is captured by energy criterion.

Finally, we briefly mention the key findings of the analytical study. First, analytical solutions show that phase transition does not take place at one point, but simultaneously in a small region. The analytical formulas for the width of PTR clear reveal the roles of the material and geometrical parameters. Second, three-region solutions capture the localized inhomogeneous deformations, indicating that the stress always stays close to Maxwell stress in the stress-elongation curve. Third, the transition from homogeneous deformations to localized inhomogeneous deformations is clearly demonstrated at nucleation of a new phase. The stress drop or rise is captured in stress-elongation curves, and the roles of parameters are disclosed. For instance, as the radius of cylinder becomes smaller than a critical value (given by material constants), there is a sharp stress drop as observed in experiments. As multiple solutions exist for certain fixed total elongation, it provides insights into the difficulties of direct numerical computations [49], e.g., mesh sensitivity and convergence difficulty.

The manuscript is arranged as follows. Section 2 briefly recalls 3D formulation and its simplification to a 1D system in three different regions (AR, MR and PTR). In Section 3, a symmetric case for a finite cylinder is considered. In Sections 4 to 6, the two-region solutions, three-region solutions and connecting solutions are constructed, either analytically or semi-analytically. Then in Section 7, the phase transition process is investigated through stress-elongation curves in displacement-controlled process. Finally, concluding remarks are made.

2 Problem formulation and simplification

As this work is a continuation of previous work in [46], for completeness this section is a brief re-account of the problem setup in [46], with necessary modifications and remarks. The previous work considers the case of an infinitely long cylinder, while this work studies the case of finite cylinder. Thus, the boundary conditions are different, and the aim now is to study the transition process from homogeneous deformations to localized inhomogeneous deformations and the associated instability. Another complication in current work is that during the transition process the number and location of interfaces (regions) are not known beforehand, so it is a free boundary problem. As we will see later that, the infinitely long case [46] can be considered as a special (and idealized) case of current one for the three-region localized inhomogeneous deformations.

2.1 Constitutive model

This subsection briefly recalls the constitutive model proposed by Rajagopal and Srinivasa in [50, 25], used in previous analytical works [41, 46]. The model is a phenomenological model with an internal variable defined as the volume fraction $\alpha$ of martensite phase. It is described by two independent functions: Helmholtz free energy $\Phi(\mathbf{F},\alpha,T)$ and rate of mechanical dissipation $\xi(\alpha)$, where $\mathbf{F}$ is deformation gradient and $T$ is the absolute temperature.

This work focuses on the isothermal case. By the first and second laws of thermodynamics, we have the following relations

[TABLE]

where $\mathbf{\Sigma}$ is the nominal stress tensor (or the transpose of first Piola-Kirchhoff stress). The specification of $\Phi$ relies on the natural configuration associated with phase state $\alpha$, defined through deformation gradient tensor

[TABLE]

where diag is a diagonal matrix in cylindrical polar coordinates, $\mathbf{G}$ is the natural configuration of martensite phase, and $s_{1}$ and $s_{2}$ (usually $\ll 1$) are two material parameters, called the specific strains normal to and along the axis of the cylinder. The elastic part of the deformation gradient $\mathbf{F}$ is $\mathbf{F}_{\alpha}=\mathbf{F}\mathbf{G}_{\alpha}^{-1}$, and the Helmholtz free energy per reference volume is given by

[TABLE]

where $\Phi_{1}$ and $\Phi_{2}$ are respectively the strain energy functions of the austenite and martensite phases, $B$ is interfacial constant [51], and $\phi_{1}$ and $\phi_{2}$ are respectively the thermal free energies of the austenite and martensite phases.

The dissipation function $\xi$ is adopted as

[TABLE]

where $\pm$ correspond to loading and unloading processes respectively, and $k^{\pm}$ and $Y^{\pm}$ are referred to as the dissipative constants. Substituting (4) into $(\ref{eq1})_{2}$ leads to the evolution law of $\alpha$

[TABLE]

Here $\dot{\alpha}=0$ means that the material response is elastic, e.g., when the material initially undergoes elastic deformation in pure austenite phase ($\alpha=0$). When phase transition occurs, the variable $\alpha$ in PTR is determined by $(\ref{eq10})_{2}$.

**Remark: **The interfacial energy $B\alpha(1-\alpha)$ in (3) leads to a non-convex energy and is closely related to the strain-softening behaviour and the associated instability. Similar terms of $\alpha$ were introduced in other works [e.g., 52, 27] to account for mixing effect, resulting in softening response of the material. The terms in dissipation (4) will also affect the strain-softening, and the constants $k^{\pm}$ have competing effects with $B$.

2.2 The 3D mechanical system

Figure 1 shows the slender cylinder (wire) with finite length under uniaxial average stress $\gamma$. Without loss of generality, the length of cylinder is set to be $L=1$, so the dimensionless radius $a$ here is interpreted as ratio of radius to length, and all the displacements below are already scaled by $L$. For a slender cylinder, we have $a\ll 1$. The reference configuration is entirely made up of austenite phase. The cylindrical polar coordinates are adopted, and $(R,\Theta,Z)$ and $(r,\theta,z)$ denote a generic point in the reference and current configurations, respectively. The general deformation can be described by

[TABLE]

where $V$ and $W$ are radial and axial displacements. The deformation gradient $\mathbf{F}$ is easily calculated by taking the derivatives, and the stress tensor $\mathbf{\Sigma}$ can be obtained from $(\ref{eq1})_{1}$.

We consider a quasi-static process for the phase transition of SMAs, as in some displacement-controlled experiments. By neglecting the body force, the mechanical field equations are

[TABLE]

where subscripts denote the components of stress along $R,\theta,Z$ directions. Next we consider the boundary conditions, which are different from [46]. First, the traction-free boundary conditions on the lateral surface are

[TABLE]

Second, on the end surfaces of the cylinder, we have

[TABLE]

where $\gamma$ is the average stress. For a force-controlled process, $\gamma$ is a given parameter, whereas for a displacement-controlled process as in many experiments, the parameter $\gamma$ is to be determined by the total displacement of cylinder.

The 3D mechanical system (7-9) for the two displacements $(V,W)$ will be coupled with the phase state variable $\alpha$, since $\mathbf{\Sigma}$ involves $\alpha$. The extra phase state $\alpha$ is determined by the evolution law $(\ref{eq10})$.

2.3 The 1D governing equations

This subsection presents the simplified 1D equations in different regions based on the preceding 3D system. The derivation is based on the series-asymptotic technique in [37, 38], which takes advantage of two small parameters: the dimensionless radius $a$ and the scale of typical axial strain. The details are given in [46], and here we present a simple modified version of the derivation and final system, with necessary important remarks.

First, we define the radial strain by $\tilde{V}=V/R$, and take the series expansion for the two strains

[TABLE]

where and hereafter subscript $Z$ denotes derivative with respect to $Z$. In the derivation, we only keep the leading order terms of the typical axial strain. Recursive relations can be derived from (7) for high-order coefficients (terms with $k\geq 1$) in (10) in terms of leading-order coefficient. They are different in three different regions: austenite region (AR, $\alpha=0$), martensite region (MR, $\alpha=1$) and phase transition region (PTR, $0<\alpha<1$). For instance, for PTR in loading process, $\alpha$ is first determined by $-\partial\Phi/\partial\alpha=A^{+}$, and then it is substituted into (7) to derive such recursive relations.

**Remark: **The series expansion (10) in terms of $R^{2}$ is due to the symmetry of the problem, since the terms with odd powers naturally disappear. The factor $2(1+k)$ in $(\ref{eq17})_{1}$ is to make all recursive relations in simple and compact form [46], otherwise the expressions of recursive relations are more complicated (see, e.g., [53]).

With recursive relations, the conditions in (8) provide the 1D governing equations for the leading-order strain vector (axial and radial strains)

[TABLE]

where superscript $T$ means transpose of a vector. By truncating the system and keeping upto $O(a^{2})$ terms, the governing equations for $\mathbf{U}_{0}$ in three regions take the form

[TABLE]

where the matrices $\mathbf{H},\bar{\mathbf{H}}$ and column vectors $\mathbf{f}_{0},\mathbf{f}_{1},\bar{\mathbf{f}}$ are given in Appendix A. The associated boundary conditions at $Z=0,1$ can be obtained from (9) by consistent truncations. The explicit form of boundary conditions depends on the region where the end point $Z=0$ or $Z=1$ is located. If both ends are in AR, we have

[TABLE]

where the row vector $\mathbf{b}^{(2)}$ is given in Appendix A. If the end point lies in other regions, either MR or PTR, the conditions should be modified as

[TABLE]

where the elastic moduli $\xi_{1},\xi_{2}$, the row vector $\bar{\mathbf{b}}^{(2)}$ and constant $\alpha_{00}$ are given in Appendix A.

Since the above simplified 1D system for $\mathbf{U}_{0}$ may consist of multiple regions, one needs to consider the interface between different regions and propose suitable connection conditions on interfaces, for which we turn to the next subsection.

**Remark: **There are many situations for the 1D system during the phase transition process. First, it may contain just one region, e.g., AR, this is the case when $\gamma$ is small with pure elastic deformations in austenite phase. Then the equation $(\ref{eq23})_{1}$ and condition (13) form the complete system, which leads to the unique solution of $\mathbf{U}_{0}$ corresponding to the homogeneous deformation. Second, for certain stress $\gamma$, the cylinder would contain multiple phase regions, e.g., AR and PTR in $(\ref{eq23})$. In this case, one needs to solve the associated system with interface conditions. Mathematically this is a free boundary problem, as the shape and position of interface are unknown and even the number of regions/interfaces is unknown beforehand.

2.4 Connection conditions on interface

During the phase transition process and for certain stress $\gamma$, there exist multiple phase regions in the cylinder. The interface between different regions is a free boundary. In general the shape of interface is unknown, and the general interfaces have been investigated in [46] and it turns out that the planar interface is favored (from an energy perspective). Planar interfaces are also adopted by others in 3D analytical studies [35, 36]. Here, we restrict ourselves to planar interfaces, for both simplicity and availability of analytical results, which provide more insights into the mechanism of phase transition and roles of parameters.

Figure 2 schematically illustrates two planar interfaces between three regions in the reference configuration. Interfaces $I_{AP}$ and $I_{PM}$ correspond to the surface $\alpha=0$ and $\alpha=1$, where the positions $Z_{0}$ and $Z_{2}$ of the interfaces are to be determined. For second-order vector equations in (LABEL:eq23) and considering the two unknown positions, one needs five connection conditions on each interface. These conditions can be obtained based on continuity of phase state $\alpha$ and deformation gradient $\mathbf{F}$ (see Section 5.2 of [46] for details).

To present the conditions and other expressions, we first define

[TABLE]

where $E$ and $\nu$ are Young’s modulus and Poisson’s ratio. Then, $\gamma^{\ast}$ is the dimensionless nominal stress, and others are combinations of material parameters.

On the interface $I_{AP}$, the five connection conditions are

[TABLE]

where superscripts $A$ and $P$ denote solutions in AR and PTR respectively. On the interface $I_{PM}$, we have

[TABLE]

where superscript $M$ denotes solutions in MR.

In short, depending on the number of interfaces in the cylinder, one needs to suitably choose the connection conditions on either $I_{AP}$ or $I_{PM}$ or both, to complete the system in (LABEL:eq23-LABEL:3eq5).

**Remark: **In this paper, interface means the boundary between two different regions (AR, PTR or MR), which is different from the phase interface in sharp interface models. In the present diffuse interface model, the transformation front is identified with a small region (PTR) rather than a sharp interface. This is similar to models with strain gradient terms [20, 21], with a parameter for the width of a smooth PTR. In fact, the present local 3D model without strain gradient terms naturally gives a width of PTR, determined by the material and geometric parameters.

3 Symmetric case and general solutions

In the previous section, we have obtained the simplified 1D system, with possible interface conditions. Since one does not know the number and locations of regions beforehand, the governing equations and boundary conditions in the system will vary. Mathematically speaking, this is a tricky free boundary problem. The instability during phase transition is reflected by the existence of multiple solutions and the transition to a more stable solution. We recall the phase transitions process in experiments and consider a symmetric case. Then, we present the general solutions in each region, with integrating constants to be determined.

In the following we take the loading process for illustration, and unloading process can be similarly studied with slight modifications. From experiments, initially the SMA wire undergoes elastic homogeneous deformations in austenite phase. This is the case when the cylinder contains only one phase region AR, and it is easy to prove that this is a trivial case that the solution is unique, corresponding to the homogeneous deformation. At certain critical stress, the nucleation of martensite phase occurs, i.e., the transition from homogeneous deformations to inhomogeneous deformations takes place. Finally during propagation of martensite band, it is in a localized inhomogeneous deformation. Clearly, the localized inhomogeneous deformation corresponds to a case with three regions AR, PTR and MR. In between, there should be transient state of inhomogeneous deformations with two regions AR and PTR, and this state could be either energetically stable or unstable.

Nucleation of martensite phase can occur either at end point or in the middle of cylinder, depending on the imperfection of material and experimental setup at two ends. To avoid the complicated interaction between phase interface and end surface [54], we will focus on the symmetric case, where the martensite band nucleates in the middle of the cylinder. In this context, we concentrate on the interval [0, 0.5], and adopt symmetric conditions at the middle surface $Z=0.5$ in place of the end condition at $Z=1$. Since the solution is symmetric about $Z=0.5$, it is natural to adopt

[TABLE]

i.e., the two leading-order strains (the axial strain $W_{0Z}$ and the radial strain $\tilde{V}_{0}$) are symmetric. Actually, this is equivalent to the symmetry of profile of the cylinder, justified in the next section. As a result, we only need to examine the left half of the cylinder with (18) and the other condition (13) or (LABEL:3eq5) at $Z=0$.

In the following, we present the general solutions in each region of (LABEL:eq23). They are constructed by finding a particular solution and studying the eigenvalues of matrix $\mathbf{H}$ or $\bar{\mathbf{H}}$. The general solutions in AR ($\alpha=0$) and MR ($\alpha=1$) are given by

[TABLE]

where all $C_{i}$ ($i=1,..,4,9,..,12$) are integrating constants and

[TABLE]

where $\tilde{Z}=Z-Z_{0}$, and the explicit expressions of $d_{1},d_{2},q_{1},q_{2}$ are given in Appendix B. In the solution, we have made a shift in $Z$ by $Z_{0}$ (any shift will be a solution), which makes it easier to determine the integrating constants in later sections. The general solution in PTR is

[TABLE]

where $C_{i}$ ($i=5,..,8$) are integrating constants and

[TABLE]

The constants $d_{3},d_{4},\bar{{q}}_{1},\bar{{q}}_{2}$ are given in Appendix B.

Based on the particular solutions ($C_{i}=0$) in (19, 21, LABEL:eq60), Figure 3 shows the 1D stress-strain (axial strain) relation for homogeneous states of the loading and unloading processes, where the Maxwell stresses $\gamma_{m}^{\pm}$ are defined by

[TABLE]

Figure 3 shows the strain-softening behaviour, with reasonable material constants in Appendix C. Each point in Figure 3 represent a solution for homogeneous deformation. These are the unique solutions when there is only one region, AR, PTR or MR at different stress $\gamma^{\ast}$. However, the homogeneous solution with pure PTR region, corresponding to the softening branch in the curve, is locally unstable from an energy point of view [45]. This part of stress-strain curve will not be observed in experiments. In the next sections, we will focus on inhomogeneous deformations where multiple regions are present, and investigate the phase transition process. For inhomogeneous deformations, we restrict ourselves to the possibly smallest number of interfaces and regions, as adding more regions and interfaces will make solutions unfavored energetically. Since the general solution is given in each region, determining the solutions for inhomogeneous deformations is equivalent to determining the associated integrating constants.

4 Two-region solutions

In this section, we study the case when the cylinder contains only two regions in $[0,0.5]$, AR and PTR from left to right (see Figure 4). We first present the analytical results for the two-region solutions for given stress $\gamma^{\ast}$, and then calculate the total elongations for the solutions, which are often measured in displacement-controlled experiments on SMAs .

4.1 Analytical results for given $\gamma^{\ast}$

In this case, the general solution $\mathbf{U}_{0}$ in $[0,0.5]$ is composed of $(\ref{3eq7})_{1}$ and (21), which contains eight integrating constants $C_{1},..,C_{8}$ and one unknown position $Z_{0}$ from the interface $I_{AP}$. The two boundary conditions in (13) at $Z=0$, two symmetric conditions in (18) at $Z=0.5$ and five connection conditions $(\ref{eq87})$ at interface $Z=Z_{0}$ can totally determine the nine unknowns.

First, the integrating constants $C_{1},..,C_{8}$ can be expressed in terms of $Z_{0},\gamma^{\ast}$ and material parameters. Since the transition region is in the middle of the cylinder, the terms containing $C_{3}$ and $C_{4}$ (with factor $e^{-d_{1}Z_{0}}$) in solutions $(\ref{3eq7})_{1}$ at the left end are exponentially small. As a result, the two boundary conditions (13) at $Z=0$ lead to $C_{1}=C_{2}=0$. Then, substituting $(\ref{3eq7})_{1}$ into conditions (16) leads to

[TABLE]

and the conditions in $(\ref{eq87})_{2,3}$ and (18) imply

[TABLE]

where $\bar{Z}_{0}=0.5-Z_{0}$ is the width of PTR in this case, and the material constants $q_{11},q_{22},\bar{q}_{11},\bar{q}_{22}$ are defined in Appendix B.

**Remark: **From general solutions of $\mathbf{U}_{0}$ together with the above results and recursive relations, we have the following equivalence

[TABLE]

Thus, the symmetric conditions of the two strains $W_{0Z}$ and $\tilde{V}_{0}$ are equivalent to the symmetric conditions of the profile (radial displacement) at all points on the middle surface. This equivalence also holds for the three-region solutions in the next section, so there is no need to distinguish the two sets of symmetric conditions.

Next, by substituting the above expressions of $C_{1},..,C_{8}$ into the last condition $(\ref{eq87})_{1}$, we obtain an equation for $\bar{Z}_{0}$ (note $\bar{Z}_{0}=0.5-Z_{0}$) in the form

[TABLE]

where the material constants $n_{11},n_{12},n_{13}$ are independent of $\gamma^{\ast}$ and are given in Appendix B. Thus for $\gamma^{\ast}<\gamma_{NM}^{\ast}:=D_{\phi}^{+}/s_{2}$, the unknown $\bar{Z}_{0}$ or interface position $Z_{0}$ is determined by the nonlinear equation in curly brackets in (27) and does not depend on the given parameter $\gamma^{\ast}$. Therefore, determining the two-region solutions has reduced to solving a nonlinear equation for $\bar{Z}_{0}$.

In order to get closed-form solutions of $\bar{Z}_{0}$ (and hence $\mathbf{U}_{0}$), we neglect the exponentially small terms of $O(e^{-d_{3}\bar{Z}_{0}})$ in (27) and obtain

[TABLE]

where $\bar{Z}_{0}^{(k)}$ is the $k^{th}$ solution for $\bar{Z}_{0}$, and $-\pi/2<\arctan x<\pi/2$. This formula is numerically verified by solving the whole system of $C_{1},..,C_{8},Z_{0}$, with chosen material constants in (63) and $a=0.03$. The first few solutions are easily obtained by Newton’s method, and the solution with smallest width $\bar{Z}_{0}$ is given in Appendix C. By comparison with the direct numerical results, the above analytical approximation (28) is very accurate, e.g.,

[TABLE]

The relative error of $\bar{Z}_{0}^{(k)}$ is almost 0 for any $k\geq 1$. For a given $\gamma^{\ast}$ in certain interval, each $\bar{Z}_{0}^{(k)}$ determines one solution of $\mathbf{U}_{0}$, called the period-k solution.

Finally, based on the result (28), we analyze the dependence of the width of PTR on material and geometrical parameters. By using the small parameter $\epsilon_{1}:=D_{k}^{+}/s_{2}^{2}$ and retaining two leading terms, we obtain

[TABLE]

Clearly, the width $\bar{Z}_{0}^{(k)}$ for any $k$ is proportional to radius $a$, and increases with $k$. The proportional dependence on geometric parameter $a$ is expected in both experiments and theoretical works [21, 23]. With $k=1$ and the condition $s_{1}=s_{2}/2$ (volume preserving condition for the natural configuration $\mathbf{G}$ in (2)), it reduces to

[TABLE]

It depends on material parameters mainly through $\epsilon_{1}$ with leading power $-1/4$, whereas the Poisson’s ratio $\nu$ only has minor effect. Thus to leading order, the width $\bar{Z}_{0}^{(1)}$ is proportional to the square root of axial transformation strain $\sqrt{s_{2}}$, and it becomes larger as $D_{k}^{+}$ decreases. From the experiments reported in [55], the width of PTR (transformation front) is closely related to the stress drop of the stress-strain curve, and tends to increase as the stress drop decreases. Here the term $D_{k}^{+}$ defined in (LABEL:eq52) is associated with stress drop or strain softening behavior, and is determined by the competition between the interfacial energy (constant $B$) and the dissipation (constant $k^{+}$). It is shown in [10] that stress drop decreases with decrease of $D_{k}^{+}$. Therefore, the qualitative feature of the analytical formula agrees well with experiments. In addition, the appearance of this natural width is due to both the strain-softening behavior and the geometric effect, since the high-dimensional effect is retained in the formulation (hidden in recursive formulas).

4.2 Period-1 solution and the total elongation

In the previous subsection, period-k solutions for given $\gamma^{\ast}$ are obtained. Actually, for a displacement-controlled process as in experiments, the parameter $\gamma^{\ast}$ is to be determined by a given total elongation based on energy criterion. One can expect that the solution with larger width $\bar{Z}_{0}$ would have larger total energy (defined in (LABEL:eq76), see also [46]). This subsection then focuses on the period-1 solution with smallest width $\bar{Z}_{0}^{(1)}$, and presents the total elongation of such solution with varying $\gamma^{\ast}$.

First, we present the solutions of the leading-order axial strain $W_{0Z}$ and phase state $\alpha$ for period-1 solution and some given $\gamma^{\ast}$. They can be represented as

[TABLE]

where $\hat{f}_{w}$ and $\hat{f}_{\alpha}$ are some shape functions independent of $\gamma^{\ast}$ and $D_{\phi}^{+}$ (recall $\gamma^{\ast}_{NM}=D_{\phi}^{+}/s_{2}$), whose expressions are omitted for brevity. The period-1 two-region solutions are obtained only for $\gamma^{\ast}$ in certain interval, since the phase state is restricted by $0<\alpha<1$, which leads to the range

[TABLE]

For the chosen material constants in (63) we have

[TABLE]

Figure 5 shows the profiles of $W_{0Z}(Z)$ and $\alpha(0,Z)$ in the interval $[0.3,0.5]$ for varying $\gamma^{\ast}$, indicated near the curves. The red curves represent the solutions in PTR, and blue curves represent solutions in AR and extend as a constant until left end $Z=0$. It is shown that the width of PTR is fixed independent of $\gamma^{\ast}$, while the maximum values of the curves at $Z=0.5$ increases with decease of $\gamma^{\ast}$. As $\gamma^{\ast}$ varies, the shape of the curves for $W_{0Z}$ or $\alpha$ in PTR remains the same, characterized by the function $\hat{f}_{w}$ or $\hat{f}_{\alpha}$.

The fixed width $\bar{Z}_{0}^{(1)}$ of PTR is an important feature of two-region solution, and its dependence on parameters has been discussed in the previous subsection. As we will see in Section 7 that the solution will transition from a homogeneous deformation (solution) to a period-1 two-region solution during phase transition. This implies that phase transition always takes place simultaneously in a fixed PTR region of width $\bar{Z}_{0}^{(1)}$ rather than at a point. Importantly, this is obtained by naturally keeping high-dimensional effect in the 3D local model instead of introducing strain gradient terms [21].

The total elongation is defined by

[TABLE]

By (LABEL:3eq52) it can be expressed in the form

[TABLE]

where $K_{2}$ only depends on material constants (not on $a$). The explicit expression of $K_{2}$ is omitted for brevity, but its approximation is given below. For the chosen material values in (63), we have $K_{2}=76.6342$. Figure 6 shows the stress-elongation curves of two-region solutions for three different radii $a$ in red lines, where the reference black line corresponds to homogeneous deformation in pure austenite phase (cf. Figure 3). Clearly, each curve is a line, passing through the peak point $(\gamma^{\ast}_{NM},\gamma^{\ast}_{NM})$. As the radius $a$ tends to 0, the red line of the two-region solution tends to the black line $\gamma^{\ast}=\Delta$. As the radius $a$ increases, the line will rotate to the right with the peak point fixed. There is a critical radius $a_{c}=1/K_{2}$ (between $0.01$ and $0.02$ here) such that the red line is vertical, which is crucial for the stability later.

To get an idea of the roles of material constants on elongation, we expand $K_{2}$ in terms of the small parameter $\epsilon_{1}$ defined above (30). With $s_{1}=s_{2}/2$, we get the two-term approximation

[TABLE]

which means $K_{2}$ depends on material parameters mainly through $\epsilon_{1}$ with leading power $-5/4$, while Poisson’s ratio $\nu$ has little impact. This implies $a_{c}\sim\epsilon_{1}^{5/4}$ to leading order. In particular, as $D_{k}^{+}$ decreases, $K_{2}$ will increase and thereby $K_{1}$ will decrease by (36) for fixed $a$, thus the line of two-region solution in Figure 6 will rotate to the right with fixed peak point.

**Remark: **As mentioned in the previous subsection, $D_{k}^{+}$ is crucial parameter associated with stress drop or strain softening behavior. As $D_{k}^{+}$ tends to 0, $K_{2}$ tends to $\infty$ and the red line of the two-region solution in Figure 6 tends to a horizontal line $\gamma^{\ast}=\gamma^{\ast}_{NM}$, together with an infinite PTR. This agrees very well with experiments in [55]. This makes sense since if there is no strain softening behaviour, there is no need to study the two-region solutions and the horizontal line is just the critical case of the stable homogeneous state with only PTR (the second branch of stress-strain curve in Figure 3 is no longer decreasing).

5 Three-region solutions

In this section, we investigate the case with three regions, i.e., AR, PTR and MR in the left half of the cylinder (see Figure 7). This corresponds to the localized inhomogeneous deformations during propagation of phase transformation front in experiments. The previous work in [46] is a special (ideal) case when the MR region is relatively large or PTR is far away from middle surface. The results here are much richer than the previous one, and the solutions are quite sensitive with small MR. We first present the semi-analytical results for the solution with given $\gamma^{\ast}$, and then compute the total elongation.

In this case, the general solution $\mathbf{U}_{0}$ in $[0,0.5]$ is composed of (19) and (21), which contains twelve integrating constants $C_{1},..,C_{12}$ and two unknown positions $Z_{0}$ and $Z_{2}$ for interfaces. The boundary conditions (13) at $Z=0$, symmetric conditions (18) at $Z=0.5$ and ten connection conditions (16) and (17) at interfaces can totally determine the unknowns.

The integrating constants $C_{1},..,C_{12}$ can be expressed by $\gamma^{\ast},Z_{0},Z_{2}$ and material constants. The constants $C_{1},..,C_{4}$ are the same as the previous section (see (24)). By the expressions $(\ref{eq87})_{2,3}$ and $(\ref{eq87_1})_{2,3}$, the constants $C_{5},..,C_{8}$ are given by

[TABLE]

where $\bar{Z}_{2}=Z_{2}-Z_{0}$ is the width of PTR in this case, and the material constants $\bar{q}_{11},\bar{q}_{22}$ are given in Appendix B. One can similarly obtain the expressions of $C_{9},..,C_{12}$, given in Appendix B.

Then for a given $\gamma^{\ast}$, the system with fourteen unknowns has reduced to a system of two unknowns $\bar{Z}_{2}$ and $\bar{Z}_{0}$ (recall $\bar{Z}_{0}=0.5-Z_{0}$), determined by the two conditions $(\ref{eq87})_{1}$ and $(\ref{eq87_1})_{1}$. First from $(\ref{eq87})_{1}$, the width $\bar{Z}_{2}$ of PTR is determined by

[TABLE]

where

[TABLE]

and $n_{11},n_{12},n_{13}$ are given in Appendix B. Numerical evidence shows that there exist solutions only for $\gamma^{\ast}$ in a small neighbourhood of $\gamma_{m}^{+}$. For the special case $\gamma^{\ast}=\gamma_{m}^{+}$, the equation (39) reduces to $\hat{n}_{1}(\bar{Z}_{2})=0$, which is exactly the infinite case in [46]. For this special case, we recall that the approximate period-k solution is given by

[TABLE]

where $-\pi<\mathrm{arccot}\,x<0$. It is verified by direct numerical results [46] that this formula is very accurate, with less than 1% error for all $k$.

For the general case $\gamma^{\ast}\neq\gamma_{m}^{+}$, there are singularities in $f_{2}(\bar{Z}_{2};\gamma^{\ast})$ due to the term $\tan(d_{4}\bar{Z}_{2}/2)$ in $\hat{m}_{1}$. There are two situations. First, one set of period-k solutions is just a small perturbation of (41), since $\gamma^{\ast}$ is close to $\gamma_{m}^{+}$. One can easily derive the regular period-k solutions

[TABLE]

Second, there are extra period-k solutions for $\bar{Z}_{2}$ near the singularities, namely,

[TABLE]

Figure 8(a) shows the curve of $f_{2}$ with $\gamma^{\ast}=0.00819$, $a=0.03$ and the material constants in (63), indicating that the first singularity is near $0.055$. The two roots near $0.055$ and $0.06$ correspond to the two period-1 solutions of $\bar{Z}_{2}$ with $k=1$ in (42) and (43). Figure 8(b) shows that the singularity disappears for the special case $\gamma^{\ast}=\gamma_{m}^{+}$.

Next by $(\ref{eq87_1})_{1}$, the equation for determining $\bar{Z}_{0}$ can be written as

[TABLE]

where $t:=\bar{Z}_{0}-\bar{Z}_{2}$ is the width of MR, and the functions $\hat{m}_{3}(t)$ and $\hat{n}_{3}(t)$ are given by

[TABLE]

It is not straightforward to obtain analytical solutions of $\bar{Z}_{0}$ from (44), instead we discuss the typical features of this equation for given $\gamma^{\ast}$ and given period-1 solution of $\bar{Z}_{2}$ (either $\bar{Z}_{2,r}^{(1)}$ in (42) or $\bar{Z}_{2,s}^{(1)}$ in (43)). Corresponding to $\bar{Z}_{2,r}^{(1)}$ in (42), one can obtain two solutions for $\bar{Z}_{0}$. The two sets of solutions for $(\bar{Z}_{0},\bar{Z}_{2})$ are verified by the numerical results of solving the whole nonlinear system of 14 equations with Newton’s method. No solution for $\bar{Z}_{0}$ is found corresponding to $\bar{Z}_{2,s}^{(1)}$ in (43). In addition, when $\gamma^{\ast}$ is extremely close to $\gamma_{m}^{+}$ (say $<10^{-8}$), more than two solutions of $\bar{Z}_{0}$ are found corresponding to $\bar{Z}_{2,r}^{(1)}$ in (42). In this case, the solutions are quite sensitive to $\gamma^{\ast}$ due to the term $(\gamma^{\ast}-\gamma_{m}^{+})\hat{n}_{3}(t)$ in $(\ref{3eq61})$. This is certainly expected, since for infinitely long cylinder in [46] we always get $\gamma^{\ast}=\gamma_{m}^{+}$, and thus any relatively large $\bar{Z}_{0}$ with $\gamma^{\ast}$ extremely close to $\gamma_{m}^{+}$ is always a solution with exponentially small error.

To sum up, for a given $\gamma^{\ast}$, the solution sets of $(\bar{Z}_{2},\bar{Z}_{0})$ can be determined by a system of two algebraic equations (39, 44), and each set can totally determine a solution for $\mathbf{U}_{0}$. For the period-1 solutions, we find that the width of PTR is always given by $\bar{Z}_{2,r}^{(1)}$ in (42). The analysis of $\bar{Z}_{2,r}^{(1)}$, including the dependence of material parameters, can be done easily. Many features are similar to those in the preceding section around (31). The main differences are that it slightly varies with $\gamma^{\ast}$ and is slightly larger than the width of PTR in two-region solutions, the $\bar{Z}_{0}^{(1)}$ in (28).

Since determination of $\bar{Z}_{0}$ from (39, 44) is very sensitive to $\gamma^{\ast}$, one could easily miss some solutions for a given $\gamma^{\ast}$, thereby missing part of stress-elongation curve. We find out that although one $\gamma^{\ast}$ corresponds to many solutions, one $\bar{Z}_{0}$ corresponds to only one period-1 solution. Therefore, from the perspective of numerical stability, it is a good strategy to control $\bar{Z}_{0}$ instead of $\gamma^{\ast}$ to find the solution set $(\gamma^{\ast},\bar{Z}_{2},\bar{Z}_{0})$ by (39, 44). In addition, the quantity $\bar{Z}_{0}$ has a one-to-one correspondence to the total elongation $\Delta$ of cylinder, making it easier to compare with curves in displacement-controlled experiments. This strategy is adopted in the computation of solutions of $\mathbf{U}_{0}$ and stress-elongation (stress-strain) curves. For two different $\bar{Z}_{0}$, Figure 9 shows the profiles of the axial strain $W_{0Z}$ in a region near $Z=0.5$. It shows that the middle smooth curve in PTR connects the low strain austenite phase and high-strain martensite phase. As $\bar{Z}_{0}$ increases (from Figure 9 to Figure 9), the PTR or transformation front gradually moves from middle to the ends.

From the solution of the axial strain $W_{0Z}$, one can calculate the total elongation of the cylinder by

[TABLE]

With material constants in (63) and $a=0.03$, Figure 10 shows the stress-elongation curve of the period-1 three-region solutions, where each blue dot is from one solution and the reference red line is from the previous two-region solutions. The stress always stays very close to the Maxwell stress $\gamma_{m}^{+}$, which explains the sensitivity to $\gamma^{\ast}$. This feature is in good agreement with the stress plateau in the curves from many experiments [8, 4, 9]. As $\Delta$ decreases to about 0.012, the width of MR gradually decreases to 0. The three-region solutions then degenerates to one solution with essentially two regions at a limit case. However, the stress-elongation curves of two-region and three-region solutions are not smoothly connected (see the gap indicated by the red circle in Figure 10), so the limit three-region solution is not included in the previous two-region solutions. We will examine the connection between the two-region and three-region solutions in the next section.

6 Connecting solutions

The mismatch of the period-1 two-region and three-region solutions is due to the interaction of interface $I_{PM}$ with the middle surface, and some intermediate solutions are missing. For the two-region solutions, the limit case is that some point on the middle surface $Z=0.5$ is on phase interface $\alpha=1$. More precisely, the center point of the middle surface first satisfies $\alpha=1$, when $\gamma^{\ast}=\gamma_{c}^{\ast}$ in (34), while other points do not. For the three-region solutions, the limit case is that the MR disappears, i.e. the whole middle surface is exactly the interface $\alpha=1$ (they coincide). In between, there are solutions that the interface condition $\alpha=1$ is satisfied from one point of the middle surface to the whole surface. In this case, both the interface conditions and symmetric conditions should be taken into account for determining such intermediate solutions. Now we adopt one interface condition and one symmetric condition

[TABLE]

which are satisfied by both limit solutions. More explicitly, the condition $(\ref{3eq63})_{1}$ is given by

[TABLE]

where the material constants $\bar{m}_{6},..,\bar{m}_{8}$ are given in Appendix B.

Following similar analysis as in Section 4, the expressions of $C_{1},..,C_{4}$ will be the same (see $(\ref{3eq11})$ and above), but $C_{5},..,C_{8}$ will be determined in terms of $\bar{Z}_{0}$ and $\gamma^{\ast}$ by $(\ref{eq87})_{2,3}$ and $(\ref{3eq63})$, whose explicit expressions are given in Appendix B. Finally the condition $(\ref{eq87})_{1}$ provides a single nonlinear algebraic equation for $\bar{Z}_{0}$

[TABLE]

where the lengthy expression is omitted. For this connecting case, the above equation makes sense only for $\gamma^{\ast}$ between two critical stresses determined by the two limit solutions in two-region and three-region cases. This small interval of $\gamma^{\ast}$ is below but not far from $\gamma^{+}_{m}$, and is independent of the radius $a$. For a given $\gamma^{\ast}$ in this small interval, there is a unique period-1 solution for $\bar{Z}_{0}$, between the two values for the width of PTR in the limit solutions. This in turn determines all other unknowns and the period-1 solution for $\mathbf{U}_{0}$. At the two critical stresses, one can recover the two limit solutions of the previous sections, thus such intermediate solutions naturally connect the two-region and three-region solutions. With $a=0.03$ and material constants in (63), Figure 10 shows the total elongation of the above connecting solutions indicated by the yellow dots, which smoothly connect the red curve and blue points for the two-region and three-region solutions.

7 The transition process

Based on results in the previous sections, we now investigate the phase transition process for the displacement-controlled loading process (unloading process is similar). The optimal solution is chosen by the energy criteria, if multiple solutions exist. We adopt the total pseudo-elastic energy as

[TABLE]

We focus on the homogeneous deformations and the inhomogeneous deformations from period-1 solutions in the previous sections, since one can easily show that period-k solution ($k\geq 2$) with larger width of PTR are not favored in terms of the above energy [46].

For the case of a relatively large radius $a>a_{c}=1/K_{2}$, where $K_{2}$ is given in (36)) and (37), the transition process is smooth and occurs gradually. Take $a=0.03$ with material constants in (63) for example, Figure 10 shows that there is a unique period-1 solution for each fixed total elongation $\Delta$ near nucleation of martensite phase. Therefore, the curve in displacement-controlled process goes exactly along the elongation curves in Figure 10, with a gradual smooth transition from homogeneous deformations to localized inhomogeneous deformations. At the critical stress or elongation $\gamma^{\ast}=\Delta=\gamma^{\ast}_{NM}$, the homogenous deformation in austenite phase is replaced by the two-region solutions with fixed PTR, and the stress drops gradually as $\Delta$ further increases. Subsequently, the two region solutions are replaced by connecting solutions and three-region solutions, as the PTR moves along the cylinder as $\Delta$ increases. The stress during propagation of PTR always stays close to Maxwell stress $\gamma_{m}^{+}$ as $\Delta$ increases. These features agree well with the experiments in [8, 9].

For the case of a relatively small radius $a<a_{c}$, the transition occurs abruptly, indicated by a jump in solutions in a snap-back bifurcation111One can also analogically think of the two cases with large $a$ and small $a$ as the supercritical and subcritical bifurcations used in dynamical systems.. Take $a=0.01$ with material constants in (63) for example, Figure 11 shows total elongation of homogeneous deformation in austenite phase and all the period-1 inhomogeneous deformations near the nucleation of martensite phase. It is clear that there exist three solutions for certain fixed elongation in $[\Delta_{1},\Delta_{2}]$, where

[TABLE]

and $K_{1}$ is given in $(\ref{3eq55})$. As the radius $a$ decreases, $\Delta_{1}$ will decrease and hence the interval $[\Delta_{1},\Delta_{2}]$ will become larger. The nucleation point lies within this interval, and will be determined by using the energy (LABEL:eq76). The energy difference is defined by

[TABLE]

where $\mathcal{W}^{+}(sol_{homo})$ is the energy with the homogeneous deformation (solution) in austenite phase which is set as a reference value, and $\mathcal{W}^{+}(sol_{inhomo})$ is the energy with any other inhomogeneous deformations. If $\Delta\mathcal{W}^{+}$ is positive, the homogeneous deformation is more stable. Otherwise if $\Delta\mathcal{W}^{+}$ is negative, the inhomogeneous deformations are more stable, and the one with minimal $\Delta\mathcal{W}^{+}$ is the optimal one. Figure 12 shows the energy differences in (52) scaled by a factor $\pi a^{2}$ for the period-1 inhomogeneous deformations. The red curve in Figure 12 depicts the energy difference for the two-region solutions corresponding to the red line in Figure 11), and the dots depict the energy difference for the connecting solutions and the three-region solutions corresponding to dots in Figure 11. By the minimal energy criterion, the solution will jump from homogeneous state in austenite phase to a localized inhomogeneous state (the limit two-region solution) at $\Delta=\Delta_{1}$, as shown in Figure 12. Then the stress goes up a bit with connecting solutions and soon goes along the stress plateau with three-region solutions, corresponding to the propagation of transformation front in experiments [8, 9].

The above quantitative analyses for large and small radius $a$ agree with the experimental observations in Figure 8 of [9]. With relatively large $a$, the transition process in experiments is smooth and stress-strain curve is in agreement with Figure 10. With relatively small radius of wire, the stress drop in experiments is very sharp, indicating an abrupt transition at one point, which is consistent with the sharp stress drop in Figure 12. This is due to the transition of solutions in a snap-back bifurcation in Figure 11. Figures 8 and 9 of [9] also show decreasing peak stress with increasing radius $a$, which does not agree with above analysis. One reason could be that the above analysis is the ideal case by minimal energy criterion, whereas in practical case the local stability, the imperfection of the material and the perturbations during experiments could influence the exact nucleation (peak) stress (see remark below). The other reason could be that the current work neglects some high-order nonlinear terms in hyperelastic energy, the importance of which will be left for future investigation.

**Remark: **From the above analysis, for a relatively small radius $a$ the stress will drop to the plateau at total elongation $\Delta_{1}$ before it reaches the peak of 1D intrinsic stress-strain curve at $\Delta_{2}$. This is the ideal case by minimal energy criterion, where the transition takes places at $\Delta_{1}$, by selecting the solution with global minimal energy. In practical experiments, the transition point would be influenced by the energy difference (barrier) of the solutions, the imperfection of the material and perturbations in the experimental setup and operations. As the homogeneous deformation is locally stable (local energy minimum) [45], the transition is very likely to occur after $\Delta_{1}$, until the perturbation is enough to overcome the energy barrier. For a relatively large $a$, in ideal case transition takes place at $\Delta_{2}$. But since the energy barrier is 0 at this transition point, perturbations or imperfections can lead to a transition well before total elongation $\Delta$ reaches $\Delta_{2}$. Therefore, experimental curves with small $a$ tend to climb further along intrinsic linear curve at austenite phase in Figure 3 and curves with large $a$ tend to concave down before reaching peak [9]. As a result, experiments show higher nucleation stress (larger $\Delta$) for smaller $a$ in [9].

The unloading process can be investigated similarly. We briefly mention the main modifications. The phase transition criterion changes to $D=-A^{-}$, and the equation in PTR changes accordingly. The 1D homogeneous unloading curve and Maxwell stress $\gamma_{m}^{-}$ modify as in Figure 3. The boundary surfaces are initially in martensite phase, i.e., the boundary condition $(\ref{3eq5})_{1}$ should be used. The two-region solutions will contain MR and PTR, and three-region solutions will have MR, PTR and AR from left to right in interval $[0,0.5]$. The typical features of stress-elongation curve will be similar, except that there is stress rise at nucleation of austenite phase instead of stress drop.

Finally, we present the stress-elongation curves for loading and unloading processes with $a=0.01$. Figure 13 shows the total elongations for homogeneous deformations in austenite and martensite phases and inhomogeneous deformations from period-1 solutions. Figure 13 shows the complete loading and unloading curves based on energy criterion, which are in good agreement with experiments in [8, 4, 9]. The four embedded figures are the enlarged parts of curves near nucleation and coalescence of austenite and martensite phases, where $\Delta_{1},\Delta_{3},\Delta_{4},\Delta_{5}$ denotes the total elongations at the four transition points.

8 Concluding remarks

Based on a 3D constitutive model, we have systematically investigated the phase transition process of SMA wire from homogeneous deformations to localized inhomogeneous deformations. The 3D system is simplified to a 1D piecewise linear system, while keeping essential high-dimensional effects. Free interfaces between different phase regions are included, and inhomogeneous deformations with two regions and three regions are studied analytically. The results clearly demonstrate important features in the transition process, such as the stress drop (rise) at nucleation of the martensite (austenite) phase and stress plateau at Maxwell stress during propagation of transformation front. The roles of material and geometrical constants are revealed. The width of PTR is proportional to radius and the dependence on material parameters is mainly through a parameter associated with strain softening behaviour. As radius decreases, the stress drop at nucleation of the martensite phase tends to be very sharp, due to transition of solutions in a snap-back bifurcation. The snap-back bifurcation sheds lights on the difficulties of direct numerical computations, since it is not an easy task to capture all the possible solutions with free interfaces in simulation and select the optimal one from these solutions. This study could provide insights into the mechanism and instability of other materials or structures with strain softening behaviour, and shed light on designing functionally graded materials to avoid instability.

The present work is restricted to the symmetric case, where phase transition takes place in the middle. When phase transition occurs at one end, there will be complicated interaction between phase interface and end surface (like the connecting solution in this paper), which will be left for future investigation. We also plan to extend the work to the investigation of geometrically and functionally graded SMAs.

Acknowledgement

The work was done at CityU and supported by a grant from the National Nature Science Foundation of China (Project No.: 11572272) and GRF grant (Project No.: CityU 11302417) from the Research Grants Council of Hong Kong SAR, China.

Appendix A Constants and expressions in section 2

The two moduli $\xi_{1}$ and $\xi_{2}$ are the first-order incremental elastic moduli of the strain energy $\Phi_{1}$ or $\Phi_{2}$ (assumed the same at first order) [46], which are related to the lamé constants [56] (or Young’s modulus $E$ and Poisson’s ratio $\nu$)

[TABLE]

The row vectors $\mathbf{b}_{k}^{(1)},\mathbf{b}^{(2)},\bar{\mathbf{b}}_{k}^{(1)},\bar{\mathbf{b}}^{(2)}$ ($k\geq 0$) are given by

[TABLE]

The matrix $\mathbf{H}$ is given by

[TABLE]

Similarly $\bar{\mathbf{H}}$ is defined in the same way by replacing $\mathbf{A},\mathbf{b}_{0}^{(1)},\mathbf{b}_{1}^{(1)},\mathbf{b}^{(2)}$ with $\bar{\mathbf{A}},\bar{\mathbf{b}}_{0}^{(1)},\bar{\mathbf{b}}_{1}^{(1)},\bar{\mathbf{b}}^{(2)}$, see also Appendix B of [46]. The vectors $\mathbf{f}_{0},\mathbf{f}_{1},\bar{\mathbf{f}}$ are given by

[TABLE]

where for loading process $\alpha_{00}$ is

[TABLE]

Appendix B Constants and expressions in solutions

The constants $d_{1},d_{2},q_{1},q_{2}$ are given by

[TABLE]

The quantities $d_{3},d_{4},\bar{q}_{1},\bar{q}_{2}$ can be expressed by the components of $\bar{\mathbf{H}}$, and with $\epsilon_{1}:=D_{k}^{\pm}/s_{2}^{2}$, we get

[TABLE]

The constants in the analytical results are given by

[TABLE]

The integrating constants $C_{9},..,C_{12}$ in three-region solutions are given by

[TABLE]

The integrating constants $C_{5},..,C_{8}$ in connecting solution are given by

[TABLE]

Appendix C Material constants and one numerical solution

The following material constants are used to show numerical results and figures

[TABLE]

One two-region solution from whole system by Newton’s method is given by

[TABLE]