Xiaojun Gu
email: gu.xiaojun@live.cn

Wael Zaki
email: wael.zaki@kustar.ac.ae

Claire Morin
email: claire.morin@emse.fr

Ziad Moumni
email: ziad.moumni@ensta-paristech.fr

Weihong Zhang
email: zhangwh@nwpu.edu.cn

Time Integration and Assessment of a Model for Shape Memory Alloys Considering Multiaxial Nonproportional Loading Cases

Keywords: shape memory alloys, nonproportional loading, multiaxial loading, phase transformation, martensite reorientation, multisurface inelasticity

The paper presents a numerical implementation of the ZM model for shape memory alloys (Zaki and Moumni, 2007a) that fully accounts for nonproportional loading and its influence on martensite reorientation and phase transformation. Detailed derivation of the time-discrete implicit integration algorithm is provided, including an explicit closed-form expression for the continuous material jacobian. The algorithm is used for finite element simulations using Abaqus, in which the model is implemented by means of a user material subroutine. Extensive validation of the model is provided against multiple sets of experimental and numerical simulation data taken from the literature.

Introduction

Shape memory alloys (SMAs) are capable of sustaining severe inelastic deformation that can be recovered by heating [START_REF] Otsuka | Shape Memory Materials[END_REF]. This behavior is explained by the ability of SMAs to undergo transformation between a higher symmetry austenite phase and a lower symmetry martensite phase. The martensite phase can be inelastically deformed by detwinning and reorientation of variants characterized by different crystallographic orientations [START_REF] Funakubo | Shape memory alloys[END_REF]. Starting from the 1970's, the use of SMAs in engineering applications has seen significant development [START_REF] Humbeeck | Non-medical applications of shape memory alloys[END_REF] and it now spans a number of important fields ranging from biomedicine to civil engineering and aeronautics [START_REF] Stöckel | Nitinol medical devices and implants[END_REF]. A key feature of many shape memory alloys is superelasticity, which refers to the capacity of these materials to undergo substantial deformation when subjected to mechanical loading at temperatures exceeding the so-called "austenite finish" temperature and to recover their undeformed shape once the load is removed. Finite element analysis of this superelastic behavior in cases of uniaxial and multiaxial proportional and nonproportional isothermal loading is the focus of this paper.

Uniaxial or multiaxial loading of shape memory alloys can induce phase transformation accompanied by the nucleation or shrinking of martensite variants [START_REF] Otsuka | Shape Memory Materials[END_REF][START_REF] Bodaghi | A phenomenological SMA model for combined axialtorsional proportional/non-proportional loading conditions[END_REF]. When the loading direction changes, preferred martensite variants begin to form at the expense of others leading to a change in orientation of the inelastic strain in the martensitic phase. Both phase change and martensite reorientation can take place at the same time in presence of multiaxial nonproportional loading, which motivates the need for SMA models capable of accounting for these processes simultaneously. Numerical integration of the constitutive equations in this case is particularly challenging because it typically involves solving highly nonlinear equations while satisfying multiple consistency conditions and intrinsic constraints on the state variables. Proper numerical integration of SMA models is of special importance in such areas as the mechanics of fracture of SMAs, where multiaxial loading conditions and severe stress gradients are experienced by the material in the neighbourhood of the crack tip, as well as in analyzing fatigue of SMA structures, where accurate simulation of the multiaxial behavior of the material is needed for numerical prediction of fatigue life using appropriate failure criteria.

In the literature, the behavior of SMAs subjected to nonproportional loading has been the subject of intense investigation. [START_REF] Sittner | Experimental study on the thermoelastic martensitic transformation in shape memory alloy polycrystal induced by combined external forces[END_REF] performed biaxial tension-torsion experiments on polycrystalline CuAlZnMn SMA samples considering rectangular and triangular stress-and straincontrolled loading cases. In both cases, the inelastic deformation of the material was found to be fully recoverable by unloading. [START_REF] Lim | Mechanical Behavior of an Ni-Ti Shape Memory Alloy Under Axial-Torsional Proportional and Nonproportional Loading[END_REF] reported several biaxial proportional and nonproportional loading experiments for thin-wall tubes of superelastic NiTi SMAs under circular normal/shear stress-and strain-controlled loading. The material response was found to exhibit tensile-compressive asymmetry when subjected to cyclic tensile-compressive loading. Moreover, force-and strain-controlled experiments were conducted at different strain rates, with and without hold times, and the results were reported including the influence of thermomechanical coupling. [START_REF] Bouvet | Mechanical Behavior of a Cu-Al-Be Shape Memory Alloy Under Multiaxial Proportional and Nonproportional Loadings[END_REF] carried out tension/internal pressure tests and biaxial compression tests using CuAlBe SMA samples where the authors also observed tensile compressive asymmetry in the material behavior.

A model was later proposed [START_REF] Bouvet | A phenomenological model for pseudoelasticity of shape memory alloys under multiaxial proportional and nonproportional loadings[END_REF]) that takes into account this asymmetry, as well as the influence of temperature, return point memory and the nonproportionality of the applied load on SMA behavior. The model was validated against experimental data obtained under general multiaxial loading cases. [START_REF] Panico | A three-dimensional phenomenological model for martensite reorientation in shape memory alloys[END_REF] proposed a macroscopic phenomenological model based on the framework of thermodynamics of irreversible processes. The model accounts for the influence of multiaxial stress states and non-proportional loading histories and was shown to reasonably agree with the experimental data in [START_REF] Sittner | Experimental study on the thermoelastic martensitic transformation in shape memory alloy polycrystal induced by combined external forces[END_REF]. [START_REF] Grabe | Path dependence and multiaxial behavior of a polycrystalline NiTi alloy within the pseudoelastic and pseudoplastic temperature regimes[END_REF] reported experimental results for polycrystalline NiTi subjected to several multiaxial loading cases in a wide temperature range in order to investigate superelasticity and the one-way shape memory effect of the material. Based on these experiments, the use of von Mises equivalent stress in deriving loading functions for phase transformation and martensite reorientation was found inappropriate as it failed to capture tensilecompressive and tensile-torsional asymmetries. [START_REF] Arghavani | A 3-D phenomenological constitutive model for shape memory alloys under multiaxial loadings[END_REF] presented a phenomenological constitutive model that accounts for key features of SMA behavior when the material is subjected to proportional as well as non-proportional loading. The model was successfully used to simulate a number of experimental results taken from the literature. [START_REF] Chemisky | Constitutive model for shape memory alloys including phase transformation, martensitic reorientation and twins accommodation[END_REF] derived a model for shape memory alloys based on earlier work by the authors. The model focuses on SMA behavior at lower stress levels and takes into account several effects associated with SMA behavior such as tensile-compressive asymmetry and temperature-driven phase transfor-mation. Simulation results are reported, without validation, at different temperatures for a case of multiaxial loading consisting of a square tensionshear loading path. Another model was developed by [START_REF] Stebner | Explicit finite element implementation of an improved three dimensional constitutive model for shape memory alloys[END_REF]. The model was numerically integrated using an explicit scheme.

Care was taken in eliminating the need for user-calibration of the numerical integration parameters and reducing the sensitivity of the integration algorithm to mass scaling for faster computation. [START_REF] Bodaghi | A phenomenological SMA model for combined axialtorsional proportional/non-proportional loading conditions[END_REF] derived a simple and robust phenomenological model for SMAs using the framework of continuum thermodynamics of irreversible processes. The model was used to simulate SMAs subjected to uniaxial and biaxial loading, including proportional and nonproportional tensile-torsional loading. Good agreement with experimental data was achieved for these loading cases. [START_REF] Auricchio | Theoretical and numerical modeling of shape memory alloys accounting for multiple phase transformations and martensite reorientation[END_REF] recently proposed a model for SMAs that accounts for several aspects of SMA behavior in presence of multiaxial loading. The model was numerically integrated using Fischer-Burmeister functions to account for the Kuhn-Tucker conditions governing the evolution of several state variables, including the volume fractions of single-and multi-variant martensites and the directional orientation of single-variant martensite. The use of the Fischer-Burmeister functions allowed the authors to dispense with the conventional search for active loading surfaces by substituting the Kuhn-Tucker inequalities with a set of equivalent nonlinear equalities. The model is validated against experimental data for several loading cases, including a case of combined tension-torsion involving a helical SMA spring. A recent model was also proposed by [START_REF] Mehrabi | A thermodynamicallyconsistent microplane model for shape memory alloys[END_REF], which consisted of a thermodynamically acceptable evolution of the microplane model for SMAs. In the spirit of the microplane theory, the authors considered different stress projection methods to resolve the three-dimensional constitutive material behavior into a combination of uniaxial models defined in specific spatial directions.

The model was then used to simulate experiments on a thin hollow SMA cylinder subjected to tension, torsion and combined tension-torsion loading cases. Only proportional loading was considered for the combined loading case Zaki and Moumni (2007b,a) and [START_REF] Moumni | Theoretical and numerical modeling of solidsolid phase change: Application to the description of the thermomechanical behavior of shape memory alloys[END_REF] derived a phenomenological model that can simultaneously account for phase transformation and martensite detwinning and reorientation. The model was improved to account for tensile-compressive asymmetry [START_REF] Zaki | Modeling tensile-compressive asymmetry for superelastic shape memory alloys[END_REF][START_REF] Zaki | An approach to modeling tensile-compressive asymmetry for martensitic shape memory alloys[END_REF], plastic yielding [START_REF] Zaki | An extension of the ZM model for shape memory alloys accounting for plastic deformation[END_REF], and thermomechanical coupling (Morin et al., 2011a) including cyclic effects (Morin et al., 2011b;[START_REF] Moumni | Cyclic behavior and energy approach to the fatigue of shape memory alloys[END_REF]Morin et al., 2011c;[START_REF] Morin | A comprehensive approach for fatigue analysis of shape memory alloys[END_REF]. However, numerical integration of the model in presence of non-proportional multiaxial loading was only considered so far in the case of martensitic SMAs (Zaki, 2012a,b). This paper presents detailed integration for the Zaki-Moumni (ZM) model for SMAs capable of undergoing combined phase transformation and martensite reorientation when subjected to nonproportional multiaxial loading.

Proper numerical integration and extensive validation of the model for such general loading cases has not been attempted before and therefore constitutes the novelty of this work. The implicit time integration procedure is implemented by means of a user material subroutine into the finite element software Abaqus.

The paper provides a brief review of the ZM models in section 2 followed by a general presentation of the boundary value problem and time-discrete equations for a SMA subjected to isothermal loading in section 3, including a closed-form expression for the continuous material Jacobian. Numerical simulations are then carried out in section 4 for a number of loading cases, many of which are taken from the literature and shown to be properly simulated by the model. A general conclusion and outlook are finally provided in section 5.

Summary of the analytical model

The ZM model for shape memory materials is developed based on the framework of generalized standard materials [START_REF] Halphen | Plastic and visco-plastic materials with generalized potential[END_REF] that was slightly extended to account for constraints on the state variables [START_REF] Moumni | Sur la modélisation du changement de phase à l'état solide[END_REF][START_REF] Moumni | Theoretical and numerical modeling of solidsolid phase change: Application to the description of the thermomechanical behavior of shape memory alloys[END_REF]. Using standard index notation for tensor representation and implied summation over repeated indices, the derivation of the model gives the following stress-strain relation:

σ ij = K ijkl (ε kl -zε ori kl ) (1)
where σ ij and ε kl are the stress and total strain tensors, z is the volume fraction of martensite and ε ori kl is the inelastic strain caused by martensite detwinning and reorientation. The elastic stiffness tensor K ijkl is a function of phase composition. It is given by the expression

K ijkl = (1 -z)K A,-1 ijkl + zK M,-1 ijkl -1 , (2)
in which K A ijkl and K M ijkl are the elastic stiffness matrices of austenite and martensite respectively.

The volume fraction of martensite cannot be less than zero or greater than one and the recoverable inelastic deformation of martensite is limited by a maximum ε 0 that depends on the material used. These constraints are mathematically expressed as z ≥ 0, 1 -z ≥ 0, and ε 0 -ε ori eq ≥ 0 (3) where ε ori eq is an equivalent orientation strain taken in this paper to be

ε ori eq = 2 3 ε ori ij ε ori ij . ( 4
)
The evolution of the dissipative variables z and ε ori ij is governed by the Kuhn-Tucker conditions

F 1 z σ ij , z, ε ori ij ≤ 0, ż ≥ 0, żF 1 z = 0, (5)
for forward transformation, in which 0 ≤ z < 1;

F 2 z σ ij , z, ε ori ij ≤ 0, -ż ≥ 0, -żF 2 z = 0, (6)
for reverse transformation, in which 0 < z ≤ 1;

F ori σ ij , z, ε ori ij ≤ 0, η ≥ 0, ηF ori = 0, (7)
for martensite detwinning and reorientation, in which η is the inelastic multiplier and the rate of orientation strain is provided by the normality rule

εori ij = η ∂F ori ∂X ij = 3 2 η X ij X vm = ηN ij . (8)
In the above equation, X vm is the von Mises equivalent of the thermodynamic force X ij , conjugate to the orientation strain ε ori ij , and

N ij = 3 2 X ij
Xvm is a vector indicating the direction of the orientation strain rate εori ij in strain space. To simplify numerical integration, martensite is assumed to be fully oriented by the applied stress as soon as the phase transformation from austenite to martensite takes place. This simplification is in accordance with results in the literature [START_REF] Patoor | Shape memory alloys, part I: general properties[END_REF]. For the ZM model used, this assumption is equivalent to considering that the stress σ rf needed for complete detwinning of martensite is less than the stress σ ms required to initiate forward phase transformation. to indicate differentiation with respect to space coordinate i, the problem consists in solving the following set of equations over Ω for t ∈ [0, T ]:

• Static equilibrium σ ij,j + f i = 0 in Ω, σ ij n j = T d i over ∂Ω T , (9)
where n j is an outward unit vector normal to the boundary ∂Ω T ,

• Kinematic boundary conditions

ξ i = ξ d i over ∂Ω ξ , (10)
where ξ i is the displacement vector,

• Strain compatibility

ε ij = 1 2 (ξ i,j + ξ j,i ) in Ω, (11)
• Constitutive equations

σ ij = K ijkl ε kl -zε ori kl in Ω (12)
where K ijkl is given by (2), and z and ε ori ij are governed by the conditions (3) to (8).

The boundary value problem above is solved using the finite element method. Time integration of the constitutive equations is done using an implicit algorithm that takes into account the possibility of simultaneous phase transformation and martensite detwinning and reorientation in presence of nonproportional multiaxial loading.

The assumption of complete orientation of martensite as soon as it forms reduces the constraint on the orientation strain to the equality ε 0 -ε ori eq = 0. The loading function F ori governing the orientation of martensite can be written in this case as

F ori = X vm -zY where X ij = s ij - 2 3ε 2 0 s kl ε ori kl ε ori ij , (13)
in which s ij is the stress deviator and Y is a material parameter that defines the stress onset of martensite detwinning in uniaxial tension. The expressions of the phase transformation functions remain unchanged with respect to the original ZM model. They are given by

F 1 z = 1 2 3 2 El M A s ij s ij + P M A σ 2 ii -C(T ) + s ij ε ori ij -(G + b)z -a(1 -z) -(α -β)z + β 2 2 3 ε ori ij ε ori ij , (14)
F 2 z = - 1 2 3 2 El M A s ij s ij + P M A σ 2 ii + C(T ) -s ij ε ori ij + (G -b)z -a(1 -z) + (α -β)z + β 2 2 3 ε ori ij ε ori ij , (15)
where El M A , P M A , α, β, a, b, G are material parameters and C(T ) is a function of temperature.

Considering a discretization of the time interval [0, T ] into N intervals, the increments of the state variables z and ε ori ij for load increment n + 1, n ∈ {0, N }, are determined by enforcing the appropriate consistency conditions assuming z and ε ori ij are known at load increment n. In the equations below, the load increments n and n + 1 are indicated in the subscript and the solution iteration k is indicated in the superscript of the relevant variable.

The incremental equations are solved using a classical predictor-corrector strategy [START_REF] Simo | Computational Inelasticity[END_REF] as follows:

• An elastic stress increment is first determined that achieves static equilibrium:

z (0) n+1 = z n and ε ori,(0) ij,n+1 = ε ori ij,n (16) σ (0) ij,n+1 = K (0) ijkl,n+1 ε kl,n+1 -z (0) n+1 ε ori,(0) kl,n+1 (17)
• The stress increment is then corrected and the state variables are updated in accordance with the consistency conditions. If F z is any of the two loading functions F 1 z and F 2 z governing phase change, these conditions must satisfy

F z σ (k) ij , z (k) , ε ori,(k) ij ≤ 0 (18) F ori σ (k) ij , z (k) , ε ori,(k) ij ≤ 0 (19)
where the subscript n + 1 is eliminated for readability and where each inequality reduces to a strict equality when the relevant process is active. The equalities can be approximated using first-order Taylor series expansions of the loading functions F z and F ori in the neighbourhood of

σ (k-1) ij , z (k-1) , ε ori,(k-1) ij
, which can then be solved iteratively. The Taylor series approximations are written explicitly as follows:

F (k-1) z + ∂F z ∂σ ij (k-1) ∆σ (k) ij + ∂F z ∂z (k-1) ∆z (k) + ∂F z ∂ε ori ij (k-1) ∆ε ori,(k) ij = 0 (20) F (k-1) ori + ∂F ori ∂σ ij (k-1) ∆σ (k) ij + ∂F ori ∂z (k-1) ∆z (k) + ∂F ori ∂ε ori ij (k-1) ∆ε ori,(k) ij = 0 (21)
where ∆v

(k) = v (k) -v (k-1)
is the increment of variable v at iteration k.

If phase transformation and martensite orientation are active at the same time, ∆z (k) and η (k) are obtained by solving the time-discrete equations ( 20) and ( 21) simultaneously, where η (k) is the discrete inelastic multiplier used to express the discrete normality rule at iteration k as

∆ε ori,(k) ij = η (k) N (k) ij . (22)
N ij in the previous equation is approximated at iteration k by its value at the previous iteration, with the initial value corresponding to an elastic prediction. Substituting ( 22) into ( 20) and ( 21) gives

F (k-1) z + ∂F z ∂σ ij (k-1) ∆σ (k) ij + ∂F z ∂z (k-1) ∆z (k) + η (k) ∂F z ∂ε ori ij (k-1) N (k) ij = 0, (23)
F (k-1) ori + ∂F ori ∂σ ij (k-1) ∆σ (k) ij + ∂F ori ∂z (k-1) ∆z (k) + η (k) ∂F ori ∂ε ori ij (k-1) N (k) ij = 0. ( 24
)
The time-discrete stress-strain relation can then be written as ∆σ

(k) ij = K (k-1) ijkl ∆ε (k) kl -∆K -1,(k) ijkl σ (k-1) ij + ∆z (k) ε ori,(k-1) kl + z (k-1) ∆ε ori,(k) kl (25)
where ∆ε

(k)
kl = 0 because ε ij is not updated in the corrector step and the increment of the elastic compliance tensor is given by

∆K -1,(k) ijkl = K M,-1 ijkl -K A,-1 ijkl ∆z (k) . ( 26
)
Combining equations ( 22), ( 25) and ( 26) gives ∆σ

(k) ij = -K (k-1) ijkl R (k-1) kl ∆z (k) + z (k-1) η (k) N (k) kl (27)
where

R (k-1) kl = K M,-1 ijkl -K A,-1 ijkl σ (k-1) ij + ε ori,(k-1) kl (28)
The increments of the martensite volume fraction ∆z (k) and the multiplier η (k) can now be calculated as

∆z (k) = A (k) ori F (k-1) ori -B (k) ori F (k-1) z A (k) z B (k) ori -B (k) z A (k) ori , (29)
η (k) = B (k) z F (k-1) z -A (k) z F (k-1) ori A (k) z B (k) ori -B (k) z A (k) ori , (30)
where

A (k) z = ∂F z ∂z (k-1) -K (k-1) ijkl ∂F z ∂σ ij (k-1) R (k-1) kl , (31)
A (k) ori = ∂F z ∂ε ori ij (k-1) N (k) ij -z (k-1) K (k-1) ijkl ∂F z ∂σ ij (k-1) N (k) kl , (32)
B (k) z = ∂F ori ∂z (k-1) -K (k-1) ijkl ∂F ori ∂σ ij (k-1) R (k-1) kl , (33)
B (k) ori = ∂F ori ∂ε ori ij (k-1) N (k) ij -z (k-1) K (k-1) ijkl ∂F ori ∂σ ij (k-1) N (k) kl , (34)
For phase transformation without martensite reorientation, the increment of martensite volume fraction at iteration k is given by the simple equation

∆z (k) = - F (k-1) z A (k) z . ( 35
)
Similarly, the inelastic multiplier at iteration k in the case of martensite reorientation without phase transformation is given by

η (k) = - F (k-1) ori B (k) ori (36)

Derivation of a continuous material Jacobian

For small volume changes, the continuous material Jacobian at the beginning of increment n + 1 is given by

J ijkl = dσ ij dε kl (37)
where the right-hand term is evaluated at the end of increment n.

The time-continuous stress-strain relation can be written in differential form as

dσ ij = K ijkl [dε kl -(R kl dz + ηzN kl )] (38)
For combined phase transformation and martensite reorientation, substituting dz and η in the above with the expressions

dz = A ori K ijkl ∂F ori ∂σ ij -B ori K ijkl ∂Fz ∂σ ij A z B ori -B z A ori dε kl , (39)
η = B z K ijkl ∂Fz ∂σ ij -A z K ijkl ∂F ori ∂σ ij A z B ori -B z A ori dε kl , (40)
leads to the following closed-form expression of the continuous material Jacobian:

J ijkl = dσ ij dε kl = K ijkl -K ijpq R pq A ori K rskl ∂F ori ∂σ rs -B ori K rskl ∂F z ∂σ rs A z B ori -B z A ori -zK ijpq N pq B z K rskl ∂F z ∂σ rs -A z K ijrs ∂F ori ∂σ rs A z B ori -B z A ori (41)
The same procedure gives the following expressions of the material Jacobian for the case of phase transformation with no reorientation of martensite and for the case of martensite reorientation without phase transformation:

J ijkl = dσ ij dε kl = K ijkl -K ijpq R pq K rskl ∂F z ∂σ rs A z (42)
and

J ijkl = dσ ij dε kl = K ijkl -zK ijpq N pq K rskl ∂F ori ∂σ rs B ori (43)
respectively.

Algorithmic setup

A "User Material" subroutine is developed for the simulation of phase transformation and martensite reorientation in the finite element analysis software Abaqus. The subroutine uses the implicit integration algorithm presented below in order to update the stress components and the state variables at the end of each load increment.

1. Read the parameters of the model, 2. Determine the elastic trial state:

• Calculate the increment of trial stress σ (0) ij,n+1 using (17),

• Calculate initial trial values for the loading functions

F (0) z,n+1 , F (0) ori,n+1 ; 3. If F (0) z,n+1 ≤ 0 and F (0)
ori,n+1 ≤ 0, the trial solution is acceptable:

• Set the stress vector equal to the trial stress vector,

• Set the state variables to their values at the beginning of the increment,

• Set the material Jacobian J ijkl equal to the elastic stiffness tensor

K ijkl ; 4. Else if F (0) z,n+1 ≤ 0 and F (0)
ori,n+1 > 0, (only martensite reorientation is active):

• Calculate the increment of the inelastic multiplier η n+1 using (30)

• Update the state variables,

• Calculate the material Jacobian J ijkl using (43);

Else if

F (0) z,n+1 > 0, F (0)
ori,n+1 ≤ 0 (only phase transformation is active):

• Calculate the increment of martensite volume fraction ∆z n+1 using (29),

• Update the state variables,

• Calculate the material Jacobian J ijkl using (42);

6. Else if F (0) z,n+1 > 0 and F (0)
ori,n+1 > 0 (martensite reorientation or phase transformation or both):

• Calculate the increment of the inelastic multiplier η n+1 and the increment of martensite volume fraction ∆z n+1 using ( 29) and (30),

• (a) If (∆z n+1 > 0 and F 1,(0) z,n+1 > 0 and η n+1 > 0) or (∆z < 0 and F 2,(0) z,n+1 > 0 and η n+1 > 0), then (simultaneous phase transformation and martensite reorientation):

-Update the stress vector and the state variables, -Calculate the material Jacobian J ijkl using Eq.( 41 -Reset the state variables to their values at the beginning of the increment, -Calculate the increment of martensite volume fraction ∆z n+1 using (29), -Update the state variables, -Calculate the material Jacobian J ijkl using (42).

Results and discussions

The problem defined in the previous section is solved by means of the finite element method for several superelastic SMA structures subjected to uniaxial and multiaxial proportional and nonproportional loading. The section is organized as follows:

• Proportional loading is considered in section 4.1, which includes numerical simulation results for pure tension, pure torsion, as well as multiaxial tension-torsion and tension-internal pressure loading that are compared to experimental data from [START_REF] Bouvet | A phenomenological model for pseudoelasticity of shape memory alloys under multiaxial proportional and nonproportional loadings[END_REF].

• A hypothetical situation involving nonproportional loading is presented in section 4.2 in order to investigate the ability of the model to account for multiaxial loading in presence of simultaneous phase transformation and martensite reorientation.

• In section 4.3, the model is used to simulate a number of numerical as well as experimental results available in the literature for SMAs subjected to nonproportional multiaxial loading. The reference data used for the validation are those reported by [START_REF] Sittner | Experimental study on the thermoelastic martensitic transformation in shape memory alloy polycrystal induced by combined external forces[END_REF], [START_REF] Bouvet | Mechanical Behavior of a Cu-Al-Be Shape Memory Alloy Under Multiaxial Proportional and Nonproportional Loadings[END_REF], [START_REF] Stebner | Explicit finite element implementation of an improved three dimensional constitutive model for shape memory alloys[END_REF], [START_REF] Panico | A three-dimensional phenomenological model for martensite reorientation in shape memory alloys[END_REF], [START_REF] Arghavani | A 3-D phenomenological constitutive model for shape memory alloys under multiaxial loadings[END_REF], [START_REF] Grabe | Path dependence and multiaxial behavior of a polycrystalline NiTi alloy within the pseudoelastic and pseudoplastic temperature regimes[END_REF]. In each case, the parameters used for the simulations using the ZM model were obtained by curve fitting to a subset of the reported data.

Simulation of proportional loading cases

The four loading cases shown in figure 1 The simulated uniaxial response of the material is shown in figures 2(a) cases, tensile and shear superelasticity is obtained with full inelastic strain recovery when the load is removed. The model is now used to simulate the behavior of a thin-wall SMA cylinder following [START_REF] Bouvet | A phenomenological model for pseudoelasticity of shape memory alloys under multiaxial proportional and nonproportional loadings[END_REF]. The cylinder is subjected to proportional multiaxial loading consisting of a combination of axial loading and internal pressure. The material parameters used for the simulation are reported in

for
A 0 f 300 K T 340 K
Table 2: material parameters obtained using data form [START_REF] Bouvet | A phenomenological model for pseudoelasticity of shape memory alloys under multiaxial proportional and nonproportional loadings[END_REF].

The loading considered in this case corresponds to the hoop versus axial stress shown in figure 5, in which a constant slope is obtained because of the proportional nature of the load.

The resulting stress-strain behavior of the material is shown in figure 6 where two superelastic stress-strain loops are shown representing the material behavior in the axial and hoop directions respectively. The simulations are in good agreement with the experimental data in both cases. In particular,

Multiaxial reorientation tests

A preliminary numerical simulation test is performed in this section in order to illustrate the ability of the model to simulate simultaneous phase transformation and martensite reorientation and to account for martensite reorientation at saturation, in which case the equivalent orientation strain reaches its maximum value ε 0 . In case the loading magnitude continues to increase beyond the saturation point, the components of the orientation strain tensor ε ori ij may continue to evolve to accommodate possible modifications in the loading direction while the equivalent orientation strain remains constant.

The parameters used in this section are those listed in table 1. The proposed simulation involves an axial loading corresponding to a maximum normal stress of 450 MPa followed by an increase in shear stress from 0 to 200 MPa as shown in figure 7(a). The maximum normal stress is such that maximum inelastic strain is achieved in tension (complete transformation of austenite to fully oriented martensite) as shown in figure 7(b). The superposition of shear stress leads to an increase in the yz strain component at the expense of the normal xx strain in such way that the total equivalent orientation strain remains constant (same figure). The evolution of the inelastic strain components is an expression of the physical process of martensite reorientation, in which preferred martensite variants form at the expense of other variants that are less-favorably oriented with respect to the applied load. From a phenomenological point of view, within the frame of the ZM model for shape memory alloys, this behavior is captured through the use of the loading function for martensite orientation F ori .

Further testing of the model is carried out using the square loading path shown in figure 8(a), which consists of the following loading steps starting from an unloaded state: Despite the total inelastic strain being recoverable by unloading, individual strain components do not vanish if the corresponding stress components are unloaded. Rather, the components of the inelastic strain tensor are scaled down in magnitude when the volume fraction of martensite decreases. This is consistent with the physics of superelastic deformation in SMAs, where shape recovery results from the transformation of martensite into austenite, not from the inelastic deformation being recovered within the martensite phase.

Validation of the model

The model is used to simulate a number of cases reported in the literature involving superelastic shape memory alloys subjected to nonproportional loading. Both stress-and strain-controlled experiments are considered using data obtained by several research groups including [START_REF] Sittner | Experimental study on the thermoelastic martensitic transformation in shape memory alloy polycrystal induced by combined external forces[END_REF], [START_REF] Bouvet | Mechanical Behavior of a Cu-Al-Be Shape Memory Alloy Under Multiaxial Proportional and Nonproportional Loadings[END_REF], [START_REF] Stebner | Explicit finite element implementation of an improved three dimensional constitutive model for shape memory alloys[END_REF], [START_REF] Panico | A three-dimensional phenomenological model for martensite reorientation in shape memory alloys[END_REF], [START_REF] Arghavani | A 3-D phenomenological constitutive model for shape memory alloys under multiaxial loadings[END_REF], [START_REF] Grabe | Path dependence and multiaxial behavior of a polycrystalline NiTi alloy within the pseudoelastic and pseudoplastic temperature regimes[END_REF].

Strain-controlled tension-shear loading

Following the recently published work of [START_REF] Stebner | Explicit finite element implementation of an improved three dimensional constitutive model for shape memory alloys[END_REF] the displacement-driven loading shown in figure 10(a). The loading involves nonproportional tensile-shear strain in which a maximum tensile strain of 0.04 is achieved. The material parameters used for the simulation, equivalent to those in [START_REF] Stebner | Explicit finite element implementation of an improved three dimensional constitutive model for shape memory alloys[END_REF], are listed in table 3. The temperature used for the simulation is artificial and is only meant to guarantee superelastic behavior of the shape memory alloy.

The displacement-controlled loading results in square-shaped shear strain versus normal strain behavior. The stress curve in figure 10(b) has a more complicated shape, which underlines the well-established nonlinearity in SMA behavior. Numerical simulations are in agreement with those obtained by [START_REF] Stebner | Explicit finite element implementation of an improved three dimensional constitutive model for shape memory alloys[END_REF]. A better fit can likely be obtained by fine-tuning the parameters of the ZM model. Some deviation is inevitable however because of intrinsic differences between the two models. simulation results of the two models include a steeper increase in the volume fraction of martensite from time step 1 to 2, which indicates that a more pronounced hardening is predicted by the ZM model leading to higher maximum value of the martensite volume fraction z at time step 2 that persists in subsequent load steps. In addition, the curve predicted by the ZM model features flat portions that indicate a delay in martensite reorientation with respect to changes in the loading direction. The delay is related to the formulation of the phase transformation function, which decreases slightly when the shear or tensile strain components decrease at the beginning of load steps 3, 4 and 8 (time steps 2-3, 3-4 and 7-8 respectively). The flat portions are indeed accompanied with a decrease of the von Mises equivalent stress, which strongly influences the phase transformation functions. Such decrease is not observed in the simulations of [START_REF] Stebner | Explicit finite element implementation of an improved three dimensional constitutive model for shape memory alloys[END_REF], which results in the absence of flat plateaus.

Strain-controlled tension-torsion loading

Experimental data reported by [START_REF] Grabe | Path dependence and multiaxial behavior of a polycrystalline NiTi alloy within the pseudoelastic and pseudoplastic temperature regimes[END_REF] shows a strong dependence of the material behavior on the loading history. The loading considered by the authors is the one shown in figure 12 The experiments by [START_REF] Grabe | Path dependence and multiaxial behavior of a polycrystalline NiTi alloy within the pseudoelastic and pseudoplastic temperature regimes[END_REF] ZM model with the parameters listed in table 4. Since the original ZM model is used, in which asymmetry in the material behavior is not considered, the simulations could not predict the tensile-compressive asymmetry seen in figure 12(b). The overall agreement with the experimental data remains satisfactory nevertheless and does produce a trend similar to the one observed in the figure 12(b), where the first load path results in shear stresses with higher magnitudes in the second and third quadrants and lower magnitudes in the first and fourth quadrants compared to the shear stresses obtained with the second load path for extremal values γ ′ = ±0.015 of the shear strain. A good fit is also achieved for the stress-strain behavior in shear shown in figure 12(d). The difference between the experimental and simulation results is more pronounced, however, for the axial stress-strain curves in figure 12(c).

Stress-controlled tension-torsion loading

Experimental data from [START_REF] Sittner | Experimental study on the thermoelastic martensitic transformation in shape memory alloy polycrystal induced by combined external forces[END_REF] is used to validate the model for the case of nonproportional stress-controlled loading in tension and torsion. The experiments were carried out on thin-walled cylindrical SMA samples. The authors used the following expressions for the equivalent stress and strain in reporting their results:

σ eq = σ 2 + τ 2 inv τ inv = C S τ, (44)
ε eq = ε 2 + γ 2 inv γ inv = γ C E , (45)
where τ inv and γ inv are energetically equivalent shear stress and strain and C S and C E are empirical coefficients. von Mises equivalent stress and strain are obtained for C S = C E = √ 3. The above expressions for the equivalent stress and strain do not appear to be frame-invariant for the choice of C S considered in the referenced paper. Since different expressions are used for σ eq and ε eq in the ZM model, two different loading cases are considered for the finite element simulations:

1. The first loading case is shown in figure 14(a), in which the normal and shear stress used for the simulation are exactly those reported in [START_REF] Sittner | Experimental study on the thermoelastic martensitic transformation in shape memory alloy polycrystal induced by combined external forces[END_REF], but the equivalent stress is not the same.

2. The second loading case is shown in figure 16(a), in which the shear stress used for the simulation produces an equivalent stress identical to the one in [START_REF] Sittner | Experimental study on the thermoelastic martensitic transformation in shape memory alloy polycrystal induced by combined external forces[END_REF].

For the first loading case, the material parameters used for the simulations are listed in table 5.

The parameters provide an acceptable fit of the uniaxial SMA response in tension and in shear as shown in figure 13. The results are in good agreement overall, with the model predictions being particularly consistent with those of [START_REF] Arghavani | A 3-D phenomenological constitutive model for shape memory alloys under multiaxial loadings[END_REF]. This is likely because both the ZM model and the model by [START_REF] Arghavani | A 3-D phenomenological constitutive model for shape memory alloys under multiaxial loadings[END_REF] propose similar loading functions and evolution rules for the inelastic deformation of SMAs. Because of the intrinsically different expressions used for the equivalent stress and strain, it is not possible to accurately fit the material behavior in tension and shear simultaneously to those reported by Sittner et al.

In order to validate the model for the case of nonproportional loading, the loading shown in figure 14 in tension and shear is shown in figures 14(c) and 14(d) respectively. In these figures, the results of numerical simulations using the ZM model are in good agreement overall with the experimental data and with the predictions of the other models considered. From figure 14(b), the material response appears to remain elastic throughout the tensile loading step and starts to become inelastic during the subsequent shear loading step where a deviation from the vertical direction is observed at an approximate shear strain of 0.005. In addition, it is interesting to observe that the predicted axial strain continues to increase during tensile unloading at constant shear stress, which is inconsistent with the experimental curve. A similar observation can also be made for the results obtained by [START_REF] Panico | A three-dimensional phenomenological model for martensite reorientation in shape memory alloys[END_REF]. The stressstrain behavior of the material in tension and in shear is shown in figures 14(c) and 14(d) respectively. In figure 14(c), it is interesting to see that the experimental curve appears to be nonlinear during the final unloading step, suggesting either an inelastic or a nonlinear elastic behavior of the material, which is not predicted by any of the models. point located on the external wall of the sample. The curves present notable differences that result from the relatively low radius-to-thickness ratio of the sample used by [START_REF] Bouvet | Mechanical Behavior of a Cu-Al-Be Shape Memory Alloy Under Multiaxial Proportional and Nonproportional Loadings[END_REF], which is 4.67. The assumption of a thin-walled cylinder may not therefore be adequate. The stress gradient across the thickness is indeed not negligible as shown in figure 20.

Conclusion

Numerical integration of the ZM model was presented in detail for the general case of SMAs undergoing phase transformation and martensite detwinning and reorientation. The proposed derivations include an analytical

3.

  Boundary value problem and algorithmic setup 3.1. Boundary value problem A superelastic SMA structure occupying a volume Ω ∈ R 3 is subjected to isothermal loading in the time interval [0, T ]. The load consists of body forces f i (x, t) acting over Ω × [0, T ] and contact forces T d i (x, t) acting over ∂Ω T × [0, T ]. Displacement is constrained over the remaining boundary ∂Ω ξ in the time interval [0, T ]. Starting from an initial configuration at time t = 0 in which the state variables are known and using the subscript notation ,i

  ); (b) Else if (∆z n+1 < 0 and F 1,(0) z,n+1 > 0 and η n+1 > 0) or (∆z n+1 > 0 and F 2,(0) z,n+1 > 0 and η n+1 > 0), then (martensite reorientation only): -Reset the state variables to their values at the beginning of the increment, -Calculate the inelastic multiplier η n+1 using (30), -Update the state variables, -Calculate the material Jacobian J ijkl using (43); (c) Else if (∆z n+1 > 0 and F 1,(0) z,n+1 > 0 and η n+1 < 0) or (∆z < 0 and F 2,(0) z,n+1 > 0 and η n+1 < 0, then (phase transformation only):

Figure 1 :

 1 Figure 1: Stress-controlled proportional loading cases.

Figure 2 :

 2 Figure 2: Simulated stress-strain behavior for pure tension and pure shear loading.

Figure 3 :

 3 Figure 3: Simulated stress-strain behavior and shear versus axial strain for tensiondominated tensile-shear loading (loading path 3).

Figure 4 :

 4 Figure 4: Simulated stress-strain behavior and shear versus axial strain for sheardominated tensile-shear loading (loading path 4).

Figure 5 :Figure 6 :

 56 Figure 5: Proportional tension-internal pressure loading.

1.Figure 7 :

 7 Figure 7: Simulation of martensite reorientation under multiaxial loading: (a) axial-shear loading path; (b) evolution of the martensite volume fraction and orientation strain.

Figure 8 :Figure 9 :

 89 Figure 8: Biaxial non-proportional loading: (a) stress-controlled axial-shear loading, (b) shear versus axial strain response.

Figure 10 :Figure 11 :

 1011 Figure10: Numerical prediction of the stress and strain behavior of the SMA and comparison with numerical simulation results from[START_REF] Stebner | Explicit finite element implementation of an improved three dimensional constitutive model for shape memory alloys[END_REF].

  (a), which consists in two strain-controlled butterfly loading paths that start and end at the same unloaded state. The maximum axial and shear strains reached in each case are ε = γ ′ = γ/ √ 3 = 0.015. Path 2 follows the loading sequence for path 1 in reverse order, such that the same intermediate states of shear and axial strain are achieved. The resulting material response in figure 12(b) shows a clear difference in terms of the shear versus axial stress curves corresponding to the two loading cases.

Figure 12 :

 12 Figure 12: Comparison of experimental and simulation results for the butterfly straincontrolled loading in figure 12(a).

  tner et al. (1995) as well as numerical simulations by[START_REF] Arghavani | A 3-D phenomenological constitutive model for shape memory alloys under multiaxial loadings[END_REF];[START_REF] Panico | A three-dimensional phenomenological model for martensite reorientation in shape memory alloys[END_REF] are reported on the same figures for comparison.

  Figure 13: Tensile and shear stress-strain curves used to determine the parameters of the ZM model.

  Figure 14: Experimental versus simulation results for the square tension-shear loading in figure 14(a).

4. 3 Figure 16 :Figure 17 :

 31617 Figure 16: Experimental versus simulation results for the loading shown in figure 16(a).

Figure 18 :

 18 Figure 18: Uniaxial tensile curve used to determine the parameters of the ZM model.

Figure 19 :

 19 Figure19: Simulation of the experimental data by[START_REF] Bouvet | Mechanical Behavior of a Cu-Al-Be Shape Memory Alloy Under Multiaxial Proportional and Nonproportional Loadings[END_REF].

Figure 20 :

 20 Figure 20: Stress variation in the radial direction for the cylinder used for numerical simulations.

Table 1 :

 1 used for the simulations are listed in table 1. Material Parameters used for the finite element simulations in section 4.1

	Parameter	Value	Parameter	Value
	E A	30340 MPa	a	5.16 MPa
	E M	18000 MPa	b	6.36 MPa
	ν	0.3	ε 0	0.04
	Y	30 MPa	G	13.17 MPa
	α	500 MPa	β	1250 MPa
	ξ	0.20 MPa	κ	4.16 MPa
	A 0 f	300 K	T	340 K

table 2 .

 2

		500
	P a ]	400
	τ [ M
	r e s s	300
	r S t
	h e a
	σ ; S	200
	r e s s
	l S t i a A x	100								Path 4: Tension
										Path 4: S hear
		0
		0	0 .	02	0	. 04	0 .	06	0 . 08	0 .1	0 .1 2	0 .1	4
						Axial Strain ε; Shear Strain γ

Table 3 :

 3 Parameters used for simulating SMA behavior in presence of nonproportional loading following[START_REF] Stebner | Explicit finite element implementation of an improved three dimensional constitutive model for shape memory alloys[END_REF].

Table 4 :

 4 Parameters used for the simulations following[START_REF] Grabe | Path dependence and multiaxial behavior of a polycrystalline NiTi alloy within the pseudoelastic and pseudoplastic temperature regimes[END_REF].

	are simulated using the

  The experimental results of Sit-

	Parameter	Value	Parameter	Value
	E A	30340 MPa	a	1.84 MPa
	E M	18000 MPa	b	1.3950 MPa
	ν	0.3	ε 0	0.02
	Y	30 MPa	G	11.46 MPa
	α	1000 MPa	β	2500 MPa
	ξ	0.116 MPa	κ	1.34 MPa
	A 0 f	300 K	T	340 K

Table 5 :

 5 Material parameters used for the simulation of nonproportional loading of SMAs following[START_REF] Sittner | Experimental study on the thermoelastic martensitic transformation in shape memory alloy polycrystal induced by combined external forces[END_REF].

Table 6 :

 6 Material parameters used for simulating the experimental data in[START_REF] Sittner | Experimental study on the thermoelastic martensitic transformation in shape memory alloy polycrystal induced by combined external forces[END_REF] for the loading case in figure16(a).

Acknowledgements Prof. Ziad Moumni would like to thank the 1000 Talent Plan and Northwestern Polytechnical University for their financial support. Dr. Wael Zaki would like to acknowledge the financial support of Khalifa University through KUIRF level 2 research fund no. 210031. Prof. Weihong Zhang would like to thank the National Natural Science Foundation of China (51221001), the National Basic Research Program of China (2012CB025904) and the 111 Project (B07050) for the financial support.


Fanny Frauziols
email: frauzi-ols@emse.fr

Jérôme Molimard

Laurent Navarro

Pierre Badel

Magalie Viallon

Rodolphe Testa

Prediction of the Biomechanical Effects of Compression Therapy by Finite Element Modeling and Ultrasound Elastography

Keywords: Medical Compression Stocking, Multi-modality, Soft tissues, Supersonic Shear Imaging

Goal: In the present study, the biomechanical response of soft tissues from the fascia cruris to the skin is studied in the human leg under elastic compression. Methods:

The distribution of elastic moduli in these tissues is measured for a volunteer at inactive and active muscle states using transient ultrasound elastography (TUSE). After registering the elasticity maps against magnetic resonance imaging (MRI) scans of the same volunteer, patient-specific finite element (FE) models are developed for the leg cross-section at inactive and active muscle states. Elastic properties obtained with TUSE are assigned at each Gauss point of the models. The response to 20 mmHg elastic compression is eventually predicted with the models. Results: Results show significantly higher elastic moduli in the fascia cruris tissue and also a significant increase of elastic moduli at active muscle state. Conclusion: This seems to have a marginal impact on pressure maps in the soft tissues of the leg predicted by the FE models. There is still an effect on the reduction of vein diameter induced by elastic compression, which is decreased at active muscle state. Significance: The discussion of the paper highlights the benefits of using elastography to reconstruct patient-specific FE models of soft tissues.

I. INTRODUCTION

A. Medical condition: causes and treatment D IFFERENT draining mechanisms prevent blood from stagnating in the legs of healthy subjects, such as the effect of muscular contraction, the activity of the sympathetic nervous system or the venous valve system [START_REF] Marieb | Anatomie et physiologie humaine[END_REF]. In case of default in one of these mechanisms, chronic venous insufficiency (CVI) may appear. Symptoms of CVI range from the feelings of "heavy leg" to ulcers. CVI affects about thirty percent of the French population [START_REF] Floury | Données sociales 1996 : la société franc ¸aise[END_REF]. Elastic compression of the limb using medical compression stockings (MCS) is a widely used treatment since it is efficient for almost all forms of CVI and can be employed at all stages of the disease [START_REF] Partsch | Chapter 10 -mechanism and effects of compression therapy[END_REF]. Fig. 1. Anatomy of the human leg -Cross-section of the human right leg. The fascia cruris makes a separation between the superficial tissues, composed by the adipose tissue, the skin and the superficial veins, and the deep tissues, composed by the muscular compartments and the deep veins. Adapted from [START_REF] Braus | Anatomie des menschen: ein lehrbuch für studierende und ärzte[END_REF].

B. Motivations for looking at superficial tissues

The benefits of MCS are not to be proven, but still, for some patients, the therapeutic goal is not reached. It has been shown that the ability of the treatment to achieve the desired medical goal is closely related to the question of transmission of pressure through the soft tissues [START_REF] Avril | Mixed experimental and numerical approach for characterizing the biomechanical response of the human leg under elastic compression[END_REF], as a result of which biomechanical models have been developed [START_REF] Dubuis | Identification of the material parameters of soft tissues in the compressed leg[END_REF], [START_REF] -Y. Rohan | Biomechanical response of varicose veins to elastic compression: a numerical study[END_REF]. These models investigated the transmission of pressure applied with compression garments [START_REF] -Y. Rohan | Biomechanical response of varicose veins to elastic compression: a numerical study[END_REF]- [START_REF] Wang | Analysis of flow and wall shear stress in the peroneal veins under external compression based on real-time mr images[END_REF]. Early studies focused on finding material properties to input in the numerical models [START_REF] Avril | Mixed experimental and numerical approach for characterizing the biomechanical response of the human leg under elastic compression[END_REF], [START_REF] Bouten | Identification des propriétés mécaniques des tissus constitutifs du mollet pour l'étude mécanique de la contention[END_REF], [START_REF] Avril | Patient-specific modeling of leg compression in the treatment of venous deficiency[END_REF]. The pressure was found to be heterogeneous inside the leg and with a great inter-subject variability. This gives a first statement on the subject-specific efficacy shown in the literature. Then, studies started to consider the deformation of veins in these models [START_REF] Dai | The effects of external compression on venous blood flow and tissue deformation in the lower leg[END_REF], [START_REF] Narracott | A validated model of calf compression and deep vessel collapse during external cuff inflation[END_REF], [START_REF] Wang | Finite element analysis of the deformation of deep veins in the lower limb under external compression[END_REF]. Recent studies [START_REF] Rohan | Prediction of the biomechanical effects of compression therapy on deep veins using finite element modelling[END_REF] have shown that the biomechanical action and the benefits of elastic compression are different for deep veins and for superficial veins (Fig. 1). Indeed, at active muscle state, the intramuscular pressure in the tissues surrounding deep veins varies by range ten times larger than the pressure applied by MCS. The direct impact of elastic compression on the caliber of deep veins is then relatively limited, especially in the standing position [START_REF] Rohan | Prediction of the biomechanical effects of compression therapy on deep veins using finite element modelling[END_REF]- [START_REF] Partsch | Narrowing of leg veins under compression demonstrated by magnetic resonance imaging (mri)[END_REF]. The situation is different for superficial veins. When the transmural pressure increases in these veins, their dilatation is only limited by their elastic properties and the elastic properties of the surrounding tissues.

C. Problem definition

The focus should be on the superficial veins network of the leg when looking at CVI treated by MCS. Also, to account for the inter-subject variability, the mechanical properties of each subject should be used. To that extent, this paper answers the two following questions:

• How do the elastic properties of soft tissues vary from the skin surface to muscle compartments, including the fascia cruris (Fig. 1)? • Does muscular activation have an impact on the elastic properties of superficial soft tissues?

D. Proposed solution by the use of a new technology

Elastography is a non-invasive way of measuring in vivo mechanical properties by inducing a static compression or transient vibrations in soft tissues [START_REF] Bensamoun | Determination of thigh muscle stiffness using magnetic resonance elastography[END_REF]. A new ultrasound based technology called Supersonic Shear Imaging (SSI) was introduced in 2004 [START_REF] Bercoff | Supersonic shear imaging: a new technique for soft tissue elasticity mapping[END_REF], for real-time soft tissue elasticity mapping. This transient elastography technique is used here for the first time to identify the material property distribution of the soft tissues from the fascia cruris to the skin in the human leg. After registering the obtained elasticity maps against magnetic resonance imaging (MRI) scans, patient-specific finite element (FE) models are developed for entire leg cross-sections at inactive and active muscle states. The response to 20 mmHg elastic compression is eventually predicted with these models.

E. General outlines of the paper

The material and methods section of this paper give a thorough explanation of the methodology for registering both imaging modalities, and for the implementation of the FE simulations. Then, in the result section, the distributions of the elastic moduli inside the leg at active and inactive muscle states are shown and the responses of the leg under 20 mmHg of elastic compression for different states of the leg are presented. Finally, these results and the limitations of this study are explained in the discussion section.

II. MATERIALS AND METHODS

A. Data acquisition and processing

A 22 year old healthy female volunteer underwent a MRI scan (Fig. 2) and SSI ultrasound elastography, following informed consent. Before acquiring ultrasound data, the subject was asked to perfom an isometric maximum voluntary contraction (MVC), and the response of a torque sensor, attached to the foot, was reported. 10% of this MVC was taken as baseline for contraction, to ensure good repeatability. An in-house program coded in LabVIEW R (National Instruments, Austin, TX, USA) was used to display the forces and give the volunteer visual feedback to ensure that the desired force was maintained when acquiring ultrasound data. Fig. 2. MRI -Fat-saturated multi-planar of the whole legs of the subject obtained from a high resolution 3D fast dual-echo sequence (Dixon [START_REF] Dixon | Simple proton spectroscopic imaging[END_REF], [START_REF] Lee | Fatty infiltration of the liver: demonstration by proton spectroscopic imaging. preliminary observations[END_REF]). Imaging parameters were as follows: repetition time (TR) 6.94 ms 2.45 ms 3.675 ms 5.69 ms, echo time (TE) 2.45 ms / 3.675 ms, flip angle 9 degrees, field of view (FOV) 330 mm x 330 mm, acquisition matrix 320 x 320, number of partitions 176, bandwidths 500/820 Hz/pixel, parallel acquisition acceleration factor 2 (grappa technique). The scan time per step was 3 min 24 sec. Subjects were imaged in supine feet first position with the arms along the body. Care was taken to avoid any compression of the leg muscle on the MR table when positioning the subject. Fig. 3. Experimental Setup -The probe is mounted on a ring, which allows collecting images around the leg. The foot is fixed through a shoe to a torque sensor. A feedback of the torque is given to the subject with the monitor. Data, including B-mode scans and SWE maps are acquired with Aixplorer R (SuperSonic Imagine, Aix-en-Provence, France).

An ultrasound elastographic device (Aixplorer R ; SuperSonic Imagine, Aix-en-Provence, France) was used to acquire the distribution of elastic moduli in the soft tissues of the leg. In the rest of the article, a distinction will be made between B-mode scans and Shear Wave Elasticity (SWE) maps. The B-mode scans (conventional echography) account for the anatomy of the leg, corresponding to the echogenicity of the medium displayed in gray scale. The SWE maps display the elastic moduli calculated from the velocity fields of shear waves inside the medium, displayed in a false color scale. An experimental setup was designed so as to hold the ultrasound probe at different positions all around the leg (Fig. 3). A special attention was paid to ensure that the plane of acquisition remains perpendicular to the tibial axis. To ensure minimal movement of the subject, her right shoe was attached to the plate. To ensure that the probe applies minimal compression to the tissues, a small gel bath was designed to place the probe at approximately half a centimeter from the skin. The probe is doing a 360 • rotation around the leg.

At each position around the leg, SWE and B-mode scans were simultaneously acquired, using a 4-15 MHz probe with 256 piezoelectric elements (SuperLinear T M SL15-4), first at inactive muscle state, second at active muscle state corresponding to 10% MVC. This required several consecutive acquisitions at different angular positions denoted β i , 1 ≤ i ≤ N (N = 37). Adjacent pairs of B-mode scans, say I i and I i-1 acquired at angles β i and β i-1 , were then stitched together using an inhouse program coded in Matlab R . After picking up a minimum of four anatomical landmarks (veins, crossing of aponeuroses) in each pair of adjacent images, the rigid transformation F 0 i relating the coordinate systems of both adjacent images was calculated by leastsquares minimization. Then, a cross-correlation factor was estimated across the overlapping area and this cross-correlation factor was optimized by updating the rigid transformation through the Nelder-Mead algorithm [START_REF] Nelder | A simplex method for function minimization[END_REF], yielding finally the rigid transformation F i . Repeating the procedure throughout 2 ≤ i ≤ N , all pairs of adjacent images were consecutively stitched together. Aligning the coordinate system of I 1 with the global coordinate system, the rigid transformation G i = F 2 × . . . × F i was applied to each image I i in order to align it within the global coordinate system. However, additional constraints had to be considered to ensure an optimal cross-correlation between I 1 and G N (I N ). Indeed, because the subject could move during the test, the images collected might have come from different transversal planes of the leg. Also, echography has a distortion coming from the width of the beam. Additionally, the stitching process was not fully automatic (points picked by user), thus, user dependent error could occur. To correct these errors, several options were considered (For details, see appendix A).

• Option 1. Find a rigid transformation F such that by applying

F × F i instead of F i to each I i , G N (I N ) = F × F 2 × . . . × F × F i × . . . × F × F N (I N )
maximizes the cross-correlation between I 1 and I N . • Option 2. Model the obtained unclosed contour of the leg as a curved beam and deform it until closure. This yields a set of transformations

F b i for each image I i ensuring that F b 2 × F 2 × . . . × F b i × F i × . . . × F b N × F N = ι with ι the identical transformation. • Option 3. Minimize a global factor: d = d(I 1 , G N (I N )) + N -1 i=1 d(G i (I i ), G i+1 (I i+1 ))
With d(I a , I b ) the objective function to be minimized representing the distance between anatomical landmarks supposed to be on top of each other after the transformation of successive images. Using the Nelder-Mead algorithm, d was minimized by updating the parameters of the G i rigid transformations, for all the transformations simultaneously.

For each angular position denoted β i , 1 ≤ i ≤ N , (N = 37), four SWE maps were acquired and averaged, yielding an elasticity image denoted J i . Then, each image was aligned within the global coordinate system by applying the rigid transformation: G i (J i )(G 1 = ι), permitting a stitched reconstruction of the whole distribution of elastic moduli.

B. FE model implementation

The MRI leg cross-section and the ultrasound reconstructions (B-mode and SWE) were registered using the bones edges (tibia and fibula) visible in both imaging modalities. From the skin segmentation of the B-mode reconstruction, the global geometry of the leg was defined. The contours and positions of the bones as well as the great and small saphenous veins [START_REF] Vin | Varices[END_REF] was added from the segmentation of the MRI leg cross-section. This geometry was used to create a 2D FE mesh, made of 4-node plane strain elements, using the meshing tools available in Abaqus R (Fig. 4). Soft tissues were represented with a Neo-Hookean hyper-elastic constitutive model. Their strain energy per unit of reference volume was defined such as [START_REF] Hibbit | Abaqus theory and user manuals version 6[END_REF]

U = C 10 ( Ī1 -3) + 1 D 1 (J el -1) 2
With Ī1 the first deviatoric strain invariant and C 10 and D 1 the material parameters reported in Table I. For each Gauss point from the fascia cruris to the skin, SWE data were available and the following material property was assigned:

C 10 = E 6
, where E was the locally measured elastic modulus using SWE. Deeper than the fascia cruris a value of C 10 = 3 kPa was assigned [START_REF] Dubuis | Identification of the material parameters of soft tissues in the compressed leg[END_REF] for the model at inactive muscle state. Regarding the model at 10% MVC, a parametric study was performed on the elastic moduli input deeper than the fascia cruris to evaluate the influence of likely evolving material properties in contracted muscles. The MCS was defined as a circle of radius equal to the radius of the sock at inactive muscle state. A class II MCS, corresponding to 15-20 mmHg external compression was used. Material properties are reported in Table I. The vein wall was modeled using five elements through the whole 1 mm thickness. The material properties were linearized following Rohan et al [START_REF] -Y. Rohan | Biomechanical response of varicose veins to elastic compression: a numerical study[END_REF]. The simulation entailed two steps. At first, the radius of the MCS was increased until the whole leg fitted in it. Then, this kinematic condition was removed and the contact between skin and MCS was activated. In order to match the experimental conditions where the subject lied supine, a blood pressure of 15 mmHg was applied on the inner surface of the vein walls. Both bones were fixed in the model. The interaction between the skin and the sock was enforced using a Coulomb frictional behavior with a skin-to-textile friction coefficient fixed at 0.3 as reported in literature [START_REF] Gerhardt | Skin-textile friction and skin elasticity in young and aged persons[END_REF]. From the Abaqus R library for contact algorithms, the penalty formulation was used in the tangential and normal directions. In the latter, an exponential contact pressure/over-closure relationship was defined. Two sets of simulations were computed:

• The biomechanical response of the leg to external compression at inactive and active (corresponding to 10% MVC) muscle states. • A parametric study to evaluate the influence of the constitutive parameters of the deep soft tissues on the response of superficial soft tissues to elastic compression.

III. RESULTS

A. B-mode assembly

A stitched reconstruction of the 37 B-mode scans without constraining the contour to be closed is shown in Fig. 5. A significant gap remains, which justifies the use of an additional constraint for the stitching approach. The results obtained with the different stitching approaches providing a closed contour are shown in Fig. 6. The sum of the differences of grayscale values in the overlapping areas were derived for the three options. Their values are 30.2 for option 1, 29.6 for option 2 and 29.4 for option 3. Option 3 is the approach which provides an optimal match between all pairs of images (Fig. 6 -c). The rigid transformations of option 3 are the ones used further to reconstruct the stitched elasticity maps for the leg at inactive and active muscle states (Fig. 7). Stitching errors were finally quantified through a criterion defined as the average distance between 10 anatomical landmarks in each pair of B-mode scans. Average errors for the entire stitched reconstruction are respectively 0.67 mm and 0.83 mm at inactive and active muscle states. The maps of elastic moduli reconstructed from SWE measurements are showed in Fig. 7 at inactive state and active muscle state corresponding to 10% MVC. They both are heterogeneous, with overall higher values at the surface (skin) and at the location of fascia cruris. Elastic moduli (Fig. 7 -b) are globally higher at active muscle state corresponding to 10% MVC than at inactive muscle state (Fig. 7 -a). A mean value of 12.4 kPa was measured for the elastic moduli (from skin to fascia cruris) at inactive muscle state and 17.2 kPa at active muscle state corresponding to 10% MVC.

B. Elastic compression of the leg

The first studied case is for the muscles at inactive muscle state (Fig. 8). Hydrostatic pressure was chosen for the analysis because it does not depend on the coordinate system and can be used to predict local fluid flows. The second studied case is for the muscles at 10% MVC. To the best of our knowledge, the literature does not report global C 10 value for 10% MVC inside the leg. To analyze the impact of the input value below the fascia cruris in our contracted model, a parametric study increasing this value up to 20 times the value inputted for the model at inactive muscle state was performed. Increasing the elastic moduli of the muscular compartments affects the hydrostatic pressure distributions derived from these simulations by accentuating the tension and compression found between the skin and the fascia cruris, until this inner C 10 overpasses a threshold value of 90 kPa. Then, the effect remains stable (Fig. 9). For all cases, the obtained maps (Fig. 8 and9) show heterogeneous distributions. The maximal pressure is located where the curvature radius is minimal (near the tibia) and conversely, minimal pressures are found on flat regions of the leg. An important effect of elastic compression is the impact on the area of the superficial veins. To study this effect, the FE model incorporated two major veins of the leg. For both veins (small and great saphenous veins), the effect of elastic compression is a slight decrease of the vein cross-sectional area (Tables II andIII). The decrease is significantly larger for the great saphenous vein than for the small saphenous vein, due to the larger curvature at the location of the great saphenous vein and due to the stiffer properties around the small saphenous vein. Also, the decrease in cross-sectional area diminishes as the elastic modulus increases inside the fascia cruris (Table III). the images. The quantitative spatial error on the stitching gave us a first quality criterion of this methodology. By assembling B-mode scans around the leg, a wide field of view was explored, comparable to the extended field-of-view [START_REF] Weng | US extended-field-of-view imaging technology[END_REF] (EFOV) software, used by clinicians to explore large areas, such as abdomen [START_REF] Kim | Extended field-of-view sonography advantages in abdominal applications[END_REF] or musculoskeletal system [START_REF] Lin | Extended field of view sonography in musculoskeletal imaging[END_REF]. However, to the best of our knowledge, extended field-of-view was never performed along with SWE.

Other imaging modalities can overcome this drawback. For example, magnetic resonance elastography has been used to explore elasticity of soft tissue in large anatomical areas such as thighs [START_REF] Bensamoun | Determination of thigh muscle stiffness using magnetic resonance elastography[END_REF] or livers [START_REF] Bensamoun | Measurement of liver stiffness with two imaging techniques: magnetic resonance elastography and ultrasound elastometry[END_REF] but due to a low spatial resolution, it is not appropriate to study sub-cutaneous tissues. Moreover, external vibrator is needed to generate low frequency shear waves compared to SSI where the ultrasonic probe remotely generates the shear wave inside the region of interest through pressure of radiation [START_REF] Bercoff | Supersonic shear imaging: a new technique for soft tissue elasticity mapping[END_REF] allowing ease of use and minimal discomfort to the subject.

To improve the SWE methodology, a different design of the probe that follows the contour of the leg may be an alternative in the future.

B. Material properties of the subcutaneous tissues

A mean value of 12.4 kPa was measured for the elastic moduli (from skin to fascia cruris) at inactive muscle state and 17.2 kPa at active muscle state corresponding to 10% MVC. These are consistent with literature, where values range from 11 to 32 kPa for fat [START_REF] Dubuis | Identification of the material parameters of soft tissues in the compressed leg[END_REF], [START_REF] Linder-Ganz | Assessment of mechanical conditions in sub-dermal tissues during sitting: a combined experimental-MRI and finite element approach[END_REF], [START_REF] Tran | In vivo characterization of the mechanical properties of human skin derived from MRI and indentation techniques[END_REF]. The two questions to be addressed in the introduction were here answered. The distribution of the elastic moduli in the superficial tissues are heterogeneous and the muscular activation tends to increase the values of the elastic moduli. The in vivo material properties of subcutaneous tissues have been an intense topic of research for the past 10 years. Measurement of skin elasticity have been performed through traditional tests such as suction [START_REF] Gerhardt | Skin-textile friction and skin elasticity in young and aged persons[END_REF], [START_REF] Hendriks | The relative contributions of different skin layers to the mechanical behavior of human skin in vivo using suction experiments[END_REF] or torsion. More recently, to account for anisotropy, micro-robot equipped with force and displacement sensors have been developed and used to demonstrate the strongly anisotropic and viscoelastic behavior of skin [START_REF] Flynn | Mechanical characterisation of in vivo human skin using a 3d force-sensitive micro-robot and finite element analysis[END_REF]. However, these techniques can only give local elasticity measures. Another approach combining finite element analysis, MRI and indentation was able to define more global mechanical parameters for the different layers of the human skin [START_REF] Tran | In vivo characterization of the mechanical properties of human skin derived from MRI and indentation techniques[END_REF] and subcutaneous tissues [START_REF] Linder-Ganz | Assessment of mechanical conditions in sub-dermal tissues during sitting: a combined experimental-MRI and finite element approach[END_REF]. All these different methods developed, although promising, cannot address the heterogeneities of mechanical properties found in subcutaneous soft tissues which can be obtained with SWE.

C. Prediction of the effects of elastic compression

Rohan et al. [START_REF] -Y. Rohan | Biomechanical response of varicose veins to elastic compression: a numerical study[END_REF] performed a computational parametric study and they showed that the biomechanical response of superficial veins to elastic compression is subject to three mechanical factors: the size of the cross-section, the position on the leg contour (local radius of curvature of the leg) and the "type" of fat of the patient. This highlighted the strong patient-specific response of the leg to external compression. By using the elastic moduli measured by SWE, the effect of heterogeneity in material property distribution and the change in elastic moduli of 10% MVC on pressure maps and veins contraction was determined in this study. This effect remains marginal, though, in the range of moduli measured for the volunteer considered here. This shows that the distribution of perivenous tissue pressure [START_REF] Gardon-Mollard | La compression médicale[END_REF] is mostly determined by the shape and the size of each subject's leg (local radii of curvature) for superficial veins such as the great saphenous veins and the small saphenous veins. The hydrostatic pressure distributions reported here confirm the idea that MCS work towards increasing perivenous tissue pressure on superficial veins consequently reducing the transmural pressure. Trans-mural pressure is also known to be related to the tension of the vein wall according to the Laplace law [START_REF] Gusic | Mechanical properties of native and ex vivo remodeled porcine saphenous veins[END_REF]. Reduction of this tension implies lower alterations in the vein wall associated with various pathologies.

Several researches [START_REF] Ramsey | Electromyographic and biomechanic analysis of anterior cruciate ligament deficiency and functional knee bracing[END_REF]- [START_REF] Kraemer | Effects of a whole body compression garment on markers of recovery after a heavy resistance workout in men and women[END_REF] have shown that muscle activation is modified by medical devices similar to MCS. Indeed, skin mechanoreceptors, provide not only exteroceptive by also proprioceptive effects, through Pacinian and Ruffini mechanoreceptors [START_REF] Johansson | Tactile sensory coding in the glabrous skin of the human hand[END_REF]. These are some of the receptors found in the subcutaneous tissue. They sustain essentially static pressure and stretch [START_REF] Martin | Modality coding in the somatic sensory system[END_REF]. In the skin and adipose tissue, inhomogeneous material properties around the locations of the mechanoreceptors may also cause significant differences in the spatial distribution of candidate relevant stimuli of the pressure. Therefore, the inhomogeneous material properties may have a stronger impact on the proprioceptive effects of elastic compression.

D. Limitations and perspectives

A question could be raised on the nature of elastic moduli given by SWE in soft tissues. Indeed, one assumption used to retrieve the elastic moduli from shear wave group velocity is that the medium is isotropic and homogenous. It is a fairly common assumption to consider adipose tissue as such [START_REF] Linder-Ganz | Assessment of mechanical conditions in sub-dermal tissues during sitting: a combined experimental-MRI and finite element approach[END_REF], [START_REF] Tran | In vivo characterization of the mechanical properties of human skin derived from MRI and indentation techniques[END_REF], [START_REF] Azar | Methods for modeling and predicting mechanical deformations of the breast under external perturbations[END_REF]. For skeletal muscular tissues, which are made of fibers, it has been previously considered as transversely isotropic [START_REF] Van Loocke | A validated model of passive muscle in compression[END_REF], [START_REF] Weiss | A constitutive model and finite element representation for transversely isotropic soft tissues[END_REF]. Here, a 2D model in transversal plane is constructed, so muscles can be considered as isotropic. Nevertheless, further investigation would be needed to check the relevance of elastic moduli given by the Aixplorer R on thin structures such as fascia cruris and skin. Apart from the assumptions used in SWE, a technical limitation constraining this methodology is not having a good access to deep tissue mechanical properties. This might be overcome by using a probe with a lower frequency that would provide useful information to study the transmission of pressure to deep veins. A higher level of muscular activation should be investigated to draw strong conclusions on the effects of muscular activation on the transmission of pressure.

V. CONCLUSION

In this paper, a new methodology coupling MRI and TUSE data was developed and implemented in a finite element model which is patient-specific in terms of geometry and material properties. This method was used to study the biomechanical response of the leg, under 20 mmHg of elastic compression. The focus was on the material properties of the superficial soft tissues from the skin to the fascia cruris, their impact on the hydrostatic pressure distribution inside the leg and the closure of the great and small saphenous veins. This opens the way to promising applications for predicting computationally the personalized action of biomechanical devices on soft tissues.

APPENDIX STITCHING ULTRASOUND IMAGES

Three options were considered to correct the problem of closure of the stitching of all the B-mode scans together.

• Option 1. Find a rigid transformation F such that by applying

F × F i instead of F i to each I i , G N (I N ) = F × F 2 × . . . × F × F i × . . . × F × F N (I N )
maximizes the cross-correlation between I 1 and I N .

Here, the error was to be systematic between all the images. The opening angle was divided by the number of images, and this correcting angle was used to correct the rotational component of the rigid body transformations, then the two other components of the rigid body transformations parameters (the two translations) were computed by maximizing the crosscorrelation between successive images. These three correcting components, the rotation and the two translations are represented by F in the equation. • Option 2. Model the obtained unclosed contour of the leg as a curved beam and deform it until closure. This yields a set of transformations F b i for each image

I i ensuring that F b 2 × F 2 × . . . × F b i × F i × . . . × F b N × F N = ι with ι the identical transformation.
Here, the error was not considered to be systematic and the idea was to use the finite element method to compute each individual error. Along the fascia cruris, a hundred points were manually picked on the unclosed stitching of the leg (Fig. 10 -a). The first and last points were taken as the same anatomical landmarks. Then, these points were used to define a curved beam on Abaqus R . Two boundary conditions were applied. The first boundary condition constrained one end of the line as fixed and the second boundary condition induced a displacement on the other end of the line to make both ends meet. Finally, the finite element simulation provided the displacements and rotations of all points (Fig. 10 -b). These displacements and rotation, F b i , were used to update each previously found rigid body transformations parameters to close the contour of the leg. Here, the goal was to find the rigid body transformations parameters globally instead of finding a correcting factor for the previous unclosed stitching. To do so, ten anatomical landmarks were picked on the image I i and the same ten anatomical landmarks were picked on the image I i+1 . This process was repeated through all the images and circled back to the first one, meaning that ten anatomical landmarks were picked on the image I N and the same ten anatomical landmarks were picked on the image I 1 (Fig. 11). Then, for each pairs of images, I i and I i+1 , the rigid body transformation parameters G i were computed and all the points picked were placed in the same global coordinate system. Following this initial estimate of G i , an objective function, d, was defined as the distance of each anatomical landmarks in the global coordinate system which are supposed to be on top of each other. Using the Nelder-Mead algorithm, d was minimized by updating the rigid body transformations parameters G i .

Fig. 4 .

 4 Fig. 4. Finite element mesh of the 2D patient-specific geometry (c) reconstructed from the segmentation of the transversal cross-section of the subject right leg acquired by MRI (a) for the bone and veins and from the ultrasound assembly (b) for the skin. It consists of continuum plane strain elements, except for skin and MCS, composed of truss elements. A relatively finer discretization is used in the vicinity of the great saphenous and small saphenous vein walls (c).

Fig. 5 .

 5 Fig. 5. Initial assembly of B-mode scans -The reconstruction is not closed (red circles).

Fig. 6 .

 6 Fig. 6. Closing contour methods comparison : a) option 1 b) option 2 c) option 3 optimization.

Fig. 7 .

 7 Fig. 7. Elastic moduli distribution in kPa at inactive muscle state (a) and at active muscle state corresponding to 10% MVC (b). The position of the fascia cruris is marked as dotted lines.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Hydrostatic pressure distribution in kPa for the leg at inactive muscle state, during 20 mmHg of elastic compression.

• Option 3 .

 3 Minimize a global factor:d = d(I 1 , G N (I N )) + N -1 i=1 d(G i (I i ), G i+1 (I i+1 ))With d(I a , I b ) the objective function to be minimized representing the distance between the anatomical landmarks placed in the global coordinate system which are supposed to be on top of each other. Using the Nelder-Mead algorithm, d was minimized by updating the parameters of the G i rigid transformations, for all the transformations simultaneously.

Fig. 10 .

 10 Fig. 10. Schematic overview of the second option of the stitching method. The contour of the aponeurosis is displayed on a stitching of B-mode scans (a). Pre-computation (green) and post-computation (blue) of the closure of the contour by finite element modeling (b).

Fig. 11 .

 11 Fig. 11. Schematic overview of the third option of the stitching method. The same anatomical landmarks are picked for each image of all the pairs of images.

TABLE III PARAMETRIC

 III STUDY. PERCENTAGE OF CONTRACTION (I.E., DECREASE OF THE CROSS-SECTIONAL AREA) OF THE GREAT AND SMALL SAPHENOUS VEINS UNDER AN ELASTIC COMPRESSION OF 15-20 MMHG FOR THE MODEL AT 10% MVC.

		Contraction (%)
		Great Saphenous Vein	Small Saphenous Vein
	a	8.1	5.7
	b	8.0	5.3
	c	7.8	4.6
	d	7.8	4.2
	e	7.7	3.9
	f	7.7	3.8
	g	7.7	3.6
	h	7.7	3.5
	i	7.7	3.3
	j	7.6	3.4
	k	7.6	3.3

ACKNOWLEDGMENT

The authors want to thank Pr. Rodolphe Le Riche for his support in developing the stitching methodology and the IVTV platform for its experimental support.

grateful to the C.N.R.S. (French Scientific National Research Center) for funding the Ph.D. grant of F.F and IFRESIS for its financial support.