A Shell Finite Element Model for Superelasticity of Shape Memory Alloys

A finite element formulation for the analysis of large strains of thinwalled shape memory alloys is briefly presented. For the shell model we use a seven-kinematic-parameter model for large deformations and rotations, which takes into account the through-the-thickness stretch and can directly incorporate a fully 3D inelastic constitutive equations. As for the constitutive model, we use a large strain isotropic formulation that is based on the multiplicative decomposition of the deformation gradient into the elastic and the transformation part and uses the transformation deformation tensor as an internal variable. Numerical examples are presented to illustrate the approach.

When a SMA in the M phase is subjected to mechanical loading, it deforms and a detwinning process of martensite variants occurs. At the macroscopic scale, this is a change of shape, while at the micro-scale, martensite variants are oriented in a more preferable way. During unloading the martensite variants do not change, therefore pseudo-plastic deformations are present at zero load. At this stage, shape recovery can be achieved by subjecting the SMA material to a temperature above the austenitic transformation finish temperature (A f ), where only A phase is stable. During the temperature rise, a phase transformation from M to A begins at the austenitic transformation start temperature (A s ) and is completed at A f . On the other hand, cooling of SMA causes martensitic transformation with formation of self-accommodated martensite variants. The transformation begins at martensitic transformation start temperature (M s ) and ends at martensitic transformation finish temperature (M f ). In contrast to heating, where a shape recovery occurs, cooling of SMA to its initial temperature does not cause any change in shape. It is worth noting that shape recovery is only possible up to a certain degree of deformation.
i) the ability to remember its original shape when a deformed SMA part is subjected to a high temperature (called the shape memory effect) and ii) the ability to withstand large strains (up to 8%) without permanent plastic deformation (known as superelasticity), see Fig. 20.1 (left). Both properties are attributed to the fact that SMAs are found in two different phases, Fig. 1 (right). The high temperature parent phase is called austenite or austenitic phase (A) and the low temperature product phase is called martensit or martensitic phase (M). The crystal structure of A is highly symmetric, which is why it can be only in one variant. On the other hand, M can be found in a large number of variants due to the lower symmetry of its crystal structure. Twinned or self accommodated martensite (M t ), which is stable at low stress state, occurs in different variants. At a stress state higher than the critical, a process of detwinning starts, which results in detwinned (or stress induced) martensite (M d ), i.e. the variant with the preferable orientation for a given stress state. In the phase diagram in Fig.1 (right), areas with stable crystal lattices and transformation regions are shown schematically. a large deformation. During unloading an endothermic reverse martensite transformation from martensite to austenite takes place and the deformation vanishes. The mechanism is not elastic because a transformation (change of the crystal lattice) takes place. Depending on the strain rate and heat exchange between the SMA material and the environment, the material may heat up or cool down during the transformation.
A number of 3D SMA material models were proposed. For implementation within the framework of the finite element (FE) method, macroscale phenomenological SMA models are a preferable choice. These models differ in various aspects, but the biggest difference is whether they are designed to solve small strain or large strain problems. The large strain models assume a multiplicative decomposition of the deformation gradient. One of the first large strain SMA models was developed by Auricchio and Taylor (1997) and numerically implemented in Auricchio (2001). Stupkiewicz and Petryk (2013) showed how to reformulate a small strain model for the finite strain regime, for superelasticity with tension-compression asymmetry and anisotropy. In Reese and Christ (2008) the deformation gradient is decomposed into elastic and transformation parts and the latter is further divided into a recoverable and a plastic part.
In this work, we apply a version of the large strain SMA model proposed in Arghavani et al (2011);Souza et al (1998); Evangelista et al (2010). The model can predict the superelastic response and the shape memory effect of polycrystalline SMA. It assumes isothermal transformations and neglects the tension-compression asymmetry as well as functional fatigue due to cyclic loading. In the framework of the finite element method, the 3D SMA material models are usually incorporated into the 3D solid finite elements. A FE implementation of an SMA model in a plate or shell finite element formulation is very rare. One of the reasons is that the standard shell theories of Kirchhoff and Reissner-Mindlin type cannot directly include 3D constitutive models because the plane stress constraint has to be enforced (and this is not a trivial task for an inelastic model, see e.g. Dujc and Brank (2012)). There are, however, the 3D-shell finite elements and the solid-shell finite elements, which are designed in such a way that they can directly use 3D constitutive equations without modification, see for example Brank et al (2002); Brank (2005) ;Brank et al (2008). In this work we rely on a 3D shell model proposed in Brank (2005). Our numerical formulation, which is presented below, can be used to simulate the nonlinear behaviour (due to mechanical loading) of (very) thin-walled shape memory alloys. They can undergo large deformations, large rotations (the formulation described in Brank and Ibrahimbegovic (2001); Ibrahimbegovic et al (2001) is used) and large strains.

Constitutive Model for SMA
In this section we revisit the 3D constitutive model for SMA that was originally developed by Souza et al (1998) for small strains and later extended to finite strains by Evangelista et al (2010); Arghavani et al (2011). Similar to the finite strain Superelasticity (also called pseudoelasticity) is exhibited for a temperature above A f , where the loading/unloading response is characterised by a nonlinear behaviour with hysteresis. During loading, a stress-induced martensite is formed during the exothermic martensite transformation. At the macroscopic level, this is observed as plasticity (see e.g. Ibrahimbegovic, 2009), a multiplicative decomposition of the deformation gradient into elastic and transformation parts is postulated: With Eq. (20.1), the initial (i.e. undeformed), intermediate and current (i.e. deformed) configurations are introduced. In this work, the total Lagrangian formulation is used to describe large deformations of a shell, which requires the derivation of constitutive equations with respect to the initial configuration. To achieve this goal, however, we also use the tensors defined at the intermediate configuration. One such tensor is C e = F T e F e , while C t = F T t F t is defined at the initial configuration. The Cauchy-Green deformation tensor and the Green-Lagrange strain tensor are: respectively, where 1 is the unit tensor. The velocity gradient tensor L = F F -1 and its symmetric part (i.e. rate of deformation) d = 1 2 (L + L T ) are also used. The dot stands for the (pseudo-)time derivative. The relation between d and the strain rate E is: According to experimental observations, the transformation in SMA is (almost) isochronic, which is expressed by det(F t ) = 1 that yields tr(d t ) = 0. The Helmholtz free energy ψ must depend on F e only through C e in order to satisfy material objectivity. It is assumed that ψ depends also on C t or yet the transformation strain tensor temperature T, and that it can be additively decomposed into elastic part ψ e and transformation part ψ t : We assume material isotropy and choose a neo-Hooke type of hyperelastic strain energy function: where μ and λ are Lamé's coefficients, and I 1 and I 3 are the first and third invariants of C e , respectively. It can be shown that these invariants equal the first and third invariants of C C -1 t , which enables writing the strain energy function (20.6) in terms of tensors from the initial configuration, i.e. ψ e (C e ) =ψ e (C,C t ). The transformation part of free energy is chosen as (after Arghavani et al, 2011;Souza et al, 1998): Here, τ M (T) = β T − T 0 provides temperature dependency of material response, h, β and T 0 are the material parameters of the SMA model, · are Macaulay brackets and · is standard tensor norm. In Eq. (20.7), a step function is introduced in order to enable a (computational) enforcement of the limit of transformation strains ε L . This is another material parameter that can be obtained experimentally from the uniaxial test as the absolute value of the maximum transformation strain.
The Clausius-Duhem inequality form of the second law of thermodynamics states: where S is the second Piola-Kirchoff stress tensor and η is entropy. Substituting (20.5) into (20.9), we obtain: It can be shown that 1 can be expanded as: For an elastic case with no change in transformation and temperature, the expression for the stress tensor follows from (20.10) and (20.11) as: For a case of a temperature change with no change in transformation, the expression for entropy is obtained from (20.10) as: It is assumed that the relations (20.12) and (20.13) are also valid in the case of transformation, which is the only case that still has to be considered. To this end, it can be shown that the derivatives in 2 can be expressed as: is a normalized transformation tensor, and is a penalty-like parameter that results from sub-differential of the indicator function Using equations (20.12), (20.13) and (20.14), the dissipation at the transformation case is obtained from (20.10) as: Some mathematical manipulations and (20.3) yield the equality X : Moreover, a decomposition of L t into a symmetric and a skew-symmetric parts, d t and w t , respectively, leads to equality P : L t = P : d t . By further using the notation K = F t X F T t and , where K plays a role of the back-stress tensor, the inequality of the dissipation at transformation can be rewritten as: It can be shown by some mathematical manipulations that Y , which plays the role of an effective stress tensor, can be expressed in terms of tensors defined at initial configuration as: (20.17) Z By choosing a transformation function f ( ) ≤ 0 and assuming that the transformation case corresponds to the stress state giving f ( ) = 0, the evolution equation for the transformation case can be obtained. We postulate that among all admissible states of transformation, we choose the one that renders maximum dissipation D t or yet minimum of −D t . Recasting the original problem into the unconstrained minimization problem can be done by using the Lagrange multiplier method and by defining the Lagrange function as: where ζ is Lagrange multiplier. In order to find minimum of L ( , ζ), the stationarity, primal feasibility, dual feasibility and complementary slackness conditions (also known as Kuhn-Tucker conditions and loading/unloading conditions) must hold: which can be replaced (in a consistent manner) by another evolution equation, which is more suitable for numerical implementation based on the total Lagrangian formulation: (20.24) This concludes a short derivation of the SMA material model used in this work. The equations that define the above SMA large strain constitutive model are: the chosen Helmholtz free energy function defined by (20.5), (20.6) and (20.7); the stress tensors that can be computed by (20.12), (20.14) and (20.18); the transformation function (20.22); the loading/unloading conditions (20.20); and the evolution equation (20.24). We will omit discussing its numerical implementation.

Seven-parameter Shell Model
In this section, we briefly present a 3D-shell model that can incorporate a 3D constitutive model without a modification. Let the position vector to the material point of the initial shell configuration be defined as: where ξ 1 and ξ 2 are curvilinear coordinates that parameterize the shell mid-surface A , ξ 3 ∈ [− h 2 , h 2 ] is a straight through-the-thickness coordinate, h is the initial shell thickness and A is the shell director vector. In what follows, we will omit showing From the stationarity condition (20.20) 1 , the evolution equation for the rate of transformation deformation d t is obtained as: In this work we choose the following transformation function (according to Arghavani et al, 2011;Souza et al, 1998), which resembles classical yield functions for metals: where D is the deviatoric part of , and R is elastic region radius, which is another parameter of the material model that has to be evaluated experimentally. Finally, combining (20.21) and (20.22) yields the following form of the evolution equation: functions and functional dependency on the convective coordinates we introduced above for the sake of brevity. The position vector to the material point of the deformed shell configuration is assumed as: where u is the mid-surface displacement vector, a is the rotated shell director that preserves the original unit length and the parameters λ and q define a constant and linear through-the-thickness stretching, respectively. The rotation of A into a is described by a singularity-free formulation that uses two large rotation parameters, i.e., a = a( θ 1 ,θ 2 ) ; we refer to e.g. Brank and Ibrahimbegovic (2001); Ibrahimbegovic et al (2001) for details. Equation (20.26) introduces seven kinematic parameters of the adopted shell model, which can be collected in a vector as Φ = {u,θ 1 ,θ 2 ,λ,q}.
The components of the Green-Lagrange strain tensor can be defined with respect to the above introduced curvilinear coordinates as: where G i is the contravariant base vector defined as G i · G k = δ i k . The covariant base vector is given as: 20.28) and δ i j is the Kronecker delta symbol. The strains E ij are polynomials up to the fourth order with respect to ξ 3 . In this work, we neglect the terms of orders three and four because, as our numerical experiments show, they have a negligible influence on the results. On the other hand, it is important that all strains have non-zero terms of order one and two. This allows the implementation of a fully 3D constitutive equations without modelling errors, which is not possible with many other shell models due to the inherent 2D character of the description of the shell kinematics.
The virtual work equation (i.e. the weak form of the equilibrium equations) is the starting point for finite element discretization. It can be written as: where G int is the virtual work of the internal forces, G ext is virtual work of the external forces acting on the shell, δE is variation of the strain tensor field that can be obtained as where is a scalar parameter, and δΦ represents the variation of the fields of seven kinematic parameters of the model. It is worth noting that the 2 nd Piola-Kirchhoff stress tensor field depends on shell kinematics as well as on the internal variable of the above described SMA constitutive model, which is C t . Here, contributes to the stiffness due to geometric effects, and is a consistent material operator, which is the fourth order tensor. After the introduction of discretization and interpolation in the framework of the finite element method, Eq. (20.29) yields a system of highly nonlinear equations (due to the arbitrariness of δΦ), where the unknowns are kinematic parameters at nodes and internal variables at Gauss integration points. The solution of such a system is based on an operator-split technology, which is applied within an incrementaliterative Newton-Raphson solution method. Namely, at each iteration the computation of the Gauss point internal variable C t (and consequently the computation of the Gauss point stresses S and the consistent material tangent operator C) is split from the computation of the nodal kinematic parameters. This is achieved by applying two sequential procedures, where the results of the first procedure are immediately used in the second procedure. Namely, the constitutive equations are first enforced at the integration points by updating the internal variables of the constitutive model. This local update is followed by the solution of the system of equations for an update of nodal kinematic parameters using the consistent tangent stiffness matrix resulting from (20.31).
The numerical examples in the next section are computed using a four node element with the assumed natural strain treatment of the transverse shear strains and through-the-thickness strains in order to avoid transverse shear and thickness lockings. The element has four integration points over the mid-surface and three through-the-thickness integration points, all of which are of the Gaussian type. In order to facilitate the implementation of the SMA constitutive model described in the previous section, a local Cartesian frame is introduced at each integration point.

Numerical Examples
The material and shell models described above were implemented into the finite element code using AceFEM (Korelc and Wriggers, 2016) by using a generator of the finite element code AceGen, see e.g. Korelc and Stupkiewicz (2014); Hudobivnik and Korelc (2016).
Three examples are presented below. The first and the second example show that our results are in in good agreement with the reference results presented in Arghavani et al (2011). The third example is somehow more demanding and illustrates the ability of the derived formulation to predict a superelastic response of a thin-walled curved structure. The material parameters typical for NiTi are applied for all considered examples: To obtain the superelastic response, the ambient temperature T is set to 37 • C. The following figures show the components of the second Piola-Kirchhoff stress tensor and the components of the Green-Lagrange strain tensor.

Square Wall Under Uniaxial Loading
A square wall with an edge length of a = 10 mm and a thickness of t = 0.01 mm is subjected to uniaxial tension and compression. The mesh (5x5 elements), the boundary conditions that allow a homogeneous uniaxial stress state over the wall and the load are shown in Fig. 20.2. The load is q 1 = λq and q 2 = 0, where q = 14 N/mm and λ is the load multiplier. Loading and unloading in tension was applied first, followed by loading and unloading in compression. Figure 20.3 shows computed uniaxial superelastic response at a Gauss point. Our results (almost) exactly match the results from reference Arghavani et al (2011). The reason for a very small discrepancy at large load levels is the use of Saint-Venant Kirchhoff strain energy function in Arghavani et al (2011), while our choice is the Neo-Hookean strain energy function.

Square Wall Under Biaxial Loading
Geometry, boundary conditions (which allow a homogenous stress state over the wall) and mesh for this example remain the same as for the previous one, see Fig. 20.2. However, the load is now q 1 = λ 1 q and q 2 = λ 2 q where q = 7 N/mm, and λ 1 and λ 2 are load multipliers. The load is therefore non-proportional, with F 1 = q 1 a and F 2 = q 2 a changing in a butterfly-like pattern as shown in Fig. 20.4. The loading was applied in five steps: (i) λ 1 and λ 2 increase from 0 to 1, (ii) λ 1 decreases to −1 at λ 2 = 1, (iii) λ 1 increases to 1 and λ 2 decreses to −1, (iv) λ 1 decreases to −1 at λ 2 = 1, and (v) λ 1 and λ 2 go to 0.
The curve in Fig. 20.5 explains how the relation between the in-plane normal strains is changing during the loading. The agreement with the solution found in Arghavani  et al (2011) is very good. We would like to point out that our load application was different from the one in Arghavani et al (2011). In our case, the mesh with edge loads was considered, while in Arghavani et al (2011) a single Gauss-point algorithm was tested under a butterfly-like stress control (with missing details how it was applied for a strain-driven update algorithm). This is the reason for a difference between our results and the reference results in Fig. 20.6, which shows the relation between the in-plane normal stresses during loading. Part of the difference may also be due to the use of different strain energy functions. It is worth noting that for the first two examples the plane stress condition is completely reproduced by the derived shell model due to the small thickness to edge ratio t/a = 10 −3 .

Compression of a Twisted Beam
Twisted beam with length L = 12 mm, width w = 2 mm and thickness t = 0.25 mm, see Fig. 20.7, is clamped at one end and subjected at the opposite end to imposed compressive axial displacements u = λu 0 , where u 0 = 1.1 mm and λ is load multiplier, which goes from 0 to 1 and back to 0. The mesh consists of 16 elements in the longitudinal and 8 elements in the transverse direction. Figure 20.8 shows deformed mesh at λ = 1, with coloured contours representing E t at the midsurface (the red colour denotes the largest transformation and the blue colour the smallest). At the final configuration, where λ is back to 0, the stresses are zero and no transformation is observed. The initial shape is fully recovered and the superelastic response is obtained. This can be nicely seen from the diagram in Fig. 20.9, which shows the reaction force at edge 1 as a function of the imposed compressive axial displacement of edge 2.

Conclusions
The aim of this work was to derive a finite element formulation that can be used for large deformation and stability analysis of curved thin-walled shape memory   alloys. It is well documented in the finite element literature that for thin shell-like structures the 3D-solid finite element models are inappropriate, because they are very likely introducing a considerable modelling error. For this purpose we have revisited a seven-parameter large deformations and large rotations shell model (that can be used for thin and thick shells), and use it as a framework for implementation of a 3D finite strain material model for shape memory alloys. Several presented numerical examples illustrate a very satisfying performance of resulting finite element formulation.
Questions on the numerical implementation of the considered SMA material model, a comparison with experimental results, further development of the finite element formulation in order to be able to perform large solution steps and to be insensitive to mesh distortions (following Lavrenčič and Brank, 2020;Brank, 2008) and the simulation of the buckling process of shape memory alloys (using Stanić and Brank, 2017) will be answered and shown in a separate publication.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.