# Constitutive laws

Here the phenomenological constitutive laws that are included in SMART+ are presented. These laws are implemented in the form of ABAQUS user materials (UMAT routines).
In the constitutive models the small strain theory assumptions are considered and thermal strains are accounted. The total strain increment $$\Delta\varepsilon^{\textrm{tot}}_{ij}$$ is additively decomposed into an elastic part $$\Delta\varepsilon^{\textrm{el}}_{ij}$$, a thermal part $$\Delta\varepsilon^{\textrm{th}}_{ij}$$ and (if it exists) an inelastic part $$\Delta\varepsilon^{\textrm{in}}_{ij}$$,
$\Delta\varepsilon^{\textrm{tot}}_{ij}=\Delta\varepsilon^{\textrm{el}}_{ij}+\Delta\varepsilon^{\textrm{th}}_{ij}+\Delta\varepsilon^{\textrm{in}}_{ij}.$
The inelastic strain can be plastic strain for elastoplastic materials, transformation strain for shape memory alloys etc. The elastic strain increment is connected with the stress increment $$\Delta\sigma_{ij}$$ through the elastic stiffness tensor $$L_{ijkl}$$
$\Delta\sigma_{ij}=L_{ijkl}~\Delta\varepsilon^{\textrm{el}}_{kl},$
while the thermal strain increment is related with the increment of temperature $$\Delta T$$ through the thermal expansion coefficients tensor $$\alpha_{ij}$$,
$\Delta\varepsilon^{\textrm{th}}_{ij}=\alpha_{ij}~\Delta T.$

The numerical schemes of the models are based on the return mapping algorithms. In these algorithms, the user provides i) the initial values of all the variables, like the total strain $$\varepsilon^{\textrm{tot(init)}}_{ij}$$, stress $$\sigma^{\textrm{init}}_{ij}$$, temperature $$T^{\textrm{init}}_{ij}$$ etc, and ii) the increments of the total strain $$\Delta\varepsilon^{\textrm{tot}}_{ij}$$ and the temperature $$\Delta T^{\textrm{init}}_{ij}$$ . The program computes a) the total final stress $$\sigma^{\textrm{fin}}_{ij}$$ and b) the tangent stiffness tensor $$L^t_{ijkl}={\partial\sigma_{ij}}/{\partial\varepsilon^{\textrm{tot}}_{kl}}$$. For all the constitutive models three options are considered:

1. 1D,
2. plane stress,
3. generalized plane strain/3D analysis.

## 1. Elasticity

Two types of elastic materials are utilized in SMART+:

• Elastic isotropic materials,
• elastic transversely isotropic materials.

### 1.1 Isotropic materials

In thermoelastic isotropic materials three parameters are required:

1. The Young modulus $$E$$,
2. the Poisson ratio $$\nu$$,
3. the coefficient of thermal expansion $$\alpha$$.

The elastic stiffness tensor and the thermal expansion coefficients tensor are written in SMART+ formalism as
$\boldsymbol{L}=\left(\begin{matrix} L_{1111} & L_{1122} & L_{1122} & 0 & 0 & 0 \\ L_{1122} & L_{1111} & L_{1122} & 0 & 0 & 0 \\ L_{1122} & L_{1122} & L_{1111} & 0 & 0 & 0 \\ 0 & 0 & 0 & L_{1212} & 0 & 0 \\ 0 & 0 & 0 & 0 & L_{1212} & 0 \\ 0 & 0 & 0 & 0 & 0 & L_{1212} \end{matrix}\right),~~~ \boldsymbol{\alpha}=\left(\begin{matrix} \alpha & 0 & 0 \\ 0 & \alpha & 0 \\ 0 & 0 & \alpha \end{matrix}\right),$
with
$L_{1111}=\frac{E(1-\nu)}{(1+\nu)(1-2\nu)},~~~L_{1122}=\frac{E\nu}{(1+\nu)(1-2\nu)},~~~L_{1212}=\frac{E}{2(1+\nu)}.$
Details on the the elastic stiffness tensor of isotropic media can be found in Lai et al 2010. The tangent stiffness tensor in this case is $$\boldsymbol{L}^t=\boldsymbol{L}$$. Moreover, the increment of the elastic strain is given by
$\Delta\varepsilon^{\textrm{el}}_{ij}=\Delta\varepsilon^{\textrm{tot}}_{ij}-\alpha\Delta T\delta_{ij},$
where $$\delta_{ij}$$ implies the Kronecker delta.

• In the 1D case only one component of stress is computed, through the relation
$\sigma^{\textrm{fin}}_{11}=\sigma^{\textrm{init}}_{11}+E\Delta\varepsilon^{\textrm{el}}_{11}.$
• In the plane stress case only three components of stress are computed, through the relations
$\left(\begin{matrix} \sigma^{\textrm{fin}}_{11} \\ \sigma^{\textrm{fin}}_{22} \\ \sigma^{\textrm{fin}}_{12} \end{matrix}\right) =\left(\begin{matrix} \sigma^{\textrm{init}}_{11} \\ \sigma^{\textrm{init}}_{22} \\ \sigma^{\textrm{init}}_{12} \end{matrix}\right)+\frac{E}{1-\nu^2} \left(\begin{matrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{matrix}\right) \left(\begin{matrix} \Delta\varepsilon^{\textrm{el}}_{11} \\ \Delta\varepsilon^{\textrm{el}}_{22} \\ 2\Delta\varepsilon^{\textrm{el}}_{12} \end{matrix}\right).$
• In the generalized plane strain/3D analysis case the stress tensor is computed through the relation
$\sigma^{\textrm{fin}}_{ij}=\sigma^{\textrm{init}}_{ij}+L_{ijkl}~\Delta\varepsilon^{\textrm{el}}_{kl}.$

### 1.2 Transversely isotropic materials

In thermoelastic transversely isotropic materials eight parameters are required:

1. The axis of transverse isotropy (1,2 or 3),
2. the axial Young modulus $$E_L$$,
3. the transverse Young modulus $$E_T$$,
4. the axial Poisson ratio $$\nu_{TL}$$,
5. the transverse Poisson ratio $$\nu_{TT}$$,
6. the axial shear modulus $$G_{LT}$$,
7. the axial coefficient of thermal expansion $$\alpha_L$$,
8. the transverse coefficient of thermal expansion $$\alpha_T$$.

When the axis of transverse isotropy is 1, the elastic stiffness tensor and the thermal expansion coefficients tensor are written in SMART+ formalism as
$\boldsymbol{L}=\left(\begin{matrix} L_{1111} & L_{1122} & L_{1122} & 0 & 0 & 0 \\ L_{1122} & L_{2222} & L_{2233} & 0 & 0 & 0 \\ L_{1122} & L_{2233} & L_{2222} & 0 & 0 & 0 \\ 0 & 0 & 0 & L_{1212} & 0 & 0 \\ 0 & 0 & 0 & 0 & L_{1212} & 0 \\ 0 & 0 & 0 & 0 & 0 & L_{2323} \end{matrix}\right),~~~ \boldsymbol{\alpha}=\left(\begin{matrix} \alpha_L & 0 & 0 \\ 0 & \alpha_T & 0 \\ 0 & 0 & \alpha_T \end{matrix}\right),$
where
$\begin{array}{c}\displaystyle{L_{1111}=\frac{E_L}{\omega}~(\nu^2_{TT}-1), ~~~L_{1122}=-\frac{E_L}{\omega}~\nu_{TL}~(1+\nu_{TT}), ~~~L_{2222}=\frac{E_L~\nu^2_{TL}-E_T}{\omega},} \\ \displaystyle{L_{2233}=-\frac{E_L~\nu^2_{TL}+E_T~\nu_{TT}}{\omega},~~~L_{1212}=G_{LT},~~~L_{2323}=\frac{E_T}{2(1+\nu_{TT})},}\\ \displaystyle{\omega=\frac{1}{E_T}(1+\nu_{TT})~(2E_L~\nu^2_{TL}+E_T~(\nu_{TT}-1)).}\end{array}$
Details on the elastic stiffness tensor of transversely isotropic media can be found in Christensen 1979. For axis of transverse isotropy being 2 or 3, the above tensors are properly rotated. The tangent stiffness tensor in this case is $$\boldsymbol{L}^t=\boldsymbol{L}$$. Moreover, the increment of the elastic strain is given by
$\Delta\varepsilon^{\textrm{el}}_{ij}=\Delta\varepsilon^{\textrm{tot}}_{ij}-\alpha_{ij}\Delta T.$

• In the 1D case only one component of stress is computed, through the relation
$\sigma^{\textrm{fin}}_{11}=\sigma^{\textrm{init}}_{11}+L_{1111}\Delta\varepsilon^{\textrm{el}}_{11}.$
• In the plane stress case only three components of stress are computed, through the relations
$\left(\begin{matrix} \sigma^{\textrm{fin}}_{11} \\ \sigma^{\textrm{fin}}_{22} \\ \sigma^{\textrm{fin}}_{12} \end{matrix}\right) =\left(\begin{matrix} \sigma^{\textrm{init}}_{11} \\ \sigma^{\textrm{init}}_{22} \\ \sigma^{\textrm{init}}_{12} \end{matrix}\right)+\boldsymbol{K} \left(\begin{matrix} \Delta\varepsilon^{\textrm{el}}_{11} \\ \Delta\varepsilon^{\textrm{el}}_{22} \\ 2\Delta\varepsilon^{\textrm{el}}_{12} \end{matrix}\right),$
with
$\boldsymbol{K}=\left(\begin{matrix} \displaystyle{L_{1111}-\frac{L_{1133}L_{3311}}{L_{3333}}} & \displaystyle{L_{1122}-\frac{L_{1133}L_{3322}}{L_{3333}}} & \displaystyle{L_{1112}-\frac{L_{1133}L_{3312}}{L_{3333}}} \\ \displaystyle{L_{2211}-\frac{L_{2233}L_{3311}}{L_{3333}}} & \displaystyle{L_{2222}-\frac{L_{2233}L_{3322}}{L_{3333}}} & \displaystyle{L_{2212}-\frac{L_{2233}L_{3312}}{L_{3333}}} \\ \displaystyle{L_{1211}-\frac{L_{1233}L_{3311}}{L_{3333}}} & \displaystyle{L_{1222}-\frac{L_{1233}L_{3322}}{L_{3333}}} & \displaystyle{L_{1212}-\frac{L_{1233}L_{3312}}{L_{3333}}} \end{matrix}\right).$
• In the generalized plane strain/3D analysis case the stress tensor is computed through the relation
$\sigma^{\textrm{fin}}_{ij}=\sigma^{\textrm{init}}_{ij}+L_{ijkl}~\Delta\varepsilon^{\textrm{el}}_{kl}.$

## 2. Elastoplasticity

The elastoplastic constitutive law implemented in SMART+ is a rate independent, isotropic, von Mises type material with exponential isotropic hardening. Six parameters are required:

1. The Young modulus $$E$$,
2. the Poisson ratio $$\nu$$,
3. the coefficient of thermal expansion $$\alpha$$,
4. the von Mises equivalent yield stress limit $$\sigma_{Y}$$,
5. the hardening parameter $$k$$,
6. the hardening exponent $$m$$.

The constitutive law is given by the set of equations
$\begin{matrix} {\sigma}_{ij}=L_{ijkl}\left({\varepsilon}^{\textrm{tot}}_{kl}-\alpha_{kl}\left(T-T^{\textrm{ref}}\right)-{\varepsilon}^{\textrm{p}}_{kl}\right),\\ \dot{\varepsilon}^{\textrm{p}}_{ij}=\dot{p}\Lambda_{ij},~~~\Lambda_{ij}=\frac{3}{2}\frac{\sigma’_{ij}}{\overline{\sigma}}, ~~~\sigma’_{ij}=\sigma_{ij}-\frac{1}{3}\sigma_{kk}\delta_{ij},~~~\overline{\sigma}=\sqrt{\frac{3}{2}\sigma’_{kl}\sigma’_{kl}},\\ \Phi=\overline{\sigma}-\sigma_{Y}-kp^m\leq 0,~~~ \dot{p}\geq0,~~~ \dot{p}~\Phi=0, \end{matrix}$
where $${\varepsilon}^{\textrm{p}}_{ij}$$ is the plastic strain tensor, $$p$$ is the plastic multiplier, $$\sigma’_{ij}$$ is the deviatoric part of the stress and $$\overline{\sigma}$$ is the von Mises equivalent stress (Lemaitre and Chaboche, 2002). Moreover, $$T^{\textrm{ref}}$$ is a reference temperature (usually the temperature at the beginning of the analysis). Finally, $$\delta_{ij}$$ denotes the Kronecker delta.

In SMART+ the elastoplastic material constitutive law is implemented using a return mapping algorithm. Two options are provided: the convex cutting plane and the closest point projection approach. The updated stress is provided for 1D, plane stress, and generalized plane strain/3D analysis according to the forms of elastic isotropic materials.

## 3. Viscoelasticity

The viscoelastic constitutive law implemented in SMART+ is based on the Zener model (spring in series with a Kelvin model). Seven parameters are required:

1. The bulk modulus, $$K_1$$, and the shear modulus, $$G_1$$, of a 3D spring attached in series with a Kelvin type spring-damper system,
2. the bulk modulus, $$K_2$$, and the shear modulus, $$G_2$$, of a 3D spring in parallel with a viscous damper,
3. the shear damping parameter, $$\eta_s$$, and the bulk damping parameter, $$\eta_b$$, of a viscous damper,
4. the coefficient of thermal expansion $$\alpha$$.

The analytical form of this constitutive law is
$\begin{array}{r}\displaystyle{\dot{\varepsilon}^{\textrm{tot}}_{ij}+\frac{G_2}{\eta_s}{\varepsilon}^{\textrm{tot}}_{ij}+C{\varepsilon}^{\textrm{tot}}_{kk}\delta_{ij}-\alpha\dot{T}\delta_{ij}-\alpha\left(\frac{G_2}{\eta_s}+3C\right)(T-T^{\textrm{ref}})\delta_{ij}} =\displaystyle{\frac{1}{2G_1}\dot{\sigma}_{ij}+A\dot{\sigma}_{kk}\delta_{ij}}\\ \displaystyle{+E\sigma_{ij}+D{\sigma}_{kk}\delta_{ij}},\end{array}$
where
$\begin{array}{rl} A=&\displaystyle{\frac{1}{3}\left(\frac{1}{3K_1}-\frac{1}{2G_1}\right)}, \\ B=&\displaystyle{\frac{1}{3}\left(\frac{1}{3\eta_b}-\frac{1}{2\eta_s}\right)},\\ C=&\displaystyle{3BK_2+\frac{1}{2\eta_s}\left(K_2-\frac{2G_2}{3}\right)},\\ D=&\displaystyle{B+C\frac{1}{3K_1}+\frac{G_2}{\eta_s}A},\\ E=&\displaystyle{\frac{1}{2\eta^s}\left(1+\frac{G_2}{G_1}\right)}, \end{array}$
$$\delta_{ij}$$ denotes the Kronecker delta, and $$T^{\textrm{ref}}$$ is a reference temperature (usually equal to the temperature at the beginning of the analysis). Introducing the backward Euler scheme
$\begin{array}{c}\varepsilon^{\textrm{tot(fin)}}_{ij}=\varepsilon^{\textrm{tot(init)}}_{ij}+\Delta\varepsilon^{\textrm{tot}}_{ij},~~~ \Delta\varepsilon^{\textrm{tot}}_{ij}=\dot{\varepsilon}^{\textrm{tot}}_{ij}\Delta t,~~~ \sigma^{\textrm{fin}}_{ij}=\sigma^{\textrm{init}}_{ij}+\Delta\sigma_{ij},~~~ \Delta\sigma_{ij}=\dot{\sigma}_{ij}\Delta t,\\ T^{\textrm{fin}}=T^{\textrm{init}}+\Delta T,~~~ \Delta T=\dot{T}\Delta t,\end{array}$
and after proper calculations, the numerical form on the Zener model is written
$\Delta\sigma_{ij}=L^{\Delta\varepsilon}_{ijkl}\Delta\varepsilon^{\textrm{tot}}_{kl}+ L^{\varepsilon}_{ijkl}\varepsilon^{\textrm{tot(init)}}_{kl}+L^{\sigma}_{ijkl}\sigma^{\textrm{init}}_{kl} +L^{\Delta T}_{ij}\Delta T+L^{T}_{ij}(T^{\textrm{init}}-T^{\textrm{ref}}),$
where
$\begin{array}{rl} F=&\displaystyle{\frac{1}{2G_1}+E\Delta t+3A+3D\Delta t}, \\ H=&\displaystyle{\frac{1}{2G_1}+E\Delta t}, \end{array}$
$\begin{array}{l} L^{\Delta\varepsilon}_{ijkl}=\displaystyle{\frac{1}{H}\left[C\Delta t-\frac{A+D\Delta t}{F}\left(1+\frac{G_2}{\eta_s}\Delta t+3C\Delta t\right)\right]3\mathcal{I}^{\textrm{vol}}_{ijkl}+\frac{1}{H}\left(1+\frac{G_2}{\eta_s}\Delta t\right)\mathcal{I}_{ijkl}},\\ L^{\varepsilon}_{ijkl}=\displaystyle{\frac{1}{H}\left[C\Delta t-\frac{A+D\Delta t}{F}\left(\frac{G_2}{\eta_s}\Delta t+3C\Delta t\right)\right]3\mathcal{I}^{\textrm{vol}}_{ijkl}+\frac{1}{H}\frac{G_2}{\eta_s}\Delta t~\mathcal{I}_{ijkl}},\\ L^{\sigma}_{ijkl}=\displaystyle{-\frac{1}{H}\left[D\Delta t-\frac{A+D\Delta t}{F}\left(E\Delta t+3D\Delta t\right)\right]3\mathcal{I}^{\textrm{vol}}_{ijkl} -\frac{1}{H}E\Delta t~\mathcal{I}_{ijkl}},\\ L^{\Delta T}_{ij}=\displaystyle{-\frac{1}{H}\alpha\left(1+\frac{G_2}{\eta_s}\Delta t+3C\Delta t\right)\left(1-3\frac{A+D\Delta t}{F}\right)\delta_{ij}},\\ L^{T}_{ij}=\displaystyle{-\frac{1}{H}\alpha\left(\frac{G_2}{\eta_s}\Delta t+3C\Delta t\right)\left(1-3\frac{A+D\Delta t}{F}\right)\delta_{ij}}, \end{array}$
and $$\mathcal{I}_{ijkl}$$ and $$\mathcal{I}^{\textrm{vol}}_{ijkl}$$ denote the fourth order symmetric identity tensor and the fourth order volumetric identity tensor respectively.

## 4. Shape memory alloys

The shape memory alloy constitutive law that is implemented in SMART+ is based on the model developed by Lagoudas et al (2012). Twenty parameters are required:

1. The Young moduli of martensite, $$E^M$$, and austenite, $$E^A$$,
2. the Poisson ratios of martensite, $$\nu^M$$, and austenite, $$\nu^A$$,
3. the coefficient of thermal expansions of martensite, $$\alpha^M$$, and austenite, $$\alpha^A$$,
4. the minimum, $$H^{\textrm{min}}$$, and the ultimate, $$H^{\textrm{sat}}$$, transformation strains,
5. the critical equivalent von Mises stress, $$\overline{\sigma}^{\textrm{crit}}$$, and the transformation strain exponent, $$k$$,
6. the four transformation temperatures of martensitic start, $$T^{M_s}$$, martensitic finish, $$T^{M_f}$$, austenitic start, $$T^{A_s}$$, and austenitic finish, $$T^{A_f}$$,
7. the forward, $$C^{M}$$, and the reverse, $$C^{A}$$, transformation slopes,
8. the four exponents, $$n_1$$, $$n_2$$, $$n_3$$, $$n_4$$, related with the smoothness of the transformation.

According to this model, the total strain $$\boldsymbol{\varepsilon}^{\textrm{tot}}$$, is decomposed in three parts, the elastic, $$\boldsymbol{\varepsilon}^{\textrm{el}}$$, the thermal $$\boldsymbol{\varepsilon}^{\textrm{th}}$$ and the transformation strain, $$\boldsymbol{\varepsilon}^{\textrm{tr}}$$,
${\varepsilon}^{\textrm{tot}}_{ij}={\varepsilon}^{\textrm{el}}_{ij}+{\varepsilon}^{\textrm{th}}_{ij}+{\varepsilon}^{\textrm{tr}}_{ij}.$
The elastic strain is connected with the stress though the linear relation
${\varepsilon}^{\textrm{el}}_{ij}={M_{ijkl}}{\sigma}_{kl}.$
$${M}_{ijkl}={L}^{-1}_{ijkl}$$ is the elastic compliance tensor of the SMA, which depends on the martensitic voloume fraction $$\xi$$ via a rule of mixtures relation
${M}_{ijkl}=\xi{M}^M_{ijkl}+(1-\xi){M}_{ijkl}^{A}.$
The elastic compliance tensors of martensite, $${M}_{ijkl}^{M}$$, and austenite, $${M}_{ijkl}^{A}$$, are computed like typical elastic isotropic materials. The thermal strain is connected with the temperature via the relation
${\varepsilon}^{\textrm{th}}_{ij}=\left(\xi{\alpha}^M_{ij}+(1-\xi){\alpha}_{ij}^{A}\right)T-{\alpha}_{ij}^{A}T^{\textrm{ref}},$
where $$T^{\textrm{ref}}$$ is a reference temperature, usually equal to the temperature at the beginning of the analysis. The coefficients of thermal expansion tensors of martensite, $${\alpha}_{ij}^{M}$$, and of austenite, $${\alpha}_{ij}^{A}$$, are computed like typical elastic isotropic materials.

With regard to the phase transformation, the term forward transformation denotes the process in which the austenite is transformed into martensite $$\left(\dot{\xi}>0\right)$$, while the term reverse transformation denotes the process in which the martensite is transformed into austenite $$\left(\dot{\xi}<0\right)$$. The transformation strain is provided through an evolution law of the form \begin{aligned} \dot{{\varepsilon}}^{\textrm{tr}}_{ij}={\Lambda}_{ij}\dot{\xi}\,, && {\Lambda}_{ij}=\left\{ \begin{array}{ll} \sqrt{\frac{3}{2}}H^{\textrm{cur}}\frac{{\sigma}_{ij}'}{\overline{\sigma}}\,, & \dot{\xi}>0\,, \\ \frac{{\varepsilon}^{t-r}_{ij}}{\xi^{t-r}}\,, & \dot{\xi}<0\,, \end{array}\right.\end{aligned} where $$\sigma'_{ij}=\sigma_{ij}-\frac{1}{3}\sigma_{kk}\delta_{ij}$$ is the deviatoric part of the stress, $$\delta_{ij}$$ is the Kronecker delta, $$\overline{\sigma}=\sqrt{\frac{3}{2}\sigma’_{kl}\sigma’_{kl}}$$ is the von Mises equivalent stress, $${\varepsilon}_{ij}^{t-r}$$ is the transformation strain tensor at the end of the forward transformation, $$\xi^{t-r}$$ is the martensitic volume fraction at the end of the forward transformation (for complete transformation loops $$\xi^{t-r}$$=1). Moreover, $$H^{\textrm{cur}}$$ is the stress dependent maximum transformation strain, provided by the formula
$H^{\textrm{cur}}=\left\{ \begin{array}{ll} H^{\textrm{min}}, & \overline{\sigma}\leq\overline{\sigma}^{\textrm{crit}}, \\ H^{\textrm{min}}+\left(H^{\textrm{sat}}-H^{\textrm{min}}\right)\left[1-\textrm{exp}\left(-k\left(\overline{\sigma}-\overline{\sigma}^{\textrm{crit}}\right)\right)\right], & \overline{\sigma}>\overline{\sigma}^{\textrm{crit}}. \end{array}\right.$
Both forward and reverse transformation processes are limited by the value of $$\xi$$, which should be between 0 and 1. For the forward and reverse transformation processes appropriate plasticity-type yield surfaces are identified $$\Phi^{\textrm{forward}}$$ and $$\Phi^{\textrm{reverse}}$$. The appropriate Kuhn-Tucker conditions are written as
$\begin{matrix} \Phi^{\textrm{forward}}=\pi-(Y+D\sigma_{ij}{\Lambda}_{ij})\leq 0,~~~ \dot{\xi}\geq0,~~~ \dot{\xi}~\Phi^{\textrm{forward}}=0, \\ \Phi^{\textrm{reverse}}=-\pi-(Y+D\sigma_{ij}{\Lambda}_{ij})\leq 0,~~~ \dot{\xi}\leq0,~~~ \dot{\xi}~\Phi^{\textrm{reverse}}=0, \end{matrix}$
where
$\begin{matrix}\begin{array}{r}\pi=\sigma_{ij}\Lambda_{ij}+\frac{1}{2}\sigma_{ij}\left({M}^M_{ijkl}-{M}^A_{ijkl}\right)\sigma_{kl}+\sigma_{ij}\left({\alpha}^M_{ij}-{\alpha}^A_{ij}\right)\left(T-T^{\textrm{ref}}\right)\\+\rho\Delta s_0 T-\rho\Delta u_0-f,\end{array}\\ f=\left\{\begin{array}{ll} \frac{1}{2}a_1\left(1+\xi^{n_1}-(1-\xi)^{n_2}\right)+a_3, & \dot{\xi}>0, \\ \frac{1}{2}a_2\left(1+\xi^{n_3}-(1-\xi)^{n_4}\right)-a_3, & \dot{\xi}<0, \end{array}\right.\end{matrix}$ and $$\rho\Delta s_0$$, $$\rho\Delta u_0$$, $$Y$$, $$D$$, $$a_1$$, $$a_2$$, $$a_3$$ are material parameters, connected with the transformation temperatures $$T^{M_s}$$, $$T^{M_f}$$, $$T^{A_s}$$, $$T^{A_f}$$, the transformation slopes $$C^{M}$$, $$C^{A}$$, and the four exponents, $$n_1$$, $$n_2$$, $$n_3$$, $$n_4$$ (for further details, see Lagoudas et al, 2012). In SMART+ the shape memory alloy constitutive law is implemented using the convex cutting plane return mapping algorithm. The updated stress is provided for 1D, plane stress, and generalized plane strain/3D analysis according to the forms of elastic isotropic materials.

## 5. Return mapping algorithms

An inelastic, rate independent, material with a plasticity-type behavior is described by a constitutive law which is given by the set of equations
$\begin{matrix} {\sigma}_{ij}=L_{ijkl}\left({\varepsilon}^{\textrm{tot}}_{ij}-{\varepsilon}^{\textrm{in}}_{ij}\right),\\ \dot{\varepsilon}^{\textrm{in}}_{ij}=\dot{p}\Lambda_{ij},~~~\Lambda_{ij}=\Lambda_{ij}(\boldsymbol{\sigma}),\\ \Phi=\Phi(\boldsymbol{\sigma},p)\leq 0,~~~ \dot{p}\geq0,~~~ \dot{p}~\Phi=0, \end{matrix}$
where a dot above a scalar, vector or tensor denotes time derivative of the scalar, the vector or the tensor respectively. The scalar $$p$$ is an inelastic multiplier, $$\Lambda_{ij}$$ is usually stress dependent and denotes direction of the inelastic strain evolution and $$\Phi$$ is a scalar function of yield-type that dictates whenever inelastic strains are generated. For the computational implementation of the above nonlinear material in a finite element framework, SMART+ utilizes the strain-driven return mapping algorithms (Ortiz and Simo, 1986). The general procedure of these algorithms is described as follows:

1. In the first step, usually called “elastic prediction”, the inelastic strains are assumed not to evolve and only generation of elastic strains are considered. Thus, during this step, the strain increment provided by the solver of the equilibrium equations (for instance, the SMART+ solver or a finite element analysis) is considered to be elastic, i.e. $$\dot{\varepsilon}_{ij}^{\textrm{tot}}=\dot{\varepsilon}_{ij}^{\textrm{el}}$$.
2. In the second step, the “inelastic correction”, the total strain is fixed, and the error in the stresses is corrected by developing inelastic strains. Thus, during this step $${\varepsilon}^{\textrm{tot}}_{ij}$$ does not evolve, while $$\dot{\varepsilon}_{ij}^{\textrm{in}}\neq 0$$.

The second step requires the numerical integration of the evolution equations for the inelastic strains and the identification of the tangent modulus that is necessary for the finite element framework. In case of temperature dependent material constitutive laws, the temperature increment is considered constant during the inelastic correction (further details about the role of temperature in a return mapping algorithm can be found in Lagoudas, 2008).

The return mapping algorithms used in SMART+ are based on a backward Euler fully implicit numerical scheme, where the value of a given quantity $$x$$ is updated from the previous time step $$n$$ to the current $$n+1$$ by the formula
$x^{(n+1)}=x^{(n)}+\Delta x^{(n+1)}.$
For a nonlinear processes, the above expression can be solved only iteratively. During such a process, the current value is updated from iteration to iteration by the expression
$x^{(n+1)(k+1)}=x^{(n+1)(k)}+\delta x^{(n+1)(k)},~~~ \delta x^{(n+1)(k)}=\Delta x^{(n+1)(k+1)}-\Delta x^{(n+1)(k)}.$
The iteration process finishes when $$\delta x^{(n+1)(k)}$$ converges to zero.

When designing a constitutive law for a finite element software like ABAQUS (UMAT routine), the main idea is that at the time increment $$n$$ the analysis is completed and the finite element code provides i) all the variables (total strain, stress, inelastic strains etc) at the time increment $$n$$ and ii) the increment of total strain $$\Delta\varepsilon^{\textrm{tot}(n+1)}_{ij}$$. The aim of a UMAT routine is to compute a) the total stress and the inelastic variables at time step $$n+1$$, and b) the tangent stiffness tensor $$L^t_{ijkl}={\partial\sigma_{ij}}/{\partial\varepsilon^{\textrm{tot}}_{kl}}$$ for the next finite element calculation.

Using the Backward Euler framework the stress-strain relation is written in incremental form at time step $$n+1$$ as
$\Delta{\sigma}_{ij}^{(n+1)}=L_{ijkl}\left(\Delta{\varepsilon}_{kl}^{\textrm{tot}(n+1)}-\Delta{\varepsilon}_{kl}^{\textrm{in}(n+1)}\right)\,.$
Since for each loading step the total strain increment is provided by the global solver and is thus known, it is the role of the inelastic prediction algorithm to find the current stress. During iterative correction, the total strain at the current time step is considered constant such that $$\delta{\varepsilon}^{\textrm{tot}(n+1)(k)}=0$$. For the scalar-valued internal variable $$p$$ holds
$p^{(n+1)(k+1)}=p^{(n+1)(k)}+\delta p^{(n+1)(k)},~~~\Delta p^{(n+1)(k+1)}=\Delta p^{(n+1)(k)}+\delta p^{(n+1)(k)}.$

From this point several approaches exist in the literature of nonlinear materials. In SMART+ two options are considered: the convex cutting plane and the closest point projection. These two methodologies have been discussed in detail by Ortiz and Simo (1986), Simo and Hughes (1998) for plastic and viscoplastic materials, and by Qidwai and Lagoudas (2000) for shape memory alloys.

### 5.1 Convex cutting plane

In the convex cutting plane the inelastic strains are updated using the following incremental formulation:
$\Delta{\varepsilon}_{ij}^{\textrm{in}(n+1)(k+1)}=\Delta p^{(n+1)(k+1)}~{{\Lambda}_{ij}^{(n+1)(k)}}.$
In this scheme the main assumption is that the direction of the inelastic strain rate remains the same between two consequent iteration steps. Thus the iterative change in stress during correction is computed by
$\delta{\sigma}_{ij}^{(n+1)(k)}=-\delta{p}^{(n+1)(k)}~{L_{ijkl}}{{\Lambda}_{kl}^{(n+1)(k)}}.$
Linearization of the constraint equation $$\Phi$$ about its roots yields
$\begin{array}{rl}\Phi^{(n+1)(k+1)}=&{\Phi^{(n+1)(k)}}+\delta {\Phi^{(n+1)(k)}} \\=&\displaystyle{\Phi^{(n+1)(k)}+\frac{\partial\Phi^{(n+1)(k)}}{\partial{\sigma}_{ij}}\delta{\sigma}_{ij}^{(n+1)(k)}+\frac{\partial\Phi^{(n+1)(k)}}{\partial p}\delta p^{(n+1)(k)}\simeq0}.\end{array}$
Combining the last two relations one obtains
$\delta {p}^{(n+1)(k)}=\frac{-{\Phi^{(n+1)(k)}}}{\mathcal{A}^{(n+1)(k)}},~~~\mathcal{A}^{(n+1)(k)}=-\frac{\partial\Phi^{(n+1)(k)}}{\partial{\sigma}_{ij}}{L_{ijkl}} {{\Lambda}_{kl}^{(n+1)(k)}}+\frac{\partial\Phi^{(n+1)(k)}}{\partial p}.$
The inelastic strains and the stresses can then be updated as follows:
$\begin{array}{rl}{\varepsilon}_{ij}^{\textrm{in}(n+1)(k+1)}=&{\varepsilon}_{ij}^{\textrm{in}(n+1)(k)}+\delta {p}^{(n+1)(k)}{{\Lambda}_{ij}^{(n+1)(k)}}, \\ {\sigma}_{ij}^{(n+1)(k+1)}=&L_{ijkl}\left({{\varepsilon}_{ij}^{\textrm{tot}(n)} +\Delta{\varepsilon}_{ij}^{\textrm{tot}(n)}-{\varepsilon}_{ij}^{\textrm{in}(n+1)(k+1)}}\right).\end{array}$
The iterative scheme continues until $${\Phi^{(n+1)(k+1)}}$$ is smaller than some tolerance.

#### 5.1.1 Continuum tangent modulus

In addition to computing the updated stress and internal variables, the global finite element solver also requires the tangent modulus $$\boldsymbol{L}^t$$ which defines the current rate of change of stress with a change in total strain. Coupled thermal and mechanical analysis also requires the thermal tangent modulus. It is explained in the literature (Simo and Hughes, 1998) that the convex cutting plane form of the return mapping algorithm utilizes the continuum tangent modulus.

In the continuum tangent modulus the scope is to identify the relation between $$\textrm{d}{\sigma}_{ij}$$ and $$\textrm{d}{\varepsilon}_{ij}$$ using the continuum description of the algorithm. When inelastic strains are not generated, $${L}_{ijkl}^t={L}_{ijkl}$$. When inelastic strains are generated, the stress-strain relation is rewritten in differential form and the evolution equations are substituted. Rewriting in terms of stress, this gives
$\textrm{d}{\sigma}_{ij}=L_{ijkl}\left(\textrm{d}{\varepsilon}^{\textrm{tot}}_{kl}-\textrm{d}{p}{\Lambda}_{kl}\right).$
The constraint equation $$\Phi$$ must be satisfied for all acceptable solutions. Taking its differential results in a relation analogous to the classical “consistency condition”
$\textrm{d}{\Phi}=\frac{\partial\Phi}{\partial\sigma_{ij}}\textrm{d}{{\sigma}_{ij}}+\frac{\partial\Phi}{\partial p}\textrm{d}{p}=0.$
Combining the last two expressions yields
$\begin{array}{c}\displaystyle{\textrm{d}{p}=-\frac{\partial\Phi}{\partial{\sigma}_{ij}}{L_{ijkl}}\textrm{d}{\varepsilon}^{\textrm{tot}}_{kl}\big/\mathcal{A}, ~~~\mathcal{A}=-\frac{\partial\Phi}{\partial{\sigma}_{ij}}{L_{ijkl}}{{\Lambda}_{kl}}+\frac{\partial\Phi}{\partial p},}\\ \displaystyle{\textrm{d}{\sigma}_{ij}=L_{ijkl}^t\textrm{d}{\varepsilon}^{\textrm{tot}}_{kl}, ~~~L_{ijkl}^t=L_{ijkl}+ \frac{1}{\mathcal{A}}{L_{ijmn}}{\Lambda}_{mn}\frac{\partial\Phi}{\partial{\sigma}_{pq}}{L_{pqkl}}.}\end{array}$

### 5.2 Closest point projection

In the closest point projection the inelastic strains are updated using the following incremental formulation:
$\Delta{\varepsilon}_{ij}^{\textrm{in}(n+1)(k)}=\Delta p^{(n+1)(k)}~{{\Lambda}_{ij}^{(n+1)(k)}}=\left(p^{(n+1)(k)}-p^{(n)}\right)~{{\Lambda}_{ij}^{(n+1)(k)}}.$
The last expression can be rewritten with the help of a residual $$\boldsymbol{R}$$ as
${R}_{ij}^{(n+1)(k)}=-\Delta{\varepsilon}_{ij}^{\textrm{in}(n+1)(k)}+\Delta{p}^{(n+1)(k)}{{\Lambda}_{ij}^{(n+1)(k)}}.$
Considering that $$\boldsymbol{\Lambda}$$ is only a function of $$\boldsymbol{\sigma}$$, the linearized form of the residual can be written as
$\begin{array}{r}{R}_{ij}^{(n+1)(k+1)}={R}_{ij}^{(n+1)(k)}+\delta{R}_{ij}^{(n+1)(k)}={R}_{ij}^{(n+1)(k)}-\delta{\varepsilon}_{ij}^{\textrm{in}(n+1)(k)} +\delta{p}^{(n+1)(k)}{{\Lambda}_{ij}^{(n+1)(k)}}\\ \displaystyle{+\Delta{p}^{(n+1)(k)}\frac{\partial{\Lambda}_{ij}^{(n+1)(k)}}{\partial\sigma_{kl}}\delta{\sigma}_{kl}^{(n+1)(k)}\simeq0}.\end{array}$
Moreover, the iterative change in stress during correction is computed by
$\delta{\sigma}_{ij}^{(n+1)(k)}=-{L_{ijkl}}{\delta{\varepsilon}_{kl}^{\textrm{in}(n+1)(k)}}.$
Finally, linearization of the constraint equation $$\Phi$$ about its roots yields
$\begin{array}{rl}\Phi^{(n+1)(k+1)}=&{\Phi^{(n+1)(k)}}+\delta {\Phi^{(n+1)(k)}}\\ =&\displaystyle{\Phi^{(n+1)(k)}+\frac{\partial\Phi^{(n+1)(k)}}{\partial{\sigma}_{ij}}\delta{\sigma}_{ij}^{(n+1)(k)}+\frac{\partial\Phi^{(n+1)(k)}}{\partial p}\delta p^{(n+1)(k)}\simeq0.}\end{array}$
Combining the last three equations and after proper calculations, the numerical algorithm yields
$\begin{array}{c}\displaystyle{\delta {p}^{(n+1)(k)}= \left(-{\Phi^{(n+1)(k)}}+\frac{\partial\Phi^{(n+1)(k)}}{\partial{\sigma}_{ij}}{\mathcal{E}^{(n+1)(k)}_{ijkl}}{{R}_{kl}^{(n+1)(k)}} \right)\big/{\mathcal{F}^{(n+1)(k)}},}\\ \displaystyle{\mathcal{E}_{ijkl}^{(n+1)(k)}=\left(L^{-1}_{ijkl}+\Delta p^{(n+1)(k)}\frac{\partial{\Lambda}_{ij}^{(n+1)(k)}}{\partial\sigma_{kl}}\right)^{-1},}\\ \displaystyle{\mathcal{F}^{(n+1)(k)}=-\frac{\partial\Phi^{(n+1)(k)}}{\partial{\sigma}_{ij}}{\mathcal{E}^{(n+1)(k)}_{ijkl}}{{\Lambda}_{kl}^{(n+1)(k)}}+ \frac{\partial\Phi^{(n+1)(k)}}{\partial p}.}\end{array}$
The inelastic strains and the stresses can then be updated as follows:
$\begin{array}{rl}{\varepsilon}_{ij}^{\textrm{in}(n+1)(k+1)}=&{\varepsilon}_{ij}^{\textrm{in}(n+1)(k)}+ L^{-1}_{ijkl}\mathcal{E}^{(n+1)(k)}_{klmn}\left({R}_{mn}^{(n+1)(k)}+\delta {p}^{(n+1)(k)}{{\Lambda}_{mn}^{(n+1)(k)}}\right), \\{\sigma}_{ij}^{(n+1)(k+1)}=&L_{ijkl}\left({{\varepsilon}_{ij}^{\textrm{tot}(n)} +\Delta{\varepsilon}_{ij}^{\textrm{tot}(n)}-{\varepsilon}_{ij}^{\textrm{in}(n+1)(k+1)}}\right),\end{array}$
The iterative scheme continues until $${\Phi^{(n+1)(k+1)}}$$ is smaller than some tolerance.

#### 5.2.1 Consistent tangent modulus

In addition to computing the updated stress and internal variables, the global finite element solver also requires the tangent modulus $$\boldsymbol{L}^t$$ which defines the current rate of change of stress with a change in total strain. Coupled thermal and mechanical analysis also requires the thermal tangent modulus. It is explained in the literature (Simo and Hughes, 1998) that the closest point projection form of the return mapping algorithm utilizes the consistent tangent modulus.

In the consistent tangent modulus the scope is to identify the relation between $$\textrm{d}{\sigma}^{(n+1)}_{ij}$$ and $$\textrm{d}{\varepsilon}^{(n+1)}_{ij}$$ by direct differentiation of the numerical algorithm. When inelastic strains are not generated, $${L}_{ijkl}^t={L}_{ijkl}$$. When inelastic strains are generated, the system of equations that describe the material behavior at time step $$n+1$$ are written in discrete form as
$\begin{array}{c}{\sigma}^{(n+1)}_{ij}=L_{ijkl}\left({\varepsilon}^{\textrm{tot}(n+1)}_{kl}-{\varepsilon}^{\textrm{in}(n+1)}_{kl}\right),\\ {\varepsilon}^{\textrm{in}(n+1)}_{ij}={\varepsilon}^{\textrm{in}(n)}_{ij}+\left(p^{(n+1)}-p^{(n)}\right)~{{\Lambda}_{ij}^{(n+1)}},\\ \Phi^{(n+1)}=0. \end{array}$
Numerical differentiation of the above equations yields
$\begin{array}{rl}\textrm{d}{\sigma}^{(n+1)}_{ij}=&L_{ijkl}\left({\textrm{d}\varepsilon}^{\textrm{tot}(n+1)}_{kl}-\textrm{d}{\varepsilon}^{\textrm{in}(n+1)}_{kl}\right),\\ {\textrm{d}\varepsilon}^{\textrm{in}(n+1)}_{ij}=&\displaystyle{\frac{\partial{\varepsilon}^{\textrm{in}(n+1)}_{ij}}{\partial{\sigma}^{(n+1)}_{kl}}\textrm{d}{\sigma}^{(n+1)}_{kl}+\frac{\partial{\varepsilon}^{\textrm{in}(n+1)}_{ij}}{\partial p^{(n+1)}}\textrm{d}p^{(n+1)}}\\=& \displaystyle{\left(p^{(n+1)}-p^{(n)}\right)~\frac{\partial{\Lambda}_{ij}^{(n+1)}}{\partial{\sigma}^{(n+1)}_{kl}}\textrm{d}{\sigma}^{(n+1)}_{kl}+{\Lambda}_{ij}^{(n+1)}\textrm{d}p^{(n+1)},}\\ \textrm{d}{\Phi}^{(n+1)}=&\displaystyle{\frac{\partial\Phi^{(n+1)}}{\partial\sigma^{(n+1)}_{ij}}\textrm{d}{{\sigma}^{(n+1)}_{ij}}+\frac{\partial\Phi^{(n+1)}}{\partial p^{(n+1)}}\textrm{d}{p}^{(n+1)}=0.}\end{array}$
Combining the last three equations and after proper calculations gives
$\begin{array}{c}\displaystyle{\textrm{d}{p}^{(n+1)}=-\frac{\partial\Phi^{(n+1)}}{\partial{\sigma}^{(n+1)}_{ij}}{\mathcal{E}^{(n+1)}_{ijkl}}\textrm{d}{\varepsilon}^{\textrm{tot}(n+1)}_{kl}\big/\mathcal{F}^{(n+1)},}\\ \displaystyle{\mathcal{E}_{ijkl}^{(n+1)}=\left(L^{-1}_{ijkl}+\left(p^{(n+1)}-p^{(n)}\right)\frac{\partial{\Lambda}_{ij}^{(n+1)}}{\partial\sigma^{(n+1)}_{kl}}\right)^{-1},}\\ \displaystyle{\mathcal{F}^{(n+1)}=-\frac{\partial\Phi^{(n+1)}}{\partial{\sigma}^{(n+1)}_{ij}}{\mathcal{E}^{(n+1)}_{ijkl}}{{\Lambda}_{kl}^{(n+1)}}+ \frac{\partial\Phi^{(n+1)}}{\partial p^{(n+1)}},}\\ \displaystyle{\textrm{d}{\sigma}^{(n+1)}_{ij}=L_{ijkl}^{t(n+1)}\textrm{d}{\varepsilon}^{\textrm{tot}(n+1)}_{kl},}\\ \displaystyle{L_{ijkl}^{t(n+1)}=\mathcal{E}^{(n+1)}_{ijkl}+ \frac{1}{\mathcal{F}^{(n+1)}} \mathcal{E}^{(n+1)}_{ijmn}\Lambda^{(n+1)}_{mn}\frac{\partial\Phi^{(n+1)}}{\partial{\sigma}^{(n+1)}_{pq}}\mathcal{E}^{(n+1)}_{pqkl}.}\end{array}$