M Peigney

J P Seguin

An incremental variational approach to coupled thermo-mechanical problems in anelastic solids. Application to shape-memory alloys

Keywords: incremental variational principles, thermo-mechanical coupling, non-smooth mechanics, shape memory alloys

We study the coupled thermo-mechanical problem that is obtained by combining generalized standard materials with Fourier's law for heat conduction. The analysis is conducted in the framework of non-smooth mechanics in order to account for possible constraints on the state variables. This allows models of damage and phase-transformation to be included in the analysis. In view of performing numerical simulations, an incremental thermo-mechanical problem and corresponding variational principles are introduced. Conditions for existence of solutions to the incremental problem are discussed and compared with the isothermal case. The numerical implementation of the proposed approach is studied in detail. In particular, it is shown that the incremental thermo-mechanical problem can be recast as a concave maximization problem and ultimately amounts to solve a sequence of linear thermal problems and purely mechanical (i.e. at a prescribed temperature field) problems. Therefore, using the proposed approach, thermo-mechanical coupling can be implemented with low additional complexity compared to the isothermal case, while still relying on a sound mathematical framework. As an application, thermo-mechanical coupling in shape memory alloys is studied. The influence of the loading strain-rate on the phase transformation and on the overall stress-strain response is investigated, as well as the influence of the thermal boundary conditions. The numerical results obtained by the pro-

Introduction

This paper focuses on coupled thermo-mechanical evolutions of dissipative solids, in the geometrically linear (small strains) setting. The framework of generalized standard materials in non-smooth mechanics is considered [START_REF] Halphen | Sur les matériaux standards généralisés[END_REF][START_REF] Moreau | Nonsmooth mechanics and applications[END_REF][START_REF] Frémond | Non Smooth Thermo-mechanics[END_REF]. In that framework, the local state of the material is described by the strain ε, the temperature θ, and an internal variable α. The constitutive laws are determined from the Helmholtz free energy w and a convex dissipation potential Φ. In its original form [START_REF] Halphen | Sur les matériaux standards généralisés[END_REF], that framework covers a wide range of elasto-plastic models, including limited and nonlinear hardening. Its extension to non-smooth mechanics has been extensively studied by [START_REF] Frémond | Non Smooth Thermo-mechanics[END_REF] and allows constraints on the internal variable α to be taken into account in a rigorous fashion. That feature is crucial for the modelling of such phenomena as damage or phase-transformation, as the internal variable in such cases is typically bounded. The thermodynamic analysis of the media considered is presented in Section 2, leading to a boundary value problem for the mechanical and thermal fields. As pointed out by [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF], the time-discretization of the thermo-mechanical evolution problem is a sensitive issue because of the coupling between mechanical and thermal equations. For instance, the Euler implicit scheme leads to an incremental thermo-mechanical problem for which existence of solutions cannot generally be ensured. This is in contrast with the isothermal case, for which the Euler implicit scheme provides a well-posed incremental problem under standard assumptions of convexity on the functions w and Φ.

One objective of this paper is to propose a sound time-discretization scheme for coupled thermo-mechanical problems, retaining some essential features displayed by the Euler scheme in the isothermal case (most notably the consistency with the rate problem and the existence of solutions). A central idea is the use of a variational formulation for the incremental problem. Incremental variational principles for dissipative solids have been the focus of a lot of attention in recent years, offering new perspectives in various topics such as finite-strains elasto-viscoplasticity [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF], homogenization [START_REF] Miehe | Strain-driven homogenization of inelastic microstructures and composites based on an incremental variational formulation[END_REF][START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I. incremental variational principles[END_REF], formation and stability of microstructures [START_REF] Ortiz | Nonconvex energy minimization and dislocation structures in ductile single crystals[END_REF][START_REF] Miehe | Analysis of material instabilities in inelastic solids by incremental energy minimization and relaxation methods: evolving deformation microstructures in finite plasticity[END_REF]. Incremental variational principles for coupled thermo-mechanical problems have been proposed by [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF] in the case where the heat flux q derives from a potential χ in (∇θ)/θ, i.e. when the heat conduction law takes the form q = -χ ((∇θ)/θ). In this article, we stick with the standard Fourier's law of heat conduction q = -K∇θ, which does not fall in the format considered by [START_REF] Yang | A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids[END_REF].

In Section 3 we introduce an incremental problem for the class of coupled thermo-mechanical problems considered, along with a corresponding variational formulation. The variational formulation of the incremental thermomechanical problem serves two purposes. First, it allows the existence of solutions to be studied, as discussed in Section 3. Second, the variational formulation leads to a convenient and efficient way of solving the incremental thermo-mechanical problem. The latter can be indeed be recast as a concave maximization problem, for which well-known algorithms are available. As detailed in Section 4, an advantage of that approach is that the solution of the thermo-mechanical problem can be obtained by solving a sequence of linear thermal problems and purely mechanical (i.e. at prescribed temperature) problems. This calls for an easy implementation in an existing finite-element code. A crucial point in the analysis lies in the introduction of an auxiliary linear problem, akin to the adjoint state used in optimal control problems [START_REF] Lions | Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles[END_REF].

As an application, the proposed method is used in Section 5 to study thermo-mechanical coupling in shape-memory alloys. The significant role of thermal effects in shape-memory alloys has notably been put forward by [START_REF] Peyroux | Thermomechanical couplings and pseudoelasticity of shape memory alloys[END_REF]; [START_REF] Chrysochoos | A one-dimensional thermomechanical modeling of phase change front propagation in a sma monocrystal[END_REF]. The solid/solid phase transformation that occurs in those materials is known to produce significant amounts of heat, associated both with recoverable latent heat effects and irreversible frictional contributions. Depending on the rate of loading and on the thermal exchange conditions, the heat produced by the phase transformation may not have time to diffuse in the body and the temperature field may become inhomogeneous. In such conditions, the overall stress-strain response becomes significantly different from its isothermal counterpart, and it is mandatory to take the thermo-mechanical coupling into account. Therefore, shape-memory alloys offer a particularly relevant application of the general methods presented in this paper. In Section 5, the influence of thermal effects on the phase-transformation and on the overall stress-strain curve is investigated in detail.

Thermo-mechanical evolutions of continuous media

Thermodynamic principles

Consider the evolution (on a time interval [0, T ]) of a continuous medium occupying a domain Ω in the reference configuration. We restrict our attention to the geometrically linear setting, defining the strain ε as ε = 1/2(∇u + ∇ T u) where u is the displacement. The first principle of thermodynamics gives

t t Ėdτ + t t Kdτ = t t P dτ + t t Qdt for all 0 ≤ t ≤ t ≤ T. (1)
In equation ( 1), K and E are respectively the kinetic and the internal energy of the system. The internal energy E can be written in the form E = Ω edx where e is the internal energy density. In the right-hand side of (1), P denotes the power of external loads, and Q is the rate of heat received by the system. The upper dot in (1) denotes left-time derivative1 . The principle of virtual power gives the relation

K = P - Ω σ : εdω (2)
where σ is the stress. Expressing Q as

Q = - ∂Ω q.ndω + Ω rdω (3)
where q is the heat flux and r a heat source, the relation (1) can be rewritten as

t t Ω ( ė -σ : ε + div q -r) dωdτ = 0 for all 0 ≤ t ≤ t ≤ T. (4)
The relation (4) also holds when replacing Ω with an arbitrary subdomain Ω ⊂ Ω. Therefore, we obtain the local equation

ė -σ : ε + div q -r = 0 a.e. in Ω × [0, T ] (5)
where the abbreviation 'a.e' stands for 'almost everywhere'. The second principle of thermodynamics gives

t t Ω
Equations (5-6) can be rewritten in terms of the Helmholtz free energy density w = e -θs as

ẇ + θ ṡ + s θ = σ : ε + r -div q a.e. in Ω × [0, T ], (7) σ : ε -s θ -ẇ ≥ 0 a.e. in Ω × [0, T ].
(8)

Mechanical constitutive laws

In the framework of standard generalized materials [START_REF] Halphen | Sur les matériaux standards généralisés[END_REF], the local state of the material is described by the strain ε, the temperature θ, and an internal variable α living in a vectorial space denoted by A. The constitutive laws are determined by the Helmholtz free energy w(ε, α, θ) and a convex dissipation potential Φ( α) according to the following relations: σ = ∂w ∂ε , (9.1) (9.4) where ∂ denotes the subdifferential operator. Recall [START_REF] Brézis | Opérateurs maximum monotones et semigroupes de contractions dans les espaces de hilbert[END_REF]) that the subdifferential ∂f of a function f : A → R is the multi-valued mapping defined by

A = - ∂w ∂α , (9.2) s = - ∂w ∂θ , (9.3) A ∈ ∂Φ( α),
∂f (x) = {τ ∈ A|f (y) -f (x) ≥ τ .(y -x) ∀y ∈ A}. ( 10
)
In the following, the dissipative behaviour is assumed to be rate-independent. In such case, the dissipation potential Φ is positively homogeneous of degree 1, i.e. satisfies Φ(λ α) = λΦ( α) for any λ ∈ R + and α ∈ A.

The property (11) obviously implies that Φ(0) = 0. Moreover, combining (11) with the convexity of Φ, it can easily be shown that Φ( α) is positive for all α. For latter reference, we note the following property that is a direct consequence of ( 10) and ( 11):

∂Φ(λ α) = ∂Φ( α) for any λ ∈ R + and α ∈ A. ( 12
)
Note that the constitutive laws (9) satisfy the inequality (8). From (9) we have indeed σ : ε -s θ -ẇ = A. α.

Since A ∈ ∂Φ( α), the definition (10) gives

Φ(0) -Φ( α) ≥ -A. α. ( 13
)
As noted above, Φ is positive and vanishes at 0. The left-hand side of ( 13) is thus non-positive, in compliance with the thermodynamical principle (8).

So far we have assumed that α is unconstrained, in the sense that α is allowed to take any value in A. In order to model such phenomena as damage or phase-transformation, it is essential to consider an extension of (9) to situations where the internal variable α is constrained, in the sense that α is required to satisfy a condition of the form α ∈ T where T is a given (usually bounded) convex subset of A. In that case, as notably detailed by [START_REF] Frémond | Non Smooth Thermo-mechanics[END_REF], equations ( 9) are modified as

σ = ∂w ∂ε , A = - ∂w ∂α , s = - ∂w ∂θ , A = A d + A r , A d ∈ ∂Φ( α), A r ∈ ∂I T (α) (14)
where I T (α) is the indicator function of T (equal to 0 if α ∈ T , and infinite otherwise). The superscripts 'd' and 'r' in (14) stand for 'dissipative' and 'reversible', respectively. Let us verify that the constitutive equations ( 14) are consistent with the thermodynamic principles. From ( 14) we find

σ : ε -s θ -ẇ = A d . α + A r . α. ( 15
)
The term A d . α can be proved to be positive using the same reasoning as in the unconstrained case. Since α(t) ∈ T for all t and A r ∈ ∂I T (α), the definition (10) implies that

I T (α(t -dt)) -I T (α(t)) ≥ A r .(α(t -dt) -α(t)) for all dt. (16)
By letting dt tend towards 0 from above, we obtain that A r . α ≥ 0 where α is the left-time derivative. The left-hand side of ( 15) is thus positive, in accordance with (8). Note from (16) that if α is time-differentiable (i.e. both left-and right-time derivatives exist and are equal), then A r . α = 0. This remark will be useful in the next section. The property A r . α = 0 also shows that A r does not contribute to the energy dissipation, which explains why A r is referred to as a 'reversible' term.

Boundary value problem

We now formulate the boundary value problem that governs quasi-static evolutions of a continuous medium submitted to a prescribed loading history. Body forces f d are applied in the domain Ω. Displacements u d are imposed on a part Γ u of the boundary Γ, and tractions T d are prescribed on Γ T = Γ-Γ u . The given functions f d , u d , T d as well as the stress and state variables (σ, ε, α) in Ω, depend on the location x and the time t. However, in order to alleviate the expressions, this dependence will be omitted in the notations.

Combining the principle of virtual power (2) with the constitutive laws (14) leads to the following set of relations

σ ∈ K σ , u ∈ K u , α ∈ K α ; σ = ∂w ∂ε , A = - ∂w ∂α , s = - ∂w ∂θ ; A = A d + A r , A d ∈ ∂Φ( α) , A r ∈ ∂I T (α); (17)
where the sets K σ , K , K α are defined by

K σ = {σ| div σ + f d = 0 in Ω; σ.n = T d on Γ T }, K u = {u|u = u d on Γ u }, K α = {α|α(x) ∈ T ∀x in Ω}.
In the isothermal case, the system (17) completely determines the evolution of the structure from a given initial state. It is important to note that the term A r in (17) has a profound impact on the behaviour of the system compared to the classical plastic case. A more detailed discussion along those lines can be found in [START_REF] Peigney | Shakedown theorems and asymptotic behaviour of solids in non-smooth mechanics[END_REF].

Note that the free energy w depends on the temperature θ, so that the temperature appears implicitly in the last two equations of (17). In the coupled thermo-mechanical case, the temperature field θ(x, t) is unknown and needs to be solved for. This is accomplished by using the energy balance equation ( 7) together with a constitutive law for heat conduction. In that regard, the most classical model of heat conduction is the Fourier's law q = -K∇θ, which we adopt in the following. To simplify the presentation, the (positive) thermal conductivity K is assumed to be independent on x. Combining the energy balance (7) with Fourier's law yields the heat equation

K∆θ -θ ṡ + A d . α + A r . α + r = 0 a.e. in Ω × [0, T ]. ( 18
)
We assume in the following that α(x, t) is time-differentiable almost everywhere in Ω × [0, T ]. As noted earlier, equation ( 16) then implies that A r . α = 0 a.e. in Ω × [0, T ]. The equation ( 18) is complemented by thermal boundary conditions defined as follows: the temperature is prescribed to take a given value θ d (x, t) on a portion Γ θ of the boundary ∂Ω. The heat flux is prescribed to take a given value q d (x, t) on a portion Γ q such that Γ q ∩ Γ θ = ∅. On Γ h = ∂Ω -Γ θ -Γ q , we consider a convection condition of the form

q.n = h(θ -θ R ) (19)
where θ R is the ambient temperature of the surrounding medium and h is a (positive) heat transfer coefficient. The temperature field θ is thus a solution of the boundary value problem

θ ∈ K θ ; -K∇θ.n = q d on Γ q , -K∇θ.n = h(θ -θ R ) on Γ h ; K∆θ -θ ṡ + A d . α + r = 0 a.e. in Ω × [0, T ]; (20)
where

K θ = {θ|θ = θ d on Γ θ }.
The mechanical variables (α, ε) have an influence on the solution θ of (20), both through the dependence of the entropy s on (α, ε) and through the dissipative term A d . α. In a similar fashion, as mentioned earlier, the solution of the mechanical problem ( 17) depends on the temperature field. The systems ( 17) and ( 20) are thus to be considered as coupled.

Variational formulation of an incremental problem

To solve a system such as ( 17)-( 20), one generally resorts to a timediscretization strategy: the time history is discretized as a sequence t 1 < t 2 < • • • < t n and one estimates successively the fields at time t i using a finite time-step problem. That problem is an approximation of the timecontinuous problem, allowing the fields (σ,u,α,θ) at t 0 + δt (with δt > 0) to be estimated from their values (σ 0 ,u 0 ,α 0 ,θ 0 ) at t 0 . Such a time-discretization scheme is required to be consistent with the rate problem, in the sense that the finite-step problem coincides with the rate problem (at least formally) as the time increment δt tends towards 0. A second natural requirement is that the finite-step problem admits some solutions. In the isothermal case, such requirements are classically satisfied by the Euler implicit scheme. The corresponding finite-time step problem is

σ ∈ K σ , u ∈ K u , α ∈ K α ; σ = ∂w ∂ε , A = - ∂w ∂α ; A = A d + A r , A d ∈ ∂Φ( α -α 0 δt ) , A r ∈ ∂I T (α). ( 21
)
The incremental problem (21) can easily be verified to be consistent with the rate problem (17). Existence of solutions can be studied by using a variational formulation attached to (21). To that purpose, introduce the functionals F e and F d 0 defined as

F e (u, α, θ) = Ω w(ε, α, θ)dω - Ω f d .u dω - Γ T T d .u da; F d 0 (α) = Ω Φ α -α 0 dω. (22)
The functional F e is the potential elastic energy of the system, while F d 0 is related to the energy dissipated on the time interval [t 0 , t 0 + δt].

The functional F 0 = F e + F d 0 is differentiable with respect to u but only subdifferentiable with respect to α. The derivative of F 0 with respect to u is denoted by ∂F 0 /∂u. The subdifferential of F 0 is denoted by ∂ α F 0 and is characterized as follows: for any

A ∈ ∂ α F 0 (u, α, θ), there exists a field A d such that A d ∈ ∂Φ(α -α 0 ) and A.(α * -α) = Ω ( ∂w ∂α + A d ).(α * -α)dω ∀α * ∈ K α .
We note that A.(α * -α) can be interpreted as the directional derivative of F 0 in the direction α * -α [START_REF] Rockafellar | Convex Analysis[END_REF].

The introduction of the functional F 0 is motivated by the following property: any solution of ( 21) is a solution of the variational problem

F ind (u, α) ∈ K u × K α and A ∈ ∂ α F 0 (u, α, θ) such that : 0 ≤ ∂F 0 ∂u .(u * -u) + A.(α * -α) ∀(u * , α * ) ∈ K u × K α . (23)
Let us justify that statement. The stationarity conditions with respect to u in (23) give the equations σ = ∂w/∂ε and σ ∈ K σ . The stationarity conditions with respect to α give

-A.(α * -α) + A d .(α * -α) ≥ 0 ∀α * ∈ T ; ( 24
)
where A = -∂w/∂α and A d ∈ ∂Φ(α -α 0 ). The equation ( 24) implies that

A -A d ∈ ∂I T (α), i.e. there exists A r ∈ ∂I T (α) such that A = A d + A r .
The variational formulation (23) means that the directional derivative of F 0 in every direction is positive. Such a condition is notably satisfied if F 0 reaches a local minimum at (u, α). In the particular case where F 0 is convex, it is known that local minima coincide with global minima, and are the only solutions of the problem ( 23). In such condition, the problem ( 23) is equivalent to min

(u,α)∈Ku×Kα F e (u, α, θ) + F d 0 (α) (25)
and admits solutions in adequate functional spaces for K u and K α . This in turn ensures the existence of solutions to the isothermal incremental problem (21).

Note that the functional F d 0 is convex. A sufficient condition for F 0 to be convex is thus that the free energy w is convex in (ε, α), which is notably satisfied by a wide range of elastoplasticity models. We note, however, that the requirement of convexity on w is not mandatory for (25) to have a solution. For instance, arguments related to quasiconvexity can be used to study the existence of solution to (25) in a more general setting (see e.g. [START_REF] Dacorogna | Direct methods in the calculus of variations[END_REF]).

Let us give some interpretation of (25). When there is no dissipation, it is well known that the solutions of the equilibrium problem minimize the potential elastic energy F e . The relation ( 25) can be interpreted as an extension of the principle of energy minimization to dissipative evolutions: the state of the system at time t 0 + δt is obtained by minimizing the 'incremental energy' F 0 . In the expression of F 0 , the dissipative contribution F d 0 depends on the state α 0 of the system at time t 0 . Therefore, the state of the system at time t 0 + δt also depends on α 0 , which reflects the path-dependence of dissipative evolutions.

We now move to the coupled thermo-mechanical problem (17)(18)(19)(20). A natural choice is to use again the Euler implicit scheme, which now reads:

σ ∈ K σ , u ∈ K u , α ∈ K α , θ ∈ K θ ; (26.1) σ = ∂w ∂ε , A = - ∂w ∂α , s = - ∂w ∂θ ; (26.2) A = A r + A d , A r ∈ ∂I T (α) , A d ∈ ∂Φ (α -α 0 )/δt ; (26.3) -K∇θ.n = q d on Γ q , -K∇θ.n = h(θ -θ R ) on Γ h ; (26.4) Kδt∆θ -θ 0 (s -s 0 ) + A d .(α -α 0 ) + rδt = 0; (26.5)
where s 0 = s(ε 0 , α 0 , θ 0 ). That incremental problem is obviously consistent with the rate problem. In particular, dividing (26.5) by δt and taking the limit δt → 0 yields the rate form of the heat equation in (20). However, existence of solutions to (26) cannot be proved in general. This is essentially due to the lack of a variational formulation for (26): in contrast with the isothermal problem ( 21), the incremental problem (26) does not correspond to the stationarity conditions of a certain functional.

In order to avoid such difficulties, consider instead the following finite time-step problem: 26) in two points: a factor θ/θ 0 is added in (27.3), and a term Kδt∇θ 0 .∇(θ 0 -θ)/θ 0 is added in the heat equation (27.5). It can be verified that -just as the more intuitive scheme (26) -the incremental problem ( 27) is a consistent time-discretization of (17-20). In particular, the heat equation in ( 20) is again recovered by dividing (27.5) by δt and taking the limit δt → 0. The crucial point here is that the extra term Kδt∇θ 0 .∇(θ 0 -θ)/θ 0 is of the second order in δt.

σ ∈ K σ , u ∈ K u , α ∈ K α , θ ∈ K θ ; (27.1) σ = ∂w ∂ε , A = - ∂w ∂α , s = - ∂w ∂θ ; (27.2) A = A r + θ θ 0 A d , A r ∈ ∂I T (α) , A d ∈ ∂Φ (α -α 0 )/δt ; (27.3) -K∇θ.n = q d on Γ q , -K∇θ.n = h(θ -θ R ) on Γ h ; (27.4) Kδt[∆θ + ∇θ 0 θ 0 .∇(θ 0 -θ)] -θ 0 (s -s 0 ) + A d .(α -α 0 ) + rδt = 0. (27.5) That problem differs from (
Motivation of the scheme ( 27) is that a variational formulation can be given. The corresponding functional is an extension of ( 22) to thermomechanics, and is given by

F(u, α, θ) = F e (u, α, θ) + F d (α, θ) + F θ (θ)
where

F d (α, θ) = Ω θ θ 0 Φ α -α 0 dω; F θ (θ) = Ω θ(s 0 + δt r θ 0 ) dω + δt Ω K(- 1 2 1 θ 0 ∇θ 2 + ( ∇θ 0 θ 0 ) 2 θ) dω -δt Γq q d θ 0 θ da -h δt 2 Γ h (θ -θ R ) 2 θ 0 da. (28)
Note that F θ is a quadratic and concave function of θ. In a way similar to the isothermal case, the functional F is differentiable with respect to u and θ (with partial derivatives denoted by ∂F/∂u and ∂F/∂θ) and only subdifferentiable with respect to α. The subdifferential of F with respect to α is the multivalued mapping characterized as follows: for any A ∈ ∂ α F(u, α, θ), there exists a field A d such that A d ∈ ∂Φ(α -α 0 ) and

A.(α * -α) = Ω ( ∂w ∂α + θ θ 0 A d ).(α * -α)dω ∀α * ∈ K α .
As detailed in Appendix A, any solution of ( 27) is a solution of the following variational problem

F ind (u, α, θ) ∈ K u × K α × K θ and A ∈ ∂ α F(u, α, θ) such that 0 ≤ ∂F ∂u .(u * -u) + ∂F ∂θ .(θ * -θ) + A.(α * -α) ∀(u * , α * , θ * ) ∈ K u × K α × K θ . (29)
As in the isothermal case, the variational formulation (29) allows the existence of solutions to (27) to be studied. Assume in particular that the free energy w is convex in (u, α) and concave in θ. In such a situation, the functional F is convex with respect to the fields (u, α) and concave with respect to θ, so that a saddle point exists (in adequate functional spaces for K u , K α and K θ , depending on the growth behaviour of w and Φ at infinity). Such a saddle point (u, α, θ) verifies (29) and therefore is a solution of the incremental problem (27).

The incremental thermo-mechanical problem ( 27) is an extension of an incremental problem introduced in Peigney ( 2006) for a simplified thermomechanical setting (in which the dissipative contribution A d . α is neglected in the heat equation ( 20)).

A maximization approach for solving the incremental problem

To solve a problem such as ( 26) or ( 27), a general strategy is to directly solve the local equations using for instance a Newton-Raphson algorithm. In such a framework, a partitioning approach is often used: the mechanical and the thermal subproblems are decoupled and solved successively until convergence, as described in Algorithm 1 below (see e.g. Auricchio and Petrini (2004) for an example of that approach). The global convergence Algorithm 1 Partitioning method for the thermo-mechanical problem (26) k ← 0 while residual > tolerance do Compute (u k+1 , α k+1 ) as the solution of (26.1-3) at θ = θ k Compute θ k+1 as the solution of (26.4-5) at (u, α) = (u k+1 , α k+1 ) k ← k + 1 end while of such methods is not ensured, and in practice one can face difficulties of convergence notably when the initial guess is not close enough to the solution.

Observe that, in the case of ( 27), such strategies ignore the variational nature of the problem at hand. As an alternative, using the variational formulation of the problem, the solution of ( 27) can be found by solving a concave maximization problem, as detailed in the following. We assume that w is convex in (u, α) and concave in θ, which ensures that (27) has a solution. More precisely, as explained at the end of section 3, solutions of ( 27) are saddle points of F, i.e solutions of the max-min problem max

θ∈K θ min (u,α)∈Ku×Kα F e (u, α, θ) + F d (α, θ) + F θ (θ). ( 30
)
The problem (30) can be rewritten as max

θ∈K θ J(θ) (31)
where the functional J is defined as

J(θ) = F θ (θ) + min (u,α)∈Ku×Kα {F e (u, α, θ) + F d (α, θ)} (32)
and can be proved to be concave (see Appendix B).

Solving the set of partial differential equations ( 27) amounts to solve the concave maximization problem (31) with respect to the temperature field. A lot of well-known algorithms can be used to solve such a maximization problem, some of them being built-in functions of scientific calculation softwares. Such algorithms (like the Broyden-Fletcher-Goldfarb-Shanno algorithm for instance) are iterative and typically require the computation of J and its gradient J (or at least of an ascent direction) at each iteration. In this regard, note from (32) that the calculation of J(θ) amounts to solve the minimization problem min

(u,α)∈Ku×Kα F e (u, α, θ) + F d (α, θ) (33)
for which the local equations (expressing the stationarity of the functional) read as

u ∈ K u , σ ∈ K σ , α ∈ K α ; σ = ∂w ∂ε , A = - ∂w ∂α ; A = A d + A r , A d ∈ θ θ 0 ∂Φ( α -α 0 δt ) , A r ∈ ∂I T (α). ( 34
)
That problem is formally identical to the isothermal problem ( 21), with a dissipation potential set equal to (θ/θ 0 )Φ. The calculation of J(θ) thus amounts to solve a incremental problem at a fixed temperature field. Special care must be taken in the calculation of J (θ) because, due to the non-differentiable nature of Φ, it is not ensured that J is differentiable everywhere. However, since J is concave, there exists a directional derivative in every direction [START_REF] Rockafellar | Convex Analysis[END_REF]. For a given θ ∈ K θ and θ such that θ = 0 on Γ θ , the directional derivative DJ(θ; θ) at θ in direction θ is defined by

DJ(θ; θ) = lim t-→0 + J(θ(t)) -J(θ) t
where θ(t) = θ + t θ ∈ K θ . As detailed in Appendix B, the directional derivative of J satisfies the property

DJ(θ, θ) ≥ ∂F ∂θ . θ (35) with ∂F ∂θ . θ = Ω [θ 0 (s 0 -s) + rδt + A d .(α -α 0 ) +Kδt( ∇θ 0 2 θ 0 - ∇θ.∇θ 0 θ 0 + ∆θ)] θ θ 0 dω -δt Γq q d + K∇θ.n θ 0 θ da -δt Γ h h(θ -θ R ) + K∇θ.n θ 0 θ da. (36)
In order to determine an ascent direction for J, consider the solution W of the following linear problem:

∆W -W = - 1 θ 0 [Kδt(∆θ + ∇θ 0 θ 0 .∇(θ 0 -θ) -θ 0 (s -s 0 ) + A d .(α -α 0 ) + rδt]; W = 0 on Γ θ ; ∇W.n = δt θ 0 (K∇θ.n + q d ) on Γ q ; ∇W.n = δt θ 0 (K∇θ.n + h(θ -θ R )) on Γ h .
(37) Substituting in the expression (36) gives

∂F ∂θ . θ = Ω (-∆W + W ) θdω - ∂Ω θ∇W.nda = Ω (∇W.∇ θ + W θ)dω.
The relation ( 35) then implies that

DJ(θ; θ) ≥ Ω (∇W.∇ θ + W θ)dω.
In particular, we have DJ(θ; W ) = Ω (W 2 + ∇W 2 ) ≥ 0, i.e. W is an ascent direction for J. Note that DJ(θ; W ) is null if and only if θ is solution of the maximization problem (31). The introduction of an auxiliary problem for determining an ascent direction is reminiscent of techniques used in optimal control theory [START_REF] Lions | Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles[END_REF]. That theory addresses minimization problems of the form min a J(s(a)) in which s(a) is solution of a set of partial differential equations parametrized by a. In the optimal control terminology, a is the control variable and s(a) is the state variable. In general, the set of partial differential equations that defines s(a) cannot be solved in closed-form, so that the dependence of J with respect to a remains implicit. The explicit determination of the gradient J is thus not straightforward. An effective method is to introduce a so-called adjoint state for expressing J . That adjoint state is usually defined as the solution of an ad hoc set of linear partial differential equations. In solid mechanics, the optimal control theory has been extensively used for inverse problems [START_REF] Bui | Fracture Mechanics, Inverse Problems and Solution[END_REF] and also proved to be useful for other classes of nonlinear problems (Peigney andStolz, 2001, 2003;[START_REF] Stolz | Optimal control approach in nonlinear mechanics[END_REF]. Note that the maximization problem (31) can be interpreted as an optimal control problem, the temperature field θ being the control variable and the mechanical fields (u, α) being the state variable. The fields (u, α) are indeed solution of the set of partial differential equations (34) in which θ

Algorithm 2 Maximization method for the thermo-mechanical problem ( 27) k ← 0 while residual > tolerance do Compute J(θ k ) by solving (34) at θ = θ k Compute W k as the solution of (37)

θ k+1 = θ k + ascent(J(θ k ), W k ) k ← k + 1 end while
acts as an external parameter. The field W defined in (37) can be interpreted as the adjoint state for that problem.

Collecting the results obtained leads to Algorithm 2 for solving the incremental thermo-mechanical problem (27). At each iteration, a mechanical problem at a fixed temperature field is solved for evaluating J, and a scalar problem is solved for determining an ascent direction W k . The results are used to feed an ascent algorithm for updating the temperature field. One of the simplest choices consists in performing a line search in the direction W k .

Algorithm 2 retains some attractive features of Algorithm 1, such as the decoupling between mechanical and thermal subproblems. The mechanical subproblems in both those algorithms have the same structure. In contrast, the thermal subproblem (37) in Algorithm 2 is linear, whereas the thermal problem (26.4-5) in Algorithm 1 generally is not. The proposed method has the advantage of relying on a sound mathematical framework, which bodes well for robustness and convergence properties.

Application to shape-memory alloys

The formulation derived so far is now applied to shape memory alloys. Such material indeed exhibit strong thermo-mechanical coupling (notably through latent heat effects) and therefore offer a relevant application for illustrating the proposed method. It is not the purpose of this paper to give a detailed presentation of shape-memory alloys (see e.g. the books by [START_REF] Otsuka | Shape memory materials[END_REF]; [START_REF] Bhattacharya | Microstructure of Martensite, oxford materials Edition[END_REF]). We simply mention that the peculiar properties of those materials stem from a solid/solid phase transformation between different crystallographic structures, known as austenite and martensite. The martensitic lattice has less symmetry than the austenitic one, which leads one to distinguish several martensitic variants, identified as individual phases.

Here we consider a micromechanical model of shape-memory alloys, for which the internal variable α = (α 1 , • • • , α n ) corresponds to the volume fractions of the n martensitic variants. Because of mass conservation in the phase-transformation process, the variable α must belong to the ndimensional tetrahedron T defined as

T = {α ∈ R n |α i ≥ 0 ∀i; n i=1 α i ≤ 1}. ( 38
)
We consider a Helmholtz free energy w and a dissipation potential Φ given by the following expressions (Abeyaratne et al., 1994;[START_REF] Govindjee | A multi-variant martensitic phase transformation model : formulation and numerical implementation[END_REF]Anand and Gurtin, 2003):

w(ε, α, θ) = 1 2 (ε - n i=1 α i ε tr i ) : L : (ε - n i=1 α i ε tr i ) + λ T θ T (θ -θ T ) n i=1 α i + c(θ -θ R -θ log( θ θ R )), (39)
Φ( α) = G + . α + + G -. α -. (40)
where x + denotes the positive vector whose component i is max(0, x i ). Similarly, for any vector x, x -is the positive vector with components max(0, -x i ). In (40), G + and G -are two given positive vectors of R n that characterize the mechanical dissipation in the model. In (39), the so-called transformation strains ε tr i are obtained from the crystallographic structure of the alloy considered. The elasticity tensor L is symmetric positive definite. The parameter λ T is the latent heat at the transformation temperature θ T , and c is the specific heat. We refer to [START_REF] Govindjee | A multi-variant martensitic phase transformation model : formulation and numerical implementation[END_REF]; [START_REF] Hackl | An upper bound to the free energy of n-variant polycrystalline shape memory alloys[END_REF]; [START_REF] Peigney | A non-convex lower bound on the effective free energy of polycrystalline shape memory alloys[END_REF]Peigney ( , 2013a,b) ,b) for more details and recent developments on micromechanical modelling of shape-memory alloys.

As detailed in Section 4, solving the thermo-mechanical problem (31) partly relies on solving the isothermal problem (21). For the material model considered here, the isothermal problem (21) specializes as

σ ∈ K σ , u ∈ K u , α ∈ K α ; σ = L : (ε - n i=1 α i ε tr i ); A = (ε tr 1 : σ, • • • , ε tr n : σ); A = A d + A r , A d ∈ ∂Φ( α -α 0 δt ) , A r ∈ ∂I T (α); ( 41
)
where the subdifferentials ∂Φ and ∂I T are characterized as follows:

A d = (A d 1 , • • • , A d n ) ∈ ∂Φ( α -α 0 δt ) ⇐⇒    A d i = G + i if α i > α i 0 = -G - i if α i < α i 0 ∈ [-G - i , G + i ] if α i = α i 0 and A r = (A r 1 , • • • , A r n ) ∈ ∂I T (α) ⇐⇒ there exists z ∈ R + and a i ∈ R + such that A r i = z -a i ; z(1 -n i=1 α i ); α i a i = 0.
Concerning the thermal equations, note from (39) that the entropy s is given by αi

) + [G + . α + + G -. α -] + r. ( 42
)
Apart from the external heat source r, the two terms on the right-hand side can be interpreted as internal heat sources due to phase transformation. They are respectively associated with the latent heat (reversible contribution) and with the mechanical dissipation (irreversible contribution). Depending on the rate of loading and on the thermal exchange conditions, the heat produced by the phase transformation may not have time to diffuse in the body and the temperature field may become inhomogeneous. In such conditions, the overall stress-strain response becomes significantly different from its isothermal counterpart. Such effects are explored in more detail in the following.

Influence of the strain rate

The first example we consider is inspired by the experiments of [START_REF] Shield | Orientation dependance of the pseudoelastic behaviour of single crystals of cu-al-ni in tension[END_REF]. We consider a monocrystalline CuAlNi strip with dimensions 38mm × 6mm × 1mm, submitted to displacement-controlled traction along the x-axis (Figure 1). The end section at x = 0 is clamped.

On the end section at x = L (L =38mm), the displacements along the e z and e y axes are set equal to 0 and the displacement u * (t) along the e x axis is prescribed as

u * (t) = t T u * max for 0 ≤ t ≤ T (2 -t T )u * max for T ≤ t ≤ 2T ( 43
)
where u * max and T are fixed. The loading strain rate ε * is defined as

ε * = | d dt ( u * L )| = u * max LT .
A null heat flux is prescribed on both end sections of the sample. A convection boundary condition of the form ( 19) is assumed on the lateral surface S lat . We consider two different thermal boundary conditions on the end sections, labeled as (i) and (ii). Boundary condition (i) consists in imposing a null heat flux on the end sections x ∈ {0, L}. Boundary condition (ii) consists in imposing θ = θ R on the two end sections. Those two boundary conditions can be interpreted as limiting cases of the convection condition ( 19), with h = 0 and h = ∞, respectively .

In the initial state, the structure is fully austenitic (i.e. α(x, 0) = 0 for all x) and in thermal equilibrium (i.e. θ(x, 0) = θ R ).

There are six martensitic variants in CuAlNi [START_REF] Otsuka | Shape memory materials[END_REF][START_REF] Bhattacharya | Comparison of the geometrically nonlinear and linear theories of martensitic transformation[END_REF]. The six transformation strains ε tr i , expressed in the (x, y, z) basis of Figure 1, take the form

ε tr i = R T .ε tr 0,i .R
where ε tr 0,i (i = 1, • • • , 6) are the reference transformation strains (expressed in the natural basis of the austenitic cubic lattice) and R is a rotation describing the orientation of the sample with respect to the austenitic lattice. The reference transformation strains ε tr 0,i are listed in Appendix C. The rotation R considered corresponds to the 'A1-T1b' sample in the nomenclature of [START_REF] Shield | Orientation dependance of the pseudoelastic behaviour of single crystals of cu-al-ni in tension[END_REF], and is given by

R =   0.925 0.380. 0 -0.380 0.925 0 0 0 1  
The thermo-mechanical response of the sample is obtained by solving the maximization problem (31) with a Broyden-Fletcher-Goldfarb-Shanno algorithm, using respectively the expressions (34) and (37) for evaluating J and an ascent direction. A finite-element method is used for discretizing (34) and (37) with respect to space. The mesh considered consists of 16×8×1 eight-node brick elements, as represented on Figure 1. The problem (37) is linear and does not present any substantial difficulty. The nonlinear problem (34) is more delicate, notably because of the constraint (38) on the internal variable α. We use a method proposed in [START_REF] Peigney | Numerical simulation of shape memory alloys structures using interior-point methods[END_REF] for solving (34) in the presence of such constraints. That method essentially consists in reformulating (34) as a linear complementarity problem (via a change of variables) and from there using a interior-point algorithm (see e.g. the books by Ye (1997) or [START_REF] Wright | Primal-dual interior-point methods[END_REF]) for a detailed presentation of linear complementarity problems and interior-point methods).

For CuAlNi, the results obtained with the boundary condition (i) and (ii) are qualitatively similar. The main difference is that fluctuations of the temperature are more pronounced for boundary condition (i). In the following we only present results corresponding to the boundary condition (i).

The evolution of the phase transformation during the loading is illustrated on Figure 1 (top), along with the evolution of the temperature field (bottom). Those numerical results correspond to h = 400W.m -2 .K -1 , θ R = 313K, ε * = 2.10 -3 s -1 . In the free energy (39), the elasticity tensor L is taken as isotropic with a Young's modulus equal to 26.7Gpa and a Poisson's ration equal to 0.25 [START_REF] Shield | Orientation dependance of the pseudoelastic behaviour of single crystals of cu-al-ni in tension[END_REF][START_REF] Govindjee | A multi-variant martensitic phase transformation model : formulation and numerical implementation[END_REF]. The transformation temperature θ T and the latent heat λ T are set equal to 277K and 46.7Mpa, respectively. Those values have been measured by [START_REF] Shield | Orientation dependance of the pseudoelastic behaviour of single crystals of cu-al-ni in tension[END_REF] using differential scanning calorimetry. The specific heat c is taken as 3.1Mpa [START_REF] Otsuka | Shape memory materials[END_REF]. The components of the dissipative parameters G + and G -take a common value G that is identified from Shield's experiments and is equal to 0.15Mpa (see [START_REF] Peigney | Numerical simulation of shape memory alloys structures using interior-point methods[END_REF] for more details).

As can be observed on Figure 1, the phase transformation initiates at the middle section of the sample, and propagates towards the end sections. In order to explain that behaviour, note from (41) that the condition for (austenite to martensite) phase transformation to occur at point x and time t reads min i σ(x, t) :

ε tr i = G + λ θ T (θ -θ T ). (44)
For small time t, the evolution is elastic (no phase transformation) and the temperature is equal to θ R at all point. In such case, the stress field σ(x, t) can be written as (t/T )σ E (x) where σ E (x) is solution of the following linear elasticity problem: Because of the clamped boundary conditions at x = 0 and x = L, the stress field σ E (x) is not homogeneous and exhibits some stress concentration near the end sections. From (44), phase transformation initiates at point x * such that min

div σ E = 0 , σ E .n = 0 on S lat ; σ E = L : (∇u E + T ∇u E )/2; u E = 0 on x = 0, u E = u max e x on x = L. (45)
i σ E (x * ) : ε tr i = sup x∈Ω (min i σ E (x) : ε tr i ). (46)
In the case of CuAlNi, that condition is met at the middle section, even though the Von Mises stress at that location is lower than near the end sections.

The thermo-mechanical coupling results in a heterogeneous temperature field, the phase transformation front acting as a moving heat source in accordance with (42). The variations of the temperature with respect to the am-bient temperature θ R get more pronounced when the strain rate ε * increases, as illustrated on Figure 2. On that Figure are represented the maximum and minimum temperatures in the sample at t = T , i.e. when the applied strain is maximum. Variations of the temperature are of the order of 10K for high strain rates, which is consistent with the order of magnitude observed experimentally [START_REF] Grabe | On the viscous and strain rate dependant behavior of polycrystalline niti[END_REF]. Observe on Figure 2 that θ max and θ min converge towards a limit as ε * → +∞. That behaviour can be explaining by noting that, for high values of ε * , conduction and convection effects become negligible compared to the latent heat effect. In such case, the local temperature θ(x, t) at point x and time t is directly correlated to the evolution of α(x, t) at the same point. More precisely, using the simplifying assumption G (θ/θ T )λ T , the integration of the heat equation ( 42) with respect to t gives

θ(x, t) θ R exp( λ T cθ T ( n i=1 α i (x, t)). ( 47
)
The temperature is thus maximum at points x where the material is fully transformed in martensite, i.e. at points where i α i (x, t) = 1. It follows that

θ max = θ R exp( λ T cθ T ). ( 48
)
For the material parameters of CuAlNi, the formula (48) gives θ max = 330.7K. The simulation results for ε = 0.5 gives θ max = 330.6K, which is in good agreement with (48). A formula analog to (48) holds for θ min : using (47) we have

θ min = θ R exp( λ T cθ T min x∈Ω ( n i=1 α i (x, T ))). (49)
However, at time t = T , phase transformation is initiated at all point in the sample, i.e. i α i (x, T ) > 0 for all x. The minimum value min x i α i (x, T ) is non zero and cannot be obtained in closed form. For ε * = 0.5, the numerical simulation gives min x i α i (x, T ) = 0.163 and θ min = 315.9K. Those values are compatible with (47).

On Figure 3 are represented the stress-strain curves obtained for several values of the applied strain rate ε * . The curve obtained for a very low strain rate ( ε * = 5.10 -4 s -1 ) coincides with the isothermal simulations presented in [START_REF] Peigney | Numerical simulation of shape memory alloys structures using interior-point methods[END_REF] and exhibits some distinctive plateaux in the stress response. As the applied strain rate ε * is increased, some hardening of the stress-strain curve is observed. As can be noticed on Figure 3, the hardening is strain-rate dependent. That behaviour is a direct consequence of thermal effects. Let us indeed emphasize that the mechanical behaviour of the model is rate-independent: the only parameters that introduce a time scale are those related to thermal effects (thermal conduction K, convection coefficient h, specific heat c).

The size of the hysteresis loop (here denoted by D) is an important parameter in some applications of shape-memory alloys, such as damping applications. There are contradictory experimental observations in the literature concerning the variation of D with respect to the applied strain-rate ε * . Some authors, as [START_REF] Shaw | On the nucleation and propagation of phase transformation fronts in a niti alloy[END_REF] observe that D increases with the applied strain rate, whereas others, as [START_REF] Dolce | Mechanical behaviour of sma for seismic applications. 2: austenite niti wires subjects in tension[END_REF], observe the opposite. In the present model, D varies in a non monotonic fashion with ε * , as can be noted on Figure 3. This is shown more clearly on Figure

Influence of the thermal boundary conditions

We now consider an example related to NiTi alloys. As detailed in Appendix C, there are twelve transformation strains to be considered for that material. The geometry and the mechanical boundary conditions are identical to those considered previously for CuAlNi, except that the dimensions of the sample are now 15mm×2mm×0.3mm.

On Figure 5 is presented the phase transformation evolution obtained from the two types of boundary conditions (i) and (ii). The thermal boundary condition (ii) results in a clearer localization of the phase transformation: for sufficiently high values of ε, phase transformation initiates at the end sections and propagates towards the middle section. That last scenario is in agreement with the numerical simulations performed by Anand and Gurtin (2003) for the same material and geometry, with a boundary condition of type (ii).

Recall that phase transformation initiates at points x satisfying the condition (46). That condition strongly depends on the transformation strains in the alloy considered. In the case of NiTi, the condition ( 46) is satisfied near the end sections, in contrast with the CuAlNi example considered previously. The time t * when phase transformation begins is given by

t * T min i σ E (x * ) : ε tr i = G. (50)
Note that ( 46) and ( 50) are independent on the thermal boundary conditions. The initiation of phase-transformation is thus identical for both boundary conditions (i) and (ii). As soon as phase transformation begins, heat is produced and the temperature rises locally. This is when the thermal boundary condition matters. Observe from (44) that high temperatures defavour phase transformation. In the vicinity of the end sections, the adiabatic boundary condition (i) tends to keep the temperature high compared to the fixed temperature boundary condition (ii). In such condition, phase transformation becomes favoured near the middle sections, where the temperature remains close to θ R . On the contrary, compared to the condition q = 0 used in (i), the boundary condition (ii) helps in preventing significant rise of the temperature near the end sections and therefore promotes phase transformation in those points. The strain-stress curves obtained for the boundary conditions (i) and (ii) are represented on Figure 6. Although the loading strain rate ( ε * = 10 -3 s -1 ) is the same for both simulations, the obtained stress-strain curves are dramatically different. The thermal boundary condition (i) results in a larger hardening of the response compared to (ii). Also note that the stress-strain curve for boundary condition (ii) is linear for ε > 4%, which indicates that phase transformation in martensite is complete for such level of strains. For both boundary conditions, the size D of the hysteresis loop is bigger than for the isothermal case (represented as dashed lines on Figure 6). In accordance with (46-50), the phase transformation begins for the same applied strain (about 0.76%) on the three responses.

The results of Figures 56show that phase transformation may strongly depend on the thermal boundary condition. Such an effect is not observed on the CuAlNi example considered previously. The reason is that the phase transformation in CuAlNi initiates at the middle section (see Figure 1), so that the thermal boundary conditions at the end sections do not have as strong an impact.

Concluding remarks

In this paper have been studied some incremental variational principles for the thermo-mechanical problem that results from the combination of generalized standard materials, non-smooth mechanics, and Fourier's law. Compared to the isothermal case, existence of solution is obtained under the additional requirement that the Helmholtz free energy w is concave with respect to the temperature. Building on those incremental variational principles, it has been shown that the incremental thermo-mechanical problem could be recast as a concave maximization problem. Although other routes are possible for solving the incremental problem, that particular way has the attractive feature of being simple to implement. Indeed the problem ultimately reduces to a sequence of linear scalar problems and purely mechanical problems, defined respectively by ( 34) and (37). Note that the structure of the linear scalar problem (37) to solve is independent on the material model considered (i.e. on the expression of the free energy w and of the dissipation potential Φ). All the specificities of the material model appear in the purely mechanical problem (34), which has the same structure as the incremental problem supplied by the Euler implicit scheme in the isothermal case. This is also where all the difficulties are concentrated, as nonlinearities, non-differentiabilities and constraints on the state variables may need to be handled. Provided one can solve that isothermal problem, the thermo-mechanical problem can be solved with low additional effort in terms of numerical implementation. For the micromechanical model of shape-memory alloys considered in Section 5, the isothermal algorithm proposed by [START_REF] Peigney | Numerical simulation of shape memory alloys structures using interior-point methods[END_REF] acted as a building block for the thermo-mechanical simulations. Such simulations are important for the design of SMA-based systems, as instance for estimating the response-time of SMA actuators, or for assessing the energy absorption capability of SMA dampers. Shape-memory alloys are obviously just an example of possible applications of the proposed approach, which for instance could be used for a wide class of plasticity models. In this paper, some mathematical aspects have been kept at a formal level. For instance, the choice of functional spaces for the different fields considered has not been discussed in detail. It would be interesting to explore those aspects in a more thorough full fashion, notably in order to study the convergence of solutions as the time increment tends towards zero.

where the expression of (∂F θ /∂θ). θ is given in (A.2). Since F d is linear with respect to θ, we obtain For CuAlNi, values of the lattice parameters are α = 0.0425, β = -0.0822, δ = 0.0194 [START_REF] Chu | Hysteresis and microstructures: a study of biaxial loading on compound twins of copperaluminumnickel single crystals[END_REF].

In the case of NiTi, austenite and martensite respectively have a cubic and monoclinic-I structure. There are 12 martensitic variants, with transformation strains given by: For nearly equiatomic NiTi alloys, the values of the lattice parameters are α = 0.0243, β = -0.0437, δ = 0.058, = 0.0427 [START_REF] Knowles | Crystallography of the martensitic transformation in equiatomic nickel-titanium[END_REF]. Abeyaratne, R., Kim, S., Knowles, J., 1994

Figure 1 :

 1 Figure 1: Phase transformation (top) and temperature distribution (bottom) in a CuAlNi rectangular sample loaded in traction.

Figure 2 :

 2 Figure 2: Maximum and minimum temperatures in the sample as a function of the loading strain rate.

Figure 3 :

 3 Figure 3: Stress-strain curves at several applied strain rates (CuAlNi).

  4: D increases for low value of ε * and then decreases after having reached a maximum. The maximum value of D as well as the corresponding value of ε * are not intrinsic properties of the material: they depend on the surrounding medium (notably through the convection coefficient h) as well as on the shape and dimensions of the sample. Provided such a behaviour is not just an artefact of the model and corresponds to a real phenomenon, it would explain the seemingly contradictory observations found in the literature.

Figure 4 :

 4 Figure 4: Size of the hysteresis loop as a function of the loading strain rate.

Figure 5 :

 5 Figure 5: Phase transformation in NiTi, for two different thermal boundary conditions at the end sections: null heat flux (top) and prescribed temperature (bottom)

Figure 6 :

 6 Figure 6: Stress-strain curves in NiTi for different thermal conditions

F

  d (α(t), θ(t)) -F d (α(t), .4-B.6) in (B.3), we obtain DJ(θ, θ) ≥ DF θ (θ, θ) + Ω (-s + A d .(α -α 0 )) θ dω (B.7) where the definition s = -∂w/∂θ and the relation Φ(α -α 0 ) = A d .(α -α 0 ) have been used. Comparing with (A.5) shows that the right-hand side of (B.7) is equal to (∂F/∂θ). θ. We thus obtain DJ(θ, θ) ≥ ∂F ∂θ . θ Appendix C. Transformation strains in CuAlNi and NiTiThe CuAlNi alloy obeys a cubic to orthorombic transformation, for which there are 6 (lattice correspondent) martensitic variants with reference transformation

  . A one-dimensional continuum model for shape-memory alloys. Int. J. Solids Struct. 31, 2229-2249. Anand, L., Gurtin, M., 2003. Thermal effects in the superelasticity of crystalline shape-memory materials. J. Mech. Phys. Solids 51 (6), 1015-1058. Auricchio, F., Petrini, L., 2004. A three-dimensional model describing stresstemperature induced solid phase transformations, part ii : thermomechanical coupling and hybrid composite applications. Int. J. Num. Meth. Eng. 61, 716-737.

In non-smooth mechanics, left-and right-time derivative of physical quantities may not be equal. In order to respect the principle of causality, the constitutive relations need to be written in terms of left-time derivatives (see e.g.[START_REF] Frémond | Non Smooth Thermo-mechanics[END_REF]).

Appendix A. Derivation of the variational principle for the incremental thermo-mechanical problem

In this section, we establish that solutions of the variational problem (29) are solutions of the incremental thermo-mechanical problem (27). Studying the stationarity conditions with respect to u and α in is similar to the isothermal case (29) and leads to (27.1-3). In the following we focus on the stationarity condition with respect to θ. For any θ such that θ = 0 on Γ θ , we have

(A.3) Moreover, since Φ is positively homogeneous of degree 1 and

(A.5) Therefore, the stationarity condition ∂F/∂θ = 0 give the relations (27.4-5).

Appendix B. Directional derivative of J

Assuming that w is convex in (u, α) and concave in θ, let us establish that the function J in (32) is concave. A first observation is that F θ , F e as well as F d are concave with respect to θ. In (32), the term min (u,α) {F e (u, α, θ) + F d (α, θ)} is the minimum of a family of concave functions (parametrized by (u, α)), and therefore is concave in θ. As a result, the function J in (32) is concave. The concavity of J implies the existence of the directional derivative

where θ(t) = θ + t θ and θ vanished on Γ θ . We prove in the following that DJ(θ; θ) satisfied the inequality (35). To that purpose, let us denote by (u(t), α(t)) the solution of the minimization problem min

From the definition of (u(t), α(t)) we obtain

3) We now examine the limits (as t → 0 + ) of the three terms in brackets in the right-hand side of (B.3). The expression (22) of F e gives F e (u(t), α(t), θ(t)) -F e (u(t), α(t), θ) t → where ∂w/∂θ is evaluated at (u(0), α(0), θ). Moreover, we have


Luc Andrea
email: luc.andrea@onera.fr

Gilles Hug

Laurent Chaput

Ab initio lattice thermal conductivity in pure and doped Half-Heusler thermoelectric materials

Half-Heusler phases are promising intermetallics for applications in thermoelectric generators. Optimization of their thermal transport properties is essential to improve their overall conversion efficiency. Our goal is to perform a theoretical evaluation of thermal transport properties directly from first-principles calculations for various pure and doped half-Heusler compounds. The electronic structures are modeled in the framework of the density functional theory (DFT). The ab initio thermal properties are deduced from harmonic and anharmonic interatomic force constants calculations using finite size displacements method. Many-body perturbation theory is used for calculating the phonon-phonon interactions which yields the lifetime of phonons as function of momentum and band index. Finally, thanks to a direct solution to the phonons Boltzmann transport equation we computed the ab initio thermal conductivities, which are found in good agreement with the experimental data.

Introduction

Improving the overall efficiency of motors and power generators is a prerequisite to evolve toward a more environmental friendly and sustainable world. The development of new thermoelectric materials (TE) for converting wasted heat into usable electric power can help this goal [START_REF] Snyder | Complex thermoelectric materials[END_REF] (see figure 1). Equally important, the availability of these advanced technologies in a medium-long term is an obligation for maintaining the European competitiveness in the aerospace industry. Identifying high efficiency TE materials is a challenge residing in the combination of the physical properties involved. Best TE should have a phonon-glass electron-crystal behaviour. Phononglass describe the necessity for a low lattice thermal conductivity as in amorphous glasses, and electron-crystal because the best electrical properties (Seebeck coefficient and electrical conductivity) have been found in crystalline semiconductors. The recent enhancement of the field of thermoelectric is mainly due to the discovering of new methods that can synthesise materials close to the ideal phonon-glass electron-crystal structure. Half-Heusler (HH) compounds is a class of materials that can help those requirements [START_REF] Poon | Electronic and Thermoelectric Properties of half-Heusler Alloys[END_REF][START_REF] Kim | High temperature thermoelectric properties of TiNiSn-based "Half-Heusler" compounds[END_REF][START_REF] Takeyuki | High-Thermoelectric Figure of Merit Realized in p-Type Half-Heusler Compounds: ZrCoSnxSb1-x. 2007 Japan[END_REF]. Furthermore, these materials are more specific to the aerospace industry which requires high temperature materials that can sustain aggressive environments (mechanical stress, corrosion, oxidation). The main problem encountered is the phonon-glass behaviour. Indeed, pure HH possess a high thermal conductivity that thermalise the system too fast. Recently, many efforts are deployed to modify the basic HH structure to reduce the thermal conductivity without affecting the electrical transport properties. Experimentally, several routes are explored to achieve this goal. One route consist in doping the material with isoelectronic substitutions [START_REF] Culp | Effect of substitutions on the thermoelectric figure of merit of "Half-Heusler" phases at 800 • C[END_REF]. It generate mass fluctuations helping the diffusion of heat trough the material. This doping is widely used experimentally but, only few theoretical work is available. Our goal is to give a clear theoretical description of the mechanisms involved in the calculation of the thermal conductivity in order to find new pathways to improve the TE efficiency. In this document, we present the basic method to model the thermal transport properties in such materials. It consists in studying the phonons (see figure 2) and their interactions which are the main actors in the thermal transport. The lattice thermal conductivity has been calculated for various pure undoped HH structures, including ZrCoSb which is in good agreement with the experimental work. A method to account for isoelectronic doping is presented and preliminary results applied on one HH structures are shown.

Ab initio modelling of the lattice thermal conductivity

Pure structures

To compute the thermal conductivity, the standard approach consists in solving the Boltzmann transport equation (BTE) describing the evolution of the lattice vibrations (the phonons) out of equilibrium, for instance when a temperature gradient is applied. The Boltzmann equation represent the evolution of distribution function which measure the phonon occupation number, and is written : When a temperature gradient is applied to a metal or a semiconductor, the charge carriers (electrons e -and holes h + , respectively for n-type and p-type semiconductors) are diffusing from the hot to the cold side. This produces a charge difference inducing an electrostatic potential between the two sides. The system evolve towards an equilibrium between the electrical repulsion of charge carriers and the chemical potential. This is known as the Seebeck effect. The produced voltage V is proportional to the temperature difference ∆T between the two sides, where the proportionality constant α is the Seebeck coefficient : V = α∆T . To maximise the Seebeck coefficient, single type of carriers is needed. Indeed, if both type of carriers are present (mixed p-type and ntype) and moving to the cold end, the Seebeck coefficient cancels out. The efficiency of TE materials is characterised by their figure of merit ZT : ZT = σS 2 T κ where σ is the electrical conductivity, S the Seebeck coefficient, T the temperature and κ the thermal conductivity. To maximise ZT , large electrical conductivity, low thermal conductivity and high Seebeck coefficient are required. Since these transport properties rely on interdepend material properties, a set of parameters need to be optimised. Real TE device are a collection of n-type and p-type semiconductors (not mixed) assembled in parallel. This figure shows a setting where both type of semiconductors are used to maximise the TE efficiency.

∂n (1) qp ∂t + ∂n (0) qp ∂T ∂T ∂ r • v qp = C(qp; n (1) qp ) + 1 2 D(qp; n (1) qp ), Thermoelectricity e - e - e - e - e - e - e - e - e - e - e - h + h

Phonons

Figure 2: In a gas, every molecules behave almost independently from each other. The temperature effect reside in the velocities of the molecules that are allowed to move freely and collide with each other. In a crystalline solid, atoms are localised around specific positions and bonded together. To retrieve the effect of the temperature, the same treatment as a for the gases is no more valid. In this case, the phonon description is used. Phonons are quantised vibrations of the atomic lattice. Instead of moving randomly, the atoms are oscillating coherently around their equilibrium positions. Those oscillations, the phonons, are carrying heat and can propagate through the material. This figure shows a linear chain of atoms where the atoms are allowed to move in the perpendicular direction. The first drawing show the atoms at their equilibrium positions. The second drawing show the oscillation of the atoms around their equilibrium positions. The red arrows show the displacement of the atoms from their equilibrium positions. For a given crystalline structure, thermal properties are retrieved via analysing the phonons and their collisions.

where

n qp ≈ n (0) qp +n (1)
qp is the phonon occupation function for a phonon with wave vector q in the branch p. n (0) qp is the equilibrium value and n [START_REF] Snyder | Complex thermoelectric materials[END_REF] qp is the first order deviation from equilibrium. T is the temperature, r the position, v qp the velocity. C(qp; n qp ) for the decay processes, both depending on the phonon lifetime τ qp . The lifetime is driven by the amount of interactions of the phonons. The solution of the BTE is usually found either by applying an iteration scheme [START_REF] Omini | An iterative approach to the phonon Boltzmann equation in the theory of thermal conductivity[END_REF] or in the relaxation time approximation which consists in fixing a constant lifetime. In this work, to retrieve the ab initio thermal conductivity, we aim at overcoming these two approximations with a direct solution to the phonon BTE [START_REF] Chaput | A direct solution to the phonon Boltzmann equation[END_REF] where the main input, the phonons lifetimes, were deduced thanks to perturbation theory and ab initio density functional theory (DFT). The phonons have a wide panel of interactions driving their lifetimes. Indeed, they can interact strongly with other phonons, with isotopes, boundaries, electrons, photons etc [START_REF] Bhattacharya | Effect of boundary scattering on the thermal conductivity of TiNiSn-based half-Heusler alloys[END_REF][START_REF] Zhu | Reduced Grain Size and Improved Thermoelectric Properties of Melt Spun (Hf,Zr)NiSn Half-Heusler Alloys[END_REF][START_REF] Chaput | Phonon-phonon interactions in transition metals[END_REF]. For the calculation of pure structures, we have considered interactions with isotopes, boundaries and between phonons. Every interactions add a contribution to the lifetime τ pure of the pure structure, this is encoded thanks to the Matthiessen's rule:

1 τ pure = i 1 τ i = 1 τ phonon/phonon + 1 τ isotopic + 1 τ boundary ,
where the phonon/phonon interaction lifetime term τ phonon/phonon is dominating over isotopic and boundary contribution τ isotopic , τ boundary at high temperature, while the boundary scattering dominate in the low temperature regime. To access to the phonon/phonon interaction contribution, the phonons are defined in the periodic electronic potential of the lattice. The Taylor expansion of the electronic potential E with respect to the displacements η from the equilibrium position gives,

E = E 0 + l α ∂E ∂η α (l) η α (l) + 1 2 ll αβ ∂ 2 E ∂η α (l)∂η β (l ) η α (l)η β (l ) + 1 6 ll l αβγ ∂ 3 E ∂η α (l)∂η β (l )∂η γ (l ) η α (l)η β (l )η γ (l ) + O(η 4 )
The sums over l, l and l run over all the lattice atomic positions and the sums over α, β and γ run over the spatial directions x, y and z. The first order vanishes because the development is made around the equilibrium positions. Expanded up to second order, the potential energy give rise to the harmonic contribution, which assume the role of a collection of independent quantum harmonic oscillators, the phonons. This allow us to compute basic thermodynamic quantities, such as constant volume specific heat or entropy. However, this approach is not sufficient to access to the thermal conductivity. Higher order expansions are needed to introduce an anharmonic contribution representing the phonon-phonon interactions. We use density functional theory to retrieve the harmonic and anharmonic contribution of the interatomic force constants with high accuracy. Harmonic computations are reasonably available without the necessity of high performance computing, the bottleneck resides in the anharmonic computations. Indeed, third order development of the potential energy drastically increase the total calculation time. From the DFT calculations, we can deduce the total lifetime τ total and solve the BTE. When the solution to the BTE is found, the total lattice conductivity κ can be deduced by,

κ = C p v 2 τ
Where C p is the heat capacity, v 2 the squared velocities and τ the relaxation time.

The presented method works well and do not demand excessive computing power if the system studied is highly symmetric. It has been applied to many well known structures such as silicium, diamond etc. The ab initio lattice thermal conductivity retrieved fits well with the experimental results.

Doped structures

When a structure is doped, the symmetry of the crystal changes. Most of the time, inserting defects causes a symmetry breaking where there are less symmetry operations that can be used to simplify the equations. A colossal computing power would be needed to access the thermal conductivity. Based on perturbation theory, we developed a method that can retrieve the effect of localised defects on the lattice thermal conductivity without high computing power. It consist in adding two perturbative terms to the hamiltonian,

H = H 0 + jl α p 2 α (jl) 2 δ 1 M (jl) + 1 2 jlj l αβ δA αβ (jl; j l )η α (jl)η β (j l )
Where H 0 is the unperturbed hamiltonian, the second term stands for mass fluctuations induced by the doping and the third term stands for the local fluctuations of the force constants. Those two contributions participate in reducing the phonons lifetimes and therefore also reducing the total lattice conductivity. The same way as before, both mass and force constant fluctuations add a contribution to the total lifetime 1 τ δmass and

1 τ δf orce , 1 τ total = i 1 τ i = 1 τ pure + 1 τ δmass + 1 τ δf orce ,
Those supplementary contributions diminish the total lifetimes and therefore reduce the lattice thermal conductivity. This has the overall effect of increasing the TE efficiency.

Results

Pure half-Heusler

The method described here above has been applied to few pure HH structures, including ZrCoSb which is widely studied. The lattice thermal conductivity with respect to the temperature is plotted in figure 3 for the ZrCoSb HH structure. The observed shape and value of κ is in well agreements with the experimental data. The phonon band structure decorated with the inverse phonon lifetimes is also shown figure 3. For the three lower phonon branches, we observe that the lifetimes are increasing. From this observation one can deduce that those phonons are the main actors in the thermal conduction. Indeed, if the lifetime of a phonon is large, its mean free path is also large. The phonon can travel long distances in the material and therefore easily transfer heat. We can also see that the three higher phonon branches possesses small lifetimes and therefore have a less important contribution to κ. Thanks to our method we are able to compute the spectral representation of all quantities involved in the calculations of the lattice thermal conductivity. This gives us useful informations about the most participating frequencies. Figure 4 shows the spectral lattice thermal conductivity as a function of frequency and temperature for the HH ZrCoSb. It shows the contribution of all frequencies. From this graph we can observe that the most participating range of frequency is located between 1 and 4 THz. This also confirm the high conductivity of the three lower phonon branches which are low frequency phonons. From this graph, one can also see a gap where the spectral lattice thermal conductivity is vanishing. It shows us that certain frequencies within that gap are not participating to the lattice thermal conductivity.

Doped half-Heusler, Ni 4 Ti 3 ZrSn 4

We applied our method to a doped HH structure. Compared to the previous calculations for pure structures, we added the contribution of the defects trough their mass and force constant fluctuations. As a first step, we studied the effects of an isoelectronic substitution to an HH structure. The isoelectronic substitution have the advantage of not disturbing much the electronic structure, in order to not deteriorate the electronic properties. Figure 5 shows the contribution of the mass fluctuations, the force constant fluctuations in the lattice conductivity compared to the pure structure. We can observe that the main contribution is due to the mass fluctuations. Which is coherent 0.01 with our description since the isoelectronic substitution do have a light influence on force constants.

Conclusion

We presented the main method to calculate the thermal properties of a crystalline structure. It has been applied to one HH structure, ZrCoSb where good agreements with the experimental work. Useful informations on frequencies contributions have been retrieved thanks the spectral representation. We also presented our technic that account for isoelectronic substitutions trough mass and force constant fluctuations. We observed that the mass fluctuation is the main effect that reduce the lattice thermal conductivity and therefore increase the overall thermoelectric efficiency for the HH Ni 4 Ti 3 ZrSn 4 structure. The next steps of our study will be concentrated on the effect of disorder and nanostructure inclusions on the lattice thermal conductivity.

Figure 1 :

 1 Figure 1: When a temperature gradient is applied to a metal or a semiconductor, the

( 1 )

 1 qp ) account for the collision processes and D(qp; n[START_REF] Snyder | Complex thermoelectric materials[END_REF]

Figure 3 :

 3 Figure 3: The first plot shows the lattice thermal conductivity of the HH ZrCoSb as a function of temperature. The phonon band structure decorated by the inverse lifetime is shown on the second plot. The error between experimental data and our method is below 5%.

Figure 4 :

 4 Figure 4: Spectral lattice thermal conductivity of the HH ZrCoSb. The z axis represents the effect of the temperature. The colours stands for differentiating the temperature.

Figure 5 :

 5 Figure 5: Lattice thermal conductivity for a doped HH structure Ni 4 Ti 3 ZrSn 4 . Every contribution are plotted separately and jointly.


A Pierret

L Schué

F Fossard

S Moldovan

O Ersen

F Ducastelle

J Barjon

A Loiseau

Tem

TEM and cathodoluminescence with nanometric spatial resolution on BN nanostructures

N/A

TEM and cathodoluminescence with nanometric spatial resolution on BN nanostructures

Pierret A. 1, 2 , Schué L. 1, 3 , Fossard F. 1 , Moldovan S. 4 , Ersen O. 4 , Ducastelle F. 1 , Barjon J. 3 , Loiseau A. 1 1 LEM, ONERA-CNRS, Châtillon, France, 2 CEA-CNRS-UJF group , 3 GEMaC, UVSQ-CNRS, Versailles, France, 4 IPCMS, Univ. Strasbourg-CNRS, Strasbourg, France Email of the presenting author: annick.loiseau@onera.fr Hexagonal boron nitride (h-BN) is a wide band gap semiconductor (6.4eV), which can be synthesized, as its carbon analog graphite, as bulk crystallites, nanotubes and nanosheets.

Investigation of their optoelectronic properties is made difficult because of the paucity of high quality samples and suitable investigation tools. These structures meet nevertheless a growing interest for deep UV LED and graphene engineering. A deeper understanding of the interplay between the structural and luminescence properties of different BN structures and how these properties can be further exploited for their characterization are therefore highly needed.

Such studies are now possible thanks to the recent development of dedicated photoluminescence (PL) and cathodoluminescence (CL) experiments running at 4K and adapted to the detection in the far UV range (up to 6eV) [1]. We can also combine various TEM techniques and CL experiments in a FEG-SEM with a spatial resolution of 3nm on the same nano-object. With these tools, we investigated the structure and luminescence of various structures, from high quality crystals [2], exfoliated nanosheets to multi-wall nanotubes [3].

As a result, BN materials present original optical properties, governed by excitonic effects in the 5.5-6eV energy range. Two kinds of excitonic luminescence have been identified and are called S and D lines [4]. As revealed from CL-TEM analyses, D lines are issued from defective areas (Fig 1), so that D/S ratio can be used as a qualification parameter of the defect densitiy [5]. This procedure has been applied to understand the first luminescence studies of few layers individual BN flakes [5].

Concerning nanotubes, CL images reveal that the luminescence in the 5.5-6eV energy range is strongly inhomogenous and oscillating. Thanks to a deep investigation combining different TEM techniques, we have shown that the tubes display a complex twisted faceted structure and that the twist period is correlated with the luminescence oscillations (Fig 2). Furthermore, we could show that excitons, responsible for the spectacular localization of the luminescence, are trapped to specific defects, twisted along with the faceting structure.

Finally, low-loss EELS providing an alternative approach to the nature of electronic excitations [6], we will show how it is an efficient tool to investigate the local structure and optical properties with an energy resolution below 100meV of different BN layers and nanotubes.

Fig. 1 :

 1 Fig. 1: (a) SEM image of a h-BN crystallite; (b), (c) Corresponding CL images recorded (b) on the main S line (S3-S4), and (c) on the main D line (D4). (d) Map of the D/S ratio. (e) CL spectra recorded in the areas #1 (grain boundary) and #2 (middle of the grain), indicated in (a).

Fig. 2 :

 2 Fig. 2: Images of a BN nanotube: (a) 3nm spatially resolved CL image recorded at 5.49 eV (226 nm); (b), (c) Corresponding TEM images in (b) bright-field mode, and (c) dark-field mode on the (100) reflection. (d) Heptagonal tube cross-section obtained by tomography experiment. (e) Structure of the tube as deduced from (b-d) images.

Acknowledgement: The research leading to these results has received funding from the European Union Seventh Framework Programme under grant agreement n°604391 Graphene Flagship. D. Golberg, T. Taniguchi and K. Watanabe from NIMS Japan are warmly acknowledged for providing samples.