Gang Xu

Bernard Mourrain
email: bernard.mourrain@inria.fr

Régis Duvigneau
email: regis.duvigneau@inria.fr

André Galligo
email: galligo@unice.fr

Constructing analysis-suitable parameterization of computational domain from CAD boundary by variational harmonic method

Keywords: isogeometric analysis, grid generation, harmonic mapping, variational method, analysis-suitable parameterization

In isogeometric anlaysis, parameterization of computational domain has great effects as mesh generation in finite element analysis. In this paper, based on the concept of harmonic mapping from the computational domain to parametric domain, a variational harmonic approach is proposed to construct analysis-suitable parameterization of computational domain from CAD boundary for 2D and 3D isogeometric applications. Different from the previous elliptic mesh generation method in finite element analysis, the proposed method focus on isogeometric version, and converts the elliptic PDE into a nonlinear optimization problem, in which a regular term is integrated into the optimization formulation to achieve more uniform and orthogonal isoparametric structure near convex (concave) parts of the boundary. Several examples are presented to show the efficiency of the proposed method in 2D and 3D isogeometric analysis.

Introduction

The isogeometric analysis (IGA for short) method proposed by Hughes et al. in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement[END_REF] is a new computational approach that offers the possibility of seamless integration between CAD and CAE. The method uses the same type of mathematical representation (spline representation), both for the geometry and for the physical solutions, and thus avoids this costly forth and back transformations. Moreover its reduces the number of parameters needed to describe the geometry, which is of particular interest for shape optimization.

Mesh generation, which generates a discrete grid of a computational domain from a given CAD object, is a key and the most time-consuming step in finite element analysis (FEA for short). It consumes about 80% of the overall design and analysis process [START_REF] Bazilevs | Isogeometric analysis using T-Splines[END_REF] in automotive, aerospace and ship industry. Parametrization of computational domain in IGA, which corresponds to the mesh generation in FEA, also has some impact on analysis result and efficiency. In particular, arbitrary refinements can be performed on the computational mesh in FEA, but in IGA if we compute with tensor product B-splines, we can only perform refinement operations in each parametric direction by knot insertion or degree elevation. Hence, parameterization of computational domain is also being important for IGA. As it is pointed by Cottrell et al. [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF], one of the most significant challenges towards isogeometric analysis is constructing analysis-suitable parameterizations from given CAD boundary representation.

In IGA, the parameterization of a computational domain is determined by control points, knot vectors and the degrees of B-spline objects. For 2D and 3D IGA problems, the knot vectors and the degree of the computational domain are determined by the given boundary curves/surfaces. That is, given boundary curves/surfaces, the quality of parameterization of computational domain is determined by the positions of inner control points. Hence, finding a good placement of the inner control points inside the computational domain, is a key issue. As the required mesh quality in FEA [START_REF] Wang | Non-self-overlapping Hermite interpolation mapping: a practical solution for structured quadrilateral meshing[END_REF][29] [START_REF] Wang | Algebraic grid generation on trimmed surface using non-self-overlapping Coons patch mapping[END_REF], a basic requirement of the resulting parameterization for IGA is that it doesn't have self-intersections, so that it is an injective map from the parametrization domain to the computational domain.

In the field of mesh generation in FEA, a general method is based on partial differential equations [START_REF] Spekreijse | Elliptic grid generation based on Laplace equations and algebraic transformations[END_REF]. The grid points are the solution of an elliptic partial differential equation system with Dirichlet boundary conditions on all boundaries. There are several advantages in elliptic mesh generation. The theory of partial differential equations guarantees that the mapping between physical and transformed regions will be one-to-one. Another important property is the inherent smoothness in the solution of elliptic systems. A disadvantage of elliptic method is that there will be some non-uniform grid elements near convex (concave) parts of the boundary. Moreover, special treatment should be considered to achieve orthogonal mesh [START_REF] Zhang | 2D nearly orthogonal mesh generation with controls on distortion function[END_REF][41] [START_REF] Wu | Orthogonal Grid Generation of an Irregular Region Using a Local Polynomial Collocation Method[END_REF].

Motivated from the elliptic grid generation method in FEA, in this paper, we study the analysissuitable parameterization problem of computational domain in isogeometric analysis. Our main contributions are:

• A nonlinear optimization framework based on harmonic mapping theory is proposed to generate a analysis-suitable parameterization of computational domain without self-intersections for 2D and 3D isogeometric applications.

• A regular term is integrated into the optimization formulation to achieve more uniform and orthogonal isoparametric structure near convex (concave) parts of the boundary.

• We test the parameterization results on 2D and 3D heat conduction problems to show the effectiveness of the proposed method.

The remainder of the paper is organized as follows. Section 2 reviews the related work in IGA.

Section 3 describes the variational harmonic method for analysis-suitable parameterization of 2D computational domain. The 3D case for volume parameterization is studied in Section 4. Some examples and comparisons based on the isogeometric heat conduction problem are presented in Section 5. Finally, we conclude this paper and outline future works in Section 6.

Related work

In this section, we will review some related works in IGA and parameterization of computational domains.

IGA was firstly proposed by Hughes et al. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement[END_REF] in 2005 to achieve the seamless integration of CAD and FEA. Since then, many researchers in the fields of computational mechanical and geometric computation were involved in this topic. The current work on isogeometric analysis can be classified into four categories: (1) application of IGA to various simulation and analysis problems [START_REF] Auricchio | A fully locking-free isogeometric approach for plane linear elasticity problems: A stream function formulation[END_REF] [START_REF] Xu | Analysis-suitable volume parameterization of multi-block computational domain in isogeometric applications[END_REF].

The topic of this paper belongs to the forth field. From the viewpoint of graphics applications, volume parameterization of 3D models has been studied in [START_REF] Li | Harmonic Volumetric Map-ping for Solid Modeling Applications[END_REF][START_REF] Xia | Parameterization of star shaped volumes using Green's functions[END_REF][START_REF] Xia | Direct product volume parameterization using harmonic fields[END_REF]. As far as we know, there are only a few work on the parametrization of computational domains from the viewpoint of isogeometric applications. Martin et al. [START_REF] Martin | Volumetric parameterization and trivariate B-spline fitting using harmonic functions[END_REF] proposed a method to fit a genus-0 triangular mesh by B-spline volume parameterization, based on discrete volumetric harmonic functions; this can be used to build computational domains for 3D IGA problems. A variational approach for constructing NURBS parameterization of swept volumes is proposed by Aigner et al [START_REF] Aigner | Swept volume parametrization for isogeometric analysis[END_REF]. Many free-form shapes in CAD systems, such as blades of turbines and propellers, are covered by this kind of volumes.

Cohen et.al. [START_REF] Cohen | Analysis-aware Modeling: Understanding Quality Considerations in Modeling for Isogeometric Analysis[END_REF] proposed the concept of analysis-aware modeling, in which the parameters of CAD models should be selected to facilitate isogeometric analysis. They also demonstrated the influence of parameterization of computational domains by several examples. Escobar et al. proposed a method to construct a trivariate T-spline volume of complex genus-zero solids for isogeometric application by using an adaptive tetrahedral meshing and mesh untangling technique [START_REF] Escobara | A new approach to solid modeling with trivariate T-spline based on mesh optimization[END_REF]. However, the obtained solid T-splines have elements with negative Jacobians at the Gauss quadrature points.

Zhang et al. proposed a robust and efficient algorithm to construct injective solid T-splines for

genus-zero geometry from a boundary triangulation [START_REF] Zhang | Solid T-spline construction from boundary representations for genus-zero geometry[END_REF]. Based on the Morse theory, a volumetric parameterization method of mesh model with arbitrary topology is proposed in [START_REF] Wang | Trivariate Solid T-spline construction from boundary triangulations with arbitrary genus topology[END_REF]. The above proposed methods demand a surface triangulation as input data. For the CAD boundary with spline representation, r-refinement method for generating optimal analysis-aware parameterization of computational domain is proposed based on shape optimization method [START_REF] Xu | Optimal analysis-aware parameterization of computational domain in 3D isogeometric analysis[END_REF] [START_REF] Xu | Parameterization of computational domain in isogeometric analysis: methods and comparison[END_REF]. However, it only works for specified analysis problems. In [START_REF] Xu | Analysis-suitable volume parameterization of multi-block computational domain in isogeometric applications[END_REF], a constraint optimization framework is proposed to obtain analysis-suitable volume parameterization of computational domain. Zhang et al. studied the construction of conformal solid T-spline from boundary T-spline representation by octree structure and boundary offset [START_REF] Zhang | Conformal solid T-spline construction from boundary T-spline representations[END_REF]. In this paper, based on the harmonic mapping theory, we propose a general method to construct analysis-suitable parameterization of computational domain from given CAD boundary information. 3. Variational harmonic method for planar parameterization of 2D computational domain

Problem statement

Consider a simply connected bounded domain S in two dimensional space with Cartesian coordinates (x; y) T . Suppose that S is bounded by four B-spline curves S(ξ, 0), S(ξ, 1), S(1, η) S(0, η). The parametric domain P of S should be a rectangular in two dimensional space with coordinates ξ, η, which is determined by the knot vector of boundary B-spline curves in ξ and η direction. The mapping from parametric space P to physical space S can be described as a B-spline surface S(ξ, η) = (x(ξ, η), y(ξ, η)) = n i=0 m j=0 N p i (ξ)N q j (η)P i,j with given four B-spline curves as boundaries. N p i (ξ) and N q j (η) are B-spline basis functions, P i,j are control points. Assume that this mapping is prescribed which maps the boundary of P one-to-one on the boundary of S. The parameterization problem of computational domain can be stated as: given four boundary B-spline curves, find the placement of inner control points such that the resulted planar B-spline surface is a good computational domain for isogeometric analysis.

Harmonic mapping and variational harmonic function

Harmonic mapping, which is a one-to-one transformation both for 2D and 3D regions, will be used in our parameterization method. Let σ : S → P be a harmonic mapping from S to P. From the theory of harmonic mapping, the inverse mapping σ -1 : P → S should be one-to-one (see Fig. 1 as an example). In this section, we will convert the harmonic conditions into some constraints the inner control points of planar B-spline parameterization should satisfy.

Suppose that σ : S → P is a harmonic mapping, that is ∆ξ(x, y) = ξ xx + ξ yy = 0 ∆η(x, y) = η xx + η yy = 0 By chain rules, we have

dξ = ξ x dx + ξ y dy = ξ x (x ξ dξ + x η dη) + ξ y (y ξ dξ + y η dη) dη = η x dx + η y dy = η x (x ξ dξ + x η dη) + η y (y ξ dξ + y η dη)
From above two equation, we have

ξ x x ξ + ξ y y ξ = η x x η + η y y η = 1
(1)

ξ x x η + ξ y y η = η x x ξ + η y y ξ = 0 (2)
By solving above two linear systems, we obtain

ξ x = y η J , ξ y = -x η J , η x = -y ξ J , η y = x ξ J where J = x ξ y η -x η y ξ
is Jacobian of transformation from P to S.

From above results and

∂ ∂x = ξ x ∂ ∂ξ , ∂ ∂y = ξ y ∂ ∂ξ , we have J∆ξ(x, y) = J ∂ ∂x ξ x + J ∂ ∂y ξ y = (Jξ x ∂ ∂ξ + Jη x ∂ ∂η )ξ x + (Jξ y ∂ ∂ξ + Jη y ∂ ∂η )ξ y = (y η ∂ ∂ξ -y ξ ∂ ∂η )( y ξ J ) + (-x η ∂ ∂ξ + x ξ ∂ ∂η )( -x η J ) = 0
After differential computation, we have

-y η Lx + x η Ly = 0 where L = (x 2 η + y 2 η ) ∂ 2 ∂ξ 2 -2(x ξ x η + y ξ y η ) ∂ 2 ∂ξη + (x 2 ξ + y 2 ξ ) ∂ 2 ∂η 2
Similarly, from J∆η(x, y) = 0, we have

y ξ Lx -x ξ Ly = 0
Then we obtain

Lx(ξ, η) = Ly(ξ, η) = 0 that is, LS(ξ, η) 2 = (Lx) 2 + (Ly) 2 = 0
This is a nonlinear system in terms of inner control points of the planar B-spline parameterization.

The final grid obtained by this method usually has non-uniform elements near convex(concave) boundary region. In order to solve this problem and achieve a grid with good orthogonality, we minimize the following energy function as in [START_REF] Xu | Parameterization of computational domain in isogeometric analysis: methods and comparison[END_REF] λ

1 ( S ξξ 2 + S ηη 2 +2 S ξη 2 ) +λ 2 ( S ξ 2 + S η 2 )dudv
Then by adding the harmonic mapping term into above energy function, we obtain a new energy function for parameterization of computational domain

LS(ξ, η) 2 +λ 1 ( S ξξ 2 + S ηη 2 +2 S ξη 2 ) + λ 2 ( S ξ 2 + S η 2 )dudv (3)
where λ 1 and λ 2 are positive weights to control the final parameterization results.

The Discrete Coons method is used to construct initial inner control points, and the Steepest Descent method is employed to minimize the energy function. The details will be described in the following two subsections.

Initial construction of inner control points

In order to solve this constraint optimization problem, an initial construction of inner control points is required. We rely on the discrete Coons method presented in [START_REF] Farin | Discrete coons patches[END_REF] to generate inner control points as initial value from boundary control points. See Fig. 2 (a) for an example.

Given the boundary control points P 0,j , P n,j , P i,0 , P i,m , i = 0, . . . , n, j = 0, . . . , m, the inner control points P i,j (i = 1, . . . , n -1, j = 1, . . . , m -1) can be constructed by the discrete Coons method as follows: Since the sum of the coefficients equals 1, the resulting inner control points lie in the convex hull of the boundary control points. For some given boundary curves, this construction may cause some selfintersections, and lead to an improper parameterization for IGA. See Fig. 2 (b) for such an example.

P i,j = (1 - i n )P 0,j + i n P n,j + (1 - j m )P i,0 + j m P i,m -[1 - i n i n ]   P 0,0 P 0,m P n,0 P n,m     1 -j m j m   (a) (b)

Optimization method

In the proposed approach, we minimize the objective function (3) , by moving inner control points of the computational domain. Therefore, we consider as optimization variables the coordinates of the inner control points and as cost function the error of the IGA solution. The optimization algorithm used for this study is a classical steepest-descent method in conjunction with a back-tracking line-search. For this exercise, the gradient of the cost function is approximated using a centered finite-differencing scheme.

Each iteration k of the optimization algorithm can be summarized as follows, starting from a point x k in the variable space:

1. Evaluation of perturbed points x k + ǫe k 2. Estimation of the gradient ∇f (x k ) by finite-difference

Define search direction d

k = -∇f (x k ) 4. Line search : find ρ such as f (x k + ρd k ) < f (x k )
These steps are carried out until a stopping criterion is satisfied.

Overview of variational harmonic method

In general, the proposed variational harmonic method can be summerized as follows:

Input: four coplanar boundary B-spline curves Output: inner control points and the corresponding planar B-spline surfaces • Construct the initial inner control points by discrete Coons method as shown in subsection 3.3;

• Solve the following optimization problem by using steepest-descent method

Min LS(ξ, η) 2 +λ 1 ( S ξξ 2 + S ηη 2 +2 S ξη 2 ) + λ 2 ( S ξ 2 + S η 2 )dudv
• Generate the corresponding planar B-spline surface S(ξ, η) as computational domain.

Variational harmonic method for volume parameterization of 3D computational domain

In this section, we will study the volume parameterization problem by variational harmonic method.

Suppose that S is a simply connected bounded domain in three dimensional space with Cartesian coordinates (x; y; z) T , and is bounded by six B-spline surfaces S(ξ, 0, ζ), S(ξ, 1, ζ), S(1, η, ζ), S(0, η, ζ), S(ξ, η, 0), S(ξ, η, 1).

The parametric domain P of S should be a cube in three dimensional space with coordinates ξ, η, ζ, which is determined by the knot vector of boundary B-spline surfaces in ξ, η and ζ direction. The mapping from parametric space P to physical space S can be described as a B-spline volume

S(ξ, η, ζ) = (x(ξ, η, ζ), y(ξ, η, ζ), z(ξ, η, ζ)) = n i=0 m j=0 l k=0 N p i (ξ) N q j (η)N r k (ζ)P i,j
,k with given six B-spline surfaces as boundaries. N p i (ξ), N q j (η) and N r k (ζ) are B-spline basis functions, P i,j,k are control points. Assume that this mapping is prescribed which maps the boundary of P one-to-one on the boundary of S. The volume parameterization problem of 3D computational domain can be stated as: given six boundary B-spline surfaces, find the placement of inner control points such that the resulted trivariate B-spline parametric volume is a good computational domain for 3D isogeometric analysis.

Harmonic mapping and variational harmonic function

Similar with the 2D case, the harmonic conditions in 3D case will be converted into some constraints the inner control points of trivariate B-spline volume parameterization should satisfy.

If the mapping σ : S → P is a harmonic mapping, we have

∆ξ(x, y, z) = ξ xx + ξ yy + ξ zz = 0 ∆η(x, y, z) = η xx + η yy + η zz = 0 ∆ζ(x, y, z) = ζ xx + ζ yy + ζ zz = 0
By chain rules, we have

dξ = ξ x dx + ξ y dy + ξ z dz = ξ x (x ξ dξ + x η dη + x ζ dζ) +ξ y (y ξ dξ + y η dη + y ζ dζ) + ξ z (z ξ dξ + z η dη + z ζ dζ) dη = η x dx + η y dy + η z dz = η x (x ξ dξ + x η dη + x ζ dζ) +η y (y ξ dξ + y η dη + y ζ dζ) + η z (z ξ dξ + z η dη + z ζ dζ) dζ = ζ x dx + ζ y dy + ζ z dz = ζ x (x ξ dξ + x η dη + x ζ dζ) +ζ y (y ξ dξ + y η dη + y ζ dζ) + ζ z (z ξ dξ + z η dη + z ζ dζ)
From above three equation, we have

ξ x x ξ + ξ y y ξ + ξ z z ξ = 1 η x x η + η y y η + η z z η = ζ x x ζ + ζ y y ζ + ζ z z ζ = 1 ξ x x η + ξ y y η + ξ z z η = η x x ξ + η y y ξ + η z z ξ = 0 ξ x x ζ + ξ y y ζ + ξ z z ζ = ζ x x ξ + ζ y y ξ + ζ z z ξ = 0 η x x ζ + η y y ζ + η z z ζ = ζ x x η + ζ y y η + ζ z z η = 0
By solving above three linear systems, we obtain

ξ x = y η z ζ -y ζ z η J , ξ y = - x η z ζ -x ζ z η J , ξ z = y η z ζ -y ζ z η J , η x = y ξ z ζ -y ζ z ξ J , η y = - x ξ z ζ -x ζ z ξ J , η z = y ξ z ζ -y ζ z ξ J , ζ x = y η z ξ -y ξ z η J , ζ y = - x η z ξ -x ξ z η J , ζ z = y η z ξ -y ξ z η J ,
where

J = x ξ y ξ z ξ x η y η z η x ζ y ζ z ζ
is Jacobian of transformation from P to S.

From above results and J∆ξ(x, y, z) = 0,J∆η(x, y, z) = 0, J∆ζ(x, y, z) = 0, the following condition can be derived,

Lx(ξ, η, ζ) = Ly(ξ, η, ζ) = Lz(ξ, η, ζ) = 0 (4)
where Eq. ( 4) can be rewritten as

L = a 11 ∂ 2 ∂ξ 2 + 2a 12 ∂ 2 ∂ξη + 2a 13 ∂ 2 ∂ξζ + a 22 ∂ 2 ∂η 2 + 2a 23 ∂ 2 ∂ηζ + a 33 ∂ 2 ∂ζ 2 ,
LS(ξ, η, ζ) 2 = (Lx) 2 + (Ly) 2 + (Lz) 2 = 0
Similarly with the 2D problem, it is also a nonlinear system in terms of inner control points of the Bspline volume parameterization. In order to achieve uniform and orthogonal iso-parametric structure near convex(concave) boundary region, we minimize the following energy function as in [START_REF] Xu | Analysis-suitable volume parameterization of multi-block computational domain in isogeometric applications[END_REF] λ

1 ( S ξξ 2 + S ηη 2 +2 S ξη 2 +2 S ηζ 2 + 2 S ξζ 2 + S ζζ 2 ) + λ 2 ( S ξ 2 + S η 2 + S ζ 2 )dudv
A new energy function for volume parameterization of computational domain can be obtained by combining the harmonic mapping term and the above energy function,

LS(ξ, η, ζ) 2 +λ 1 ( S ξξ 2 + S ηη 2 +2 S ξη 2 +2 S ηζ 2 +2 S ξζ 2 + S ζζ 2 ) +λ 2 ( S ξ 2 + S η 2 + S ζ 2 )dudv (5)
where λ 1 and λ 2 are positive weights to control the final parameterization results.

Initial construction of inner control points

The initial construction of inner control points for volume parameterization is based on the discrete Coons volume generation method presented in [START_REF] Farin | Discrete coons patches[END_REF] from given boundary control points.

Suppose that the given boundary surfaces are B-spline surfaces and that the opposite boundary B-spline surfaces have the same degree, number of control points and knot vectors. The interior control points P i,j,k can be constructed as linear combinations of points described by the following formulas:

P i,j,k = (1 -i/l)P 0,j,k + i/lP l,j,k + (1 -j/m)P i,0,k +j/mP i,m,k + (1 -k/n)P i,j,0 + k/nP i,j,n -[1 -i/l, i/l]   P 0,0,k P 0,m,k P l,0,k P l,m,k     1 -j/m j/m   -[1 -j/m, j/m]   P i,0,0 P i,0,n P i,m,0 P i,m,n     1 -k/n k/n   -[1 -k/n, k/n]
  P 0,j,0 P l,j,0 P 0,j,n P l,j,n N p i (ξ)N q j (η)N r k (ζ)P i,j,k .

    1 -i/l i/l   +(1 -k/n)   [1 -i/l, i/l]
where N p i (ξ), N q j (η) and N r k (ζ) are B-spline basis function with knot vectors given by boundary surfaces.

Overview

The proposed variational harmonic method for volume parameterization can be summerized as follows:

Input: six boundary B-spline surfaces Output: inner control points and the corresponding B-spline volumes

• Construct the initial inner control points by discrete Coons method as shown in subsection 4.2;

• Solve the following optimization problem by using steepest-descent method presented in ;

Min LS(ξ, η, ζ) 2 +λ 1 ( S ξξ 2 + S ηη 2 +2 S ξη 2 +2 S ηζ 2 +2 S ξζ 2 + S ζζ 2 ) +λ 2 ( S ξ 2 + S η 2 + S ζ 2 )dudv
• Generate the corresponding B-spline volume S(ξ, η, ζ) as computational domain.

Examples and comparison

In this section, we aim at presenting several examples to show the efficiency of the proposed method.

We will also give a comparison study between the original parameterization and the final parameterization constructed by the proposed method based on the 2D and 3D heat conduction problem.

Test model -heat conduction problem

For ease of presentation, we consider the second order elliptic PDE with homogeneous Dirichlet boundary condition as an illustrative model problem :

-∆U (x) = f (x) in Ω U (x) = 0 on ∂Ω (6)
where x are the Cartesian coordinates, Ω is a Lipschitz domain with boundary ∂Ω, f (x) ∈ L 2 (Ω) : Ω → R is a given source term, and U (x) : Ω → R is the unknown solution.

Starting from a planar B-spline surface (or trivariate B-spline parametric volume) as computational domain, a general framework of an isogeometric solver for 2D and 3D heat conduction problem (6) has been implemented as a plugin in the AXEL 1 platform, yielding a B-spline surface (or B-spline volume) as solution field. Additional details concerning the isogeometric solver of problem ( 6) can be found in [START_REF] Duvigneau | An introduction to isogeometric analysis with application to thermal conduction[END_REF]. The proposed variational harmonic method is implemented as a part of the isogeometric toolbox of the project EXCITING 2 .

In this paper, we test the different parameterizations of computational domains for 2D heat conduction problem [START_REF] Bazilevs | Isogeometric analysis using T-Splines[END_REF] with source term

f (x, y) = 4π 2 9 sin( πx 3 ) sin( πy 3
).

and 3D problem with source term

f (x, y, z) = π 2 3 sin( πx 3 ) sin( πy 3 ) sin( πz 3 ). (8)
For problems with unknown exact solution U , suppose that U h is the approximation solution obtained by isogeometric method, then the discrete error e = U -U h . We employ a posteriori error assessment proposed in [START_REF] Xu | Parameterization of computational domain in isogeometric analysis: methods and comparison[END_REF] to compare the initial and final parameterization of the computational domain. Example # Deg. # Con. #Iter.

Conclusion

Parameterization of computational domain is a key point in isogeometric analysis. In this paper, we propose a variational harmonic method to parameterize the 2D and 3D computational domain from the theory of harmonic mapping. The resulted parameterization is injective and has high quality near convex(concave) parts of the boundary. Examples and comparison are presented to show that the proposed methods can produce analysis-suitable parameterization of computational domain for isogeometric analysis.

As part of the future work, we will generalize the proposed method to the computational domain with complex shape, which are more important in practice.

Fig. 1 .

 1 Fig. 1. The mapping σ from physical domain S to parametric domain P.

Fig. 2 .

 2 Fig. 2. Examples of B-spline surface constructed by discrete Coons method: (a) example without self-intersection ; (b) example with self-intersections. Green control points are inner control points to be constructed.

and a 11 =

 11 (S ξ , S ξ ), a 12 = (S ξ , S η ), a 13 = (S ξ , S ζ ), a 22 = (S η , S η ), a 23 = (S η , S ζ ), a 33 = (S ζ , S ζ ).

  0,0 P 0,m,0 P l,0,0 P l,m,0 -i/l, i/l]   P 0,0,n P 0,m,n P l,0,n P l,m,n

Fig. 3 .Fig. 4 .

 34 Fig. 3. 2D example I.

Fig. 5 .

 5 Fig. 5. 3D example I.

Fig. 6 .

 6 Fig. 6. 3D example II.

Table 1 :

 1 Quantitative data in Fig.3, Fig.4, Fig.5and Fig.6. # deg.: degree of B-spline parameterization; # con.: number of control points; # iter.: number of optimization iterations.

  Fig.3 3 36 4 Fig.4 4 64 7 Fig.5 3 125 7 Fig.6 4 216 11

  a 11 = a 22 a 33 -a 2 23 , a 12 = a 13 a 23 -a 12 a 33 , a 13 = a 12 a 23 -a 13 a 22 , a 22 = a 11 a 33 -a 2 13 , a 23 = a 13 a 12 -a 11 a 23 , a 33 = a 11 a 22 -a 2 12 ,

Acknowledgment

The authors are supported by the 7th Framework Program of the European Union, project SCP8-218536 "EXCITING". The first author is partially supported by the National Nature Science Foundation of China (Nos. 61004117, 61272390, 61211130103,61003193), the Defense Industrial Technology Development Program(A3920110002), the ScientificResearch Foundationforthe Returned Overseas Chinese Scholars from State Education Ministry, and the Scientific Research Starting Foundation of Hangzhou Dianzi University (No. KYS055611029).

The approximation error e from (9) also has a B-spline form. Some h-refinement operation should be performed to achieve more accurate results for above problem. Though it is much more expensive, it can be used as an error assessment method to show the effectiveness of the proposed construction method of computational domain.

The first 2D example is shown in Fig. 3. Fig. 3 (a) presents the given boundary B-spline curves. Another 2D example is illustrated in Fig. 4. The given boundary B-spline curves is shown in Fig. 4 (a). Fig. 4(b) presents the initial parametrization of computational domain constructed by discrete Coons method. There are some self-intersections on the initial parameterization as shown by the isoparametric curves in Fig. 4(d). Fig. 4(c) shows the optimized parametrization of computational domain constructed by variational harmonic method.The optimized parameterization has no self-intersection as illustrated in Fig.

4(e).

Fig. 5 shows a 3D example, which is drawn partly to illustrate the interior information of the volume.

The given boundary B-spline surfaces and curves are shown in Fig. 5 (a) and 5 (b). Fig. 5 (c) presents the initial volume parametrization of computational domain constructed by discrete Coons method. There are some self-intersections on the initial parameterization. Fig. 5 (d) shows the final volume parameterization of computational domain without self-intersections constructed by variational harmonic method. To illustrate the quality of the parameterization, the iso-parametric surfaces of the resulted B-spline volume are presented in Fig. 5 (e) and 5 (f). Fig. 5 (g) and 5 (h) presents the corresponding solution field obtained from the volume parameterization in Fig. 5 (c) and Fig. 5 (d) . Obviously, the solution field in Fig. 5 (h) is more smooth than the solution field in Fig. 5 (g), which indicates that the final volume parameterization is better than the initial one. Another 3D example with complex geometry is shown in Fig. 6.

Quantitative data of the examples presented in Fig. 3, Fig. 4, Fig. 5 and Fig. 6 are summarized in Table .1. Overall, the optimized parameterization obtained by the variational harmonic method has better quality, and can achieve better simulation results than the initial Coons parameterization.


H Louche
email: h.louche@univ-montp2.fr

P Schlosser

D Favier
email: denis.favier@grenoble-inp.fr

L Orgéas

Heat source processing for localized deformation with non-constant thermal conductivity. Application to superelastic tensile tests of NiTi Shape Memory Alloys

Keywords: Kinematic and thermal field measurements, Infrared thermography and heat source processing, Localized phase transformation, NiTi shape memory alloy, In situ DSC

This paper presents a methodology for heat source estimation when thermal conductivity is not constant. Both thermal and kinematic field measurements are necessary. The method requires spatiotemporal synchronization of these fields and a possible method to achieve this is presented. Temperature and heat source fields can then be properly estimated at every material points in the reference or deformed configurations of the tested sample.

The proposed method also highlights the importance of precise determination of thermophysical properties in heat estimations.

The method is applied for heat source estimation during a superelastic tensile test of a NiTi shape memory alloy displaying both stress-induced phase transformation and deformation localization. Spatiotemporal change in thermal conductivity occurs during deformation of this material. Typical heat source distribution results observed during mechanical tests are given. Heat sources are integrated over time, at each material point, for local estimation of heat energy associated with stress-induced phase transformation. Such a measurement can be

qualified as a first step towards an in situ local DSC during a mechanical test. It is shown that the spatiotemporal distributions of these energy fields correlate well with strain field distributions during superelastic localized deformation.

INTRODUCTION

Thermal field measurements are valuable tools for studying thermomechanical couplings which are often associated with deformation of materials. The intensity and spatiotemporal distribution of temperature variations depend on the material behavior (elasticity, plasticity, viscosity, damage, phase transformation, etc.) but also on the different characteristics of the thermal problem (geometry, thermal properties, boundary conditions, etc.).

Previous work [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localization[END_REF] highlighted the interest of studying local heat sources that cause temperature variations, e.g. when studying Lüders bands and necking in steels [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localization[END_REF][START_REF] Louche | Thermal and dissipative effects accompanying Lüders band propagation[END_REF], Portevin Le Châtelier bands in AlMg alloys [START_REF] Louche | Thermal observations associated with the Portevin-le Chatelier effect in an Al-Mg alloy[END_REF], plasticity in Al oligocrystals [START_REF] Saai | Experimental and numerical study of the thermomechanical behavior of al bi-crystal in tension using full field measurements and micromechanical modeling[END_REF] and stress-induced phase transformation in Nickel-Titanium (NiTi) shape memory alloys [START_REF] Schlosser | Image processing to estimate the heat sources related to phase transformations during tensile tests of NiTi tubes[END_REF][START_REF] Favier | Homogeneous and heterogeneous deformation mechanisms in an austenitic polycrystalline Ti-50.8at.%Ni thin tube under tension. Investigation via temperature and strain fields measurements[END_REF]. In all of these studies, the thermal conductivity of the tested sample was assumed to be constant during deformation. However, this hypothesis no longer applies for some materials and/or types of loading.

The first objective of this paper is to present a methodology to obtain better heat source estimations when thermal conductivity is not constant. This method requires kinematic field measurements and spatiotemporal synchronization of thermal and kinematic full fields measured during a mechanical test.

The second objective is to apply the above methodology to the analysis of superelastic tensile tests in the case of nickel-titanium (NiTi) shape memory alloys (SMA). When superelastic deformation involves reversible stress-induced phase transformation, from austenitic phase (A) to martensitic phase (M), non-constant thermal conductivity occurs. NiTi SMAs have been extensively used in many engineering applications, especially in biomedical products (including stents and catheters) where the superelasticity of these alloys is required [START_REF] Favier | Influence of manufacturing methods on the homogeneity and properties of Nitinol tubular stents[END_REF]. In spite of the ever-increasing interest in SMA modeling, there are still many unsolved problems with the predictability of the thermomechanical responses of this material [START_REF] Orgéas | Stress state effect on mechanical behaviour of shape memory alloys: experimental characterisation and modeling[END_REF]. One of the main limitations is the deformation localization that appears systematically during the stress plateau in tensile tests [START_REF] Sittner | On the origin of Lüders-like deformation of Ni-Ti shape memory alloys[END_REF][START_REF] Brinson | Stress-induced transformation behavior of a polycrystalline NiTi shape memory alloy: micro and macromechanical investigations via in situ optical microscopy[END_REF][START_REF] Daly | Stress-induced martensitic phase transformation in thin sheets of nitinol[END_REF][START_REF] Pieczyska | Investigation of nucleation and propagation of phase transitions in NiTi SMA[END_REF][START_REF] Miyazaki | Lüders-like deformation observed in the transformation pseudoelasticity of a TiNi alloy[END_REF][START_REF] Shaw | Thermomechanical aspects of NiTi[END_REF][START_REF] Shaw | On the nucleation and propagation of phase transformation fronts in a NiTi alloy[END_REF][START_REF] Orgéas | Stress-induced martensitic transformation of a NiTi alloy in isothermal shear, tension and compression[END_REF]. This localization may disappear for other mechanical loadings such as compression, shear or equibiaxial tension [START_REF] Orgéas | Stress-induced martensitic transformation of a NiTi alloy in isothermal shear, tension and compression[END_REF][START_REF] Manach | Shear and Tensile Thermomechanical Behaviour of Near Equiatomic NiTi Alloy[END_REF][START_REF] Grolleau | Assessment of tensioncompression asymmetry of NiTi using circular bulge testing of thin plates[END_REF] or tensiontorsion combined tests [START_REF] Sun | Phase transformation in superelastic NiTi polycrystalline micro-tubes under tension and torsion from localisation to homogeneous deformation[END_REF]. It may also be observed during the ferroelastic tensile deformation of polycrystalline NiTi SMAs [START_REF] Liu | Lüders-like deformation associated with martensite reorientation in niti[END_REF]. This localization has been studied using qualitative optical observations [START_REF] Miyazaki | Lüders-like deformation observed in the transformation pseudoelasticity of a TiNi alloy[END_REF][START_REF] Shaw | On the nucleation and propagation of phase transformation fronts in a NiTi alloy[END_REF], multiple extensometers [START_REF] Shaw | Thermomechanical aspects of NiTi[END_REF][START_REF] Orgéas | Stress-induced martensitic transformation of a NiTi alloy in isothermal shear, tension and compression[END_REF], strain [START_REF] Favier | Homogeneous and heterogeneous deformation mechanisms in an austenitic polycrystalline Ti-50.8at.%Ni thin tube under tension. Investigation via temperature and strain fields measurements[END_REF][START_REF] Daly | Stress-induced martensitic phase transformation in thin sheets of nitinol[END_REF][START_REF] Grolleau | Assessment of tensioncompression asymmetry of NiTi using circular bulge testing of thin plates[END_REF], temperature [START_REF] Favier | Homogeneous and heterogeneous deformation mechanisms in an austenitic polycrystalline Ti-50.8at.%Ni thin tube under tension. Investigation via temperature and strain fields measurements[END_REF][START_REF] Pieczyska | Investigation of nucleation and propagation of phase transitions in NiTi SMA[END_REF][START_REF] Shaw | On the nucleation and propagation of phase transformation fronts in a NiTi alloy[END_REF] full-field measurements and recently by in situ EBSD investigations [START_REF] Mao | EBSD studies of the stress-induced B2-B19' martensitic transformation in NiTi tubes under uniaxial tension and compression[END_REF]. Temperature measurements obtained by infrared thermography was previously used in [START_REF] Favier | Homogeneous and heterogeneous deformation mechanisms in an austenitic polycrystalline Ti-50.8at.%Ni thin tube under tension. Investigation via temperature and strain fields measurements[END_REF] to estimate heat sources during the homogeneous stages of a tensile test on NiTi tubes. Such heat sources, linked to the local latent heat of the phase transformations involved in SMAs, brought new and precious infor-mation on deformation mechanisms in SMAs. When phase transformation is localized, the method to estimate localized heat sources is different; see [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localization[END_REF] for thin sheets and [START_REF] Schlosser | Image processing to estimate the heat sources related to phase transformations during tensile tests of NiTi tubes[END_REF] for thin tubes.

Thermophysical parameters, such as the specific heat C and the thermal conductivity k of the material, are necessary to estimate the heat sources. Few published studies have been attempted to monitor these thermophysical parameters in NiTi materials. For example, C and k were assumed to be constant in our previous studies carried out with polycrystaline NiTi SMAs [START_REF] Schlosser | Image processing to estimate the heat sources related to phase transformations during tensile tests of NiTi tubes[END_REF][START_REF] Favier | Homogeneous and heterogeneous deformation mechanisms in an austenitic polycrystalline Ti-50.8at.%Ni thin tube under tension. Investigation via temperature and strain fields measurements[END_REF]. However, some experimental studies (e.g. see [START_REF] Faulkner | Experimental determination of thermal and electrical properties of Ni-Ti shape memory wires[END_REF]) have shown substantial differences in C and k for austenite and martensite. The differences seem to be greater for the thermal conductivity. By working with a polycrystalline NiTi wire, Faulkner et al. [START_REF] Faulkner | Experimental determination of thermal and electrical properties of Ni-Ti shape memory wires[END_REF] obtained a ratio of two between the thermal conductivity of austenitic k A and martensite phases k M , i.e. k A =28 W m -1 K -1 and k M =14 W m -1 K -1 . Specific heat values given in the literature are also scattered, but in our study we decided to take a mean specific heat for the two phases for the sake of simplicity: C=450 J kg -1 K -1 . However, different k values (through different models) are proposed in order to take the possible variations in k with the volume fraction of martensite into account. Information on phase fraction changes during loading is thus obtained from local strain levels, measured via digital image correlation (DIC) applied to visible images recorded during the test.

Simultaneous strain field (from DIC applied to visible images) and thermal field measurements (obtained from infrared thermography) during mechanical testing were first reported in [START_REF] Wattrisse | Dissipative mechanisms during the necking of a semicrystalline polymer[END_REF] and [START_REF] Louche | Experimental energetic balance associated to the deformation of an aluminum multicrystal and monocrystal sheet[END_REF]. In the first work, kinematic fields were used to locally compute the true stress and to estimate the convective term in the heat diffusion equation. In the second paper, temperature, heat source and axial strain fields were qualitatively compared during a tensile test on an Al multicrystal. More recently, three other studies [START_REF] Chrysochoos | Fields of stored energy associated with localized necking of steel[END_REF][START_REF] Bodelot | Experimental setup for fully coupled kinematic and thermal measurements at the microstructure scale of an aisi 316L steel[END_REF][START_REF] Pottier | Study on the use of motion compensation techniques to determine heat sources. Application to large deformations on cracked rubber specimens[END_REF] proposed a major improvement in the synchronization of kinematics and thermal fields. They offer the possibility of computing the temperature field in the reference configuration, i.e. to monitor the temperature of material points, even in the case of in-plane displacements of these points. In [START_REF] Chrysochoos | Fields of stored energy associated with localized necking of steel[END_REF], the two cameras each monitored one side of a thin plate steel sample. Heat sources were estimated from the thermal fields in the reference configuration. They were then integrated over time in order to achieve energy balances. An interesting optical set-up was proposed in [START_REF] Bodelot | Experimental setup for fully coupled kinematic and thermal measurements at the microstructure scale of an aisi 316L steel[END_REF] involving the use of a dichroic mirror, thus allowing two-cameras monitoring of the same surface of the sample. In these two latter research works, the temperature fields were projected in the reference configuration through displacement fields measured by DIC. A last study [START_REF] Pottier | Study on the use of motion compensation techniques to determine heat sources. Application to large deformations on cracked rubber specimens[END_REF] proposed a motion compensation technique based on detection of reflective spots in visible and infrared images. In the present study, the same surface of the sample was simultaneously observed by a visible and an infrared camera. Thanks to a compatible coating, it was then possi-ble to compute kinematics from DIC and temperature fields from radiations coming from the surface during a tensile test.

After a description of the experimental setup (section 2), the adopted spatiotemporal synchronization of the thermal and displacement fields is presented (section 3). This processing allows us time-course variations in temperature for each material point of the sheet sample to be determined. A method to estimate heat sources in the case of different thermal conductivities between the austenite and martensite phases is then presented (section 4). Finally, typical results showing heat source and energy distributions during mechanical tests are given and compared with strain fields (section 5).

EXPERIMENTAL SET-UP

Heat source fields are estimated during the deformation of a specimen. Different types of loading can be studied. In many works, as here, quasi-static tensile tests are conducted on sheet samples at room temperature. In this study, a fast multidetector infrared camera (CEDIP Jade III MW, 14 bits, 320 x 240 pixels, thermal resolution 20 mK) was used to measure the temperature field T(M, t) in spatial points M of spatial position x  at current time t. The infrared frame rate was f IR =5 Hz. The spatial resolution (pixel size) of the thermal fields was estimated at 0.22 mm for the conducted test. The displacement field ! u(M, t) at the sample surface was obtained using a digital visible camera (Hamamatsu, 12 bits, 1280 x 1024 pixels) and DIC processing software (7D software [START_REF] Vacher | Bidimensional deformation measurement using digital images[END_REF]). The frame rate for the visible images was f v = 0.5 Hz. The grid size was 10 x 10 pixels and the spatial resolution (pixel size) was estimated to be 0.065 mm. The resolution of the kinematic fields (grid step) was then 0.65 mm (see [START_REF] Louche | Kinematic fields and acoustic emission observations associated with the Portevin Le Châtelier effect on an AlMg alloy[END_REF] and Figure 1b). The Green-Lagrange strain fields were estimated from the displacement field ! u and can be represented in the referential (Lagrangian description) or deformed (Eulerian description) configu- rations.

The sample surface was coated (spray projection) with high emissivity black paint (commercial black paint spray) with a white texture on top of it, as shown in Figure 1a. As shown in [START_REF] Schlosser | Infuence of thermal and mechanical aspects on deformation behaviour of NiTi alloys[END_REF], this coating is compatible with both the infrared and visible observations and its measured emissivity ε was 0.98, i.e. very close to black body emissivity (ε = 1). This difference was neglected and ε =1 was used in the following. The optical axis of the visible camera was normal to this surface and the infrared camera was slightly oriented (at 15°) in the normal direction. Because the thermal radiation of a black body is not influenced by disorientation of the camera below 45°, our setup allows quantitative measurement of thermal fields [START_REF] Schlosser | Infuence of thermal and mechanical aspects on deformation behaviour of NiTi alloys[END_REF].

As in [START_REF] Favier | Homogeneous and heterogeneous deformation mechanisms in an austenitic polycrystalline Ti-50.8at.%Ni thin tube under tension. Investigation via temperature and strain fields measurements[END_REF], a mean nominal axial strain ! = !l l 0 was measured from the displacement fields using the relative displacement Δl of two material points A and B located at the extremities of the gauge length (Figures 1b,1c), with l 0 =58 mm being the initial distance between these two points.

SPATIOTEMPORAL SYNCHRONIZATION OF KINEMATICS AND THERMAL FIELDS

The method proposed in the present paper is the same type as that proposed in [START_REF] Chrysochoos | Fields of stored energy associated with localized necking of steel[END_REF][START_REF] Bodelot | Experimental setup for fully coupled kinematic and thermal measurements at the microstructure scale of an aisi 316L steel[END_REF]. The main characteristics of the methods used in these two references are recalled at the end of the Introduction section. However, our approach differed slightly from these studies since the two cameras directly monitored the same sample surface.

First, the temporal synchronization necessary to get the signals at the same time is presented. Once this temporal synchronization is achieved, spatial synchronization is also done to obtain the temperature field at each material point. This synchronization will be presented in the second subsection.

Temporal synchronization

During the experiments, mechanical data (including force, crosshead displacement), visible images together with infrared images were recorded. The three acquisitions have their own time base, defined by the initial time and frequency, as plotted schematically in Figure 2. It is then often difficult to achieve, with good accuracy, temporal correspondence between the kinematic field, the temperature field and the usual mechanical data (including force). Moreover, such a situation is enhanced when the acquisition frequencies differ from each other and/or transient phenomena are observed.

The method used in this work manages, with a homemade software, the data recorded from the three different time bases. The mechanical data and times of the visible images (i.e. of the infrared ones, respectively) were projected onto the infrared time base (i.e. onto the visible time base, respectively) using linear interpolations between the temporal data.

Spatial synchronization

The kinematic fields (displacements, strains) were obtained on the reference configuration (Lagrangian or material description) or in the current configuration (Eulerian or spatial description). During such measurements, the fixed infrared camera observes the temperatures of spatial positions whereas the material points can move towards the fixed observation area (see Figure 5). The temperature fields are then given as function of spatial coordinates: T(x,y,t), where x and y are the spatial coordinates of the material point P (with P having integer pixel values) at time t.

Therefore, three coordinate systems CS (or virtual grids) are available:

-CS-V: Visible image coordinate system, typically with a resolution of 1280 x 1024 pixels. In this coordinate system the initial and current positions of the material point are noted (X V ,Y V ) and (x V ,y V ), respectively.

-CS-DIC: DIC coordinate system given by the DIC software. Using a virtual grid on the visible image, every n = 10 pixels, the resolution in both directions is divided by a factor n. The maximum resolution is then about 128 x 102 pixels. In this coordinate system the initial and current positions of the material point are noted (X DIC ,Y DIC ) and (x DIC ,y DIC ), respectively.

-CS-IR: The infrared image coordinate system, with a maximum resolution of 320 x 240 pixels. In this coordinate system, the initial and current positions of the material point are noted (X IR ,Y IR ) and (x IR ,y IR ), respectively.

In the following, the subscripts 'V, DIC or IR' are omitted when the CS are indicated. The spatial synchronization can be achieved in any of the three coordinate systems. Due to the lower resolution of the kinematic grid compared to the infrared one, CS-DIC was chosen in this work. In that case, it enables the operator to know the temperature variation T(P,t) of a material point P, which is not possible with only one infrared camera. This is the same as having small thermocouples glued on the surface at the same position as the virtual DIC extensometers.

Spatial synchronization is useful if the sample is deformed (especially in case of large or localized deformations) or if it moves toward the camera (e.g. rigid body motion at the beginning of a mechanical test). Figure 3 shows the positions of selected material points A, B and P on a tensile sample prior to and after deformation (8% axial strain). The material point A was close to the lower grip. As this grip was fixed during the test, the vertical displacement of A was small. The point B was close to the moving grip. Its typical vertical displacement could go up to (u x ) maxi = 23 pixels (in CS-IR), which was no longer negligible during 8% tensile strain. Moreover, at the beginning of loading, a rigid body motion of about u x = 4 pixels and u y = -3 pixels occurred.

With the proposed spatial synchronization, it was possible:

(i) to calculate temperature variations of material points and not of spatial points. Without spatial synchronization, the temperature variation (noted θ 0 ) was estimated at each spatial point (pixel x IR ,y IR of the IR camera) in CS-IR using:

θ 0 (x IR ,y IR ,t)=T(x IR ,y IR ,t)-T 0 (x IR ,y IR ). ( 1
)
Initial measured temperature field on the sample may be heterogenous (T 0 (x IR ,y IR )≠ constant) either if a temperature gradient (often less than 0.5 K) is observed on the sample (due to the boundary conditions and the gripping system) or if the emissivity of the sample surface is not constant. So, if during the measurement there is displacement of the material points (rigid body motion or deformation), estimating the temperature variation using θ 0 can induce substantial errors. Instead, the temperature variation (noted θ) should be estimated for a material point P through:

θ(P,t)=Τ(P,t)-T(P,0)=T(x IR ,y IR ,t)-T 0 (X IR ,Y IR ), ( 2
)
where X IR and Y IR are the initial coordinates of the material point P.

(ii) to obtain heat sources at a material point and the corresponding energies (by integrating over time heat sources) in each material point (more details are given in section 5).

Spatial synchronization method

In order to allow spatial synchronization, both cameras need to be fixed during the measurement. The two coordinate systems CS-V and CS-IR, initially independent, were spatially linked by choosing reference points at the same initial time t 0 = 0, on the visible and infrared images.

The initial position (X V ,Y V ) on the visible image of any material point P can then be linked to its position (X IR ,Y IR ) on the infrared image using a rotation of angle α, a translation vector ! T = (T x ,T y ) , and a homothetic transformation κ to change the length scale.

A transformation matrix v2ir is built to go from the visible to the infrared description with:

X IR Y IR 1 ! " # # # $ % & & & = ! cos(") sin(") 0 'sin(") cos(") 0 ! " # $ % & X V Y V 1 ! " # # # $ % & & & + T x T y 1 ! " # # # $ % & & & = v2ir X V Y V 1 ! " # # # $ % & & & , (3)
and:

v2ir = ! cos(") ! sin(") T x !! sin(") ! cos(") T y 0 0 1 " # $ $ $ % & ' ' ' . (4)
Conversely, the matrix ir2v enables one to obtain the position in the visible system from the infrared one: ir2v = v2ir !1 . Note that the new coordinates after operation (3) are generally not integer values, but fractions of pixels. The four parameters of the transformation (α, κ,T x ,T y ) are identified from the initial positions of characteristic points selected in the initial visible and infrared images. A minimum of two characteristic points is necessary. However, the positions of such points are only known, in the two images, with one pixel accuracy. To reduce this uncertainty, the user can choose as many characteristic points as possible. Average values of the unknown parameters are then obtained from the set of selected points. An example is presented in Figure 4, where three characteristic points were easily chosen in the two initial images (CS-V and CS-IR) via a steel ruler stuck on the sample surface. In order to obtain the temperature of a material point P, at current time t, the spatial coordinates (x IR ,y IR ,t) of this point P in the infrared coordinate system must be calculated. This is achieved through the following steps (Figure 5):

1. The position at time t 0 of the material point P, located in CS-IR (Figure 5b), is computed in CS-DIC (Figure 5a) using ir2v transformation.

2. At time t, the position of the material point P(t) is given by the DIC software, in CS-DIC (Figure 5c).

3. The current position in CS-IR of P(t) is then obtained using v2ir transformation (Figure 5d).

4. Finally, the temperature of P(t) is estimated with T(x IR ,y IR ,t) in CS-IR.

This procedure can be used: (i) to express, in the initial CS-IR, the temperature of a material point over the time, or (ii) to project the temperature field in the initial CS-DIC. Spatial interpolation, bilinear or bicubic (not detailed here) are necessary to obtain the temperature in sub-pixel positions. Figure 6 shows such a situation in a common case where the resolution in CS-DIC is lower than in CS-IR.

Example of synchronization

One example of spatiotemporal synchronization was chosen to illustrate the differences between the two expressions of the temperature variation (θ and θ 0 ) and also to show the capacity of the method to monitor the temperature of a material point. This example was built from a real quasi-static tensile test conducted on a NiTi tube in quasi-isothermal testing conditions [START_REF] Schlosser | Infuence of thermal and mechanical aspects on deformation behaviour of NiTi alloys[END_REF] (see Figures 7a,7b).

During this test, a deformation localization propagated from the left to the right side of the sample, associated with a phase transformation front. The instantaneous strain field E xx , plotted in Figure 7c in Lagrangian description, shows two regions of very different strain values, with a boundary around material point B. The temperature variation fields (θ and θ 0 ) are plotted in Figures 7d and7e at the same time in CS-DIC, but not with the same description: Lagrangian for θ (i.e. on the initial configuration) and Eulerian for θ 0 (i.e. on the current configuration at time t). The thin thermal band, associated with the exothermic phase transformation, is well located at the material point B (Figure 7d), at the same position as the strain jump (Figure 7c). With the field θ 0 (x,y,t), obtained from equation ( 1) and presented in the Eulerian description (Figure 7e), the position and intensity of the thermal band are modified. The intensity change is enhanced in this example because of the non-uniformity of the initial thermal image T 0 and the radial rigid displacement observed at the beginning of the loading.

The time course of strain (Figure 7f) and temperature variations (Figures 7g and7h) at points A and B show the interest of monitoring the temperature of material points. The temperature variation peaks (Figure 7g), associated with deformation localization phenomena, coincide with the strain jumps in the case of the Lagrangian description θ(X,Y,t). In contrast, a time lag (Figure 7h) is observed in the Eulerian plot θ 0 (x,y,t), where x and y are the coordinates of the pixel coinciding with points A and B at the initial time. As noted before, a slight difference in amplitude, induced by the approximated evaluation (equation ( 1)) of the temperature variation, can also be observed in this figure.

HEAT SOURCE FIELD ESTIMATIONS

Heat diffusion equation

The method to estimate heat source fields on the basis of thermal fields was extensively explained in [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localization[END_REF]. At the current time t and at a spatial point of a domain Ω, located by its spatial position vector ! x , the heat diffusion equation model is:

! C( !T !t + ! v.grad " ! """ T) + div(-kgrad " ! """ T ) = s e + s i , " ! x,t # $ %[0,t f ] ( 5
)
where T is the absolute temperature, ! v the velocity vector, ρ the mass density, C the specific heat, k the heat conduction tensor, s e the heat source of external origin (e.g. radiation), s i the heat source of internal origin (e.g. thermoelastic effects, intrinsic dissipation, latent heat rate). The last term is dominant in the case of stressinduced transformation, as illustrated in Section 5.

In the previous studies [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localization[END_REF][START_REF] Schlosser | Image processing to estimate the heat sources related to phase transformations during tensile tests of NiTi tubes[END_REF], the following hypotheses were put forward: (i) the convective term ! v.!grad " ! """ T was considered as a negligible quantity, (ii) C was kept constant during deformation, (iii) k was considered as isotropic and constant, equal to k. Moreover, two other hypotheses were necessary to compute the heat sources:

the sample was thin so that (iv) s i and (v) T were assumed to be homogeneous throughout its thickness.

Taking the previous hypotheses into account and integrating equation ( 5) in the thickness direction leads to the following bi-dimensional heat diffusion equation:

! C( !" !t + " # )-k" 2 " = s i , # ! x,t $ % 2 D &[0,t f ] (6)
where T 0 is the initial thermal equilibrium state, θ=T-T 0 is the temperature variation such that θ≪ T 0 , ! 2 !

is the Laplacian calculated in the 2D mid-plane Ω 2D of the thin sample (sheets [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localization[END_REF] or tube [START_REF] Schlosser | Image processing to estimate the heat sources related to phase transformations during tensile tests of NiTi tubes[END_REF]), and where τ is a time constant modeling heat losses from the back and front surfaces of the thin sample. Finally note that assumption (i) can be validated depending on the testing material and conditions. In the example shown in Section 5, we checked that it was relevant (see below).

Heat source estimations from thermal field measurements

Taking the room temperature T 0 into account, the infrared camera gives the temperature variation fields θ at each spatial point (pixel) of the infrared image CS-IR. The default coordinate system for the thermal fields will be CS-IR. The coordinate system is indicated if this is not the case.

The internal heat sources s i (x,y,t) are obtained through an estimation of the left hand side of equation ( 6). The method, based on image processing, is explained and discussed in [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localization[END_REF][START_REF] Schlosser | Image processing to estimate the heat sources related to phase transformations during tensile tests of NiTi tubes[END_REF]. In these previous works, all the material coefficients involved in equation ( 6) were assumed to be constant throughout the test, whatever the state of the material.

In the following, we consider a case in which the thermal conductivity is not constant. As illustrated in Part 5, this is the case during superelastic deformation of a shape memory alloy, which induces a phase transformation between austenite A and martensite M. The conduction term of equations ( 5) and ( 6) for a thin sample is then:

div(-k grad ! " !!!! ! ) = -div(k grad ! " !!!! ! ) = -k ! 2 ! -grad ! " !!!! 2 k.grad 2 ! " !!!!! ! (7)
where the thermal conductivity is still assumed to be isotropic. Four types of thermal conductivity models can then be proposed:

Model 1: constant thermal conductivity throughout the test, independent of the position and the material state:

k=constant. ( 8
)
Model 2: spatially homogeneous thermal conductivity, but varying with time k=k(t). This is the case when the deforming material is a mixture of two phases A and M, with a martensite fraction f M evolving uniformly at the macro scale over the sample gauge length during its deformation. As a first rough approximation, the thermal conductivity of the two-phase material is written:

k=k(t)=k(f M (t))= f M (t) k M + (1-f M (t) )k A . ( 9
)
Therefore, in equation ( 6), the term -kΔ2θ must be reconsidered as:

-k(f M (t))! 2 ! (P,t) = -( f M (t)k M + (1" f M (t))k A )! 2 ! (P,t). ( 10
)
Model 3: spatially heterogeneous thermal conductivity, k=k(P,t), with values being either k A or k M . As presented in part 5, this is the case during the superelastic tensile test of an NiTi alloy for which localized deformation is associated with localized phase transformation. At each material point, in a first approach, two states are possible: either a low strain value associated with an almost fully austenitic A state or a high strain value associated with an almost fully martensitic M state:

k = k A if E xx < E AM k = k M if E xx ! E AM ! " # , ( 11
)
where E xx is the local axial component of the Green-Lagrange strain tensor measured by DIC and E AM is a fixed threshold indicating the A to M phase transition (see Section 5). The value given by equation [START_REF] Daly | Stress-induced martensitic phase transformation in thin sheets of nitinol[END_REF] will be used in equation ( 6) to locally compute the associated heat source.

Model 4: heterogeneous thermal conductivity, k= k(f M (P,t)). The local value of the martensite fraction f M must be known to be able to evaluate k . As recalled in [START_REF] Favier | Homogeneous and heterogeneous deformation mechanisms in an austenitic polycrystalline Ti-50.8at.%Ni thin tube under tension. Investigation via temperature and strain fields measurements[END_REF], this fraction can be computed from the specific phase transformation heat sources ! q tr = s i ! according to:

f M (P, t) = ! f M (P, u)du = ! q tr (u) !H A"M tensile 0 t # 0 t # du , (12)
where !H A"M tensile stands for the specific enthalpy for a complete martensitic transformation (f M = 1).

Equation ( 9), with space dependence of f M , and equation ( 7) are then used.

ANALYSIS OF LOCALIZED SUPERELASTIC TENSILE TESTS OF NITI SHAPE MEMORY

ALLOYS

The proposed heat source estimation method is now applied for analysis of the superelastic deformation of a NiTi shape memory alloy. In this paper, we decided to present only one tensile test, representative of several tests conducted on NiTi sheets. In the present work, only heat source estimations using Models 1 and 3 of Part 4 are presented. Model 2 was taken into account in [START_REF] Schlosser | Infuence of thermal and mechanical aspects on deformation behaviour of NiTi alloys[END_REF] in the case of an homogeneous superelastic shear test. Experimental results using Models 2 and 4 will be presented in a forthcoming paper.

Material and experimental procedure

The used material was a commercial NiTi alloy produced in China by the PeierTech company with a composition of Ti-50.6 at.% Ni. The material was annealed for 1 h at 1023 K and then aged at 773 K for 30 min. A specimen for differential scanning calorimetry (DSC) analysis was cut from the sheet using a low-speed diamond cut-off wheel. DSC measurements are shown in Figure 8 and exhibited standard AàM and MàA phase transformations upon cooling and heating, respectively. The transformation temperatures were determined to be M s =253 K, M f =230 K, A s =256 K and A f =272 K, and the enthalpies were ΔH A-M = 13 J/g, ΔH M-A =15 J/g. The initial microstructure of the sample was thus fully austenitic at room temperature.

The presented tensile test was performed on one rectangular sheet sample of gauge length L 0 =59 mm, 7.14 mm width and 0.5 mm thickness. It was conducted at room temperature (T 0 =298 K) and at constant crosshead velocity ! U / L 0 = 4.9 10 !4 s !1 . The experimental set-up was described in part 2.

Fig. 8 DSC measurements of the studied NiTi plate.

Heat source estimation using constant thermal conductivity (Model 1)

In order to analyze this tensile test, strain, temperature variation and heat source fields were synchronized using the method presented in Part 3. The synchronization was performed at the infrared frequency f IR on the CS-DIC (pixel size: 0.65 mm). The temperature variation fields used in equation ( 6) were computed with equation ( 2). The heat sources were calculated in the current CS-IR, then projected in the reference CS-DIC.

The nominal stress-strain curve is shown in Figure 9a. Figures 9c to 9f present the Green-Lagrange axial strain E xx , the temperature variations θ and the heat source s fields at four different times, as reported in Figure 9a. The area corresponding to the CS-DIC is indicated in a scheme in Figure 9b.

Model 1 was first used to estimate heat source fields, with the arbitrary choice k=k M . In that case, the heat sources are noted s M . At this step, these first results will be considered as qualitative estimations of the heat sources.

The test started with an almost linear initial variation in the stress, where positive temperature variation was observed [START_REF] Favier | Homogeneous and heterogeneous deformation mechanisms in an austenitic polycrystalline Ti-50.8at.%Ni thin tube under tension. Investigation via temperature and strain fields measurements[END_REF]. Then the stress-strain pattern was curving around time a. At this moment, the three fields plotted in Figure 9c were almost homogeneous, despite the fact that weak thermal gradients were noted with the typical associated colorbar. A first deformation localization was then observed at the beginning of the plateau (time b) on E xx and s M fields. This localization took the shape of a thin band, noted B1 in Figure 9d. The band angle was about 59° to the tensile axis. At the same time, due to heat diffusion, weak thermal heterogeneities were observed in the thermal field θ. At the middle of the plateau (time c), localization was very important (Figure 9e)

and seems very clear considering heat source fields. Strain field monitoring revealed two regions, the first one with small strain values (dark red color) lower than 0.01 and the second one with high strain values (clear yellow color) higher than 0.06. The strain field did not give the current position of the transformation. This was not the case for the heat source field, which showed the zone where substantial phase transformation occurred. At the end of the loading (time d, Figure 9f), the strain and heat sources were also localized. As shown by the strain fields, about 30% of the sample was not transformed and was still mostly in the austenitic phase. Temperature variations increased up to 11K but, as noted before, were weakly localized. Conversely, the heat sources were strongly localized in thin-phase transformation bands with, qualitatively, the same characteristics (including width, amplitude).

In order to obtain a global view during loading, a spatiotemporal representation of these fields is plotted in but was hidden by the grips. After time b, the middle band B1 became wider and the two bands B2 and B3 that initiated outside the observation area were also widening and therefore visible. Finally, three other bands (B4, B5 and B6) appeared in the middle of "cold areas" of the sample. The bands near the top grip B2, B4, B5 and the middle band B1 had the same orientation (≈59° to the loading axis). The bottom bands B3 and B6 were oriented in the opposite direction (≈-59°). Unlike NiTi tubes, where helical propagation was observed [START_REF] Schlosser | Image processing to estimate the heat sources related to phase transformations during tensile tests of NiTi tubes[END_REF][START_REF] Favier | Homogeneous and heterogeneous deformation mechanisms in an austenitic polycrystalline Ti-50.8at.%Ni thin tube under tension. Investigation via temperature and strain fields measurements[END_REF][START_REF] Schlosser | Infuence of thermal and mechanical aspects on deformation behaviour of NiTi alloys[END_REF], the bands were stopped on the left and right edges of the sheet sample and could only grow wider along the vertical axis.

As noted before, the strain fields were still strongly localized at the end of loading (time d) since the test was stopped before the end of the stress plateau. The unloading stress plateau was not observed in this test because the material was not superelastic at this temperature. A remaining strain of ε=1.8% was noted at the end of the test.

Localization could easily be detected using the temperature field at the very beginning of localization (e.g. temperature field at time b in Figure 10b). Soon after, due to conduction phenomena, temperature fields were difficult to analyze (see θ fields in Figures 9d,10b). The spatiotemporal representation of the heat sources (Fig- ure 10c) is more contrasted, showing the continual propagation (spreading) of the active fronts of the bands, with about the same characteristics: width, amplitude and spreading velocity. The heat source bands were localized at the same positions as the strain jumps. When no strain jump was observed, the heat source was weak. θ(X,Y 0 ,t) and c) the heat sources s M (X,Y 0 ,t). The two superimposed curves are the nominal stress σ and the local axial strain ε (refer to Figure 9a to get their amplitude values).

In order to show the influence of the value selected for thermal conductivity in Model 1 to estimate the heat source, heat sources s A and s M , respectively associated with k=k A and k=k M , were plotted for two material points A and B (Figure 11). The heat sources s A , obtained with higher k values give upper bound values and s M are lower bound values of the actual heat sources. Material point A was chosen in the middle of the sample, where the first band appeared; whereas B was chosen in an area that remained outside of the localization bands (see position in Figure 9d). Their strain patterns, E xx (A) and E xx (B), are plotted in Figures 11a and11b and the heat sources in Figures 11c to 11f.

For the two points, the strain pattern was the same before time a (Figures 11a,11b). Then, between times a and b, the strain increased quickly at point A, while it stayed almost constant at B. Just before time b, E xx (A) increased suddenly when the first deformation localization band appeared (Figure 9c). After time b and until time d, strain E xx (A) increased in the localized transformation area but with a smaller strain rate, from E xx (A)=0.055 to 0.07. At the same time, E xx (B) was still constant, about 0.01, in the outside region of the band. The heat source variations for the two points, plotted in Figures 11c and11d, confirm the previous strain observations. At point A, a sudden and intense peak appeared around time b, when the first transformation band appeared, then propagated in the direction of the mobile grip (see band B1 in Figure 10c). A second smooth peak appeared at the same place just before time c. It corresponds to the widening of the band in the other direction.

The heat sources at point B were positive but very weak, but an increase was observed at the end of loading when the transformation was close to point B.

The close-ups shown in Figures 11e and11f highlight the sensitivity of the heat source estimation with the thermal conductivity value. In areas where the heat sources were strongly localized (e.g. around point A at time b here) the laplacian term in equation ( 6) was important in the heat source estimation, compared to the other ones.

Hence, the heat sources are strongly dependent on k, which explains the marked difference in amplitude between s A and s M (Figure 11e). Some differences were also observed at the beginning of unloading at point A, where the material is supposed to be initially in the martensite phase (Figure 11f). Heat sources s A appeared to be exothermic whereas they are almost zero and then negative when calculated (s M ) with k=k M .

Heat source estimation using non-constant thermal conductivity (Model 3)

The two Figures 11a and 11b between the less deformed and highly deformed states is very short. In this interval, the phase changed rapidly from almost fully austenitic to mostly martensitic. At the end of this short jump, the strain value (about 5.5%)

indicates that the phase transformation was not complete. As reported by [START_REF] Miyazaki | Lüders-like deformation observed in the transformation pseudoelasticity of a TiNi alloy[END_REF], complete phase transformation would lead to nearly 8% strain. Hence, the threshold E AM associated with Model 3 (equation ( 10)) was arbitrarily set at a mean value of 0.04. As a first approximation, the material was considered to be fully martensitic above

Energy estimation (application to the case of thermal conductivity Model 3)

From s(P,t) it is possible, for each material point P(x,y,t) of initial coordinates (X,Y), to estimate the specific heat energies q(P,t) locally released during the tensile test by the following time integration:

q(X,Y, t) = 1 ! s(X,Y, u)du 0 t !
Such a local quantity can be considered as an in-situ local DSC of the material. The energy profiles, together with the axial strain and heat source profiles along the central column of the CS-DIC (see Figure 10) are reported in Figure 13 at the same time as. At the end of loading, time d, the energy profiles reveal two typical levels, associated with the same localization area of the sample where high strains have also been measured. The upper energy values, about 27 J g -1 , are associated with exothermic enthalpy of the martensite phase transformation !H A"M tensile . As noted before, when analyzing strain intensities, the transformation at time d was not complete.

However, the intensity (27 J g -1 ) is still much higher than typical literature values for such alloys obtained by temperature induced transformation (ΔH tr = 13 J g -1 in the DSC reported in Figure 8). Such results were already obtained by [START_REF] Favier | Restoration by rapid overheating of thermally stabilised martensite of NiTi shape memory alloys[END_REF][START_REF] Liu | Stability of ageing-induced multiple-stage transformation behaviour in a Ti-Ni 50.15% Alloy[END_REF] but during free stress reverse transformation (MàA) consecutive to an initial mechanical loading.

To show the influence of spatial synchronization in the energy estimation, the energy q 0 obtained from spatial heat sources s(x,y,t) in equation ( 13) was plotted at time d in Figure 13d. This energy profile is noisier than q and the positions of the energy upper levels do not at all correspond to the upper strain levels. profiles of heat energy q, plotted in Lagrangian configuration and obtained by a time integration of s(X,Y,t). The heat energy, denoted q 0 , plotted in black dotted lines, was obtained by temporal integration of s(x,y,t).

From Figures 10, 11 and 13, it is possible to quickly verify hypothesis (i) on the weakness of the convective term used in section 4.1 to establish eq. 6. The strain jump in the band is about 0.045 in about 10 s (Figure 11), so the mean axial strain rate in the band is about 0.0045 s -1 . The bandwidth is close to 4 pixels (Figure 13b), i.e.

2.6 mm. The axial velocity associated with the localized deformation in the band is then V=0.0045*2.6=0,0117 mm s -1 . The temperature increase in the band is about 2K (Figure 10b), such that !T !t " 2 /10 = 0.2 K s -1 and the projection of grad

! " !!!! T in the X direction is --1 mm 0.77 K = 2/2.6 T grad ≈ . The term ! v.grad
" ! """" T , estimated through 0.00117*0.77=0.009 K s -1 , is then about 20-fold lower than !T !t , which is in close agreement with hypothesis (i).

Finally, through a spatiotemporal representation, Figure 14 gives a synthetic view of the axial strain field E xx and the local DSC heat q. These two experimental fields obtained by two independent set-ups (visible and infrared cameras) are very similar. The positions of the phase transformation given indirectly by E xx and q are located at the same place in these two types of representation. Due to the difficulty of taking the thermal conductivity coefficient into account, the results of the local DSC proposed in this work can be considered as a semiquantitative result. However, the confrontation with the strain field is very encouraging and indirectly gives a kind of validation of the proposed method. Moreover, if the transformation strain Δε tr and the enthalpy !H tr tensile for a complete transformation are known, these two fields can also give an estimation of the martensite phase fraction [START_REF] Favier | Homogeneous and heterogeneous deformation mechanisms in an austenitic polycrystalline Ti-50.8at.%Ni thin tube under tension. Investigation via temperature and strain fields measurements[END_REF][START_REF] Louche | Kinematic fields and acoustic emission observations associated with the Portevin Le Châtelier effect on an AlMg alloy[END_REF][START_REF] Schlosser | Infuence of thermal and mechanical aspects on deformation behaviour of NiTi alloys[END_REF].

Heat source and energy estimations: comparison of Model 1/Model 3 thermal conductivities

The differences observed between models 1 and 3 in terms of heat source s and energy estimations q can be quantified from Figure 15, in which the results at time d obtained with Model 3 for thermal conductivity are compared with those obtained with Model 1. The different plots confirm that s and q values, as estimated from Model 3 for k, are located between the lower and upper bounds obtained with Model 1 for k, i.e. respectively k=k M =8 Wm -1 K -1 and k=k A =18 Wm -1 K -1 .

There are substantial differences between the two models, indicating the key role played by the thermal conductivity values of the two phases, for heat source or energy estimations during localized phase transformations.

At the end of loading (time d), the phase transformation energy varied from about 20 to about 30 J g -1 , which are bounds for the mean value (27 J g -1 ) given in the previous section with Model 3. Finally, this study highlighted the necessity of increasing the accuracy on thermal conductivity data for this kind of material. 11)).

CONCLUSION

A methodology to obtain heat source estimations, when thermal conductivity is not constant, was proposed in this paper. It was applied to estimate heat sources during a localized superelastic tensile test of an NiTi shape memory alloy displaying both stress-induced phase transformation and localized deformation behavior.

The same side of an NiTi SMA sheet sample was observed with a visible and an infrared camera. A specific spatiotemporal method was first presented in order to obtain, at the same time and at the same material point, kinematic (displacement, strain) and thermal (temperature, heat sources) data. This synchronization: (i) allowed temperature monitoring at a material point which can move during deformation, (ii) improved the temperature variation measurement, and (iii) also improved the estimation of heat energy by temporal integration of the heat sources at each material point.

As the thermal conductivity of martensite (k=k A ) and austenite phases (k=k M ) can be different, four models of k have been proposed and discussed. Two of them were chosen and used in two methods to estimate, from the thermal fields, heat sources during the test. Coefficient k was supposed to be constant in the first method and depended locally on the phase type in the second one.

Kinematic, thermal and heat source fields have been compared during the stress plateau of a superelastic tensile test. As in [START_REF] Brinson | Stress-induced transformation behavior of a polycrystalline NiTi shape memory alloy: micro and macromechanical investigations via in situ optical microscopy[END_REF][START_REF] Daly | Stress-induced martensitic phase transformation in thin sheets of nitinol[END_REF][START_REF] Pieczyska | Investigation of nucleation and propagation of phase transitions in NiTi SMA[END_REF][START_REF] Miyazaki | Lüders-like deformation observed in the transformation pseudoelasticity of a TiNi alloy[END_REF][START_REF] Shaw | Thermomechanical aspects of NiTi[END_REF][START_REF] Shaw | On the nucleation and propagation of phase transformation fronts in a NiTi alloy[END_REF], they showed substantial localization when considering the shape of the bands propagating along the sample. The axial strain fields had low (about 0.5%) and high (about 7%) values in areas where austenite and martensite phases were respectively predominant. The strain intensity indicates that martensitic transformation was not complete at the end of the loading. The temperature fields were very diffuse and hard to directly interpret. On the other hand, the heat sources clearly exhibited, at each current time, the area of the sample where the phase transformation occurred. These areas appeared to coincide with the transition area between low and high strain values.

Finally, heat sources were integrated at each material point over the time course in order to obtain a local estimation of heat energy associated with the phase transformation. Such measurement can be qualified as an in situ local DSC during a mechanical test. Due to the difficulty of modeling thermal conductivity, the result obtained can be considered as semi-quantitative. The scattering of specific heat values was not taken into account.

The same kind of technique, taking the appropriate value for each phase, could be applied to process thermal data, as carried out for thermal conductivity. Studies are under way to measure this material property for the two phases of NiTi materials.

However, the spatiotemporal distribution of the energy fields was in close agreement with the strain field distribution. A forthcoming study will be carried out to improve the heat source estimation, in particular through better determination of the heat capacities and thermal conductivities, and to validate it, in the case of heterogeneous thermal conductivity.

  As a complement to usual mechanical data (force F and cross-head displacement U), kinematic and thermal fields are measured by visible and infrared cameras. An example of such an experimental set-up is given in Figure 1.

Fig. 1

 1 Fig. 1 Experimental set-up. a) Position of the infrared and visible cameras monitoring the same surface of the sheet sample. b) Close-up of the DIC grid showing the resolution of the kinematic fields (grid step = 0.65 mm). c) DIC grid (in green), position of the two extreme points A and B. The distances L 0 and l 0 are the gauge length and the distance between A and B, respectively.

Fig. 2

 2 Fig. 2 Example of the different characteristics of the three used time bases.

Fig. 3

 3 Fig. 3 Schematic view of the sample during a tensile test in an initial configuration (a) and in the current configuration (b). The initial and current coordinates of material points A, B, P are (X,Y) and (x,y), respectively. c) Example of an experimental temporal variation pattern of horizontal (u y ) and vertical (u x ) displacements (in CS-IR) of points A and B. At the beginning, both points move together, due to rigid body motion.

Fig. 4

 4 Fig. 4 Visible image (a) and infrared image (b) of a tube sample and a steel ruler prior the test. A grey colormap was used. Three characteristic points (A, B and C) were selected for each image on the ruler. The two images were rotated by 180° due to the position of the two cameras.

Fig. 5

 5 Fig. 5 Schematic view of the spatial synchronization. a) Visible and b) infrared image at initial time t 0 =0. The initial positions of a material point P(0) in the two images are related by the transformation ir2v or v2ir . c)

Fig. 6

 6 Fig. 6 Scheme illustrating the spatial interpolation. The resolutions of the CS-DIC and CS-IR are indicated and show a higher resolution for the infrared grid. The initial and current grids, in CS-DIC, are projected in CS-IR using v2ir transformation. In the two cases, spatial interpolation (bi-linear or bi-cubic) was used to estimate the

Fig. 7

 7 Fig. 7 Illustration of spatiotemporal synchronization in the case of a real tensile test conducted on a NiTi tube sample in quasi isothermal conditions. Position of the tube (cross hatched) and of the processing area (red background color) at a) initial (free) and b) current (loading) time t. Material points A and B are also indicated in the two situations. c) and d) Green Lagrange axial strain field E xx and temperature variation field at current time t=37s, observed in CS-DIC, in Lagrangian description. e) Temperature variation field at the same time, observed in CS-DIC, in Eulerian description. Time-course patterns of f) the axial strain and g) the temperature variation

Figure 10 .

 10 Figure 10. It represents temporal changes in an axial profile of the Green-Lagrange axial strain E xx , the temperature variation θ and the heat source s M . The grey dotted line and the superimposed black continuous line respectively represent the time-course pattern of the nominal stress and the nominal axial strain ε.

Fig. 9

 9 Fig. 9 Comparison of axial strains E xx , temperature variations θ and heat sources fields s M at several times, plotted in reference CS-DIC. a) Nominal stress strain curve and position of the selected times during loading. b) Scheme of the area corresponding to CS-DIC. c) Green Lagrange axial strain E xx (X,Y,t), temperature variation θ(X,Y,t) and heat source s M (X,Y,t) fields at the same time a, indicated in the first Figure. d) to f) same fields at times b, c, d

Fig. 10

 10 Fig. 10 Spatiotemporal representations along the tensile axis (indicated by the green dotted line of the fixed Y 0 coordinate in the left schemes) of: a) the Green Lagrange axial strain E xx (X,Y 0 ,t), b) the temperature variation

Fig. 11

 11 Fig. 11 Temporal variations in Green Lagrange axial strain E xx at a) material point A and b) material point B; see positions of points A and B in Figure 9d. c) to f) Temporal variations in heat sources s A and s M , computed with the thermal conductivity k A and k M.

  have been merged into Figure 12a in order to analyze the temporal pattern of E xx for material points A and B. This representation shows, for material point A, that the time interval (denoted D)

E

  AM and k=k M was used to compute the heat sources within this strain range (now denoted s). The two regions E xx <E AM and E xx >E AM are indicated in Figures 12b on a spatiotemporal representation along the vertical axis of the sample.

Fig. 12 a

 12 Fig. 12 a) Temporal variations in Green Lagrange axial strain E xx (A) and E xx (B), for material points A and B. The threshold E AM and the length D of the transformation in A are also indicated. b) Background image showing the spatiotemporal pattern of E xx values plotted, in Lagrangian configuration, with a two-level colorbar defined by the E AM =0.04 threshold. The two superimposed curves are the nominal stress σ and the local axial strain ε (refer to Figure 9a to get their amplitude values). The two horizontal dotted green lines indicate the spatial position of points A and B.

Fig. 13 a

 13 Fig. 13 a) Nominal stress-strain curve and position of selected times during loading. b) Axial profiles, at the selected times, of the Green Lagrange axial strain E xx , plotted in Lagrangian configuration. c) Axial profiles of heat sources s estimated with Model 3 of thermal conductivity, plotted in the same Lagrangian configuration. d) Axial

Fig. 14

 14 Fig. 14 Spatiotemporal representation along the tensile axis, indicated by the green dotted line of the fixed Y 0 coordinate in the left schemes, of: a) the Green Lagrange axial strain E xx ,(X,Y0,t) and b) of the heat energy q(X,Y,t). The two curves superimposed in these figures are the nominal stress σ and the local axial strain ε (refer to Figure 9a to get their amplitude values).

Fig. 15

 15 Fig. 15 Comparison at time d (see Figure 13a) of heat source (a) and energy (b) estimations from the two thermal conductivity models: Model 1 (with k=k M or k=k A ) and Model 3 (eq. (11)).