Giulia Scalet

Michaël Peigney

A robust and efficient radial return algorithm based on incremental energy minimization for the 3D Souza-Auricchio model for shape memory alloys

Keywords: shape memory alloys, numerical scheme, constitutive modeling, incremental energy minimization, radial return

The present paper focuses on the numerical simulation of quasi-static problems involving shape memory alloy (SMA) structures or components. Phenomenological constitutive models formulated within the continuum thermodynamics with internal variable framework describe phase transformation in a SMA by introducing a suitable set of internal variables, which may be constrained to satisfy a set of inequalities. The numerical treatment of such constraints, together with the presence of non-smooth functions and/or complementary conditions in the model formulation, is not an easy task and strongly influences the numerical convergence, algorithm robustness, and computational times. The aim of this paper is to propose a novel stateupdate procedure for the three-dimensional phenomenological model known as the Souza-Auricchio model. The proposed radial return algorithm, relying on an incremental energy minimization approach, allows for an easy implementation of model equations and internal constraints and avoids the use of regularization parameters for the treatment of non-smooth functions. Several numerical simulations assess the noticeable efficiency, robustness, and performance of the proposed approach, while comparisons with a classical algorithm proposed in the literature show the reduced computational times.

Introduction

Shape memory alloys (SMAs) are smart materials with the ability to recover large strains after thermal cycling or mechanical loading/unloading [START_REF] Lagoudas | Shape Memory Alloys: Modeling and Engineering Applications[END_REF]. Thanks to these unique properties, denoted as shape memory effect and pseudoelasticity, SMA materials are employed in numerous commercial applications, ranging from the biomedical, aerospace, automotive, and earthquake engineering fields to recent technologies as additive manufacturing and active textiles [START_REF] Jani | A review of shape memory alloy research, applications and opportunities[END_REF].

Motivated by the ever-increasing use of SMA-based components and devices, a lot of research has been dedicated to the modeling, simulation, and experimental testing of these alloys.

Several constitutive models have been proposed in the literature for describing the complex behavior of SMAs; generally, they are categorized in macroscopic, microscopic, and micro-macro models. A description of all the modeling approaches is out of the present scope; the reader is referred to [START_REF] Cisse | A review of constitutive models and modeling techniques for shape memory alloys[END_REF] for a recent review on the topic. Among the others, phenomenological macroscopic models formulated within the framework of continuum thermodynamics with internal variables attract large engineering interest. Starting from pioneering models involving only scalar internal variables [START_REF] Tanaka | A thermomechanical description of materials with internal variables in the process of phase transitions[END_REF], such models have actually reached high levels of accuracy and refinement (Auricchio et al., 2014a;[START_REF] Lagoudas | Constitutive model for the numerical analysis of phase transformation in polycrystalline shape memory alloys[END_REF][START_REF] Zaki | A three-dimensional model of the thermomechanical behavior of shape memory alloys[END_REF]. Most of the models are in fact able to describe several mechanisms and effects characterizing SMA behavior, e.g., martensite reorientation [START_REF] Arghavani | A 3-d phenomenological constitutive model for shape memory alloys under multiaxial loadings[END_REF][START_REF] Pan | Multi-axial behavior of shapememory alloys undergoing martensitic reorientation and detwinning[END_REF][START_REF] Popov | A 3d constitutive model for shape memory alloys incorporating pseudoelasticity and detwinning of self-accommodated martensite[END_REF], twinned-detwinned martensites [START_REF] Auricchio | A new flexible 3d macroscopic model for shape memory alloys[END_REF], thermomechanical coupling [START_REF] Auricchio | Gradient structures for the thermomechanics of shape-memory materials[END_REF], R-phase and anisotropy [START_REF] Sedlák | Thermomechanical model for niti-based shape memory alloys including r-phase and material anisotropy under multi-axial loadings[END_REF], asymmetric behavior in tension and compression (Auricchio et al., 2009b), transformation-induced plasticity [START_REF] Hartl | Constitutive modeling and structural analysis considering simultaneous phase transformation and plastic yield in shape memory alloys[END_REF][START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF], viscoplasticity [START_REF] Chemisky | A constitutive model for cyclic actuation of high-temperature shape memory alloys[END_REF], internal damage [START_REF] Hartl | Three-dimensional constitutive model considering transformation-induced damage and resulting fatigue failure in shape memory alloys[END_REF], functional fatigue effects [START_REF] Barrera | Macroscopic modeling of functional fatigue in shape memory alloys[END_REF], microstructure-dependent inelasticity [START_REF] Grandi | A phenomenological model for microstructure-dependent inelasticity in shape-memory alloys[END_REF], and the two-way shape memory effect [START_REF] Lexcellent | The two way shape memory effect of shape memory alloys: an experimental study and a phenomenological model[END_REF].

One of the major concerns associated to constitutive modeling is the need of suitable numerical algorithms to treat non-smooth functions and/or local constraints deriving from the set of nonlinear constitutive equations describing material behavior. The presence of equality/inequality constraints in the evolution problem and/or on internal variables as well as the introduc-tion of non-smooth functions in the model definition, often treated through regularized terms, may prevent numerical convergence and may influence numerical robustness, accuracy, and computational times. Therefore, a robust and accurate numerical implementation is of extreme importance for the design of novel structures [START_REF] Peraza-Hernandez | Design and optimization of a shape memory alloy-based self-folding sheet[END_REF] and for the development of novel technologies [START_REF] Meisel | A procedure for creating actuated joints via embedding shape memory alloys in polyjet 3d printing[END_REF]. The implementation generally consists first in the time-integration of the system of constitutive equations and then in a state-update procedure to derive the updated variables.

Several models are available in a suitable form to perform simulations of complex SMA-based geometries subjected to general thermomechanical loading paths in both finite element [START_REF] Arghavani | A 3-d phenomenological constitutive model for shape memory alloys under multiaxial loadings[END_REF][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] Zaki | Time integration of a model for martensite detwinning and reorientation under non-proportional loading using Lagrange multipliers[END_REF] and isogeometric analysis [START_REF] Auricchio | Innovative and efficient stent flexibility simulations based on isogeometric analysis[END_REF][START_REF] Dhote | 3d coupled thermomechanical phase-field modeling of shape memory alloy dynamics via isogeometric analysis[END_REF] frameworks. The state-update procedures generally adopted to treat SMA constitutive equations are based on return-map schemes, incremental energy minimization approaches, or algorithms for mathematical programming.

A wide class of SMA models implements a return-map-like procedure in an implicit framework [START_REF] Auricchio | A three-dimensional model describing stresstemperature induced solid phase transformations: solution algorithm and boundary value problems[END_REF][START_REF] Qidwai | Numerical implementation of a shape memory alloy thermomechanical constitutive model using return mapping algorithms[END_REF][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][START_REF] Hartl | Constitutive modeling and structural analysis considering simultaneous phase transformation and plastic yield in shape memory alloys[END_REF], even if some works have recently proposed implementations in an explicit environment [START_REF] Stebner | Explicit finite element implementation of an improved three dimensional constitutive model for shape memory alloys[END_REF][START_REF] Scalet | Efficiency and effectiveness of implicit and explicit approaches for the analysis of shape-memory alloy bodies[END_REF]. The work by [START_REF] Popov | A 3d constitutive model for shape memory alloys incorporating pseudoelasticity and detwinning of self-accommodated martensite[END_REF] applies an extension of the closest-point projection algorithm to treat a SMA model incorporating single-and multi-variant martensites. The advantage in using such a procedure is granted by its good numerical performance and well-established numerical properties.

To eliminate the need for a predictor-corrector-type scheme and to omit an active set search, which may become elaborate in the presence of coupled evolution equations and internal constraints, [START_REF] Kiefer | Implementation of numerical integration schemes for the simulation of magnetic sma constitutive response[END_REF] presented two alternative algorithms for the integration of the constitutive equations for magnetic SMAs, namely, the classical return-map scheme and a Fischer-Burmeister-based algorithm. The work demonstrates a greater numerical efficiency of the Fischer-Burmeister-based algorithm, compared to the classical return-mapping. In the context of SMA micromechanical modeling, it is noteworthy to cite the contributions by Bartel andHackl (2009, 2010); [START_REF] Bartel | Thermodynamic and relaxationbased modeling of the interaction between martensitic phase transformations and plasticity[END_REF], who also employed the Fischer-Burmeister complementary function. Recently, Auricchio et al. (2014a) have proposed the use of the Fischer-Burmeister function to treat the inequality constraints on the internal variables and the Kuhn-Tucker complementary inequality conditions deriving from a plasticity-like model formulation. A regularized Fischer-Burmeister function is adopted to treat the non-differentiability of the function when both the arguments are equal to zero.

An alternative approach to return-map procedures is the class of variational methods, which relies on an incremental energy minimization approach. [START_REF] Sedlák | Thermomechanical model for niti-based shape memory alloys including r-phase and material anisotropy under multi-axial loadings[END_REF] applied the Nelder-Mead minimization algorithm to solve the derived problem for SMAs and introduced a regularization energy to assure the fulfillment of constraints on internal variables. [START_REF] Stupkiewicz | A robust model of pseudoelasticity in shape memory alloys[END_REF] presented a pseudoelastic model within the incremental energy minimization framework and proposed an unified augmented Lagrangian treatment of both constitutive constraints and non-smooth dissipation function.

Another approach consists in rewriting the problem as a mathematical programming problem which is solved using general optimization methods. [START_REF] Peigney | Numerical simulation of shape memory alloys structures using interior-point methods[END_REF] applied an incremental variational approach to SMAs and reformulated the incremental problem as a linear complementarity problem. The advantage of the formulation is to solve simultaneously the equilibrium equations and the constitutive laws, taking the internal constraints into account. The obtained formulation leads to a simple and efficient numerical algorithm, solved using interior-point methods.

In the present paper, we focus on the three-dimensional phenomenological model introduced by [START_REF] Souza | Three-dimensional model for solids undergoing stress-induced phase transformations[END_REF], and then treated and generalized in [START_REF] Auricchio | A three-dimensional model describing stresstemperature induced solid phase transformations: solution algorithm and boundary value problems[END_REF]; [START_REF] Evangelista | Phenomenological 3d and 1d consistent models for shape-memory alloy materials[END_REF] (hence the denomination Souza-Auricchio model, in the following). The model has received a large attention thanks to several advantageous features [START_REF] Grandi | The souza-auricchio model for shape-memory alloys[END_REF], such as a simple plasticity-like formulation, few parameters, and a robust numerical implementation, which have allowed its application to the simulation of a wide range of devices (Auricchio et al., 2014b). Several works have focused on the numerical implementation of this model, to treat the evolution of the tensorial internal variable (i.e., the transformation strain), subjected to a saturation constraint. [START_REF] Souza | Three-dimensional model for solids undergoing stress-induced phase transformations[END_REF] adopted a return-map scheme and used a criterion for the nucleation of the product phase, while [START_REF] Auricchio | A three-dimensional model describing stresstemperature induced solid phase transformations: solution algorithm and boundary value problems[END_REF] used a regularized parameter to treat the case of vanishing transformation strain and Jähne (2012) proposed an explicit integration scheme. Recently, [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] have adopted an incremental energy minimization approach for the solution of the constitutive equations, without introducing any regularization, and have shown the effect of the regularized term on material response. The use of a regularization parameter strongly influences the predicted response at the cost of accuracy. Arghavani et al. (2011b) proposed an integration algorithm for a finite-strain extension of the Souza-Auricchio model, based on a logarithmic mapping and a nucleation-completion criterion, while Arghavani et al. (2011a) proposed an improved alternative constitutive model for a finite-strain extension of the Souza-Auricchio model, expressed in terms of only symmetric tensors.

The present paper aims to propose a novel solution algorithm for the Souza-Auricchio model. 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 evolution of the transformation strain in a finite time step incrementally minimizes a convex functional, given by the sum of the elastic energy and the dissipation function. The treatment is then based on a radial-return algorithm and a standard Newton-Raphson scheme is adopted to solve the single scalar equation in both the unsaturated and saturated phase transformation cases. The solution algorithm as well as the initial guess for the resultant nonlinear equation are also discussed. The formulation is here restricted to the framework of infinitesimal strain theory; however, the present approach can be applied to geometric nonlinear problems.

The choice of using this approach is due to its suitable variational structure which facilitates the treatment of internal constraints and allows for an efficient numerical implementation. The advantages of the proposed algorithm are its equation simplicity, easy implementation, the possibility of avoiding regularized terms in both energy/dissipation definition and transformation strain norm, and, overall, the reduced computational times. The actual nonlinear system consists of only one scalar equation for both the unsaturated and phase transformation cases. Further simplifications are possible for low-temperature cases, i.e. when temperature is lower than the equilibrium temperature T * . This results in noticeable low computational times and costs, which are of great important in engineering simulations.

To test the performance, robustness, and efficiency of the proposed algorithm several finite element analyses are presented. The simulations range from classical uniaxial and biaxial tests to more complex multiaxial problems, reproducing the pseudoelastic and shape memory properties. Moreover, the algorithm in [START_REF] Auricchio | A three-dimensional model describing stresstemperature induced solid phase transformations: solution algorithm and boundary value problems[END_REF] is implemented to provide a comparison in terms of model response, computational times, as well as iteration number and to highlight the advantages in using the proposed radial return algorithm.

The paper is organized as follows. Section 2 briefly reviews the continuum formulation of the Souza-Auricchio model, together with the related boundary value problem. Section 3 presents the equations in the time-discrete framework and Section 4 describes the proposed algorithmic scheme. Then, Section 5 presents the results of several numerical simulations. Finally, conclusions are given in Section 6.

Continuum model formulation

Constitutive laws

We briefly reviews the time-continuous formulation of the Souza-Auricchio model. The reader is referred to [START_REF] Souza | Three-dimensional model for solids undergoing stress-induced phase transformations[END_REF] for details.

The assumed control variables are the total strain ε and the absolute temperature T , while the internal variable is the transformation strain e tr . According to experimental evidences [START_REF] Wayman | Engineering Aspects of Shape Memory Alloys[END_REF]), e tr is assumed to be traceless and satisfies the saturation constraint:

e tr ≤ ε L , (1)
where ε L is a material parameter corresponding to the scalar maximum transformation strain reached at the end of the phase transformation during an uniaxial test and • denotes the Euclidean norm.

The Helmoltz free energy function Ψ = Ψ(ε, T, e tr ) is expressed as follows:

Ψ = 1 2 κ θ 2 + G e -e tr 2 + τ M e tr + 1 2 h e tr 2 + I ε L (e tr ) . (2)
Here, θ and e are the volumetric and deviatoric part of ε, respectively; different expressions can be chosen for τ M : in this case we assume τ M = β T -T * , where β is a positive parameter related to the dependence of the critical stress on temperature, T * is a reference temperature, and • indicates the positive part function; κ and G are the bulk and shear modulus, respectively; h defines the phase transformation hardening. The indicator function

I ε L (e tr ) = 0 if e tr ≤ ε L +∞ otherwise (3)
is introduced to satisfy the transformation strain constraint (1).

Following standard arguments [START_REF] Gurtin | An Introduction to Continuum Mechanics[END_REF][START_REF] Coleman | The termodynamics of elastic mmaterial with heat conduction and viscosity[END_REF], the constitutive equations are:

p = ∂Ψ ∂θ = κθ , s = ∂Ψ ∂e = 2G(e -e tr ) , X = - ∂Ψ ∂e tr = s -τ M ∂ e tr -he tr -γ e tr e tr , (4)
where p and s are the volumetric and deviatoric part of σ, respectively, and X is the thermodynamic stress-like quantity associated to e tr . If e tr = 0, the subdifferential ∂ e tr in ( 4) is given by ∂ e tr = e tr e tr .

If e tr = 0, we have

∂ e tr = {τ | tr τ = 0; τ ≤ 1}.
Note that ∂ e tr is not single-valued when e tr = 0. The variable γ in (4) results from the indicator function subdifferential ∂I ε L (e tr ) and it is defined as follows:

γ = 0 if e tr < ε L ≥ 0 if e tr = ε L , (5)
with ∂I ε L (e tr ) = γe tr / e tr . To describe phase transformation and inelasticity evolution, a Mises-type limit function F = F (X) is introduced in the following form:

F = X -R Y , (6)
where R Y is a positive material parameter corresponding to the radius of the elastic domain in the deviatoric stress space. The evolution equation for e tr takes the form:

ėtr = λ ∂F ∂X = λ X X , ( 7
)
where λ is the non-negative consistency parameter. The model is then completed by the classical Kuhn-Tucker conditions: λ ≥ 0, F ≤ 0, λF = 0 . (8)

Boundary value problem

Consider a three-dimensional continuum occupying a domain Ω and subjected to a prescribed loading history. The continuum is assumed to be at mechanical and thermal equilibrium at each time. In such a situation, the temperature field is homogeneous at each time and takes a prescribed value T (t). The stress field σ satisfies the equilibrium equation

div σ + f d (t) = 0 in Ω, σ.n = T d (t) on ∂Ω T , (9)
where f d (t) are body forces in Ω and T d (t) are tractions specified on a part ∂Ω T of the boundary ∂Ω. The strain field ε derives from a kinematically admissible displacement field u, i.e.,

ε = 1 2 (∇u + ∇ T u) in Ω, u = u d (t) on ∂Ω u , (10)
where u d (t) are given displacements on ∂Ω u = ∂Ω -∂Ω T . The functions f d , T d , and u d may also depend on the position x, but that dependence is omitted in the notations. Together, the local constitutive relations (4-8) and the field equations (9-10) define a boundary value problem that governs the evolution of the SMA continuum. In general, a discretization both in space and time is necessary to solve that boundary value problem, as discussed next.

Discrete model formulation

The time-discretization consists in introducing a finite number of time instants t 0 < • • • < t N and estimating the solution at each time t i using a timemarching approach. To do so, the crucial step is to introduce an incremental problem for estimating the values (u n+1 , e tr n+1 ) at t = t n+1 , assuming the values (u n , e tr n ) at t = t n known. A widely used time-discretization scheme is the Euler implicit scheme, which in the present context corresponds the incremental problem defined as follows:

u n+1 ∈ K u (t n+1 ), σ n+1 ∈ K σ (t n+1 ), e tr n+1 ∈ K ε L , (11.1) p n+1 = κθ n+1 , s n+1 = 2G(e n+1 -e tr n+1
) , (11.2)

X n+1 = s n+1 -τ M,n+1 ∂ e tr n+1 -he tr n+1 -γ n+1 e tr n+1 e tr n+1 , (11.3)
e tr n+1 = e tr n + ∆λ

X n+1 X n+1 , (11.4) γ n+1 ≥ 0, e tr n+1 ≤ ε L , γ n+1 ( e tr n+1 -ε L ) = 0, (11.5) ∆λ ≥ 0, X n+1 ≤ R Y , ∆λ(R Y -X n+1 ) = 0. (11.6)
In (11.1), the sets of admissible fields are defined by

K ε L = {e tr | e tr ≤ ε L for all x ∈ Ω}, K u (t) = {u|u = u d (t) on ∂Ω u }, K σ (t) = {σ|div σ + f d (t) = 0 in Ω; σ.n = T d (t) on ∂Ω T }. (12)
The incremental problem (11) being nonlinear, the issues of existence and uniqueness of solutions need to be addressed. In that regard, it is useful to observe that a variational formulation is attached to the incremental problem (11). Consider indeed the functional F(u, e tr ) defined as the sum F(u, e tr ) = F e (u, e tr ) + F d (e tr ) of the elastic potential energy

F e (u, e tr ) = Ω Ψ(ε, T (t n+1 ), e tr )dω - Ω f d .u dω - ∂Ω T T d .u da (13)
and the dissipative term

F d (e tr ) = R Y Ω e tr -e tr n dω. (14)
Using standard results in convex analysis [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], it can be proved (see Appendix A) that the solutions to (11) are minimizers of F, i.e. the incremental problem ( 11) is formally equivalent to the minimization problem inf

(u,e tr )∈Ku(t n+1 )×Kε L F(u, e tr ) ( 15
)
Variational formulations governing the evolution of generalized standard materials have been largely used in the literature in case of convex potentials and have been also extended to non-convex potentials to predict microstructure formation [START_REF] Ortiz | Nonconvex energy minimization and dislocation structures in ductile single crystals[END_REF][START_REF] Miehe | Homogenization of inelastic materials at finite strains based on incremental variational principles. application to the texture analysis of polycrystals[END_REF][START_REF] Carstensen | Non-convex potentials and microstructures in finite-strain plasticity[END_REF]. In the present case observe that the function F is convex in (u, e tr ), with quadratic growth at infinity. It follows that minimizers of F do exist, provided adequate functional spaces are chosen for K u and K ε L . For a three-dimensional continuum, such a choice of appropriate functional spaces is a subtle mathematical issue that is beyond the scope of this paper. The situation simplifies once space discretization is performed: in practice, the spaces K u and K ε L are indeed taken as finite-dimensional, using, for instance, the finite element method. In such case, the aforementioned properties of F ensures that the (space-discretized) incremental problem does admit some solutions. If h > 0 (which is assumed in the following) then F is actually strictly convex, which ensures that the minimizer is unique.

Several strategies are available for solving the minimization problem (15). A possible approach consists in rewriting the minimization problem (15) in the nested form inf u∈Ku(t n+1 )

{ inf

e tr ∈Kε L F(u, e tr )}. ( 16
)
Observing that the variable e tr in ( 16) is free from any differential constraint, we have inf

e tr ∈Kε L F(u, e tr ) = Ω Ψ n+1 (ε) dω - Ω f d .u dω - ∂Ω T T d .u da. ( 17
)
where the function

Ψ n+1 (ε) = inf e tr : e tr ≤ε L Ψ(ε, T (t n+1 ), e tr ) + R Y e tr -e tr n ( 18
)
is obtained by solving a local nonlinear problem. Using (17), Eq. ( 16) can be rewritten as

inf u∈Ku(t n+1 ) Ω Ψ n+1 (ε) dω - Ω f d .u dω - ∂Ω T T d .u da (19)
The functional that appears in ( 19) is similar to the functional F e in ( 13), the only difference being that the Helmoltz free energy Ψ in ( 2) is replaced by the function Ψ n+1 . Solving (19) amounts to solve the equilibirum problem of an elastic nonlinear continuum with a free energy density Ψ n+1 , i.e., obeying the constitutive law σ = Ψ n+1 (ε). Note that Ψ n+1 is path-dependent since it depends on the state e tr n at the previous time instant. From a numerical point of view, Eq. ( 19) is a global nonlinear minimization problem that can be solved using a descend algorithm or a Newton method. Such algorithms make use of the first and second derivatives of the function to minimize, which in the present case amounts to calculate Ψ n+1 (ε) and the consistent tangent operator Ψ n+1 (ε). In the following we propose an efficient algorithm for calculating those quantities.

A radial return algorithm

In the present section, we focus on the local nonlinear problem (18) which, as explained before, is the central building block for solving the incremental structural problem (11). To alleviate the expressions, the term τ M,n+1 is denoted by τ M in the following. Replacing Ψ with its expression (2), the minimization problem in (18) can be rewritten as (24)

From a numerical standpoint, there are two possible approaches for solving (20). The first approach consists in constructing a minimizing sequence, using for instance a linear search algorithm. Since we are in the most favorable situation where the function to minimize is convex, convergence of such algorithms is ensured under mild conditions. The second approach consists in solving the nonlinear equation ( 21), using for instance a Newton-Raphson procedure. Generally speaking, the convergence of the Newton-Raphson procedure is conditioned by the choice of the initial guess. In the case where the initial guess is close enough to the solution, the Newton-Raphson procedure converges very fast -usually faster than linear search algorithms or the like -and it is hard to beat. In the following, we favor the second approach, i.e., solving the optimality equation ( 21). The reason is that Eq. ( 21) can be reduced to a scalar equation, as detailed next. For scalar problems, the issue of finding a good initial guess simplifies considerably, which allows one to take full advantage of the Newton-Raphson procedure.

Check for elastic evolution or fully reversed phase transformation

To solve ( 21), the first step consists in checking whether e tr n+1 = e tr n or e tr n+1 = 0. This is accomplished using the expressions (23-24). The obtained conditions take different expressions whether e tr n = 0 or e tr n = 0: • Case e tr n = 0: We have e tr n+1 = e tr n (elastic evolution) if and only if

a -τ M e tr n e tr n ≤ R Y . (25)
We have e tr n+1 = 0 (fully reversed phase transformation) if and only if

2G e + R Y e tr n e tr n ≤ τ M . (26)
• Case e tr n = 0: We have e tr n+1 = 0 if and only if

a ≤ R Y + τ M . (27)
In (26-27), a is the trial elastic state, defined by

a = 2G (e -e tr n ). ( 28
Eq. ( 31) shows that the tensor e tr n+1 -e tr n is parallel to b. Taking the norm of (31) and using the fact that G , τ M , x, R Y are positive, we find

e tr n+1 -e tr n = b -R Y τ M /x + 2G
.

Substituting in (31), we finally arrive at

e tr n+1 = e tr n + b b b -R Y τ M /x + 2G . ( 33
)
Eq. ( 33) shows that e tr n+1 is entirely determined by the positive scalar x. We are thus left with the issue of finding x. This is accomplished by taking the norm of (33), leading to the consistency equation

x = e tr n + b b b -R Y τ M /x + 2G (34)
in which x is the only unknown. Eq. ( 34) can be rewritten in the form

f (x) = 0 where f (x) = c -τ M -2G x (35)
and

c = 2G e -R Y b b . ( 36
)
The equation f (x) = 0 has a unique solution in ]0, +∞[. More precisely, we have lim

x→0 + f (x) > 0; f (x max ) ≤ 0 (37)
where Using Eq. ( 37) and the continuity of f on [0, +∞[, a direct application of the intermediate value theorem shows that the solution of f (x) = 0 lies in the bounded interval ]0, x max ]. This is a useful information for solving the equation f (x) = 0 numerically. In that regard, another useful ingredient is the derivative f (x). A direct calculation gives: This case corresponds to saturated phase transformation. Provided τ M = 0, we can use a change of variable to put (38) in the same form as (30). Introducing indeed the variable x defined by

x max = e + (R Y -τ M )/2G .
f (x) = - τ M x 2
τ M x = τ M ε L + γ n+1 ,
we can see that Eq. ( 38) is formally identical to Eq. ( 30), except that x is now submitted to the constraint 0 ≤ x ≤ ε L . Eq. ( 33) still holds, but the consistency equation for determining x now reads as

ε L = e tr n + b b b -R Y τ M /x + 2G . ( 39
)
Eq. ( 39) can be rewritten as:

g(x) = 0 where g(x) = c -ε L τ M x -2G ε L ( 40
)
and c is defined as in (36). The equation g(x) = 0 has a unique solution in ]0, ε L ]. We have:

g (x) = - τ M x 2 R Y b c c : e tr n - (b : e tr n ) b 2 b + ε L τ M x 2 .
Remark: The change of variable detailed above can be used only if τ M > 0. The case τ M = 0 can be treated by replacing τ M /x with γ n+1 in the previous developments. More precisely, starting from (38) and performing manipulations similar to those used in Sect. 4.2, we obtain, in the case τ M = 0,

e tr n+1 = e tr n + b b b -R Y γ n+1 + 2G ( 41
)
where b = a -γ n+1 e tr n . Enforcing that e tr n+1 = ε L in (41) leads to the consistency equation

h(γ n+1 ) = 0 where h(γ n+1 ) = c -ε L γ n+1 -2G ε L and c = 2G e -R Y b / b .

Summary

Algorithm 1 summarizes the proposed solution procedure. The input variables are the transformation strain e tr n at the previous time instant t n and the deviatoric strain e = Ge/(G + h/2) at the current time instant t n+1 . The output variable is the transformation strain e tr n+1 at t n+1 . The stress σ n+1 at time t n+1 is given by

σ n+1 = C : (ε -e tr n+1 )
where C is the elasticity tensor. Let 1 and I be the second-order and fourthorder identity tensors, respectively. Setting J = (1 ⊗ 1)/3 and K = I -J, we note for later reference that

C = 3κJ + 2GK, σ n+1 = C : ε -2Ge tr n+1 ( 42
)
where the second equality results from the fact that e tr n+1 is traceless. Algorithm 1 follows a two-step procedure that is typical of return mapping strategies and similar to the elastic prediction -plastic correction decomposition used in plasticity [START_REF] Simo | Computational Inelasticity[END_REF].

Step 1 consists in checking for special cases of evolutions -namely elastic evolution or fully reversed phase transformation. This is accomplished by checking whether the input data (e tr n , e ) verify certain explicit inequalities. If those inequalities are not satisfied, then the algorithm moves to Step 2 in which a nonlinar problem is solved.

Such a two-step structure is commonly found in the algorithms for the Souza-Auricchio model, discussed in Section 1. A distinctive feature of Algorithm 1 is that the nonlinear problem involves a single scalar parameter, whereas the other mentioned algorithms use 5, 6, or 7 scalar parameters. An expected outcome is a significant gain in computational costs. It is worth highlighting that the computational cost for the factorization is O(n 3 /3), where n is the matrix size [START_REF] Quarteroni | Numerical Models for Differential Problems[END_REF], and the computation of the jacobian for a system of equations increases the computational time. Another related benefit is a gain in robustness. As mentioned earlier, the Newton-Raphson procedure -which is often used for solving nonlinear problemsis known to perform very fast provided the initial guess is close enough to the solution. Failing to have a good initial guess may lead to difficulties in convergence. The issue of finding such a good initial guess simplifies considerably when the problem depends on a single scalar variable; for instance, a bisection method can be used, which consists only in iterative function evaluations that do not increase excessively the computational time.

Since it involves a single scalar unknown in Step 2, Algorithm 1 is similar in spirit to the radial return algorithm originally proposed by [START_REF] Wilkins | Calculation of elastic-plastic flow[END_REF] in perfect plasticity. This is especially clear when τ M = 0, i.e., at low temperature. In that case, there is a simple geometric interpretation of Step 2 in Algorithm 1 (see Fig. 1): Starting from the trial elastic state a = 2G (e -e tr n ), a correction along the 'radial direction' a such that 2G (e -e tr n+1 ) = R Y (in the case of unsaturated phase transformation) or e tr n+1 = ε L (in the case of saturated phase transformation). The situation is more complicated when τ M > 0. For simplicity we only discuss the case of unsaturated phase transformation. Note from (26) that the ball with center τ M e tr n / e tr n and radius R Y defines the elastic domain at time t n . If the trial elastic state a is outside of that domain, a correction is performed in such fashion that a n+1 = 2G (e -e tr n+1 ) falls on the boundary of the elasticity domain at time t n+1 , i.e., satisfies a n+1 -τ M e tr n+1 / e tr n+1 = R Y . When τ M = 0, the correction a n+1 -a is no longer parallel to a, as can be seen from Eq. ( 32). The situation is represented in Fig. 2. The increment e tr n+1 -e tr n is parallel to the correction a n+1 -a. Also note that a, a n+1 are aligned with the center τ M e tr n+1 / e tr n+1 of the elasticity domain at time t n+1 . All those geometric properties are direct consequences of Eq. (30).

Input ;

if e tr n+1 > ε L then Saturated phase transformation if τ M > 0 then Solve g(x) = 0; γ n+1 ← τ M /x; else Solve h(γ n+1 ) = 0; end b ← a -γ n+1

Tangent stiffness operator

Using the expressions given in Sects. 4.1-4.3, the tangent stiffness operator dσ n+1 /dε can be calculated explicitly. In the structure of Algorithm 1 there appear 4 distinct regimes: elastic evolution (e tr n+1 = e tr n ), fully reversed phase transformation (e tr n+1 = 0), unsaturated phase transformation (0 < e tr n+1 < ε L with e tr n+1 = e tr n ) and saturated phase transformation ( e tr n+1 = ε L with e tr n+1 = e tr n ). The operator dσ n+1 /dε takes a different expressions in each one of those regimes. In the case e tr n+1 = e tr n or e tr n+1 = 0, we have dσ n+1 /dε = C, where C is the elasticity tensor of the material, see Eq. ( 42). The expressions get more involved in the remaining cases. More precisely, we have • in the case 0 < e tr n+1 < ε L (with e tr n+1 = e tr n ):

dσ n+1 dε = C -4G 2 x( b -R Y ) b c K + X 11 ẽtr n+1 ⊗ ẽtr n+1 + X 22 b ⊗ b +X 12 (ẽ tr n+1 ⊗ b + b ⊗ ẽtr n+1 ) (43)
• in the case e tr n+1 = ε L (with e tr n+1 = e tr n ):

dσ n+1 dε = C-4G 2 ε L ( b -R Y ) b c (K-ẽ tr n+1 ⊗ẽ tr n+1 + R Y b -R Y (1 -α 2 ) b ⊥ ⊗b ⊥ ) (44) R Y • 0 • a a n+1 •τ M
X 11 = ( b -R Y )τ M X , X 22 = R Y X ( b c b -R Y -τ M ), X 12 = αR Y τ M X , (45)
where

X = b ( c -τ M ) + (1 -α 2 )R Y τ M .
The derivation of Eqs. (43-44) is detailed in Appendix B. We emphasize the fact that expressions (43-44) are fully explicit and do not involve any matrix inversion, thus making for an efficient numerical computation of the tangent stiffness operator.

Numerical simulations

This section presents several numerical simulations to show the performance of the proposed radial return algorithm. The conducted finite element analyses include simple uniaxial and biaxial tests as well as two complex simulations, reproducing both pseudoelasticity and shape memory effect. The algorithm proposed in [START_REF] Auricchio | A three-dimensional model describing stresstemperature induced solid phase transformations: solution algorithm and boundary value problems[END_REF] (in the following, denoted as classical for simplicity) has been implemented to provide a comparison. The adopted material parameters are reported in Table 1.

The proposed algorithm has been implemented within a user-defined material subroutine (UMAT) of the finite element software ABAQUS/Standard (Abaqus, 2010). We use the package AceGen [START_REF] Korelc | Multi-language and multi-environment generation of nonlinear finite element codes[END_REF] of the symbolic software Mathematica to generate the UMAT, according to the methodology proposed in [START_REF] Boatti | Development of an effective and user-friendly numerical framework for the simulation of complex smart material components and devices[END_REF].

Table 1: Material parameters adopted in all the numerical simulations (Auricchio et al., 2009a). Regularization δ is used only in the classical algorithm.

Symbol Value Unit

E 53,000 MPa ν 0.33 - h 1,000 MPa ε L 0.056 - T * 243 K β 6.1 MPa/K R Y 100 MPa δ 10 -6 -

Isothermal tests

We start by considering several three-dimensional isothermal tests. First, we simulate uniaxial tension-compression tests on a single 8-node hexahedral element, under displacement control and prescribed homogeneous temperature. Figure 3(a) shows the applied boundary conditions. In particular, a displacement is applied and increased up to 0.07 mm and three temperature values of 350,255,and 200 K are considered. Figures 4(a), 4(c), and 4(e) reports the stress-strain curves predicted by the proposed algorithm for time steps, respectively, of 0.002, 0.02, and 0.2 s, corresponding to 1000, 100, and 10 increments. The curves demonstrate A qualitative agreement between the two algorithms is evident; however, the predicted responses differ at the beginning of the forward phase transformation and at the end of the reverse phase transformation. This is due to the presence of the regularized parameter δ in the classical algorithm, which affects model response, especially when the norm of the transformation strain e tr tends to zero. As a further comparison, we consider the uniaxial tensile test at constant temperature of 350 K and time step of 0.02 s. Figures 5(a) and (b) report the number of iterations required by the (local) Newton-Raphson scheme in the unsaturated and saturated cases (see Sections 4.2 and 4.3) versus the controlled strain component 33 , for the two algorithms during loading and unloading, respectively. The local iteration number is the same for all the integration points. The values 33 = 6.72 • 10 -2 and 33 = 7.00 • 10 -2 correspond to the saturated case for both the algorithms. Moreover, it is noticed that the classical algorithm starts the phase transformation earlier during loading (i.e., at 33 = 2.80 • 10 -3 ) and ends the phase transformation later during unloading (i.e., again at 33 = 2.80 • 10 -3 ). This increases the computational time for the classical algorithm, since it has to solve a nonlinear system of equations, while the proposed algorithm has to perform an elastic check and the iteration number is equal to zero in Figure 5. It can be noticed that the number of iterations is higher for the classical algorithm in the unsaturated case, thus resulting in longer computational times.

We now consider biaxial tests on a single 8-node hexahedral element, under displacement control and prescribed homogeneous temperature. Figure 3(b) shows the applied boundary conditions. The loading consists of a butterfly-shaped history at a constant temperature of 400 K. A maximum displacement of 0.1 mm is applied, as shown in Figure 6. To demonstrate the robustness of the proposed solution algorithm, Figure 7(a) reports the curves in terms of the non-zero stress components predicted by the proposed algorithm for time steps, respectively, of 0.004, 0.04, and 0.4 s, corresponding to 1000, 100, and 10 increments. Figure 7(b) shows a good matching between the curves obtained with the proposed and the classical algorithm. Differences between the two algorithms are clearer in Figure 8, where the transformation strain norm e tr is plotted. The classical algorithm does not reach zero values of the norm under zero stress and strain states (see time instants of 2 and 4 s).

Finally, we conclude this subsection with the simulation of a square plate with a circular hole, subjected to displacement-control loading at a constant high temperature. The plate has an edge length of 100 mm, a thickness of 2 mm, and a central circular hole of radius 10 mm. According to the symmetry of the problem, we model only one quarter of the plate by applying appropriate boundary conditions. The mesh is composed of 8,712 8-node hexahedral elements and 13,467 nodes, as shown in Figure 9; a mesh refinement has been performed to choose the appropriate mesh. A displacement u Y is applied at the top side of the plate at a temperature of 400 K (see Figure 9). A maximum value of 5 mm is applied and then set back to zero, as reported in Figure 10(a). The simulation is performed using variable time steps between 10 -6 and 0.1 s.

The deformed mesh under the maximum displacement is reported in Figure 10(b). Figure 11 shows the reaction force of the bottom side of the plate versus the applied displacement diagram; a comparison with the results obtained with the classical algorithm is also provided. As it can be observed, the pseudoelastic response of the plate is correctly reproduced by both the algorithms, but differences are evident mostly at the beginning of forward transformation and at the end of reverse transformation, i.e., around displacements of 0.2 and 1.3 mm. This is clear from Figure 12 which represents the contour plot of the norm of the transformation strain e tr at different time instants. Main discrepancies in the distribution are evident at the beginning of the analysis at 0.125 s and corresponding displacement of 1.25 mm.

To assess the computational efficiency of the proposed algorithm, Figures 13(a) and 13(b) report, respectively, the loading time increments and the total global iterations during the analysis. The results are compared to those obtained with the classical algorithm. As it can be observed, the classical algorithm requires smaller time increments to converge, especially in correspondence of phase transformation. On the contrary, the proposed algorithm needs a larger number of iterations; this can be due to the use of larger time increments. Table 2 reports the total increments to complete the analysis for the two algorithms; it is worth to highlight the lower number of time increments for the proposed solution scheme. As a further comparison, we also report the wall-clock times necessary to complete the analysis in Table 3. The values clearly show the lower times required by the proposed algorithm.

Thermal-cycling tests

We now consider several three-dimensional thermal-cycling tests.

First we reproduce the shape memory effect through a uniaxial test on a single 8-node hexahedral element under force-control and constant low temperature, followed by a heating cycle. The applied boundary conditions are those considered for the tension-compression uniaxial tests (see Figure 3(a)), except that a force is now applied on the top surface of the element. Initially, a temperature of 200 K is prescribed and a force of 400 N is applied and then set back to zero. Then, the temperature is increased up to 400 K at zero stress. Figure 14 reports the stress-strain-temperature curves predicted

P ROP OS E D

A n a l y s i s t i me compares the curves obtained with the proposed and the classical algorithm.

The differences are evident at around 260 K, when phase transformation starts.

We now simulate a uniaxial test on a 8-node hexahedral element under force-control and constant high temperature, followed by a thermal cycle. Initially, a temperature of 303 K is prescribed and a force of 300 N is applied and kept constant. Then, the temperature is decreased up to 121.2 K and finally increased up to 454.5 K. Figure 16(a) reports the strain-temperature curves predicted by the proposed algorithm for time steps, respectively, of 0.001, 0.01, and 0.1 s, corresponding to 17000, 1700, and 170 increments. Figures 16(b) compares the curves obtained with the proposed and the classical algorithm.

Finally, we conclude with the simulation of a helical spring actuator, subjected first to a tensile force at constant high temperature and then to a thermal cycling at constant applied force. The spring 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 8node hexahedral elements and 4486 nodes, as shown in Figure17(a); a mesh refinement has been performed to choose the appropriate mesh. An axial tensile force F of 20 N is applied at one section of the spring at a temperature of 423 K, while the other section is fixed (see Figure 17(a)). 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. A thermal cycling is applied at constant force F , by cooling the spring up to 169.2 K and then heating up to the initial temperature (see Figure 18(a)). The simulation is performed using variable time step increments between 10 -6 and 0.01 s. The deformed mesh at 169.2 K is reported in Figure 17(b), while Figure 12 represents the related contour plot of the norm of the transformation strain e tr .

Figure 18(b) shows the displacement of the free end of the spring versus temperature; a comparison with the results obtained with the classical algorithm is also provided. As it can be observed, the two models present differences at the beginning of forward transformation and at the end of reverse transformation.

To assess the computational efficiency of the proposed algorithm, Figures 20(a) and 20(b) report, respectively, the loading time increments and the total global iterations during the analysis. The results are compared to those obtained with the classical algorithm. As it can be observed, the classical algorithm requires smaller time increments to converge, especially in correspondence of phase transformation. On the contrary, the proposed algorithm proceeds with a larger number of iterations; this is attributed to the use of larger time increments. Table 2 reports the total increments to complete the analysis for the two formulations; it is worth to highlight the lower number of time increments for the proposed solution scheme. To provide a further comparison, we also report the wall-clock time necessary to complete the analysis in Table 3, which is considerably lower (about 1/4) for the proposed algorithm.

Conclusions

This paper has presented a new algorithmic scheme for the Souza-Auricchio model. The main idea relies in the exploitation of the variational structure of model equations, based on an incremental energy minimization approach. This results in a nonsmooth optimization problem which, using convex analysis, can be solved without any kind of regularization. This is a welcome feature in practice: as illustrated in Section 5, adding a regularization biases the structural response, especially at the onset of phase transformation.

The proposed approach leads to a radial return scheme that is simple to implement. That scheme has the key feature of being based on the solution of a single scalar equation. The presented approach has been verified to be efficient in practice through several numerical finite element simulations, leading to low computational costs. Such an advantage can be exploited for the computer-based design of complex SMA-based devices.

We note that the proposed approach can be adapted to SMA constitutive models involving tensorial internal variables and constraints. On a related note, the easy implementation offers the possibility of extending it to more complex models, e.g., [START_REF] Auricchio | A three-dimensional model describing stress-induces solid phase transformation with residual plasticity[END_REF].

In this paper, thermal equilibrium was assumed to hold at each time, so that the temperature was acting as a given external parameter. That assumption is valid for small loading rates. When the loading rate increases, thermomechanical effects -most notably through latent heat and thermal diffusion effects -may become significant. Handling such effects can be achieved by combining the proposed approach with incremental variational principles proposed by [START_REF] Peigney | An incremental variational approach to coupled thermo-mechanical problems in anelastic solids. Application to shapememory alloys[END_REF] for coupled thermomechanical problems. where A is the symmetric fourth-order tensor defined as: The tensor A can be inverted in closed form to give:

A = I - R Y b b ⊥ ⊗ b ⊥ .
A -1 = I + Y b ⊥ ⊗ b ⊥

  inf e tr G e tr 2 -2G e tr : e + τ M e tr + R Y e tre tr n + I ε L (e tr ) (20) where G = G + h/2 and e = (G/G )e. The function to minimize in (20) is strictly convex in e tr and therefore admits a unique solution. That solution, denoted by e tr n+1 in the following, is characterized by the optimality equation 0 ∈ 2G (e tr n+1 -e ) + τ M ∂ e tr n+1 + R Y ∂ e tr n+1 -e tr n + γ n+1 γ n+1 ≥ 0, e tr n+1 ≤ ε L , γ n+1 ( e tr n+1 -ε L ) = 0, (22) ∂ e tr n+1 -e tr n = e tr n+1 -e tr n e tr n+1 -e tr n if e tr n+1 = e tr n {τ : tr τ = 0, τ ≤ 1} if e tr n+1 = e tr n , {τ : tr τ = 0, τ ≤ 1} if e tr n+1 = 0.

  ) 4.2. Unsaturated phase transformation If none of the inequalities (25-27) is satisfied, then e tr n+1 / ∈ {0, e tr n } and Eq. (21) becomes 0 = 2G (e tr n+1 -e ) + τ M x = e tr n+1 and γ n+1 satisfying (22). To solve (29), first consider the equation 0 = 2G (e tr n+1 -e ) + τ M drop the constraint e tr ≤ ε L . Eq. (30) can written as 0 = 2G (e tr n+1 -e tr n ) -a + τ M

  In the case τ M = 0, the equation f (x) = 0 can be solved directly. From Eqs. (33-34) we have indeed x = e tr n + phase transformation If the value e tr n+1 provided by (33-34) satisfies the inequality e tr n+1 ≤ ε L , then the solution to problem (20) is obtained. Otherwise, we need to solve 0 = 2G (e tr n+1 -e ) + τ M γ n+1 ≥ 0, e tr n+1 = ε L .

Figure 1 :

 1 Figure 1: Radial return for τ M = 0, case of unsaturated phase transformation.

Figure 2 :

 2 Figure 2: Radial return for τ M > 0, case of unsaturated phase transformation.

Figure 3 :

 3 Figure 3: Applied boundary conditions for the (a) uniaxial tension-compression and (b) biaxial tests under displacement-control and constant prescribed temperature.

Figure 4 :

 4 Figure 4: Uniaxial tests under displacement-control and constant prescribed temperatures of 350, 255, and 200 K. (a)-(c)-(e) Stress-strain curves obtained with the proposed algorithm for different time steps ∆t. (b)-(d)-(f) Comparison between the proposed and the classical algorithm.

Figure 5 :

 5 Figure 5: Uniaxial test under displacement-control and constant prescribed temperatures of 350 K. Comparison between the number of iterations required by the (local) Newton-Raphson scheme at the integration point level in the unsaturated and saturated case for the classical and proposed algorithm. (a) Loading and (b) unloading cases.

Figure 6 :Figure 7 :

 67 Figure 6: Biaxial test under displacement-control and constant prescribed temperature of 400 K. The applied butterfly-shaped history is represented.

Figure 8 :

 8 Figure 8: Biaxial test under displacement-control and constant prescribed temperature of 400 K. Trend of the norm of the transformation strain e tr obtained with the proposed and classical algorithm.

Figure 9 :

 9 Figure9: Pseudoelastic test on a square plate with a circular hole. Adopted mesh and applied boundary conditions.

Figure 10 :Figure 11 :

 1011 Figure 10: Pseudoelastic test on a square plate with a circular hole. (a) Loading history in terms of the applied displacement u Y and (b) deformed mesh at the maximum applied displacement.

Figure 12 :Figure 13 :

 1213 Figure 12: Pseudoelastic test on a square plate with a circular hole. Contour plot of the norm of the transformation strain e tr at three different time instants. Comparison between the proposed and the classical algorithm.

Figure 15 :

 15 Figure 15: Uniaxial test under force-control and constant low temperature, followed by a heating cycle. Comparison between the proposed and the classical algorithm. (a) Stressstrain and (b) strain-temperature curves.

FFFigure 17 :

 17 Figure 17: Thermal cycling constant tensile axial force a spring. (a) Adopted mesh and applied boundary conditions. (b) Deformed spring after cooling at 169.2 K and constant force of 20 N.

Figure 18 :

 18 Figure 18: Thermal cycling at constant tensile axial force on a helical spring. (a) Loading history. (b) Curves in terms of displacement of the free end of the spring and prescribed temperature. Comparison between the proposed and the classical algorithm.

Figure 19 :Figure 20 :

 1920 Figure 19: Contour plot of the norm of the transformation strain e tr at 169.2 K, obtained with the proposed and classical algorithm for the thermal cycling at constant tensile axial force on the helical spring.

  The first inequality in (37) follows from the assumption that none of the inequalities (25-27) are satisfied. To obtain the second inequality in (37), observe from (36) that c ≤ 2G e + R Y .

Substituting in (35) gives f (x) ≤ 2G ( e -x) + R Y -τ M and the result follows.

Table 2 :

 2 Total time increments to complete the analysis. Comparison between the proposed and the classical algorithm.

	Analysis	Proposed Classical
	Pseudoelastic plate with hole 128	174
	Spring actuator	230	812

Table 3 :

 3 Wall-clock times [s] to complete the performed analyses with one CPU. Comparison between the proposed and the classical algorithm.

	Analysis	Proposed Classical
	Pseudoelastic plate with hole 2781	3923
	Spring actuator	3507	15053
	by the proposed algorithm for time steps, respectively, of 0.003, 0.03, and 0.3
	s, corresponding to 1000, 100, and 10 increments. Figures 15(a) and 15(b)

  Figure14: Uniaxial test under force-control and constant low temperature, followed by a heating cycle. Stress-strain-temperature curves obtained with the proposed algorithm for different time steps ∆t.

						400
				33 [MPa]	200 300		t = 0.003 s
								t = 0.03 s
				Stress		0 100		t = 0.3 s
						200
				250 Temperature T [K] 300	350	400	0	0.01	0.02 Strain 33 [-] 0.04 0.05 0.03	0.06
		500
			Proposed
			Classical
		400
	33 [MPa]	200 300
	Stress	100
		0
	0 -100	0.01	0.02	0.03	0.04	0.05	0.06
				Strain 33 [-]

Preprint submitted to European Journal of Mechanics -A/SolidsDecember 20, 2016

Acknowledgment

The authors would like to acknowledge Prof. F. Auricchio and Prof. A. Constantinescu for their advice. This work is funded by Fondazione Cariplo, grant number 2013-1779 and by the French National Research Agency (project Fast3D-ANR-11-BS09-012-01).

Appendix A. Variational formulation of the Euler implicit scheme Let (u n+1 , e tr n+1 ) be a solution to the minimization problem (15). We show in the following that (u n+1 , e tr n+1 ) necessarily satisfies the set of equations (11) that define the Euler implicit scheme. A solution (u n+1 , e tr n+1 ) to the minimization problem (15) necessarily satisfies the optimality condition

where ∂F is the subdifferential of F [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]. Observe that F is differentiable with respect to u. Eq. (A.1) can be rewritten as 0 = ∂ u F(u n+1 , e tr n+1 ), (A.2)

where ∂ u F is the derivative of the F with respect to u, and ∂ e tr F is the subdifferential of F with respect to e tr . Reproducing the classical reasoning used in linear elasticity, Eq. (A.2) leads to Eqs. (11.1-2). Eq. (A.3) needs to be treated more carefully because F is only subdifferentiable with respect to e tr . Since e tr is free from differential constraints, first observe that (A.3) gives the local equation

where Φ(e tr ) = R Y e tre tr n . Further note from (2) that Ψ can be written as Ψ = Ψ 0 + Ψ 1 + I ε L with Ψ 0 = (1/2)κ θ 2 + G ee tr 2 + (1/2) h e tr 2 and Ψ 1 = τ M e tr . Since all the functions Φ, Ψ 0 , Ψ 1 , I ε L are proper, lower-semicontinuous and convex, we have [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])

The expressions of ∂ e tr Ψ 1 and ∂ e tr I ε L have been given in Sect. 2.1. We have Using (A.5), Eq. (A.6) can be seen to be equivalent to (11.4)(11.6). This completes the proof that the stationarity conditions (A.1) coincide with the local equations ( 11) of the Euler implicit scheme.

Appendix B. Derivation of the tangent stiffness operator

From (42) we have

.

Replacing in (B.9) and recalling that b ⊥ : ẽtr n+1 = 0, we find:

Substituting in (B.2) gives the expression (44) for the tangent stiffness operator.


Zhenlian

Ma Peng

Yewen Zhang
email: yewen.zhang@tongji.edu.cn

Stéphane Holé
email: stephane.hole@tongji.edu.cn

Feihu Zheng

Min Gu

Zhenlian An

Analysis and interpretation of photo-stimulated discharge spectrum for polypropylene films under different electric and geometrical conditions

Keywords: photo-stimulated discharge, space charge, trap depth, polypropylene, photoemission

Analysis and interpretation of photo-stimulated discharge spectrum for polypropylene films under different electric and geometrical conditions

Introduction

Space charge accumulation affects the electrical properties of the polymeric insulation by distorting the electric field in the bulk. This can produce degradations and even premature breakdown for high-voltage applications [1][2][3][4][5]. Therefore, the suppression of charge buildup in insulating polymers has been a research focus for decades [6,7]. Modification of charge traps has been demonstrated to be an effective way to modulate the space charge formation and transportation in insulating polymeric materials. It was suggested that shallow and deep traps are formed respectively by the existence of physical defects (e.g. lamella-amorphous interface) and chemical defects (e.g. additives and byproducts) in polymers [8]. Though many kinds of techniques, e.g. thermal methods and optical methods, have been developed to measure the trap depths, experimental results remain scattered [9][10][11][12].

Among those methods, photo-stimulated discharge (PSD) method is a direct and powerful technique to obtain the charge trap depths in polymeric materials [10,[13][14][15][START_REF] He | ICPADM'09[END_REF]. Comparing to thermally stimulated discharge (TSD) method, PSD method detects much deeper traps without destroying the samples. In recent years, important advances have been made in understanding charge traps by PSD method. For example, Gonzalez et al demonstrated that only charges trapped in the deepest traps are responsible for the piezoelectric behavior by selective discharge of cellular polypropylene films [START_REF] Mellinger | [END_REF]. Yue et al investigated the evolution of deep charge traps in polyimide under low-energy proton irradiation and found that the decrease of PSD intensities is due to the irradiation-induced degradation of carbonyl groups [18]. Based on PSD theory about release of charges from the mono-and multi-level of traps, Zhu et al proposed two modified data processing methods for more accurately obtaining the information of trapped charges [19].

However, PSD spectrum is induced not only by the de-trapping process of injected space charges but also by the photoemission from electrodes and polymer materials in practical experiments [START_REF] Tetsuji | 6 th ISE[END_REF][START_REF] Murata | [END_REF]. The photoemission process from the electrode is influenced by both the work function of the electrodes and the surface states of the polymers [22][23][24]. Photoemission from polymeric materials in UV wavelength is from the surface states of the sample [START_REF] Murata | [END_REF][25][26][27]. Because the work function of the illuminated electrode and the dielectric surface states are usually in the same photon energy range in PSD experiments, a comprehensive understanding of the photoemission mechanism from the illuminated electrode and dielectric surface is of great significance for correctly interpreting PSD spectra. However, up to now, photoemission mechanisms have seldom been studied in PSD measurements. When Takai et al. measured the PSD spectra for studying traps in -irradiated polyethylene, they observed a broad peak around 240 nm for which no explanations were given [15]. In our previous studies, some current peaks at short wavelengths (240 nm and 300 nm) cannot be interpretted by charge de-trapping mechanisms [28]. A strong photocurrent in the ultraviolet (UV) region has been observed and was suggested to originate from the photoelectric effect of the electrode, but no further evidences have been given [13,[START_REF] Kan | presented at ISEI[END_REF].

In this paper, photocurrent spectra for polypropylene films under different electric and geometrical conditions are studied. The samples with specifically deposited finger electrode are adopted to find out the origins of the photocurrent at different wavelengths and to confirm the electric conditions for generating these photocurrents. In order to identify the wavelength region of charge de-trapping current in the PSD spectra, a sample metalized with semi-transparent solid electrode is also tested during and after voltage applications and then their corresponding photocurrent spectra are compared. has been sputtered or evaporated in high vacuum on both sides of the sample as 30 mm diameter electrodes. Au electrode is semitransparent (20 nm thickness) and Al electrode is opaque (180 nm thickness). One kind of electrode has a circular shape (solid electrode) and the other kind has a finger-like shape (finger electrode). The samples with solid or finger electrodes are illustrated in figure 1. The samples are sandwiched by two plate copper electrodes of 30 mm diameter, one having a ring shape of 20 mm inner diameter in order to allow monochromatic light illuminating on the center of the sample. For short in the following, samples are named after their electrode combinations. For example, the sample with finger and solid Au respectively as front and back electrodes is called finger Au -solid Au sample.

Experiment

The samples are illuminated with monochromatic light generated by a 450 W Xe arc lamp (Osram XBO 450/4) resolved by a monochromator (Dongwoo Optron DM-201i). The photocurrent is measured with a Keithley 6514 electrometer while the illuminating light scans from longer to shorter wavelengths. The bandwidth of the monochromatic light is 20 nm and the scanning rate is 24 nm/min. The details of the measurement setup can be found in previous publications [19,[START_REF] Ma | [END_REF]. The light illumination side is the front electrode and the photocurrent from the front electrode to the back electrode is taken as the positive current. The charging duration before the scanning measurement is 1 hour. Measurements are carried out in air in a double screen dark room at room temperature and with a relative humidity of about 55%. Because the signal-to-noise ratio (SNR) of the output voltage from available positive voltage supply is not as good as that from negative voltage supply, we have used negative voltage supply as much as possible unless explicitly noted otherwise. The polarity reversal of the electric field in the sample is performed by connecting the negative voltage supply either to the front electrode or to the back electrode. The electrometer is then connected to the grounded electrode.

Results and discussion

The normalized current of a charged sample (photocurrent generated per incident photon on the sample surface) is used by many researchers to characterize the amount of charges de-trapped from their traps when illuminated by photons of a particular wavelength [START_REF] Mellinger | [END_REF][START_REF] Suzuoki | [END_REF]. However, our purpose is to find the possible processes involved in the experiment and to further distinguish the physical mechanisms from each other for different photocurrent peaks.

Besides, the spectral energy distribution of Xe lamp has the same feature among different scanning measurements.

Therefore, the original photocurrents under different electrical conditions can be directly compared to find the different causes of the peaks in the photocurrent spectra.

Photocurrent spectra for samples with different electrode combinations

Figure 2 shows the photocurrent spectra for samples with different electrode combinations, with the scanning wavelength in the UV region from 350 nm to 200 nm. It can be seen that there is no obvious photocurrent of solid Au -solid Au sample when -100 V is applied to either front or back electrodes. However, the photocurrent of finger Au -solid Au sample is different since almost no photocurrent is observed when -100 V is applied to the back electrode whereas a large photocurrent peak is detected when -100 V applied to the front electrode. The photocurrent peak features an extremum positioned at about 224 nm and a cutoff wavelength at about 260 nm.

Notice that considering the difference between finger and solid electrodes, the photocurrents must be generated from the edges of the finger electrode. Since accurate estimation of the peak cutoff is very difficult, the cutoff wavelength is here defined by the intersection of the linear interpolation of photocurrent falling edge with the baseline. In order to explain the large peak in the photocurrent spectra, the experiment is also carried out on finger Alsolid Au sample. The position of the peak is shifted to the longer wavelengths compared with that from finger Ausolid Au sample. The spectrum has its extremum and cutoff wavelength at about 250 nm and 290 nm respectively.

The observed shift seems to be related to the different work functions of the electrode. Because the peak extremum is associated with both light source spectral energy distribution and quantum yield by light wavelength, the cutoff wavelength of the current peak is used to represent the relationship between the photocurrent response and the electrode work function.

There is a positive correlation between the peak cutoff wavelength and the electrode work function as shown in table 1. We thus speculate that the peak results from electron emission via photoelectric effect from the front finger electrode. However, the differences between the work function and the peak cutoffs of electrode Au and Al are also obvious, both being 0.33 eV. This difference may be caused by the monochromatic light bandwidth and the local electric field effect near the electrode edge in this experiment. The local electric field effect is verified by the experiments under different voltages in section 3.3. There are two possible routes for the current flow from front electrode to the chamber wall or to the polymer material [33] as shown in the inset of figure 3. To check the route of the photocurrent, the measurement is carried out twice. For the first measurement, the electrometer is connected to the back electrode and -100 V voltage is applied to the front electrode. For the second measurement, the electrometer is connected to the front electrode and +100 V voltage is applied to the back electrode. As the electrode field in the sample is same in the two cases, if there is part of the current flowing to the chamber wall, the amplitude of the current in the first measurement should be smaller than that in the second measurement. The results are shown in figure 3. It can be seen that the shape of photocurrent spectra is almost the same and the amplitude of the current in the first measurement is even larger than that in the second measurement. This shows that, if the current from the electrode to the chamber wall exists, it is not the main contribution to the measured current under the UV light illumination. Additional experiments show that amplitude slightly reduces from one measurement to the following. This phenomenon can be explained by the electric field shielding effect as demonstrated in the next section.

Evolution of photocurrent spectra for long time UV illumination measurement

In order to find the reason for the amplitude decrease phenomenon in two successive measurements, a new finger Au -solid Au sample is used to carry out long time illumination experiments. The negative voltage -100 V is applied to the front electrode. Experimental procedures are as shown in figure 4. First, the sample is illuminated by monochromatic light at 224 nm for one hour without voltage application (between t0 and t1 as shown in figure 4). Second, the photocurrent spectrum is measured after one hour voltage application (denoted as "the 1st

scanning", between t2 and t3). Third, the sample is illuminated under 224 nm monochromatic light with the voltage on for about 50 minutes. Then the photocurrent spectrum is measured for the second time after the voltage being applied for another one hour (denoted as "the 2nd scanning", between t5 and t6). Last, the sample is thermally treated under short-circuit at 353 K for 4 hours and then the photocurrent spectra are measured for the third time after one hour voltage application (denoted as "the 3rd scanning", between t8 and t9). The experimental results are shown in figure 5. The photocurrent gradually decreased to 1 pA with the monochromatic illumination at 224 nm for 50 minutes, as shown in the inset of figure 5. A first assumption is that 50 minutes illumination at 224 nm may affect the film electrical property so as to induce the photocurrent decrease. Because the photocurrent spectrum for the 1st scanning is almost the same as the spectrum for the sample without enduring long time illumination (as shown in figure 3), this assumption is excluded. As mentioned above, the photoemission electrons from the electrode move toward the film. Then some charges can accumulate in polypropylene near the finger edge under the light illumination. Accordingly, the applied electric field near the edge is counter balanced by the local electric field from the accumulated charges. When the total electric field goes under a threshold value, the photoemission process can be tremendously weakened. Therefore, we assume that space charge shielding effect causes the decrease of the photocurrent under long time illumination. It can be seen that, even though there is one hour duration between the third step and the 2nd scanning, the current detected at 224 nm remains very low, about 1 pA. This indicates that the accumulated photoemission charges are very stable at room temperature. However, the photocurrent spectrum is partially recovered after the 4 hours heat treatment on the sample. It is reasonable to assume that accumulated charges move away from their positon, thus the amplitude of the shielding field weakens accordingly. Consequently, photoemissions from the electrode can be observed again in the 3rd scanning.

Figure 5. Evolution of photocurrent spectra due to long time monochromatic illumination at 224 nm when -100 V is applied to the front electrode. The inset shows the decrease of photocurrent with time under the illumination.

Comparison of photocurrent spectra under different voltages

Based on the above analysis, space charge accumulation occurs in polypropylene film under UV illumination but their positions are still to elucidate. Moreover, charge emission from electrode edge to polymer has also been reported without UV illumination by many researchers by using surface potential measurement [34,35]. Therefore, electric field emission from electrode can contribute along with photoemission to the measured results. In order to further understand the experimental results, experiments in a wide range of voltage levels (5 V, 15 V, 30 V, 100 V, 300 V, 500 V and 700 V) are carried out. Finger Au -solid Au samples are used in these experiments. First, a negative voltage is applied to the back electrode of the sample for one hour and then the photocurrent spectrum is measured. Second, the photocurrent spectrum is measured again after short-circuiting the sample 15 min in the dark. This experimental procedure is repeated when negative voltage is applied to the front electrode. For each measurement at different voltage levels, a new sample is used to avoid the effect of the remnant charges on the experimental results.

Figure 6(a) shows the photocurrent spectra when voltage from -5 V to -700 V is applied to the back electrode of the sample. Almost no photocurrent response is observed when the voltage amplitude is smaller than 100 V and the photocurrent gradually increases with the voltage amplitude above 300 V. Moreover, it can also be observed that the amplitude of each spectrum increases with the increase of photon energy up to about 210 nm. However, it seems that the photocurrent features a much steeper increase when the wavelength is shorter than about 250 nm.

Thus, the photocurrent below and above 250 nm must be involved in different charge migration mechanisms.

Considering the polarity of the measured current, two possibilities, i.e. photoemission from the back electrode and photoemission from polypropylene surface between the finger electrodes, can cause the photocurrent spectrum in the region below 250 nm. First, the current is generated by the photoemission from the back electrode. But the peak cutoff of the photocurrent is not at 260 nm as expected, which may be caused by the optical absorption spectrum of the polypropylene film. An inversed sample with finger electrode at the back side (solid Au -finger Au sample) is used to test that former possibility. For comparison, the back plate copper electrode used in this experiment also has a ring shape. The experimental results are shown in the inset of figure 6(a). It can be seen that the cutoff of the peak is the same as that measured when the front electrode is fingered. Therefore, it can be assumed that the photoemission from the surface states of polypropylene between the finger electrodes accounts for the photocurrent in the wavelength region below 250 nm. The surface states are formed due to the surface features of polymer film like the topological roughness, residues, and polymer oxidation etc [36]. Murata has also reported photoemission phenomena from surface states of polymers when he studied the effect of surface charging on surface states by using UV illumination measurement [25]. The associated mechanisms of photocurrent spectra above 250 nm will be discussed later. The short-circuit photocurrent spectra measured after the removal of the voltage are as shown in figure 6(b).

An obvious peak with a cutoff wavelength from about 259 nm to 266 nm is visible from the photocurrent spectra for samples that have been pre-charged under voltage above 300 V amplitude. There is no response in the photocurrent spectra for samples under voltage below 100 V amplitude. It can be noticed that the cutoff wavelength is almost the same than the one obtained under -100 V applied to the front electrode (as shown in figure 3). It can be speculated that positive charges accumulate in polypropylene near the finger electrode edge during the 1st step and thus a local positive electric field builds up. Therefore, photoelectrons are generated from the electrode edge if the local electric field is strong enough. However, the position of photocurrent extremum is about 240 nm which is 16 nm larger than that measured under -100 V applied to the front electrode. One possibility is that the peak with lower amplitude may be more easily influenced by the photoemission electrons from polymer which forms a reverse current, and the measured photocurrent is the photoemission current from electrode superimposed by the photoemission current from the polymer. But the spectrum measured under -30 V applied to the front electrode has almost the same amplitude, and has no wavelength shift as shown in figure 6(c).

Another possibility is that accumulated positive charges can be partially recombined by the electrons from the electrode during the scanning process. Accordingly, local electric field near the electrode edge weakens and so does the photoemission current amplitude at the shorter wavelengths. This results in a measured photocurrent peak with a shift toward longer wavelengths compared to photocurrents measured under -100 V. A question arises as to what physical mechanism accounts for the generation of the accumulated charges? They can come from field emission from electrode due to one hour high voltage application [34] and/or from the photoemission from the polymer during the UV illumination at the beginning of the procedure. We have carried out experiments on the sample without UV scanning measurement during the voltage application but a similar photocurrent peak is also observed during the short-circuit measurements. This indicates that field emission from finger electrode can also produce amounts of charges which finally generate a strong local electric field.

Figure 6(c) shows the photoemission current when a voltage from -5 V to -700 V is applied to the front electrode of the sample. The current peak slowly increases with the voltage amplitude, while the extremum position firstly slightly shifts to shorter wavelengths when the voltage goes from -30 V to -300 V and then shifts inversely to longer wavelengths when the voltage amplitude increases further. If electrons were generated just by the photoelectric effect, there should be a saturation of the photocurrent with the increase of applied voltage.

However, the inset of figure 6(c) does not show a saturation. Therefore, a field-assisted photoelectric emission mechanism can be expected, especially when higher voltage amplitude is applied [37]. In addition, the existence of field-assisted photoemission can illustrate the longer-wavelength-direction shift of the photocurrent peak when the voltage amplitude increases from 300 V to 700 V (cutoffs from 257 nm to 265 nm). Indeed, electrons can be emitted from the electrode by photons with less energy under higher electric field. It can be noticed that the cutoff wavelengths under different voltages are in the same range as that in figure6 (b). However, the shorter-wavelength-direction shift of the photocurrent peak at lower voltage from -30V to -300V is more difficult to understand. Here, a tentative explanation is given on the basis of charge shielding effect. Under very low voltages (-5 V to -30 V), the charge shielding effect is negligible because almost no accumulated charges exist near the electrode edge. Thus the peak position is not influenced. When the voltage is increased (-100 V to -300 V), charges accumulate in polypropylene film and produce a large shielding field compared with the applied electric field. Thus charges reduce the photoemission from the electrode. Therefore, the photocurrent peak between -100 V and -300 V is slightly shifted to shorter wavelengths. Therefore, the photocurrent at the wavelength shorter than 250 nm in these spectra results from the photoemission from polypropylene. A very broad peak can be observed above 250 nm, which may be associated with the photoconduction current. It has to be pointed out that the current above 250 nm is so low (about 0.4 pA) that it is difficult to notice it when measuring slightly large currents (e.g. figure 6(b) and (c)). For avoiding the disturbance from finger electrode edge, the current in the region above 250 nm is studied further by using solid Au -solid Au sample in section 3.4. It is quite obvious that the strong field emission and the photoemission from the electrode come from the electric field distortion at the finger electrode edge. Ansoft Maxwell, which is a powerful tool for electromagnetic field simulation, is used to model the electric field distribution. Space charge accumulation in the film or near the finger electrode edge is not taken into account in the modeling process to simplify the numerical simulations. Here, the results are focused in the region indicated in the top left of figure 1. Figure 7(a) and (b) show the simulation when 100 V is applied to the back electrode in the experiment. Two regions exhibit strong electric field distortion near point 1 and point 2 as shown in figure 7(a). Line 1 and line 2 are parallel to the polypropylene surface and electrode edges and pass to point 1 and point 2, respectively. To have a clear view of the field distribution, the field strength along line 1 and line 2 is plotted in figure 7(b). It can be seen that the field strength at point 1 and point 2 is much higher than the average field strength in the film bulk (110 7 V/m). At 24 nm distance from both point 1 and 2, the electric field strength is still about 10 times higher than the average field strength. The photoemission process from polypropylene probably happens near point 1 and the field emission and the photoemission from the electrode can happen both at point 1 and point 2.

Because the band gap of polymeric insulating material is usually above 8 eV [8], photoemission due to interband transition of polypropylene is excluded. Based on the above discussions, it can be assumed that the is due to either the photoemission from electrodes or the de-trapping process of trapped charges, which will be discussed below. As shown in figure 8, the obvious photoemission current peak is absent in all spectra when solid Ausolid Au samples are used. Almost no photocurrent peaks are observed in the measurements during or after the application of 100 V (curves 1 and 2). This indicates that no charges have accumulated in the film. However, two photoconduction current peaks at 273 nm and 320 nm are observed when -700 V is applied (curves 3 to 6).

Photocurrent spectra for solid Au -solid Au sample

Because surface states are tremendously diminished when solid electrodes are used [38], the two current peaks should be related to accumulated charges in the bulk of the film. It has been reported that homo space charges accumulate in polypropylene film measured by thermal pulse (TP) method [39]. Therefore, there also exists field shielding effect from the bulk charges when the sample is under high voltage, though there is no large electric field distortion without finger electrodes. Accordingly, the electric field strength can then be greatly depressed at the interface between the electrode and the polypropylene film.

Based on the above discussion, it is likely that photoemission from the electrode is impossible to occur. We then assume that the photocurrent peaks at 320 nm and 273 nm are both caused by the de-trapping process of the photo-stimulated charges. Because there is a discharge of the sample after the voltage removal, the amplitude of the photo-stimulated de-trapping current is smaller when the sample is under short-circuit than when it is under applied voltage as shown in figure 8. Because there is a positive correlation between the amplitude of the charge de-trapping current and the amount of space charges stored in the film, the ratio between the de-trapping current amplitudes at 273 nm and 320 nm in the spectra can be seen as the relative quantity index of charges in deeper and shallower traps in a sample. The ratio between the two peak amplitudes for one sample should be the same whenever the sample is either under applied voltage or in short-circuit. Even considering different lifetimes of charges for deeper and shallower traps, the ratio of the spectra in the short-circuit measurements (curves 4 and 6)

should be larger because of the relative faster decay rate of charges in 320 nm traps [START_REF] Ma | [END_REF]. But the experimental results in figure 8 do not confirm this view. The only explanation is that photoelectrons from electrodes partially contribute to the peak at 273 nm. Of course, the photoemission from polymer is also possible, but it is not considered here because the wavelength region below 250 nm has not obvious photocurrent peak. In conclusion, the current peak at 320 nm is certainly due to charge de-trapping mechanism and the current peak at 273 nm is certainly due to the combined contributions of charge de-trapping and photoemission from electrode.

Conclusion

The photocurrent spectra for samples under different electric and geometrical conditions have been compared in this paper. When a local strong positive electric field exists near the finger electrode, photoemission from finger electrode is observed. The generated electrons migrate to the polypropylene film and are partially captured by the surface states of polymer to form a stable charge accumulation near the finger electrode edge. The continued accumulation of charges can gradually increase the shielding field and finally inhibit the subsequent photoemission process. When a local large negative electric field is formed near the finger electrode, photoemission from surface states of polypropylene can be observed, with a cutoff wavelength at about 250 nm.

An obvious photoconduction of polypropylene film can be observed in the wavelength region above 250 nm under high electric field. The measured current spectrum includes photoemission current from the electrode and the charge de-trapping current in polypropylene film. The photoconduction current peaks at 273 nm and 320 nm are mainly induced by de-trapping charges and the peak at 273 nm is also partially contributed by the photoelectrons from the electrode.

Figure 1 .

 1 Figure 1. Sample geometry with either solid or finger electrodes. The electric field distribution in the indicated region is shown in figure 7.

Figure 2 .

 2 Figure 2. Comparison of photocurrent spectra for solid Au -solid Au, finger Au -solid Au and finger Al -solid Au samples when -100 V is applied to the front or the back electrode. The photocurrent spectra for solid Au -solid Au sample are vertically shifted to avoid overlapping.

Table 1 .Figure 3 .

 13 Figure 3. Comparison of photocurrent spectra when -100 V is applied to the front electrode and +100 V is applied to the back electrode. The inset shows the two possibilities of current routes in the experiment.

Figure 4 .

 4 Figure 4. Procedure for the measurement when -100 V is applied to the front electrode of the sample. Time t0 and t3: monochromatic light on; Time t1 and t7: voltage on; Time t2, t5 and t8: light scanning started; Time t4, t6: light off; Time t9: voltage off and light off.

Figure 6 .

 6 Figure 6. Comparison of photocurrent spectra for samples (a) when a negative voltage is applied to the back electrode, (b) after the removal of a negative voltage from the back electrode, (c) when a negative voltage is applied to the front electrode, and (d) after the removal of the negative voltage from the front electrode. The inset of figure 6 (a) shows the photocurrent spectra for solid Au -finger Au sample, and the inset of figure 6 (c) shows the dependence of peak amplitude on the applied voltage.

Figure 6 (

 6 Figure 6(d) shows similar photocurrent spectra as those shown in figure6(a). A steep increase can also be observed when the wavelength is below about 250 nm. As in figure6(b), negative charges can also have accumulated in polypropylene film. This forms similar local negative electric field near the finger electrode edge.

Figure 7 .

 7 Figure 7. Electric field distribution near the finger electrode edge (a) in two dimensions and (b) along line 1 and line 2 when the finger Au -solid Au sample is under voltage 100 V. Point 1 and point 2 are the electrode tips in the paper plane and line 1 and line 2 are two lines passing through the two points.

  electrons trapped in the surface states of polypropylene can contribute to the photoemission current in the wavelength region below 250 nm only if a local strong negative electric field is formed near the edge of finger electrode. Photoemission from the electrode happens under a local strong positive electric field near the edge of finger electrode and is affected by both the electrode work function and the local electric field. The charged surface states region can extend to a certain distance along the film surface from the finger electrode edge because of the migration of charges under their own electric field. If the photoelectron accumulation locates in the bulk of the film, the current above the wavelength 250 nm in figure 6(d) should be much larger than that measured in the experiment. Hence, most accumulated charges should just be trapped in the surface states of polypropylene near finger electrode edges, which accounts for the formation of large shielding field. The photocurrent above 250 nm

Figure 8 .

 8 Figure 8. Photo-stimulated discharge spectra for the solid Au -solid Au sample.

Figure 1 .

 1 Figure 1. Sample geometry with either solid or finger electrodes. The electric field distribution in the indicated region is shown in figure 7.

Figure 2 .

 2 Figure 2. Comparison of photocurrent spectra for solid Au -solid Au, finger Au -solid Au and finger Al -solid Au samples when -100 V is applied to the front or the back electrode. The photocurrent spectra for solid Au -solid Au sample are vertically shifted to avoid overlapping.

Figure 3 .

 3 Figure 3. Comparison of photocurrent spectra when -100 V is applied to the front electrode and +100 V is applied to the back electrode. The inset shows the two possibilities of current routes in the experiment.

Figure 4 .

 4 Figure 4. Procedure for the measurement when -100 V is applied to the front electrode of the sample. Time t0 and t3: monochromatic light on; Time t1 and t7: voltage on; Time t2, t5 and t8: light scanning started; Time t4, t6: light off; Time t9: voltage off and light off.

Figure 5 .

 5 Figure 5. Evolution of photocurrent spectra due to long time monochromatic illumination at 224 nm when -100 V is applied to the front electrode. The inset shows the decrease of photocurrent with time under the illumination.

Figure 6 .

 6 Figure 6. Comparison of photocurrent spectra for samples (a) when a negative voltage is applied to the back electrode, (b) after the removal of a negative voltage from the back electrode, (c) when a negative voltage is applied to the front electrode, and (d) after the removal of the negative voltage from the front electrode. The inset of figure6 (a) shows the photocurrent spectra for solid Au -finger Au sample, and the inset of figure6 (c) shows the dependence of peak amplitude on the applied voltage.

Figure 7 .

 7 Figure 7. Electric field distribution near the finger electrode edge (a) in two dimensions and (b) along line 1 and line 2 when the finger Au -solid Au sample is under voltage 100 V. Point 1 and point 2 are the electrode tips in the paper plane and line 1 and line 2 are two lines passing through the two points.

Figure 8 .

 8 Figure 8. Photo-stimulated discharge spectra for the solid Au -solid Au sample.

Figure 2 .

 2 Figure 2.

Figure 6 .

 6 Figure 6.

Figure 8 .

 8 Figure 8.

Table 1 .

 1 Relationship between work function of the electrode material and photocurrent peak cutoff.

	electrode material work function (eV) [32] cutoff (nm)/energy (eV)
	Au	5.1	260 / 4.77
	Al	4.28	314 / 3.95

Acknowledgments

We gratefully acknowledge the financial support from the National Natural Science Foundation of China (NSFC Nos. 51477118 and 51277133) and PhD Programs Foundation of Ministry of Education of China (No. 20130072110046).