.. _models-shape_memory_alloy:

************************
Shape Memory Alloy Model
************************

Theory
======

The shape memory alloy (SMA) model is used to describe the thermomechanical response of intermetallics (e.g. NiTi, NiTiCu, NiTiPd, NiTiPt) that can undergo a reversible, 
diffusionless, solid-to-solid martensitic transformation.  Specifically, the materials have a high-symmetry (typically cubic) austenitic crystallographic structure 
at high temperature and/or low stress.  At lower-temperatures and/or high stress the crystallographic structure is transformed to a lower symmetry (typically orthorhombic or 
monoclinic) martensitic phase.  The change in structure and symmetry may be taken advantage of to produce large inelastic strains of :math:`\approx` 1-8%.  Importantly, 
this class of materials differentiates itself from TRIP steels in that the transformation is *reversible* and a variety of thermomechanical loading paths 
have been conceived of to take advantage of this behavior.  A notable application of these materials is as an actuator in smart, morphing structures.

Phenomenologically, the macroscopic behavior of SMAs is typically discussed in effective stress-temperature space via a *phase diagram* 
like in :numref:`fig-sma-phaseDiagram`.  The four lines denoted :math:`M_s,~M_f,~A_s,` and :math:`A_f` indicate the martensitic start, martensitic finish, 
austenitic start, and austenitic finish transformation surfaces.  *Forward transformation* 
(from an austenitic to a martensitic state) is described by the martensitic start and finish surfaces.  
Specifically, the former refers to the thermomechanical conditions at which transformation 
will initiate while the latter corresponds to complete transformation.  The difference between the two surfaces is associated 
with internal hardening effects due to microstructure (i.e. texture, back stresses).  Transformation from martensite to austenite 
is referred to as *reverse* and is characterized by the austenitic start and finish surfaces.  Detailed 
discussion of the crystallography and phenomenology may be found in [:footcite:`mat:ref:lagBook`, :footcite:`mat:ref:otsukaRen`] (In the martensitic configuration, 
the crystallographic structure can either self-accommodate in a *twinned* configuration producing no macroscopic inelastic strain or 
an internal or external stress field may be used to *detwin* the microstructure thereby producing the desired inelastic strain.  
For simplicity, this distinction is bypassed in this brief text and the interested reader should consult the referenced works.).

.. _fig-sma-phaseDiagram:

.. figure:: ../../_static/figures/phaseDiagram.png
   :align: center
   :scale: 25 %

   Representative phase diagram of shape memory alloys highlighting characteristic loading paths (:math:`\left(A\right)` and :math:`\left(B\right)`), transformation surfaces, and phases.

Two responses characteristic of SMAs may also be represented via the phase diagram.  These are the *actuation response* and the *pseudoelastic* 
(often referred to as *superelastic* in the literature) responses.  The first (actuation) is indicated by path :math:`A` in :numref:`fig-sma-phaseDiagram`.  In this case, 
a mechanical bias load is applied to the SMA and the material is then thermally cycled through forward and reverse transformation.  The resulting transformation first produces 
and then removes the large transformation strains of SMAs and is commonly used for (surprisingly) actuation applications.  At higher temperatures (:math:`T>A_f`), 
mechanical loading may be used induce forward and, upon unloading, reverse transformation as indicated in path :math:`B` 
of :numref:`fig-sma-phaseDiagram`.  Through such a cycle, 
a distinctive flag shape in the stress strain response is observed through which large amounts of energy may be dissipated while producing no permanent 
deformations.  As such, this loading path is often considered for vibration isolation or damping applications.    

In LAME, the response of SMAs is described by the phenomenological model of Lagoudas and coworkers [:footcite:`mat:ref:lagoudas:2012`].  This model was motivated 
by actuator applications and it describes the inelastic deformation associated with martensitic transformation through two internal state variables -- 
the scalar martensitic volume fraction, :math:`\xi`, and tensorial transformation strain tensor, :math:`\varepsilon^{\text{tr}}_{ij}`.  Before proceeding it should be noted 
that the *structural* response of SMA specimens and components exhibit a rate dependency associated with the strong thermomechanical coupling 
of SMAs.  Specifically, the transformation process gives off/absorbs large amounts of energy via the *latent heat of transformation*.  The rate dependence 
observed is a result of the characteristic time scale associated with thermal transport of this heat.  In pure mechanical analyses (like Sierra/SM), this 
means quasistatics loadings are typically considered (a strain rate of :math:`\approx 1x10^{-4}` and/or heating/cooling rate of :math:`\approx 2^{\circ}C/\text{min}`).  Formulations 
accounting for the full coupling have been developed but require more complex implementations.

To begin, the model assumes an additive decomposition of the total, elastic, thermal, and transformation deformation (strain) rates respectively denoted 
by :math:`D_{ij},~D^{\text{e}}_{ij},~D^{\text{th}}_{ij}` and :math:`D^{\text{tr}}_{ij}` producing a total deformation rate of the form,

.. math::

   D_{ij}=D^{\text{e}}_{ij}+D^{\text{th}}_{ij}+D^{\text{tr}}_{ij}.
  


With respect to the thermoelastic deformations, it is noted that the different crystallographic phases have different thermoelastic constants.  
Previous studies have demonstrated that a rule of mixtures on the compliance and other material properties of the form,

.. math::

   \mathbb{S}_{ijkl}=\mathbb{S}_{ijkl}^A+\xi\left(\mathbb{S}_{ijkl}^M-\mathbb{S}_{ijkl}^A\right)=\mathbb{S}_{ijkl}^A+\xi\Delta \mathbb{S}_{ijkl}, \\
   \mathbb{\alpha}_{ij}=\alpha^A\delta_{ij}+\xi\left(\alpha^M\delta_{ij}-\alpha^A\delta_{ij}\right)=\alpha^A\delta_{ij}+\xi\Delta\alpha\delta_{ij},
 

 
in which :math:`\mathbb{S}_{ijkl}` and :math:`\alpha_{ij}` are the current effective compliance and coefficient of thermal expansion and the 
superscripts :math:`A` and :math:`M` denote thermoelastic properties in the austenitic and martensitic configuration.  The symbol :math:`\Delta` is used 
to indicate the difference in a property between the martensitic and austenitic phases while :math:`\delta_{ij}` is the Kronecker delta.  Isotropy is assumed 
for all these properties and the compliances are determined via the definition of elastic moduli and Poisson's ratio of the two phases -- 
:math:`E^A`, :math:`E^M`, :math:`\nu^M`, and :math:`\nu^M`.  The two Poisson ratios are often the same and take typical values for metals (:math:`\nu^A\approx\nu^M\approx 0.3`) 
while the elastic moduli can differ by a factor of more than two.  For instance the austenitic modulus of NiTi is typically given as :math:`\approx` 70 GPa while 
the martensitic one is :math:`\approx` 30 GPa (Given the lower symmetry of the martensitic phase the determination of an isotropic elastic modulus 
can vary with characterization methodology.  In this case, the apparent elastic modulus measured from macroscopic thermoelastic tests should be used.).  
Importantly, this difference means that the thermoelastic properties and corresponding deformations vary with transformation.  As such, 
the corresponding rates of deformation are given as,

.. math::

   D^{\text{e}}_{ij}=\dot{\xi}\Delta \mathbb{S}_{ijkl}\sigma_{kl}+\mathbb{S}_{ijkl}\stackrel{\circ}{\sigma}_{kl}, \\
   D^{\text{th}}_{ij}=\dot{\xi}\Delta\alpha\delta_{ij}\left(\theta-\theta_0\right)+\alpha\delta_{ij}\dot{\theta},

where :math:`\theta` and :math:`\theta_0` are the current and reference temperature and :math:`\sigma_{ij}` is the symmetric Cauchy stress.  
Note, in using the SMA model a temperature field *must* be defined.  The stress 
rate may then be shown to be,

.. math::

   \stackrel{\circ}{\sigma}_{ij}=\mathbb{C}_{ijkl}\left(D_{kl}-\alpha\delta_{kl}\dot{\theta}
                                -\dot{\xi}\left(\Delta\mathbb{S}_{klmn}\sigma_{mn}+\Delta\alpha\delta_{kl}\left(\theta-\theta_0\right)\right)
                                -D^{\text{tr}}_{kl}\right),



with :math:`\mathbb{C}_{ijkl}` being the current stiffness tensor defined as :math:`\mathbb{C}_{ijkl}=\mathbb{S}^{-1}_{ijkl}`.

To describe the transformation strain evolution, it is assumed that these deformations evolve with (and only with) the martensitic 
volume fraction, :math:`\xi`.  The corresponding flow rule is given as,

.. math::

   D_{ij}^{\text{tr}}=\dot{\xi}\Lambda_{ij},



and :math:`\Lambda_{ij}` is the transformation direction tensor assumed to be of the form,

.. math::
   :label: sma-eq:lambda

   \Lambda_{ij}=\left\{\begin{array}{cc}H^{\text{cur}}\left(\bar{\sigma}_{vM}\right)\frac{3}{2}\frac{s_{ij}}{\bar{\sigma}_{vM}} & \dot{\xi}\geq 0 \\
                                        \frac{\varepsilon^{\text{tr}-rev}_{ij}}{\xi^{rev}} & \dot{\xi} < 0 \end{array} \right. .



In :eq:`sma-eq:lambda`, :math:`H^{\text{cur}}` is the transformation strain magnitude that is dependent on the von Mises effective stress, :math:`\bar{\sigma}_{vM}`, 
and :math:`s_{ij}` is the deviatoric stress.  With forward transformation defined in this way, it is assumed that deformation is shear-based and follows 
a :math:`J_2` like flow direction.  For reverse transformation (:math:`\dot{\xi}<0`), the postulated form is utilized to ensure complete recovery of transformation 
strains with martensitic volume fraction.  In other words, all transformation strain components are zero-valued at :math:`\xi=0`.  Without enforcing this condition 
in this way, non-proportional loading paths could be constructed producing a non-zero transformation strain when the material is austenitic.  The transformation strain 
at load reversal, :math:`\varepsilon^{\text{tr}-rev}_{ij}`, and martensitic volume fraction at load reversal, :math:`\xi^{rev}`, are then tracked (via the implementation) 
and used for this purpose. 

The transformation strain magnitude, :math:`H^{\text{cur}}`, is a function of the von Mises effective stress (:math:`\bar{\sigma}_{vM}`) and is introduced to incorporate 
detwinning effects without introducing an additional internal state variable complicating the model.  Specifically, at low stress values, this function 
returns a minimum value.  If the microstructure is self-accommodated this value will be zero.  A decaying exponential is used such that as the 
stress increases the value of the strain magnitude becomes that of the maximum value incorporating both crystallographic and texture effects.  
The given functional form is,

.. math::

   H^{\text{cur}}=\left\{\begin{array}{cc} H_{\text{min}} & \bar{\sigma}_{vM} \leq \sigma_{\text{crit}} \\
                                           H_{\text{min}} + \left(H_{\text{sat}}-H_{\text{min}}\right)\left(1-\exp\left(-k\left(\bar{\sigma}_{vM}
                                           -\sigma_{\text{crit}}\right)\right)\right)
                                           & \bar{\sigma}_{vM} > \sigma_{\text{crit}} \end{array} \right. ,
  


where :math:`H_{\text{min}},~H_{\text{sat}},~k,` and :math:`\sigma_{\text{crit}}` are model parameters giving the minimum transformation strain magnitude, 
maximum transformation strain magnitude, exponential fitting parameter governing the transition zone, and critical stress values (in some ways 
analogous to the detwinning stress).

The evolution of martensitic transformation process is governed by a transformation function serving an analogous role to the yield function in plasticity.  
This function is given by,

.. math::

   f\left(\sigma_{ij},\theta,\xi\right)=\pm\phi\left(\sigma_{ij},\theta,\xi\right)-\bar{\sigma}\left(\sigma_{ij}\right),



with :math:`\phi` begin the thermodynamic driving force for transformation and :math:`\bar{\sigma}` the critical value.  The :math:`\pm` is used to denote either forward (:math:`+`) 
or reverse (:math:`-`) transformation.  This transformation function and the associated forms are derived from continuum thermodynamic considerations and the 
details of that process are neglected here for brevity but may be found in [:footcite:`mat:ref:lagoudas:2012`].  The functional forms of these variables are given as,

.. math::

   \phi\left(\sigma_{ij},\theta,\xi\right) & = \sigma_{ij}\Lambda_{ij}+\frac{1}{2}\sigma_{ij}\Delta\mathbb{S}_{ijkl}\sigma_{kl}+\sigma_{ij}\Delta\alpha\delta_{ij}\left(\theta-
                       \theta_0\right)+\rho\Delta s_0\theta-\rho\Delta u_0-f^t\left(\xi\right),\nonumber \\
   \bar{\sigma}\left(\sigma_{ij}\right) & = \sigma_0+D\sigma_{ij}\Lambda_{ij}, 


 
in which :math:`\rho\Delta s_0` and :math:`\rho\Delta u_0` are the differences in reference entropy and internal energy of the two phases, :math:`D` is a calibration 
parameter intended to capture variations in dissipation with stress, and :math:`f^t\left(\xi\right)` 
is the hardening function.  With respect to this latter term, empirical observations were used to arrive at a postulated form of,

.. math::

   f^t\left(\xi\right)=\left\{\begin{array}{cc}\frac{1}{2}a_1\left(1+\xi^{n_1}-\left(1-\xi\right)^{n_2}\right)+a_3 & \dot{\xi}\geq 0 \\
                                               \frac{1}{2}a_2\left(1+\xi^{n_3}-\left(1-\xi\right)^{n_4}\right)-a_3 & \dot{\xi}<0\end{array}\right. ,



with :math:`a_1,~a_2,` and :math:`a_3` being fitting parameters and :math:`n_1,~n_2,~n_3,` and :math:`n_4` are exponents fit to match the smooth transformation from elastic 
to inelastic deformations at the start of forward, end of forward, start of reverse, and end of reverse transformation respectively.

Before proceeding, one final note should be given in regards to calibration.  Specifically, some of the model parameters just listed 
(:math:`a_1,~a_2,~a_3,~D,~\sigma_0,~\rho\Delta s_0` and :math:`\rho\Delta u_0`) are not easily identified or conceptualized in terms of common thermomechanical 
experiments.  Some easily identifiable parameters (:math:`M_s,~M_f,~A_s`, and :math:`A_f`), however, are not evident in the theoretical formulation.  
Conditions associated with these terms and some physical constraints may be used to determine the model parameters in terms of these 
more accessible properties.  These relations are,

.. math::
   :label: eqn:sma-theory:rhoS

   \rho\Delta s_0=\frac{-2\left(C^MC^A\right)\left[H^{\text{cur}}\left(\sigma\right)+\sigma\frac{\partial H^{\text{cur}}}{\partial\sigma}
                        +\sigma\left(\frac{1}{E^M}-\frac{1}{E^A}\right)\right]}{C^M+C^A}|_{\sigma=\sigma^*}, \\
   D = \frac{\left(C^M-C^A\right)\left[H^{\text{cur}}\left(\sigma\right)+\sigma\frac{\partial H^{\text{cur}}}{\partial\sigma}
             +\sigma\left(\frac{1}{E^M}-\frac{1}{E^A}\right)\right]}{\left(C^M+C^A\right)\left[H^{\text{cur}}\left(\sigma\right) +
             \sigma\frac{\partial H^{\text{cur}}}{\partial\sigma}\right]}|_{\sigma=\sigma^*}, \\
   a_1=\rho\Delta s_0\left(M_f-M_s\right),~\qquad~a_2=\rho\Delta s_0\left(A_s-A_f\right), \\
   a_3=-\frac{a_1}{4}\left(1+\frac{1}{n_1+1}-\frac{1}{n_2+1}\right)+\frac{a_2}{4}\left(1+\frac{1}{n_3+1}-\frac{1}{n_4+1}\right), \\
   \rho\Delta u_0=\frac{\rho\Delta s_0}{2}\left(M_s+A_{f}\right),~\qquad~\sigma_0=\frac{\rho\Delta s_0}{2}\left(M_s-A_f\right)-a_3,
 


in which :math:`\sigma^*` is the scalar stress measure in which the calibration is performed at.  For additional discussion on the 
characterization of SMAs and calibration of this model, the user is referred to [:footcite:`mat:ref:hartl:2010:1`, :footcite:`mat:ref:hartl:2010:2`].

Implementation
==============

Similar to the various plasticity models in LAMÉ, an elastic predictor-inelastic corrector approach is used to 
perform the stress updating routine.  Unlike the other models, however, in the shape memory alloy routine a convex cutting 
plane (CCP) return mapping algorithm (RMA) is used in lieu of the closest point projection.  This difference essentially simplifies 
the integration of flow rule and the corresponding problem at the cost of some algorithmic stability.  
Prior studies [:footcite:`mat:ref:qidwai:00`] have shown that this implementation is sufficient for convergence in most cases while 
providing a substantial savings in cost.  The specific implementation used here is that of [:footcite:`mat:ref:lagoudas:2012`].

To compute an elastic trial state, a trial stress is determined assuming purely thermoelastic deformations such that,

.. math::

   T^{tr}_{ij} = \mathbb{C}_{ijkl}\left(\xi^n\right)\left(d_{kl}\Delta t - \alpha_{kl}\left(\xi^n\right)\Delta \theta\right), ~\qquad~ \Delta \theta = \theta^{n+1}-\theta^n.



In this case, it is assumed that the temperature fields are known at :math:`t_{n+1}` and :math:`t_n` (denoted :math:`\theta^{n+1}` and :math:`\theta^n`, respectively) and 
the thermoelastic properties are computed using the martensitic volume fraction at the previous time step :math:`\xi^n`.  At this stage, 
a perturbation stress (:math:`T_{ij}^{per}=T_{ij}^n+\beta\left(T^{tr}_{ij}-T^n_{ij}\right)` with :math:`\beta<<1`) is computed and used 
to determine local variations of the thermodynamic driving force, :math:`\phi`.  This is necessary to determine the direction of transformation 
(forward or reverse).  Using the full trial stress to this end can produce spurious results in some thermally-driven cases.  The trial 
yield function value is then computed as,

.. math::

   f^{tr}=f\left(T_{ij}^{tr},\theta^{n+1},\xi^n\right)=\pm\phi\left(T_{ij}^{tr},\theta^{n+1},\xi^n\right)-\bar{\sigma}\left(T_{ij}^{tr}\right).


If :math:`f^{tr}<0`, no nonlinear deformation occurs and the trial solution is accepted as the material state at :math:`t=t_{n+1}`.  When this 
condition is not satisfied, the CCP-RMA routine is used to correct the trial state and return it to the yield surface.

To perform the inelastic correction, the Newton-Raphson method is iteratively used to update the material state (:math:`T_{ij}` and :math:`\xi`) 
until convergence is achieved.  Denoting the current and next iteration by :math:`\left(k\right)` and :math:`\left(k+1\right)`, respectively, produces 
updating expressions of the form,

.. math::

   T_{ij}^{\left(k+1\right)} & = T_{ij}^{\left(k\right)}+\Delta T_{ij}, \nonumber \\
   \xi^{\left(k+1\right)} & = \xi^{\left(k\right)}+\Delta \xi,



with :math:`\xi^{\left(0\right)}` and :math:`T_{ij}^{\left(0\right)}` initialized to :math:`\xi^{n}` and :math:`T_{ij}^{tr}`, respectively.  The key difference 
between the CCP and closest point projection (CPP) methods is associated with how the inelastic strain flow rules are integrated.  
In the former method, an explicit evaluation of the flow direction is utilized while the latter is associated with a fully implicit 
expression.  For the CPP algorithms, this implicit expression means the flow rule must be solved in a nonlinear system of equations 
with the consistency equation.  Relaxing this assumption via the CCP method, however, produces an explicitly evaluated flow rule of,

.. math::

   \varepsilon^{t\left(k+1\right)}_{ij}=\varepsilon^{t\left(k\right)}_{ij}+\Delta \xi\Lambda_{ij}^{\left(k\right)}.
  


Importantly, this means that the only nonlinear equation to be solved is the scalar consistency equation (:math:`f=0`) which can be linearized 
such that, 

.. math::

   \Delta \xi = -\frac{f^{\left(k\right)}}{\frac{\partial f^{\left(k\right)}}{\partial\xi}-\frac{\partial f^{\left(k\right)}}{\partial T_{ij}}
              \mathbb{C}_{ijkl}\left(\xi^{\left(k\right)}\right)\left(\Delta\mathbb{S}_{klmn}T^{\left(k\right)}_{mn}+
              \Delta\alpha_{kl}\Delta\theta+\Lambda^{\left(k\right)}_{kl}\right)},



and the stress increment is then found as,

.. math::

   \Delta T^{\left(k\right)}=-\mathbb{C}_{ijkl}\left(\xi^{\left(n\right)}\right)\left(\Delta\mathbb{S}_{klmn}T^{\left(k\right)}_{mn} + 
                                                    \Delta\alpha_{kl}\Delta\theta+ \Lambda^{\left(k\right)}_{kl}\right)\Delta\xi.


Verification
============

The shape memory alloy model is verified through a series of thermomechanical loadings.  
The material properties and model parameters for these investigations are given in :numref:`tab-sma-verProps`.  
These properties correspond to those given in Table 3.4 in [:footcite:`mat:ref:lagBook`] with all :math:`n'`s assumed to be 1 
and setting :math:`E^M=E^A`.

.. _tab-sma-verProps:

.. csv-table:: The material and model parameters for the shape memory alloy model used during verification test.
   :align: center
   :delim: &

   :math:`E^A`      & 55 GPa & :math:`E^M`  & 55 GPa
   :math:`\nu^A`    & 0.33 & :math:`\nu^M` & 0.33
   :math:`\alpha^A` & 22.0\ :math:`x10^{-6}\frac{1}{\text{K}}` & :math:`\alpha^M` & 22.0\ :math:`x10^{-6}\frac{1}{\text{K}}`
   :math:`M_s`    & 245 K & :math:`A_s` & 270 K
   :math:`M_f` & 230 K  & :math:`A_f` & 280 K
   :math:`C^M` & 7.4 :math:`\frac{\text{MPa}}{\text{K}}` & :math:`C^A` & 7.4 :math:`\frac{\text{MPa}}{\text{K}}`
   :math:`H_{\text{min}}` & 0.056 & :math:`H_{\text{sat}}` & 0.056

It should also be clear that because :math:`H_{\text{min}}=H_{\text{sat}}` the model response is independent of the values of :math:`\sigma^{\text{crit}}` and 
:math:`k`.  For convenience, values of :math:`k=1.0x10^{6}` and :math:`\sigma^{\text{crit}}=0` will be used.  Additionally, :math:`\sigma^*` will be taken to be zero although 
inspection of :eq:`eqn:sma-theory:rhoS` and consideration of the relative magnitudes of the transformation strain and the difference 
in elastic strain similarly indicates an invariance in the model response to this parameter with constant :math:`H^{\text{cur}}`.  The default prestrain 
values are also utilized such that the SMA is initially austenitic.

Uniaxial Stress -- Pseudoelasticity
-----------------------------------

First, the isothermal (:math:`\theta>A_f`) pseudoelastic response through a uniaxial stress loading is explored.  
Importantly, the simplifications and model parameters described above (:math:`E^A=E^M=E`, :math:`H^{\text{cur}}\left(\bar{\sigma}_{vM}\right)=H`, 
:math:`C^A=C^M=C`, :math:`n_i=1`) allow for a simple analytical description of the pseudoelastic response (essentially trilinear).  For instance, 
given the constant slopes of the transformation surfaces, the stresses needed to induce or complete transformation are simply given 
by,

.. math::

   \sigma^{\beta}\left(\theta\right)=C\left(\theta-\beta\right),~\qquad \beta=M_s,~M_f,~A_s,~A_f,

where :math:`\sigma^{M_s}\left(\theta\right)` is the stress needed to start forward transformation at temperature, :math:`\theta`.  Given a stress value, the strain 
and material state may be completely determined by knowing the martensitic volume fraction, :math:`\xi`.  Specifically, the axial (here taken to be the :math:`1` 
direction) strain is simply,

.. math::

   \varepsilon_{11}=\frac{\sigma}{E}+\xi H,
  
and the lateral strains are 

.. math::

   \varepsilon_{22}=\varepsilon_{33}=-\nu\frac{\sigma}{E}-\frac{1}{2}\xi H,

in which the fact that the transformation strain tensor is deviatoric is being leveraged.  
The martensitic volume fraction may then simply be found by noting that :math:`f=0` 
during transformation.  Therefore, for forward transformation,

.. math::
   :label: eq:sma-ver-xi

   \xi=\left\{\begin{array}{cc} 0 & \sigma \leq \sigma^{M_s} \\
                           \frac{1}{a_1}\left(\sigma H+\rho\Delta s_0\theta-\rho\Delta u_0 - a_3-\sigma_0\right) & \sigma^{M_s} < \sigma < \sigma^{M_f} \\
                            1 & \sigma \geq \sigma^{M_f} \end{array} \right . .

A comparable expression is easily determined for reverse transformation.

 
The results of this simple analytical expression and those determined by Adagio are presented in :numref:`fig-ver-pseudo`  
for three different temperatures.  :numref:`fig-ver-pseudo`\ (a) presents the stress-strain response under these conditions while 
:numref:`fig-ver-pseudo`\ (b) presents the evolution of the martensitic volume fraction. 

.. _fig-ver-pseudo:

.. subfigure:: AB
   :subcaptions: below
   :align: center

   .. image:: ../../_static/figures/pseudoStressStrain.png
      :alt: Axial stress-strain
      :scale: 25 %

   .. image:: ../../_static/figures/pseudoMVF.png
      :alt: MVF
      :scale: 25 %

   Axial stress-strain response (a) and martensitic volume fraction (b), :math:`\xi`, 
   evolution determined analytically and via adagio for three different ambient temperatures :math:`\theta=300,~320` and :math:`340` K.

Constant Stress Actuation
-------------------------

To consider thermally driven transformation, the constant stress actuation response is investigated.  In such a loading, a mechanical load 
is applied at high temperature (:math:`\theta>A_f`) and held constant while the specimen is cooled through forward transformation and then heated 
back to its initial state.  Given the aforementioned simplifications to the model parameters, the analytical response is determined 
in a very similar fashion to that of pseudoelasticity.  In this instance, critical temperatures needed for transformation are first determined by 

.. math::

   \beta^{\sigma}=\beta + \frac{\sigma}{C},~\qquad \beta=M_s,~M_f,~A_s,~A_f,

with :math:`M_s^{\sigma}` being the temperature needed to start forward transformation at an effective stress, :math:`\sigma`.  The zero-stress value is 
:math:`M_s`.  Similarly, the axial and lateral strains may be adjusted as,

.. math::

   \varepsilon_{11} & = \frac{\sigma}{E}+\xi H+\alpha\left(\theta-\theta_0\right), \\
   \varepsilon_{22} = \varepsilon_{33} & = -\nu\frac{\sigma}{E}-\frac{1}{2}\xi H + \alpha\left(\theta-\theta_0\right).

The martensitic volume fraction is found through relations :eq:`eq:sma-ver-xi` albeit with the piecewise intervals defined in 
terms of temperature (*e.g* :math:`\sigma^{Ms}<\sigma<\sigma^{M_f}~\leftrightarrow~M_f^{\sigma}<\theta<M_s^{\sigma}`).  Results 
for the axial strain-temperature, lateral strain-temperature, and martensitic volume fraction-temperature as determined 
analytically and via adagio are presented below in :numref:`fig-ver-actuation`\ (a), :numref:`fig-ver-actuation`\ (b), 
and :numref:`fig-ver-actuation`\ (c), respectively.

.. _fig-ver-actuation:

.. subfigure:: ABC
   :subcaptions: below
   :align: center

   .. image:: ../../_static/figures/actuationStrain.png
      :alt: Axial stress-strain
      :scale: 25 %

   .. image:: ../../_static/figures/actuationLatStrain.png
      :alt: Laterial Strain
      :scale: 25 %

   .. image:: ../../_static/figures/actuationMVF.png
      :alt: MVF
      :scale: 25 %

   Axial stress-strain (a), lateral strain (b), and martensitic volume fraction (c), :math:`\xi`, evolution as a function of temperature as determined analytically and numerically (Sierra/SM).  Results are presented for three different applied bias stresses :math:`\sigma=100,~200` and :math:`300` MPa.

User Guide
==========

.. code-block:: sierrainput

   BEGIN PARAMETERS FOR MODEL SHAPE_MEMORY_ALLOY
     #
     # Elastic constants
     #
     YOUNGS MODULUS = <real>
     POISSONS RATIO = <real>
     SHEAR MODULUS  = <real>
     BULK MODULUS   = <real>
     LAMBDA         = <real>
     TWO MU         = <real>
     #
     # Thermoelastic properties of two crystallographic phases
     #
     ELASTIC MODULUS AUSTENITE  = <real>
     POISSON RATIO AUSTENITE    = <real>
     CTE AUSTENITE              = <real>
     ELASTIC MODULUS MARTENSITE = <real>
     POISSON RATIO MARTENSITE   = <real>
     CTE MARTENSITE             = <real>
     #
     # Phase diagram parameters
     #
     MARTENSITE START                  = <real>
     MARTENSITE FINISH                 = <real>
     AUSTENITE START                   = <real>
     AUSTENITE FINISH                  = <real>
     STRESS INFLUENCE COEFF MARTENSITE = <real>
     STRESS INFLUENCE COEFF AUSTENITE  = <real>
     #
     # Transformation strain magnitude parameters
     #
     H_MIN          = <real>
     H_SAT          = <real>
     KT             = <real>
     SIGMA_CRITICAL = <real>
     #
     # Calibration parameters
     #
     N1         = <real>
     N2         = <real>
     N3         = <real>
     N4         = <real>
     SIGMA STAR = <real>
     T0         = <real>
     #
     # Initial phase conditions
     #
     XI0                 = <real> (0.0)
     PRESTRAIN_DIRECTION = <int> (0)
     PRESTRAIN_MAGNITUDE = <real> (0.0)
     #
   END [PARAMETERS FOR MODEL SHAPE_MEMORY_ALLOY]


In the command blocks that define the Shape Memory Alloy model:


-  Although the thermoelastic constants of the phases are defined separately, the definition of these constants in this form is necessary for the global solver.  Typical values of the phases should be applied.
- The isotropic elastic moduli of the austenitic (:math:`E^A`) and martensitic phases (:math:`E^M`) are defined with the ``ELASTIC MODULUS AUSTENITE`` and ``ELASTIC MODULUS MARTENSITE`` command lines, respectively.  Note, alternative elastic constants (e.g. bulk or shear moduli) may not be used.
- The isotropic Poisson's ratio of the austenitic (:math:`\nu^A`) and martensitic phases (:math:`\nu^M`) are defined with the ``POISSON RATIO AUSTENITE`` and ``POISSON RATIO MARTENSITE`` command lines, respectively.  Note, alternative elastic constants (e.g. lame constant) may not be used.
- The isotropic coefficient of thermal expansion of the austenitic (:math:`\alpha^A`) and martensitic phases (:math:`\alpha^M`) are defined with the ``CTE AUSTENITE`` and ``CTE MARTENSITE`` command lines, respectively.  Note, given the phase and history dependence of the material thermal expansion, the use of artificial or thermal strain functions may not lead to desired results.  The use of these constants in encouraged instead.
- The zero stress, smooth transformation temperatures corresponding to the start and end of forward transformation (martensitic start :math:`M_s` and finish :math:`M_f`, respectively) and start and end of reverse transformation (austenitic start :math:`A_s` and finish :math:`A_f`, respectively) are given by the (in order) command lines ``MARTENSITE START``, ``MARTENSITE FINISH``, ``AUSTENITE START``, and ``AUSTENITE FINISH``. 
- The stress influence coefficients giving the slope of the forward and reverse transformation surfaces (:math:`C^M` and :math:`C^A`, respectively) are given by the ``STRESS INFLUENCE COEFF MARTENSITE`` and ``STRESS INFLUENCE COEFF AUSTENITE``, respectively.
- The stress dependence of the transformation strain magnitude requires four coefficients.  These are the minimum transformation strain magnitude (:math:`H_{\text{min}}`), the saturation (or maximum) magnitude (:math:`H_{\text{sat}}`), exponential fitting coefficient (:math:`k`), and critical effective stress value below which the magnitude is minimum (:math:`\sigma_{\text{crit}}`).  These parameters are defined via the ``H_MIN``, ``H_SAT``, ``KT``, and ``SIGMA_CRITICAL`` command lines, respectively.
- The smooth hardening fitting constants :math:`n_1,~n_2,~n_3,` and :math:`n_4` correspond to the degree of smoothness (essentially how gradual the transformation is) of the martensitic start, martensitic finish, austenitic start, and austenitic finish transformation surfaces.  They are given by the ``N1``, ``N2``, ``N3``, and ``N4`` command lines, respectively, and should take values :math:`0<n_i\leq 1`.
- The stress level of transformation at which calibration is performed is denoted by :math:`\sigma^*` and given by the command line ``SIGMA STAR``.  For thermally induced transformation this corresponds to the bias stress level while in pseudo-elastic loadings it corresponds to the stress level at which the material is roughly evenly split between martensite and austenite.
- The zero-strain reference temperature is denoted :math:`\theta_0` and prescribed via the ``T0`` command line.
- Three optional parameters describing the initial state of the material may be input.  These parameters are intended for the case in which the material is initially martensite to allow for initial heating and transformation recovery.  The first is the initial martensitic volume fraction, :math:`\xi\left(t=0\right)`, input via the ``XI0`` command line.  If this parameter is not specified the default value is :math:`0.0` representative of an austenitic material.  A value between :math:`0.0` and :math:`1.0` may be entered to initialize the material to partially (or fully) martensitic.  Corresponding initial transformation strains may be entered via the ``PRESTRAIN_DIRECTION`` :math:`\left(n^{\text{ps}}\right)` and ``PRESTRAIN_MAGNITUDE`` :math:`\left(||\varepsilon^{\text{tr}}_{ij}\left(t=0\right)||\right)` commands.  The first (an integer between one and three) gives the direction of transformation (in global Cartesian space) and the magnitude of the inelastic strain in that direction is given by a fraction (between 0 and 1) of :math:`H_{\text{sat}}` via the second ``PRESTRAIN_MAGNITUDE`` line.  As the transformation strain tensor is deviatoric, the other two directions are specified by preserving that the tensor be trace less.  Note, the ``PRESTRAIN_DIRECTION`` and ``PRESTRAIN_MAGNITUDE`` cannot be specified without a non-zero ``XI0`` definition.
- Thermal strains functions and commands in Sierra should not be used in conjunction with the ``shape_memory_alloy`` model.  
- In order for this material model to function a temperature must be defined on every block that uses the model.

Output variables available for this model are listed in :numref:`out-tab-smastvar`.

.. _out-tab-smastvar:

.. csv-table:: State Variables for SHAPE MEMORY ALLOY Model
   :align: center
   :delim: &
   :header: Name, Description

   ``MVF``  & martensitic volume fraction, :math:`\xi` 
   ``TransStrain`` & transformation strain tensor, :math:`\varepsilon_{ij}^{\text{tr}}` 

.. raw::
   html

   <hr>

.. footbibliography::
