M Peigney

G Scalet

F Auricchio

A time-integration algorithm for a 3D constitutive model for SMAs including permanent inelasticity and degradation effects

Keywords: Shape-memory alloys, incremental energy minimization, permanent inelasticity

Components based on shape-memory alloys are often subjected to several loading cycles that result in substantial alteration of material behavior. In such a framework, accurate models as well as robust and efficient numerical approaches become essential to allow for the simulation of complex devices. The present paper focuses on the numerical simulation of quasi-static problems involving shape memory alloy (SMA) structures or components subjected to multiple loading-unloading cycles. A novel state-update procedure for a three-dimensional phenomenological model able to describe the saturation of permanent inelasticity, including degradation effects, is here proposed. The algorithm, being of the predictor-corrector type and relying on an incremental energy minimization approach, is based on elastic checks, closed-form solutions of polynomial equations, and nonlinear scalar equations solved through a combination of Newton-Raphson and bisection methods. This allows for an easy implementation of model equations and to avoid the use of regularization parameters for the treatment of non-smooth functions. Numerical results assess the good performances of the proposed approach in predicting both pseudoelastic and shape-memory material behavior under cyclic loading as well as algorithm robustness.

Introduction

Shape-memory alloys (SMAs) are metallic alloys possessing the unique properties known as pseudoelasticity (PE) and shape-memory effect (SME). The material is in fact able to recover the original shape through a phase transformation caused by the imposition of a stress (i.e., PE) and/or temperature (i.e., SME) field. Such unique thermo-mechanical properties make SMAs an effective material for several innovative technological applications in the biomedical up to the mechanical field [START_REF] Jani | A review of shape memory alloy research, applications and opportunities[END_REF].

Cyclic loading is one important feature of many of these applications, no matter whether they exploit mechanical or thermal recovery (i.e., PE and one/two-way SME, respectively) [START_REF] Eggeler | Structural and functional fatigue of NiTi shape memory alloys[END_REF]. Typical examples of cyclic loading are the pulsatile blood pressure, that is applied to cardiovascular devices as stents or aortic valves, temperature cycling in actuation components as robotic grippers or thermal valves, or force cycles in damping applications.

A factor that limits the service life of SMA-based applications subjected to cyclic loading is fatigue, both in terms of material integrity (i.e., structural fatigue) and of the change of functional properties and reversibility (i.e., functional fatigue) [START_REF] Robertson | Mechanical fatigue and fracture of Nitinol[END_REF]. The thermo-mechanical response of SMA materials under cyclic loading is however more complex than the response of classical metals, due to the occurrence of phase transformation and plastic deformation, which can lead to different physical situations [START_REF] Feng | Shakedown analysis of shape memory alloy structures[END_REF]. Experimental evidences [START_REF] Van Humbeek | Cycling effects, fatigue and degradation of shape memory alloys[END_REF][START_REF] Kang | Ratcheting deformation of superelastic and shape memory niti alloys[END_REF][START_REF] Miyazaki | Fatigue life of Ti-50 at.% Ni and Ti-40Ni-10Cu (at.%) shape memory alloy wires[END_REF][START_REF] Tobushi | Low-cycle fatigue of TiNi shape memory alloy and formulation of fatigue life[END_REF][START_REF] Wang | Phase transformation behavior of pseudoelastic NiTi shape memory alloys under large strain[END_REF][START_REF] Arrigoni | Cyclic effects in shape-memory alloys: a one-dimensional continuum model[END_REF][START_REF] Owusu-Danquah | On the cyclic stability of the thermomechanical behavior of NiTi shape memory cylindrical actuators[END_REF] have reported that thermal cycling in one-way SME applications suffers a decrease in the exploitable displacement. On the other hand, mechanical cycling in PE components determines an increasing level of permanent deformation, that saturates on a stable value after a certain number of cycles, shifts the hysteresis loop downward, lowers its height and width, and decreases the level of dissipated energy. Such physical evidences originate from the combination of residual martensitic phase and transformationinduced plasticity, that is the formation of microscopic plastic deformation during the stress-induced transformation.

Moreover, such effects are present not only in the widely-used SMAs based on Nickel-Titanium, but also in other types of SMAs, and recent studies have also investigated the behavior of additive-manufactured SMAs [START_REF] Reginald | Multi-scale shape memory effect recovery in niti alloys additive manufactured by selective laser melting and laser directed energy deposition[END_REF].

For these reasons, both understanding the underlying processes and incorporating them in constitutive modeling are of utmost importance to effectively predict material response and to support the design of SMA components.

Several models taking into account the inelastic strain build-up due to not-completed reverse phase transformation or/and plasticity, its accumulation during cyclic loads, and degradation effects, are available from the literature; see, e.g., [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF][START_REF] Jiang | A constitutive model for isothermal pseudoelasticity coupled with plasticity[END_REF][START_REF] Barrera | Macroscopic modeling of functional fatigue in shape memory alloys[END_REF][START_REF] Hartl | Three-dimensional modeling and numerical analysis of rate-dependent irrecoverable deformation in shape memory alloys[END_REF][START_REF] Saint-Sulpice | A 3D super-elastic model for shape memory alloys taking into account progressive strain under cyclic loadings[END_REF][START_REF] Zaki | A 3D model of the cyclic thermomechanical behavior of shape memory alloys[END_REF][START_REF] Peng | A constitutive model for transformation, reorientation and plastic deformation of shape memory alloys[END_REF][START_REF] Bo | Thermomechanical modeling of polycrystalline smas under cyclic loading. Part III: evolution of plastic strains and two-way shape memory effect[END_REF][START_REF] Cheikh | A review of constitutive models and modeling techniques for shape memory alloys[END_REF][START_REF] Waimann | Modeling the cyclic behavior of shape memory alloys[END_REF][START_REF] Paiva | A constitutive model for shape memory alloys considering tensile-compressive asymmetry and plasticity[END_REF][START_REF] Petrini | The role of inelastic deformations in the mechanical response of endovascular shape memory alloy devices[END_REF][START_REF] Govindjee | A shape memory alloy model for uranium-niobium accounting for plasticity[END_REF][START_REF] Lagoudas | Modeling of transformation-induced plasticity and its effect on the behavior of porous shape memory alloys. Part I: constitutive model for fully dense SMAs[END_REF]. Recently, such models have been used in connection with shakedown theorems [START_REF] Peigney | Shakedown theorems and asymptotic behaviour of solids in non-smooth mechanics[END_REF][START_REF] Peigney | On shakedown of shape memory alloys structures[END_REF] and fatigue approaches [START_REF] Gu | Shakedown based model for high-cycle fatigue of shape memory alloys[END_REF][START_REF] Auricchio | A shakedown analysis of high cycle fatigue of shape memory alloys[END_REF].

We shall here focus on the three-dimensional model by Auricchio et al. [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF], later generalized in [START_REF] Barrera | Macroscopic modeling of functional fatigue in shape memory alloys[END_REF], which is capable of describing permanent inelastic effects in both pseudo-elastic and shape-memory behaviors with a low number of physical parameters.

In the modeling framework, it is important to provide a robust and efficient numerical approach to treat model equations and to allow for the simulation of complex devices. Model equations generally involve numerous tensorial and scalar internal variables, subjected to constraints, and include evolution equations in order to describe several physical effects and transformations. Therefore, the numerical implementation in this case is particularly challenging.

In general, the state-update procedures adopted to treat SMA constitutive equations are mostly based on returnmap schemes, e.g., [START_REF] Lagoudas | Constitutive model for the numerical analysis of phase transformation in polycrystalline shape memory alloys[END_REF][START_REF] Zaki | Time integration of a model for martensite detwinning and reorientation under non-proportional loading using Lagrange multipliers[END_REF], while only in the last years the attention towards incremental energy minimization approaches [START_REF] Sedlák | Thermomechanical model for NiTi-based shape memory alloys including R-phase and material anisotropy under multi-axial loadings[END_REF][START_REF] Stupkiewicz | A robust model of pseudoelasticity in shape memory alloys[END_REF][START_REF] Artioli | An incremental energy minimization state update algorithm for 3D phenomenological internal-variable sma constitutive models based on isotropic flow potentials[END_REF][START_REF] Scalet | A robust and efficient radial return algorithm based on incremental energy minimization for the 3D Souza-Auricchio model for shape memory alloys[END_REF] or algorithms for mathematical programming [START_REF] Peigney | Numerical simulation of shape memory alloys structures using interior-point methods[END_REF][START_REF] Peigney | An incremental variational approach to coupled thermo-mechanical problems in anelastic solids. Application to shape-memory alloys[END_REF] is being increased. So far, the solution of the model [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF][START_REF] Barrera | Macroscopic modeling of functional fatigue in shape memory alloys[END_REF] has been performed by means of an elastic-predictor inelastic-corrector return map procedure with a δ-regularized version to control the smoothness of the norm regularization. Since the model includes two tensorial internal variables to describe material behavior, i.e., the transformation and permanent inelastic strains, the return-mapping algorithm involves 10-12 scalar parameters. A standard Newton-Raphson scheme has been adopted in [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF][START_REF] Barrera | Macroscopic modeling of functional fatigue in shape memory alloys[END_REF] to solve the nonlinear system of equations in both the unsaturated and saturated cases.

The model has been tested on uniaxial and biaxial tests in the Matlab environment, but its investigation in a threedimensional finite element (FE) framework is lacking. The large number of scalar parameters may in fact increase computational costs and cause trouble of convergence when using Newton-Raphson procedures.

The aim of the present work is to propose a new time-integration algorithm for the numerical implementation of the model described in [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF][START_REF] Barrera | Macroscopic modeling of functional fatigue in shape memory alloys[END_REF]. As it will be demonstrated, the proposed algorithm can be readily integrated in a finite-element code for solving boundary value problems.

Among the several numerical approaches cited above, the proposed algorithm belongs to the class of variational methods relying on an incremental energy minimization approach. The idea of applying such an approach to SMAs stems from previous works, e.g., [START_REF] Sedlák | Thermomechanical model for NiTi-based shape memory alloys including R-phase and material anisotropy under multi-axial loadings[END_REF][START_REF] Stupkiewicz | A robust model of pseudoelasticity in shape memory alloys[END_REF][START_REF] Artioli | An incremental energy minimization state update algorithm for 3D phenomenological internal-variable sma constitutive models based on isotropic flow potentials[END_REF][START_REF] Scalet | A robust and efficient radial return algorithm based on incremental energy minimization for the 3D Souza-Auricchio model for shape memory alloys[END_REF], and it is here applied to the model described in [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF][START_REF] Barrera | Macroscopic modeling of functional fatigue in shape memory alloys[END_REF]. The incremental energy minimization approach has been successfully applied in [START_REF] Scalet | A robust and efficient radial return algorithm based on incremental energy minimization for the 3D Souza-Auricchio model for shape memory alloys[END_REF] to the original model [START_REF] Souza | Three-dimensional model for solids undergoing stress-induced phase transformations[END_REF], in which the transformation strain is the only internal variable and permanent inelastic strains are not taken into account. The algorithm developed in [START_REF] Scalet | A robust and efficient radial return algorithm based on incremental energy minimization for the 3D Souza-Auricchio model for shape memory alloys[END_REF] is here extended to the model under consideration, taking into account both transformation-induced strains and permanent inelasticity. In the present case, two tensorial internal variables (i.e., the transformation strain and the permanent inelastic strain) are introduced and their evolution in a finite time step incrementally minimizes a convex functional, given by the sum of the free-energy energy and the dissipation functional. The proposed algorithm is based on elastic checks, closed-form solutions of polynomial equations, and nonlinear scalar equations solved through a combination of Newton-Raphson and bisection methods. This allows to avoid the difficulties mentioned above, when using Newton-Raphson procedures. This suitable variational structure facilitates the treatment of internal constraints and allows for an efficient numerical implementation. Other advantages of the proposed algorithm are its easy implementation and, overall, the possibility of avoiding regularized terms in both energy/dissipation definition and norms, that may affect material response as well as numerical convergence.

To test the performance and robustness of the proposed algorithm several FE analyses are presented. The simulations range from classical uniaxial tests to more complex representative problems, involving both pseudoelasticity and shape-memory effect. The complex problems involve the three-dimensional analyses of a stent strut and of a spring actuator, given the importance of such devices under cyclic loading conditions [START_REF] Pelton | Fatigue and durability of Nitinol stents[END_REF][START_REF] Elahinia | Shape memory alloy actuators: Design, fabrication and experimental evaluation[END_REF]. The results show that the model implemented with the proposed algorithm is able to catch material response for several sets of materials parameters and different time steps.

The paper is organized as follows. Section 2 briefly reviews the continuum equations of the model under investigation. Section 3 presents the equations in the time-discrete framework and describes the proposed algorithmic scheme.

Then, Section 4 presents the results of several numerical simulations. Finally, conclusions are given in Section 5.

Model equations

This section briefly recalls main model continuum equations in the small strain regime, as presented in [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF][START_REF] Barrera | Macroscopic modeling of functional fatigue in shape memory alloys[END_REF].

The model assumes the total strain ε and the absolute temperature T as control variables, while the transformation strain e tr and the permanent inelastic strain q as internal ones. Both e tr and q are symmetric trace-free second order tensors. Specifically, the transformation strain e tr describes the strain associated to the austenite-martensite phase transformation and the permanent inelastic strain q gives a measure of the part of e tr that cannot be recovered when unloading to a zero stress state, since e tr has no fully reversible evolution. The transformation strain e tr is required to satisfy the constraint:

e tr ≤ ε L (1)
ε L being a material parameter corresponding to the maximum transformation strain reached at the end of the transformation during a uniaxial test. The norm • in Eq. ( 1) is the Euclidean norm, as defined by e tr = √ e tr : e tr where : denotes the contraction with respect to the last two indices (e.g., a : b = ij a ij b ji ).

The Helmholtz free-energy density function is expressed as

Ψ = Ψ 0 + I ε L (2)
where

Ψ 0 (ε, e tr , q) = 1 2 K θ 2 + G e -e tr 2 + τ M e tr -q + 1 2 H e tr 2 + 1 2 h q 2 -Ae tr : q (3)
and

I ε L (e tr ) =      0 if e tr ≤ ε L +∞ otherwise.
In [START_REF] Robertson | Mechanical fatigue and fracture of Nitinol[END_REF], θ and e are the volumetric and the deviatoric part of ε; K and G are, respectively, the bulk and the shear modulus; τ M = β T -T 0 , where β is a material parameter related to the dependence of the critical stress on the temperature, T 0 is the temperature below which only martensite phase is stable, and • is the positive part function;

H, h, and A define, respectively, the hardening of the phase transformation, the saturation of the permanent inelastic strain evolution, and model degradation. The energy term I ε L (e tr ) is the indicator function associated with the constraint [START_REF] Eggeler | Structural and functional fatigue of NiTi shape memory alloys[END_REF].

For later reference, we note that the energy in Eq. ( 2) is strictly convex provided that:

hH -A 2 > 0. (4)
In the following, the condition ( 4) is assumed to be satisfied.

The dissipation function originally considered in ref. [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF] is defined as:

Φ( εtr , q) = R Y max( εtr , κ q ) (5)
where R Y and κ are non-negative material parameters1 . We note that Eq. ( 5) can be written as follows:

Φ( εtr , q) = R Y ( εtr , q) κ,∞
where:

(v 1 , v 2 ) κ,∞ = max( v 1 , κ v 2 )
is a weighted supremum norm. As detailed by Barrera et al. (15), other choices of norms can be made, leading to other expressions of (rate-independent) dissipation functions. In particular, the supremum norm • κ,∞ in Eq. ( 5) could be replaced by the weighted taxicab norm:

(v 1 , v 2 ) κ,1 = v 1 + κ v 2
or the weighted Euclidean norm:

(v 1 , v 2 ) κ,2 = v 1 2 + κ 2 v 2 2 .
As explained in ref. [START_REF] Barrera | Macroscopic modeling of functional fatigue in shape memory alloys[END_REF], the norms • κ,1 and • κ,2 lead to results that are more consistent with experiments than the norm • κ,∞ . In the following, we choose the norm • κ,2 which seems to be the simplest one to handle for three-dimensional numerical implementation. To alleviate the notations, the norm • κ,2 is denoted by • κ from now on. The dissipation function is thus assumed to be of the form:

Φ = R Y ( ėtr , q) κ = R Y ėtr 2 + κ 2 q 2 . ( 6
)
Following standards arguments, the stress-strain relation is obtained by differentiating the free energy function Ψ with respect to the strain ε, yielding:

p = ∂Ψ ∂θ = Kθ, s = ∂Ψ ∂e = 2G(e -e tr ) (7)
where p and s are the hydrostatic and deviatoric part of the stress σ, respectively. Similarly, the thermodynamic forces (X, Q) associated to the internal variables (e tr , q) are usually defined by the relations X = -∂Ψ/∂e tr and Q = -∂Ψ/∂q. In the present case, however, special care must be taken because Ψ is only subdifferentiable in (e tr , q).

In such case, the usual definition needs to be amended as -(X, Q) ∈ ∂Ψ where ∂ denotes the subdifferential operator with respect to (e tr , q). It follows that -(X, Q) ∈ (-s + He tr -Aq, hq -Ae tr ) + τ M ∂ e trq + ∂I ε L (e tr ).

The reader is referred to, e.g., [START_REF] Brézis | Opérateurs maximum monotones et semigroupes de contractions dans les espaces de Hilbert[END_REF][START_REF] Rockafellar | Convex analysis[END_REF] for an in-depth introduction to subdifferentials and related tools in convex analysis. We simply recall here that the subdifferential of a convex function F (e tr , q) is the multi-valued operator ∂F defined by: ∂F (e tr , q) = {(a, b) : F (ẽ tr , q) ≥ F (e tr , q) + a : (ẽ tre tr ) + b : (q -q) ∀(ẽ tr , q)}.

In particular, we have:

∂ e tr -q =     
(e trq, qe tr ) e trq if e tr = q

(τ , -τ ) : tr τ = 0, τ ≤ 1 if e tr = q (9)
and

∂I ε L (e tr ) =              (0, 0) if e tr < ε L (γe tr , 0) : γ ≥ 0 if e tr = ε L ∅ if e tr > ε L (10)
Note that ∂ e trq is multi-valued when e tr = q. Similarly, ∂I ε L (e tr ) is multi-valued when e tr = ε L .

The evolution equation for (e tr , q) is determined by the dissipation function Φ as:

(X, Q) ∈ ∂Φ( εtr , q). (11)
Following the framework of standard generalized materials [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF], Eq. ( 11) respects the second law of thermodynamics for any choice of positive, convex dissipation function that vanishes at the origin. Using expression (6) of the dissipation function yields:

∂Φ( εtr , q) =        R Y εtr , κ 2 q ( εtr , q) κ if ( εtr , q) = (0, 0) C if ( εtr , q) = (0, 0) (12)
with:

C = {(τ 1 , τ 2 ) : tr τ 1 = tr τ 2 = 0, τ 1 2 + 1 κ 2 τ 2 2 ≤ R 2 Y }. (13)
Relation [START_REF] Owusu-Danquah | On the cyclic stability of the thermomechanical behavior of NiTi shape memory cylindrical actuators[END_REF] can be rewritten in a more familiar form by noting that it is equivalent to:

( εtr , q) ∈ ∂Φ * (X, Q) (14)
where Φ * is the Legendre transform of the dissipation function Φ, as defined by Φ * (X, Q) = sup ( εtr , q) X : εtr + Q :

q -Φ( εtr , q). Using expression [START_REF] Kang | Ratcheting deformation of superelastic and shape memory niti alloys[END_REF], it can be calculated that Φ * is the indicator function of the domain C in [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF]. It follows that Eq. ( 14) becomes:

( εtr , q) = λ X X , 1 κ 2 Q Q (15)
with conditions:

λ ≥ 0, X 2 + 1 κ 2 Q 2 -R 2 Y ≤ 0, λ( X 2 + 1 κ 2 Q 2 -R 2 Y ) = 0. (16)
Eqs. ( 15) and ( 16) correspond to a normality flow rule for the variables (e tr , q). Correspondingly, the domain C introduced in Eq. ( 13) can be interpreted as the elasticity domain of the material. In the space of (X, Q) variables, the domain C has an ellipsoidal shape with axis R Y and κR Y .

For later reference, we note that Eqs. ( 8) and ( 11) can be combined as:

(0, 0) ∈ (-s + He tr -Aq, hq -Ae tr ) + τ M ∂ e trq + ∂I ε L (e tr ) + ∂Φ( εtr , q). (17)

Incremental algorithm

This section discusses the time-discretization of the constitutive laws [START_REF] Miyazaki | Fatigue life of Ti-50 at.% Ni and Ti-40Ni-10Cu (at.%) shape memory alloy wires[END_REF] and [START_REF] Saint-Sulpice | A 3D super-elastic model for shape memory alloys taking into account progressive strain under cyclic loadings[END_REF]. Time discretization consists in introducing a finite number of time instants t 0 < • • • < t N and estimating the state at each time instant t n in a time-marching approach. Let p n , s n , e tr n , q n be respectively the hydrostatic stress, deviatoric stress, transformation strain, and permanent inelastic strain at time t n . We focus on the central issue of estimating the state (p n+1 , s n+1 , e tr n+1 , q n+1 ) at current time t n+1 , assuming that : (i) the control variables at current time t n+1 (i.e, the total strain ε n+1 and the temperature T n+1 ) are prescribed, (ii) the state (e tr n , q n ) at previous time t n is known and satisfies the constraint e tr n ≤ ε L . A natural way of performing the state update is to discretize Eq. ( 17) using an implicit Euler scheme as:

(0, 0) ∈ (-s n+1 +He tr n+1 -Aq n+1 , hq n+1 -Ae tr n+1 )+τ M,n+1 ∂ e tr n+1 -q n+1 +∂I ε L (e tr n+1 )+∂Φ( e tr n+1 -e tr n t n+1 -t n , q n+1 -q n t n+1 -t n )

with

p n+1 = Kθ n+1 , s n+1 = 2G(e n+1 -e tr n+1 ), τ M,n+1 = β T n+1 -T 0 . (19)
In Eq. ( 19), θ n+1 and e n+1 are the volumetric and deviatoric part of ε n+1 , respectively. To alleviate the notations, the scalar τ M,n+1 will be denoted by τ M in the following. Eliminating s n+1 between Eqs. ( 18) and [START_REF] Peng | A constitutive model for transformation, reorientation and plastic deformation of shape memory alloys[END_REF] and noting that ∂Φ is positively homogeneous of degree 0, we obtain the equation:

(0, 0) ∈ (2G e tr n+1 -a, hq n+1 ) + τ M ∂ e tr n+1 -q n+1 -A(q n+1 , e tr n+1 ) + ∂I ε L (e tr n+1 ) + ∂Φ(e tr n+1 -e tr n , q n+1 -q n ) (

where G = G + H/2 and a = 2Ge n+1 . In Eq. ( 20), the unknowns are the two deviatoric tensors (e tr n+1 , q n+1 ).

In a finite-element framework, Eq. ( 20) typically needs to be solved at each Gauss point. A computationally efficient algorithm for solving Eq. ( 20) is thus essential. Aside from nonlinearity, there are two main difficulties in solving Eq. [START_REF] Bo | Thermomechanical modeling of polycrystalline smas under cyclic loading. Part III: evolution of plastic strains and two-way shape memory effect[END_REF]. The first one is that several terms in Eq. ( 20) are not differentiable, meaning that those terms are multi-valued for certain values of (e tr n+1 , q n+1 ). The second difficulty lies in the dimensionality of the problem: compared (for instance) to the Souza-Auricchio model [START_REF] Souza | Three-dimensional model for solids undergoing stress-induced phase transformations[END_REF][START_REF] Auricchio | A three-dimensional model describing stress-temperature induced solid phase transformations: solution algorithm and boundary value problems[END_REF], there are now two internal variables (e tr , q) instead of one (namely, e tr ), hence the dimensionality jumps from 5 to 10.

For our purpose, it is important to note that a variational formulation is attached to Eq. [START_REF] Bo | Thermomechanical modeling of polycrystalline smas under cyclic loading. Part III: evolution of plastic strains and two-way shape memory effect[END_REF]. Consider indeed the convex function F(e tr , q) defined by F(e tr , q) = Ψ(ε n+1 , e tr , q) + Φ(e tre tr n , qq n )

where Ψ is the Helmoltz energy function and Φ is the dissipation potential, as introduced in ( 2) and ( 6) respectively.

Since Ψ and Φ are proper, lower-semicontinuous and convex, the function F is also proper, lower-semicontinuous and convex [START_REF] Brézis | Opérateurs maximum monotones et semigroupes de contractions dans les espaces de Hilbert[END_REF][START_REF] Rockafellar | Convex analysis[END_REF]. Hence the solutions to the minimization problem inf (e tr , q) F(e tr , q)

are characterized by the optimality condition (0, 0) ∈ ∂F(e tr , q). (

Moreover, we have ∂F(e tr , q) = ∂Ψ(ε n+1 , e tr , q) + ∂Φ(e tre tr n , qq n ) where ∂Ψ is the subdifferential of Ψ with respect to the variables (e tr , q). From expression (2) we find ∂F(e tr , q) = (2G e tra, hq) + τ M ∂ e trq -A(q, e tr ) + ∂I ε L (e tr ) + ∂Φ(e tre tr n , qq n ) (

Any (e tr n+1 , q n+1 ) satisfying Eq. ( 20) clearly verifies Eqs ( 22) and ( 23), i.e. is a solution to the minimization problem [START_REF] Cheikh | A review of constitutive models and modeling techniques for shape memory alloys[END_REF]. Conversely, any solution to [START_REF] Cheikh | A review of constitutive models and modeling techniques for shape memory alloys[END_REF] verifies Eq. ( 20). The formulations [START_REF] Bo | Thermomechanical modeling of polycrystalline smas under cyclic loading. Part III: evolution of plastic strains and two-way shape memory effect[END_REF] and ( 21) are thus equivalent.

Eq. ( 21) can be interpreted as an incremental energy minimization problem. The state update equation ( 20) of the Euler implicit scheme can be interpreted as the optimality condition in the minimization problem [START_REF] Cheikh | A review of constitutive models and modeling techniques for shape memory alloys[END_REF].

The variational formulation [START_REF] Cheikh | A review of constitutive models and modeling techniques for shape memory alloys[END_REF] notably allows one to justify the existence and uniqueness of the updated state (e tr n+1 , q n+1 ) introduced in Eq. [START_REF] Bo | Thermomechanical modeling of polycrystalline smas under cyclic loading. Part III: evolution of plastic strains and two-way shape memory effect[END_REF]. Under the requirement (4), the function F is indeed strictly convex, with infinite growth at infinity, and therefore admits a unique minimizer. The solution (e tr n+1 , q n+1 ) to ( 20) is thus uniquely defined. In the following, we will take advantage of the variational formulation [START_REF] Cheikh | A review of constitutive models and modeling techniques for shape memory alloys[END_REF] for determining (e tr n+1 , q n+1 ) in a robust and consistent manner. The overall strategy that we propose breaks down into three steps. The first step (see Section 3.1) consists in checking whether (e tr n+1 , q n+1 ) = (e tr n , q n ). If not, the second step (see Section 3.2) consists in solving the minimization problem obtained by ignoring the term I ε L (e tr ) in Eq. ( 2). If the state (ẽ tr , q) obtained in such fashion satisfies the constraint ẽtr ≤ ε L , then we have (e tr n+1 , q n+1 ) = (ẽ tr , q). Otherwise, we move to the third and final step (see Section 3.3) which consists in solving Eq. ( 22) with respect to pairs (e tr n+1 , q n+1 ) that saturates the constraint, i.e., such that e tr n+1 = ε L .

Elastic evolution

In this section we detail the procedure for checking whether (e tr n+1 , q n+1 ) = (e tr n , q n ) i.e., whether the incremental evolution is elastic. From Eqs.( 12) and ( 20), the condition for the evolution to be elastic is:

(b, -c) ∈ τ M ∂ e tr n -q n + ∂I ε L (e tr ) + C (24)
where:

b = a -2G e tr n + Aq n , c = hq n -Ae tr n .
Detailing condition [START_REF] Petrini | The role of inelastic deformations in the mechanical response of endovascular shape memory alloy devices[END_REF] any further requires to distinguish between different cases, depending on the value of (e tr n , q n ).

3.1.1

Case e tr n = q n with e tr n < ε L

We first consider the situation where e tr n = q n and e tr n < ε L . Using Eqs. ( 9) and [START_REF] Arrigoni | Cyclic effects in shape-memory alloys: a one-dimensional continuum model[END_REF], condition [START_REF] Petrini | The role of inelastic deformations in the mechanical response of endovascular shape memory alloy devices[END_REF] becomes:

(u, v) ∈ C (25)
where u and v are defined by

u = b -τ M e tr n -q n e tr n -q n , v = -c + τ M e tr n -q n e tr n -q n .
From the expression of C in Eq. ( 13), condition (25) can be rewritten as:

κ 2 u 2 + v 2 ≤ κ 2 R 2 Y . (26)
3.1.2 Case e tr n = q n with e tr n = ε L

We now consider the situation where e tr n = q n and e tr n = ε L . In such case, by Eqs. ( 9) and ( 10), condition ( 24) is satisfied if and only if (u -γe tr n , v) ∈ C for some γ ≥ 0. This is equivalent to requiring that:

inf γ≥0 κ 2 u -γe tr n 2 + v 2 ≤ κ 2 R 2 Y . (27)
Let u ⊥ = u -(u : e tr n )e tr n /ε 2 L be the projection of u on the orthogonal of e tr n , so that:

u -γe tr n 2 = (u : e tr n )/ε L -γε L 2 + u ⊥ 2 .
It follows that:

inf γ≥0 u -γe tr n 2 = u - u : e tr n ε 2 L e tr n 2 =      u ⊥ 2 if u : e tr n ≥ 0, u 2 if u : e tr n ≤ 0.
Substituting in Eq. ( 27), we obtain that the condition for (e tr n , q n ) to be the solution to problem (21) reads as:

κ 2 u - u : e tr n ε 2 L e tr n 2 + v 2 ≤ κ 2 R 2 Y . (28)

3.1.3

Case e tr n = q n with e tr n < ε L

In the case e tr n = q n with e tr n < ε L , condition ( 24) is satisfied if and only if:

(b -τ , -c + τ ) ∈ C (29)
for some deviatoric tensor τ such that τ ≤ τ M . Condition ( 29) can be rewritten as :

inf τ : τ ≤τ M κ 2 b -τ 2 + c -τ 2 ≤ κ 2 R 2 Y . (30)
Using expansion:

κ 2 b -τ 2 + c -τ 2 = (κ 2 + 1) τ 2 -2τ : (κ 2 b + c) + κ 2 b 2 + c 2
as well as the Cauchy-Schwarz inequality, it can be easily seen that the tensor τ reaching the infimum in problem [START_REF] Auricchio | A shakedown analysis of high cycle fatigue of shape memory alloys[END_REF] is positively collinear to κ 2 b + c. Condition (30) can thus be simplified as:

F 1 ≤ κ 2 R 2 Y ( 31
)
where

F 1 = min 0≤t≤τ M (κ 2 + 1)t 2 -2t κ 2 b + c + κ 2 b 2 + c 2 . ( 32
)
Solving the quadratic minimization problem defining F 1 leads to the following expressions:

F 1 =      κ 2 1 + κ 2 b -c 2 if κ 2 b + c ≤ (κ 2 + 1)τ M , (κ 2 + 1)τ 2 M -2τ M κ 2 b + c + κ 2 b 2 + c 2 otherwise.

3.1.4

Case e tr n = q n with e tr n = ε L

We finally consider the situation where e tr n = q n and e tr n = ε L . In that case, condition ( 24) is satisfied if and only if

(b -τ -γe tr n , -c + τ ) ∈ C (33)
for some γ ≥ 0 and some deviatoric tensor τ such that τ ≤ τ M . Condition (33) can be rewritten as:

F 2 ≤ κ 2 R 2 Y ( 34
)
where

F 2 = inf γ ≥ 0, τ : τ ≤ τ M κ 2 b -τ -γe tr n 2 + c -τ 2 .
The quadratic minimization problem defining F 2 can be solved in closed form, although the expressions are more involved than those obtained for F 1 in Eq. [START_REF] Zaki | Time integration of a model for martensite detwinning and reorientation under non-proportional loading using Lagrange multipliers[END_REF]. We set:

b / / = a : e tr n ε L + (A -2G )ε L , c / / = (h -A)ε L , d / / = τ M κ 2 b / / + c / / κ 2 b + c , B = a 2 -a : e tr n ε L 2 .
Table 1: Conditions for elastic evolution.

e tr n < ε L e tr n = ε L e tr n = q n κ 2 u 2 + v 2 ≤ κ 2 R 2 Y κ 2 u - u : e tr n ε 2 L e tr n 2 + v 2 ≤ κ 2 R 2 Y e tr n = q n F 1 ≤ κ 2 R 2 Y F 2 ≤ κ 2 R 2 Y
Omitting the detail of the calculations, F 2 is given by the following expressions:

-if κ 2 b + c ≥ (κ 2 + 1)τ M : F 2 =              (κ 2 + 1)τ 2 M -2τ M κ 2 b + c + κ 2 b 2 + c 2 if b / / ≤ d / / κ 2 1 + κ 2 B if b / / ≥ d / / and c 2 / / + κ 2 1 + κ 2 2 B ≤ τ 2 M τ 2 M x 2 -1 c 2 / / κ 2 + x 2 + (c / / -x) 2 if b / / ≥ d / / and c 2 / / + κ 2 1 + κ 2 2 B ≥ τ 2 M -if κ 2 b + c ≤ (κ 2 + 1)τ M : F 2 =              - 1 1 + κ 2 κ 2 b + c 2 + κ 2 b 2 + c 2 if b / / ≤ c / / κ 2 1 + κ 2 B if b / / ≥ c / / and c 2 / / + κ 2 1 + κ 2 2 B ≤ τ 2 M τ 2 M x 2 -1 c 2 / / κ 2 + x 2 + (c / / -x) 2 if b / / ≥ c / / and c 2 / / + κ 2 1 + κ 2 2 B ≥ τ 2 M
The scalar x that appears in the above expressions is obtained by solving the polynomial equation

(x 2 -τ 2 M )(c / / + κ 2 x) 2 + κ 4 Bx 2 = 0 (35) on [-τ M , min(τ M , b / / )].
Since the polynomial in Eq. ( 35) is of degree 4, its roots can be obtained in closed form.

Summary

The conditions for the incremental evolution to be elastic are summarized in Table 1. If those conditions are satisfied, then the solution (e tr n+1 , q n+1 ) to problem [START_REF] Bo | Thermomechanical modeling of polycrystalline smas under cyclic loading. Part III: evolution of plastic strains and two-way shape memory effect[END_REF] is simply obtained as (e tr n+1 , q n+1 ) = (e tr n , q n ). We emphasize that all the calculations needed for checking the conditions in Table 1 can be done in closed form, without resorting to a nonlinear solver, as in the case of ref. [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF]. In the most complex case ( e tr n = ε L with e tr n = q n and c 2 / / +

κ 2 1 + κ 2 2 B ≥ τ 2 M
), a polynomial equation of degree 4 needs to be solved.

Unsaturated phase transformation

If the conditions reported in Table 1 are not satisfied, then (e tr n+1 , q n+1 ) = (e tr n , q n ), i.e., phase transformation occurs.

To carry out the state update in that case, we first consider the incremental energy minimization problem obtained by dropping the term I ε L (e tr ) in [START_REF] Eggeler | Structural and functional fatigue of NiTi shape memory alloys[END_REF]. The corresponding minimization problem can be written as inf (e tr , q) F 0 (e tr , q)

where F 0 (e tr , q) = Ψ 0 (ε n+1 , e tr , q) + Φ(e tre tr n , qq n ) (37

)
and Ψ 0 is defined as in Eq. ( 3). For later reference, we note that problem (36) can be written in a more explicit form as inf (e tr , q) G e tr 2 -e tr : a + τ M e trq + 1 2 h q 2 -Ae tr : q + Φ(e tre tr n , qq n ).

Since F 0 is strictly convex and grows to infinity as (e tr , q) tends to infinity, problem (36) admits a unique solution that we denote by (ẽ tr , q). The latter is characterized by the optimality condition (0, 0) ∈ F 0 (ẽ tr , q), i.e.

(0, 0) ∈ (2G ẽtra, hq) + τ M ∂ ẽtrq -A(q, ẽtr ) + ∂Φ(ẽ tre tr n , qq n ).

(39)

Case ẽtr = q

In order to find (ẽ tr , q), the strategy that we propose consists in first carrying out the optimization with respect to pairs (e tr , q) such that e tr = q. We thus consider the problem

inf q F 0 (q, q). (40)
The minimization problem ( 41) is easier to solve than problem [START_REF] Scalet | A robust and efficient radial return algorithm based on incremental energy minimization for the 3D Souza-Auricchio model for shape memory alloys[END_REF], because it involves only one tensorial unknown, instead of two. From expression (37) of F 0 , (40) can be rewritten as

inf q 1 2 G q 2 -q : a + Φ(q -e tr n , q -q n ) (41)
where

G = 2G + h -2A.
Let q * be the (unique) solution to problem [START_REF] Elahinia | Shape memory alloy actuators: Design, fabrication and experimental evaluation[END_REF]. If e tr n = q n , then q * is characterized by the optimality condition:

0 = G q * -a + R Y y ((1 + κ 2 )q * -e tr n -κ 2 q n ) (42)
with

y = (q * -e tr n , q * -q n ) κ . (43)
From Eq. ( 42) we obtain :

q * = ya + R Y (e tr n + κ 2 q n ) yG + R Y (1 + κ 2 ) . ( 44
)
Substituting [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF] in [START_REF] Rockafellar | Convex analysis[END_REF] gives an equation in which y is the only unknown. After some manipulation, that equation is found to read as

a 0 + a 1 y + a 2 y 2 + a 3 y 3 + a 4 y 4 = 0 (45)
with:

a 0 = κ 2 R 2 Y (1 + κ 2 ) e tr n -q n 2 , a 1 = 2R Y G κ 2 e tr n -q n 2 , a 2 = (G e tr n -a, G q n -a) 2 κ -R 2 Y (1 + κ 2 ) 2 , a 3 = -2R Y (1 + κ 2 )G , a 4 = -G 2 .
Since the polynomial equation ( 45) is of degree 4, the value of y can be obtained in closed form. Substituting the result in Eq. ( 44) gives the value of q * . The solution q * to problem [START_REF] Elahinia | Shape memory alloy actuators: Design, fabrication and experimental evaluation[END_REF] being found, we proceed to check if (ẽ tr , q) = (q * , q * ), i.e., if (q * , q * ) is the solution to problem [START_REF] Scalet | A robust and efficient radial return algorithm based on incremental energy minimization for the 3D Souza-Auricchio model for shape memory alloys[END_REF]. Eq. [START_REF] Souza | Three-dimensional model for solids undergoing stress-induced phase transformations[END_REF] shows that it happens to be the case if:

(h -A)q * + R Y y κ 2 (q * -q n ) ≤ τ M . (46)
Remark: in the case e tr n = q n , one first needs to check whether:

G q n -a ≤ R Y 1 + κ 2 . ( 47
)
If condition ( 47) is satisfied, then q * = q n so that (q * , q * ) = (e tr n , q n ). We already know from Section 3.1 that (e tr n , q n ) is not the solution to problem [START_REF] Cheikh | A review of constitutive models and modeling techniques for shape memory alloys[END_REF]. It follows that (q * , q * ) = (e tr n , q n ) is not the solution to problem [START_REF] Scalet | A robust and efficient radial return algorithm based on incremental energy minimization for the 3D Souza-Auricchio model for shape memory alloys[END_REF] neither.

If condition (47) is not satisfied, q * is obtained by expressions ( 44)-( 45), and condition [START_REF] Vrahatis | From linear to nonlinear iterative methods[END_REF] for checking if (ẽ tr , q) = (q * , q * ) still applies.

Case ẽtr = q

If (q * , q * ) does not provide the solution to problem [START_REF] Scalet | A robust and efficient radial return algorithm based on incremental energy minimization for the 3D Souza-Auricchio model for shape memory alloys[END_REF], then we necessarily have ẽtr = q and the optimality conditions (39) become:

       0 = 2G ẽtr -a +τ M ẽtr - q x -Aq + R Y ẽtr -e tr n y 0 = hq +τ M q -ẽtr x -Aẽ tr + κ 2 R Y q -q n y (48)
where x = ẽtrq and y = (ẽ tre tr n , qq n ) κ . For given x and y, Eq. ( 48) can be viewed as a linear system in (ẽ tr , q) and put in matrix form:

   2G + x + y -x -A -x -A x + h + κ 2 y       ẽtr q    =    a + y e tr n κ 2 y q n    (49)
where x = τ M /x and y = R Y /y. Let M be the 2×2 matrix that appears in the left-hand side of system [START_REF] Boatti | Development of an effective and user-friendly numerical framework for the simulation of complex smart material components and devices[END_REF]. We have:

det M = 2G h -A 2 + x G + x y (1 + κ 2 ) + y (2κ 2 G + h) + κ 2 y 2 . ( 50
)
Using condition (4), it can easily be checked that 2G h -A 2 > 0 so that det M > 0 for any positive (x, y). System (49) can be thus inverted to give:

   ẽtr q    = 1 det M    x + h + κ 2 y x + A x + A 2G + x + y       a + y e tr n κ 2 y q n    . (51)
Through relations (51), (ẽ tr , q) are expressed as explicit functions of the two unknown scalars (x, y). It remains to find the value of (x, y). To do so, we note from system (51) that:

ẽtr -q = 1 det M d (52) with d = (h + κ 2 y -A)(a + y e tr n ) + (A -2G -y )κ 2 y q n .
A crucial observation is that tensor d is independent on x. We can thus use Eq. ( 52) to obtain x as a function of y.

More precisely, taking the norm of Eq. (52) gives:

x det M = d .
Using Eq. ( 50), we find:

x = d -τ M (G + y (1 + κ 2 )) 2G h -A 2 + y (2κ 2 G + h) + κ 2 y 2 (53)
and substituting in system (51) yields:

   ẽtr q    = 1 Y     (h + κ 2 y ) - τ M d (h -A + κ 2 y ) 2 A + τ M d (A -h -κ 2 y )(A -2G -y ) A + τ M d (A -h -κ 2 y )(A -2G -y ) (2G + y ) - τ M d (A -2G -y ) 2        a + y e tr n κ 2 y q n    (54)
where

Y = 2G h -A 2 + y (2κ 2 G + h) + κ 2 y 2 .
We are now left with the issue of finding the scalar y. This is accomplished by solving the equation:

f (y) = 0 (55)
where:

f (y) = (ẽ tr -e tr n , q -q n ) κ -y
in which (ẽ tr , q) are expressed as functions of y via system (54). Since (ẽ tr , q) is uniquely defined, the equation f (y) = 0 has a unique solution in ]0, +∞[. In practice, the nonlinear equation (55) needs to be solved numerically.

It can be verified that f (y) → -∞ as y → +∞ and that f (y) converges towards a non negative value as y → 0. In practical computations, those properties can be useful for initializing a bisection algorithm.

Saturated phase transformation

If the values (ẽ tr , q) found in Section 3.2 are such that ẽtr ≤ ε L , then (ẽ tr , q) satisfies the optimality condition [START_REF] Bo | Thermomechanical modeling of polycrystalline smas under cyclic loading. Part III: evolution of plastic strains and two-way shape memory effect[END_REF] so that (e tr n+1 , q n+1 ) = (ẽ tr , q). Otherwise, the solution (e tr n+1 , q n+1 ) to problem (21) necessarily saturates the constraint, i.e., verifies e tr n+1 = ε L .

Case e tr n+1 = q n+1

As in Section 3.2.1, we first carry out the optimization with respect to pairs (e tr , q) such that e tr = q. The strategy consists in determining the (unique) solution q * to the minimization problem

inf q F(q, q) (56)
and check whether (e tr n+1 , q n+1 ) = (q * , q * ). Problem (56) can equivalently be rewritten as

inf q: q ≤ε L 1 2 G q 2 -q : a + Φ(q -e tr n , q -q n ). (57)
Observe that if the value q * introduced in Section 3.2.1 satisfies q * ≤ ε L , then q * = q * . In that case, we already know from Section 3.2 that (e tr n+1 , q n+1 ) = (q * , q * ). In the following, we examine the situation where q * > ε L . In that case, q * necessarily saturates the constraint in (57) , i.e.

q * = ε L .

In the special case where e tr n = q n with e tr n = ε L , it is possible that q * = q n . This occurs if:

√ B ≤ R Y 1 + κ 2 . ( 58
)
In that case, as we know from Section 3.1 that (e tr n , q n ) is not the solution to problem (21), we have again (e tr n+1 , q n+1 ) = (q * , q * ).

Except in the very special case where e tr n = q n , e tr n = ε L , and inequality (58) are satisfied, the optimality condition in (57) becomes

0 = -a + R Y y ((1 + κ 2 )q * -e tr n -κ 2 q n ) + γ * q * ( 59
)
with γ * ≥ 0 and y = (q * -e tr n , q * -q n ) κ .

(60)

It follows that:

q * = ε L ya + R Y (e tr n + κ 2 q n ) ya + R Y (e tr n + κ 2 q n ) . ( 61
)
Subsituting ( 61) in (60) leads to the equation

g(y) = 0 (62)
where

g(y) = y 2 + 2ε L ya : u 0 + R Y u 0 2 ya + R Y u 0 -ε 2 L (1 + κ 2 ) -e tr n 2 -κ 2 q n 2 and u 0 = e tr n + κ 2 q n .
In practice, Eq. (62) needs to be solved numerically. Substituting the obtained value for y in (61) gives q * . Having found the solution q * to problem (56), we check whether (q * , q * ) happens to be the solution to problem [START_REF] Cheikh | A review of constitutive models and modeling techniques for shape memory alloys[END_REF]. Using

Eqs. ( 10)- [START_REF] Bo | Thermomechanical modeling of polycrystalline smas under cyclic loading. Part III: evolution of plastic strains and two-way shape memory effect[END_REF], we obtain that (q * , q * ) is the solution to problem [START_REF] Cheikh | A review of constitutive models and modeling techniques for shape memory alloys[END_REF] if:

(h -A)q * + R Y y κ 2 (q * -q n ) ≤ τ M , a + R Y y u 0 ≥ (G + R Y y (1 + κ 2 ))ε L . (63)
If condition (63) is satisfied, then (e tr n+1 , q n+1 ) = (q * , q * ).

Case e tr n+1 = q n+1

If the procedure described in Section 3.3.1 does not the provide the solution to problem [START_REF] Cheikh | A review of constitutive models and modeling techniques for shape memory alloys[END_REF], then e tr n+1 = q n+1 and e tr n+1 = ε L . Hence the optimality condition (20) becomes:

       0 = 2G e tr n+1 -a +τ M e tr n+1 -q n+1 x -Aq n+1 + R Y e tr n+1 -e tr n y +γe tr n+1 0 = hq n+1 +τ M q n+1 -e tr n+1 x -Ae tr n+1 + κ 2 R Y q n+1 -q n y (64)
with

x = e tr n+1 -q n+1 , y = (e tr n+1 -e tr n , q n+1 -q n ) κ , γ ≥ 0.

The solution to problem (21) can be obtained by solving the nonlinear problem:

h(γ) = 0 (65)
where:

h(γ) = e tr (γ) -ε L
and (e tr (γ), q(γ)) denotes the solution to the unconstrained problem (parameterized by γ).

inf (e tr ,q) 1 2 γ e tr 2 + G e tr 2 -e tr : a + τ M e trq + 1 2 h q 2 -Ae tr : q + Φ(e tre tr n , qq n ).

Problem ( 66) is formally identical to problem [START_REF] Peigney | An incremental variational approach to coupled thermo-mechanical problems in anelastic solids. Application to shape-memory alloys[END_REF], the only difference being that the term G e tr 2 in Eq. ( 38) is replaced by ( 1 2 γ + G ) e tr 2 . Consequently, the method presented in Section 3.2 can be directly used for solving problem (66). The updated state is obtained as (e tr n+1 , q n+1 ) = (e tr (γ), q(γ)), where γ is the solution to problem (65).

Summary

The pseudocode of the proposed algorithm in summarized below in Algorithm 1. The input variables are the internal variables (e tr n , q n ) at time t n and the control variables (e n+1 , T n+1 ) at current time t n+1 . The output is the state variables (e tr n+1 , q n+1 ) at time t n+1 , from which the stress can be deduced using Eq. ( 19). The proposed algorithm results from a careful analysis of the incremental energy minimization problem [START_REF] Cheikh | A review of constitutive models and modeling techniques for shape memory alloys[END_REF] and delivers the solution of the time-discretized problem [START_REF] Bo | Thermomechanical modeling of polycrystalline smas under cyclic loading. Part III: evolution of plastic strains and two-way shape memory effect[END_REF] in all cases. The presented algorithm can be readily implemented in a FE code for solving three-dimensional boundary value problems, as will be demonstrated in Section 4.

The overall structure of the algorithm is of the predictor-corrector type. One first checks (through the conditions in Table 1) whether the elastic guess happens to give the solution. If not, the value of (e tr n+1 , q n+1 ) is updated so as to satisfy the optimality condition [START_REF] Bo | Thermomechanical modeling of polycrystalline smas under cyclic loading. Part III: evolution of plastic strains and two-way shape memory effect[END_REF]. That updating procedure proceeds in a two-step fashion by distinguishing between the cases of unsaturated phase transformation and saturated phase transformation. In general, iterative solvers are needed for solving the nonlinear equations (55), (62), (65) that arise in the updating procedure. We emphasize that all those nonlinear equations are scalar and can thus be solved in a very robust fashion by using (for instance) a combination of bisection and Newton methods (see, e.g., ref. [START_REF] Vrahatis | From linear to nonlinear iterative methods[END_REF][START_REF] Mejzlik | A bisection method to find all solutions of a system of nonlinear equations[END_REF]).

We recall that the idea, the structure, and the formalism of elastic/unsaturated/saturated evolution stem from previous works, see, e.g., [START_REF] Sedlák | Thermomechanical model for NiTi-based shape memory alloys including R-phase and material anisotropy under multi-axial loadings[END_REF][START_REF] Stupkiewicz | A robust model of pseudoelasticity in shape memory alloys[END_REF][START_REF] Artioli | An incremental energy minimization state update algorithm for 3D phenomenological internal-variable sma constitutive models based on isotropic flow potentials[END_REF][START_REF] Scalet | A robust and efficient radial return algorithm based on incremental energy minimization for the 3D Souza-Auricchio model for shape memory alloys[END_REF]. The present algorithm can be in fact interpreted as an extension of a radial return algorithm proposed in ref. [START_REF] Scalet | A robust and efficient radial return algorithm based on incremental energy minimization for the 3D Souza-Auricchio model for shape memory alloys[END_REF] for the original Souza-Auricchio model in which e tr is the only internal variable.

The algorithm in ref. [START_REF] Scalet | A robust and efficient radial return algorithm based on incremental energy minimization for the 3D Souza-Auricchio model for shape memory alloys[END_REF] involves a single scalar parameter, whereas more conventional return-mapping algorithms (see, e.g., ( 45)) involve 5-7 parameters (namely, the components of e tr , the plastic multiplier, and a Lagrange multiplier associated with the constraint e tr ≤ ε L ). For the model considered in this paper -which includes the permanent inelastic strain q as an additional internal variable -it can be expected that conventional return-mapping algorithms would involve 10-12 scalar parameters, as reported in ref. [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF]. This may increase computational costs and cause trouble of convergence when using Newton-Raphson procedures. Such difficulties are avoided by the presented algorithm, since it has the distinctive property of involving only scalar nonlinear equations.

Algorithm 1 Pseudocode of the proposed algorithm 1: a ← 2Ge n+1 , u 0 ← e tr n + κ 2 q n 2: if the conditions in Table 1 are verified then (e tr n+1 , q n+1 ) ← (e tr n , q n ) Elastic evolution 4: else

5: if e tr n = q n and G q n -a ≤ R Y √ 1 + κ 2 then 6:
q * ← e tr Calculate y by solving the polynomial equation ( 45)

10:

q * ← ya + R Y u 0 yG + R Y (1 + κ 2 ) 11: if (h -A)q * + R Y y κ 2 (q * -q n ) ≤ τ M then 12:
(ẽ tr , q) ← (q * , q * ) 13:

go to line 18 Calculate y by solving f (y) = 0 in Eq. (55)

17:

Calculate (ẽ tr , q) by Eq. (54)

18:

if ẽtr ≤ ε L then 19:

(e tr n+1 , q n+1 ) ← (ẽ tr , q) Unsaturated phase transformation if ( q * ≤ ε L ) or (e tr n = q n and e tr n = ε L and

√ B ≤ R Y √ 1 + κ 2 ) then 22:
go to line 30 Calculate y by solving g(y) = 0 in Eq. ( 62)

25: q * ← ε L ya + R Y u 0 ya + R Y u 0 26: if (h -A)q * + R Y y κ 2 (q * -q n ) ≤ τ M and a + R Y y u 0 ≥ (G + R Y y (1 + κ 2 ))ε L then 27:
(e tr n+1 , q n+1 ) ← (q * , q * ) 28:

go to line 35 Calculate γ by solving h(γ) = 0 in Eq. (65)

31:

(e tr n+1 , q n+1 ) ← (e tr (γ), q(γ)) Saturated phase transformation

Tangent stiffness operator

The tangent stiffness operator dσ/dε can be obtained in a fully explicit form. Let 1 and I be the second-order and fourth order identity tensors, respectively. Setting J = (1 ⊗ 1)/3 and K = I -J, we have from Eq. ( 7):

dσ dε = C -4G 2 de tr n+1 da K (67)
where C = 3K I + 2G K is the elasticity tensor. In the case of an elastic evolution, we obviously have de tr n+1 /da = 0, so that dσ/dε = C. In the following, we give the expression of de tr n+1 /da = 0 and dσ/dε in the case where phase transformation occurs. In accordance with previous developments, different expressions are obtained depending on the values of (e tr n+1 , q n+1 ). We refer to APPENDIX A for the derivation of the results presented next.

Unsaturated phase transformation ( e tr n+1 < ε L )

Let us first consider the situation where e tr n+1 = q n+1 . As detailed in APPENDIX A, we have :

de tr n+1 da = X K + y X y 2 -y X U 0 2 U 0 ⊗ U 0 (68)
where:

U 0 = (1 + κ 2 )q n+1 -e tr n -κ 2 q n , X = (G + y (1 + κ 2 )) -1 . ( 69
)
We recall that y = (q n+1 -e tr n , q n+1 -q n ) κ and y = R Y /y. Substituting Eq. (68) into Eq. (67) gives:

dσ dε = C -4G 2 X K + y X y 2 -y X U 0 2 U 0 ⊗ U 0 . ( 70
)
The expression of the operator de tr n+1 /da gets more involved in the situation where e tr n+1 = q n+1 . Introducing the tensors U = (h + κ 2 y -A)(e tr n+1 -q n+1 ),

V = (x + h + κ 2 y )(e tr n+1 -e tr n ) + κ 2 (x + A)(q n+1 -q n ), V = (x + A)(e tr n+1 -e tr n ) + κ 2 (2G + x + y )(q n+1 -q n ), (71)
and the scalar quantities

X 11 = y y 2 V : (e tr n+1 -e tr n ) + κ 2 V : (q n+1 -q n ) , X 22 = x (G + (1 + κ 2 )y ), X 12 = √ x y xy (V -V ) : (e tr n+1 -q n+1 ), (72)
it can be calculated (see APPENDIX A) that:

de tr n+1 da = x + h + κ 2 y D K + Y 11 U ⊗ U + Y 22 V ⊗ V + Y 12 (U ⊗ V + V ⊗ U ) (73)
where:

Y 11 = 1 Z x x 2 (1 - X 11 D ) , Y 22 = 1 Z y y 2 (1 - X 22 D ) , Y 12 = 1 Z √ x y xy X 12 D and D = (2G + x + y )(x + h + κ 2 y ) -(x + A) 2 , Z = (D -X 22 )(D -X 11 ) -X 2 12 . (74)
Expression (73) holds if e tr n+1 < ε L and e tr n+1 = q n+1 . In such case, we obtain from Eq. (67) that:

dσ dε = C -4G 2 x + h + κ 2 y D K + Y 11 U ⊗ U + Y 22 V ⊗ V + Y 12 (U ⊗ V + V ⊗ U ) . ( 75
)
3.5.2 Saturated phase transformation ( e tr n+1 = ε L )

In the case e tr n+1 = q n+1 with e tr n+1 = ε L , we obtain:

de tr n+1 da = X P + y X y 2 -y X U ⊥ 0 2 U ⊥ 0 ⊗ U ⊥ 0 ( 76
)
where X = ε L /( a + y (e tr n + κ 2 q n )), P = K -(e tr n+1 ⊗ e tr n+1 )/ε 2 L and U ⊥ 0 = PU 0 . In Eq. ( 76), U 0 is defined as in Eq. ( 69). The operator P can be interpreted as the projector on the orthogonal to e tr n+1 in the deviatoric space. It follows from Eq. ( 76) that:

dσ dε = C -4G 2 X P + y X y 2 -y X U ⊥ 0 2 U ⊥ 0 ⊗ U ⊥ 0 . ( 77
)
In the case e tr n+1 = q n+1 with e tr n+1 = ε L , the operator de tr n+1 /da takes the form:

de tr n+1 da = x + h + κ 2 y D P + Ỹ11 U ⊥ ⊗ U ⊥ + Ỹ22 V ⊥ ⊗ V ⊥ + Ỹ12 (U ⊥ ⊗ V ⊥ + V ⊥ ⊗ U ⊥ ) (78)
where U ⊥ = PU , V ⊥ = PV and (U , V ) are defined as in Eq. ( 71). The scalar Ỹ11 , Ỹ22 , Ỹ12 in (78) are defined by:

Ỹ11 = 1 Z x x 2 (1 - X11 D ) , Ỹ22 = 1 Z y y 2 (1 - X22 D ) , Ỹ12 = 1 Z √ x y xy X12 D
where

D = (2G + γ + x + y )(x + h + κ 2 y ) -(x + A) 2 , Z = ( D -X22 )( D -X11 ) -X2 12 and X11 = y y 2 V ⊥ : (e tr n+1 -e tr n ) + κ 2 W : (q n+1 -q n ) , X22 = x x 2 (G + γ + (κ 2 + 1)y ) Pq n+1 2 + D x + k + κ 2 y (e tr n+1 : (q n+1 -e tr n+1 )) 2 ε 2 L , X12 = √ x y xy (V ⊥ -W ) : (e tr n+1 -q n+1 ). (79)
The tensor W that appears in Eq. ( 79) is given by:

W = PV + D x + h + κ 2 y κ 2 ε 2
L (e tr n+1 : (q n+1 -q n ))e tr n+1 .

From Eqs. ( 67) and (78) we finally obtain:

dσ dε = C -4G 2 x + h + κ 2 y D P + Ỹ11 U ⊥ ⊗ U ⊥ + Ỹ22 V ⊥ ⊗ V ⊥ + Ỹ12 (U ⊥ ⊗ V ⊥ + V ⊥ ⊗ U ⊥ ) . ( 80
)
To close this section, we emphasize that the various expressions obtained for the tangent stiffness operator are fully explicit and do not involve any matrix inversion, thus making for an efficient and accurate numerical evaluation of dσ/dε.

Numerical results

This section presents the results obtained by the numerical implementation of the proposed algorithm in a FE framework.

Table 2: Adopted model parameters taken from ref. [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF].

Parameter We implemented the algorithm within a user-defined material subroutine (UMAT) of the FE software ABAQUS/-Standard [START_REF]Abaqus. Abaqus, analysis user's manual[END_REF]. Following the methodology proposed in ref. [START_REF] Boatti | Development of an effective and user-friendly numerical framework for the simulation of complex smart material components and devices[END_REF], we exploited the capabilities of the package AceGen (50) of the symbolic software Mathematica to generate the UMAT.

In the pseudocode of the algorithm as presented in Sect. 3.4, there appear several equality conditions of the form u = v, where u and v are two scalar quantities. In the numerical implementation, such equality conditions have been replaced by inequality conditions of the form |u -v| ≤ tol where tol is an absolute tolerance parameter set to 10 -6 .

Indeed, round-off errors in practical computations almost always prevent floating-point variables to be exactly equal.

Since the solution (e tr n+1 , q n+1 ) to the minimization problem [START_REF] Cheikh | A review of constitutive models and modeling techniques for shape memory alloys[END_REF] depends continuously on the data (e n+1 , e tr n , q n ), the exact value of the tolerance parameter tol is expected not to impact the robustness of the algorithm.

Numerical experiments aim to verify model implementation as well as to evaluate algorithm performances and robustness. We perform several tests involving both PE and SME, of increasing complexity, starting from simple uniaxial tests to more complex three-dimensional boundary-value problems on real devices.

We consider the three sets of material properties proposed in ref. [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF] and specified in Table 2, in order to verify the main features of the model, compared to the results in ref. [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF], as well as to show the role played by each parameter.

Uniaxial tests

We first simulate several pseudoelastic tests on a single 8-node hexahedral element, under force control and prescribed homogeneous temperature. In particular, we apply a pressure varying cyclically between a maximum and a minimum value at a fixed temperature.

The first pseudoelastic test involves 50 tension cycles with permanent inelasticity. The total analysis time is 10 s and each cycle lasts 0.2 s. We adopt material parameters from Set 1 (see Table 2), where both parameters h and A are equal to zero, while κ is different from zero. The applied pressure varies between 0 and 500 MPa and temperature is set equal to 298 K. Figure 1(a) shows the stress-strain curve, where it can be observed that the non-zero parameter κ gives rise to a permanent inelasticity phenomenon. This is also clear by observing Figures 1(c)-(d), representing the evolution of the internal variables e tr and q in time. The second pseudoelastic test involves 50 tension cycles with saturating permanent inelasticity. The total analysis time is 10 s and each cycle lasts 0.2 s. We adopt material parameters from Set 2 (see Table 2), where A is equal to zero, while h and κ are different from zero. The applied pressure varies between 0 and 500 MPa and temperature is set equal to 298 K. Figure 2(a) shows the stress-strain curve; in such a case, the non-zero parameter h gives rise to the saturation of permanent inelasticity, for which the permanent strain saturates and does not exceed the threshold.

The evolution of the internal variables e tr and q in time reported in Figures 2(c)-(d) is similar to the previous test. The third pseudoelastic test involves 50 tension cycles with saturating permanent inelasticity, including degradation effect. The total analysis time is 10 s and each cycle lasts 0.2 s. We adopt material parameters from Set 3 (see Table 2), where all the parameters are different from zero. The applied pressure varies between 0 and 500 MPa and temperature is set equal to 298 K. Figure 3(a) shows the stress-strain curve; in such a case, the non-zero parameter A couples the two internal variables and results in shifting down the loops. The evolution of the internal variables e tr and q in time, reported in Figures 3(c)-(d), is again similar to previous tests. Figure 3(b) reports the strain rate over time for time steps of 0.001 and 0.01 s, corresponding to 200 and 20 increments per tension cycle.

Set 3 has been used also to perform the fourth pseudoelastic test, simulating 10 tension cycles followed by 10 compression cycles. Figure 4 shows the stress-strain curve, where the applied pressure varies between -500 and 500

MPa.

It is worth highlighting that all the reported results are in perfect agreement with the results reported in ref. [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF]. We then simulate two shape-memory tests on a single 8-node hexahedral element. For all these tests we adopt material parameters from Set 3 (see Table 2).

The first shape-memory test consists of 25 tension cycles at constant low temperature, each one followed by heating strain recovery. The total analysis time is 10 s and each cycle lasts 0.4 s. Initially, a pressure of 500 MPa is applied at a low temperature of 200 K. Then, the temperature is increased up to 400 K at zero applied pressure. Figure 5(a)

shows the stress-strain curve. An inelastic effect is activated, so that we observe only a partial shape recovery (see Figure 5(c)). Also for this test, algorithm robustness is verified (see Figure 5(b)).

The second shape-memory test is performed under force control and constant low temperature, followed by cyclic heating. Initially, a temperature of 200 K is prescribed and a pressure of 250 MPa is applied. Then, the temperature is increased up to 500 K at constant applied pressure. The temperature is then varied cyclically between 200 and 500 K. The total analysis time is 186.6 s and each thermal cycle lasts 3.2 s. We adopt material parameters from Set 3 (see Table 2). Figure 6(a) shows the strain-temperature curve. It is observed that the saturation is reached at the first thermal cycle; therefore, subsequent cycles determine only a shift of the curve due to the presence of permanent deformation (see Figure 6(b)). The test has been repeated by adopting different time steps, respectively, of 0.001 and 0.01 s, corresponding to 320 and 32 increments per thermal cycle. Only the case with time steps of 0.001 s is shown; however, also for this test robustness has been verified.

Pseudoelastic stent strut

The use of SMAs to manufacture cardiovascular stents is increasing, since pseudoelasticity enables the essential selfdeployment behavior. Once implanted in the patient, stents are subjected to million of pulsatile cycles and their fatigue resistance is therefore a topic of wide interesting in the medical, industrial, and scientific communities. From We propose here to investigate the behavior of a representative stent strut subjected to cycling loading under two strain-controlled conditions.

The strut is obtained from the planar CAD representation of a stent geometrically resembling a Bard ViVEXX carotid stent (C.R. Bard Angiomed GmbH & Co., Germany). Figure 7 shows the Bard ViVEXX carotid stent geometry and the adopted mesh of the stent strut, consisting of 7872 eight-node hexahedral finite elements and 11480 nodes. A mesh refinement has been performed to determine the appropriate model.

We adopt material parameters from Set 3 (see Table 2). We apply a tensile cyclic displacement u to one side of the strut, while the other side is fully constrained (see Figure 7). All other surfaces are traction free. A temperature of 310 K is prescribed. The first loading history consists in applying a displacement u of 1 mm and then in varying it between 1 and 0.5 mm for 30 cycles. The second history consists in the application of a displacement of 1.5 mm and then in varying it between 1.25 and 1 mm for 12 cycles. The total analysis time is 25 s. For both the analyses we adopt an adaptive time step varying between 10 -6 and 0.05 s. Figures 8(a)-(b) show the applied displacement u versus the total reaction force of the fixed side plot for the two loading histories, respectively. As observed, the non-zero parameter A couples the two internal variables and results in shifting down the loops, up to a saturation curve. Figures 9(a) and (b) report, respectively, the loading time increments and the total global iterations during the analysis for the second loading history. As it can be observed, small time increments and a high number of total global iterations are required only during the loading phase, which ends at an analysis time of 1 s with a final applied displacement of 1.5 mm. Then, the trend repeats during cycling between 1.25 and 1 mm from 2 s to 25 s.

Helical spring

We conclude this section with the simulation of a helical spring, which can be used as thermal actuator and as element for passive vibration isolation. Therefore, spring components are generally tested in both pseudoelastic and shape-memory regime. The spring under investigation has a wire diameter of 1.0 mm, an external coil diameter of 6.0 mm, a pitch of 2.5 mm, 2 active coils, and an initial length of 5.0 mm. The mesh is composed of 3712 eight-node hexahedral finite elements and 4486 nodes, as shown in Figure 10; a mesh refinement has been performed to choose the appropriate mesh. We adopt material parameters from Set 3 (see Table 2). One end of the spring is fixed (see Figure 10).

For the pseudoelastic tests, we perform multiple tension cycles by applying an axial tensile force F at one section of the spring at a constant temperature of 298 K (see Figure 10). All the nodes of the section where the force is applied are constrained against the two translations in the directions orthogonal to the axial one. The force F varies cyclically between zero and a maximum value. Three maximum values are considered of 20, 40, and 100 N and for each of them we perform 10, 50, and 50 cycles, respectively. The total analysis time is 20 s. All the simulations are performed using variable time step increments between 10 -6 and 0.05 s. Figures 11(a)-(c) show the applied force-axial displacement output for the three pseudoelastic loading histories.

For the shape-memory tests, we perform 19 thermal cycles at constant force F of 5 N. Temperature varies between 208.6 and 298 K. The total analysis time is 40 s. All the simulations are performed using variable time step increments between 10 -6 and 0.5 s. Figure 11(d) reports the axial displacement versus temperature plot.

In both pseudoelastic and shape-memory tests, the non-zero parameter A couples the two internal variables and results in shifting down the loops, up to a saturation curve. Figures 12(c) and (d) report the values for the shape-memory test. As it can be observed, for the pseudoelastic test, small time increments and a high number of total global iterations are required only during the loading phase, which ends at an analysis time of 1 s with a final applied force of 100 N. For the shape-memory tests larger time increments are sufficient for global convergence. For both tests, the trend repeats during cycling. This paper has presented a new algorithmic scheme for a three-dimensional model taking into account permanent inelastic phenomena and degradation effects due to cyclic loading in SMAs. A variational structure of model equations, based on an incremental energy minimization approach, has been presented. The resulting algorithm is free from any kind of regularization, which may affect the structural response, and is simple to implement. The key feature of the algorithm is that it relies only on scalar nonlinear equations, thus avoiding to solve simultaneously for 10-12 scalar unknowns as could be expected in a conventional return-mapping algorithm. Such simplifications are achieved thanks to the incremental energy minimization structure and ultimately trace back to the isotropic nature of the considered constitutive model. Several numerical FE simulations have demonstrated the correctness of the scheme in predicting material response, its ability to solve complex boundary values problems, as well as its robustness.

Although our algorithm is specific to the choice (6) for the dissipation potential, there is no major difficulty in adapting it to the dissipation potentials ( εtr , q) → R Y ( εtr , q) κ,∞ and ( εtr , q) → R Y ( εtr , q) κ,1 mentioned in Sect. 2. In such cases, however, the elasticity domain would have corners so additional subcases are expected to appear in the algorithm. Among other further possible developments, it would also be interesting to extend the present approach to geometrically nonlinear problems.

Financial disclosure

None reported.

A.1.2 Case e tr n+1 = q n+1

In the case of unsaturated phase transformation with e tr n+1 = q n+1 , differentiating the optimality conditions [START_REF]Abaqus. Abaqus, analysis user's manual[END_REF] yields We recall that x = e tr n+1 -q n+1 , y = (e tr n+1 -e tr n , q n+1 -q n ) κ , x = τ M /x and y = R Y /y. From the definition of x, we have by differentiation: dx = (e tr n+1 -q n+1 ) x : (de tr n+1 -dq n+1 ) (A7) which can be rewritten as: where U and V are the tensors introduced in Eq. ( 71), and D is defined by Eq. (74). Substituting in Eq. (A11) and using Eq. (A10) leads to the expression (73) of the de tr n+1 /da, from which the tangent stiffness operator directly follows.

A.2 Saturated phase transformation

In the case of saturated phase transformation, we have e tr The property (A13) will proved to be useful in the following.

A.2.1 Case e tr n+1 = q n+1

Assuming e tr n = q n , differentiating (59) gives: da = (γ + y (1 + κ 2 ))de tr n+1 -y y U 0 dy + e tr n+1 dγ. which correspond to Eq. (76).

A.2.2 Case e tr n+1 = q n+1

In the case of saturated phase transformation with e tr n+1 = q n+1 , the optimality conditions read as

Figure 1 :

 1 Figure 1: First uniaxial pseudoelastic test: 50 tension cycles with permanent inelasticity. (a) Axial stress-axial strain curve obtained for ∆t = 0.001 s. (b) Axial strain-analysis time curve for different time steps ∆t. (c) Axial components of e tr and q versus analysis time plot for ∆t = 0.001 s. (d) Axial components of e tr and q during the first loading cycle.

Figure 1 (

 1 b) reports the strain rate over time for different adopted time steps, respectively, of 0.001 and 0.01 s, corresponding to 200 and 20 increments per tension cycle. The curves demonstrate the robustness of the proposed solution algorithm.

Figure 2 :

 2 Figure 2: Second uniaxial pseudoelastic test: 50 tension cycles with saturating permanent inelasticity. (a) Axial stressaxial strain curve obtained for ∆t = 0.001 s. (b) Axial strain-analysis time curve for different time steps ∆t. (c) Axial components of e tr and q versus analysis time plot for ∆t = 0.001 s. (d) Axial components of e and q during the first loading cycle.

Figure 2 (

 2 Figure 2(b) shows again the robustness of the algorithm, by reporting the strain rate over time for time steps of 0.001 and 0.01 s, corresponding to 200 and 20 increments per tension cycle.

Figure 3 :

 3 Figure 3: Third uniaxial pseudoelastic test: 50 tension cycles with saturating permanent inelasticity, including degradation effect. (a) Axial stress-axial strain curve obtained for ∆t = 0.001 s. (b) Axial strain-analysis time curve for different time steps ∆t. (c) Axial components of e tr and q versus analysis time plot for ∆t = 0.001 s. (d) Axial components of e tr and q plot during the first loading cycle.

Figure 4 :

 4 Figure4: Fourth uniaxial pseudoelastic test: 10 tension cycles followed by 10 compression cycles with saturating permanent inelasticity, including degradation effect. Axial stress-axial strain curve obtained for ∆t = 0.001 s.

Figure 5 :

 5 Figure 5: First uniaxial shape-memory test: 25 tension cycles, each one followed by heating strain recovery, with saturating permanent inelasticity, including degradation effect. (a) Axial stress-axial strain curve obtained for ∆t = 0.001 s. (b) Axial strain-analysis time curve for different time steps ∆t. (c) Axial components of e tr and q during the first loading cycle.

Figure 6 :

 6 Figure 6: Second uniaxial shape-memory test: multiple thermal cycles under constant applied pressure, with saturating permanent inelasticity, including degradation effect. (a) Axial strain-temperature curve obtained for ∆t = 0.001 s. (b) Axial component of e tr and of q versus analysis time plot for ∆t = 0.001 s.

  u

Figure 7 :

 7 Figure 7: SMA cardiovascular stent strut: initial geometry, adopted mesh, and boundary conditions.

Figure 8 :

 8 Figure8: SMA cardiovascular stent strut: applied displacement vs. reaction force of the fixed side for two different loading histories: (a) first a 1 mm displacement is applied and then it varies between 1 and 0.5 mm; (b) first a 1.5 mm displacement is applied and then it varies between 1.25 and 1 mm.

Figure 9 :

 9 Figure 9: SMA cardiovascular stent strut. Analysis of the case where a 1.5 mm displacement is applied and then varied between 1.25 and 1 mm. Trends of (a) the loading time increments and (b) the total global iterations during the analysis.

FixedFFigure 10 :

 10 Figure 10: SMA spring: initial geometry, adopted mesh, and boundary conditions.

Figure 11 :

 11 Figure 11: Pseudoelastic tests on SMA spring: displacement of the free end vs. applied force of (a) 20 N, (b) 40 N, and (c) 100 N. Shape-memory tests on SMA spring: (d) temperature vs. displacement of the free end.

Figure 12 :

 12 Figure 12: SMA spring. Trends of the loading time increments and the total global iterations during the analysis for the (a)-(b) pseudoelastic test with applied force varying between 0 and 100 N and (c)-(d) the shape-memory test.

  x + y )de tr n+1 -(x + A)dq n+1 -(x + A)de tr n+1 (x + h + κ 2 y )dq n+1

2 Ux x 2 U

 22 the definition of y we have by differentiation:* ⊗ U * -y y 2 V * ⊗ V * ) : * ⊗U * -y y 2 V * ⊗V * can be inverted in closed-form. Setting Ũ = M -1 : U * and Ṽ = M -1 : V * ,it can indeed be verified that:M -x x 2 U * ⊗ U * -y y 2 V * ⊗ V * -1 = M -1 2G -y )(e tr n+1 -q n+1 ) x )(e tr n+1 -e tr n ) + κ 2 (2G + x + y )(q n+1 -q n )   

n+1 2 =

 2 ε L hence by differentiation: de tr n+1 : e tr n+1 = 0. (A12) Let P = K -(e tr n+1 ⊗ e tr n+1 )/ε 2 L be the projection on the orthogonal to e tr n+1 in the deviatoric space. Eq. (A12) implies that Pde tr n+1 = e tr n+1 . (A13)

-e tr n κ 2

 2 (q n+1 -q n ) γ + x + y )de tr n+1 -(x + A)dq n+1 -(x + A)de tr n+1 (x + h + κ 2 y )dq n+1 to those used in Section A.1.2, Eq. (A18) can be rewritten as:whereN = Mx x 2 U * ⊗ U * -y y 2 V * ⊗ V * with (U * , V *) are defined as in Eqs. (A8)-(A9). Let P be the linear operator defined by: that PN P = P MPx x 2 PU * ⊗ PU * -y y 2 PV * ⊗ PV * . Relation (A22) can be inverted by a relation similar to (A11), leading to expression (78) for de tr n+1 /da.

Note that our notation κ corresponds to 1/κ in the paper by Auricchio et al.[START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF].

To conclude the discussion, Figures12(a) and (b) report, respectively, the loading time increments and the total global iterations during the analysis for the pseudoelastic test with applied force varying between 0 and 100 N, while

Conflict of interest

The authors declare no potential conflict of interests.

APPENDIX A Derivation of the tangent stiffness operator

This section reports the derivation of the tangent stiffness operator introduced in Section 3.5. Accordingly to previous results and discussion, the derivation is performed by distinguishing between the cases of unsaturated and saturated phase transformation. Let us assume that e tr n = q n . Setting U 0 = (1 + κ 2 )q n+1 -e tr n -κ 2 q n , differentiating [START_REF] Brézis | Opérateurs maximum monotones et semigroupes de contractions dans les espaces de Hilbert[END_REF] gives:

A.1 Unsaturated phase transformation

where we recall that y = R Y /y. The variation dy can be obtained by differentiation of the relation y = (q n+1e tr n , q n+1 -q n ) κ as:

Hence,

Therefore:

The expressions (A3) and (A4) can be verified to remain valid when e tr n = q n .

Note that the value of the Lagrange multiplier γ can be obtained from ( 59) and (61), yielding γ + y (1 + κ 2 ) = ( a + y (e tr n + κ 2 q n ))/ε L . Using the expression (A1) of dy, we obtain:

with X = ε L /( a + y (e tr n + κ 2 q n )). In view of (A13), Eq. (A14) can be rewritten as

Applying P to both sides of Eq. (A15) gives

Eq. (A16) can be inverted in a way similar to (A2), yielding: