Modeling large deformations of thin-walled SMA structures by shell ﬁnite elements

Many shape memory alloy (SMA) applications, such as biomedical devices, electromechan- ical actuators, and elastocaloric cooling devices, are based on thin-walled ﬂat or shell-like structures. An advanced design of such structures requires the development of an eﬃcient and accurate numerical tool for simulations of very thin and curved SMA structures that may experience large deformations and even buckling upon thermo-mechanical loading. So far, ﬁnite element models for ﬁnite strain deformations of SMA structures have been based on 3D solid formulations, which are relatively ineﬃcient for solving (thin) shell problems. In this paper, we present a ﬁnite element model for the analysis of shape memory alloy shells. Our model is based on a 7-parameter, large-rotation, one-director shell formulation, which takes into account a fully three-dimensional form of the constitutive equations for the isothermal transformations of isotropic superelasticity, as well as the shape-memory effect in a simpliﬁed way. In fact, we present three 4-node shell ﬁnite elements for SMAs. Two of them use the assumed natural strain concepts for the transverse shear strains, through-the-thickness normal strain, and membrane strains. The third element is a combination of the assumed natural strain and the enhanced assumed strain concepts, applied to satisfy the zero through-the-thickness-normal-stress condition for thin geometries to a high degree of accuracy. After a detailed description of the SMA ﬁnite element models for shells in the ﬁrst part of the paper, numerical examples in the second part illustrate the approach. Compared to 3D solid SMA formulations, our results show excellent accu- racy, even with a signiﬁcantly reduced number of degrees of freedom, which consequently translates into a reduction in the computational time.


Introduction
Shape memory alloys (SMAs) belong to a group of smart materials that are characterized by a change in their properties when acted upon by external stimuli (stress and temperature for SMAs). Typically, SMAs are divided into Ni-Ti-based, Cubased, and Fe-based groups, among which Ni-Ti is the most widely used in a large number of applications [1] . In medicine,  Ni-Ti is used for stents, surgical tools and bone implants due to its biocompatibility and comparable material properties with those of bone [2][3][4] . In robotics, SMAs are used as noise-free, small-size actuators with a large force-to-weight ratio [5] . In aeronautics, SMAs are employed for actuators, vibration-damping elements, seals, deployment mechanisms [6,7] , and morphing wings [8] . In civil engineering, SMAs are used as seismic-response-control and damping elements [9] . Applications in energy and process engineering are also under intensive development [10,11] , especially in elastocaloric cooling technology [12] , due to its great potential as an alternative to vapor-compression refrigeration technology, which is relatively inefficient and still uses environmentally harmful refrigerants [13] .
SMAs have two special properties, illustrated in Fig. 1 (left). First, the shape-memory effect (SME) is the ability to recover an initial shape after being deformed and subsequently subjected to high temperature, and second, superelasticity (SE) is the ability to recover from large deformations (up to 8%) without any permanent deformation. Both properties are attributed to the fact that SMAs are found in two different phases, see Fig. 1 (right): a low-temperature product phase called martensite (M) and a high-temperature parent phase called austenite (A). Because the crystal structure of A is highly symmetric, it exists in only one variant, while M can be found in many variants due to the lower symmetry of its crystal structure. Twinned or self-accommodated martensite M t , which is stable at low stresses, occurs in different variants. At a critical stress σ s the detwinned (or stress-induced) martensite M d (the variant with the preferable orientation for a given stress state) occurs. Figure 2 (left) shows the stress-strain-temperature diagram of the SME and the corresponding crystal lattices. When M is subjected to mechanical loading at low temperature, it deforms and the detwinning process occurs. On the macro-scale, this is a change of shape, while on the micro-scale, martensite variants occur that are preferable for the given stress state.
During unloading, martensite variants do not change; therefore, pseudo-plastic deformation ε SME is present at zero loads.
Shape recovery occurs by subjecting the material to a temperature above the austenitic finish temperature A f , where only A is stable. During a temperature increase, the phase transformation from M to A begins at the austenitic transformation start temperature A s and finishes at A f . On the other hand, cooling of a SMA causes a martensitic transformation (MT) with the formation of self-accommodated martensite variants M t . The transformation begins at the martensitic transformation's start temperature M s and ends at the martensitic transformation's finish temperature M f . In contrast to heating, cooling to its initial temperature does not cause any shape change (except the volume change due to thermal contraction). Shape recovery is possible if the deformation is not too large, otherwise, permanent plastic deformation occurs, as shown in Fig. 2 .
Superelasticity (also called pseudoelasticity) is exhibited by a SMA when it is mechanically loaded and unloaded above the temperature A f . A hysteresis occurs, as shown in Fig. 2 (right): while loading, stress-induced martensite M d is formed during an exothermic MT (which is observed at the macroscopic level as a large deformation), and while unloading, an endothermic reverse MT from M to A takes place and the deformation vanishes. The mechanism is not elastic because the transformation, i.e., the change of the crystal lattice, occurs. The reversible deformation is limited by the maximum transformation strain ε t . Depending on the strain rate and the heat exchange between the SMA and the environment, the SMA can heat up or cool down during the transformation, which is called the elastocaloric effect.
Numerous SMA material models have been developed, and these can be roughly grouped into three categories. Micromodels focus on phenomena such as grain nucleation and growth, detwinning, and reorientation [14,15] . Micro-macromodels combine microscopic features with macro-scale thermodynamics [16,17] , while macro-models use thermodynamic variables to describe the macro-observed phenomena. Many of the macro-models describe superelasticity and the shape memory effect for small strains [18,19] . Some include tension-compression asymmetry based on a modified transformation limit function [20,21] and the difference between austenite and martensite material properties by a nonlinear order parameter [22] . To model detwinning and reorientation, a separate treatment of self-accommodated and stress-induced martensite is required, as well as a criterion for reorientation [23][24][25][26] , while the effect of cyclic loading can be considered using an internal variable representing the accumulated martensite [27][28][29] . It is worth mentioning that not only has isothermal superelastic behavior been studied: temperature changes during a transformation were modeled by treating the released latent heat as a volume heat source [30][31][32] .
The main motivation for our work stems from the need for an accurate prediction of the response of mechanically loaded SMA structures, which could be applied in elastocaloric devices. Obviously, thin-walled structures are preferable as they enable the rapid (efficient) transfer of heat [54] . Currently, the macroscopic SMA material models have been implemented almost exclusively into 3D solid finite elements, which can be very inefficient for the modeling of thin structures. To the best of our knowledge, there have been only a few SMA formulations devised for structural elements; e.g., in [55] for plates, in [56,57] for membranes, in [58] for solid-beams, and in [59] for a small strain corotational shell formulation. There have been no attempts to combine a finite strain SMA material model with a shell finite element formulation, although it is well known that solid finite elements can be very inaccurate when representing the bending of thin-walled structures. For this reason, it seems rational to focus on the development of shell finite element formulations for the simulation of thin and curved SMA structures, especially because there exist thin-walled formulations that can incorporate fully 3D constitutive models. Let us mention 7-parameter (7-p) shells and solid-shells (e.g. [60] ), the latter being rotation-free that simplify computer implementation.
The main focus of this work is combining the finite strain SMA material model with the 7-p shell finite element model. To this end, three 7-p, 4-node, shell finite element formulations were designed with the extensive use of the assumed natural strain (ANS) concepts that were originally developed for 5-parameter shell formulations. In particular, our 7-p formulation that is called 2ANS uses the ANS concepts for the transverse-shear strains [61] and the through-the-thickness stretching [62] . The 7-p formulation that is called 3ANS additionally uses the ANS concept for the membrane strains, in particular, the approach that was proposed in [63] and extensively tested in [63][64][65][66] , is adopted. The membrane ANS makes the 3ANS formulation low-sensitive to mesh distortion. The third 7-p shell formulation is a combination of the ANS and the enhanced assumed-strain (EAS) concepts designed to better enforce the zero-through-the-thickness-normal-stress condition for thin SMAs, with respect to the 2ANS and 3ANS formulations. It must be emphasized that the resulting 7-p shell finite element formulations are a much better choice than the solid finite element formulations for the numerical simulations of very thin and curved SMA structures that can undergo large deformations, large rotations, and buckling, an understanding of which is an interest for some SMA applications, e.g., seismic applications and elastocaloric cooling technology.
The rest of the text is organized in the following way. Section 2 revises the finite strain SMA material model from [35,36] . The focus is on an effective solution procedure for local systems of equations when performing an update of the internal variables. In Section 3 , the 7-p shell finite element formulations are presented. In Section 4 , numerical examples illustrate the capability of the derived formulation, and the conclusions are drawn in Section 5 .

Continuous framework
In this section we revise the 3D constitutive model for a SMA that was originally presented by Souza et al. [50] for small strains and later extended to finite strains by Evangelista et al. [35] and Arghavani et al. [36] .

Model development
For a deformable body, the deformation gradient F , det F = J > 0 , can be decomposed with polar decomposition as F = R U = V R , where U and V are the symmetric positive definite right and left stretch tensors, respectively, and R is a proper orthogonal rotation tensor. The right Cauchy-Green deformation tensor and the Green-Lagrange strain tensor are: where 1 is the identity tensor. Similar to the finite-strain plasticity [67] , a multiplicative decomposition of F into elastic F e and transformation F t parts is postulated: and Eqs. (1) 1 , (2) , and (3) 1 can be used to relate C e with C as: It turns out to be useful to define the velocity gradient tensor L = ˙ F F −1 and the symmetric and skew-symmetric parts of L that are the rate of the deformation tensor d = 1 2 ( L + L T ) and the vorticity tensor w = 1 2 ( L − L T ) . With respect to the dot denoting the time derivative, it is worth noting that for time-independent material models, such as the one considered here, the notion "time" is used for the loading parameter that linearly increases/decreases the level of the external loading. The relation between d and the time derivative of the Green-Lagrange tensor ˙ E (i.e., the strain rate) can be obtained by using L as (see Eq. (A.1) in Appendix): According to the experimental results, the transformation is (almost) isochronic, which is expressed by det F t = 1 . Taking the time derivative gives tr d t = 0 (see Eq. The Helmholtz free energy ψ must depend on F e only through the elastic right Cauchy-Green deformation tensor C e , because of the material objectivity. It is assumed that ψ can be additively split into the elastic ψ e and transformation parts ψ t , and that the latter depends on the transformation Green-Lagrange deformation tensor E t = 1 2 ( C t − 1 ) (which plays the role of the internal variable) and temperature T : Moreover, let ψ be valid for both the austenite and martensite phases. For an isotropic SMA, the hyperelastic strain-energy function ψ e ( C e ) can be expressed by the invariants of C e as ψ e ( C e ) = ψ e (I C e , II C e , I I I C e ) . Let us choose the following strainenergy function of the neo-Hookean type: where μ and λ are Lams coefficients. Because II I C e = det C e = det ( C C −1 t ) , see Eq. (4) , and , the strain-energy function ( Eqs. (7) and (6) ) can be rewritten in terms of the tensors from the reference configuration The transformation part ψ t of the Helmholz free energy is chosen after [36] as: where τ M (T ) provides the temperature dependency of the material response h is the material-hardening parameter due to the transformation, β is a material parameter related to the dependence of the critical stress on the temperature, and T 0 = M f . Furthermore, the Macaulay bracket used here is defined as x = x −| x | 2 , and the norm in Eq. (9) is given as (11) handles the constraint on the transformation strain norm where ε L is another parameter of the material model linked to the maximum transformation strain norm at the end of the phase transformation. The indicator function states that the violation of inequality (12) is inadmissible. The Clausius-Duhem inequality of the second law of thermodynamics states the following: where S is the second Piola-Kirchhoff stress tensor and η is the entropy. By inserting Eq. (6) into Eq. (13) , we obtain: Let us elaborate on the expression labeled by ➀. The time derivative of ˙ C e can be obtained by using Eq. (4) and the relation S : For an elastic isothermal case, where D = 0 and there is no change in transformation, the expression for the stress tensor and entropy follow from Eq. (14) and Eq. (16) as: Because of Eq. (4) , we can conclude that It is assumed that the relations (17) are also valid in the case of a transformation that remains to be considered. To this end, let us introduce the back-stress tensor X = ∂ ψ t E t and examine the part of the expression (14) , labelled with ➁ (see Eq. (A.5) in Appendix): where N = E t || E t || and γ is a sub-differential of the indicator function ∂I (|| E t || ) : The back-stress tensor X depends on the temperature (for T > M f ) and the transformation strain tensor. The parameter γ controls the end of the phase transformation and ensures the elastic response of the martensite after a full transformation. By using Eqs. (16) , (17) and (19) in (14) , the inelastic dissipation for the case of transformation can be expressed as: The product X : ˙ E t in (21) can be further manipulated by using (5) and the relation tr ( A B ) = tr ( B A ) to obtain X : ˙ The decomposition L t = d t + w t and the symmetry of P result in P : w t = 0 (because w t is skew-symmetric), therefore the following holds P : L t = P : d t . These allow us to rewrite the dissipation inequality (21) as: By further introducing a transformation function f ( Z ) ≤ 0 (with its proposal left open for now) and assuming that the transformation case corresponds to the stress state giving f ( Z ) = 0 , the evolution equation for the transformation case can be obtained. We postulate that among all the admissible states of transformation, we choose the one that renders the maximum dissipation D t or the minimum of −D t . Recasting the original problem into the unconstrained minimization problem is possible by using the method of Lagrangian multipliers with the Lagrangian function defined as: where ˙ ζ is the Lagrange multiplier. In order to find the minimum of L ( Z , ˙ ζ ) , the stationarity, primal feasibility, dual feasibility and complementary slackness conditions must hold: From the stationarity condition (24) 1 , the evolution equation for the rate of transformation deformation d t is obtained: Let us choose the following expression for the transformation function (after [35,36] ), which is similar to the classic yield functions for metals: where Z D = ( P − K ) D is a deviatoric part of Z and the material parameter R defines the radius of the elastic domain and has a similar role as the yield stress in plasticity. By further using the relation for the symmetric tensor ∂ || A || A = A || A || , Eq. (25) yields the following evolution equation: which completes the list of the basic equations for the considered SMA material model.

Equations in the reference configuration
To have the equations of the SMA material model given with respect to the reference configuration, the evolution Eq. (27) must be recast in terms of C and E t = 1 2 ( C t − 1 ) (or yet C t ), which is the internal variable. To this end, we use the procedure from Eqs. (5) and (27) to obtain the equality: In order to facilitate the expression (28) , let us first investigate F T where the non-symmetric stress like tensor Y was introduced as: that finally leads to: Taking advantage of the symmetry of Y D C t , the norm || Therefore, in the view of the above, the evolution Eq. (28) for ˙ C t can be written as: from where the definition of the non-symmetric tensor A is obvious. The transformation function (26) and Kuhn-Tucker conditions (24) can also be rewritten in terms of Y as: The temperature effects in the material are linked with the back stress X and therefore with the tensor N , see Eq. (19) , which is the consequence of the choice for the transformation part of the free energy (9) . The temperature effects should also be taken into account in the elastic case; however, N is undefined for E t = 0 , which calls for an approximation of N for the elastic case. To this end, we can consider Eq. (34) during nucleation (i.e., at beginning of the forward phase transformation), which yields (assuming X = 0 ) which indicates that during nucleation E t evolves in the direction of ( C S el ) D . Here, S el = S ( C t = 1 ) . Following [36] , we adopt the direction (37) for all the elastic states. Thus, N is redefined for computational purposes as:

Discrete framework
In the finite-element inelasticity, the solution in time is advanced through a sequence of discrete time points, t 0 , ..., t n , t n +1 , ..., t end . Because of the single-step integration method, it suffices to focus on the generic increment from t n to t n +1 . In our case, the central problem is: given the values of the state variables at t n , i.e., C (t n ) = C n and C t (t n ) = C t,n , find their (converged) values at t n +1 , i.e. C n +1 and C t,n +1 , by taking into account the constraint (12) . To that end, the implicit methods to solve the evolution equation for C t,n +1 and the equilibrium equations for C n +1 are combined with the Newton-Raphson iterative method. The operator-split method is used to simplify the computation of the central problem. The update of C t,n +1 for a given value of C (i ) n +1 is followed by the update of the latter for the given value of the former. The update of C t,n +1 takes place at Gauss integration points, while the update of C (i ) n +1 requires a consistent tangent operator and an assembly of the finite-element contributions into a global, non-linear system of equations, with the nodal generalized displacements as the unknowns. The index (i ) denotes an iteration related to the global system.

Computations at the Gauss point
To integrate the evolution Eq. (34) , we use the exponential mapping approximation, see, e.g., [68] , which algorithmically conserves the inelastic volume (i.e., preserves det C t = det F t = 1 ) as required by the SMA model. For a given system of the first-order differential equations ˙ y = G ( y ) y , the exponential approximation leads to y n +1 = exp ((t n +1 − t n ) G n +1 ) y n . By applying this formula to Eq. (34) , we obtain: We compute the asymmetric matrix exponent in Eq. (40) as proposed in [69,70] . A Gauss-point update starts with the evaluation of the trial state, which is obtained by keeping the transformation flow frozen, i.e., C tr t,n +1 = C t,n . With these values and given C (i ) n +1 , the trial stresses S tr n +1 and Y tr n +1 are computed. The trial state is a natural way to check the loading/unloading conditions. The elastic loading, unloading (or neutral loading) are simply verified if the nucleation/completion criterion (given below) during the trial state is satisfied, which enforces ζ tr n +1 = 0 . In that case, the step is elastic. In the opposite case, inelastic loading is the result, and we need to compute the updated values of C t,n +1 and the parameters ζ n +1 and γ n +1 , which will enforce the consistency. The computation of the admissible state variables is given in terms of a set of non-linear equations with C t,n +1 , ζ n +1 and γ n +1 as the unknowns.
For the inelastic step, there are two possibilities: the unsaturated state where || E t || < ε L , and the saturated state where || E t || = ε L . As a consequence, two sets of equations are generated. For the unsaturated state, the following system (labelled hereinafter as PT1) needs to be solved, where the last equation is simply the enforcement of the parameter γ n +1 to 0. For the saturated state, Eqs. (41) and (42) This set of equations is labelled hereinafter as PT2. Let us write PT1 and PT2 in a more compact form as: and the vector of unknowns as The local Newton-Raphson method (with (k ) denoting an iteration) is used to solve PT1 and PT2 as: with the update and repeated iterations until the convergence || h (k ) n +1 || < tol is reached, with tol = 10 −8 for the numerical examples presented below.
The PT1 system (41) - (43) handles the transitions between the elastic and transformation steps, which are nucleation (i.e., || E t,n || = 0 and || E t,n +1 || > 0 ) and completion (i.e., || E t,n || > 0 and || E t,n +1 || = 0 ). For these steps, the tensor N is defined as given in Eq. (38) , in accordance with the proposal in [36] . Therefore, N n +1 and Y D n +1 are computed as: The trial value f tr of the transformation function (35) , which is computed with Eq. (50) , the primal feasibility Kuhn-Tucker condition, and k τ > 0 give the following nucleation condition (NC): When the condition (51) is violated, the transition from the elastic to the transformation state takes place in the current step. Similarly, with k τ < 0 , the completion condition (CC) is obtained: When the condition (52) is conformed, the transition from the transformation to the elastic state takes place. We note that the left-hand side of the final equations in (51) and (52) is always positive, as is the right-hand side of the NC in the inequality (51) . However, the right-hand side of the CC in Eq. (52) can be negative for τ < R , which, however, only occurs when not dealing with a superelastic regime.

Solving nonlinear systems at the Gauss point
It is suggested in [36] that PT1 is always solved first for an inelastic step. If the PT1 solution, labelled as h PT1 n +1 , is inadmissible because of the violation of condition (12) , it is suggested to proceed by PT2 and compute h PT2 n +1 with h PT1 n +1 as the initial guess. It is our experience, however, that this procedure does not always work well.
Our numerical experiments show that much better convergence properties are obtained for the inelastic step using the algorithm presented in Fig. 3 . It has four options: (i) For γ n = 0 , PT1 is solved first (with h n as the initial guess). (ii) If the PT1 solution violates the condition (12) , i.e., if || E t,n +1 || > ε L is computed, it proceeds with PT2 and takes h PT1 n +1 as the initial guess. (iii) For γ n > 0 , PT2 is solved first, however, the 7th component of the initial guess h n , i.e., ζ n , is multiplied by a factor k . Our numerical tests suggest that k = 2 is an effective choice. (iv) If the PT2 solution violates the condition (20) , i.e., if it computes γ n +1 < 0 , we proceed with PT1 and the modified initial guess h n that enforces γ n = 0 .
Implementation of the line-search method, which prevents too large increments h (k ) in Eq. (47) , also improves considerably the convergence properties. The line search introduces the parameter α LS that controls the size of the increment as: where 0 < α LS ≤ 1 . Usually, see [41,71] , α (1) LS = 1 at the beginning of each iteration. If the inequality is computed, α LS is reduced by a factor 0 < β < 1 (usually β = 0 . 5 ) as α ( j+1) The loop iterates as long as the condition (54) holds. By lowering the norm of Q n +1 , the convergence is improved in each iteration of the local Newton-Raphson loop, which, however, might require a large number of iterations.
For the saturated state, the largest contribution to || h (k ) n +1 || is the change of the parameter γ . Therefore, it seems reasonable at first glance to define α LS with respect to the increment γ . However, γ → ∞ during loading (and γ → 0 at unloading) makes the limitation of || h (k ) n +1 || troublesome. A better alternative is to define a limit value for the Lagrange multiplier ζ > 0 and relate it with α LS as The choice (55) does not guarantee descending in each local Newton-Raphson iteration because the inequality (54)

Seven-parameter shell model and finite-element formulations
We briefly present a 7-p finite-rotation-shell model. For more details we refer to, e.g., [62,[72][73][74] . The extensible-director kinematic hypothesis is assumed, along with the quadratic through-the-thickness variation of displacements. We adopt the standard notation for the indices: small Greek letters for indices from 1 , 2 and small Latin letters for indices from 1 , 2 , 3 .

7-p shell model
The shell is embedded in 3D space with a fixed orthonormal basis { e i } . It is described as a surface with an extensibledirector field, with the position vector to the material point in the initial configuration given as: where ξ i are the convected curvilinear coordinates. Hereinafter, we will omit writing the functions' and functionals' arguments for brevity. Thus, ϕ 0 describes the shell mid-surface M , D is the normal-to-the-M vector field of unit length that is called the shell director, h is the initial shell thickness, and A is the domain for the mid-surface parametrization. The position vector to the material point in the deformed (current) configuration is assumed to be: where u is the mid-surface displacement vector field, λ is the scalar through-the-thickness-stretching field, and q is the hierarchical scalar field related with the quadratic variation of displacements in the through-the-thickness direction. Moreover, a is the rotated shell director field in the deformed configuration, see, e.g., [75][76][77] for details. With (56) and (57) , the covariant base vectors for the initial and deformed configurations are where d ,α = λ ,α a + λa ,α and f ,α = q ,α a + q a ,α . The dual base vectors G i are defined through the orthogonal condition G i · G j = δ i j , where δ i j is the Kronecker delta symbol. The Green-Lagrange strain tensor for the considered shell model is defined as: with the covariant components of the strain tensor given as functions of ξ 3 , see [73] : Following [73] , we truncate (60) after the quadratic term, thus introducing the reduction The explicit forms of the covariant components of the reduced strain tensor E = E i j G i G j are: We note that in the computer code the covariant components (62) -(64) are transformed to the Cartesian components using the local Cartesian basis defined for each integration point. The weak form of the equilibrium equations at time t n +1 is (index n + 1 is omitted hereinafter for brevity): is the variation of the reduced strain tensor.
After the spatial discretization of the mid-surface M with n el non-overlapping finite elements, such that M ≈ n el e =1 A e , and after using the appropriate interpolation for generalized displacements and their admissible variations, the functional (65) becomes an assembly of the finite-element contributions, with the mesh-related nodal values for generalized displacements and the Gauss-point-related values for the transformation tensor C t as the unknowns: Here, A denotes the finite-element-assembly operator, e.g., [67] , and G e int and G e ext denote the element's approximation of the first and second terms on the right-hand side of Eq. (65) , respectively: The current values of the stress tensor are obtained using the algorithm presented in Section 2.2 , which is one of the two steps of the operator-split method. For the other step, where the solution A n el e =1 e is computed for a given loading, the linearization of (67) is needed, which can be written as Lin Here, C = ∂ S / ∂ E is the fourth-order tensor, called the consistent inelastic tangent operator, which is computed already in the first step, and

Three 7-p quadrilateral finite-element formulations
We describe the assumed natural strain (ANS) and the enhanced assumed strain (EAS) interpolations that are applied for three 7-p quadrilateral finite elements used to compute the numerical examples.

Basic interpolations
The initial mid-surface and the initial shell director over the quadrilateral 7-p shell finite element are interpolated as The deformed configuration of the element is approximated with: Note that the interpolations for the shell director (71) and (72) demand applying the rotations and rotational parameters only at the nodes of the finite element.

ANS enhancements
To avoid the transverse shear, the membrane, and the through-the-thickness lockings, and to reduce the sensitivity of the quadrilateral shell element to geometric distortion, the ANS interpolations are used.
For the membrane covariant components of the strain tensor, we adopt the ANS interpolation proposed by Kulikov and Plotnikova [63] for the five-parameter (5-p) quadrilateral shell element: where A, B, C, D and E define the points where the strains are evaluated with Eqs. (61) -(64) , see Fig. 4 . The a factors that appear in the above expressions measure the mesh distortion, e.g. [63][64][65][66] . For our element called 3ANS, they were computed for the last converged configuration [65] , and for our element called EAS, they were computed for the initial configuration.
For the transverse shear covariant components of the strain tensor, we adopt the ANS interpolation proposed by Dvorkin and Bathe [61] for the quadrilateral based on the 5-p shell model: Equations (61)

EAS enhancements
The performance of the element can also be improved by application of the EAS concept, which relies on the enhancement of the covariant components of the strain tensor by adding the enhancing components. The in-plane covariant components of the strain tensor, H αβ , K αβ and L αβ , can be additively enhanced by covariant components of the enhancing strain tensor denoted as ˜ H αβ , ˜ K αβ and ˜ L αβ . The latter can be collected in vectors as and interpolated over the element as where the interpolation was chosen based on experience with 5-p shell quadrilaterals, e.g., [65,66,78] , as In (80) , α H αβ , α K αβ and α L αβ are vectors of the additional EAS parameters of the element, with each vector having five parameters. The EAS enhancement of the in-plane covariant strains can therefore be written as E αβ = H αβ + ξ 3 K αβ + ξ 3 2 L αβ . Also, the through-the-thickness normal covariant components of the strain tensor, H 33 and K 33 , can be additively enhanced by the EAS strains. We use the following interpolations where TH = ξ η ξη , (83) and the vectors α H 33 and α K 33 consist of three additional EAS parameters each. We note that in the computations, the components of the strain and stress tensors with respect to the local Cartesian basis are used in the Gauss points. The transformation of the covariant generalized-displacement-based, ANS, and EAS strain components is performed in the standard manner, see, e.g., [66] . The variational framework for the EAS enhancements is a modified Hu-Washizu functional, which is without the stress tensor (e.g., [78,79] ), and with the strain tensor defined as E EAS = E + E . Its stationary point provides two equations: among which (85) is used for condensation of the additional EAS parameters at the element level. As shown for the 5-p shell model quadrilateral [65,80,81] , combining the ANS and EAS concepts as E EAS * = E ANS + E leads to very efficient formulations. In that case, (84) and (85) The interpolations for the three finite elements applied in the section with numerical examples are listed in Table 1 , where "Disp." stands for standard, generalized-displacement-based interpolation, and the number in parenthesis denotes the number of EAS parameters in the EAS formulations.

Numerical examples
The computer codes for the above-described material and shell models were produced using the finite-element code generator AceGen, see [69,70] . They were implemented into the finite-element computer code AceFEM [82] that was used to compute the numerical examples given below. Three sets of material parameters are used: Set 1, see (88) , is taken from Arghavani et al. [36] , Set 2, see (89) , is taken from Evangelista et al. [35] , and Set 3, see (90) , is calibrated to match the experiment in [83] . Set 1: Set 2: To obtain the superelastic response, the ambient temperature T is set to 37 • C for Set 1, to 11.85 • C for Set 2, and to 27.5 • C for Set 3. The numerical examples were computed with all three formulations (2ANS, 3ANS and EAS) using a 2 × 2 × 3 Gaussian quadrature, i.e., with three integration points in the thickness direction. For examples with identical results for all formulations, the label "2ANS" is used in the figures below. Where the results are identical for the 2ANS and 3ANS, but not for the EAS, the labels "2ANS" and "EAS" are used. The components of the second Piola-Kirchhoff stress tensor and the components of the Green-Lagrange strain tensor are shown when presenting the response at the integration point.

Tension of a square element
A square element with an edge length L = 10 mm and a thickness h = 0.01 mm is subjected to tension. The supports are applied in such a way that the homogeneous stress state is obtained (see Fig. 5 (left)). The edge 3 is subjected to an imposed displacement u x = u max λ, where u max = 1 mm and λ is a load multiplier. Set 1 of the material parameters is used. Loading and unloading with a total of 320 steps is considered, with λ increasing from 0 to 1 and decreasing back to 0.
The stress-strain relation at the integration point is shown in Fig. 5 (right). Comparing the curve with that from Arghavani et al. [36] (labelled as "Ref."), only a small discrepancy is observed for large stresses due to the use of different strainenergy functions. Namely, we use the Neo-Hokean strain-energy function, while in [36] the Saint-Venant Kirchoff model was applied. Moreover, the reference curve from Arghavani et al. [36] was obtained by directly solving local equations that appear at the Gauss point for the plane-stress case, while our results are obtained by applying a load to a single finite element. Characteristic points are marked on the curve to show how the variables ζ and γ are changing during loading and unloading. In the beginning, the response is elastic, therefore ζ = γ = 0 . During the forward martensitic transformation, the PT1 system is being solved and ζ is increasing. When the transformation is completed, the saturated state is present, and therefore the PT2 system gives an admissible converged solution with the constant ζ and an increasing γ . The higher the load, the greater the value for γ . At the beginning of the unloading, the first few steps are elastic, because the Kuhn-Tucker condition for the limit function is not violated ( f tr < 0 ). With further unloading, ζ remains constant and γ is decreasing towards 0. At the beginning of the reverse transformation, γ = 0 and the unsaturated state is present. During the reverse transformation, ζ is increasing, while γ remains at 0 and the PT1 system is being solved. After completion of the reverse transformation, the elastic response is obtained, i.e., ζ and γ are constant.

Shear of a square element
A square element with the edge length L = 10 mm and thickness h = 0.01 mm is subjected to a shear load, see Fig. 6 (left). The displacement u x = λ u max is imposed on edge 2, where u max = 1.5 mm and λ is the load multiplier. Loading and unloading in a total of 80 steps is considered by increasing λ from 0 to 1 and decreasing it back to 0. The Set 1 of the material parameters is used. With the boundary conditions shown in Fig. 6 (left), a homogeneous stress state is obtained. We note that for a mesh of several elements, this is no more the case, because this kind of loading does not produce a pure shear state in a nonlinear neo-Hookean solid. Nevertheless, we use a single-element test in order to compare our results with the pure shear test reported in [36] . We can see from Fig. 6 (right) that our curve matches perfectly with the curve "Ref." from Arghavani et al. [36] . It is worth mentioning that the reference curve is obtained by directly solving the local equations that appear at the Gauss point for the pure-shear case, while our results are obtained by applying a load on a single finite element.

Square wall in tension and compression
A square wall with the edge length L = 10 mm and thickness h = 0.01 mm is subjected to non-proportional loading. The wall is supported along the edges 1 and 2 in such a way that a homogeneous stress state is obtained (see Fig. 7 (left)). A structured mesh of 5 × 5 elements is used, along with the Set 1 of material parameters.
The loads are q 1 = λ 1 q 1 ,max and q 2 = λ 2 q 2 ,max , where λ 1 and λ 2 are the load multipliers, and q 1 ,max = 7 . 524 N/mm and q 2 ,max = 7 . 559 N/mm. The values of the load multipliers are adjusted load-step-wise, as shown in Fig. 7 (right), so that the butterfly-like stress pattern from Fig. 8 (left) is obtained, with the maximum stresses S 11 = S 22 = 700 MPa. This allows us to compare our results with those from Arghavani et al. [36] in Fig. 8 (right), where the normal in-plane strains are shown. Our results fit well with the reference results, which were obtained with the Gauss-point algorithm for the butterfly-like stress control shown in Fig. 8 (left).

Beam in tension
A straight beam of length L = 100 mm, width w = 20 mm, and thickness h = 1 mm is subjected on edge 2, see Fig. 9 , to a tensile load q 1 = q λ, where q = 400 N/mm and λ is the load multiplier (after [35] ). The mesh is 10 × 1 elements, and the Set 2 material parameters are used. The loading and unloading is considered, with λ increasing from 0 to 1 and decreasing back to 0. Four cases of boundary conditions are studied. The first three cases have the edge 1 supported as u x = u y = u z = 0 .   Additionally, for the first case (C1), the thickness change (i.e., the 6th and the 7th nodal degrees of freedom) is restricted to zero over the whole mesh. For the second case (C2), the thickness change is restricted to zero only on edge 1, while for the third case (C3), no additional boundary conditions are applied. For the fourth case (C4), the point with coordinates (0,0,0) is restrained as u x = u y = u z = 0 , and the rest of the edge 1 is restrained as u x = u z = 0 (giving a homogeneous stress state). Thus, C1 is the plane-strain case, and C3 and C4 are the plane-stress cases.
The results show that for the plane-strain case C1, E 33 = 0 for 2ANS and 3ANS, and E 33 ≈ 10 −3 for EAS. For the planestress cases (C3 and C4), S 33 ≈ 0 for EAS, which, however, is not the case for 2ANS and 3ANS. For C2, the influence of the edge support is much larger for 2ANS and 3ANS than for EAS. Figure 10 shows the axial force at edge 1 versus the axial displacement of the point with initial coordinates (L, w, 0) . Interestingly, the difference between the computed stresses and strains for the three formulations reflects only slightly in the load-displacement curves for C1, C3 and C4, and it is considerable only for C2. The curves are compared with those of Evangelista et al. [35] , who used a mesh of 10 × 1 2D elements that were claimed as plane stress, and C3 boundary conditions. It can be observed that none of our curves matches with that from Evangelista et al. [35] . It seems that the results from Evangelista et al. [35] are neither plane stress nor plane strain, but, unfortunately, a description of the applied 2D element and the corresponding SMA formulation are not provided in [35] . Part of the difference between our results and those of [35] can be attributed to our use of the nucleation and completion conditions, while [35] used a regularized norm for || E t || , which means that the material is treated as inelastic from the start of the analysis. This makes a sharp transition at the beginning of the transformation (when loading) and at the end of transformation (when unloading) in our formulations. Figure 11 shows distribution of non-zero stresses through the thickness of EAS element: at the near-support-point marked with P in Fig. 9 , for C3, at λ = 1, and for 3, 5 and 7 through-the-thickness integration points.
Small values of S 33 (in the range of 10 −3 S 11 ) and S 13 (in the range of 10 −2 S 11 ) are observed, with integral of S 13 over the thickness yielding zero. It seems that our 7-p formulations exhibit slight 3D effects (expressed by non-zero S 33 and S 13 Fig. 11. Through the thickness stress distribution for C3 and EAS element, at point P. in this particular case) in regions where localized 3D effects may be expected. It also seems that such 3D effects are more pronounced when (volume preserving) transition takes place. However, this has a negligible influence on the global (stress resultant) level of the response. For example, C3 and C4 load-displacement curves in Fig. 10 are practically identical, and for C4, the only non-zero stress is S 11 .

Curved beam
A curved beam with thickness h = 1 mm (see Fig. 12 ) is analyzed (after [35] ) for four cases of boundary conditions (C1-C4), wherein the first three cases are the same as those in the previous example, while for C4 the point with coordinates ( R o , 0,0) is restrained as u x = u y = u z = 0 , and the rest of the edge 1 is restrained as u y = u z = 0 . On edge 2, the line load q 1 = q λ is applied, where q = 60 N/mm and λ is the load multiplier, which is increased from 0 to 1 and decreased back to 0. Set 2 of the material parameters and a regular mesh with 40 elements in the circumferential and 10 elements in the radial direction are used.
In Fig. 12 , the amount of transformation in the middle surface is shown. Figure 13 shows the relation between the applied force and the displacement in the −x direction of the point with the initial coordinates (−R o , 0 , 0) . The curves are compared with those from reference [35] , where a mesh of 40 × 10 2D elements was used for a small-strain SMA, labelled as "Ref. SS", and for a large-strain SMA, labelled as "Ref. FS". As for the previous example, none of our curves matches the one from Evangelista et al. [35] , indicating again that the 2D element in [35] is neither plane stress nor plane strain. The curves C2, C3 and C4 are similar, while the plane-strain curve C1 shows the stiffest response. For this example, there is no significant difference between C2 and C3 (and C4), in contrast with the previous example. In general, the EAS formulation is softer than the 2ANS and 3ANS. It is worth noting that for the 2D tests presented in [36] (i.e., the first three examples), the results of our formulations match perfectly with those from Arghavani et al. [36] . On the other hand, for the 2D tests presented in [35] (i.e., this one   and the previous one), the results of our formulations do not match with those from Evangelista et al. [35] . We attribute the difference to the unknown 2D formulation in [35] . As shown in Section 4.8 for a 3D example, our results match perfectly with those from Evangelista et al. [35] .

Cook's membrane
A Cook's membrane of length L = 48 mm and thickness h = 1 mm is considered (see Fig. 14 (left)). Edge 1 is fully clamped, i.e., all 7 degrees of freedom are set to 0. Load q 1 = q λ is applied on edge 2, where q = 625 N/mm and λ is the load multiplier that goes from 0 to 1 and back to 0. A regular mesh of 16 × 16 elements and Set 1 of the material parameters are used. Figure 14 (right) presents the relationship between the vertical reaction force on edge 1 and the displacement in the −y direction for the middle point on edge 2. Figure 15 shows the initial and deformed meshes, along with the transformation strain norm || E t || for the middle surface of 2ANS. The stress state through the thickness is not homogeneous due to the imposed boundary conditions on edge 1. Small differences between the three formulations are observable, and like previous examples, the EAS response is slightly softer.

Bending of a cylindrical panel
A cylindrical panel with radius R = 5 mm, angle θ = π / 4 , width w = 2 mm, and thickness h = 0.5 mm is subjected to a bending moment m 1 (see Fig. 16 ). Edge 1 is fully clamped, i.e., all 7 degrees of freedom are set to 0. Load m 1 = m λ is applied on edge 2, where m = 300 Nmm/mm and λ is the load multiplier. A mesh with 32 elements in the bending direction and 16 elements in width direction is used. Loading and unloading is considered, where λ is changed from 0 to 1 and decreased back to 0. Set 1 of the material parameters is used. This is a pure bending test, while all the previous examples were pure membrane tests. Figure 17 (left) shows the superelastic response of the cylindrical panel by presenting the relationship between rotation (around −y axis) of the middle point on edge 2 and the reaction moment on edge 1. In order to show the ability of our algorithm to find solutions for fully transformed material, the value of the applied load is high. A solution is only possible to obtain by using the line-search method, when solving local equations at the Gauss points (see Section 2.2.2 ), with ζ limit = 0 . 005 . Without applying the line search, the algorithm finds converged solutions only for the applied moments that are less than 200 Nmm. The superelastic response is shown in Fig. 17 (right), where displacements in −x and −z of the middle point on edge 2 are shown. The initial and deformed meshes for λ = 0 . 2 and λ = 1 are shown in Fig. 18 , along with the amount of transformation. The "upper" and "lower" surfaces are located at ±0 . 5 3 / 5 h from the middle surface. The transformation at the "upper" and "lower" surfaces is not completely symmetric with respect to the middle surface due to slight curvature deformation in the shorter direction.

Helical spring
A helical spring (see Fig. 19 (left)) of inner radius R 1 = 2 . 5 mm, outer radius R 2 = 3 . 5 mm, step s = 1 . 2566 mm, and thickness h = 0 . 5 mm is considered after [35] . Edge 1 is fully clamped, i.e., all degrees of freedom are 0, while edge 2 is subjected to load q 1 = q λ, where q = 4 . 14 N/mm and λ is the load multiplier. A regular mesh with 25 elements in the circumferential and 6 elements in the radial direction is used. Following [35] , the mechanical and temperature loading are applied in three steps, where the duration of each step is 1 unit of dimensionless pseudo-time (see Fig. 19 (right)). The loading procedure is defined in the following way: (i) increase of λ from 0 to 1 at a constant ambient temperature of 11.85 • C, (ii) an increase of the temperature from 11.85 • C to 526.85 • C at λ = 1 , and (iii) a decrease of the temperature back to 11.85 • C at λ = 1 . Set 2 of the material parameters is used.
The relation between the average stress on edge 1, denoted as Q , and the displacement in the −z direction of the middle point on edge 2 is shown in Fig. 20 (left). The average stress is obtained by dividing the vertical reaction force by the cross-section area (R 2 − R 1 ) h = 0 . 5 mm 2 . Our results are compared with those from Evangelista et al. [35] , where a mesh of 25 × 5 × 5 8-node solid elements was used. Despite the small discrepancies, which might be because of the differences between the solid and shell formulations, meshes, and hyperelastic functions, our curves match well with the reference curves. For 2ANS, the maximum displacement at the end of the first loading step is almost identical with the reference curve, while at the end of the second step the difference in the displacements is 6.4%. The EAS formulation is softer than the 2ANS; therefore, the maximum displacement at the end of the first loading step is higher; however, the displacement change during the second loading step is similar to that of the reference curve. Fig. 20 (right) shows the displacement of the middle point on edge 2 during the loading procedure. At the end of the third step, when the temperature is back to the initial value, the displacement caused by the temperature change is also recovered. This is also seen in Fig. 20 (right), where the initial and deformed shapes at the end of each step for the 2ANS formulation are shown. Figure 21 presents distribution of stresses across the helical spring at the end of step (i) for EAS element. Stresses at the levels of the three through-the-thickness integration points are shown.
The order of S 33 is rather high, roughly 1 / 10 of the order of other stresses. This is because the shell is not very thin, h/R 2 = 1 / 7 , and due to slight 3D effects in more stressed part of the spring.

Half sphere
A half sphere of radius r = 10 mm, with a hole on top, angle θ = 2 π / 5 , and thickness h = 0.4 mm, is subjected to 2 point forces, see Fig. 22 .
This is a standard benchmark test for shells, because of its response that couples the membrane and bending effects. Double symmetry is considered; therefore, only a quarter of the half sphere is analysed by taking into account the appropriate symmetry boundary conditions. In order to prevent rigid-body motion, the displacement of point P 0 with the coor-  dinates (r sin (π / 10) , 0 , r cos (π / 10)) is restrained in the z direction. At points P 1 = (r, 0 , 0) and P 2 = (0 , r, 0) , forces F x = F λ and F y = F λ are applied, respectively, where F = 100 N and λ is the load multiplier.
Four cases with fine and coarse meshes are considered. The first and second cases (C1 and C2) have regular and distorted meshes with 4 × 4 elements, respectively. The third and fourth cases (C3 and C4) utilize regular and distorted meshes with 16 × 16 elements, respectively. The distorted meshes have the ratio L 1 /L 2 = 16 . Loading and unloading is considered by increasing λ from 0 to 1 and back to 0, and Set 1 of the material parameters is used. Figure 23 shows the relationship between the force and the displacement of points P 1 and P 2 in the x and −y directions, respectively. For finer meshes, the results for 2ANS and 3ANS and for both meshes are almost identical; therefore, only one curve is shown to represent all four analyses. Similarly, the difference between the regular and distorted finer meshes for the EAS are negligible. For the coarse meshes, the 3ANS results are much closer to the fine-mesh results than the 2ANS results, while the EAS results are the closest to the fine-mesh result. For the coarse regular (R) mesh, the difference between the 2ANS, 3ANS and EAS results is large. For the coarse distorted (D) mesh, the difference in the results between the 2AND and 3ANS is smaller. However, the EAS results are different. This is because of the way the distortion is applied (the mesh is finer in the vicinity of the forces). This example shows the advantage of 3ANS over 2ANS for shell problems where the interaction between the bending and membrane effects is present. The 3ANS and EAS give much better results for coarse and distorted meshes. In Table 2 , u x of point P 1 and u y of point P 2 are collected for a particular load multiplier for all elements and meshes. Figure 24 shows the initial and deformed fine regular mesh at λ = 1 . The amount of transformation is shown as well, with the "upper" and "lower" surfaces defined in the same way, with respect to the thickness, as in the Example 4.7 above. As expected, the transformation is not homogeneous through the thickness and a larger level of transformation is present in areas near the application of the forces.

Twisted beam
Similar to example in [84] a beam of length L = 12 mm and for 2 π twisted along its axis, with width w = 2 mm and thickness h = 0.2 mm (see Fig. 25 (lower)) is analyzed. It has one edge clamped (i.e., all 7 degrees of freedom are 0 in this

Table 2
Displacements u x (point P 1 ) and u y (point P 2 ), in [ mm ] , for a particular value of load multiplier. case). At the opposite edge, the axial displacement u = u max λ is imposed, where u max = 1.1 mm and λ is load multiplier. Set 1 of material parameters is used for regular and distorted meshes, with 32 elements in longitudinal and 16 elements in the transverse direction. For the latter mesh, the ratio r = L 1 /L 2 = 1 . 1 is applied (see Fig. 25 (upper)). For the regular mesh, the loading and unloading are considered through the displacement control, with increasing λ from 0 to 1 and decreasing it back to 0 in order to obtain a superelastic response. However, for the distorted mesh, the symmetry with respect to the beam axis is lost, which induces small geometric imperfections that trigger instability, and the arc length method [85] is  used, which does not allow the control of loading and unloading by prescribing λ, because the latter becomes unknown that is computed in the solution process. Figure 26 shows relationship between the reaction force on edge 1 and the displacement of the mid-point on edge 2, both in the −x direction. It can be seen that for the regular mesh, the loading and unloading generate superelastic responses without any instability. A deformed regular mesh is shown in Fig. 25 (lower) for λ = 0 . 65 , with colours indicating the amount of transformation in the mid-surface (we can see the symmetry of the transformation with respect to the midaxis). This is in contrast with the results for the distorted mesh, when buckling occurs, for the 2ANS and 3ANS elements at approximately 0.5 mm and for the EAS element at approximately 0.4 mm. The arc-length starts with unloading at some point, but the algorithm is not able to find much of the subsequent solution path. Figure 25 (upper) shows the deformed distorted mesh at the maximum displacement. This example illustrates how small imperfections (induced here simply with a mesh distortion) generate very different results for a buckling-sensitive SMA.

Tube in tension
A cylinder of length L = 30 mm, mid-diameter d = 4.326 mm and thickness h = 0.348 mm is subjected to axial tension, see Fig. 27 (left). On both edges, all degrees of freedom are set to 0, except for axial displacement at edge 2, which is u = u max λ, where u max = 1 . 1 mm and λ is load multiplier. Loading and unloading is considered by increasing λ from 0 to 1 and back to 0. Regular mesh with 40 × 40 elements is used. The main purpose of this example is identification of material parameters by comparing our results with experimental data from Helm and Haupt [83] , and use these parameters to compute the next example. The calibration yielded Set 3 data given in (90) for the ambient temperature T = 27 . 5 • C , which was roughly the average temperature during the experiment. Figure 27 (right) shows our results (with Set 3 data), experimental results [83] and numerical solution [41] with 3D solid finite element. Axial stress is σ = F /A , where F is axial reaction force and A is initial cross section area. Axial strain is ε = l /l , where l is displacement difference between the nodes initially at z = 15 mm ± 0 . 75 mm and l is the initial distance between these nodes (i.e. 1.5 mm). Some discrepancies are present between our curves and the experiment curve ("Ref. Exp"), while differences between our curves and the curve from Reese and Christ [41] ("Ref. Sim") are minor. At the beginning of the martensite transformation and at the end of the reverse transformation our curves are non-smooth because of nucleation and completion condition.

Tube in tension and torsion
The cylinder from the previous example is subjected to tension and torsion, after [83] and [41] , see Fig. 28 . Edge 1 is clamped (i.e. all degrees of freedom are 0), while edge 2 is subjected to axial stretching u z = u max λ 1 followed by torsion φ z = φ max λ 2 , where λ 1 and λ 2 are load multipliers. Torsional rotation is applied by imposing u x and u y due to φ z . The thickness change and the director rotations (i.e. 4th to 7th degrees of freedom) are set to 0 on edge 2. The u max = 0 . 903 and φ max = 0 . 722 are such that 3% normal and shear deformations are reached, with normal deformation defined as ε = L/L (where L is axial displacement on edge 2) and shear deformation defined as γ / √ 3 , where γ = (φ z d) / (2 L ) . Loading consists of 4 steps: (i) λ 1 is increasing from 0 to 1 at λ 2 = 0 , (ii) λ 2 is increasing from 0 to 1 at λ 1 = 1 , (iii) λ 1 is decreasing from 1 to 0 at λ 2 = 1 , and (iv) λ 2 is decreasing from 1 to 0 at λ 1 = 0 . Regular mesh with 40 × 40 elements and Set 3 material parameters are used. Figure 29 (left) compares our results with experimental results [83] ("Ref. Exp") and simulation from Reese and Christ [41] ("Ref. Sim"), in which 3D solid elements were used. In Fig. 29 (left), σ = F /A , where F is axial reaction force and A is initial cross section area, and τ = (Md/ 2) / (π / 32 (D 4 in )) , where M is torsional moment and D in and D out are the inner and the outer diameter of the tube, respectively. The match of our curve to the experiment and the reference simulation is acceptable, although not perfect. The main reason for the differences between the simulations (2ANS and "Ref. Sim") are different constitutive equations used. However, both simulations lack the compression-tension asymmetry, which may be the reason for the deviations from the experiment. Another significant difference between the experiment and simulations is temperature change of the material. In [83] , temperature changes are measured during the experiments, indicating non-isothermal conditions, whereas in simulations isothermal loading is considered. Figure 29 (right) shows distribution of the transformation, which is similar as in [41] . Figure 30 compares normal and shear stress-strain curves with experiment [83] . Normal stress-strain curve matches well the experiment, especially in the first and the second loading step, while the  matching for the shear stress-strain curve is worse. We assume that the main reason is temperature change of the sample as a consequence of transformation during the experiment. Namely, the temperature increase at loading steepens the slope of the curve during transformation.

Conclusions
In this work we have derived, for the first time, a shell finite element formulation for the simulation of deformations of thin-walled shape-memory alloys under mechanical and thermal loadings. Such formulations are much more suited to thinwalled and curved SMA structures than solid finite element formulations, which can have difficulties correctly representing the bending and buckling of the shell-like structures. As for the SMA material model, we revised the finite strain model presented in [35] and [36] . Novel approaches related to its efficient finite element implementation are described. As for the shell finite element formulations, we have applied a seven-parameter shell model and three corresponding finite elements that extensively use an assumed natural strain and enhanced assumed strain enhancements. Although only one formulation can fulfill the condition of the zero-through-the-thickness normal stress, and the other two can do that only approximately, the corresponding load-displacement curves differ only slightly. Also, although the formulation with the enhanced assumed strain enhancements can approximate better the incompressibility of the transformation strains, its results do not differ much from the formulations with the assumed natural strain formulations.
Numerical examples show that for an excellent match for the results obtained with the solid SMA formulation, the number of elements needed for the shell SMA model is less than 25% of the number of elements needed for the solid SMA model. We believe that for thin-walled SMAs the same (or even superior) accuracy of the results can be obtained with shell formulations in comparison with solid formulations, for a considerably reduced number of elements in the mesh, and for a significant reduction in the computational time. Moreover, two of the applied shell finite element formulations, i.e. 3ANS and EAS, are designed to be low-sensitive with respect to mesh distortion, which can further reduce the necessary number of elements in the mesh.
Although the numerical simulations offer the rapid generation of results, great attention is required for the determination of the input-material parameters to obtain reliable results. The process of determination for the SMA material parameters from experiments is trivial for small strain formulations, whereas for finite strain formulations more caution is needed, since the finite strain material models include a differentiation between the (reference, intermediate and current) configurations. Therefore, our future work will involve finite strain SMA material parameter determinations as part of developing opportunities for an enhancement of the current SMA shell formulation, e.g., a smooth transition of the thermomechanical response at the beginning and the end of the phase transformation and a consideration of the compression-tension asymmetry.